PSO粒子群算法(鸟群算法)计算二元函数极值(C语言实现matlab工具箱实现)
Posted 铖铖的花嫁
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了PSO粒子群算法(鸟群算法)计算二元函数极值(C语言实现matlab工具箱实现)相关的知识,希望对你有一定的参考价值。
算法讲解和感悟
PSO算法是经典的智能优化算法,在数学建模等比赛中非常常用,求解时的效果不错。对于智能优化算法,个人倾向于matlab实现,因为计算起来非常方便。但是这次因为老师的要求,准备C语言实现PSO算法求二元函数极值。说实话,第一次听到老师的要求时吓了一跳,因为从来没想过尝试过C语言实现这个算法,感觉这操作多少是有点阴间了。不过静下心来写写发现其实还是很容易的,可能是因为粒子群优化算法比较简单。下面是一个不太标准的流程图,个人认为认真读代码3分钟就能看懂(可以看C语言代码中的PSO_search)。
这里浅浅讲一下自己对PSO算法的理解,个人认为PSO算法实际上就是有方向的暴力搜索,一开始瞎搜,然后粒子们逐渐朝表现最好的粒子的方向靠拢着搜索,最后会收敛于某个点。
粒子群优化算法里面有粒子的加速度和速度的最大最小值,个人认为这个东西有点像深度学习里的学习率,用来控制收敛的速度。大了可能在最优值附近摇摆不定,小了可能陷入局部最优出不来。
实例
待求极值函数如下:
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试试结果。先浅浅画俩图:
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:
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])
运行结果如下:
经过证明,我们C语言的实现是没什么大问题的。
以上是关于PSO粒子群算法(鸟群算法)计算二元函数极值(C语言实现matlab工具箱实现)的主要内容,如果未能解决你的问题,请参考以下文章
MATLAB教程案例11基于PSO粒子群优化算法的函数极值计算matlab仿真及其他应用