C语言实现粒子群算法(PSO)一
最近在温习C语言，看的书是《C primer Plus》，忽然想起来以前在参加数学建模的时候，用过的一些智能算法，比如遗传算法、粒子群算法、蚁群算法等等。当时是使用MATLAB来实现的，而且有些MATLAB自带了工具箱，当时有些只是利用工具箱求最优解问题，没有自己动手亲自去实现一遍，现在都忘的差不多了。我觉得那样层次实在是很浅，没有真正理解算法的核心思想。本着“纸上得来终觉浅，绝知此事要躬行”的态度，我决定现在重新复习一遍算法，然后手工用C语言重新实现一遍。说做就做，我第一个实现的算法是相对来说比较简单的粒子群算法(与遗传算法等相比，至少我自己觉得实现要简单一些)。
首先简单介绍一下启发式算法和智能算法。粒子群算法、遗传算法等都是从传统的搜索算法演变而来的启发式算法。启发式算法（heuristic algorithm)是相对于最优化算法提出的。一个问题的最优算法求得该问题每个实例的最优解。启发式算法可以这样定义：一个基于直观或经验构造的算法，在可接受的花费（指计算时间和空间）下给出待解决组合优化问题每一个实例的一个可行解，该可行解与最优解的偏离程度一般不能被预计，但是通常情况下启发式算法可以给出接近最优解的不错的解，但是无法保证每次它都可以得到很好的近似解。启发式算法中有一类被称之为智能算法，所谓"智能"二字，指的是这种算法是通过模仿大自然中的某种生物或者模拟某种现象而抽象得到的算法，比如遗传算法就是模拟自然界生物自然选择，优胜劣汰，适者生存而得到的进化算法，粒子群是源于对于鸟类捕食行为的研究，而模拟退火算法则是根据物理学中固体物质的退火过程抽象得到的优化算法。智能算法兴起于上个世纪80年代左右，之后就一直发展迅速，除了传统的智能算法之外，近几年又涌现出了一些新的算法比如鱼群算法、蜂群算法等。
言归正传，下面来介绍今天的主角：粒子群算法。粒子群算法的基本原理如下（参考《MATLAB智能算法30个案例分析》）：
假设在一个D维的搜索空间中，由n个粒子组成的种群X=(X1,X2,..,Xn),其中第i个粒子表示为一个D维的向量Xi=(xi1,xi2,xiD),代表第i个粒子在D维搜索空间中的位置，亦代表问题的一个潜在解。根据目标函数即可以计算出每个粒子位置Xi对应的适应度值。第i个粒子的速度为Vi = (Vi1,Vi2,...,ViD),其个体极值为Pi=(Pi1,Pi2,...,PiD),种群的群体极值为Pg=(Pg1,Pg2,...,PgD)。在每次迭代的过程中，粒子通过个体极值和群体极值更新自身的速度和位置，即：
Vid(k+1)=w*Vid(k)+c1*r1*(Pid(k)-Xid(k))+c2*r2*(Pgd(k)-Xid(k))
Xid(k+1) = Xid(k) + Vid(k+1)
其中w为惯性权重，如果不考虑可以默认为1,后面还会再详细讨论w对于PSO的影响。d=1,2,..,D;i=1,2,...,n;k为当前迭代次数;Vid为粒子的速度;c1和c2是非负的常数，称为加速度因子;r1和r2是分布于[0,1]之间的随机数。为了防止粒子的盲目搜索，一般建议将其位置和速度限制在一定的区间内。
下面是我用C语言实现的求一个二元函数最大值的粒子群算法：
1 /* 2 * 使用C语言实现粒子群算法(PSO) 3 * 参考自《MATLAB智能算法30个案例分析》 4 * update: 16/12/3 5 * 本例的寻优非线性函数为 6 * f(x,y) = sin(sqrt(x^2+y^2))/(sqrt(x^2+y^2)) + exp((cos(2*PI*x)+cos(2*PI*y))/2) - 2.71289 7 * 该函数有很多局部极大值点，而极限位置为(0,0),在(0,0)附近取得极大值 8 */ 9 #include<stdio.h> 10 #include<stdlib.h> 11 #include<math.h> 12 #include<time.h> 13 #define c1 1.49445 //加速度因子一般是根据大量实验所得 14 #define c2 1.49445 15 #define maxgen 300 // 迭代次数 16 #define sizepop 20 // 种群规模 17 #define popmax 2 // 个体最大取值 18 #define popmin -2 // 个体最小取值 19 #define Vmax 0.5 // 速度最大值 20 #define Vmin -0.5 //速度最小值 21 #define dim 2 // 粒子的维数 22 #define PI 3.1415926 //圆周率 23 24 double pop[sizepop][dim]; // 定义种群数组 25 double V[sizepop][dim]; // 定义种群速度数组 26 double fitness[sizepop]; // 定义种群的适应度数组 27 double result[maxgen]; //定义存放每次迭代种群最优值的数组 28 double pbest[sizepop][dim]; // 个体极值的位置 29 double gbest[dim]; //群体极值的位置 30 double fitnesspbest[sizepop]; //个体极值适应度的值 31 double fitnessgbest; // 群体极值适应度值 32 double genbest[maxgen][dim]; //每一代最优值取值粒子 33 34 //适应度函数 35 double func(double * arr) 36 { 37 double x = *arr; //x 的值 38 double y = *(arr+1); //y的值 39 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; 40 return fitness; 41 42 } 43 // 种群初始化 44 void pop_init(void) 45 { 46 for(int i=0;i<sizepop;i++) 47 { 48 for(int j=0;j<dim;j++) 49 { 50 pop[i][j] = (((double)rand())/RAND_MAX-0.5)*4; //-2到2之间的随机数 51 V[i][j] = ((double)rand())/RAND_MAX-0.5; //-0.5到0.5之间 52 } 53 fitness[i] = func(pop[i]); //计算适应度函数值 54 } 55 } 56 // max()函数定义 57 double * max(double * fit,int size) 58 { 59 int index = 0; // 初始化序号 60 double max = *fit; // 初始化最大值为数组第一个元素 61 static double best_fit_index[2]; 62 for(int i=1;i<size;i++) 63 { 64 if(*(fit+i) > max) 65 max = *(fit+i); 66 index = i; 67 } 68 best_fit_index[0] = index; 69 best_fit_index[1] = max; 70 return best_fit_index; 71 72 } 73 // 迭代寻优 74 void PSO_func(void) 75 { 76 pop_init(); 77 double * best_fit_index; // 用于存放群体极值和其位置(序号) 78 best_fit_index = max(fitness,sizepop); //求群体极值 79 int index = (int)(*best_fit_index); 80 // 群体极值位置 81 for(int i=0;i<dim;i++) 82 { 83 gbest[i] = pop[index][i]; 84 } 85 // 个体极值位置 86 for(int i=0;i<sizepop;i++) 87 { 88 for(int j=0;j<dim;j++) 89 { 90 pbest[i][j] = pop[i][j]; 91 } 92 } 93 // 个体极值适应度值 94 for(int i=0;i<sizepop;i++) 95 { 96 fitnesspbest[i] = fitness[i]; 97 } 98 //群体极值适应度值 99 double bestfitness = *(best_fit_index+1); 100 fitnessgbest = bestfitness; 101 102 //迭代寻优 103 for(int i=0;i<maxgen;i++) 104 { 105 for(int j=0;j<sizepop;j++) 106 { 107 //速度更新及粒子更新 108 for(int k=0;k<dim;k++) 109 { 110 // 速度更新 111 double rand1 = (double)rand()/RAND_MAX; //0到1之间的随机数 112 double rand2 = (double)rand()/RAND_MAX; 113 V[j][k] = V[j][k] + c1*rand1*(pbest[j][k]-pop[j][k]) + c2*rand2*(gbest[k]-pop[j][k]); 114 if(V[j][k] > Vmax) 115 V[j][k] = Vmax; 116 if(V[j][k] < Vmin) 117 V[j][k] = Vmin; 118 // 粒子更新 119 pop[j][k] = pop[j][k] + V[j][k]; 120 if(pop[j][k] > popmax) 121 pop[j][k] = popmax; 122 if(pop[j][k] < popmin) 123 pop[j][k] = popmin; 124 } 125 fitness[j] = func(pop[j]); //新粒子的适应度值 126 } 127 for(int j=0;j<sizepop;j++) 128 { 129 // 个体极值更新 130 if(fitness[j] > fitnesspbest[j]) 131 { 132 for(int k=0;k<dim;k++) 133 { 134 pbest[j][k] = pop[j][k]; 135 } 136 fitnesspbest[j] = fitness[j]; 137 } 138 // 群体极值更新 139 if(fitness[j] > fitnessgbest) 140 { 141 for(int k=0;k<dim;k++) 142 gbest[k] = pop[j][k]; 143 fitnessgbest = fitness[j]; 144 } 145 } 146 for(int k=0;k<dim;k++) 147 { 148 genbest[i][k] = gbest[k]; // 每一代最优值取值粒子位置记录 149 } 150 result[i] = fitnessgbest; // 每代的最优值记录到数组 151 } 152 } 153 154 // 主函数 155 int main(void) 156 { 157 clock_t start,finish; //程序开始和结束时间 158 start = clock(); //开始计时 159 srand((unsigned)time(NULL)); // 初始化随机数种子 160 PSO_func(); 161 double * best_arr; 162 best_arr = max(result,maxgen); 163 int best_gen_number = *best_arr; // 最优值所处的代数 164 double best = *(best_arr+1); //最优值 165 printf("迭代了%d次，在第%d次取到最优值，最优值为:%lf.\n",maxgen,best_gen_number+1,best); 166 printf("取到最优值的位置为(%lf,%lf).\n",genbest[best_gen_number][0],genbest[best_gen_number][1]); 167 finish = clock(); //结束时间 168 double duration = (double)(finish - start)/CLOCKS_PER_SEC; // 程序运行时间 169 printf("程序运行耗时:%lf\n",duration); 170 return 0; 171 }
我运行C采用的是Ubuntu16 下的gcc编译器，运行结果截图如下：
多次运行结果差不多，基本每次都可以很接近最优解。而且发现C语言运行时间要远快于MATLAB实现（我记得MATLAB要用好几秒，这里就不贴MATLAB代码进行运行时间对比了），只需要耗时0.004秒左右。这里只讨论了基本的粒子群算法，后面一篇我还会对于粒子群的参数w进行详细的讨论，讨论不同的w参数的取法对于粒子群寻优能力的影响。
