粒子群_Matlab程序代碼

1、粒子群優(yōu)化算法(panicle swarm optimization，PSO)是kennedy和Eberhart在研究鳥類和魚類的群體行為基礎(chǔ)上于1995年提出的一種群智能算法，其思想米源予人工生命和演化計算理論，模仿鳥群飛行覓食行為，通過鳥集體協(xié)作使群體達(dá)到最優(yōu)。1.粒子群算法的原理PSO中，每個優(yōu)化問題的解看作搜索空間中的一只鳥(即粒子)，所有的粒子都有一個被優(yōu)化的函數(shù)決定的適應(yīng)值，并且有一個速度決定它們飛翔的方向和速率，粒子們追隨當(dāng)前的最優(yōu)粒子在解空間中搜索。算法首先初始化一群隨機(jī)粒子，然后通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個“極值”即個體極值和全局極值來更新自己的速度與

2、位置。在D維目標(biāo)搜索空間中，由種群數(shù)為m的粒子組成粒子群，其中第f個粒子在第d維的位置為Xid，其飛行速度為Vid，該粒子當(dāng)前搜索到的最優(yōu)位置為Pid(goodvalue)和整個粒子群當(dāng)前的最優(yōu)位置Pgd(bestvalue)。每維的速度與位置更新公式如下W為慣性權(quán)重，C1和C2為學(xué)習(xí)因子，rand（）0,1范圍內(nèi)變化的隨機(jī)數(shù)。2.參數(shù)介紹與設(shè)置（1）ww是保持粒子運動慣性的參數(shù)，能使種群擴(kuò)展搜索空間，獲得較好的求解效果。較大的w有利于群體在更大的范圍內(nèi)進(jìn)行搜索。而較小的w能夠保證群體收斂到最優(yōu)位置，所以w的選擇及在迭代中的變化對搜索能力和跳出局優(yōu)能力具有重要影響。一般將w設(shè)定為0.8左右。（

3、1）加速因子c1和c2c1和c2用于調(diào)整粒子自身經(jīng)驗和社會經(jīng)驗在其運動中的作用，表示將每個粒子拉向pbest和gbest位置的隨機(jī)加速項的權(quán)重，低的值允許粒子在被拉回前可以在目標(biāo)區(qū)域外徘徊, 而高的值則導(dǎo)致粒子突然沖向或越過目標(biāo)區(qū)域。如果c1=0, 則粒子沒有認(rèn)知能力。在粒子的相互作用下,能到達(dá)新的搜索空間, 但容易陷入局部極值點; 如果c2=0, 粒子間沒有社會信息共享, 其算法變成一個多起點的隨機(jī)搜索; 如果c1=c2=0, 粒子將一直以當(dāng)前的速度飛行, 直到到達(dá)邊界。恰當(dāng)?shù)剡x擇c1與c2能較好的得到最優(yōu)解，一般都設(shè)定為2。（2）最大速度vmaxvmax決定當(dāng)前位置與最好位置之間的區(qū)域的分

4、辨率（或精度）。如果vmax太高, 微?？赡軙w躍最優(yōu)解; 如果vmax太小, 則微粒易陷入局部最優(yōu)。引入慣性權(quán)重w可消除對vmax的需要, 因為兩者的作用都是維護(hù)全局和局部搜索能力的平衡。當(dāng)vmax增加時, 可通過減小w來平衡搜索, 而w的減小可使得所需的迭代次數(shù)變小。但常常將設(shè)定一個最大的vmax。較好地選擇和調(diào)整參數(shù)能夠加大搜索能力，避免PSO的早熟現(xiàn)象。3.算法步驟及設(shè)計PSO的算法框架如下所述，圖1給出了PS0的算法流程。(1)初始化所有的個體(粒子)，初始化它們的速度和位置，并且將個體的歷史最優(yōu)goodvalue設(shè)為當(dāng)前位置，而群體中最優(yōu)的個體歷史位置作為當(dāng)前的全局最優(yōu)bestva

5、lue。(2)在當(dāng)代的進(jìn)化中，計算各個粒子的適應(yīng)度函數(shù)值。(3)如果該粒子當(dāng)前的適應(yīng)度函數(shù)值比其歷史最優(yōu)值要好，那么歷史最優(yōu)將會被當(dāng)前位置所替代。（4）如果該粒子的歷史最優(yōu)比全局最優(yōu)要好，那么全局最優(yōu)將會被該粒子的歷史最優(yōu)所替代。(5)對每個粒子按照公式(1)和公式(2)對速度和位置進(jìn)行更新。(6)進(jìn)化代數(shù)增加1，如果還沒有到達(dá)結(jié)束條件，轉(zhuǎn)到(2)步，否則輸出全局最優(yōu)bestvalue并結(jié)束。開始初始化粒子群計算每個粒子新位置的適應(yīng)值根據(jù)粒子適應(yīng)值更新個體最優(yōu)值和全局最優(yōu)值根據(jù)公式更新粒子速度和位置N是否滿足終止條件Y輸出最優(yōu)解結(jié)束圖一 PSO算法流程4.程序設(shè)計與仿真（1）程序設(shè)計說明采用M

6、ATLAB軟件編程，程序粒子種群個數(shù)設(shè)定為20個，搜索范圍為全局版本，終止條件是粒子都到達(dá)最優(yōu)點。（2）函數(shù)，其圖形如圖2所示，其在x=0.5，y=3.處有最小值f=-7.25.運行程序大約迭代220步，就可以讓所有的粒子到達(dá)最優(yōu)點。圖3中坐標(biāo)圖從左到右從上到下分別為初始化種群的位置、k=50的粒子位置、k=100的粒子位置及迭代結(jié)束粒子的位置.從圖中可見能夠搜索到最優(yōu)點。圖2 f1三維圖形 圖3 仿真粒子位置圖（2）Rastrigin函數(shù)，當(dāng)n=2時，為二維圖形如圖4所示，在x1=0，x2=0.處有最大值f=0；在程序中取相反數(shù)求其最小值，迭代k=670步左右將可以將所有的粒子推向一點最優(yōu)點

7、。取前十個粒子的位置如圖5，可知可以快速地得到滿意的最優(yōu)解，且能夠跳出局優(yōu)。圖4 f2函數(shù)n=2三維圖形圖5 粒子位置坐標(biāo)（3）rosenbrock函數(shù)：當(dāng)n=2時，函數(shù)為二維，當(dāng)x1=x2=1時有最小值f=0，其圖形如圖6，可知在最優(yōu)點附近，函數(shù)值下降很慢，利用將所有粒子推向一點終止條件，大約迭代k=1200次左右結(jié)束?？梢詮淖鴺?biāo)圖中（圖7）看到粒子的在迭代過程中位置。圖6 f3函數(shù)n=2三維圖形圖7 粒子坐標(biāo)位置當(dāng)n=5時，即是5維的rosenbrock函數(shù)，當(dāng)x1=x2=x3=x4=x5=1時有最小值f=0.采用粒子推向一點作為終止條件，運行程序迭代10萬次（人為設(shè)定）后，全局最優(yōu)值已到

8、了理想值（f=0）但此時粒子的位置還是相離較遠(yuǎn)，且運行時間較長（大約一分鐘），考慮程序目的就是找到最優(yōu)值，而不是讓所有粒子到達(dá)一點，因此采用迭代次數(shù)作為程序終止條件。依次增大設(shè)定的迭代次數(shù)觀察運行結(jié)果，如下表1表1 取不同迭代次數(shù)最優(yōu)點迭代次數(shù)k全局最優(yōu)點位置x1 x2 x3 x4 x5最優(yōu)解10000.0794 0.2819 0.5254 0.7267 0.85300.841530000.3718 0.6126 0.7832 0.8851 0.94070.215150000.5751 0.7583 0.8710 0.9328 0.96450.0816100000.7064 0.8413 0.

9、9165 0.9572 0.97840.0348500000.9937 0.9969 0.9984 0.9992 0.99960.00000.9999 1.0000 1.0000 1.0000 1.00000.0000由表中可知，隨著維數(shù)的增加，雖然迭代過程也是收斂的，但尋找最優(yōu)點很難，前已述及，PSO的搜索能力與參數(shù)設(shè)置有關(guān)，且全局版本相對局部版本有較差的跳出局優(yōu)的能力。為此若要增加搜索能力、提高迭代次數(shù)需要考慮改進(jìn)參數(shù)和更換局部版本搜索方式。（5）程序代碼（1）主函數(shù)：clearclchold offpsize=20;%粒子個數(shù)的設(shè)置pd=5;%粒子的維數(shù)lz=zeros(psize,pd

10、);for i=1:psize %粒子群隨機(jī)生成。psize行pd列，pd維。 lz(i,:)=rand(1,pd).*20-10;%【-1010之間】endlv=lz;lzfigure=lz;%畫初始圖用goodvalue=lz;%每個粒子自己歷史最好值初始化,psize行pd列vmax=20;%速度上限c1=2;c2=2;%學(xué)習(xí)因子w=0.729;%隨機(jī)因子和慣性因子。ri=0;bestvalue=zeros(1,pd);%全局最好值初始化，1行pd列for j=1:pd bestvalue(1,j)=goodvalue(1,j);endfnew=zeros(1,psize);for j=

11、1:psize fnew(j)=fpso(lz(1,:); endfold=fnew;flagstop=0;%終止標(biāo)志k=0;%迭代次數(shù)記錄f0=fpso(bestvalue);%適應(yīng)值初始化（默認(rèn)）while flagstop=0 %*.1.* for i=1:psize%適應(yīng)值比較，更新各自歷史最好值（位置） fnew(i)=fpso(lz(i,:);%記錄每次每個粒子的適應(yīng)值，便于以后設(shè)置終止條件 if fnew(i)fold(i) fold(i)=fnew(i);%fold記錄每個粒子的最好歷史值 for j=1:pd goodvalue(i,j)=lz(i,j); end end e

12、nd%*.2.* for i=1:psize%每個粒子歷史最好值相比較，更新全局最好值(及各維) f1=fold(i); if f1f0 f0=f1; for j=1:pd bestvalue(1,j)=goodvalue(i,j); end end end %*粒子趨一點終止條件*% % flagstop0=max(abs(fold);%比較當(dāng)次的所有粒子的適應(yīng)值， % flagstop1=min(abs(fold);%若它們之間差別較小，則可停止。 % if (flagstop0-flagstop1)1000 %迭代次數(shù)過大，強(qiáng)制終止，以防陷入死循環(huán) flagstop=1; end %*

13、%*速度和位置更新* for j=1:pd for i=1:psize lv(i,j)=w.*lv(i,j)+(c1*rand(1).*(goodvalue(i,j)-lz(i,j). +(c2*rand(1).*(bestvalue(1,j)-lz(i,j);%更新速度 if lv(i,j)vmax%速度最大值約束 lv(i,j)=vmax; end lz(i,j)=lz(i,j)+lv(i,j);%更新粒子位置 end end %*% % ri=ri+1; % rivalue(ri)=f0; % if ri=5000 % ri=0; % flagstop0=max(abs(rivalue)

14、; % flagstop1=min(abs(rivalue); % if (flagstop0-flagstop1)0.001 % flagstop=1; % end% end %*畫圖*% if k=0 %畫出初始粒子的位置 subplot(2,2,1); t=1:psize; plot(lzfigure(t,1),lzfigure(t,2),*); axis(-10,10,-10,10); end for j=2:3%分別畫出迭代300次和600次粒子的位置 if k=300*j subplot(2,2,j); t=1:psize; % gtext(k=600); plot(lz(t,1)

15、,lz(t,2),*); axis(0,10,0,10); end end if flagstop=1 %畫出終止時粒子的位置 subplot(2,2,4); t=1:psize; plot(lz(t,1),lz(t,2),*); axis(0,10,0,10); end%*% k=k+1;endfprintf(總迭代次數(shù)：k=%ld。n,k-1);fprintf(每個粒子的歷史最好值：n);for i=1:psize fprintf(粒子%d： %6.4f,%6.4fn,i,goodvalue(i,1),goodvalue(i,2);endfprintf(全局最好值：%6.4f,%6.4f.

16、n,bestvalue(1,1),bestvalue(1,2);fprintf(最優(yōu)解(最大值)：%6.4f.n,fpso(bestvalue(1,:);for j=1:20 fprintf(粒子%d位置：,j) ; fprintf(x=%6.4f，y=%6.4f,lz(j,1),lz(j,2); fprintf(n) ;end（2）調(diào)用函數(shù)：function fvalue=fpso(x)%fvalue=(x(1)+2)*(x(1)-3)+(x(2)-4)*(x(2)-2);%f1pd=5;%粒子的維數(shù)%*f3*fvalue=0; for r=1:(pd-1) f0=100*(x(r)-x(r+1)2)2+(x(r+1)-1)2; fvalue=fvalue+f0; end %* %fvalue=(x(1)-2)*(x(1)-2.5)+(x(2)-4)