優化演算法—粒子群演算法(PSO)

04-07

優化演算法—粒子群演算法(PSO)

一、粒子群演算法的概述

粒子群演算法(PSO)屬於群智能演算法的一種，是通過模擬鳥群捕食行為設計的。假設區域里就只有一塊食物（即通常優化問題中所講的最優解），鳥群的任務是找到這個食物源。鳥群在整個搜尋的過程中，通過相互傳遞各自的信息，讓其他的鳥知道自己的位置，通過這樣的協作，來判斷自己找到的是不是最優解，同時也將最優解的信息傳遞給整個鳥群，最終，整個鳥群都能聚集在食物源周圍，即我們所說的找到了最優解，即問題收斂。

二、粒子群演算法的流程

粒子群演算法通過設計一種無質量的粒子來模擬鳥群中的鳥，粒子僅具有兩個屬性：速度和位置，速度代表移動的快慢，位置代表移動的方向。每個粒子在搜索空間中單獨的搜尋最優解，並將其記為當前個體極值，並將個體極值與整個粒子群里的其他粒子共享，找到最優的那個個體極值作為整個粒子群的當前全局最優解，粒子群中的所有粒子根據自己找到的當前個體極值和整個粒子群共享的當前全局最優解來調整自己的速度和位置。粒子群演算法的思想相對比較簡單，主要分為：1、初始化粒子群；2、評價粒子，即計算適應值；3、尋找個體極值；4、尋找全局最優解；5、修改粒子的速度和位置。下面是程序的流程圖：

（PSO流程）

下面我們具體解釋下流程圖裡面的每一個步驟：

1、初始化

首先，我們需要設置最大的速度區間，防止超出最大的區間。位置信息即為整個搜索空間，我們在速度區間和搜索空間上隨機初始化速度和位置。設置群體規模。

2、個體極值與全局最優解

個體極值為每個粒子找到的歷史上最優的位置信息，並從這些個體歷史最優解中找到一個全局最優解，並與歷史最優解比較，選出最佳的作為當前的歷史最優解。

3、更新速度和位置的公式

更新公式為：

其中，稱為慣性因子，稱為加速常數，一般取。表示區間上的隨機數。表示第個變數的個體極值的第d維。表示全局最優解的第d維。

4、終止條件

有兩種終止條件可以選擇，一是最大代數：；二是相鄰兩代之間的偏差在一個指定的範圍內即停止。我們在實驗中選擇第一種。

三、實驗

我們選擇的測試函數是：Griewank。其基本形式如下：

圖像為：

（Griewank函數圖像）

在實驗中我們選擇的維數是20；MATLAB程序代碼如下：

主程序：

[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片

c1=2;%學習因子

c2=2;%學習因子

Dimension=20;

Size=30;

Tmax=500;

Velocity_max=1200;%粒子最大速度

F_n=2;%測試函數名

Fun_Ub=600;%函數上下界

Fun_Lb=-600;

Position=zeros(Dimension,Size);%粒子位置

Velocity=zeros(Dimension,Size);%粒子速度

Vmax(1:Dimension)=Velocity_max;%粒子速度上下界

Vmin(1:Dimension)=-Velocity_max;

Xmax(1:Dimension)=Fun_Ub;%粒子位置上下界，即函數自變數的上下界

Xmin(1:Dimension)=Fun_Lb;

[Position,Velocity]=Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin);

Pbest_position=Position;%粒子的歷史最優位置，初始值為粒子的起始位置，存儲每個粒子的歷史最優位置

Gbest_position=zeros(Dimension,1);%全局最優的那個粒子所在位置，初始值認為是第1個粒子

for j=1:Size

Pos=Position(:,j);%取第j列，即第j個粒子的位置

fz(j)=Fitness_Function(Pos,F_n,Dimension);%計算第j個粒子的適應值

end 數據分析師培訓

[Gbest_Fitness,I]=min(fz);%求出所有適應值中最小的那個適應值，並獲得該粒子的位置

Gbest_position=Position(:,I);%取最小適應值的那個粒子的位置，即I列

for itrtn=1:Tmax

time(itrtn)=itrtn;

Weight=1;

r1=rand(1);

r2=rand(1);

for i=1:Size

Velocity(:,i)=Weight*Velocity(:,i)+c1*r1*(Pbest_position(:,i)-Position(:,i))+c2*r2*(Gbest_position-Position(:,i));

end

%限制速度邊界

for i=1:Size

for row=1:Dimension

if Velocity(row,i)>Vmax(row)

Veloctity(row,i)=Vmax(row);

elseif Velocity(row,i)<Vmin(row)

Veloctity(row,i)=Vmin(row);

else

end

Position=Position+Velocity;

%限制位置邊界

for i=1:Size

for row=1:Dimension

if Position(row,i)>Xmax(row)

Position(row,i)=Xmax(row);

elseif Position(row,i)<Xmin(row)

Position(row,i)=Xmin(row);

else

end

for j=1:Size

P_position=Position(:,j);%取一個粒子的位置

fitness_p(j)=Fitness_Function(P_position,F_n,Dimension);

if fitness_p(j)< fz(j) %粒子的適應值比運動之前的適應值要好，更新原來的適應值

Pbest_position(:,j)=Position(:,j);

fz(j)=fitness_p(j);

end

if fitness_p(j)<Gbest_Fitness

Gbest_Fitness=fitness_p(j);

end

[Gbest_Fitness_new,I]=min(fz);%更新後的所有粒子的適應值，取最小的那個，並求出其編號

Best_fitness(itrtn)=Gbest_Fitness_new; %記錄每一代的最好適應值

Gbest_position=Pbest_position(:,I);%最好適應值對應的個體所在位置

end

plot(time,Best_fitness);

xlabel(迭代的次數);ylabel(適應度值P_g);

初始化：

[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片

function [Position,Velocity] = Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin)

for i=1:Dimension

Position(i,:)=Xmin(i)+(Xmax(i)-Xmin(i))*rand(1,Size); % 產生合理範圍內的隨機位置，rand(1,Size)用於產生一行Size個隨機數

Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);

end

適應值計算：

[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片

function Fitness=Fitness_Function(Pos,F_n,Dimension)

switch F_n

case 1

Func_Sphere=Pos(:)*Pos(:);

Fitness=Func_Sphere;

case 2

res1=Pos(:)*Pos(:)/4000;

res2=1;

for row=1:Dimension

res2=res2*cos(Pos(row)/sqrt(row));

end

Func_Griewank=res1-res2+1;

Fitness=Func_Griewank;

end

最終的收斂曲線：

（收斂曲線）