粒子群优化算法

粒子群优化算法属于群智能(swarm intelligence)优化算法。群智能分两种，一种是粒群优化，另一种是蚁群优化。

群智能概念

       假设你和你的朋友正在寻宝，每个人有个探测器，这个探测器可以知道宝藏到探测器的距离。你们一群人在找，每个人都可以把信息共享出去，就跟打dota时你可以有你队友的视野，你可以知道其他所有人距离宝藏的距离，这样，你看谁离宝藏最近，就向谁靠近，这样会使你发现宝藏的机会变大，而且，这种方法比你单人找要快的多。

       这是一个群行为(swarm behavior)的简单实例，群中各个体交互作用，使用一个比单一个体更有效的方法求解全局目标。可以把群(swarm)定义为某种交互作用的组织或Agent之结构集合，在群智能计算研究中，群的个体组织包括蚂蚁，白蚁，蜜蜂，黄蜂，鱼群，鸟群等。在这些群体中，个体在结构上是很简单的，而它们的集体行为却可能变得相当复杂。研究人员发现，蚂蚁在鸟巢和食物之间的运输路线，不管一开始多随机，最后蚂蚁总能找到一条最短路径。

粒群优化概念

粒群优化(particle swarm optimization,PSO)算法是一种基于群体搜索的算法，它建立在模拟鸟群社会的基础上。粒群概念的最初含义是通过图形来模拟鸟群优美和不可预测的舞蹈动作，发现鸟群支配同步飞行和以最佳队形突然改变飞行方向并重新编队的能力。这个概念已经被包含在一个简单有效的优化算法中。

在粒群优化中，被称为“粒子”(particle)的个体通过超维搜索空间“流动”。粒子在搜索空间中的位置变化是以个体成功地超过其他个体的社会心理意向为基础的，因此，群中粒子的变化是受其邻近粒子（个体）的经验或知识影响的。一个粒子的搜索行为受到群中其他粒子的搜索行为的影响。由此可见，粒群优化是一种共生合作算法。

算法描述

先通过一个形象的场景来描述一下：5只鸟觅食，每个鸟都知道自己与食物的距离，并将此信息与其他鸟共享。

一开始，5只鸟分散在不同的地方，假设没只鸟每秒钟更新自己的速度和方向，问题是怎么更新呢？

每只鸟记下自己离食物最近的位置，称为pbest，pbest0,pbest1,..分别表示5只鸟的pbest，从这里面选一个gbest，组里最好的。

每过一秒钟，鸟更新自己的速度v(此处为矢量)，

v_new  = v_old + c1*rand()*(pbest-pcurrent) +c2*rand()*(gbest-pcurrent)

c1,c2一般为2，rand()产生0~1的随机数。

然后以这个速度飞行1秒，再次更新，最终离食物越来越近。

 以下伪码摘自百度百科

　程序的伪代码如下

　　For each particle

　　____Initialize particle

　　END

　　Do

　　____For each particle

　　________Calculate fitness value

　　________If the fitness value is better than the best fitness value (pBest) in history

　　____________set current value as the new pBest

　　____End

　　____Choose the particle with the best fitness value of all the particles as the gBest

　　____For each particle

　　________Calculate particle velocity according equation (a)

　　________Update particle position according equation (b)

　　____End

　　While maximum iterations or minimum error criteria is not attained

　　在每一维粒子的速度都会被限制在一个最大速度Vmax，如果某一维更新后的速度超过用户设定的Vmax，那么这一维的速度就被限定为Vmax。

程序实例

计算f(x) = x*x - 20x + 100 的最小值及取最小值时的x

1. #include <iostream>  
2. #include <cmath>  
3. #include <cstdlib>  
4. using namespace std;  
5.  
6. #define C1 2  
7. #define C2 2  
8. #define VMAX 5.0  
9. #define MAX_ITERATIONS 100  
10. float rand01()  
11. {  
12.         return (float) (rand()/(double)RAND_MAX);  
13. }  
14. struct particle{  
15.         float current;  
16.         float pbest;  
17. };  
18.  
19. float fitness(float x)  
20. {  
21.         return x*x - 20*x + 100;  
22. }  
23.  
24. float gbest = 10000;  
25. struct particle p[5];  
26. float v[5] = {0};  
27.  
28. void init_particles()  
29. {  
30.         int i;  
31.         for(i = 0; i < 5; i++)  
32.         {  
33.                 p[i].current = -2+i;  
34.                 p[i].pbest = p[i].current;  
35.         }  
36. }  
37. void find_gbest()  
38. {  
39.         int i;  
40.         for(i = 0; i < 5; i++)  
41.         {  
42.                         if(fitness(gbest) > fitness(p[i].current))  
43.                                 gbest = p[i].current;  
44.         }  
45. }  
46. void adjust_v()  
47. {  
48.         int i ;  
49.         for(i = 0; i < 5; i++)  
50.         {  
51.                 v[i] = v[i] + C1*rand01()*(p[i].pbest - p[i].current) + C2*rand01()*(gbest - p[i].current);  
52.                 if(v[i] > VMAX)  
53.                         v[i] = VMAX;  
54.         }  
55. }  
56. void pso()  
57. {  
58.         int i,iter_num;  
59.         iter_num = 1;  
60.         while(iter_num < MAX_ITERATIONS)  
61.         {  
62.                 /*for(i = 0; i < 5; i++)  
63.                 {  
64.                         cout <<"p"<<i<<":current "<<p[i].current<<" pbest "<<p[i].pbest<<endl;  
65.                 }  
66.                 cout <<"gbest:"<<gbest<<endl;  
67.                 cout <<endl;  
68.                 getchar();*/ 
69.                 for(i = 0; i < 5; i++)  
70.                 {  
71.                         if(fitness(p[i].current) < fitness(p[i].pbest))  
72.                                 p[i].pbest = p[i].current;  
73.                 }  
74.                 find_gbest();  
75.                 adjust_v();  
76.                 for(i = 0; i < 5; i++)  
77.                         p[i].current += v[i];  
78.                 iter_num ++;  
79.         }  
80. }  
81. int main()  
82. {  
83.  
84.         init_particles();  
85.         pso();  
86.         printf("After %d iterations,gbest is %f\n",MAX_ITERATIONS,gbest);  
87.         return 0;  
88. }  

运行结果

下面是每次迭代后的结果，可以看到，经过5次迭代，结果就很好了。

1. After 1 iterations  
2. p0:current -2 pbest -2  
3. p1:current -1 pbest -1  
4. p2:current 0 pbest 0  
5. p3:current 1 pbest 1  
6. p4:current 2 pbest 2  
7. gbest:10000  
8.  
9.  
10. After 2 iterations  
11. p0:current 1.15506 pbest -2  
12. p1:current 3.79064 pbest -1  
13. p2:current 0.790205 pbest 0  
14. p3:current 2.53646 pbest 1  
15. p4:current 2 pbest 2  
16. gbest:2  
17.  
18.  
19. After 3 iterations  
20. p0:current 6.15506 pbest 1.15506  
21. p1:current 8.58128 pbest 3.79064  
22. p2:current 5.79021 pbest 0.790205  
23. p3:current 5.87216 pbest 2.53646  
24. p4:current 4.17373 pbest 2  
25. gbest:3.79064  
26.  
27.  
28. After 4 iterations  
29. p0:current 11.1551 pbest 6.15506  
30. p1:current 13.3719 pbest 8.58128  
31. p2:current 10.7902 pbest 5.79021  
32. p3:current 9.79741 pbest 5.87216  
33. p4:current 8.27141 pbest 4.17373  
34. gbest:8.58128  
35.  
36.  
37. After 5 iterations  
38. p0:current 13.8766 pbest 11.1551  
39. p1:current 10.1764 pbest 8.58128  
40. p2:current 14.7492 pbest 10.7902  
41. p3:current 13.7227 pbest 9.79741  
42. p4:current 13.2714 pbest 8.27141  
43. gbest:9.79741  
44.  
45.  
46. After 6 iterations  
47. p0:current 8.03327 pbest 11.1551  
48. p1:current 6.98078 pbest 10.1764  
49. p2:current 13.2414 pbest 10.7902  
50. p3:current 4.78856 pbest 9.79741  
51. p4:current 11.6974 pbest 8.27141  
52. gbest:10.1764  
53.  
54.  
55. After 7 iterations  
56. p0:current 5.84287 pbest 11.1551  
57. p1:current 9.25245 pbest 10.1764  
58. p2:current 5.23059 pbest 10.7902  
59. p3:current -3.28694 pbest 9.79741  
60. p4:current 9.93147 pbest 11.6974  
61. gbest:10.1764  

 本文转自nxlhero 51CTO博客，原文链接：http://blog.51cto.com/nxlhero/734212，如需转载请自行联系原作者