粒子群算法
粒子群算法，也稱粒子群優化算法或鳥群覓食算法（Particle Swarm Optimization），縮寫為 PSO， 是近年來由J. Kennedy和R. C. Eberhart等開發的一種新的進化算法(Evolutionary Algorithm - EA)。PSO 算法屬於進化算法的一種，和模擬退火算法相似，它也是從隨機解出發，通過迭代尋找最優解，它也是通過適應度來評價解的品質，但它比遺傳算法規則更為簡單，它沒有遺傳算法的“交叉”(Crossover) 和“變異”(Mutation) 操作，它通過追隨當前搜索到的最優值來尋找全局最優。這種算法以其實現容易、精度高、收斂快等優點引起了學術界的重視，並且在解決實際問題中展示了其優越性。粒子群算法是一種並行算法。
1算法介紹
粒子群優化算法的基本思想是通過群體中個體之間的協作和信息共享來尋找最優解。
PSO的優勢在於簡單容易實現並且沒有許多參數的調節。目前已被廣泛應用於函數優化、神經網絡訓練、模糊系統控制以及其他遺傳算法的應用領域。
1.1問題提出
設想這樣一個場景：一群鳥在隨機的搜索食物。 在這個區域里只有一塊食物，所有的鳥都不知道食物在哪。但是它們知道自己當前的位置距離食物還有多遠。 那么找到食物的最優策略是什么？最簡單有效的就是搜尋目前離食物最近的鳥的周圍區域。
1.2問題抽象
鳥被抽象為沒有質量和體積的微粒(點)，並延伸到N維空間，粒子在N維空間的位置表示為矢量Xi＝(x1，x2，…，xN)，飛行速度表示為矢量Vi＝(v1，v2，…，vN)。每個粒子都有一個由目標函數決定的適應值(fitness value)，並且知道自己到目前為止發現的最好位置(pbest)和現在的位置Xi。這個可以看作是粒子自己的飛行經驗。除此之外，每個粒子還知道到目前為止整個群體中所有粒子發現的最好位置(gbest)(gbest是pbest中的最好值)。這個可以看作是粒子同伴的經驗。粒子就是通過自己的經驗和同伴中最好的經驗來決定下一步的運動。
1.3算法描述
PSO初始化為一群隨機粒子(隨機解)。然后通過迭代找到最優解。在每一次的迭代中，粒子通過跟蹤兩個“極值”(pbest,gbest)來更新自己。
在找到這兩個最優值后，粒子通過下面的公式來更新自己的速度和位置。
i＝1，2，…，M，M是該群體中粒子的總數；Vi 是粒子的速度； pbest和gbest如前定義； rand()是介於（0、1）之間的隨機數； xi是粒子的當前位置。c1和c2是學習因子，通常取c1＝ c2＝2 在每一維，粒子都有一個最大限制速度Vmax，如果某一維的速度超過設定的Vmax ，那么這一維的速度就被限定為Vmax 。（ Vmax >0）以上面兩個公式為基礎，形成了后來PSO 的標准形式。
1.4PSO算法流程
Step1：初始化一群微粒(群體規模為m)，包括隨機位置和速度；
Step2：評價每個微粒的適應度；
Step3：對每個微粒，將其適應值與其經過的最好位置 pbest作比較，如果較好，則將其作為當前的最好位置pbest;
Step4：對每個微粒，將其適應值與其經過的最好位置 gbest作比較，如果較好，則將其作為當前的最好位置gbest;
Step5：根據(2)、(3)式調整微粒速度和位置；
Step6：未達到結束條件則轉Step2。
迭代終止條件根據具體問題一般選為最大迭代次數Gk或和微粒群迄今為止搜索到的最優位置滿足預定最小適應閾值。
1.5參數分析
方程中pbest和gbest分別表示微粒群的局部和 全局最優位置，當c1＝0時，則粒子沒有了認知能力， 變為只有社會的模型(social-only):
被稱為全局PSO算法.粒子有擴展搜索空間的能力，具有 較快的收斂速度，但由於缺少局部搜索，對於復雜問題 比標准PSO 更易陷入局部最優。
當c2＝0時，則粒子之間沒有社會信息，模型變為 只有認知(cognition-only)模型：
被稱為局部PSO算法。由於個體之間沒有信息的 交流，整個群體相當於多個粒子進行盲目的隨機 搜索，收斂速度慢，因而得到最優解的可能性小。
群體規模m 一般取20～40，對較難或特定類別的問題 可以取到100～200。
最大速度Vmax決定當前位置與最好位置之間的區域的 分辨率(或精度)。如果太快，則粒子有可能越過極小 點;如果太慢，則粒子不能在局部極小點之外進行足 夠的探索，會陷入到局部極值區域內。這種限制可以 達到防止計算溢出、決定問題空間搜索的粒度的目的。
學習因子c1和c2代表將每個粒子推向pbest 和gbest位置的統計加速項的權值。較低的值允許粒子 在被拉回之前可以在目標區域外徘徊，較高的值導致粒 子突然地沖向或越過目標區域。
2算法模型
目標函數：
約束條件：
3算法實現流程
4算法代碼：
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define pop_size 20//種群規模
#define dim 2//粒子維度
#define PI 3.1415926 //圓周率
#define iteration_max 300//迭代次數
#define c2 2 //學習因子
#define c1 2 //學習因子
#define vmax 0.5//速度最大值
#define vmin -0.5//速度最小值
#define popmax 2//個體位置的最大值
#define popmin -2//個體位置的最小值
double pop[pop_size][dim];//種群位置的數組
double speed[pop_size][dim];//種群速度
double fitness[pop_size];//種群適應度數組
double gbest[dim];//群體極值的位置
double pbest[pop_size][dim];//個體極值的位置
double fitnesspbest[pop_size];//個體極值適應度值
double fitnessgbest;//群體極值適應度值
double iteration_best[iteration_max][dim]; //每一代最優值取值粒子
double result[iteration_max]; //定義存放每次迭代種群最優值的數組
double fitness_func(double *arr){
double x=*arr;
double y=*(arr+1);
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;
//printf("%lf\n",fitness);
return fitness;
}
void init()
{
for(int i=0;i<pop_size;i++)
{
for(int j=0;j<dim;j++)
{
pop[i][j] = (((double)rand())/RAND_MAX-0.5)*4;//-2到2之間的隨機數
speed[i][j] = ((double)rand())/RAND_MAX-0.5; //-0.5到0.5之間
}
fitness[i]=fitness_func(pop[i]);
//printf("fit:%lf\n",fitness[i]);
}
}
double *max(double *fitness,int size)
{
int index=0;
double max=*fitness;
static double best_fit_index[2];
for(int i=0;i<size;i++)
{
if(*(fitness+i)>max)
{
max=*(fitness+i);
index=i;
}
}
best_fit_index[0]=index;
best_fit_index[1]=max;
return best_fit_index;
}
void pso(){
init();
double * best_fit_index; // 用於存放群體極值和其位置(序號)
best_fit_index = max(fitness,pop_size); //求群體適應度最大值
int index=(int)(*best_fit_index);
printf("best:%lf\n",*(best_fit_index+1));
//群體極值的位置
for(int i=0;i<dim;i++)
{
gbest[i]=pop[index][i];
}
//個體極值的位置
for(int i=0;i<pop_size;i++)
{
for(int j=0;j<dim;j++)
{
pbest[i][j]=pop[i][j];
}
}
//個體極值適應度值
for(int i=0;i<pop_size;i++)
{
fitnesspbest[i]=fitness[i];
}
//群體極值適應度值
fitnessgbest=*(best_fit_index+1);
//迭代尋優
for(int i=0;i<iteration_max;i++)
{
for(int j=0;j<pop_size;j++)
{
for(int k=0;k<dim;k++)
{
//更新速度
double rand1=(double)rand()/RAND_MAX; //0到1之間的隨機數
double rand2=(double)rand()/RAND_MAX; //0到1之間的隨機數
speed[j][k]=speed[j][k]+c1*rand1*(pbest[j][k]-pop[j][k])+c2*rand2*(gbest[k]-pop[j][k]);
if(speed[j][k]>vmax)
{
speed[j][k]=vmax;
}
if(speed[j][k]<vmin)
{
speed[j][k]=vmin;
}
//更新粒子的位置
pop[j][k]=pop[j][k]+speed[j][k];
if(pop[j][k]>popmax)
{
pop[j][k]=popmax;
}
if(pop[j][k]<popmin)
{
pop[j][k]=popmin;
}
}
//計算新粒子的適應度值
fitness[j]=fitness_func(pop[j]);
}
for(int j=0;j<pop_size;j++)
{
//更新個體極值的位置
if(fitness[j]>fitnesspbest[j])
{
for(int k=0;k<dim;k++)
{
pbest[j][k]=pop[j][k];
}
fitnesspbest[j]=fitness[j];
}
//更新群體極值的位置
if(fitness[j]>fitnessgbest)
{
for(int k=0;k<dim;k++)
{
gbest[k]=pop[j][k];
}
fitnessgbest=fitness[j];
}
}
for(int k=0;k<dim;k++)
{
iteration_best[i][k] = gbest[k]; // 每一代最優值取值粒子位置記錄
}
result[i] = fitnessgbest; // 每代的最優值記錄到數組
printf("iteration %d:\n",i);
printf("Optimal value particle position: x=%lf, y=%lf\n",iteration_best[i][0],iteration_best[i][1]);
printf("Optimal fitness=%lf\n",result[i]);
}
printf("--------------------------------------------------------\n");
printf("final result:\n");
printf("Optimal value particle position: x=%lf, y=%lf\n",gbest[0],gbest[1]);
printf("Optimal value=%lf\n",fitnessgbest);
}
int main()
{
int i=0;
srand((unsigned int)time(NULL));
pso();
return 0;
}
5結果
可以看出函數的最優解為1.005236。
6結果分析
求解函數的圖形：
取第0、49、99、149、199、249、299代的20個粒子的位置做分析，可以看出粒子在迭代中進行收斂，最終在x=0，y=0處取得最優解。（圖中橫縱坐標分別為x、y）：
gnuplot繪圖代碼：
set xrange [-2:2] set yrange [-2:2] set size 1,1 set multiplot layout 2,2 set size 0.5,0.5 set origin 0,0.5 set title "0 iteration" plot '/home/unubtu/data0.dat' with points lc rgb 'blue' set origin 0.5,0.5 set title "49 iteration" plot '/home/unubtu/data1.dat' with points lc rgb 'blue' set origin 0,0 set title "99 iteration" plot '/home/unubtu/data2.dat' with points lc rgb 'blue' set origin 0.5,0 set title "149 iteration" plot '/home/unubtu/data3.dat' with points lc rgb 'blue' set xrange [-2:2] set yrange [-2:2] set size 1,1 set multiplot layout 2,2 set size 0.5,0.5 set origin 0,0.5 set title "199 iteration" plot '/home/unubtu/data4.dat' with points lc rgb 'blue' set origin 0.5,0.5 set title "249 iteration" plot '/home/unubtu/data5.dat' with points lc rgb 'blue' set origin 0,0 set title "299 iteration" plot '/home/unubtu/data6.dat' with points lc rgb 'blue' set isosamples 50,50 set pm3d set xlabel ‘x’ set ylabel ‘y’ set zlabel ‘f(x,y)’ splot [-2:2] [-2:2] sin(sqrt(x**2+y**2))/(sqrt(x**2+y**2))+exp((cos(2*pi*x)+cos(2*pi*y))/2)-2.71289