最近在学习粒子群算法，看了很多资料都有点摸不清头脑，直到看了一篇博客中超级简洁的粒子群C++实现代码，才明白粒子群算法的原理，真心感谢博主，在此贴出博主的博客地址：
http://blog.sina.com.cn/s/blog_4ed027020100c9lx.html
在前面的文章中，我们提到的梯度下降法、牛顿法、高斯-牛顿法，以及LM算法等都属于最优化算法，这些最优化算法的共同点是优化一组可能的解，使该解尽可能接近真实解。粒子群算法（也称为PSO算法）也属于最优化算法，但该算法与上述提到算法的主要区别在于，该算法同时优化多组可能的解，最后在多组可能的解中选择最接近真实解的一组作为最终解。每一组可能的解就是一个粒子，因此多组可能解组成了粒子群，这就是该算法名称的由来。下面将讲解我对粒子群算法的理解，并将该算法应用于一个简单例子的最优化求解。
1. 粒子群算法的核心思路
粒子群算法被提出的灵感来源于鸟群觅食。如下图，鸟群觅食过程中，每只鸟沿着各个方向飞行去寻找食物，每只鸟儿都能记住到目前为止自己在飞行过程中最接近食物的位置，同时每只鸟儿之间也有信息共享，它们会比较到目前为止各自与食物之间的最近距离，从各自的最近距离中，选择并记忆整体的一个最近距离位置。由于每只鸟儿都是随机往各个方向飞行的，各自的最近距离位置与整体最近距离位置不断被更新，也即它们记忆中的最近位置越来越接近食物，当它们飞行到达的位置足够多之后，它们记忆的整体最近位置也就达到了食物的位置。
抽象成数学模型，每只鸟儿就是一个粒子，食物的位置也就是问题的最优解，鸟儿与食物的距离也即当前粒子的目标函数值。比如要求z=x*x+y*y的最优解，设置粒子数为6个，如下图所示，相当于各个粒子在曲面上滚动，每个粒子i（0≤i≤5）记住自己滚动过程中（迭代）的最低位置（位置越低z值越小）Zi（0≤i≤5），同时不同粒子之间也有交流，它们会比较Z0~Z5，从中取得一个最小值作为当前的全局最低位置。
2. 粒子群算法的实现
假设目标函数有n个输入参数：f(x1, x2, ..., xn)。
对于每个粒子，它包含的元素：
(1) 当前时刻位置：(x1, x2, ..., xn)
(2) 当前时刻的目标函数值（也称为适应度）：f(x1, x2, ..., xn)
(3) 该粒子的历史最优位置：(x1_pbest, x2_pbest, ..., xn_pbest)
(4) 该粒子的历史最优目标函数值：f_pbest = f(x1_pbest, x2_pbest, ..., xn_pbest)
对于全部粒子，它们共有的元素：
(1) 粒子总数：num
(2) 迭代总次数：cnt。粒子群算法也是一个迭代的过程，需要多次迭代才能获取到理想的最优解。
(3) 全局最优位置：(x1_gbest, x2_gbest, ..., xn_gbest)
(4) 全局最优目标函数值：f_gbest = f(x1_gbest, x2_gbest, ..., xn_gbest)
(5) 位置随机化的上下限：xmin，xmax。迭代开始的时候以及迭代的过程中，均需要对粒子的位置进行随机分布，需要设置随机分布的上下限，不然随机分布偏离得太远，会严重影响优化结果。
(6) 速度的上下限：Vmin，Vmax。迭代过程中，速度也具有一定的随机性，需要限制速度的大小在一定范围内，不然如果速度值太大也会严重影响优化结果。
(7) 速度计算参数：c1、c2。通常取1.0~1.8的值。
粒子位置的更新公式：
参考上述提到的博客，对于每一个粒子的每一个位置参数xi，其位置更新公式如下，其中randf为0~1的随机数，c1与c2通常取1.0~1.8之间的固定值。
为了加快收敛速度，根据速度具有惯性的原理，后来人们提出了在速度计算中增加惯性，也即加上上一轮迭代时该位置参数的速度。如下式，其中w为0~1的参数，通常随着迭代次数的增加，逐渐减小w，因为越到后面可能解就越接近真实解，迭代收敛就越慢，所以需要减小w来放慢速度，否则容易错过最优解。
C++代码实现（本代码参考了上述提到的博客的代码）：
定义全局参数：
#define DATA_SIZE 500
const int NUM = 300; //粒子总数
const float c1 = 1.65; //速度计算参数1
const float c2 = 1.65; //速度计算参数2
const float xmin = -150; //位置下限
const float xmax = 150; //位置上限const float vmin = -250.5; //速度下限
const float vmax = 250.5; //速度上线
定义粒子结构体：
struct particle
{float x[DATA_SIZE]; //该粒子的当前时刻位置float bestx[DATA_SIZE];//该粒子的历史最优位置float f; //该粒子的当前适应度float bestf; //该粒子的历史最优适应度
}swarm[NUM]; //总共有num个粒子
获取0~1的随机数：
#define randf ((rand()%10000+rand()%10000*10000)/100000000.0) //产生0-1随机浮点数
目标函数，显然目标函数的真实解为(0, 1, 2，..., dim-1)：
float f1(float x[], int dim)
{float z = 0;for (int i = 0; i < dim; i++){z += (i+20)*(x[i] - i)*(x[i] - i);}return z;
}
迭代优化过程：
/*
best为全局最优位置
DIM为目标函数的输入参数总个数，也即每个粒子的位置参数总个数
*/
void PSO(float *best, int DIM)
{for (int i = 0; i < DIM; i++)//初始化全局最优{best[i] = randf*(xmax - xmin) + xmin; //随机生成xmin~xmax的随机数}//初始化全局最优目标函数值为一个非常大的值double gbestf = 100000000000000.0; for (int i = 0; i < NUM; i++) //初始化粒子群{particle* p1 = &swarm[i]; //获取每一个粒子结构体的地址for (int j = 0; j < DIM; j++){//随机生成xmin~xmax的随机数作为该粒子的位置参数p1->x[j] = randf*(xmax - xmin) + xmin;}p1->f = f1(p1->x, DIM); //计算当前时刻的目标函数值p1->bestf = 100000000000000.0; //初始化每个粒子的历史最优目标函数值为一个非常大的值}float *V = (float *)calloc(DIM, sizeof(float));const int cnt = 2000; //总共2000次迭代float w = 0.0025/(cnt-1); //设置w初值，后面逐渐减小for (int t = 0; t < cnt; t++) //cnt次迭代{for (int i = 0; i < NUM; i++) //NUM个粒子{particle* p1 = &swarm[i];for (int j = 0; j < DIM; j++) //计算速度，并更新位置{//w*(cnt-1-t)随着t的增加逐渐减小V[j] = w*(cnt-1-t)*V[j] + c1*randf*(p1->bestx[j] - p1->x[j]) + c2*randf*(best[j] - p1->x[j]);//V[j] = c1*randf*(p1->bestx[j] - p1->x[j]) + c2*randf*(best[j] - p1->x[j]);//限制速度的上下限V[j] = (V[j] < vmin) ? vmin : ((V[j] > vmax) ? vmax : V[j]);p1->x[j] = p1->x[j] + V[j]; //更新位置参数}p1->f = f1(p1->x, DIM); //计算当前粒子的当前时刻目标函数值//如果当前粒子的当前时刻目标函数值小于其历史最优值，则替换该粒子的历史最优if (p1->f < p1->bestf) {for (int j = 0; j < DIM; j++){p1->bestx[j] = p1->x[j];}p1->bestf = p1->f;}//如果当前粒子的历史最优值小于全局最优值，则替换全局最优值if (p1->bestf < gbestf) {for (int j = 0; j < DIM; j++){best[j] = p1->bestx[j];}//这里有点不理解，为什么替换全局最优值之后要重新随机分布该粒子的位置参数？//如果没有这部分代码，整体优化效果就会很差for (int j = 0; j < DIM; j++) {p1->x[j] = randf*(xmax - xmin) + xmin;}gbestf = p1->bestf;printf("t = %d, gbestf = %lf\n", t, gbestf);}}}free(V);
}
主函数：
void main(void)
{int DIM = 400;//维数float gbestx[DATA_SIZE];//全局最优位置PSO(gbestx, DIM);printf("全局最优位置:\n");for(int i = 0; i < DIM; i++){printf("%f ", gbestx[i]);if (i > 0 && i % 7 == 0)printf("\n");}printf("\n");
}
运行上述代码，得到结果如下，可知获取的最优解还是非常接近真实解的。
计算速度时，分别加入惯性与不加入惯性，记录每轮迭代的全局最优目标函数值，如下图所示，可以看到加入惯性之后全局最优目标函数值减小得更快。
