C語言實現標準PSO演算法

簡介以及基本PSO的實現：http://blog.csdn.net/fovwin/article/category/1256709

相對於基本PSO，標準PSO加入了慣性權重係數W，W代表者粒子群對全域空間的搜尋能力和局部收斂速度的權衡。

也就是說，若慣性權重W越大，速度更新式子的第一項所佔的比重，即速度慣性項，比較大。粒子就會比較不受約束，可以衝來衝去。不受世俗的約束，也就不容易變“俗”（陷入局部極小點出不來）；若慣性權重W越小，速度更新式子的後兩項所佔的比重，即個體曆史極值和全域最優解，比較大。粒子被自己的經驗和社會的世俗所約束了，聽“媽媽”的話，當然人生走的順暢，但是你一般在你20歲就知道你後輩子的路，即容易陷入局部極小點出不來——過早收斂早熟。

所以一般在一開始一般設定較大的W值，讓它充分搜尋全域空間（就像小時候不聽話，有時候老爸老媽說了這樣不行，我們也會自己試了再說，呵呵），到粒子迭代的後半階段，設定W較小，可以提高收斂速度（後半輩子大部分按經驗活著）。

本文的W按照如下設定：

/*慣性權重係數*/#define W_START 1.4#define W_END0.4#define W_FUN(W_START-(W_START-W_END)*pow((double)cur_n/ITE_N,2))

當然，我們可以有的別的設定公式，但是大體呈現下降的趨勢，也可以用模糊控制規則（下下步實現目標）來設定～～～

至於作用效果麼，比基本的能更快早到最優解～～～不過木有pp，1.X版再說～～～

相對與Verson 0.0版減少好幾個函數，第一版初級階段，主要為了學習原理，Version 1.0版代碼更加簡潔，將pso.c和pso.h獨立出來，好以後擴充，通過define來設定你的待最佳化函數（下一步嘗試多種詭異的函數，有的還蠻漂亮的），但是還要完善，但是加入了擴充適應度函數的架構～～～將初始化第一次計算P，X和之後的合并為一體～～～

main.c

#include "pso.h"#include "stdio.h"int main(int argc, const char *argv[]){cur_n=0;RandInitofSwarm();             //初始化粒子群while(cur_n++ != ITE_N){UpdatePandGbest();    //更新個體曆史最優解P和全域最優解GBestUpdateofVandX();        //速度和位置更新，即飛翔}getchar();return 0;}

pso.h

#ifndef _PSO_H_#define _PSO_H_/*各種適應度函數選擇，要用哪個，就設定為1,但只能有一個為1*/#define TEST 0  #define GOLDSTEIN_PRICE0#define SCHAFFER 1#define HANSEN0#define NEEDLE0#define Dim 2//粒子維度#define PNum 20//種群規模#define ITE_N  50  //最大迭代次數int cur_n;//當前迭代次數/*慣性權重函數*/#define W_START 1.4#define W_END0.4#define W_FUN(W_START-(W_START-W_END)*pow((double)cur_n/ITE_N,2))/*個體和種群結構體*/struct PARTICLE{double X[Dim];double P[Dim];double V[Dim];double Fitness;} particle;struct SWARM{struct PARTICLE Particle[PNum];int GBestIndex;double GBest[Dim];double C1;double C2;double Xup[Dim];double Xdown[Dim];double Vmax[Dim];} swarm;/*是的，只要三個就好了，更好理解點感覺*/void RandInitofSwarm(void);void UpdateofVandX(void);void UpdatePandGbest(void);#endif

pso.c

#include "stdio.h"#include "stdlib.h"#include "time.h"#include "math.h"#include "pso.h"//初始化種群void RandInitofSwarm(void){int i, j;        //學習因子C1,C2swarm.C1 = 2.0;swarm.C2 = 2.0;for(j = 0; j < Dim; j++){swarm.Xdown[j] = -10;    //搜尋空間範圍swarm.Xup[j] = 10;swarm.Vmax[j] = 0.1;       //粒子飛翔速度最大值}srand((unsigned)time(NULL));for(i = 0; i < PNum; i++){for(j = 0; j < Dim; j++){swarm.Particle[i].X[j] = rand() / (double)RAND_MAX * (swarm.Xup[j] - swarm.Xdown[j]) + swarm.Xdown[j];//-Xdown~Xupswarm.Particle[i].V[j] = rand() / (double)RAND_MAX * swarm.Vmax[j] * 2 - swarm.Vmax[j];//-Vmax~Vmax}}}/*計算適應度函數的適應度，待進一步完善*/static double ComputAFitness(double X[]){/*OKmin-f(0,0)=3*/#if TESTreturn X[0]*X[0]+X[1]*X[1]+3;#endif/*min-f(0,-1)=3.x,y-(-2,2)*/#if GOLDSTEIN_PRICEreturn (1+pow(X[0]+X[1]+1,2)*(19-14*X[0]+3*pow(X[0],2)-14*X[1]+6*X[0]*X[1]+ 3*pow(X[1],2))*(30+pow((2*X[0]-3*X[1]),2)*\(18-32*X[0]+12*pow(X[0],2)+48*X[1]-36*X[0]*X[1]+ 27*pow(X[1],2))));#endif/*min-f(0,0)=-1.x,y-(-4,4)*/#if SCHAFFER//函數:Schaffer's F6return -1*(0.5-(sin(sqrt(X[0]*X[0]+X[1]*X[1]))*\sin(sqrt(X[0]*X[0]+X[1]*X[1]))-0.5)/pow(( 1+0.001*(X[0]*X[0]+X[1]*X[1])),2));#endif/*此函數局部極值點有760個。min-f(x,y)=-176.541793.x,y-(-10,10)(-7.589893,-7.708314),(-7.589893,-1.425128)(-7.5898893,4.858057),(-1.306708,-7.708314)(-1.306707,-1.425128),(-1.306708,4.858057)(4.976478,-7.708314),(4.976478,-1.425128)(4.976478,4.858057)*/#if HANSENint i;double temp1=0,temp2=0;double hansenf=0;for (i=1;i <= 5;i++){temp1+=i*cos((i-1)*X[0]+i);temp2+=i*cos((i-1)*X[1]+i);}hansenf=-1*temp1*temp2;return hansenf;#endif/*min-f(0,0)=-3600,x,y-(-5.12,5.12)該問題的全域最優解被最差解包圍,局部極值點：(-5.12,5.12),(-5.12,-5.12)(5.12,-5.12),(5.12,5.12)f=-2748.78*/#if NEEDLEreturn -1*pow((3/(0.05+pow(X[0]-1,2)+pow(X[1]-1,2))),2);#endif}//update  V and Xvoid UpdateofVandX(void){int i, j;srand((unsigned)time(NULL));for(i = 0; i < PNum; i++){for(j = 0; j < Dim; j++)swarm.Particle[i].V[j] = W_FUN * swarm.Particle[i].V[j] +rand() / (double)RAND_MAX * swarm.C1 * (swarm.Particle[i].P[j] - swarm.Particle[i].X[j]) +rand() / (double)RAND_MAX * swarm.C2 * (swarm.GBest[j] - swarm.Particle[i].X[j]);for(j = 0; j < Dim; j++){if(swarm.Particle[i].V[j] > swarm.Vmax[j])swarm.Particle[i].V[j] = swarm.Vmax[j];if(swarm.Particle[i].V[j] < -swarm.Vmax[j])swarm.Particle[i].V[j] = -swarm.Vmax[j];}for(j = 0; j < Dim; j++){swarm.Particle[i].X[j] += swarm.Particle[i].V[j];if(swarm.Particle[i].X[j] > swarm.Xup[j])swarm.Particle[i].X[j] = swarm.Xup[j];if(swarm.Particle[i].X[j] < swarm.Xdown[j])swarm.Particle[i].X[j] = swarm.Xdown[j];}}}/*更新個體極值P和全域極值GBest*/void UpdatePandGbest(void){int i, j;//update of P if the X is bigger than current Pfor (i = 0; i < PNum; i++){if (swarm.Particle[i].Fitness < ComputAFitness(swarm.Particle[i].P)){for(j = 0; j < Dim; j++){swarm.Particle[i].P[j] = swarm.Particle[i].X[j];}}}//update of GBestfor(i = 0; i < PNum; i++)if(ComputAFitness(swarm.Particle[i].P) < ComputAFitness(swarm.Particle[swarm.GBestIndex].P))swarm.GBestIndex = i;for(j = 0; j < Dim; j++){swarm.GBest[j] = swarm.Particle[swarm.GBestIndex].P[j];}printf("The %dth iteraction.\n",cur_n);printf("GBestIndex:%d \n",swarm.GBestIndex );printf("GBest:" );for(j=0;j<Dim;j++){printf("%.4f ,",swarm.GBest[j]);}printf("\n" );printf("Fitness of GBest: %f \n\n",ComputAFitness(swarm.Particle[swarm.GBestIndex].P));}

Version 1.0 標準PSO