MATLAB-粒子群算法-最详细原理详解及代码分析

量子计算9 · 发表于 2022-3-27 12:17

  本文主要分为四个部分，对MSTLAB中的粒子群优化算法进行最详细的解释与代码展示。
  这可能是知乎上最全面的相关文章，码字不易，请大家动动小手，点个赞哈哈哈~
目录

一、概念篇
二、思想及公式解析篇
三、例题解析篇
四、代码篇

一、概念篇

中文名：粒子群优化算法
英文名：Particle Swarm Optimization（PSO）
  一群鸟在寻找食物时，根据个体探索与群体信息共享找到食物最多的地方——这便是粒子群优化算法的基本思想来源。
  粒子群优化算法是一种新的群智能算法，是在自然界中鸟群觅食行为的启发下提出的。粒子群优化算法的基本思想是通过群体中个体之间的协作与信息共享来寻找最优解。由于该算法结构简单、需要调节的参数少等优点，现已被应用于参数优化、神经网络等领域。

二、思想及公式解析篇

1、思想解析

  粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度与位置。其中，速度代表粒子移动的快慢，位置代表粒子移动的方向。个体层面，每个粒子在搜索空间中单独的搜寻最优解，并不断记录与更新最优解，即pbest。全局层面，所有粒子的最优解与其他粒子互相共享，所有粒子中最优的个体最优解记为全局最优解，该解也会不断记录与更新，即gbest。
种群中的所有粒子根据自己找到的个体最优解与整个种群的全局最优解来调整自己的速度与位置。
  以下分别是一维与二维的粒子群优化算法更新过程。

2、公式解析

  粒子群优化算法的公式是对这句话的数学表达：
种群中的所有粒子根据自己找到的个体最优解与整个种群的全局最优解来调整自己的速度与位置。

$v_{i}=v_{i}+c_{1}\times rand()\times (pbest_{i}-x_{i})+c_{2}\times rand()\times (gbest_{i}-x_{i})$       (1)

$x_{i}=x_{i}+v_{i}$                                                                                                                (2)
  公式(1)、(2)分别是更新速度与位置的公式。其中，i=1,2,...,N, N是粒子个数。
$v_{i}$ 是粒子的速度， $x_{i}$ 是粒子的位置， $rand()$ 是介于(0, 1)之间的随机数， $c_{1}$ 与 $c_{2}$ 是学习因子。
  在公式(1)中，第一部分称为记忆项，表示上一次的速度对本次更新的影响；第二部分称为自身认知项，表示从当前点指向粒子个体历史最优解的一个矢量，即粒子自身经验对本次更新的影响；第三部分称为群体认知项，表示从当前点指向全局历史最优解的一个矢量，即种群经验对本次更新的影响。
  公式(1)、(2)统称为粒子群优化算法的标准形式。

$v_{i}=\omega \times v_{i}+c_{1}\times rand()\times (pbest_{i}-x_{i})+c_{2}\times rand()\times (gbest_{i}-x_{i})$ (3)
  公式(3)中的 $\omega$ 称为惯性因子。
惯性因子较大，全局寻优能力强，局部寻优能力弱；
惯性因子较小，全局寻优能力弱，局部寻优能力强。
  由于粒子群优化算法初始化为一群遍布在取值范围内的随机粒子，分布范围广，因此更新初期应使惯性因子取值较大，全局寻优能力强。在更新后期，粒子逐渐收敛到一个或多个最优解，仅需在有限解中做出判断，因此更新后期应使惯性因子取值较小，局部寻优能力强。因此介绍目前使用较多的线性递减权值策略：

$\omega = (\omega_{ini}-\omega_{end})(G_{k}-g)/G_{k}+\omega_{end}$                                                                   (4)
  在公式(4)中， $\omega_{ini}$ 表示初始值， $\omega_{end}$ 表示最终更新值。
  公式(3)、(2)称为标准粒子群优化算法。

三、例题解析篇

  已知函数 $f(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}$ ， $-10\leq x_{1},x_{2}\leq10$ 求f的最小值。

四、代码篇

1、清空运行环境

clc;
clear;
close all;2、参数初始化

C1 = 2;             % 加速系数
C2 = 2;
PopG = 100;          % 迭代次数
PopSize = 20;       % 种群规模
D = 2;             % 维度
Popmax = 1;          % 种群上下边界值
Popmin = 0;          %
Vmax = 0.15 * Popmax; % 速度限定，通常在搜索范围的10%~20%均可，这里取15%
Vmin = -Vmax;       % Vmax的相反数
Wini = 0.9;          % W为惯性因子omega，Wini为初始值
Wend = 0.4;          % Wend为终止值
K = 51;             % 适应度曲线K次叠加平均ps：由于粒子群优化算法属于一种随机算法，每一次运行程序的结果、适应度曲线会有所区别，为降低这种随机性，取运行51次程序的适应度曲线的平均结果作为最终的适应度曲线。
3、第零代种群粒子初始化

for j = 1: PopSize                   % 将popsize个粒子依次初始化
Pop(j, :) = Popmax * rand(1, D);  % 第i个粒子初始位置
V(j, :) = Vmax * rand(1, D);    % 第i个粒子初始速度
[fitness(j), ~] = Ackley(Pop(j, :)); % 第i个粒子适应度
endps：1、Ackley函数为适应度评价函数，其图像如封面所示。
2、所谓第零代，可以理解为进入循环的第一代做准备。
4、第零代种群寻优

[Bestfitness, Bestindex] = min(fitness); % 第一代种群最佳适应度
PBest = Pop;                            % 个体历史最佳位置
GBest = Pop(Bestindex, :);             % 全局历史最佳位置
fitnesspbest = fitness;                % 个体历史最佳适应度
fitnessgbest = Bestfitness;             % 全局历史最佳适应度5、更新

for k = 1: K
fprintf(&#39;第%d次运行程序\n&#39;, k);
for i = 1: PopG       % 每一代循环
      for j = 1: PopSize % 每个粒子循环
         %%% 速度更新
         W = (Wini - Wend) * (PopG - i) / PopG + Wend; % 线性递减惯性因子
         V(j,:) = W*V(j,:) + C1*rand*(PBest(j,:)-Pop(j,:)) + C2*rand*(GBest-Pop(j,:));
         %%% 越界处理
         V(j,find(V(j,:) > Vmax)) = Vmax;
         V(j,find(V(j,:) < Vmin)) = Vmin;
         %%% 位置更新
         Pop(j,:) = Pop(j,:) + V(j,:);
         %%% 越界处理
         Pop(j,find(Pop(j,:) > Popmax)) = Popmax;
         Pop(j,find(Pop(j,:) < Popmin)) = Popmin;
         %%% 计算适应度值
         [fitness(j), ~] = ActionPredict(Pop(j, :));
         %%% 个体历史最优更新
         if fitness(j) < fitnesspbest(j)
            PBest(j, :) = Pop(j, :);    % 第i代种群第j个粒子历史最佳位置更新
            fitnesspbest(j) = fitness(j); % 第i代种群第j个粒子历史最佳适应度更新
         end
      end
      %%% 种群历史最优更新
      Popgbest = find(fitnesspbest < fitnessgbest);
      % 第i代种群所有粒子最佳历史适应度与全局历史最佳适应度比较
      if isempty(Popgbest) == 0
         if length(Popgbest) ~= 1
            % 如果多个粒子历史最佳适应度都优于全局历史最佳适应度，随机挑选
            Popgbest = randsrc(1, 1, Popgbest);
         end
         GBest = Pop(Popgbest, :);          % 第i代种群全局历史最佳位置更新
         fitnessgbest = fitnesspbest(Popgbest);% 第i代种群全局历史最佳适应度更新
      end
      result(k, i) = fitnessgbest;  % 储存历代全局历史最佳适应度值
end
%%% 历代最佳位置粒子不做记录，以最后一次迭代作为历代最佳位置粒子，填入GBEST
GBEST(k, :) = GBest;
end
[~, INDEX] = min(result(:, PopG));  % 比较所有运行程序中，最后一次迭代的fitnessgbest
FITNESS_RESULT = mean(result);    % 适应度曲线取K次叠加平均，避免随机性
GBEST_RESULT = GBEST(INDEX, :);    % 全局最佳粒子位置取K次中最优6、绘图

semilogy(FITNESS_RESULT);
grid on;
%%% 坐标含义说明
title(&#39;适应度曲线 &#39;);
xlabel(&#39;进化代数&#39;)
ylabel(&#39;适应度&#39;);如有错误，还请大家多多指正~
最后，希望大家在阅读完之后，能够点个赞，谢谢大家鸭~

参考出处
文章的第二部分参考自~青萍之末~的CSDN文章《最优化算法之粒子群算法(PSO)》
文章的第四部分在koko可可：Matlab标准粒子群优化算法代码基础上进行修改与添加

		自动登录	找回密码
密码			立即注册

MATLAB-粒子群算法-最详细原理详解及代码分析

本帖子中包含更多资源