算法篇----簡易版的粒子群（PSO）實作（Matlab語言）

前言

    粒子群算法實作起來并不是很難，算法思想可以參加我上一篇博文，不多說了。好了，Matlab版的粒子群走起。

1 定義變量

    粒子群算法有很多參數，做實驗的時候會糾結在參數問題上，這裡就随機設定了。有時候，參數好壞，是成敗關鍵。沒有修過Matlab語言的朋友不用擔心，之後會把C或Java版的也總結一下。知識嘛，多多益善。Let's  go ！

    根據公式，需要c1 c2 r1 r2 w M 

    代碼如下：

c1 = 2;%學習因子

c2 = 2;%學習因子

w = 0.7;%慣性權重

D = 10;%搜尋次元

M = 40;%種群大小

%r1、r2 由randn生成，不聲明了

2 初始化種群

    這裡是随機生成在不同次元上的粒子，就是說在D個方向上撒上M個粒子。

    代碼如下：

for i=1:M

for j=1:D

x(i,j)=randn; %随機初始化位置

v(i,j)=randn; %随機初始化速度

end

end

3 計算粒子适應度

    通過調用适應度函數，計算适應度值。在此處，初始化個體的适應度值p(i)和全局最優的适應度值gbest。

    代碼如下：

for i=1:M

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

gbest=x(1,:); %gbest為全局最優适應度值

4 更新個體最優值

    将适應度函數計算的結果和目前的個體适應度值相比較，将适應度值較小的位置儲存下來。

    代碼如下：

for i = 1:M

if fitness(x(i,:),D) < p(i)

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

end

5 更新全局最優值

    上面已經将全局最優位置初始化，隻需要将它和目前個體的最優位置相比較即可，位置代入适應度函數，取最小的值對應的位置，儲存該位置即可。

    代碼如下：

if p(i) < fitness(pg,D)

gbest=y(i,:);

end

6 更新粒子的速度和位置

    根據公式來寫即可，對于不同的粒子群算法這個地方會需要改動。

    代碼如下：

for i=1:M

v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));

x(i,:)=x(i,:)+v(i,:);

end

7 适應度函數

    适應度函數最好另外建一個Script檔案，把它放進去。适應度函數可以根據要研究的問題做相應的改動，下面列出一個函數，僅供測試。

    代碼如下:

function result = fitness(X,D)

sum = 0.0;

for i = 1:D

sum = sum + X(i)^2;

end

result = sum;

8  完整代碼

    到目前，應該有兩個m檔案，一個是适應度函數的，一個是粒子群的。不過，還需要在粒子群的檔案中加入一些代碼，主要是調整代碼順序，疊代次數和精度問題，具體參見代碼即可。

    完整代碼如下：

 clear all;

clc;

format long;

%------給定初始化條件----------------------------------------------

c1=2; %學習因子1

c2=2; %學習因子2

w=0.7; %慣性權重

MaxDT=1000; %最大疊代次數

D=10; %搜尋空間維數（未知數個數）

M=40; %初始化群體個體數目

eps=10^(-6); %設定精度(在已知最小值時候用)

%------初始化種群的個體(可以在這裡限定位置和速度的範圍)------------

for i=1:M

for j=1:D

x(i,j)=randn; %随機初始化位置

v(i,j)=randn; %随機初始化速度

end

end

%------先計算各個粒子的适應度，并初始化p(i)和gbest--------------------

for i=1:M

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

gbest=x(1,:); %gbest為全局最優

for i=2:M

if fitness(x(i,M:),D) < fitness(gbest,D)

gbest=x(i,:);

end

end

%------進入主要循環，按照公式依次疊代，直到滿足精度要求------------

for t=1:MaxDT

for i=1:M

v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(gbest-x(i,:));

x(i,:)=x(i,:)+v(i,:);

if fitness(x(i,:),D)<p(i)

p(i)=fitness(x(i,:),D);

y(i,:)=x(i,:);

end

if p(i)<fitness(gbest,D)

gbest=y(i,:);

end

end

end

%------顯示計算結果

disp('*************************************************************')

disp('函數的全局最優位置為：')

Solution=gbest'

disp('最後得到的優化極值為：')

Result=fitness(gbest,D)

disp('*************************************************************')

9 實驗結果

如圖：