粒子群演算法(8)—混合粒子群演算法的實現

混合粒子群演算法將全域粒子群演算法與局部粒子群演算法結合，其速度更新採用公式

其中G（k+1）是全域版本的速度更新公式，而L（k+1）是局部版本的速度更新公式，混合粒子群演算法採用H（k+1）的公式。

位置更新公式

因為是局部版本與全域版本相結合，所以，粒子群的初始化函數應該與局部版本的相同，這裡就不列出了，參看粒子群演算法(7）中的LocalInitSwarm函數。

關鍵還是混合粒子群演算法的單步更新函數，函數名為HybridStepPso

代碼如下：

function [ParSwarm,OptSwarm]=HybridStepPso(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)<br />% 功能描述：混合粒子群演算法。將全域版本與局部版本相混合。基本的粒子群演算法的單步更新位置,速度的演算法<br />%<br />%[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)<br />%<br />%演算法思想：全域版本的速度更新公式：vq(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pq-w)<br />%pg為個體曆史最優，pq為全域最優<br />%局部版本的速度更新公式：vl(k+1)=w*v(k)+c1*r1*(pg-w)+c2*r2*(pl-w) pl為鄰域最優<br />%現在速度更新公式vh=n*vq+(1-n)*vl;n屬於0到1的一個數<br />% 輸入參數：ParSwarm:粒子群矩陣，包含粒子的位置，速度與當前的目標函數值<br />%輸入參數：OptSwarm：包含粒子群個體最優解與全域最優解的矩陣<br />%輸入參數：ParticleScope:一個粒子在運算中各維的範圍；<br />% 輸入參數：AdaptFunc：適應度函數<br />%輸入參數：LoopCount：迭代的總次數<br />%輸入參數：CurCount：當前迭代的次數<br />%<br />% 傳回值：含意同輸入的同名參數<br />%<br />%用法：[ParSwarm,OptSwarm]=LocalStepPsoByCircle(ParSwarm,OptSwarm,AdaptFunc,ParticleScope,MaxW,MinW,LoopCount,CurCount)<br />%<br />% 異常：首先保證該檔案在Matlab的搜尋路徑中，然後查看相關的提示資訊。<br />%<br />%編製人：XXX<br />%編製時間：2010.5.6<br />%參考文獻：XXX<br />%參考文獻：XXX<br />%<br />%修改記錄<br />%----------------------------------------------------------------<br />%2010.5.6<br />%修改人：XXX<br />% 添加2*unifrnd(0,1).*SubTract1(row,:)中的unifrnd(0,1)隨機數，使效能大為提高<br />%修改混合因子，從小到大，開始從CurCount/LoopCount開始<br />%修改C1=2.05,C2=2.05<br />%修改速度的範圍時區間每維範圍的0.5(一半)<br />%參照基於MATLAB的粒子群最佳化演算法程式設計<br />%<br />% 總體評價：使用這個版本的調節係數，效果比較好<br />%</p><p>%容錯控制<br />if nargin~=8<br /> error('輸入的參數個數錯誤。')<br />end<br />if nargout~=2<br /> error('輸出的個數太少，不能保證迴圈迭代。')<br />end</p><p>%開始單步更新的操作</p><p>%*********************************************<br />%***** 更改下面的代碼，可以更改慣性因子的變化*****<br />%---------------------------------------------------------------------<br />% 線形遞減策略<br />w=MaxW-CurCount*((MaxW-MinW)/LoopCount);<br />%---------------------------------------------------------------------<br />%w 固定不變策略<br />%w=0.7;<br />%---------------------------------------------------------------------<br />% 參考文獻：陳貴敏，賈建援，韓琪，粒子群最佳化演算法的慣性權值遞減策略研究，西安交通大學學報，2006，1<br />%w 非線形遞減，以凹函數遞減<br />%w=(MaxW-MinW)*(CurCount/LoopCount)^2+(MinW-MaxW)*(2*CurCount/LoopCount)+MaxW;<br />%---------------------------------------------------------------------<br />%w 非線形遞減，以凹函數遞減<br />%w=MinW*(MaxW/MinW)^(1/(1+10*CurCount/LoopCount));<br />%*****更改上面的代碼，可以更改慣性因子的變化*****<br />%*********************************************<br />%更改下面代碼可以改變混合因子的取值<br />%-----------------------------------------------<br />Hybrid=CurCount/LoopCount;<br />%-----------------------------------------------</p><p>% 得到粒子群群體大小以及一個粒子維數的資訊<br />[ParRow,ParCol]=size(ParSwarm);</p><p>%得到粒子的維數<br />ParCol=(ParCol-1)/2;<br />GlobleSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身曆史最優解位置減去粒子當前位置<br />LocalSubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%粒子自身曆史最優解位置減去粒子當前位置<br />LocalSubTract2=OptSwarm(ParRow+1:end-1,:)-ParSwarm(:,1:ParCol);%粒子鄰域最優解位置減去粒子當前位置<br />%*********************************************<br />%***** 更改下面的代碼，可以更改c1,c2的變化*****<br />c1=2.05;<br />c2=2.05;<br />%---------------------------------------------------------------------<br />%con=1;<br />%c1=4-exp(-con*abs(mean(ParSwarm(:,2*ParCol+1))-AdaptFunc(OptSwarm(ParRow+1,:))));<br />%c2=4-c1;<br />%----------------------------------------------------------------------<br />%***** 更改上面的代碼，可以更改c1,c2的變化*****<br />%*********************************************<br />for row=1:ParRow<br /> GlobleSubTract2=OptSwarm(ParRow*2+1,:)-ParSwarm(row,1:ParCol);%全域最優的位置減去每個粒子當前的位置<br /> LocalTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*LocalSubTract1(row,:)+c2*unifrnd(0,1).*LocalSubTract2(row,:);<br /> GlobleTempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*GlobleSubTract1(row,:)+c2*unifrnd(0,1).*GlobleSubTract2;<br /> TempV=Hybrid.*GlobleTempV+(1-Hybrid).*LocalTempV;<br /> %限制速度的代碼<br /> for h=1:ParCol<br /> if TempV(:,h)>ParticleScope(h,2)/2.0<br /> TempV(:,h)=ParticleScope(h,2)/2.0;<br /> end<br /> if TempV(:,h)<-ParticleScope(h,2)/2.0<br /> TempV(:,h)=(-ParticleScope(h,2)+1e-10)/2.0;<br /> %加1e-10防止適應度函數被零除<br /> end<br /> end </p><p> % 更新速度<br /> ParSwarm(row,ParCol+1:2*ParCol)=TempV;</p><p> %*********************************************<br /> %***** 更改下面的代碼，可以更改約束因子的變化*****<br /> %---------------------------------------------------------------------<br /> %a=1;<br /> %---------------------------------------------------------------------<br /> a=0.729;<br /> %***** 更改上面的代碼，可以更改約束因子的變化*****<br /> %*********************************************</p><p> % 限制位置的範圍<br /> TempPos=ParSwarm(row,1:ParCol)+a*TempV;<br /> for h=1:ParCol<br /> if TempPos(:,h)>ParticleScope(h,2)<br /> TempPos(:,h)=ParticleScope(h,2);<br /> end<br /> if TempPos(:,h)<=ParticleScope(h,1)<br /> TempPos(:,h)=ParticleScope(h,1)+1e-10;<br /> end<br /> end</p><p> %更新位置<br /> ParSwarm(row,1:ParCol)=TempPos;</p><p> % 計算每個粒子的新的適應度值<br /> ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));<br /> if ParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))<br /> OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);<br /> end<br />end<br />%for迴圈結束<br />%確定鄰域的範圍<br />linyurange=fix(ParRow/2);<br />%確定當前迭代的鄰域範圍<br />jiange=ceil(LoopCount/linyurange);<br />linyu=ceil(CurCount/jiange);<br /> for row=1:ParRow<br /> if row-linyu>0&&row+linyu<=ParRow<br /> tempM =[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:row+linyu,:)];</p><p> [maxValue,linyurow]=max(tempM(:,2*ParCol+1));<br /> if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))<br /> OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);<br /> end<br /> else<br /> if row-linyu<=0<br /> %該行上面的部分突出了邊界，下面絕對不會突破邊界<br /> if row==1<br /> tempM=[ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];<br /> [maxValue,linyurow]=max(tempM(:,2*ParCol+1));<br /> if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))<br /> OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);<br /> end<br /> else<br /> tempM=[ParSwarm(1:row-1,:);ParSwarm(ParRow+row-linyu:end,:);ParSwarm(row+1:row+linyu,:)];<br /> [maxValue,linyurow]=max(tempM(:,2*ParCol+1));<br /> if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))<br /> OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);<br /> end<br /> end<br /> else<br /> %該行下面的部分突出了邊界，上面絕對不會突破邊界<br /> if row==ParRow<br /> tempM=[ParSwarm(ParRow-linyu:row-1,:);ParSwarm(1:linyu,:)];<br /> [maxValue,linyurow]=max(tempM(:,2*ParCol+1));<br /> if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))<br /> OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);<br /> end<br /> else<br /> tempM=[ParSwarm(row-linyu:row-1,:);ParSwarm(row+1:end,:);ParSwarm(1:linyu-(ParRow-row),:)];<br /> [maxValue,linyurow]=max(tempM(:,2*ParCol+1));<br /> if maxValue>AdaptFunc(OptSwarm(ParRow+row,:))<br /> OptSwarm(ParRow+row,:)=tempM(linyurow,1:ParCol);<br /> end<br /> end<br /> end<br /> end<br /> end%for<br /> %尋找適應度函數值最大的解在矩陣中的位置(行數)，進行全域最優的改變<br />[maxValue,row]=max(ParSwarm(:,2*ParCol+1));<br />if AdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow*2+1,:))<br /> OptSwarm(ParRow*2+1,:)=ParSwarm(row,1:ParCol);<br />end</p><p>

注意代碼的91行到96行，這幾行就是混合粒子群速度更新公式，其他部分基本與前面的實現一樣。

最後還是一個把這兩個函數組裝在一起的函數，同樣採用LocalPsoProcessByCircle函數，詳細見粒子群演算法(7）的內容，最後還是給出一個應用執行個體。

Scope=[-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10;-10 10];//粒子每維的限制範圍<br />qun=20;//粒子群的規模<br />lizi=10;//每個粒子的維數<br />[Result,OnLine,OffLine,MinMaxMeanAdapt,BestofStep]=LocalPsoProcessByCircle(qun,lizi,Scope,@localinitswarm,@Hybridsteppso,@Rastrigin,0,0,1000,0);

注意：在這個LocalPsoProcessByCircle函數中，使用HybridStepPso作為單步更新的函數，其餘基本與局部粒子群演算法相同。

經過本人的實際測試，運行條件相同，最好的是局部版本的PSO，混合的PSO並不像有些文獻上說的那麼好，也許是我實現的不對，如果有那個大俠實現的效果更好，可以給我聯絡，我們可以共用代碼。

同時也希望那些磚家、叫獸們共用你們的效果非常好的代碼。

本人已經實現了一個PSO的工具箱，不過效果不好，本人水平低劣，又需要的可以聯絡我。

不知道CSDN能不能做連結下載，如果可以，請告訴我，我做個連結，大家可以隨便下載，共同交流。