MATLAB粒子群優化演算法（PSO）

作者：凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/

一、介紹

• 粒子群優化演算法（Particle Swarm Optimization Algorithm）是一種群智能演算法，為了尋求全局最優。群體迭代，粒子在解空間追隨最優的粒子進行搜索。
• 粒子群演算法的思想源於對鳥群捕食行為的研究，模擬鳥集群飛行覓食的行為，鳥之間通過集體的協作使群體達到最優目的，是一種基於Swarm Intelligence的優化方法。
• 馬良教授在他的著作《蟻群優化演算法》一書的前言中寫到：“自然界的蟻群、鳥群、魚群、羊群、牛群、蜂群等，其實時時刻刻都在給予我們以某種啟示，只不過我們常常忽略了大自然對我們的最大恩賜！”
• 設想這樣一個場景：一群鳥在隨機搜索食物，已知
• （1）在這塊區域里只有一塊食物;  （2）所有的鳥都不知道食物在哪裡;  （3）但它們能感受到當前的位置離食物還有多遠.  
• 那麼:找到食物的最優策略是什麼呢？ 
• 搜尋目前離食物最近的鳥的周圍區域，根據自己飛行的經驗判斷食物的所在。

• PSO的基礎: 信息的社會共用 

二、演算法思路

1. 每個尋優的問題解都被想像成一隻鳥，稱為“粒子”。所有粒子都在一個D維空間進行搜索。
2. 所有的粒子都由一個fitness function 確定適應值以判斷目前的位置好壞。
3. 每一個粒子必須賦予記憶功能，能記住所搜尋到的最佳位置。
4. 每一個粒子還有一個速度以決定飛行的距離和方向。這個速度根據它本身的飛行經驗以及同伴的飛行經驗進行動態調整。

三、演算法流程與偽代碼

• c1: weight of local information
• c2: weight of global information
• p: particle's position
• v: path direction
• rand: random variable, 0<rand<1
• number of particles: 10~50
• c1+c2=4
• v=v+c1*rand*(pbest-p)+c2*rand*(gbest-p);中右邊第一項v: inertia, c1*rand*(pbest-p): personal influence, c2*rand*(gbest-p): social influence.

四、MATLAB程式

粒子群優化演算法求解函數y=x+x*sin(x)在[-10, 10]的最大值.

pso.m

function [best_x, best_y]=pso()
%find the max value of a function
%Author: kailugaji https://www.cnblogs.com/kailugaji/
posMax=10;posMin=-10;  %range of feasible solution
x=posMin:0.1:posMax;
y=my_fun(x);  %object function
plot(x,y);hold on; %plot function
popsize=30;  %number of particles
max_iter=100;  %maximum number of iterations
position=rand(popsize,1);
position=position*20-ones(popsize,1)*10;   %particle's position
vMax=1;vMin=-1;  %range of velocity
velosity=2*rand(popsize,1)-1;  %path direction
pbest=position;  %current optimum
gbest=position(1);  %global optimum
wBegin=1;wEnd=0.05;
c1=2;  %weight of local information
c2=2;  %weight of global information
for i=1:max_iter
    position_new=position+velosity;  %change direction of each particle
    velosity_new=(wEnd+(wBegin-wEnd)*(100-i)/100).*velosity+c1*rand().*(position-pbest.*ones(popsize,1))+c2*rand().*(position-gbest.*ones(popsize,1));  %inertia+personal influence+social influence
    for j=1:popsize
        if(position_new(j)>posMax)
            position_new(j)=posMax;
        end
        if(position_new(j)<posMin)
            position_new(j)=posMin;
        end
        if(velosity_new(j)>vMax)
            velosity_new(j)=vMax;
        end
        if(velosity_new(j)<vMin)
            velosity_new(j)=vMin;
        end
         if(my_fun(position_new(j))>my_fun(pbest(j)))
            pbest(j)=position_new(j);
        end
        if(my_fun(position_new(j))>my_fun(gbest))
            gbest=position_new(j);
        end
    end
    position=position_new;
    velosity=velosity_new;
end
plot(gbest,my_fun(gbest),'or');
best_x=gbest;
best_y=my_fun(gbest);

my_fun.m

function y=my_fun(x)
y=x+x.*sin(x);

五、實驗結果

>> [best_x, best_y]=pso()
best_x =
    8.0051
best_y =
   15.9190

當x=8.0051時，粒子群優化演算法終止迭代，此時獲得的最大值為y=15.9190。

六、參考文獻

• 粒子群優化演算法
• 優化 | 粒子群演算法介紹