PSO
最近需要用上一點最優化相關的理論，特地去查了些PSO算法相關資料，在此記錄下學習筆記，附上程序代碼。基礎知識參考知乎大佬文章，寫得很棒！
傳送門
背景
起源：1995年，受到鳥群覓食行為的規律性啟發，James Kennedy和Russell Eberhart建立了一個簡化算法模型，經過多年改進最終形成了粒子群優化算法(Particle Swarm Optimization, PSO) ，也可稱為粒子群算法[1] 。
特點：粒子群算法具有收斂速度快、參數少、算法簡單易實現的優點（對高維度優化問題，比遺傳算法更快收斂於最優解）。
基本思想：粒子群算法的思想源於對鳥群覓食行為的研究，鳥群通過集體的信息共享使群體找到最優的目的地。
算法原理
| 鳥群覓食 | 粒子群算法 |
|---|---|
| 鳥 | 粒子 |
| 森林 | 求解空間 |
| 食物的量 | 目標函數值 |
| 每只鳥所處的位置 | 空間中的一個解(粒子位置) |
| 食物量最多的位置 | 全局最優解 |
-
每個粒子包含兩個屬性
- 速度：下一時刻粒子移動的方向與距離。
- 位置：所求問題粒子的當前解。
-
算法的6個重要參數
在\(D\)維空間中，有\(n\)個粒子。
-
第\(i(i=1,2,\cdots,n)\)個粒子的位置
\[\begin{align} X_{i}=\left [ x_{i1},x_{i2},\cdots,x_{iD}\right ] \end{align} \] -
第\(i(i=1,2,\cdots,n)\)個粒子的速度（粒子移動方向和大小）
\[\begin{align} V_{i}=\left [ v_{i1},v_{i2},\cdots,v_{iD}\right ] \end{align} \] -
第\(i(i=1,2,\cdots,n)\)個粒子的搜索到的最優位置（個體最優解）
\[\begin{align} P^{pbest}_{i}=\left [ p_{i1},p_{i2},\cdots,p_{iD}\right ] \end{align} \] -
群體搜索到的最優位置（群體最優解）為
\[\begin{align} P^{gbest}=\left [ p^{gbest}_{1},p^{gbest}_{2},\cdots,p^{gbest}_{D}\right ] \end{align} \] -
第 \(i(i=1,2,\cdots,n)\) 個粒子搜索到的最優位置的適應值（優化目標函數的值)為
\[\begin{align} f^p_i:\text{粒子歷史最優適應值} \end{align} \] -
群體搜索到的最優位置的適應值為
\[\begin{align} f^g:\text{群體歷史最優適應值} \end{align} \]
-
-
算法流程圖
-
速度更新公式
\[\begin{align} V^{k+1}_{i}=\omega \cdot V^{k}_{i}+c_1r_1(P^{pbest}_{i}-X^{k}_{i})+c_2r_2(P^{gbest}_{i}-X^{k}_{i}) \end{align} \]\(r_1,r_2\)為區間\([0,1]\)內的隨機數，\(\omega\)為權重因子，\(c_1,c_2\)分別為個體、群體學習因子
-
第一項：慣性部分
由慣性權重和粒子自身速度構成，表示粒子對先前自身運動狀態的信任。
-
第二項：認知部分
表示粒子本身的思考，即粒子自己經驗的部分，可理解為粒子當前位置與自身歷史最優位置之間的距離和方向。
-
第三項：社會部分
表示粒子之間的信息共享與合作，即來源於群體中其他優秀粒子的經驗，可理解為粒子當前位置與群體歷史最優位置之間的距離和方向
\[\text{更新速度方向}=\text{慣性方向 + 個體最優方向 + 群體最優方向} \] -
-
位置更新公式
\[\begin{align} X^{k+1}_{i}=X^{k}_{i}+V^{k}_{i} \end{align} \]
算法的參數解釋
-
粒子群規模\(n\)
一個正整數，推薦取值范圍：[20,1000]，簡單問題一般取20~40，較難或特定類別的問題可以取100~200。較小的種群規模容易陷入局部最優；較大的種群規模可以提高收斂性，更快找到全局最優解，但是相應地每次迭代的計算量也會增大；當種群規模增大至一定水平時，再增大將不再有顯著的作用。
-
粒子維度\(D\)
粒子搜索的空間維數即為自變量的個數。
-
迭代次數\(K\)
推薦取值范圍：[50,100]，典型取值：60、70、100
-
慣性權重\(\omega\)
慣性權重\(\omega\)描述粒子上一代速度對當前代速度的影響。\(\omega\)值較大，全局尋優能力強，局部尋優能力弱；反之，則局部尋優能力強。當問題空間較大時，為了在搜索速度和搜索精度之間達到平衡，通常做法是使算法在前期有較高的全局搜索能力以得到合適的種子，而在后期有較高的局部搜索能力以提高收斂精度。所以\(\omega\)不宜為一個固定的常數。
當\(\omega=1\) 時，退化成基本粒子群算法，當 \(\omega=0\) 時，失去對粒子本身經驗的思考。推薦取值范圍：$$[0.4,2]$$，典型取值：0.9、1.2、1.5、1.8
在解決實際優化問題時，往往希望先采用全局搜索，使搜索空間快速收斂於某一區域，然后采用局部精細搜索以獲得高精度的解。線性變化策略：隨着迭代次數的增加，慣性權重\(\omega\)不斷減小，從而使得粒子群算法在初期具有較強的全局收斂能力，在后期具有較強的局部收斂能力。
\[\omega= \omega_{max}-(\omega_{max}-\omega_{min})\frac{k}{k_{max}} \]\(\omega_{max},\omega_{min}\)為最大，最小慣性權重值，\(k,k_{max}\)為當前迭代次數與最大迭代次數。
-
學習因子\(c_1,c_2\)
也稱為加速系數或加速因子（這兩個稱呼更加形象地表達了這兩個系數的作用）
\(c_1=0\)：無私型粒子群算法，“只有社會，沒有自我”，迅速喪失群體多樣性，易陷入局部最優而無法跳出。
\(c_2=0\)：自我認知型粒子群算法，“只有自我，沒有社會”，完全沒有信息的社會共享，導致算法收斂速度緩慢 。
\(c_1,c_2\)都不為0，稱為完全型粒子群算法，更容易保持收斂速度和搜索效果的均衡，是較好的選擇。 推薦取值范圍：[0,4]；
典型取值：\(c_1=c_2=2、c_1=1.6,c_2=1.8、c_1=1.6,c_2=2\)針對不同的問題有不同的取值，一般通過試湊來調整這兩個值。
粒子群算法條件限制
位置限制
限制粒子搜索的空間，即自變量的取值范圍，對於無約束問題此處可以省略。
速度限制
為了平衡算法探索能力與開發能力，需要設定一個合理的速度范圍，限制粒子的最大速度。
\(V_m\)較大時，粒子飛行速度快，探索能力強，但粒子容易飛過最優解
\(V_m\)較小時，飛行速度慢，開發能力強，但是收斂速度慢，且容易陷入局部最優解
\(V_m\)一般設為每維變量變化范圍的10％~20％
優化停止准則
一般有兩種：
-
最大迭代步數
-
可接受的滿意解：上一次迭代后最優解的適應值與本次迭代后最優解的適應值之差小於某個值后停止優化
初始化：
如果粒子的初始化范圍選擇得好的話可以縮短優化的收斂時間，我們需要根據具體的問題進行分析。如果根據我們的經驗判斷出最優解一定在某個范圍內，則就在這個范圍內初始化粒子；如果無法確定，則以粒子的取值邊界作為初始化范圍。
MATLAB編程實現
實例效果
采用面相對象的編程思路，將完整的\(PSO\)算法以及數據成員封裝到類中，\(demo.m\)實例以及\(PSO.m\)類代碼見后續
\(Demo1D.m\)以一維函數\(f(x)=xsin(x) \cdot cos(2x)-2xsin(x)\)求\([0,20]\)范圍內最小值的PSO算法實例，迭代80時，結果達到收斂。
算法特性
- 支持多維空間變量
- 支持設定位置與速度許可范圍（for each dimension）
- 支持常系數與線性壓縮權重參數
對於多優化目標還需進一步改進。
\(demo1D.m\)
clc;clear all;
f= @(x)(x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x)); % 函數表達式
pso=PSO(f,50,1);
pso.SetRangeLim([-10,40],[-4,4]); %設置位置與速度范圍
%pso.SetOmega("const",0.2);
pso.SetOmega("linear",[0.5,2]); %設置權重Omega因子
pso.solve(200); %迭代求解100步
%% Plot
figure(111)
for i=1:1:pso.k_max
ezplot(f,[-10,0.01,40]);
hold on
x=pso.record_X_V{i,1}(:,1);
plot(x,f(x),'o'); title(num2str(i,"istep k = %d"))
hold off
pause(0.1)
end
\(demo2D.m\)
clc;clear all;
%%
f= @(x)x(:,1).^2+x(:,2).^2;
pso=PSO(f,50,2);
pso.SetRangeLim([-10,-10;10,10],[-2,-2;2,2]); %設置位置與速度范圍
pso.SetOmega("linear",[0.4,2]); %設置權重Omega因子
pso.Set_C_factor(1,1); %設置學習因子
pso.solve(200); %迭代求解100步
%% plot
figure(112)
[X,Y] = meshgrid(-10:0.5:10,-10:0.5:10);
Z=X.^2+Y.^2;
line=[0:1:10].^2;
%for i=1:1:pso.k_max
for i=160
contour(X,Y,Z,line,'ShowText','on');
hold on
x=pso.record_X_V{i,1}(:,1);
y=pso.record_X_V{i,1}(:,2);
plot(x,y,'o'); title(num2str(i,"istep k = %d"))
hold off
pause(0.1)
end
\(PSO.m \quad class\)
classdef PSO < handle
%PSO 算法類
properties
%% 目標函數
func;
%%
N=50; %種群個數，默認值50
D=1; %變量空間維數
k_max=100;%最大迭代步數,默認值為100
%%
omega; %權重因子/范圍
flag_omega="const";%權重形式
c1=1; %個體學習因子
c2=1; %群體學習因子
%% 粒子屬性
X; %位置
V; %速度
flag_lim=false;
lim_X; %位置限制范圍
lim_V; %速度限制范圍
%% 記錄
X_pbest; %粒子個體最佳位置
X_gbest; %粒子群體最佳位置
func_pbest; %每個個體的歷史最佳適應度
func_gbest; %種群歷史最佳適應度
%%
k; %當前迭代步
%%
record_func; %記錄位置與速度
record_X_V;
flag_plot=false;
end
methods
function self = PSO(func,N,D)
%PSO 構造函數
% fun:目標函數
% N :粒子數目
% D :空間維數
self.func=func;
self.N=N;
self.D=D;
self.X=zeros(self.N,self.D);
self.V=zeros(self.N,self.D);
self.X_pbest=zeros(self.N,self.D);
self.X_gbest=zeros(1,self.D);
self.func_pbest=ones(self.N,1)*1.0e10;
self.func_gbest=1.0e10;
fprintf("創建PSO Model完成\n");
end
function solve(self,k_max)
%求解
self.k=1;
self.record_func=zeros(k_max,self.D+1);%初始化記錄容器
self.record_X_V=cell(k_max,1);
self.k_max=k_max;
self.Initial(); %初始化
while(self.k<=self.k_max)
self.Update_p_g_best(); %更新群體
self.Update_V_X(); %更新速度、位置
self.k=self.k+1;
self.Report(10);
end
end
function SetRangeLim(self,limX,limV)
%設置空間位置與速度限制范圍
% limX(1,i):第i維數據最小值
% limX(2,i):第i維數據最大值
self.flag_lim=true;
self.lim_X=ones(2,self.D);
self.lim_X(1,:)=limX(1,:);
self.lim_X(2,:)=limX(2,:);
self.lim_V=ones(2,self.D);
self.lim_V(1,:)=limV(1,:);
self.lim_V(2,:)=limV(2,:);
end
function Set_C_factor(self,c1,c2)
%設置學習因子
% c1:個體學習因子
% c2:群體學習因子
self.c1=c1;
self.c2=c2;
end
function Initial(self)
%初始化粒子位置
if self.flag_lim
for d=1:1:self.D
for i=1:1:self.N
self.X(i,:)=self.lim_X(1,d)+(self.lim_X(2,d)-self.lim_X(1,d))*...
rand(1,self.D);
self.V(i,:)=self.lim_V(1,d)+(self.lim_V(2,d)-self.lim_V(1,d))*...
rand(1,self.D);
end
end
else
for i=1:1:self.N
self.X(i,:)= rand(1,self.D);
self.V(i,:)= rand(1,self.D);
end
end
fprintf("PSO Model初始化完成\n");
end
function SetOmega(self,flag,omega)
self.flag_omega=flag;
self.omega=omega;
end
function Update_V_X(self)
%更新速度矢量，與當前位置
% flag == const 常系數權重 omega為常值
% flag == linear 線性權重 omega為范圍
if self.flag_omega == "const"
omega_=self.omega;
elseif self.flag_omega == "linear"
omega_=self.omega(2)-(self.omega(2)-self.omega(1))*(self.k/self.k_max);
end
for i=1:1:self.N
self.V(i,:)=self.V(i,:)*omega_+...
self.c1*bsxfun(@times,self.X_pbest(i,:)-self.X(i,:),rand(1,self.D))+...
self.c2*bsxfun(@times,self.X_gbest-self.X(i,:),rand(1,self.D));
end
if self.flag_lim
% 限制速度范圍
for i=1:1:self.D
tmpV=self.V(:,i);
tmpV(tmpV<self.lim_V(1,i))=self.lim_V(1,i);
tmpV(tmpV>self.lim_V(2,i))=self.lim_V(2,i);
self.V(:,i)=tmpV;
end
end
self.X=self.X+self.V; %g更新位置
if self.flag_lim
% 限制位置范圍
for i=1:1:self.D
tmpX=self.X(:,i);
tmpX(tmpX<self.lim_X(1,i))=self.lim_X(1,i);
tmpX(tmpX>self.lim_X(2,i))=self.lim_X(2,i);
self.X(:,i)=tmpX;
end
end
self.record_X_V{self.k,1}=[self.X,self.V];
end
function Update_p_g_best(self)
% 更新X_pbest，X_gbest
func_x=self.func(self.X); % 個體當前適應度(函數值)
for i=1:1:self.N
if self.func_pbest(i) > func_x(i)
self.func_pbest(i)=func_x(i);
self.X_pbest(i,:)=self.X(i,:);
end
end
if self.func_gbest > min(self.func_pbest)
[self.func_gbest, nmin] = min(self.func_pbest); % 更新群體歷史最佳適應度
self.X_gbest = self.X_pbest(nmin, :); % 更新群體歷史最佳位置
end
self.record_func(self.k,:)=[self.X_gbest,self.func_gbest];
end
function Report(self,step)
if(mod(self.k,step)==0)
fprintf("\nCurrent istep k = %d\n",self.k);
fprintf("X = %f func= %f\n",self.X_gbest,self.func_gbest);
fprintf("Residual(func) = %f \n",abs(self.record_func(self.k-1,2)-self.record_func(self.k-2,2)));
end
end
end
end
參考文獻
[1] Kennedy J . Particle swarm optimization[J]. Proc. of 1995 IEEE Int. Conf. Neural Networks, (Perth, Australia), Nov. 27-Dec. 2011, 4(8):1942-1948 vol.4.