先看兩個概述：

1.

2.

好了，進入主題：

PSO演算法是基於群體智慧理論的優化演算法,群體中的粒子在每次迭代搜尋的過程中,通過跟蹤群體2個極值:粒子本身所找到的最優解Pbest和群體找到的最優解Gbest來動態調整自己位置和速度[5, 6],完成對問題尋優,對於如下的函式優化問題

maxf(x1,x2,…,xn)

s. t  R1j≤xj≤R2j,   j=1,2,…,n(1)

其計算步驟如下[1]:

(1)對粒子群中粒子的位置和速度進行隨機初始化。設定群體的規模為N,則隨機生成如下矩陣

x11x12…x1nv11v12…v1n

x21x22…x2nv21v22…v2n

… … … … … … … …

xN1xN2…xNnvN1vN2…vNn

, (2)

其中{xij, i=1,2,…,N, j=1,2,…,n}表示群體中i粒子的位置為j,vij是對應它的速度,二者均為區間[R1j,R2j]上均勻分佈的隨機數。

(2)計算每個粒子的適應度(目標函式值)。

(3)計算粒子所經歷的最好位置pbesti(t)=(xi1,xi2,…,xim),也就是粒子所經歷過的具有最好適應度的位置,由下面兩個式子確定

pbesti(t+1)= pbesti(t)　如果f(x1(t+1),x2(t+1),…,xn(t+1))<pbesti(t)),

pbesti(t+1)= 

計算群體中所有粒子經歷過最好位置,即全域性最好位置

Gbesti(t)=maxf(Pbest1(t)), f(Pbest2(t)),…, f(PbestN(t)) ,    (4)

(4)依據下式對粒子的速度和位置進行進化

vij(t+1)=vij+c1r1(Pbesti(t)-xij(t))+c2r2(Gbesti(t)-xij(t)),   (5)

其中c1,2為加速度常數(學習速率), r1,2為[0,1]均勻分佈的隨機數。

(5)判斷結束條件,目標函式的適應度達到足夠好或者進化到預先設定的代數,否則返回步驟(2),繼續進行

3用MATLAB實現粒子群演算法

MATLAB是MathWorks公司的產品,是一個功能強大的科學與工程計算軟體,集成了計算、視覺化和程式編制等功能,其使用直觀、簡潔並符合人們思維習慣的程式碼給使用者提供了一個友好的開發環境[9, 10],利用MATLAB矩陣運算的強大功能編寫粒子群演算法程式有很好的優勢。

3. 1粒子群的初始化

initpop函式的功能是實現群體的初始化, popsize表示群體的大小, dimsize表示粒子維數,它由變數的維數決定, xmin是變數的下限, xmax是變數的上限。

function pop=initpop(popsize, dimsize)

pop=unifrnd(xmin, xmax, popsize, 2* dimSize);

3. 2計算粒子的適應度和確定群體的Pbest和Gbest

Calobjvalue函式的功能是計算目標函式的適應度,其公式是採用本文的示例模擬,使用者可以根據不同的優化問題予以修改。

function objvalue=calobjvalue(pop)

for i=1: popsize

obfuct1=sin(sqrt(pop(,i 1) ^2+pop(,i 2) ^2)) ^2-0. 5;

obfuct2=(1. 0+0. 001* (sqrt(x(,i 1) ^2+x(,i 2) ^2)) ^2);

objvalue(,i 1)=0. 5+obfuct1/obfuct2;

end

Pbest=pop(:, 1: dimsize);

[ggBest, xindex]=max(objvalue);

xtemp=pop(xindex, 1: dimsize);

Gbest=xtemp;

3. 3粒子速度和位置的更新

粒子速度和位置的更新是基於式(5)與(6),並生成新的粒子群體.

function pop2=renew(pop)

for t=1: popsize

for dimIndex=1: dimsize

w=wcmax-(wcmax-wcmin)* (generation/maxgeneration);

sub1=Pbest(t, dimIndex)-pop(t, dimIndex);

sub2=Gbest(1, dimIndex)-pop(t, dimIndex);

tempV=w* pop(t, dimszie+dimIndex)+c1* unifrnd(0, 1)* sub1+c2* unifrnd(0, 1)* sub2;

if tempV>speedmax

pop(t, dimszie+dimIndex)=speedmax;

elseif tempV<(-speedmax)

pop(t, dimsize+dimIndex)=-speedmax;

else

pop(t, dimszie+dimIndex)=tempV;

end

tempposition=pop(t, dimIndex)+pop(t, dimsize+dimIndex);

if tempposition>xmax

pop(t, dimIndex)=xmax;

elseif tempposition<xmin

pop(t, dimIndex)=xmin;

else

pop(t, dimIndex)=tempposition;

end

end

3. 4粒子Pbest和Gbest的更新

粒子在進化過程中依據其適應度,調節個體最好位置Pbest和群體最好位置Gbest

function [Pbest,Gbest]=regulate(pop)

for i=1: popsize

obfuct1=sin(sqrt(x(,i 1) ^2+x(,i 2) ^2)) ^2-0. 5;

obfuct2=(1. 0+0. 001* (sqrt(x(,i 1) ^2+x(,i 2) ^2)) ^2);

objvalue(,i 1)=0. 5+obfuct1/obfuct2;

obfuct3=sin(sqrt(pBest(,i 1) ^2+pBest(,i 2) ^2)) ^2-0. 5;

obfuct4=(1. 0+0. 001* (sqrt(pBest(,i 1) ^2+pBest(,i 2) ^2)) ^2);

pvaluer(,i 1)=0. 5+obfuct3/obfuct4;

end

obfuct1=sin(sqrt(gBest(1) ^2+gBest(2) ^2)) ^2-0. 5;

obfuct2=(1. 0+0. 001* (sqrt(gBest(1) ^2+gBest(2) ^2)) ^2);

objvaluetemp=0. 5+obfuct1/obfuct2;

for i=1: popsize

if objvaluer(,i 1)<pvaluer(,i 1)

PBest(,i 1: dimszie)=pop(,i 1: dimszie);

end

if objvaluer(,i 1)<objvaluertemp

GBest=pop(,i 1: dimszie);

xtemp=pop(,i 1: dimsize);

end

end