一、實驗題目

計算如下二元函數的最小值：

（其中自變量x、y的范圍均為[-50, 50]）

• 用matlab代碼實現。
• 代碼必須能一鍵運行。
• 最后輸出x,y和z的最優值及收斂到最優值所需的迭代步數。
• 算法關鍵參數需要注釋清楚（如權重因子、學習因子、速度變化范圍等）。

二、算法原理

粒子群優化算法求最優解

1、問題描述

D維空間中，有N個粒子；

??? 粒子i位置：xi=(xi1,xi2,…xiD)，將xi代入適應函數f(xi)求適應值；

??? 粒子i速度：vi=(vi1,vi2,…viD)

??? 粒子i個體經歷過的最好位置：pbesti=(pi1,pi2,…piD)

??? 種群所經歷過的最好位置：gbest=(g1,g2,…gD)

通常，在第d（1≤d≤D）維的位置變化范圍限定在內，速度變化范圍限定在?內（即在迭代中若?超出了邊界值，則該維的速度或位置被限制為該維最大速度或邊界位置）

2、算法公式

粒子i的第d維速度更新公式：

???

粒子i的第d維位置更新公式：

???—第k次迭代粒子i飛行速度矢量的第d維分量

—第k次迭代粒子i位置矢量的第d維分量

??? c1,c2—加速度常數，調節學習最大步長

??? r1,r2—兩個隨機函數，取值范圍[0，1]，以增加搜索隨機性

??? w —慣性權重，非負數，調節對解空間的搜索范圍

粒子速度更新公式包含三部分：

?? 第一部分為粒子先前的速度

?? 第二部分為“認知”部分，表示粒子本身的思考，可理解為粒子i當前位置與自己最好位置之間的距離。

?? 第三部分為“社會”部分，表示粒子間的信息共享與合作，可理解為粒子i當前位置與群體最好位置之間的距離。

3、算法流程

（1）Initial：

初始化粒子群體（群體規模為n），包括隨機位置和速度。

設定參數運動范圍，設定學習因子c1、c2，最大進化代數G，kg表示當前的進化代數。在一個D維參數的搜索解空間中，粒子組成的種群規模大小為Size，每個粒子代表解空間的一個候選解，其中第i（1≤i≤Size）個粒子在整個解空間的位置表示為Xi，速度表示為Vi。第i個粒子從初始到當前迭代次數搜索產生的最優解，個體極值Pi，整個種群目前的最優解為BestS。隨機產生Size個粒子，隨機產生初始種群的位置矩陣和速度矩陣。

（2）Evaluation：

?根據fitness function ，評價每個粒子的適應度。

個體評價（適應度評價）：將各個粒子初始位置作為個體極值，計算群體中各個粒子的初始適應值f（Xi），并求出種群最優位置。

（3）Update the speed and position:

更新粒子的速度和位置，產生新種群，并對粒子的速度和位置進行越界檢查，為避免算法陷入局部最優解，加入一個局部自適應變異算子進行調整。

其中，kg=1，2，…，G，i=1，2，…，Size，r1和r2為0到1的隨機數，c1為局部學習因子，c2為全局學習因子，一般取c2大些。

(4)Find the Pbest：

?? 對每個粒子，將其當前適應值與其個體歷史最佳位置（pbest）對應的適應值做比較，如果當前的適應值更高，則將用當前位置更新歷史最佳位置pbest。

(5）Find the Gbest：

對每個粒子，將其當前適應值與全局最佳位置（gbest）對應的適應值做比較，如果當前的適應值更高，則將用當前粒子的位置更新全局最佳位置gbest。

（6）Update the Velocity：

?根據公式更新每個粒子的速度與位置。如未滿足結束條件，則返回步驟2。通常算法達到最大迭代次數???? 或者最佳適應度值的增量小于某個給定的閾值時算法停止。

三、實驗代碼與結果

1、實驗代碼

function main%由于原代碼求取的最大值，最小值的話直接將目標函數取成相反數即可 %很顯然，結果只能無限逼近 3與5 跟算法原理有關(類似遺傳算法）clear all; close all; %（1）初始化粒子群算法參數 min=-50;max=50;%粒子位置范圍 Vmax=25;Vmin=-25;%粒子運動速度范圍 c1=1.3;c2=1.7; %學習因子[0，4]wmin=0.20;wmax=0.90;%慣性權重 G=400; % 最大迭代次數 Size=100; %初始化群體個體數目for i=1:Gw(i)=wmax-((wmax-wmin)/G)*i; %隨著優化進行，應降低自身權重 end for i=1:Sizefor j=1:2x(i,j)=min+(max-min)*rand(1); %隨機初始化位置v(i,j)=Vmin +(Vmax-Vmin)*rand(1); %隨機初始化速度end end%（2）計算各個粒子的適應度，并初始化Pi、plocal和最優個體BestS for i=1:Sizep(i)=func(x(i,:));y(i,:)=x(i,:);if i==1plocal(i,:)=evaluate_localbest(x(Size,:),x(i,:),x(i+1,:));elseif i==Sizeplocal(i,:)=evaluate_localbest(x(i-1,:),x(i,:),x(1,:));else plocal(i,:)=evaluate_localbest(x(i-1,:),x(i,:),x(i+1,:));end endBestS=x(1,:);%初始化最優個體BestS for i=2:Sizeif func(x(i,:))>func(BestS)BestS=x(i,:);end end%（3）進入主循環 for kg=1:Gfor i=1:SizeM=1;if M==1v(i,:)=w(kg)*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(plocal(i,:)-x(i,:));%局部尋優：加權，實現速度的更新elseif M==2v(i,:)=w(kg)*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(BestS-x(i,:)); %全局尋優：加權，實現速度的更新endfor j=1:2 %檢查速度是否越界if v(i,j)<Vminv(i,j)=Vmin;elseif x(i,j)>Vmaxv(i,j)=Vmax;endendx(i,:)=x(i,:)+v(i,:)*1; %實現位置的更新for j=1:2 %檢查位置是否越界if x(i,j)<minx(i,j)=min;elseif x(i,j)>maxx(i,j)=max;endend %自適應變異,避免粒子群算法陷入局部最優if rand>0.60k=ceil(2*rand);x(i,k)=min+(max-min)*rand(1);end %（4）判斷和更新if i==1plocal(i,:)=evaluate_localbest(x(Size,:),x(i,:),x(i+1,:));elseif i==Sizeplocal(i,:)=evaluate_localbest(x(i-1,:),x(i,:),x(1,:));else plocal(i,:)=evaluate_localbest(x(i-1,:),x(i,:),x(i+1,:));endif func(x(i,:))>p(i) %判斷當此時的位置是否為最優的情況，當不滿足時繼續更新p(i)=func(x(i,:));y(i,:)=x(i,:);endif p(i)>func(BestS)BestS=y(i,:);endend Best_value(kg)=func(BestS); end figure(1); kg=1:G; plot(kg,-Best_value,'r','linewidth',2); xlabel('generations');ylabel('Fitness function');% display('Best Sample=');disp(BestS); display('X=');disp(BestS(1)); display('Y=');disp(BestS(2)); % display('Biggest value=');disp(Best_value(G)); display('Z=');disp(Best_value(G)); for i=1:Gif Best_value(i)==Best_value(end)disp('收斂到最優值需要的迭代步數');disp(i);break;end endendfunction f = func(x) f=-(0.5*(x(1)-3)^2+0.2*(x(2)-5)^2-0.1); endfunction f =evaluate_localbest(x1,x2,x3)%求解粒子環形鄰域中的局部最優個體 K0=[x1;x2;x3]; K1=[func(x1),func(x2),func(x3)]; [maxvalue index]=max(K1); plocalbest=K0(index,:); f=plocalbest; end

2、實驗結果

X=3.0022Y=5.0068Z=0.1000收斂到最優值需要的迭代步數379

PS:該代碼有智能控制劉金錕第四版更改，更加詳細的內容請查看原作。

總結

以上是生活随笔為你收集整理的粒子群优化算法求解函数最值的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。