PSO（粒子群）演算法學習筆記

阿新 • • 發佈：2021-01-27

PSO（粒子群）演算法學習筆記

PSO演算法的概述

PSO演算法是一種全域性優化的演算法，它模擬的是鳥群或者魚群的一種彼此共享資訊去搜尋食物的過程。
PSO演算法與遺傳演算法(GA)類似，但是其少了GA中“交叉”，“變異”的操作，因此總體來看PSO演算法操作更加簡單，更加容易實現
PSO演算法可以理解為使用一群粒子來模擬鳥群個體的變化，每一個粒子在每一輪都會進行下一位置的搜尋，下一輪粒子飛行的方向由當前時刻的方向以及自身搜尋到的最優解方向和目前搜尋到的全域性最優解方向決定，粒子飛行的位置與當前位置與下一輪搜尋的方向共同決定。

PSO演算法最核心的兩個公式如下：

V k n e x t = w V k n o w + C 1 R 1 ( P b e s t k − S k n o w ) + C 2 R 2 ( g b e s t − S k n o w ) V^{next}_k=wV_k^{now}+C_1R_1(P_{best}^{k}-S_k^{now})+C_2R_2(g_{best}-S^{now}_k) Vknext=wVknow+C1R1(Pbestk−Sknow)+C2R2(gbest−Sknow)

S k n e x t = S k n o w + V k n e x t S_k^{next}=S_k^{now}+V_k^{next}

Sknext=Sknow+Vknext

公式變數說明：
V k n e x t V_k^{next} Vknext代表了下一輪次搜尋的方向，
V k n o w V_k^{now} Vknow代表了當先搜尋的方向資訊。
S k n o w S^{now}_k Sknow代表的是當前第K個粒子搜尋的位置資訊， S k n e x t S_k^{next} Sknext代表的第K個下一次搜尋到位置資訊。 P b e s t k P_{best}^{k} Pbestk代表的是當前第K個粒子搜尋到的個體極值的位置資訊
g b e s t g_{best} gbest代表了當前全域性最優的位置資訊。

C 1 C_1 C1是個體學習因子
C 2 C_2 C2是社會學習因子
w w w是方向向量的慣性權重值。
公式原理說明：第K個物件擁有自身搜尋的最優解資訊 P b e s t k P^k_{best} Pbestk,這是從其以往的經驗中 得出來的資訊；擁有目前為止的全域性最優解 g b e s t g_{best} gbest的資訊，這是粒子之間共享的資訊，為了從當前位置得到下一位置資訊 S k n e x t S_k^{next} Sknext，粒子必須知道心得方向資訊 V k n e x t V_k^{next} Vknext,此時PSO演算法要求粒子必須在一方面維持當前方向 V k n o w V_k^{now} Vknow的慣性，另一方面要考慮指向物件最優解的方向 P b e s t k − S k n o w P_{best}^{k}-S_k^{now} Pbestk−Sknow ,以及指向全域性最優解的方向 g b e s t − S k n o w g_{best}-S^{now}_k gbest−Sknow 。這三者綜合可以得到公式（1），而公式（2）描述的是粒子位置資訊的改變方式，即每一個粒子在原先位置的基礎上，加上新的方向資訊，得到新的位置資訊

PSO演算法的步驟

首先利用隨機數初始化所有粒子的初始位置S與初始方向向量V，並初始化相關引數C1,C2,W等
設定成本函式，並以成本函式最小（大）作為全域性變數的搜尋目標
進行一輪搜尋，記錄每一個粒子歷史搜尋的最優值，最優位置
結束一輪搜尋後，比較當前輪次最小（大）成本與全域性最小（大）成本，滿足條件則替換
根據公式（1）（2）迭代得到下一輪次的搜尋方向與搜尋位置
重複執行（3），直達達到迭代上限或者得到滿意的解

一個栗子

這裡採用的是Branin函式（該函式常被用於二維的測試函式），Branin函式定義如下：

在這裡插入圖片描述
Branin函式的全域性最優解為0.397887358，有三個全域性最優解
這裡我們使用PSO演算法求解，程式碼如下：

clf;
PARAMETER=2;
v_max=5*ones(1,PARAMETER);   % 速度限制 上限
v_min=-5*ones(1,PARAMETER);  % 速度限制  下限
up_scope=[15 10];      % x y 座標範圍限制  上限
down_scope=[-15,-10];  % x y 座標的範圍限制  下限
C1=0.5;% 個體學習因子
C2=0.5;% 社會學習因子
w=0.6; %方向向量的慣性權重值
MAX_ITER=800;  % 最大迭代次數
PARTICLE_NUMBER=30;
%% 初始化位置 方向
for i=1:PARTICLE_NUMBER
    S(i,:)=(-down_scope+up_scope).*(rand(PARAMETER,1))'+down_scope;  % 位置
    V(i,:)=(v_max-v_min).*(rand(PARAMETER,1))'+v_min;   % 方向
end
%物件最優解
cost_personal=(1e10.*ones(PARTICLE_NUMBER,1));
pos_personal=S;
%全域性最優解
cost_global=100;
pos_global=[];
% 最終結果
TRUE=0.397887358;
err=1000;
iter=0;
while (iter<=MAX_ITER && err>1e-7)
    iter=iter+1; % 迭代加一
    %  一輪中各個個體對應的代價最小值  位置值
    for i=1:PARTICLE_NUMBER
        cost(i,1)=TestBR(S(i,:));
        if cost(i,1)<cost_personal(i,1)
            cost_personal(i,1)=cost(i,1);   % 個體代價更新
            pos_personal(i,:)=S(i,:);  % 位置資訊更新
        end
    end
    [val,position]=min(cost_personal);
    if (val<cost_global)
        cost_global=val;
        pos_global=pos_personal(position,:);
    end
    V=w*V+C1*rand*(pos_personal(:,:)-S(:,:))+C2*rand*(ones(30,1)*(pos_global)-S(:,:));
    S=S+V;
    err=abs(cost_global-TRUE);
    fee(iter)=val;

    figure(1);
    plot(pos_personal(:,1),pos_personal(:,2),".r",'Markersize',15);
    axis([down_scope(1) up_scope(1) down_scope(2) up_scope(2)]);
    set(get(gca, 'XLabel'), 'String', 'x座標');
    set(get(gca, 'YLabel'), 'String', 'Y座標');
    set(get(gca, 'Title'), 'String', 'PSO演算法尋優');
    pause(0.1);
end
disp("min");
disp(cost_global);
disp ("min location");
disp(pos_global);
figure(2);
plot(1:length(fee),fee,'b');
set(get(gca, 'XLabel'), 'String', '迭代次數');
set(get(gca, 'YLabel'), 'String', '當前輪次搜尋最優值');
hold on ;

function y=TestBR(x)  %的時刻
    x1=x(1);   x2=x(2);
    y=(x2-5.1/(4*pi^2)*x1^2 +5/pi*x1-6)^2+10*(1-1/(8*pi))*cos(x1)+10;
end