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參數的取法對於粒子群尋優能力的影響。
C語言實現粒子群算法(PSO)一