目錄

粒子群優化演算法概述

• 粒子群優化(PSO, particle swarm optimization)演算法是計算智慧領域，除了蟻群演算法，魚群演算法之外的一種群體智慧的優化演算法，該演算法最早由Kennedy和Eberhart在1995年提出的，該演算法源自對鳥類捕食問題的研究。

• PSO演算法首先在可行解空間中初始化一群粒子，每個粒子都代表極值優化問題的一個潛在最優解，用位置、速度和適應度值三項指標表示該粒子特徵。

• 粒子在解空間中運動，通過跟蹤個體極值Pbest和群體極值Gbest更新個體位置，個體極值Pbest是指個體所經歷位置中計算得到的適應度值最優位置，群體極值Gbest是指種群中的所有粒子搜尋到的適應度最優位置。

• 粒子每更新一次位置，就計算一次適應度值，並且通過比較新粒子的適應度值和個體極值、群體極值的適應度值更新個體極值Pbest和群體極值Gbest位置。

在每一次迭代過程中，粒子通過個體極值和群體極值更新自身的速度和位置，更新公式如下:

上式中V是速度，以當前的速度加上兩個修正項：該個體與當前行進路徑中最優個體的差距、與群體最優值的偏差。c1和c2是係數，r1、r2是隨機數。

X是位置，用當前位置加上速度。根據相鄰時間間隔預設為一個時間單位，故速度和距離可以直接相加。

PSO演算法步驟

步驟中的關鍵點提示：

1）初始化、適應度函式計算和遺傳演算法很相似。

2）群體極值是好找的。但個體極值：第一步時，每個個體的值都是個體極值，第二步後才開始有真正的個體極值的概念。

3）終止條件：比如說，達到多少次迭代次數，相鄰兩次誤差小於一定值，等等。也可以多種終止條件混合使用。

PSO（粒子群優化演算法）與GA（遺傳演算法）對比

• 相同點:

1）種群隨機初始化，上面也提到了。

2）適應度函式值與目標最優解之間：都有一個對映關係

• 不同點:

1）PSO演算法沒有選擇、交叉、變異等操作運算元。取而代之的是個體極值、群體極值來實現逐步優化的功能。GA的相關操作含義參見《遺傳演算法原理簡介及其MATLAB實踐》。

2）PSO有記憶的功能：

在優化過程中參考到了上一步的極值情況（劃橫線的部分）。按此公式計算出新的粒子位置時，若新粒子的適應度函式還不如之前的好，則這個優化方法會幫助優化程序回到之前的位置，體現了一個記憶的效果。矩形內是權重係數。

3）資訊共享機制不同，遺傳演算法是互相共享資訊，整個種群的移動是比較均勻地向最優區域移動，而在PSO中，只有gBest或lBest給出資訊給其他粒子，屬於單向的資訊流動，整個搜尋更新過程是跟隨當前最優解的過程。因此， 在一般情況下，PSO的收斂速度更快。

簡單直白地理解：GA在變異過程中可能從比較好的情況又變成不好的情況，不是持續收斂的過程，所以耗時會更長些。

PSO的MATLAB實現

MATLAB2014以上有自帶的PSO的工具箱函式。

這裡我們自己寫程式碼來實現一下PSO：

【例】利用PSO尋找極值點：y = sin(10*pi*x) / x;

%% I. 清空環境
clc
clear all

%% II. 繪製目標函式曲線圖
x = 1:0.01:2;
y = sin(10*pi*x) ./ x;
figure
plot(x, y)
hold on

%% III. 引數初始化
c1 = 1.49445;
c2 = 1.49445;

maxgen = 50;   % 進化次數  
sizepop = 10;   %種群規模

Vmax = 0.5;   %速度的範圍，超過則用邊界值。
Vmin = -0.5;  
popmax = 2;   %個體的變化範圍
popmin = 1;

%% IV. 產生初始粒子和速度
for i = 1:sizepop
    % 隨機產生一個種群
    pop(i,:) = (rands(1) + 1) / 2 + 1;    %初始種群，rands產生(-1,1)，調整到(1,2)
    V(i,:) = 0.5 * rands(1);  %初始化速度
    % 計算適應度
    fitness(i) = fun(pop(i,:));   
end

%% V. 個體極值和群體極值
[bestfitness bestindex] = max(fitness);
zbest = pop(bestindex,:);   %全域性最佳
gbest = pop;    %個體最佳
fitnessgbest = fitness;   %個體最佳適應度值
fitnesszbest = bestfitness;   %全域性最佳適應度值

%% VI. 迭代尋優
for i = 1:maxgen
    
    for j = 1:sizepop
        % 速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        V(j,find(V(j,:)>Vmax)) = Vmax;
        V(j,find(V(j,:)<Vmin)) = Vmin;
        
        % 種群更新
        pop(j,:) = pop(j,:) + V(j,:);
        pop(j,find(pop(j,:)>popmax)) = popmax;
        pop(j,find(pop(j,:)<popmin)) = popmin;
        
        % 適應度值更新
        fitness(j) = fun(pop(j,:)); 
    end
    
    for j = 1:sizepop    
        % 個體最優更新
        if fitness(j) > fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        
        % 群體最優更新
        if fitness(j) > fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end 
    yy(i) = fitnesszbest;          
end

%% VII. 輸出結果並繪圖
[fitnesszbest zbest]
plot(zbest, fitnesszbest,'r*')

figure
plot(yy)
title('最優個體適應度','fontsize',12);
xlabel('進化代數','fontsize',12);ylabel('適應度','fontsize',12);