目錄

• 1.算法
• • 1.1.原理
  • 1.2.性能比較
  • 1.3.步驟
• 2.代碼
• • 2.1.源碼及注釋
  • 2.2.執行與效果

1.算法

1.1.原理

\qquad建議沒接觸過粒子群算法的朋友先看較為基礎的全局粒子群算法原理及介紹，以下博文鏈接有詳細的講解、代碼及其應用舉例：
【Simulink】粒子群算法（PSO）整定PID參數（附代碼和講解）
\qquad這里就介紹一下全局粒子群算法和混合粒子群算法的區別。全局粒子群算法（General PSO）將粒子速度矢量影響因子分為了三類：慣性因子、自身最優因子、社會因子。每個粒子的初速度定為0，即$v_0=0$，第$j$個粒子（$1\le j\le m$）的下一次迭代的速度$v^{(j)}$由以下三部分組成：
$$v^{(j)}=w?v_0+c_1?rand?(P^{(j)}?X^{(j)})+c_2?rand?(P_G?X^{(j)})$$
注：rand是(0,1)的隨機數，$v_0$代表上一次粒子的速度。
第一部分為自身慣性因子，因為下一次的迭代次數保留了上一次的速度信息；
第二個部分為自身最優因子，$P^{(j)}$為第$j$個因子在之前所有迭代次數中自適應度最高的位置，可以理解為歷史上自身離最優解最近的位置。
第三部分為社會因子，$P_G$為種群在之前所有迭代次數中自適應度最高的位置，可以理解為歷史上所有粒子離最優解最近的位置中的最近的一個。
\qquad而混合粒子群算法（Hybrid PSO）與全局粒子群算法唯一的區別就是將社會因子分解為了兩部分——局部社會因子和全局社會因子。即$$v^{(j)}=w?v_0+c_1?rand?(P^{(j)}?X^{(j)})+c_2[q?rand?(P_G?X^{(j)})+(1?q)?rand?(P_L^{(j)}?X^{(j)})]$$
一般0＜q＜1,即兩部分社會因子都為正。與上面不同的是$P_L^{(j)}$是第$j$個粒子的局部最優解的坐標。
\qquad全局社會因子使得粒子的速度矢量有向（目前的）全局最優解趨向的趨勢，而局部社會因子則使速度矢量包含向粒子的（目前的）局部最優解趨向的趨勢。全局最優解好理解，局部最優解就是該粒子一定范圍內粒子中的最優解，一般這個范圍是固定的，而包含的粒子數量（如果為1，則是粒子本身）不確定。
\qquad全局粒子群算法收斂速度快，但是對于多極值問題容易收斂到一個局部最優解；混合粒子群算法由于要計算局部最優解的原因，收斂速度慢，但可以很好地解決全局粒子群算法收斂到局部最優解的問題。

1.2.性能比較

\qquad對粒子群算法感興趣的朋友可以讀一下面這兩部分解釋。

1.為什么全局粒子群算法容易收斂到局部最優解呢？
\qquad我們不妨假設僅有3個粒子，目前有2個極小值點A和B，假設第一次迭代（位置隨機，初速度為0），一個粒子a恰好就位于其中一個極小值點A附近，另2個粒子b和c彼此較近，但離B較遠，函數值$f(A)<f(C)<f(B)$。初速度為0，慣性因子為0。又因為a就是當前全局最優點，自身最優因子也是本身（因為迭代次數為1），所以后兩項都是0，a不會動。粒子b和c的歷史最優因子就是它現在的位置，而社會因子使它向粒子a移動，并且每移動一次，它的自身最優因子都比上次的好，即不會對它的目前的速度造成任何影響。

\qquad最終兩個粒子都會收斂到極值點A。即使極值點B的目標函數值更小一些，但由于極值點A先被粒子群“找到”，因此就像一個極小值點黑洞一樣，把所有的粒子都“吸”了過去。所以簡單地來說，全局最優解就像一個吸引其他粒子的“黑洞”，如果在找到第二個“黑洞”之前所有粒子都被吸引到了第一個“黑洞”的周圍，那就只能找到第一個局部最優解了。

123

2.那為什么混合粒子群算法就可以避免這個問題呢？
\qquad我們不妨假設一個和上文一樣的場景。目標函數只有2個極小值點A和B，第一次迭代，粒子a就位于A的附近，粒子b和c不位于B的附近，但彼此離得較近，離A較遠。而且目標函數值都比a要大。初速度為0，因此三者的慣性因子為0。a的自身最優因子就是它本身，同時也是全局最優解，因此a不會動。而b、c的自身最優因子也是本身，然而全局社會最優因子是a，即二者都有向a偏移的趨勢。

\qquad現在考慮局部社會最優因子，b和c距離較小，會受到局部社會因子的影響；而二者都離a較遠，不會受到a的影響（同樣a也不會受到b和c的影響）。假設c的函數值比b小，因此b是c的局部最優解（同時也是它自身的）。在局部社會因子的作用下，b有向c移動的趨勢；c和a均不受影響（c和a的局部最優解都是本身）。假設局部最優因子的權重合適，b應該向極小值B移動，c向a點移動。因為b有慣性因子，很可能在第二次迭代的過程中，離極小值B的距離會比c要短，滿足$f(B)<f(C)<f(A)$。那么在接下來的幾步迭代中，三者就會順利收斂到極小值B附近，而不是更大的一個極小值A，進而找到全局最優解。

1234

1.3.步驟

\qquad步驟和全局粒子群算法基本一致，下面將不一樣的步驟標紅，看過全局粒子群算法那篇博文的讀者可以重點關注一下。
【step0】確定參數維度$N$、參數范圍（即每只鳥的初始位置），確定慣性系數$c_1,c_2,w$，確定種群規模m，迭代次數n以及局部因子作用半徑$R$。

每個粒子的初始位置是隨機的，設輸入向量$x=(x_1,x_2,...,x_N{)}^T$必須滿足參數范圍限制為：
$$?\begin{array}{l}{x}_{min}^{(1)}<x_1<x_{max}^{(1)}\\ x_{min}^{(2)}<x_2<x_{max}^{(2)}\\ ...\\ x_{min}^{(N)}<x_1<x_{max}^{(N)}\end{array}$$
記$X_{min}=(x_{min}^{(1)},x_{min}^{(2)},...x_{min}^{(N)}),X_{max}=(x_{max}^{(1)},x_{max}^{(2)},...x_{max}^{(N)})$
則輸入向量$x$應滿足$X_{min}<x<X_{max}$
【step1】計算每個粒子的速度
每個粒子的初速度定為0，即$v_0=0$，第$j$個粒子（$1\le j\le m$）的下一次迭代的速度$v^{(j)}$由三部分組成：
$$v^{(j)}=w?v_0+c_1?rand?(P^{(j)}?X^{(j)})+c_2[q?rand?(P_G?X^{(j)})+(1?q)?rand?(P_L^{(j)}?X^{(j)})]$$

注：rand是(0,1)的隨機數，$v_0$代表上一次粒子的速度。
第一部分為自身慣性因子，因為下一次的迭代次數保留了上一次的速度信息；
第二個部分為自身最優因子，$P^{(j)}$為第$j$個因子在之前所有迭代次數中自適應度最高的位置，可以理解為歷史上自身離食物最近的位置。
第三部分為社會因子，由全局社會因子和局部社會因子組成。

$P_G$為全局最優解，是種群在之前所有迭代次數中自適應度最高的位置，可以理解為歷史上所有粒子離食物最近的位置中的最近的一個。
一般情況下，取$w=1,c_1=c_2=2$，當種群規模較大時，可以讓$w$的值隨迭代次數減小以增加收斂速度。
$P_L^{(j)}$為局部最優解，每個粒子都有各自的局部最優解，可以理解為以該粒子為中心，$R$為半徑的超球體內包含的所有粒子中自適應度最高的粒子的位置。q是全局社會因子的占比，比例越高局部社會因子所占的權重越小。$R$稱為局部因子作用半徑，和粒子分布平均密度有關。

【step2】按照step1的速度公式計算下一次的速度，并計算下一次的粒子位置。對于第$j$個粒子，第$k+1$次迭代（第$k+1$代）的位置$X_{k+1}^{(j)}$與第$k$次迭代的位置$X_K^{(j)}$與速度$v_k^{(k+1)}$關系為：
$$X_{k+1}^{(j)}=X_k^{(j)}+v_k^{(j+1)}?dt$$
其中$dt$是仿真間隔，一般情況下可以取1，如果希望仿真得慢一點，搜索平滑一點，可以適當減小$dt$。
【step3】計算每個粒子的自適應度$F_{k+1}^{(j)}$，為了計算出step1公式中的$P^{(j)}、P_G$和$P_L^{(j)}$。為了方便起見，我們記前k次計算得到了的$P_G$為$P_G^{(k)}$，則第k+1次迭代中我們將適應度最高的粒子位置記為$P_G^{(k+1)}$，則最終的$P_G$為：
$$P_G=?\begin{array}{l}{P}_G^{(k)}，F(P_G^{(k)})>F(P_G^{(k+1)})\\ P_G^{(k+1)}，F(P_G^{(k)})<F(P_G^{(k+1)})\end{array}$$
同樣，我們記前k次計算得到的第$j$個粒子的位置為$P_k^{(j)}$，第k+1次計算得到的第$j$個粒子的位置為$P_{k+1}^{(j)}$，則最終的$P^{(j)}$為：
$$P^{(j)}=?\begin{array}{l}{P}_k^{(j)}，F(P_k^{(j)})>F(P_{k+1}^{(j)})\\ P_{k+1}^{(j)}，F(P_k^{(j)})<F(P_{k+1}^{(j)})\end{array}$$
我們計第k次迭代的$P_L^{(j)}$為$P_{L_k}^{(j)}$，第$k$次迭代的第$i$粒子的位置為$P_k^{(i)}$則
$$P_{L_{k+1}}^{(j)}=avg[min(F(P_k^{(i)}))],||P_k^{(i)}?P_k^{(j)}||\le R$$
即所有距離第$j$個粒子不超過R的粒子中函數值最小的那個。
\qquad注意$P_L^{(j)}$不需要保存為數組，也不需要和上次的值比較，它的值依賴于每次迭代粒子的位置，也會跟著來粒子位置分布的變化為變化。

【step4】更新每個粒子的信息。
由于我們在step2的位置迭代中已經更新過粒子的位置信息，在step1中的速度迭代中已經更新過粒子的速度信息，而在step3中又更新了每個粒子的歷史最優位置$P^{(j)}$及種群最優位置$P_G$，因此實際上如果僅需要知道最優解的情況下我們不需要做這一步。
但如果需要作出參數隨迭代次數的改變的圖的話，每次迭代產生最優解$P_G^{(k)}$及最優適應度$F(P_G^{(k)})$需要用數組保存下來。

2.代碼

2.1.源碼及注釋

\qquad混合粒子群算法與全局粒子群算法的唯一區別就在一多出一個局部社會因子，計算局部社會因子就需要計算每個粒子為中心，R為半徑的“局部”內自適應度最高的粒子。因此判斷粒子與粒子間的距離成為了一大難題。將這部分代碼寫成并行形式，可以大大降低算法時間復雜度和空間復雜度。
\qquad我們定義一個鄰接矩陣
$$D=?\begin{array}{cccc}{d}_{11}& d_{12}& ...& d_{1n}\\ d_{21}& d_{22}& ...& d_{2n}\\ ?& ?& ?& ?\\ d_{n1}& d_{n2}& ...& d_{nn}\end{array}?$$
其中$d_{ij}$代表第$i$個粒子距離第$j$個粒子的距離，$d_{ii}=0$。
每次迭代過程中都計算鄰接矩陣$D$，計算第$j$粒子的局部社會因子時即參考$D$矩陣的第$j$行（或列）的數據即可。
計算鄰接矩陣$D$的函數書寫如下：(Distance_Matrix.m)

function DM=Distance_Matrix(X) m=size(X,2);%樣本數 n=size(X,1);%維數 DM=zeros(n,n);%距離矩陣 for t=1:mDif=ones(n,1)*X(t,:)-X;%坐標差矩陣Dis=sum(Dif.^2,2);%第t個樣本距離其他所有樣本的距離（矩陣）DM(t,:)=Dis; end

注意，這里X是樣本矩陣，樣本是以行向量的形式組成矩陣的。

主程序代碼如下所示：(Hybrid_POS.m)

function [Pg,fmin]=Hybrid_PSO(f,dimension,n,m,xmax,xmin,vmax,vmin,w,c1,c2,R,p,eta,ifdraw) %混合粒子群算法,f為目標函數,dimension為維度,n為代數,m為種群規模 %w,c1,c2分別為慣性權重因子，歷史因素因子，社會因素因子 %位置限幅為[xmin,xmax]，速度限幅為[vmin,vmax] %R為局部粒子群算法作用范圍 %p為局部粒子群占的比重（Max=1） %eta為w的衰減因子，范圍[0,1]，0代表不衰減，1代表epoch全跑完w=0 %ifdraw: True時為繪制粒子位置變化圖，False時不繪制 % 以下四行為位置限幅和速度限幅，對應的四個參數均為行向量 Xmax=ones(m,1)*xmax; Xmin=ones(m,1)*xmin; Vmax=ones(m,1)*vmax; Vmin=ones(m,1)*vmin; Savef=zeros(n+1,1);%記錄自適應度F的數組 if ifdrawSaveData=zeros(m,dimension,10);%記錄11代種群的位置，只需結果時不需要此數組 end dt=0.3;%位移仿真間隔 v=zeros(m,dimension);%初始速度為0 X=(Xmax-Xmin).*rand(m,dimension)+Xmin;%初始位置滿足(-xmax,xmax)內均勻分布 P=X;%P為每個粒子每代的最優位置 Local_Pg=zeros(m,dimension);%局部最優解矩陣 last_f=f(X); [fmin,min_i]=min(last_f);%Pg為所有代中的最優位置 Pg=X(min_i,:);%獲取這個最優位置的坐標 Savef(1)=fmin;N=0; for iter=1:nDM=Distance_Matrix(X);%距離矩陣for i=1:mXR_i=find(DM(i,:)<R);[~,min_i]=min(last_f(XR_i));%按行篩選小于R的點Local_Pg(i,:)=X(XR_i(min_i),:);%局部最優解的坐標endrandM = rand(m,dimension);v=w*v+c1*rand*(P-X)+c2*randM.*((1-p)*(ones(m,1)*Pg-X)+p*(Local_Pg-X));v=(v>Vmax).*Vmax+(v>=Vmin & v<=Vmax).*v+(v<Vmin).*Vmin;X=X+v*dt;X=(X>Xmax).*Xmax+(X>=Xmin & X<=Xmax).*X+(X<Xmin).*Xmin;new_f=f(X);%新的目標函數值update_j=find(new_f<last_f);P(update_j,:)=X(update_j,:);%修正每個粒子的歷史最優值[new_fmin,min_i]=min(new_f);new_Pg=X(min_i,:);Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;last_f=new_f;%保存當前的函數值fmin=min(new_fmin,fmin);%更新函數最小值Savef(iter)=fmin;if mod(iter,floor(n/10))==0%每10代記錄一次種群參數N=N+1;if ifdrawSaveData(:,:,N)=X;endendw=w-iter/n*eta*w; end if ifdrawfor j=1:10figure(j)plot(SaveData(:,1,j),SaveData(:,2,j),'o');xlim([xmin(1),xmax(1)])ylim([xmin(2),xmax(2)])end end figure % 另起一個圖 plot(Savef','b-') title('混合粒子群優化') disp(Pg) disp(fmin) end

注：代碼注釋部分是繪制數次迭代粒子分布位置，要解除需要一起解除（Ctlr+T）,否則會報錯！

2.2.執行與效果

\qquad我們首先來嘗試使用RosenBrock函數（香蕉函數）對其進行試驗。香蕉函數的極小值點為(1,1)，極小值為0。書寫為匿名函數如下：

>> f=@(x)((1-x(:,1)).^2+100*(x(:,2)-x(:,1).^2).^2);

完整的測試例程M文件：

close all clear f=@(x)((1-x(:,1)).^2+100*(x(:,2)-x(:,1).^2).^2); Pg = Hybrid_PSO(f,2,200,50,[10,10],[-10,-10],[3,3],[-3,-3],2,1,1,3,0.3,0.8,true); disp(['Pg:',num2str(Pg)]) disp(['fmin:',num2str(f(Pg))])

運行窗口輸出

>> start_HPSO1.0000 1.00005.9764e-12Pg:1 1 fmin:5.9764e-12

注：這里的冒號不可以省略，因為是樣本為多維行向量組成的矩陣，所以函數也應返回同維向量而不是一個值。

我們設定一個稍大的范圍，速度限幅設置為位置限幅約1/3的大小，局部社會因子取0.3，局部社會因子影響半徑R取位置限幅的1/6。仿真命令及結果如下：
目標函數的下降趨勢如下：

粒子群的變化趨勢如下圖：

①②③

可以發現粒子群算法順利收斂到了(1,1)附近。

希望本文對您有幫助，謝謝閱讀！

總結

以上是生活随笔為你收集整理的【MATLAB】混合粒子群算法原理、代码及详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。