粒子群算法
粒子群算法，也称粒子群优化算法或鸟群觅食算法（Particle Swarm Optimization），缩写为 PSO， 是近年来由J. Kennedy和R. C. Eberhart等开发的一种新的进化算法(Evolutionary Algorithm - EA)。PSO 算法属于进化算法的一种，和模拟退火算法相似，它也是从随机解出发，通过迭代寻找最优解，它也是通过适应度来评价解的品质，但它比遗传算法规则更为简单，它没有遗传算法的“交叉”(Crossover) 和“变异”(Mutation) 操作，它通过追随当前搜索到的最优值来寻找全局最优。这种算法以其实现容易、精度高、收敛快等优点引起了学术界的重视，并且在解决实际问题中展示了其优越性。粒子群算法是一种并行算法。
1算法介绍
粒子群优化算法的基本思想是通过群体中个体之间的协作和信息共享来寻找最优解。
PSO的优势在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。
1.1问题提出
设想这样一个场景：一群鸟在随机的搜索食物。 在这个区域里只有一块食物，所有的鸟都不知道食物在哪。但是它们知道自己当前的位置距离食物还有多远。 那么找到食物的最优策略是什么？最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
1.2问题抽象
鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)。这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
1.3算法描述
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest,gbest)来更新自己。
在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。
i＝1，2，…，M，M是该群体中粒子的总数；Vi 是粒子的速度； pbest和gbest如前定义； rand()是介于（0、1）之间的随机数； xi是粒子的当前位置。c1和c2是学习因子，通常取c1＝ c2＝2 在每一维，粒子都有一个最大限制速度Vmax，如果某一维的速度超过设定的Vmax ，那么这一维的速度就被限定为Vmax 。（ Vmax >0）以上面两个公式为基础，形成了后来PSO 的标准形式。
1.4PSO算法流程
Step1：初始化一群微粒(群体规模为m)，包括随机位置和速度；
Step2：评价每个微粒的适应度；
Step3：对每个微粒，将其适应值与其经过的最好位置 pbest作比较，如果较好，则将其作为当前的最好位置pbest;
Step4：对每个微粒，将其适应值与其经过的最好位置 gbest作比较，如果较好，则将其作为当前的最好位置gbest;
Step5：根据(2)、(3)式调整微粒速度和位置；
Step6：未达到结束条件则转Step2。
迭代终止条件根据具体问题一般选为最大迭代次数Gk或和微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。
1.5参数分析
方程中pbest和gbest分别表示微粒群的局部和 全局最优位置，当c1＝0时，则粒子没有了认知能力， 变为只有社会的模型(social-only):
被称为全局PSO算法.粒子有扩展搜索空间的能力，具有 较快的收敛速度，但由于缺少局部搜索，对于复杂问题 比标准PSO 更易陷入局部最优。
当c2＝0时，则粒子之间没有社会信息，模型变为 只有认知(cognition-only)模型：
被称为局部PSO算法。由于个体之间没有信息的 交流，整个群体相当于多个粒子进行盲目的随机 搜索，收敛速度慢，因而得到最优解的可能性小。
群体规模m 一般取20～40，对较难或特定类别的问题 可以取到100～200。
最大速度Vmax决定当前位置与最好位置之间的区域的 分辨率(或精度)。如果太快，则粒子有可能越过极小 点;如果太慢，则粒子不能在局部极小点之外进行足 够的探索，会陷入到局部极值区域内。这种限制可以 达到防止计算溢出、决定问题空间搜索的粒度的目的。
学习因子c1和c2代表将每个粒子推向pbest 和gbest位置的统计加速项的权值。较低的值允许粒子 在被拉回之前可以在目标区域外徘徊，较高的值导致粒 子突然地冲向或越过目标区域。
2算法模型
目标函数：
约束条件：
3算法实现流程
4算法代码：
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define pop_size 20//种群规模
#define dim 2//粒子维度
#define PI 3.1415926 //圆周率
#define iteration_max 300//迭代次数
#define c2 2 //学习因子
#define c1 2 //学习因子
#define vmax 0.5//速度最大值
#define vmin -0.5//速度最小值
#define popmax 2//个体位置的最大值
#define popmin -2//个体位置的最小值
double pop[pop_size][dim];//种群位置的数组
double speed[pop_size][dim];//种群速度
double fitness[pop_size];//种群适应度数组
double gbest[dim];//群体极值的位置
double pbest[pop_size][dim];//个体极值的位置
double fitnesspbest[pop_size];//个体极值适应度值
double fitnessgbest;//群体极值适应度值
double iteration_best[iteration_max][dim]; //每一代最优值取值粒子
double result[iteration_max]; //定义存放每次迭代种群最优值的数组
double fitness_func(double *arr){
double x=*arr;
double y=*(arr+1);
double fitness = sin(sqrt(x*x+y*y))/(sqrt(x*x+y*y))+exp((cos(2*PI*x)+cos(2*PI*y))/2)-2.71289;
//printf("%lf\n",fitness);
return fitness;
}
void init()
{
for(int i=0;i<pop_size;i++)
{
for(int j=0;j<dim;j++)
{
pop[i][j] = (((double)rand())/RAND_MAX-0.5)*4;//-2到2之间的随机数
speed[i][j] = ((double)rand())/RAND_MAX-0.5; //-0.5到0.5之间
}
fitness[i]=fitness_func(pop[i]);
//printf("fit:%lf\n",fitness[i]);
}
}
double *max(double *fitness,int size)
{
int index=0;
double max=*fitness;
static double best_fit_index[2];
for(int i=0;i<size;i++)
{
if(*(fitness+i)>max)
{
max=*(fitness+i);
index=i;
}
}
best_fit_index[0]=index;
best_fit_index[1]=max;
return best_fit_index;
}
void pso(){
init();
double * best_fit_index; // 用于存放群体极值和其位置(序号)
best_fit_index = max(fitness,pop_size); //求群体适应度最大值
int index=(int)(*best_fit_index);
printf("best:%lf\n",*(best_fit_index+1));
//群体极值的位置
for(int i=0;i<dim;i++)
{
gbest[i]=pop[index][i];
}
//个体极值的位置
for(int i=0;i<pop_size;i++)
{
for(int j=0;j<dim;j++)
{
pbest[i][j]=pop[i][j];
}
}
//个体极值适应度值
for(int i=0;i<pop_size;i++)
{
fitnesspbest[i]=fitness[i];
}
//群体极值适应度值
fitnessgbest=*(best_fit_index+1);
//迭代寻优
for(int i=0;i<iteration_max;i++)
{
for(int j=0;j<pop_size;j++)
{
for(int k=0;k<dim;k++)
{
//更新速度
double rand1=(double)rand()/RAND_MAX; //0到1之间的随机数
double rand2=(double)rand()/RAND_MAX; //0到1之间的随机数
speed[j][k]=speed[j][k]+c1*rand1*(pbest[j][k]-pop[j][k])+c2*rand2*(gbest[k]-pop[j][k]);
if(speed[j][k]>vmax)
{
speed[j][k]=vmax;
}
if(speed[j][k]<vmin)
{
speed[j][k]=vmin;
}
//更新粒子的位置
pop[j][k]=pop[j][k]+speed[j][k];
if(pop[j][k]>popmax)
{
pop[j][k]=popmax;
}
if(pop[j][k]<popmin)
{
pop[j][k]=popmin;
}
}
//计算新粒子的适应度值
fitness[j]=fitness_func(pop[j]);
}
for(int j=0;j<pop_size;j++)
{
//更新个体极值的位置
if(fitness[j]>fitnesspbest[j])
{
for(int k=0;k<dim;k++)
{
pbest[j][k]=pop[j][k];
}
fitnesspbest[j]=fitness[j];
}
//更新群体极值的位置
if(fitness[j]>fitnessgbest)
{
for(int k=0;k<dim;k++)
{
gbest[k]=pop[j][k];
}
fitnessgbest=fitness[j];
}
}
for(int k=0;k<dim;k++)
{
iteration_best[i][k] = gbest[k]; // 每一代最优值取值粒子位置记录
}
result[i] = fitnessgbest; // 每代的最优值记录到数组
printf("iteration %d:\n",i);
printf("Optimal value particle position: x=%lf, y=%lf\n",iteration_best[i][0],iteration_best[i][1]);
printf("Optimal fitness=%lf\n",result[i]);
}
printf("--------------------------------------------------------\n");
printf("final result:\n");
printf("Optimal value particle position: x=%lf, y=%lf\n",gbest[0],gbest[1]);
printf("Optimal value=%lf\n",fitnessgbest);
}
int main()
{
int i=0;
srand((unsigned int)time(NULL));
pso();
return 0;
}
5结果
可以看出函数的最优解为1.005236。
6结果分析
求解函数的图形：
取第0、49、99、149、199、249、299代的20个粒子的位置做分析，可以看出粒子在迭代中进行收敛，最终在x=0，y=0处取得最优解。（图中横纵坐标分别为x、y）：
gnuplot绘图代码：
set xrange [-2:2] set yrange [-2:2] set size 1,1 set multiplot layout 2,2 set size 0.5,0.5 set origin 0,0.5 set title "0 iteration" plot '/home/unubtu/data0.dat' with points lc rgb 'blue' set origin 0.5,0.5 set title "49 iteration" plot '/home/unubtu/data1.dat' with points lc rgb 'blue' set origin 0,0 set title "99 iteration" plot '/home/unubtu/data2.dat' with points lc rgb 'blue' set origin 0.5,0 set title "149 iteration" plot '/home/unubtu/data3.dat' with points lc rgb 'blue' set xrange [-2:2] set yrange [-2:2] set size 1,1 set multiplot layout 2,2 set size 0.5,0.5 set origin 0,0.5 set title "199 iteration" plot '/home/unubtu/data4.dat' with points lc rgb 'blue' set origin 0.5,0.5 set title "249 iteration" plot '/home/unubtu/data5.dat' with points lc rgb 'blue' set origin 0,0 set title "299 iteration" plot '/home/unubtu/data6.dat' with points lc rgb 'blue' set isosamples 50,50 set pm3d set xlabel ‘x’ set ylabel ‘y’ set zlabel ‘f(x,y)’ splot [-2:2] [-2:2] sin(sqrt(x**2+y**2))/(sqrt(x**2+y**2))+exp((cos(2*pi*x)+cos(2*pi*y))/2)-2.71289