數學建模：現代優化算法之粒子群算法

《數學建模算法與應用》在現代優化算法中沒有提及，依照組內進度要求，所以需要在此節補充有關“粒子群算法”的內容。

---

數學建模

• 數學建模：現代優化算法之粒子群算法
• 前言
• 一、算法簡介
• 二、算法步驟
• • 1.初始化
  • 2.更新規則
• 三、PSO算法的流程圖和偽代碼
• 四、PSO算法舉例
• 五、PSO算法的demo
• 總結
• 參考內容：

---

前言

先來看一段視頻：Matlab 粒子群算法PSO實例學習（有代碼和詳細注釋）

來自 “百科”： 粒子群優化算法(Particle Swarm optimization,PSO)又翻譯為粒子群算法、微粒群算法、或微粒群優化算法。是通過模擬鳥群覓食行為而發展起來的一種基于群體協作的隨機搜索算法。通常認為它是群集智能 (Swarm intelligence, SI) 的一種。它可以被納入多主體優化系統(Multiagent Optimization System, MAOS)。粒子群優化算法是由Eberhart博士和kennedy博士發明。

---

一、算法簡介

參考CSDN博主「lx青萍之末」的原創文章
最優化算法之粒子群算法（PSO）

二、算法步驟

1.初始化

PSO初始化為一群隨機粒子(隨機解)，然后通過迭代找到最優解，在每一次迭代中，粒子通過跟蹤兩個“極值”來更新自己。第一個就是粒子本身所找到的最優解，這個解叫做個體極值pBest，另一個極值是整個種群找到的最優解，這個極值是全局極值gBest。另外也可以不用整個種群而只是用其中一部分最優粒子的鄰居，那么在所有鄰居中的極值就是局部極值。

2.更新規則

公式(1)的第一部分稱為【記憶項】，表示上次速度大小和方向的影響；公式(1)的第二部分稱為【自身認知項】，是從當前點指向粒子自身最好點的一個矢量，表示粒子的動作來源于自己經驗的部分；公式(1)的第三部分稱為【群體認知項】，是一個從當前點指向種群最好點的矢量，反映了粒子間的協同合作和知識共享。粒子就是通過自己的經驗和同伴中最好的經驗來決定下一步的運動。以上面兩個公式為基礎，形成了PSO的標準形式。

三、PSO算法的流程圖和偽代碼

四、PSO算法舉例

【例一】

【解】

【例二】 本部分參考CSDN博主「nightmare_dimple」的原創文章
原文鏈接：https://blog.csdn.net/nightmare_dimple/article/details/74331679

對于函數 f=xsin(x)cos(2x)-2xsin(3x) ，求其在區間[0,20]上該函數的最大值。
【解】

clc;clear;close all; %% 初始化種群 f= @(x)x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x); % 函數表達式 figure(1);ezplot(f,[0,0.01,20]); N = 50; % 初始種群個數 d = 1; % 空間維數 ger = 100; % 最大迭代次數 limit = [0, 20]; % 設置位置參數限制 vlimit = [-1, 1]; % 設置速度限制 w = 0.8; % 慣性權重 c1 = 0.5; % 自我學習因子 c2 = 0.5; % 群體學習因子 for i = 1:dx = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, d);%初始種群的位置 end v = rand(N, d); % 初始種群的速度 xm = x; % 每個個體的歷史最佳位置 ym = zeros(1, d); % 種群的歷史最佳位置 fxm = zeros(N, 1); % 每個個體的歷史最佳適應度 fym = -inf; % 種群歷史最佳適應度 hold on plot(xm, f(xm), 'ro');title('初始狀態圖'); figure(2) %% 群體更新 iter = 1; record = zeros(ger, 1); % 記錄器 while iter <= gerfx = f(x) ; % 個體當前適應度 for i = 1:N if fxm(i) < fx(i)fxm(i) = fx(i); % 更新個體歷史最佳適應度xm(i,:) = x(i,:); % 更新個體歷史最佳位置end end if fym < max(fxm)[fym, nmax] = max(fxm); % 更新群體歷史最佳適應度ym = xm(nmax, :); % 更新群體歷史最佳位置endv = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(ym, N, 1) - x);% 速度更新% 邊界速度處理v(v > vlimit(2)) = vlimit(2);v(v < vlimit(1)) = vlimit(1);x = x + v;% 位置更新% 邊界位置處理x(x > limit(2)) = limit(2);x(x < limit(1)) = limit(1);record(iter) = fym;%最大值記錄 % x0 = 0 : 0.01 : 20; % plot(x0, f(x0), 'b-', x, f(x), 'ro');title('狀態位置變化') % pause(0.1)iter = iter+1; end figure(3);plot(record);title('收斂過程') x0 = 0 : 0.01 : 20; figure(4);plot(x0, f(x0), 'b-', x, f(x), 'ro');title('最終狀態位置') disp(['最大值：',num2str(fym)]); disp(['變量取值：',num2str(ym)]);

結果如下：

由上圖可以看出算法已成功找出了最優解，其最優解為18.3014，而其最大值為32.1462。
如果想看粒子群算法中粒子的搜索過程的可以將代碼中注釋掉的三行代碼放上去。效果是這樣的：

【例三】
版權聲明：本部分參考CSDN博主「簫韻」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/qq_41053605/article/details/88924287

（1）構建目標函數。這里使用只有兩個參數的函數，這樣便于畫出三維圖型。

function y=A11_01(x) y=x(1)^2+x(2)^2-x(1)*x(2)-10*x(1)-4*x(2)+60;

首先從直觀上看看這個函數：

x1=-15:1:15; x2=-15:1:15; [x1,x2]=meshgrid(x1,x2); y=x1.^2+x2.^2-x1.*x2-10.*x1-4.*x2+60; mesh(x1,x2,y);

圖 1 目標函數的三維網格圖

（2）整體實現代碼

clear clc%繪制原圖 圖1目標函數的三維網格圖 x1=-15:1:15; x2=-15:1:15; [x1,x2]=meshgrid(x1,x2); y=x1.^2+x2.^2-x1.*x2-10.*x1-4.*x2+60; mesh(x1,x2,y); hold on;%%預設參數 n=100; %粒子群的規模 d=2; %變量個數 c1=2; c2=2; w=0.9;%權重一般設為0.9 K=50; %迭代次數%%分布粒子的取值范圍及速度的取值范圍 x=-10+20*rand(n,d); %x在[-10,10]中取值 v=-5+10*rand(n,d); %v在[-5,5]中取值%%計算適應度 fit=zeros(n,1); for j=1:nfit(j)=A11_01(x(j,:)); end pbest=x; ind=find(min(fit)==fit); gbest=x(ind,:); h=scatter3(x(:,1),x(:,2),fit,'o'); %圖2 粒子的初始分布圖%%更新速度與位置 for i=1:Kfor m=1:nv(m,:)=w*v(m,:) + c1*rand*(pbest(m,:)-x(m,:)) + c2*rand*(gbest-x(m,:));%rand是[0,1]隨機數v(m,find(v(m,:)<-5))=-5;%這里發現速度小于-5時取-5v(m,find(v(m,:)>5))=5;%這里發現速度大于5時取5x(m,:)=x(m,:)+0.5*v(m,:);x(m,find(x(m,:)<-10))=-10;%這里發現位置小于-10時取-10x(m,find(x(m,:)>10))=10;%這里發現位置大于10時取10%重新計算適應度fit(m)=A11_01(x(m,:));if x(m,:)<A11_01(pbest(m,:))pbest(m,:)=x(m,:);endif A11_01(pbest(m,:))<A11_01(gbest)gbest=pbest(m,:);endendfitnessbest(i)=A11_01(gbest);pause(0.01); %為了直觀，每更新一次，暫停0.01秒h.XData=x(:,1);h.YData=x(:,2);h.ZData=fit; end % hold off; % plot(fitnessbest); % xlabel('迭代次數');

（3）粒子的初始分布圖

圖2 粒子的初始分布圖

（4）粒子群移動圖

圖 3 粒子群移動圖

（5）迭代過程

圖 4 迭代過程圖
從圖中可以看出，迭代次數基本上在13,14的時候就達到最優，即目標函數取得最小值為8.

五、PSO算法的demo

#include <iostream> #include <vector> #include <cmath> #include <map> #include <algorithm> #include <random> #include <ctime> #include <Eigen/Dense> using namespace Eigen; using namespace std;const int dim = 1;//維數 const int p_num = 10;//粒子數量 const int iters = 100;//迭代次數 const int inf = 999999;//極大值 const double pi = 3.1415; //定義粒子的位置和速度的范圍 const double v_max = 4; const double v_min = -2; const double pos_max = 2; const double pos_min = -1; //定義位置向量和速度向量 vector<double> pos; vector<double> spd; //定義粒子的歷史最優位置和全局最優位置 vector<double> p_best; double g_best; //使用eigen庫定義函數值矩陣和位置矩陣 Matrix<double, iters, p_num> f_test; Matrix<double, iters, p_num> pos_mat;//定義適應度函數 double fun_test(double x) {double res = x * x + 1;return res; }//初始化粒子群的位置和速度 void init() {//矩陣中所有元素初始化為極大值f_test.fill(inf);pos_mat.fill(inf);//生成范圍隨機數static std::mt19937 rng;static std::uniform_real_distribution<double> distribution1(-1, 2);static std::uniform_real_distribution<double> distribution2(-2, 4);for (int i = 0; i < p_num; ++i){pos.push_back(distribution1(rng));spd.push_back(distribution2(rng));}vector<double> vec;for (int i = 0; i < p_num; ++i){auto temp = fun_test(pos[i]);//計算函數值//初始化函數值矩陣和位置矩陣f_test(0, i) = temp;pos_mat(0, i) = pos[i];p_best.push_back(pos[i]);//初始化粒子的歷史最優位置}std::ptrdiff_t minRow, minCol;f_test.row(0).minCoeff(&minRow, &minCol);//返回函數值矩陣第一行中極小值對應的位置g_best = pos_mat(minRow, minCol);//初始化全局最優位置 }void PSO() {static std::mt19937 rng;static std::uniform_real_distribution<double> distribution(0, 1);for (int step = 1; step < iters; ++step){for (int i = 0; i < p_num; ++i){//更新速度向量和位置向量spd[i] = 0.5 * spd[i] + 2 * distribution(rng) * (p_best[i] - pos[i]) +2 * distribution(rng) * (g_best - pos[i]);pos[i] = pos[i] + spd[i];//如果越界則取邊界值if (spd[i] < -2 || spd[i] > 4)spd[i] = 4;if (pos[i] < -1 || pos[i] > 2)pos[i] = -1;//更新位置矩陣pos_mat(step, i) = pos[i];}//更新函數值矩陣for (int i = 0; i < p_num; ++i){auto temp = fun_test(pos[i]);f_test(step, i) = temp;}for (int i = 0; i < p_num; ++i){MatrixXd temp_test;temp_test = f_test.col(i);//取函數值矩陣的每一列std::ptrdiff_t minRow, minCol;temp_test.minCoeff(&minRow, &minCol);//獲取每一列的極小值對應的位置p_best[i] = pos_mat(minRow, i);//獲取每一列的極小值，即每個粒子的歷史最優位置}g_best = *min_element(p_best.begin(), p_best.end());//獲取全局最優位置}cout << fun_test(g_best); }int main() {init();PSO();system("pause");return 0; }

總結

參考內容：

版權聲明：本文部分為CSDN博主「lx青萍之末」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/daaikuaichuan/article/details/81382794

————————————————
版權聲明：本文部分為CSDN博主「nightmare_dimple」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/nightmare_dimple/article/details/74331679

總結

以上是生活随笔為你收集整理的数学建模：现代优化算法之粒子群算法的全部內容，希望文章能夠幫你解決所遇到的問題。

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