一、簡介

粒子群優化(PSO)是一種基於群體智慧的數值優化演算法，由社會心理學家James Kennedy和電氣工程師Russell Eberhart於1995年提出。自PSO誕生以來，它在許多方面都得到了改進，這一部分將介紹基本的粒子群優化演算法原理和過程。
1.1 粒子群優化
粒子群優化(PSO)是一種群智慧演算法，其靈感來自於鳥類的群集或魚群學習，用於解決許多科學和工程領域中出現的非線性、非凸性或組合優化問題。

1.1.1 演算法思想
許多鳥類都是群居性的，並由各種原因形成不同的鳥群。鳥群可能大小不同，出現在不同的季節，甚至可能由群體中可以很好合作的不同物種組成。更多的眼睛和耳朵意味著有更多的及時發現食物和捕食者的機會。鳥群在許多方面對其成員的生存總是有益的：
覓食：社會生物學家E.O.Wilson說，至少在理論上，群體中的個體成員可以從其他成員在尋找食物過程中的發現和先前的經驗中獲益[1]。如果一群鳥的食物來源是相同的，那麼某些種類的鳥就會以一種非競爭的方式聚集在一起。這樣，更多的鳥類就能利用其他鳥類對食物位置的發現。
抵禦捕食者：鳥群在保護自己免受捕食者侵害方面有很多優勢。
更多的耳朵和眼睛意味著更多的機會發現捕食者或任何其他潛在的危險；
一群鳥可能會通過圍攻或敏捷的飛行來迷惑或壓制捕食者；
在群體中，互相間的警告可以減少任何一隻鳥的危險。
空氣動力學：當鳥類成群飛行時，它們經常把自己排成特定的形狀或隊形。鳥群中鳥的數量不同，每隻鳥煽動翅膀時產生不同的氣流，這都會導致變化的風型，這些隊形會充分利用不同的分型，從而使得飛行中的鳥類能夠以最節能的方式利用周圍的空氣。
粒子群演算法的發展需要模擬鳥群的一些優點，然而，為了瞭解群體智慧和粒子群優化的一個重要性質，值得提一下是鳥群的一些缺點。當鳥類成群結隊時，也會給它們帶來一些風險。更多的耳朵和眼睛意味著更多的翅膀和嘴，這導致更多的噪音和運動。在這種情況下，更多的捕食者可以定位鳥群，對鳥類造成持續的威脅。一個更大的群體也會需要更多的食物，這導致更多食物競爭，從而可能淘汰群體中一些較弱的鳥類。這裡需要指出的是，PSO並沒有模擬鳥類群體行為的缺點，因此，在搜尋過程中不允許殺死任何個體，而在遺傳演算法中，一些較弱的個體會消亡。在PSO中，所有的個體都將存活，並在整個搜尋過程中努力讓自己變得更強大。在粒子群演算法中，潛在解的改進是合作的結果，而在進化演算法中則是因為競爭。這個概念使得群體智慧不同於進化演算法。簡而言之，在進化演算法中，每一次迭代都有一個新的種群進化，而在群智慧演算法中，每一代都有個體使自己變得更好。個體的身份不會隨著迭代而改變。Mataric[2]給出了以下鳥群規則：

安全漫遊：鳥類飛行時，不存在相互間或與障礙物間的碰撞；
分散：每隻鳥都會與其他鳥保持一個最小的距離；
聚合：每隻鳥也會與其他鳥保持一個最大的距離；
歸巢：所有的鳥類都有可能找到食物來源或巢穴。
在設計粒子群演算法時，並沒有採用這四種規則來模擬鳥類的群體行為。在Kennedy和Eberhart開發的基本粒子群優化模型中，對agent的運動不遵循安全漫遊和分散規則。換句話說，在基本粒子群優化演算法的運動過程中，允許粒子群優化演算法中的代理儘可能地靠近彼此。而聚合和歸巢在粒子群優化模型中是有效的。在粒子群演算法中，代理必須在特定的區域內飛行，以便與任何其他代理保持最大距離。這就相當於在整個過程中，搜尋始終停留在搜尋空間的邊界內或邊界處。第四個規則，歸巢意味著組中的任何代理都可以達到全域性最優。

在PSO模型的發展過程中，Kennedy和Eberhart提出了五個判斷一組代理是否是群體的基本原則：

就近原則：代理群體應該能夠進行簡單的空間和時間計算；
質量原則：代理群體能夠對環境中的質量因素作出反應；
多響應原則：代理群體不應在過於狹窄的通道從事活動；
穩定性原則：代理群體不能每次環境變化時就改變其行為模式；
適應性原則：計算代價不大時，代理群體可以改變其行為模式。

1.1.2 粒子群優化過程
考慮到這五個原則，Kennedy和Eberhart開發了一個用於函式優化的PSO模型。在粒子群演算法中，採用隨機搜尋的方法，利用群體智慧進行求解。換句話說，粒子群演算法是一種群智慧搜尋演算法。這個搜尋是由一組隨機生成的可能解來完成的。這種可能解的集合稱為群，每個可能解都稱為粒子。
在粒子群優化演算法中，粒子的搜尋受到兩種學習方式的影響。每一個粒子都在向其他粒子學習，同時也在運動過程中學習自己的經驗。向他人學習可以稱為社會學習，而從自身經驗中學習可以稱為認知學習。由於社會學習的結果，粒子在它的記憶中儲存了群中所有粒子訪問的最佳解，我們稱之為gbest。通過認知學習，粒子在它的記憶中儲存了迄今為止它自己訪問過的最佳解，稱為pbest。

任何粒子的方向和大小的變化都是由一個叫做速度的因素決定的，速度是位置相對於時間的變化率。對於PSO，迭代的是時間。這樣，對於粒子群演算法，速度可以定義為位置相對於迭代的變化率。由於迭代計數器單位增加，速度v的維數與位置x相同。

對於D維搜尋空間，在時間步t下群體中的第ith個粒子由D維向量x i t = ( x i 1 t , ⋯   , x i D t ) T x_i^t = {(x_{i1}^t, \cdots ,x_{iD}^t)T}xit​=(xi1t​,⋯,xiDt​)T來表示，其速度由另一個D維向量v i t = ( v i 1 t , ⋯   , v i D t ) T v_i^t = {(v_{i1}^t, \cdots ,v_{iD}^t)T}vit​=(vi1t​,⋯,viDt​)T表示。第ith個粒子訪問過的最優解位置用p i t = ( p i 1 t , ⋯   , p i D t ) T p_i^t = {\left( {p_{i1}^t, \cdots ,p_{iD}^t} \right)^T}pit​=(pi1t​,⋯,piDt​)T表示，群體中最優粒子的索引為“g”。第ith個粒子的速度和位置分別由下式進行更新：
v i d t + 1 = v i d t + c 1 r 1 ( p i d t − x i d t ) + c 2 r 2 ( p g d t − x i d t ) (1) v_{id}^{t + 1} = v_{id}^t + {c_1}{r_1}\left( {p_{id}^t - x_{id}^t} \right) + {c_2}{r_2}\left( {p_{gd}^t - x_{id}^t} \right)\tag 1vidt+1​=vidt​+c1​r1​(pidt​−xidt​)+c2​r2​(pgdt​−xidt​)(1)

x i d t + 1 = x i d t + v i d t + 1 (2) x_{id}^{t + 1} = x_{id}^t + v_{id}^{t + 1}\tag 2xidt+1​=xidt​+vidt+1​(2)

其中d=1,2,…,D為維度，i=1,2,…,S為粒子索引，S是群體大小。c1和c2為常數，分別稱為認知和社交縮放參數，或簡單地稱為加速係數。r1和r2是滿足均勻分佈[0,1]之間的隨機數。上面兩個式子均是對每個粒子的每個維度進行單獨更新，問題空間中不同維度之間唯一的聯絡是通過目標函式引入的，也就是當前所找到的最好位置gbest和pbest[3]。PSO的演算法流程如下：
1.1.3 解讀更新等式
速度更新等式（1）的右側包括三部分3：
前一時間的速度v，可以認為是一動量項，用於儲存之前的運動方向，其目的是防止粒子劇烈地改變方向。
第二項是認知或自我部分，通過這一項，粒子的當前位置會向其自己的最好位置移動，這樣在整個搜尋過程中，粒子會記住自己的最佳位置，從而避免自己四處遊蕩。這裡需要注意的是，pidt-xidt是一個方向從xidt到pidt的向量，從而將當前位置向粒子的最佳位置吸引，兩者的順序不能改變，否則當前位置會遠離最佳位置。
第三項是社交部分，負責通過群體共享資訊。通過該項，粒子向群體中最優的個體移動，即每個個體向群體中的其他個體學習。同樣兩者應該是pgbt-xidt。
可以看出，認知尺度引數c1調節的是粒子在其最佳位置方向上的最大步長，而社交尺度引數c2調節的是全域性最優粒子方向上的最大步長。圖2給出了粒子在二維空間中運動的典型幾何圖形。

圖2 粒子群優化過程中粒子移動的幾何說明

從更新方程可以看出，Kennedy和Eberhart的PSO設計遵循了PSO的五個基本原則。在粒子群優化過程中，在d維空間中對一系列時間步進行計算。在任何時間步，種群都遵循gbest和pbest的指導方向，即種群對質量因素作出反應，從而遵循質量原則。由於速度更新方程中有均布隨機數r1和r2，在pbest和gbest之間的當前位置隨機分配，這證明了響應原理的多樣性。在粒子群優化過程中，只有當粒子群從gbest中接收到較好的資訊時，才會發生隨機運動，從而證明了粒子群優化過程的穩定性原則。種群在gbest變化時發生變化，因此遵循適應性原則。
1.2 粒子群優化中的引數
任何基於種群的演算法的收斂速度和尋優能力都受其引數選擇的影響。通常，由於這些演算法的引數高度依賴於問題引數，因此不可能對這些演算法的引數設定給出一般性的建議。但是，已有的理論和/或實驗研究，給出了引數值的一般範圍。與其他基於種群的搜尋演算法類似，由於在搜尋過程中存在隨機因素r1和r2，因此通用PSO的引數調整一直是一項具有挑戰性的任務。PSO的基礎版本只需要很少的引數。本章只討論了[4]中介紹的PSO基礎版本的引數。

一個基本的引數是群體規模，它通常是根據問題中決策變數的數量和問題的複雜性經驗地設定的。一般建議20-50個粒子。

另一個引數是縮放因子c1和c2。如前所述，這些引數決定了下一個迭代中粒子的步長。也就是說，c1和c2決定了粒子的速度。在PSO的基礎版本中，選擇c1=c2=2。在這種情況下，粒子s速度的增加是不受控制的，這有利於更快的收斂速度，但不利於更好地利用搜索空間。如果我們令c1=c2>0，那麼粒子會吸引到pbest和gbest的平均值。c1>c2設定有利於多模態問題，而c2>c1有利於單模態問題。在搜尋過程中，c1和c2的值越小，粒子軌跡越平滑，而c1和c2的值越大，粒子運動越劇烈，加速度越大。研究人員也提出了自適應加速係數[5]。
停止準則不僅是粒子群演算法的引數，也是任何基於種群的元啟發式演算法的引數。常用的停止準則通常基於函式評估或迭代的最大次數，該次數與演算法所花費的時間成正比。一個更有效的停止準則是基於演算法的搜尋能力，如果一個演算法在一定的迭代次數內沒有顯著地改進解，那麼應該停止搜尋。

二、原始碼

%function [gBest]=PSO()%PSO主函式    
pso_size=100;%群體規模
c1=0.5;c2=0.5;%學習因子
w_max=0.8;%最大權重
w_min=0.4;%最小權重
w=w_max;
Pb=10000000;%適應度值 
Pb1=10000000;
Pb2=10000000;
dimens=1;%待優化問題的維數   
run_max=100;%迭代次數上限
X=zeros(1,pso_size);
V=zeros(1,pso_size);
Xb=zeros(1,pso_size);
Xb1=zeros(1,pso_size);
Xb2=zeros(1,pso_size);
Yb=zeros(1,pso_size);
Yb1=zeros(1,pso_size);
Yb2=zeros(1,pso_size);
P=zeros(1,8000);
PL=zeros(1,pso_size);
Pp=0;Pp1=0;Pp2=0;
dpBest=zeros(1,pso_size);
dgBest=0;
gBest=0;
s1=zeros(1,run_max);
s2=zeros(1,run_max);
q=0;
L=0.05;%!!!!!!精度
t=1;% !!!!!!!!!!P(t)中t初值為1
for i=1:pso_size%初始化粒子速度位置
    X(i)=-5+12*rand;
    X1(i)=X(i);
    X2(i)=X(i);
    Xb(i)=X(i);
    Xb1(i)=X(i);
    Xb2(i)=X(i);
    V(i)=(-12+24*rand)/10;
end
for i=1:pso_size   
    Yb1(i)=f1(X(i));%計算f1中個體適應度值,步驟3
    if Yb1(i)<Pb1
        Pb1=Yb1(i);%f1中全域性極值  步驟3,4
        Pp1=X(i);
    end
end
for i=1:pso_size
    Yb2(i)=f2(X(i));%計算f2中個體適應度值
    if Yb2(i)<Pb2
        Pb2=Yb2(i);%f2中全域性極值    步驟3,4
        Pp2=X(i);
    end
end  
gBest=(Pp1+Pp2)/2;%  步驟5
Pp=gBest;
dgBest=abs(Pp1-Pp2);%計算全域性最優值的距離
for i=1:pso_size
    dpBest(i)=abs(Xb1(i)-Xb2(i));%計算各粒子間的距離
    if dpBest(i)<dgBest 
       q=rand;
       if q>0.5
           Xb(i)=Xb1(i);
       else
           Xb(i)=Xb2(i);
       end %個體極值隨機選取    
    else
        Xb(i)=(Xb1(i)+Xb2(i))/2;%取兩函式個體極值的均值
    end
    PL(i)=(Yb1(i)+Yb2(i))/2;%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if PL(i)<L
    
    P(t)=X(i);
    t=t+1;
    end%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end
t%P
for count=1:run_max     %迴圈迭代
    for k=1:pso_size
            V(k)=w*V(k)+c1*rand*(Xb(k)-X(k))+c2*rand*(Pp-X(k));
            if abs(V(k))>12     %判斷是否超出Vmax
                if V(k)>0
                    V(k)=12;
                else
                    V(k)=-12;
                end
            end
             X(k)=X(k)+V(k);
             X1(k)=X(k);
             X2(k)=X(k);
             temp1=f1(X1(k)); %$$$$在f1中計算個體適應度值$$$$
             if temp1<Yb1(k)
                 Yb1(k)=temp1;%$$$$判斷是否更改個體極值$$$$
                 Xb1(k)=X1(k);
             end
             if Yb1(k)<Pb1
                 Pb1=Yb1(k);%$$$$得到新的f1中全域性極值Pb1$$$$
                 Pp1=X1(k);
             end
             temp2=f2(X2(k));  %$$$$在f2中計算個體適應度值$$$$
             if temp2<Yb2(k)
                 Yb2(k)=temp2;%$$$$判斷是否更改個體極值$$$$
                 Xb2(k)=X2(k);
             end
             if Yb2(k)<Pb2
                 Pb2=Yb2(k);  %$$$$
                 Pp2=X2(k);   %$$$$
             end
             gBest=(Pp1+Pp2)/2;    %得到全域性最優
              Pp=gBest;    %全域性極值,
             dgBest=abs(Pp1-Pp2);   %$$$$計算全域性最優值的距離
             dpBest(k)=abs(Xb1(k)-Xb2(k));
             if dpBest(k)<dgBest    %$$$$比較個體極值的距離$$$$
                 q=rand;
             if q>0.5
              Xb(k)=Xb1(k);
            else
             Xb(k)=Xb2(k);
             end    %$$$$個體極值隨機選取  
        else
            Xb(k)=(Xb1(k)+Xb2(k))/2;   %$$$$取兩函式個體極值的均值，得到目標均衡適應度值$$$$
        end
            PL(k)=(Yb1(k)+Yb2(k))/2;
          %flag=0;
        if PL(k)<L%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            a=0;
            for e=1:t-1
                if P(e)==X(k)
                a=1;
                end
            end
            if a~=1
                    P(t)=X(k);
                    t=t+1;
            end
        end
     %   if PL(k)<L        
      %      for e=1:t
       %         if P(t)不等於X(k)
    %P(t)=X(k);
    %t=t+1;
%end

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423