粒子群算法

一. PSO算法的基本理論

粒子群（PSO）算法是依托群鳥覓食的模型尋找最優解。群體特征基于以下3個簡單的行為準則：

• 沖突避免：群體在一定空間移動，個體有自己的移動意志，但不能影響影響其他個體移動，避免碰撞和争執。
• 速度比對：個體必須配合中心移動速度，不管在方向、距離和速率上都必須互相配合。
• 群體中心：個體将會向群體中心移動，配合群體中心向目标前進。

群鳥覓食其實是一個最佳決策的過程， 與人類決策的過程相似。Boyd和Re chars on探索了人類的決策過程，并提出了個體學習和文化傳遞的概念。根據他們的研究成果，人們在決策過程中常常會綜合兩種重要的資訊：

• 第一種是自己的經驗，即他們根據以前自己的嘗試和經曆，已經積累了一定的經驗，知道怎樣的狀态會比較好。
• 第二種是他人的經驗，即從周圍人的行為擷取知識，從中知道哪些選擇是正面的，哪些選擇是消極的。

同理，群鳥在覓食的過程中，每隻鳥的初始狀态都處于随機位置，且飛翔的方向也是随機的。每隻鳥都不知道食物在哪裡，但是随着時間的推移，這些初始處于随機位置的鳥類通過群内互相學習、資訊共享和個體不斷積累字覓食物的經驗，自發組織積聚成一個群落，并逐漸朝唯一的目标-—食物前進。每隻鳥能夠通過一定經驗和資訊估計目前所處的位置對于能尋找到食物有多大的價值，即多大的适應值；每隻鳥能夠記住自己所找到的最好位置，稱之為局部最優(pbest) 。此外， 還能記住群鳥中所有個體所能找到的最好位置， 稱為全局最優(gbest) ， 整個鳥群的覓食中心都趨向全局最優移動， 這在生物學上稱之為“同步效應”。通過鳥群覓食的位置不斷移動，即不斷疊代，可以使鳥群朝食物步步進逼。

在群鳥覓食模型中，每個個體都可以被看成一個粒子，則鳥群可以被看成一個粒子群。假設在一個 D D D維的目标搜尋空間中，有 m m m個粒子組成一個群體，其中第 i i i個粒子( i = 1 , 2 , ⋯   , m i= 1,2,\cdots,m i=1,2,⋯,m)，位置表示為 X i X_i Xi​ = ( x i 1 , ⋯   , x i D x^1_{i},\cdots,x^D_{i} xi1​,⋯,xiD​)，即第 i i i個粒子在 D D D維搜尋空間中的位置為 X i X_i Xi​。換言之，每個粒子的位置就是一個潛在解，将 X i X_i Xi​代入目标函數就可以計算出其适應值，根據适應值的大小衡量其優劣。粒子個體經曆過的最好位置記為 P i = ( p i 1 , ⋯   , p i D ) P_i = (p^1_{i},\cdots,p^D_{i}) Pi​=(pi1​,⋯,piD​)，整個群體所有粒子經曆過的最好位置記為 P g = ( p g 1 , ⋯   , p g D ) P_g = (p^1_{g},\cdots,p^D_{g}) Pg​=(pg1​,⋯,pgD​)。粒子 i i i的速度記為 V i = ( v i 1 , ⋯   , v i D ) V_i = (v^1_{i},\cdots,v^D_{i}) Vi​=(vi1​,⋯,viD​)。

粒子群算法采用下列公式對粒子所在的位置不斷更新(機關時間1)：

v i d = w v i d + c 1 r 1 ( p i d − x i d ) + c 2 r 2 ( p g d − x i d ) (1) {v^d_{i}=wv^d_{i}+c_1r_1(p^d_i-x^d_i)+c_2r_2(p_g^d-x_i^d) \tag{1}} vid​=wvid​+c1​r1​(pid​−xid​)+c2​r2​(pgd​−xid​)(1) x i d = x i d + α v i d (2) x_i^d=x_i^d+\alpha v_i^d \tag{2} xid​=xid​+αvid​(2)

• 自 我 認 知 ： c 1 r 1 ( p i d − x i d ) 自我認知：c_1r_1(p^d_i-x^d_i) 自我認知：c1​r1​(pid​−xid​)
• 社 會 經 驗 ： c 2 r 2 ( p g d − x i d ) 社會經驗：c_2r_2(p_g^d-x_i^d) 社會經驗：c2​r2​(pgd​−xid​)
• i = 1 , 2 , ⋯   , m i = 1,2,\cdots,m i=1,2,⋯,m
• d = 1 , 2 , ⋯   , D d = 1,2,\cdots,D d=1,2,⋯,D
• w w w是非負數，稱為慣性因子
• 加速度常數 c 1 、 c 2 c_1、c_2 c1​、c2​是非負數
• r 1 、 r 2 r_1、r_2 r1​、r2​是[0,1]範圍内變換的随機數
• α \alpha α稱為限制因子，目的是控制速度的權重

v i d ∈ [ − v m a x d ， v m a x d ] v_i^d\in[-v^d_{max}，v^d_{max}] vid​∈[−vmaxd​，vmaxd​]範圍内，即粒子 i i i的飛行速度 V i V_i Vi​被一個最大速度 V m a x = ( v m a x 1 , ⋯   , v m a x D ) V_{max}=(v_{max}^1,\cdots,v_{max}^D) Vmax​=(vmax1​,⋯,vmaxD​)所限制。如果目前時刻粒子在某維的速度 v i d v_i^d vid​更新後超過該維的最大飛翔速度 v m a x d v_{max}^d vmaxd​則目前時刻該維的速度被限制在 v m a x d v_{max}^d vmaxd​。 v m a x v_{max} vmax​ 為常數，可以根據不同的優化問題設定。

疊代終止條件根據具體問題設定，一般達到預定最大疊代次數或粒子群且前為止搜尋到的最優位置滿足目标函數的最小允許誤差。

二. PSO算法程式設計流程

• 步驟1：初始化粒子群：速度、位置、慣性因子、加速常數、最大疊代次數、算法終止的最小允許誤差。
• 步驟2：評價每個粒子的初始适應值。
• 步驟3：将初始适應值作為目前每個粒子的局部最優值，并将各适應值對應的位置作為每個粒子的局部最優值所在的位置。
• 步驟4：将最佳初始适應值作為目前全局最優值，并将最佳适應值對應的位置作為全局最優值所在的位置。
• 步驟5：依據式(1)更新每個粒子目前的飛翔速度。
• 步驟6：對每個粒子的飛翔速度進行限幅處理，使之不能超過設定的最大飛翔速度。
• 步驟7：依據式(2)更新每個粒子目前所在的位置。
• 步驟8：比較目前每個粒子的适應值是否比曆史局部最優值好，如果好，則将目前粒子适應值作為粒子的局部最優值，其對應的位置作為每個粒子的局部最優值所在的位置。
• 步驟9：在目前群中找出全局最優值，并将目前全局最優值對應的位置作為粒子群的全局最優值所在的位置。
• 步驟10：重複步驟5~9，直到滿足設定的最小誤差或最大疊代次數。
• 步驟11：輸出粒子群的全局最優值和其對應的位置以及每個粒子的局部最優值和其對應的位置。

三. MATLAB程式設計實作

function main()
clc;clear all; close all;
tic;                                    % 程式運作計時
E0 = 0.001;                             % 允許誤差
MaxNum = 100;                           % 粒子最大疊代次數
narvs = 1;                              % 目标函數的自變量個數
particlesize = 30;                      % 粒子群規模
c1 = 2;                                 % 個體經驗學習因子
c2 = 2;                                 % 社會經驗學習因子
w =0.6;                                 % 慣性因子
vmax = 0.8;                             % 粒子的最大飛翔速度,rand(m,n)産生m*n的随機矩陣，範圍(0,1)
x = -5 + 10 * rand(particlesize, narvs);% 粒子所在的位置 (rand産生的大小為0,1),規模是 粒子群數和參數需求數 設定了x的取值範圍[-5,5] 
v = 2*rand(particlesize,narvs);         % 粒子的飛翔速度,生成每個粒子的飛翔速度，由于是隻有一個變量，是以速度是一維的

% 用inline定義适應度函數以便将子函數檔案與主程式檔案放到一起
% 目标函數是:y = 1+(2.1*(1- x + 2*x.^2).*exp(-x.^2 / 2)) # 與Python不同的是，這裡必須要寫成.*
% .^之類的，因為定義不同
%inline指令定義适應度函數如下：

fitness = inline('1/(1+(2.1*(1 - x + 2 * x.^2).*exp(-x.^2/2)))','x'); % 這裡求倒數，還在分母上加了個1，確定不會出現分母為0的情況，轉為求最小值位置

% inline函數定義可以大大降低程式運作速度

for i= 1:particlesize
    for j = 1:narvs
        f(i) = fitness(x(i,j));     %i是粒子個數，j是次元就是變量個數，x(i,j)就是表示的x的值
    end
end 
% 完成了對每一個粒子，在各自位置上的适應值
% 粒子開始學習
personalbest_x=x;         % 用于存儲對于每一個粒子最佳經曆點的x值
personalbest_faval=f;     % 同時存儲對于每一個粒子的最佳經曆點的數值，用于更改 
[globalbest_faval,i] = min(personalbest_faval); % min函數傳回的第一個是最小值，還有一個就是最小值的下标，這裡就是告訴了是在哪個粒子上
globalbest_x = personalbest_x(i,:);   % 這個是必定是全局最優點的位置
k = 1; % 開始疊代計數
while k <= MaxNum   % 當疊代次數達到設定的最大值的時候，就不要再進行疊代了
    for i = 1:particlesize   % 對于每一個粒子
         for j = 1:narvs
            f(i) = fitness(x(i,j)); % 得到每個粒子的目前位置 在函數上的适應值 
         end
        if f(i) < personalbest_faval(i)   % 如果這個值是小于個人最優解的位置的時候，就更新,我們經過轉換，是以隻用考慮求最小值的情況
            personalbest_faval(i) = f(i); % 将第i個粒子的個人最優解設定為全局最優解
            personalbest_x(i,:) = x(i,:); % 同時更改最優解的位置
        end
    end 
   [globalbest_faval,i] = min(personalbest_faval); 
    globalbest_x = personalbest_x(i,:); % 更新全局 全局資訊由個體資訊描述組成

    for i = 1:particlesize              %更新粒子群裡每個個體的最新位置
        v(i,:) = w*v(i,:) + c1*rand*(personalbest_x(i,:) - x(i,:)) + c2*rand*(globalbest_x-x(i,:)); % 由個人和全局的最佳資訊資料進行更新移動速度
        % 上面中rand會随機生成一個rand(0,1)然後會随機的降低學習因子的比例
        for j = 1:narvs   % 這個個循環确定了每個自變量上的速度，有沒有超過對應的最大值
            if v(i,j) > vmax;
                v(i,j) = vmax;
            else if v(i,j) < -vmax;
                v(i,j) = -vmax;
            end
        end 
        x(i,:) = x(i,:) + v(i,:); % 通過速度更新位置
 end
    if abs(globalbest_faval) < E0,break,end         
    k = k + 1;
end
Value1 = 1/globalbest_faval - 1; % 還記得上面做了一個加1，求倒數的操作麼?
Value1 = num2str(Value1);
%strcat指令可以實作字元的組合輸出
disp(strcat('the maximum value',' = ', Value1)); % 主要是在這進行了展示
%輸出最大值所在的橫坐标的位置
Value2 = globalbest_x;    % 得到了全局最優解的位置
Value2 = num2str(Value2);
disp(strcat('the maximum x = ', Value2));

% 繪圖
x = -5:0.01:5;
y = 2.1*(1 - x + 2 * x.^2).*exp(-x.^2/2);
plot(x,y,'m--','linewidth',3); % m表示的是粉紅色,-是表示的是連續的曲線線
hold on;
plot(globalbest_x, 1/globalbest_faval-1,'kp','linewidth',4);
legend('目标函數','搜尋到的最大值');
xlabel('x'); % 給x軸貼标簽
ylabel('y'); % 給y軸貼标簽
grid on;
end
           

四. 算法舉例