一、粒子群演算法簡介

1 引言
自然界中的鳥群和魚群的群體行為一直是科學家的研究興趣所在。生物學家Craig Reynolds在1987年提出了一個非常有影響的鳥群聚集模型，在他的模擬中，每一個個體都遵循：避免與鄰域個體相撞：匹配鄰域個體的速度；飛向鳥群中心，且整個群體飛向目標。模擬中僅利用上面三條簡單的規則，就可以非常接近地模擬出鳥群飛行的現象。1990年， 生物學家Frank Heppner也提出了鳥類模型， 它的不同之處在於：鳥類被吸引飛到棲息地。在模擬中，一開始每一隻鳥都沒有特定的飛行目標，只是使用簡單的規則確定自己的飛行方向和飛行速度，當有一隻鳥飛到棲息地時，它周圍的鳥也會跟著飛向棲息地，最終整個鳥群都會落在棲息地。
1995年， 美國社會心理學家James Kennedy和電氣工程師RussellEberhart共同提出了粒子群演算法(ParticleS warm Optimization， PSO) ， 該演算法的提出是受對鳥類群體行為進行建模與模擬的研究結果的啟發。他們的模型和模擬演算法主要對Frank Heppner的模型進行了修正，以使粒子飛向解空間並在最優解處降落。粒子群演算法一經提出，由於其演算法簡單，容易實現，立刻引起了進化計算領域學者們的廣泛關注， 形成一個研究熱點。2001年出版的J.Kennedy與R.Eberhart合著的《群體智慧》將群體智慧的影響進一步擴大[] ， 隨後關於粒子群優化演算法的研究報告和研究成果大量湧現，繼而掀起了國內外研究熱潮[2-7]。
粒子群優化演算法來源於鳥類群體活動的規律性，進而利用群體智慧建立一個簡化的模型。它模擬鳥類的覓食行為，將求解問題的搜尋空間比作鳥類的飛行空間，將每隻鳥抽象成一個沒有質量和體積的粒
子，用它來表徵問題的一個可能解，將尋找問題最優解的過程看成鳥類尋找食物的過程，進而求解複雜的優化問題。粒子群優化演算法與其他進化演算法一樣，也是基於“種群”和“進化”的概念，通過個體間
的協作與競爭，實現複雜空間最優解的搜尋。同時，它又不像其他進化演算法那樣對個體進行交叉、變異、選擇等進化運算元操作，而是將群體中的個體看作在l維搜尋空間中沒有質量和體積的粒子，每個粒子以一定的速度在解空間運動， 並向自身歷史最佳位置P best和鄰域歷史最佳位置g best聚集， 實現對候選解的進化。粒子群演算法具有很好的生物社會背景而易於理解，由於引數少而容易實現，對非線性、多峰問題均具有較強的全域性搜尋能力，在科學研究與工程實踐中得到了廣泛關注。目前，該演算法已廣泛應用於函式優化、神經網路訓練、模式分類、模糊控制等領域。

2 粒子群演算法理論
2.1粒子群演算法描述
鳥類在捕食過程中，鳥群成員可以通過個體之間的資訊交流與共享獲得其他成員的發現與飛行經歷。在食物源零星分佈並且不可預測的條件下，這種協作機制所帶來的優勢是決定性的，遠遠大於對食物
的競爭所引起的劣勢。粒子群演算法受鳥類捕食行為的啟發並對這種行為進行模仿，將優化問題的搜尋空間類比於鳥類的飛行空間，將每隻鳥抽象為一個粒子，粒子無質量、無體積，用以表徵問題的一個可行解，優化問題所要搜尋到的最優解則等同於鳥類尋找的食物源。粒子群演算法為每個粒子制定了與鳥類運動類似的簡單行為規則，使整個粒子群的運動表現出與鳥類捕食相似的特性，從而可以求解複雜的優化問題。
粒子群演算法的資訊共享機制可以解釋為一種共生合作的行為，即每個粒子都在不停地進行搜尋，並且其搜尋行為在不同程度上受到群體中其他個體的影響[8]，同時這些粒子還具備對所經歷最佳位置的記
憶能力，即其搜尋行為在受其他個體影響的同時還受到自身經驗的引導。基於獨特的搜尋機制，粒子群演算法首先生成初始種群，即在可行解空間和速度空間隨機初始化粒子的速度與位置，其中粒子的位置用於表徵問題的可行解，然後通過種群間粒子個體的合作與競爭來求解優化問題。
2.2粒子群演算法建模

粒子群演算法首先在給定的解空間中隨機初始化粒子群，待優化問題的變數數決定了解空間的維數。每個粒子有了初始位置與初始速度，然後通過迭代尋優。在每一次迭代中，每個粒子通過跟蹤兩個“極值”來更新自己在解空間中的空間位置與飛行速度：一個極值就是單個粒子本身在迭代過程中找到的最優解粒子，這個粒子叫作個體極值：另一個極值是種群所有粒子在迭代過程中所找到的最優解粒子，這個粒子是全域性極值。上述的方法叫作全域性粒子群演算法。如果不用種群所有粒子而只用其中一部分作為該粒子的鄰居粒子，那麼在所有鄰居粒子中的極值就是區域性極值，該方法稱為區域性粒子群演算法。

2.3粒子群演算法的特點
粒子群演算法本質是一種隨機搜尋演算法，它是一種新興的智慧優化技術。該演算法能以較大概率收斂於全域性最優解。實踐證明，它適合在動態、多目標優化環境中尋優，與傳統優化演算法相比，具有較快的計
算速度和更好的全域性搜尋能力。
(1)粒子群演算法是基於群智慧理論的優化演算法，通過群體中粒子間的合作與競爭產生的群體智慧指導優化搜尋。與其他演算法相比，粒子群演算法是一種高效的並行搜尋演算法。
(2)粒子群演算法與遺傳演算法都是隨機初始化種群，使用適應值來評價個體的優劣程度和進行一定的隨機搜尋。但粒子群演算法根據自己的速度來決定搜尋，沒有遺傳演算法的交叉與變異。與進化演算法相比，粒子群演算法保留了基於種群的全域性搜尋策略，但是其採用的速度-位移模型操作簡單，避免了複雜的遺傳操作。
(3)由於每個粒子在演算法結束時仍保持其個體極值，即粒子群演算法除了可以找到問題的最優解外，還會得到若干較好的次優解，因此將粒子群演算法用於排程和決策問題可以給出多種有意義的方案。
(4)粒子群演算法特有的記憶使其可以動態地跟蹤當前搜尋情況並調整其搜尋策略。另外，粒子群演算法對種群的大小不敏感，即使種群數目下降時，效能下降也不是很大。

3 粒子群演算法種類
3.1基本粒子群演算法

3.2標準粒子群演算法
引入研究粒子群演算法經常用到的兩個概念：一是“探索”，指粒子在一定程度上離開原先的搜尋軌跡，向新的方向進行搜尋，體現了一種向未知區域開拓的能力，類似於全域性搜尋；二是“開發”，指粒子在一定程度上繼續在原先的搜尋軌跡上進行更細一步的搜尋，主要指對探索過程中所搜尋到的區域進行更進一步的搜尋。探索是偏離原來的尋優軌跡去尋找一個更好的解，探索能力是一個演算法的全域性搜尋能力。開發是利用一個好的解，繼續原來的尋優軌跡去搜索更好的解，它是演算法的區域性搜尋能力。如何確定區域性搜尋能力和全域性搜尋能力的比例， 對一個問題的求解過程很重要。1998年， Shi Yuhui等人提出了帶有慣性權重的改進粒子群演算法[10]，由於該演算法能夠保證較好的收斂效果，所以被預設為標準粒子群演算法。其進化過程為：

在式(6.7)中，第一部分表示粒子先前的速度，用於保證演算法的全域性收斂效能；第二部分、第三部分則使演算法具有區域性收斂能力。可以看出，式(6.7)中慣性權重w表示在多大程度上保留原來的速度：W
較大，則全域性收斂能力較強，區域性收斂能力較弱；w較小，則區域性收斂能力較強，全域性收斂能力較弱。
當w=1時，式(6.7)與式(6.5)完全一樣，表明帶慣性權重的粒子群演算法是基本粒子群演算法的擴充套件。實驗結果表明：w在0.8~1.2之間時，粒子群演算法有更快的收斂速度；而當w>1.2時，演算法則容易陷入區域性極值。
另外，在搜尋過程中可以對w進行動態調整：在演算法開始時，可給w賦予較大正值，隨著搜尋的進行，可以線性地使w逐漸減小，這樣可以保證在演算法開始時，各粒子能夠以較大的速度步長在全域性範圍內探
測到較好的區域；而在搜尋後期，較小的w值則保證粒子能夠在極值點周圍做精細的搜尋，從而使演算法有較大的概率向全域性最優解位置收斂。對w進行調整，可以權衡全域性搜尋和區域性搜尋能力。目前，採用較多的動態慣性權重值是Shi提出的線性遞減權值策略， 其表示式如下：

3.3壓縮因子粒子群演算法
Clerc等人提出利用約束因子來控制系統行為的最終收斂[11] ， 該方法可以有效搜尋不同的區域，並且能得到高質量的解。壓縮因子法的速度更新公式為：

實驗結果表明：與使用慣性權重的粒子群優化演算法相比，使用具
有約束因子的粒子群演算法具有更快的收斂速度。
3.4離散粒子群演算法
基本的粒子群演算法是在連續域中搜索函式極值的有力工具。繼基本粒子群演算法之後， Kennedy和Eberhart又提出了一種離散二進位制版的粒子群演算法[12]。在此離散粒子群方法中，將離散問題空間對映到連續粒子運動空間，並適當修改粒子群演算法來求解，在計算上仍保留經典粒子群演算法速度-位置更新運算規則。粒子在狀態空間的取值和變化只限於0和1兩個值， 而速度的每一維vi y代表位置每一位xi取值為1的可能性。因此， 在連續粒子群中的vij更新公式依然保持不變， 但是P best和：best只在[0， 1] 內取值。其位置更新等式表示如下：

4 粒子群演算法流程
粒子群演算法基於“種群”和“進化”的概念，通過個體間的協作與競爭，實現複雜空間最優解的搜尋[13]，其流程如下：
(1)初始化粒子群，包括群體規模N，每個粒子的位置x；和速度Vio
(2) 計算每個粒子的適應度值fit[i] 。
(3) 對每個粒子， 用它的適應度值fit[門和個體極值P best(i)比較。如果fit[i] <P best(i) ， 則用fit[i] 替換掉P best(i) 。
(4) 對每個粒子， 用它的適應度值fit[i] 和全域性極值g best比較。如果fit[i] < 8 best， 則用fit[i] 替換g best。
(5)迭代更新粒子的速度v；和位置xj。
(6)進行邊界條件處理。
(7)判斷演算法終止條件是否滿足：若是，則結束演算法並輸出優化結果；否則返回步驟(2)。
粒子群演算法的運算流程如圖6.1所示。

5 關鍵引數說明
在粒子群優化演算法中，控制引數的選擇能夠影響演算法的效能和效率；如何選擇合適的控制引數使演算法效能最佳，是一個複雜的優化問題。在實際的優化問題中，通常根據使用者的經驗來選取控制引數。
粒子群演算法的控制引數主要包括：粒子種群規模N，慣性權重w，加速係數c和c， 最大速度Via x， 停止準則， 鄰域結構的設定， 邊界條件處理策略等[14]，
粒子種群規模N
粒子種群大小的選擇視具體問題而定，但是一般設定粒子數為20~50。對於大部分的問題10個粒子，已經可以取得很好的結果：不過對於比較難的問題或者特定型別的問題，粒子的數量可以取到100或
200。另外，粒子數目越大，演算法搜尋的空間範圍就越大，也就更容易發現全域性最優解；當然，演算法執行的時間也越長。
慣性權重w
慣性權重w是標準粒子群演算法中非常重要的控制引數，可以用來控制演算法的開發和探索能力。慣性權重的大小表示了對粒子當前速度繼承的多少。當慣性權重值較大時，全域性尋優能力較強，區域性尋優能力
較弱：當慣性權重值較小時，全域性尋優能力較弱，區域性尋優能力較強。慣性權重的選擇通常有固定權重和時變權重。固定權重就是選擇常數作為慣性權重值，在進化過程中其值保持不變，一般取值為
[0.8，1.2]：時變權重則是設定某一變化區間，在進化過程中按照某種方式逐步減小慣性權重。時變權重的選擇包括變化範圍和遞減率。固定的慣性權重可以使粒子保持相同的探索和開發能力，而時變權重可以使粒子在進化的不同階段擁有不同的探索和開發能力。
加速常數c1和c2
加速常數c和c 2分別調節向P best和g best方向飛行的最大步長， 它們分別決定粒子個體經驗和群體經驗對粒子執行軌跡的影響，反映粒子群之間的資訊交流。如果cr=c2=0，則粒子將以當前的飛行速度飛到邊界。此時，粒子僅能搜尋有限的區域，所以難以找到最優解。如果q=0，則為“社會”模型，粒子缺乏認知能力，而只有群體經驗，它的收斂速度較快，但容易陷入區域性最優；如果oy=0，則為“認知”模
型，沒有社會的共享資訊，個體之間沒有資訊的互動，所以找到最優解的概率較小，一個規模為D的群體等價於運行了N個各行其是的粒子。因此一般設定c1=C2，通常可以取c1=cg=1.5。這樣，個體經驗和群體經驗就有了同樣重要的影響力，使得最後的最優解更精確。
粒子的最大速度vmax
粒子的速度在空間中的每一維上都有一個最大速度限制值vd max，用來對粒子的速度進行鉗制， 使速度控制在範圍[-Vimax， +va max] 內，這決定問題空間搜尋的力度， 該值一般由使用者自己設定。Vmax是一個非常重要的引數，如果該值太大，則粒子們也許會飛過優秀區域：而如果該值太小，則粒子們可能無法對區域性最優區域以外的區域進行充分的探測。它們可能會陷入區域性最優，而無法移動足夠遠的距離而跳出區域性最優， 達到空間中更佳的位置。研究者指出， 設定Vmax和調整慣性權重的作用是等效的， 所以!max一般用於對種群的初始化進行設定， 即將vmax設定為每維變數的變化範圍， 而不再對最大速度進行細緻的選擇和調節。
停止準則
最大迭代次數、計算精度或最優解的最大停滯步數▲t(或可以接受的滿意解)通常認為是停止準則，即演算法的終止條件。根據具體的優化問題，停止準則的設定需同時兼顧演算法的求解時間、優化質量和
搜尋效率等多方面效能。
鄰域結構的設定
全域性版本的粒子群演算法將整個群體作為粒子的鄰域，具有收斂速度快的優點，但有時演算法會陷入區域性最優。區域性版本的粒子群演算法將位置相近的個體作為粒子的鄰域，收斂速度較慢，不易陷入區域性最優
值。實際應用中，可先採用全域性粒子群演算法尋找最優解的方向，即得到大致的結果，然後採用區域性粒子群演算法在最優點附近進行精細搜尋。
邊界條件處理
當某一維或若干維的位置或速度超過設定值時，採用邊界條件處理策略可將粒子的位置限制在可行搜尋空間內，這樣能避免種群的膨脹與發散，也能避免粒子大範圍地盲目搜尋，從而提高了搜尋效率。
具體的方法有很多種， 比如通過設定最大位置限制Xmax和最大速度限制Vmax， 當超過最大位置或最大速度時， 在範圍內隨機產生一個數值代替，或者將其設定為最大值，即邊界吸收。

二、案例及完整原始碼

1 案例

2 完整程式碼

%%%%%%%%%%%%%%%%%粒子群演算法求函式極值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %清除所有變數
close all;              %清圖
clc;                    %清屏
N=100;                  %群體粒子個數
D=10;                   %粒子維數
T=200;                  %最大迭代次數
c1=1.5;                 %學習因子1
c2=1.5;                 %學習因子2
w=0.8;                  %慣性權重
Xmax=20;                %位置最大值
Xmin=-20;               %位置最小值
Vmax=10;                %速度最大值
Vmin=-10;               %速度最小值
%%%%%%%%%%%%%%%%初始化種群個體（限定位置和速度）%%%%%%%%%%%%%%%%
x=rand(N,D) * (Xmax-Xmin)+Xmin;
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%%%%%%%%%%%%%%%%%%初始化個體最優位置和最優值%%%%%%%%%%%%%%%%%%%
p=x;
pbest=ones(N,1);
for i=1:N
    pbest(i)=func1(x(i,:));
end
%%%%%%%%%%%%%%%%%%%初始化全域性最優位置和最優值%%%%%%%%%%%%%%%%%%
g=ones(1,D);
gbest=inf;
for i=1:N
    if(pbest(i)<gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T);
%%%%%%%%%%%按照公式依次迭代直到滿足精度或者迭代次數%%%%%%%%%%%%%
for i=1:T
    for j=1:N
        %%%%%%%%%%%%%%更新個體最優位置和最優值%%%%%%%%%%%%%%%%%
        if (func1(x(j,:))<pbest(j))
            p(j,:)=x(j,:);
            pbest(j)=func1(x(j,:));
        end
        %%%%%%%%%%%%%%%%更新全域性最優位置和最優值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);
            gbest=pbest(j);
        end
        %%%%%%%%%%%%%%%%%跟新位置和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));
        x(j,:)=x(j,:)+v(j,:);
        %%%%%%%%%%%%%%%%%%%%邊界條件處理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin;
            end
            if (x(j,ii)>Xmax)  |  (x(j,ii)< Xmin)
                x(j,ii)=rand * (Xmax-Xmin)+Xmin;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%記錄歷代全域性最優值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
g;                         %最優個體         
gb(end);                   %最優值
figure
plot(gb)
xlabel('迭代次數');
ylabel('適應度值');
title('適應度進化曲線')

%%%%%%%%%%%%%%%%%%%%%%%適應度函式%%%%%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;

三、解題過程及執行結果

1 解題過程

2 執行結果

四、matlab版本及參考文獻

1 matlab版本
2014a

2 參考文獻
[1] 包子陽,餘繼周,楊杉.智慧優化演算法及其MATLAB例項（第2版）[M].電子工業出版社，2016.
[2]張巖,吳水根.MATLAB優化演算法原始碼[M].清華大學出版社，2017.