文章目录
- 算法
- 实例
- 求解极值 C 语言实现
- 求解极值 matlab 语言实现
- 画图
- PSO 求解
- 结论
PSO算法是经典的智能优化算法,在数学建模等比赛中非常常用,求解时的效果不错。对于智能优化算法,个人倾向于matlab实现,因为计算起来非常方便。但是这次因为老师的要求,准备C语言实现PSO算法求二元函数极值。说实话,第一次听到老师的要求时吓了一跳,因为从来没想过尝试过C语言实现这个算法,感觉这操作多少是有点阴间了。不过静下心来写写发现其实还是很容易的,可能是因为粒子群优化算法比较简单。
算法
下面是一个不太标准的流程图,个人认为认真读代码3分钟就能看懂(可以看C语言代码中的PSO_search
)。
这里浅浅讲一下自己对PSO算法的理解,个人认为PSO算法实际上就是有方向的暴力搜索,一开始瞎搜,然后粒子们逐渐朝表现最好的粒子的方向靠拢着搜索,最后会收敛于某个点。
粒子群优化算法里面有粒子的加速度和速度的最大最小值,个人认为这个东西有点像深度学习里的学习率,用来控制收敛的速度。大了可能在最优值附近摇摆不定,小了可能陷入局部最优出不来。
实例
待求极值函数如下:
求解极值 C 语言实现
C语言代码如下,主要的流程在 PSO_search
函数里写的比较清楚;两个题目可以在 fitting
函数里面切换,修改的时候记得改 #define
里的最大最小值。
#include<bits/stdc++.h>
using namespace std;
// 种群大小
#define population_size 300
// 维度=> 二元的,二维的
#define dimension 2
// 迭代次数
#define max_epoch 15
// 最大值最小值
#define population_max 10
#define population_min 0
// 定义速度
#define velocity_max 5
#define velocity_min -5
// 定义维数
#define dimension 2
// 网上说取这个很好
#define c1 1.49445
#define c2 1.49445
// 鸟群
double population[population_size][dimension];
// 鸟群速度
double velocity[population_size][dimension];
// 适应度,其实就是每个粒子对应的解
double fitness[population_size];
// 个体极值的位置
double pbest[population_size][dimension];
// 个体极值适应度的值
double fitnesspbest[population_size];
// 群体极值的位置
double gbest[dimension];
// 群体极值适应度值
double fitnessgbest;
// 存每次迭代种群后储存本次迭代最优值的数组=>存的是适应度
double result[max_epoch];
// 存种群每次迭代后最优值(适应度)时的x,y => 当然不是二维也可以,不一定时x, y
double fitnessepoch[population_size][dimension];
// 不让用库,只能手写它俩
double max(double x1, double x2){
if(x1>x2)
return x1;
else
return x2;
}
double min(double x1, double x2){
if(x1<x2)
return x1;
else
return x2;
}
// 适应度函数,其实就是目标函数
double fitting(double x, double y){
// 第一题:
// double fitness = sin(x)/x+sin(y)/y;
// 第二题
double fitness = (6.452*(x+0.125*y)*pow((cos(x)-cos(2*y)),2))/sqrt(0.8+pow((x-4.2),2)+2*pow((y-7),2))+3.226*y;
return fitness;
}
// 初始化函数,一开始还是随机初始化
void initialization(void){
for(int i = 0; i < population_size; i++){
for(int j = 0; j < dimension; j++){
// 初始化第一批粒子的位置和速度
population[i][j] = ((double)rand()/(RAND_MAX))*(population_max-population_min)+population_min;
velocity[i][j] = ((double)rand()/(RAND_MAX))*(velocity_max-velocity_min)+velocity_min;
}
// 顺便算出来第一批粒子的适应度函数
fitness[i] = fitting(population[i][0], population[i][1]);
// printf("%lf\t%lf\tfitness:%lf\n", population[i][0], population[i][1], fitness[i]);
}
}
// 作用其实就是比较所有散点对应的适应度,输出最大的适应度和对应的位置
double* max_in_community(double fit[], int size){
static double result[2];
result[0] = *fit;
result[1] = 0;
for(int i = 0; i < size; i++)
if(fit[i]>result[0]){
result[0] = fit[i];
result[1] = i;
}
// printf("%lf\n%lf\n",result[0], result[1]);
return result;
}
void PSO_search(void){
//——————————Step1——————————//
// 随机初始化种群
initialization();
//先初始化下第一轮的变量
double *best_fitness;
// 该轮所有鸟的最大值和index
best_fitness = max_in_community(fitness, population_size);
// 存群体的极值在哪里
for(int i = 0; i<dimension; i++)
gbest[i] = population[int(best_fitness[1])][i];
//——————————Step2——————————//
// 个体极值位置
for(int i = 0; i < population_size; i++)
for(int j = 0; j < dimension; j++)
pbest[i][j] = population[i][j];
// 个体极值
for(int i = 0; i < population_size; i++)
fitnesspbest[i] = fitness[i];
// 群体极值
// printf("%lf\n%lf\n", best_fitness[0], best_fitness[1]);
fitnessgbest = best_fitness[0];
//——————————Step3——————————//
// 开始迭代
for(int i=0; i<max_epoch; i++){
for(int j = 0; j<population_size; j++){
for (int k = 0; k<dimension; k++){
// 更新速度
// 这个是PSO的公式
double r1 = (double)rand()/RAND_MAX;
double r2 = (double)rand()/RAND_MAX;
velocity[j][k] = velocity[j][k] + c1*r1*(pbest[j][k]-population[j][k]) + c2*r2*(gbest[k]-population[j][k]);
// 防止速度超出范围
velocity[j][k] = min(velocity[j][k], velocity_max);
velocity[j][k] = max(velocity[j][k], velocity_min);
// 更新位置
population[j][k] += velocity[j][k];
population[j][k] = min(population[j][k], population_max);
population[j][k] = max(population[j][k], population_min);
}
// printf("(%lf,%lf)\n", population[j][0], population[j][1]);
// 计算更新后的fitness
fitness[j] = fitting(population[j][0], population[j][1]);
}
// 对于每一次迭代,都要更新极值(个体 群体)
for(int j = 0; j<population_size; j++){
if(fitness[j] > fitnesspbest[j]){
for(int k=0;k<dimension;k++)
pbest[j][k] = population[j][k];
fitnesspbest[j] = fitness[j];
}
// 对于每一次迭代,同样也需要更新群体极值
if(fitness[j] > fitnessgbest){
for(int k=0;k<dimension;k++)
gbest[k] = population[j][k];
fitnessgbest = fitness[j];
}
}
// 迭代完了, 记录下每次的结果
for(int k = 0; k<dimension; k++)
fitnessepoch[i][k] = gbest[k];
// 记录一下
result[i] = fitnessgbest;
printf("第%d次迭代:(%lf,%lf)%lf\n", i, fitnessepoch[i][0], fitnessepoch[i][1], result[i]);
}
}
int main(){
//程序开始和结束时间
clock_t begin, end;
begin = clock(); //开始计时
srand((unsigned)time(NULL)); // 初始化随机数种子
// 搜索
PSO_search();
double * solve;
solve = max_in_community(result,max_epoch);
printf("迭代了%d次,在第%d次收敛到最优值,最优解为:%lf\n", max_epoch, solve[1], solve[0]);
// 这次迭代中的最优解中找位置
printf("取到最优值的位置为(%lf,%lf).\n", fitnessepoch[int(solve[1])][0], fitnessepoch[int(solve[1])][1]);
end = clock();
printf("计算总耗时:%lf秒", end - begin);
return 0;
}
求解极值 matlab 语言实现
画图
不是很确定最后的结果对不对,拿matlab试试结果。先浅浅画俩图:
fsurf(@(x,y)sin(x)/x+sin(y)/y)
fsurf(@(x,y)(6.452*(x+0.125*y)*(cos(x)-cos(2*y))^2/(0.8+(x-4.2)^2+2*(y-7)^2)^(1/2))+3.226*y)
例题1:
例题2:
PSO 求解
matlab自带PSO求解工具计算极值,因为这个东西只能算最小值,所以对原函数取负号,得到的最小值再取负就是我们要的最大值:
[x,fval] = particleswarm(@(x)-(sin(x(1))/x(1)+sin(x(2))/x(2)),2,[-10,-10],[10,10])
[x,fval] = particleswarm(@(x)-(6.452*(x(1)+0.125*x(2))*(cos(x(1))-cos(2*x(2)))^2/(0.8+(x(1)-4.2)^2+2*(x(2)-7)^2)^(1/2)+3.226*x(2)),2,[0,0],[10,10])
运行结果如下:
例题1
例题2
结论
经过证明,我们C语言的实现是没什么大问题的。