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

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

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

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

圖像為：

（Griewank函數(shù)圖像）
在實(shí)驗(yàn)中我們選擇的維數(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;%測(cè)試函數(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)位置，初始值為粒子的起始位置，存儲(chǔ)每個(gè)粒子的歷史最優(yōu)位置  
Gbest_position=zeros(Dimension,1);%全局最優(yōu)的那個(gè)粒子所在位置，初始值認(rèn)為是第1個(gè)粒子  
 
for j=1:Size  
    Pos=Position(:,j);%取第j列，即第j個(gè)粒子的位置  
    fz(j)=Fitness_Function(Pos,F_n,Dimension);%計(jì)算第j個(gè)粒子的適應(yīng)值  
end  數(shù)據(jù)分析師培訓(xùn)
[Gbest_Fitness,I]=min(fz);%求出所有適應(yīng)值中最小的那個(gè)適應(yīng)值，并獲得該粒子的位置  
Gbest_position=Position(:,I);%取最小適應(yīng)值的那個(gè)粒子的位置，即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)';%取一個(gè)粒子的位置  
     fitness_p(j)=Fitness_Function(P_position,F_n,Dimension);  
     if fitness_p(j)< fz(j) %粒子的適應(yīng)值比運(yùn)動(dòng)之前的適應(yīng)值要好，更新原來(lái)的適應(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)值，取最小的那個(gè)，并求出其編號(hào)  
   Best_fitness(itrtn)=Gbest_Fitness_new; %記錄每一代的最好適應(yīng)值  
   Gbest_position=Pbest_position(:,I);%最好適應(yīng)值對(duì)應(yīng)的個(gè)體所在位置  
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個(gè)隨機(jī)數(shù)  
      Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);  
  end  
end  

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

最終的收斂曲線：

（收斂曲線）