優(yōu)化算法—粒子群算法(PSO)

一、粒子群算法的概述
    粒子群算法(PSO)屬于群智能算法的一種，是通過模擬鳥群捕食行為設(shè)計的。假設(shè)區(qū)域里就只有一塊食物（即通常優(yōu)化問題中所講的最優(yōu)解），鳥群的任務(wù)是找到這個食物源。鳥群在整個搜尋的過程中，通過相互傳遞各自的信息，讓其他的鳥知道自己的位置，通過這樣的協(xié)作，來判斷自己找到的是不是最優(yōu)解，同時也將最優(yōu)解的信息傳遞給整個鳥群，最終，整個鳥群都能聚集在食物源周圍，即我們所說的找到了最優(yōu)解，即問題收斂。
二、粒子群算法的流程
    粒子群算法通過設(shè)計一種無質(zhì)量的粒子來模擬鳥群中的鳥，粒子僅具有兩個屬性：速度和位置，速度代表移動的快慢，位置代表移動的方向。每個粒子在搜索空間中單獨的搜尋最優(yōu)解，并將其記為當(dāng)前個體極值，并將個體極值與整個粒子群里的其他粒子共享，找到最優(yōu)的那個個體極值作為整個粒子群的當(dāng)前全局最優(yōu)解，粒子群中的所有粒子根據(jù)自己找到的當(dāng)前個體極值和整個粒子群共享的當(dāng)前全局最優(yōu)解來調(diào)整自己的速度和位置。粒子群算法的思想相對比較簡單，主要分為：1、初始化粒子群；2、評價粒子，即計算適應(yīng)值；3、尋找個體極值；4、尋找全局最優(yōu)解；5、修改粒子的速度和位置。下面是程序的流程圖：

（PSO流程）
下面我們具體解釋下流程圖里面的每一個步驟：
1、初始化
   首先，我們需要設(shè)置最大的速度區(qū)間，防止超出最大的區(qū)間。位置信息即為整個搜索空間，我們在速度區(qū)間和搜索空間上隨機(jī)初始化速度和位置。設(shè)置群體規(guī)模。
2、個體極值與全局最優(yōu)解
   個體極值為每個粒子找到的歷史上最優(yōu)的位置信息，并從這些個體歷史最優(yōu)解中找到一個全局最優(yōu)解，并與歷史最優(yōu)解比較，選出最佳的作為當(dāng)前的歷史最優(yōu)解。
3、更新速度和位置的公式
   更新公式為：

其中，稱為慣性因子，稱為加速常數(shù)，一般取。表示區(qū)間上的隨機(jī)數(shù)。表示第個變量的個體極值的第d維。表示全局最優(yōu)解的第d維。
4、終止條件
有兩種終止條件可以選擇，一是最大代數(shù)：；二是相鄰兩代之間的偏差在一個指定的范圍內(nèi)即停止。我們在實驗中選擇第一種。
三、實驗
    我們選擇的測試函數(shù)是：Griewank。其基本形式如下：

圖像為：

（Griewank函數(shù)圖像）
在實驗中我們選擇的維數(shù)是20；MATLAB程序代碼如下：
主程序：
[plain] view plain copy 在CODE上查看代碼片派生到我的代碼片
c1=2;%學(xué)習(xí)因子  
c2=2;%學(xué)習(xí)因子  
Dimension=20;  
Size=30;  
Tmax=500;  
Velocity_max=1200;%粒子最大速度  
 
F_n=2;%測試函數(shù)名  
 
Fun_Ub=600;%函數(shù)上下界  
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;%粒子位置上下界，即函數(shù)自變量的上下界  
Xmin(1:Dimension)=Fun_Lb;  
[Position,Velocity]=Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin);  
 
Pbest_position=Position;%粒子的歷史最優(yōu)位置，初始值為粒子的起始位置，存儲每個粒子的歷史最優(yōu)位置  
Gbest_position=zeros(Dimension,1);%全局最優(yōu)的那個粒子所在位置，初始值認(rèn)為是第1個粒子  
 
for j=1:Size  
    Pos=Position(:,j);%取第j列，即第j個粒子的位置  
    fz(j)=Fitness_Function(Pos,F_n,Dimension);%計算第j個粒子的適應(yīng)值  
end  數(shù)據(jù)分析師培訓(xùn)
[Gbest_Fitness,I]=min(fz);%求出所有適應(yīng)值中最小的那個適應(yīng)值，并獲得該粒子的位置  
Gbest_position=Position(:,I);%取最小適應(yīng)值的那個粒子的位置，即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  
    end  
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  
    end  
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) %粒子的適應(yīng)值比運(yùn)動之前的適應(yīng)值要好，更新原來的適應(yīng)值  
         Pbest_position(:,j)=Position(:,j);  
         fz(j)=fitness_p(j);  
     end  
     if fitness_p(j)<Gbest_Fitness  
         Gbest_Fitness=fitness_p(j);  
     end  
  end  
  [Gbest_Fitness_new,I]=min(fz);%更新后的所有粒子的適應(yīng)值，取最小的那個，并求出其編號  
   Best_fitness(itrtn)=Gbest_Fitness_new; %記錄每一代的最好適應(yīng)值  
   Gbest_position=Pbest_position(:,I);%最好適應(yīng)值對應(yīng)的個體所在位置  
end  
plot(time,Best_fitness);  
xlabel('迭代的次數(shù)');ylabel('適應(yīng)度值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); % 產(chǎn)生合理范圍內(nèi)的隨機(jī)位置，rand(1,Size)用于產(chǎn)生一行Size個隨機(jī)數(shù)  
      Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);  
  end  
end  

適應(yīng)值計算：
[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  

最終的收斂曲線：

（收斂曲線）