最近在学习粒子群算法，看了很多资料都有点摸不着头脑。 在某博客中看到超简洁的粒子群c实现代码才知道粒子群算法的原理。 衷心感谢博主。 在此处发布博客的博客地址：
http://blog.sina.com.cn/s/blog_4ed027020100c9lx.html
在前面的文章中，梯度下降法、mdl法、善良的皮卡丘-mdl法、LM算法等属于优化算法。 这些优化算法的共同点是优化一组可能的解，使其尽可能接近真正的解。 粒子群算法(也称为PSO算法)也是一种优化算法，但该算法与上述算法的主要区别在于该算法同时优化了多个可能解，最后在多个可能解中选择最接近真实解作为最终解由于每个可能的解都是一个粒子，多个可能的解构成粒子群是算法名称的由来。 说明我对粒子群算法的理解，并将该算法应用于简单例子的优化求解。
1. 粒子群算法的核心思路
粒子群算法被提出的灵感来源于鸟群觅食。 如下图所示，在鸟群觅食的时候，各鸟向各个方向飞来觅食。 每只鸟都记得至今为止自己在飞行中最接近诱饵的位置。 另外，各鸟之间也有信息共享，比较与各自食物的最近距离，从各自的最近距离中，选择整体的一个最近距离的位置进行记忆。 每只鸟都随机地向各个方向飞行，所以各自的最近位置和整体最近位置都在不断更新。 也就是说，记忆中的最近位置接近食物，飞行充分到达后，记忆的整体最近位置也到达食物的位置。
抽象为数学模型，每只鸟是一个粒子，食物位置是问题的最优解，鸟与食物的距离是当前粒子的目标函数值。 例如，求出z=x*x y*y的最佳解，将粒子数设为6个时，如下图所示，相当于各粒子在曲面上滚动，各粒子I(0I5 )在自身滚动过程中的(反复)的最低位置)位置越低，z值越小
2. 粒子群算法的实现
假设目标函数有n个输入参数f(x1，x2，xn )。
对于每个粒子，它包含的元素
(1)当前时间位置(x1，x2，xn ) ) ) )。
)2)目前的目标函数值(也称为适应度(f ) x1，x2，xn ) )。
)3)该粒子的历史最佳位置(x1_pbest，x2_pbest，xn_pbest ) ) ) ) ) ) ) ) ) 652 )
)4)该粒子的历史最佳目标函数值(f_pbest=f(x1_pbest，x2_pbest，xn_pbest ) ) ) ) ) ) ) )。
所有粒子通用的元素：
(1)粒子总数) num
)2)总迭代次数： cnt。 粒子群算法也是一个迭代的过程，需要多次迭代才能得到理想的最优解。
(3)全局最佳位置(x1_gbest，x2_gbest，xn_gbest ) ) ) ) ) ) ) )全局最佳位置) ) ) ) 652 )
(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++代码实现（本代码参考了上述提到的博客的代码）：
定义全局参数：
# definedata _ size 500 const intnum=300； //粒子总数constfloat
c1 = 1.65; //速度计算参数1const float c2 = 1.65; //速度计算参数2const float xmin = -150; //位置下限const float xmax = 150; //位置上限const float vmin = -250.5; //速度下限const float vmax = 250.5; //速度上线定义粒子结构体：
获取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 = %lfn", 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");}运行上述代码，得到结果如下，可知获取的最优解还是非常接近真实解的。
计算速度时，分别加入惯性与不加入惯性，记录每轮迭代的全局最优目标函数值，如下图所示，可以看到加入惯性之后全局最优目标函数值减小得更快。