最新第6弹！带你轻松理解MATLAB中的粒子群算法及其改良版本，持续分享深度解析

一、经典粒子群PSO算法

1 思想来源

粒子群优化(Particle Swarm Optimization,PSO) 作为进化计算的一个分支，是由 Eberhart 和 Kennedy 于 1995 提出的一种搜索算法。该算法在提出时是为了优化非线性函数，亦即求解非线性函数在某个求解范围内的最优解。

粒子群优化是一种模拟自然界生物的活动以及群体智能的随机搜索算法，该算法吸取了人工生命(Artificial Life)、鸟群觅食(Birds Flocking)、鱼群学习(Fish Schooling) 和 群理论(Swarm Theory) 的思想，通过对以上群体智能行为的抽象模仿，从而缔造出粒子群优化 。从这我们也可以看出，粒子群优化，其本质就是一种 群体智能优化 ，该优化算法就是通过模拟生物群体在觅食或者躲避天敌时展现的群体智能行为而产生的。

2 基本原理

2.1 群体智能(Swarm Intelligence)

在自然界中，鸟类是通过什么养的机制发现猎物的呢？或许，我们很少见到鸟群协同寻找食物，但是在动物世界中，我们经常可以看到庞大的鱼群聚集在一起躲避被猎杀。以及狼群和猎狗等等，它们是如何发现猎物的呢？

当我们不断剖析其中原因，会发现，狼群以及猎狗在发现猎物之后会通过声音向同伴传递信息，使得同伴不断向猎物所在地不断聚集。而对于其他的亦是如此，虽然我们无法理解鸟群以及鱼群是如何进行协同配合来躲避天敌或是寻找食物，但是在它们之间一定是存在某种可供交流的机制，使得各个个体之间的信息、经验能够互享，随着群体中经验的累加，就会展现出一种超越个体的 群体智能 ，这也是我们为什么说人多力量大，其实人类也是具有群体智能的，并且我们老祖宗很早之前就已经得出这个结论。

2.2 算法原理

在上文我们已经分析了群体智能的基本原理，现在最关键的，便是如何将其要素进行抽象，最后变成能够进行编程的规则要求。

在之前的场景中已经提出，一群分散寻找食物的狼或者猎狗，我们以狼群为例，一开始并不知道猎物的所在，但是会有一个间接的机制会让狼群知道猎物的所在，比如猎物留下的足迹以及气味等。于是狼就会随着这些痕迹不断分散搜索，寻找气味最浓的位置方向，并通过声音或其他信息交流的方式通知其他狼，以此不断调整整个群体的位置，保证不断向猎物逼近，以至于最后发现猎物并捕食。这样，狼群中的每个个体都有了一个指导方向，总体的方向是猎物的方向，在这个过程中，所有个体不断地配合，根据自身的经验以及整个群体的经验调整自己的速度以及跟踪方向。

在很多的文献当中，更多的是以鸟群为例，但是在我个人理解看来，对于鸟群觅食，我们见的少之又少，对于理解这些概念或许有些不太适合。而对于狼群或者猎狗，却是再合适不过。此外，在上面的小标题中，我写的是 群体智能 ，而不是粒子群智能，这也是我对于该算法的考虑。在我看来，它的重心，亦即算法核心在于群体协同所产生超乎个体的智能，所以我更愿意将它称之为 群体智能 ，至于该优化算法的名字中的 粒子 一词，其实是将群体中的个体抽象为一个点，一个带有位置速度等性质的点，我们称之为 粒子。

下面通过一个表格以狼群捕食和粒子群优化算法的基本定义进行对照

在智能算法领域，却有人提出了狼群算法，所以有必要进行解释，在该文只是以狼群的行为帮助大家更好的理解粒子群优化算法，所以在这里进行注明。

3 粒子群优化算法基本流程

3.1 基本流程

粒子群优化(PSO)要求每个个体(粒子)在进化的过程中维护两个向量，分别是速度向量v i = [ v i 1 , v i 2 , ⋯   , v i D ] v_i=[v_i^1,v_i^2,\dotsb,v_i^D]v i =[v i1 ,v i2 ,⋯,v iD ]和位置向量x i = [ x i 1 , x i 2 , ⋯   , x i D ] x_i=[x_i^1,x_i^2,\dotsb,x_i^D]x i =[x 1 ,x i2 ,⋯,x iD ]，其中i ii是指粒子的编号，D DD是求解问题空间的维数(通俗一点就是未知数的个数)。粒子的速度决定了自身运动的方向和速率，而其所处的位置则是解空间中的一个可行解，然后通过特定的评价方式以及评价函数评价该解的优劣。

3.1.1 两个主要变量

自身历史最优(p B e s t pBestpBest) 在该算法中，每个粒子还需要维护自身的历史最优位置

pBest

也就是每次粒子到达一个位置，都会通过指定的评价函数以及评价方式对该位置的解进行评价，如若此时得到的解的质量要优于之前的历史最优解，则更新该粒子的历史最优位置。所以在该算法不断优化的过程中，该粒子的最优位置也就在不断地改变，p B e s t pBestpBest 向量也就在不断地更新

全局最优(g B e s t gBestgBest) 除了需要进行局部最优的维护，另一个就是全局最优位置的维护，使用g B e s t gBestgBest来表示。有了局部最优就相当于狼群中的每只狼都有自己搜索得到最好的结果，然后传递给同伴，也就是整个狼群。狼群通过分析所有的局部最优结果，得到一个全局的最优结果。所以全局最优就是在局部最优内再挑选出其内最优的结果。这个全局最优将会对所有粒子都会有引导作用，使得整个种群往全局最优位置收敛。

PSO 的算法步骤如下：

种群初始化：初始化种群中各粒子的速度和位置，分别按照公式(1)和公式(2)。并将个体历史最优p B e s t pBestpBest设为当前值，将种群最优的个体作为当前g B e s t gBestgBest值。

迭代搜索，计算每次迭代后的粒子位置的适应度函数值(也就是评判粒子的质量)。

判断当前粒子的质量是否优于历史最优，如是，则将历史最优值设置为当前值；如不是，则不进行操作。

判断当前全局最优值是否优于历史的全局最优值，如是，则更新全局最优值；否则，不进行操作。

对粒子进行速度和位置的更新,更新方式如公式( 3 ) (3)(3)，公式( 4 ) (4)(4)。 如果没有达到结束条件，转到步骤(2)，否则输出 gBest并结束。

PSO 粒子依靠速度和位置在二维空间中更新示意图如下所示:

在公式(3)中，c 1 , c 2 c_1,c_2c 1 ,c 2 是加速系数(Acceleration Coefficients,也称学习因子)，在原论文中建议取固定值2.0。r a n d 1 d 和 r a n d 2 d rand_1^d和rand_2^drand 1d 和rand 2d是两个在区间[0,1]上的随机数。在更新过程中，需要对速度进行限制，速度一般限制在解的取值范围的10%~20%内。对于公式(4)，位置更新必须是合法的，所以在每次进行更新之后都要检查更新后的位置是否在问题空间中，如不在，则需要进行必要的修正，修正的方法可以是重新随机设定或者限定在边界。

4. 粒子群优化模型介绍

粒子群优化中，越过最优的个体会被拉回去，所有粒子都保留了良好解决方案的知识。此外粒子群优化是有记忆的，在运行过程中会把历史最优粒子的位置进行记录，并且在后面的运行过程中不断对后面的粒子进行引导。 在《A New Optimizer Using Particle Swarm Theory》 一文中，作者 Russell Eberhart 和 James Kennedy 进一步论述了粒子群优化算法，在该文中描述了使用粒子群方法对非线性函数的优化。 讨论并比较了两种范式的实现，包括最近开发的面向本地的范式，描述了这两种范例的基准测试。在这其中提出的两种范式便是粒子群优化算法(Particle Swarm Optimization,PSO)的两种形式，但是一般我们更容易见到其中一种，也就是GBEST 模型，在第一次描述PSO算法中，就是以GBEST为例进行举例。除却 GBEST模型，在该文中还提出另外一种模型，那就是LBEST 模型。下面，我们将细致地介绍这两种模型。

4.1 GBEST模型

在之前已经介绍了粒子群优化算法的思想来源以及定义的对照表，所以想必大家对gbest 是有印象的。在第一篇文章中，我们将其称之为全局最优解。标准的“ GBEST”粒子群算法非常简单，它是开发的粒子群优化的原始形式。

4.1.1 算法执行流程

Step1：初始化种群，也即粒子数N，同时初始化每个粒子的速度和初始位置，维度为解的维度。

Step2：初始化pbest以及gbest，并给其赋初始值，可选择取值边界。定义最大迭代次数以及起始迭代值。这个设置就需要根据语言而定，比如MATLAB，其下标是从1开始，所以一般起始迭代我们设置为1.

Step3：对每一个粒子位置进行适应性评估，倘若当前位置优于历史最优位置，则更新该粒子的历史最优解。

Step4：从所有粒子的历史最优解获取适应性最优的那个位置，并与当前全局最优位置gbest适应性进行比较，如果适应性大于当前位置，则更新全局最优位置。

Step5：按照公式更新粒子速度，更新速度公式为：v i d = v i d + c 1 × r a n d i d × ( p B e s t i d − x i d ) + c 2 × r a n d i d × ( g B e s t d − x i d ) v_i^d=v_i^d+c1\times rand_i^d\times (pBest_i^d-x_i^d)+c2\times rand_i^d\times (gBest^d-x_i^d)v

d =v i +c1×rand id ×(pBest d −x id )+c2×randid ×(gBest d −x id )Step6：判断粒子速度是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。

Step7：按照公式更新粒子位置，更新粒子位置公式为：x i d = x i d + v i d x_i^d=x_i^d+v_i^dx id =x id +v dStep8：判断粒子位置是否超过取值范围，将超过取值范围的重新赋值（可以值为边界值，或者重新随机生成）。

Step9：判断迭代次数是否大于最大迭代次数，如不大于，则跳至Step3。否则输出最优解

4.2 LBEST 模型

除了GBEST模型之外，在这之后还提出了LBEST模型。在这种模型中，粒子仅具有自己的信息以及邻居的信息，而不是整个团队的全局最优信息。也就是说粒子没有了全局最优导向粒子，在这个过程中只会接收按照排序位置相邻的粒子的信息。粒子朝着由“pbest”和“ lbest”定义的点移动，而不是朝着“pbest”和“gbest”位置进行移动。这样的模型，粒子的多样性得到了很大程度上的维护，但是当粒子的相邻位置适应度较差时，粒子的适应值也会变差，同时，优于粒子的分布过于分散，导致粒子即使发现较好位置，也没有足够的粒子进行开发。可能会使其整体解的质量下降。 那么LBEST模型是如何进行更新的呢？在邻域= 2模型中，粒子（i）将其误差值与粒子（i-1）和粒子（i + 1）进行比较。 在这里，我们使用领域=2的模型进行描述。它的基本步骤和GBEST类似，只不过在更新的时候略有变化，同时也不再需要gbest变量来存储全局最优。

4.2.1 算法执行流程

Step1：初始化种群，也即粒子数N，同时初始化每个粒子的速度和初始位置，维度为解的维度。

Step2：初始化pbest，并给其赋初始值，可选择取值边界。定义最大迭代次数以及起始迭代值。这个设置就需要根据语言而定，比如MATLAB，其下标是从1开始，所以一般起始迭代我们设置为1.

Step3：对每一个粒子位置进行适应性评估，倘若当前位置优于历史最优位置，则更新该粒子的历史最优解。

Step4：按照公式更新粒子速度，更新速度公式为：v i d = v i d + c 1 × r a n d i d × ( p B e s t i d − x i d ) + c 2 × r a n d i − 1 d × ( x i − 1 d − x i d ) + c 3 × r a n d i + 1 d × ( x i + 1 d − x i d ) v_i^d=v_i^d+c1\times rand_i^d\times (pBest_i^d-x_i^d)+c2\times rand_{i-1}^d\times (x_{i-1}^d-x_i^d)+c3\times rand_{i+1}^d\times (x_{i+1}^d-x_i^d)v

id =v id +c1×rand d ×(pBest id −x id )+c2×rand i−1

d ×(x i−1d −xid )+c3×rand i+1d ×(x +1d −x d );在这里需要注意，我们采取首尾相连的措施，也就是第一个粒子的前一个粒子是最后一个粒子，同理，最后一个粒子的后一个粒子就是第一个粒子。

Step5：判断粒子速度是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。

Step6：按照公式更新粒子位置，更新粒子位置公式为：x i d = x i d + v i d x_i^d=x_i^d+v_i^dx id =x id +vid

Step7：判断粒子位置是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。

Step8：判断迭代次数是否大于最大迭代次数，如不大于，则跳至Step3。否则输出最优解

4.2.2 非常重要

在邻域模型中，存在着很多邻域结构：例如二邻域、三邻域等多邻域结构，但是这些类似结构的编程方法都类似，无需多言。

但是在邻域模型中有一个非常重要的点，那就是对邻域的理解，或者说定义。一般说来可以分为两种：

第一种：逻辑上的邻域。在粒子群进行搜索的时候，每个粒子都会有一个位置下标，就如同我们给每个粒子进行编号一样，并且该编号不会随着后续程序的运行而发生改变。在该基础上实现多邻域粒子群优化代码会更加简单。但是会存在问题，那就是所谓的邻域粒子除了逻辑上相邻之外，在物理上并不相邻，因此在物理上或许并不能称为邻域。

第二种：物理上的邻域。从上面就可以很快明白，物理上的邻域就是在搜索过程中在粒子 ii 身边的粒子，因此在这种邻域结构中，就需要额外计算个粒子的相邻粒子，再根据计算的结果进行邻域更新。因此第二种邻域结构更符合正常的思维，但是要付出更大的空间以及时间。

4.3 ωω模型

Yuhui Shi 和 Russell Eberhart在1998年发表的一篇名叫《 A Modified Particle Swarm Optimizer》中对之前提出的GBEST模型进行进一步的分析，再该论文中，他们指出：该模型中的公式(1)是由三部分组成：第一部分是粒子的先前速度；第二和第三部分是有助于粒子速度变化的部分。通过公式(1)我们可以看到，如果没有c 1 × r a n d i d × ( p B e s t i d − x i d ) c1\times rand_i^d\times (pBest_i^d-x_i^d)c1×rand id ×(pBest id −x id )和c 2 × r a n d i d × ( g B e s t d − x i d ) c2\times rand_i^d\times(gBest^dx_i^d)c2×rand id ×(gBest d −x id )这两部分，那么粒子就会在速度所指向的方向直线来回运动，无法去探索更广泛的区域。除非解就在它们运行的直线上，还必须是在粒子所能到达的点的位置，否则无法寻找到最优解。

但是当没有第一部分时，粒子的飞行速度将只会受到局部最优位置以及全局最优位置的影响。当该是全局最优的时候，那么该粒子就会停在原地，直到其他粒子找到一个更优的解。

所以从上面的描述我们可以看出，当粒子速度更新公式中缺少第一部分的时候，那么该粒子群展示的是局部搜索能力：“因为此时所有粒子都会围绕在全局最优处”。而当保留第一部分的时候，它其实是扩大了粒子群搜索空间，使得他们更有能力搜索到新的区域。

因此，通过添加第一部分，更有可能具有全局搜索能力。局部搜索和全局搜索都有助于解决某些类型的问题。 全局搜索和局部搜索之间存在权衡对于不同的问题，局部搜索能力和全局搜索能力之间应该有不同的平衡。 考虑到这一点，将惯性权重w引入等式（1）中，如等式（2）所示。 这个 w 起到平衡全局搜索和局部搜索的作用，它可以是正常数，甚至是时间的正线性或非线性函数。

v i d = ω v i d + c 1 × r a n d i d × ( p B e s t i d − x i d ) + c 2 × r a n d i d × ( g B e s t d − x i d ) (2) v_i^d=\omega v_i^d+c1\times rand_i^d\times (pBest_i^d-x_i^d)+c2\times rand_i^d\times (gBest^d-x_i^d) \tag{2}vid =ωvid +c1×rand d ×(pBest i −x id )+c2×rand id ×(gBest d −x id )(2)同时，作者在论文也做了关测试，从测试的数据可以发现：“范围[0.9,1.21] 是使用ω \omegaω得到效果最好的区域。 ωω在这个范围内的 PSO 将有更大的概率在合理的迭代次数找到全局最优值。 在该文的测试中，发现当ωω=1.05时，PSO在迭代运行30次时都找到了全局最优值。 当ωω在[0.9,1.2] 的范围内取值时，发现当ω \omegaω = 0.9时，PSO 找到全局最优值的所需的平均迭代次数是最少的。”

原文中也指出最优的方式 ωω 一般是在[0.9,0.5]范围内线性递减

4.3.1 算法执行流程

该算法和 GBSET 模型的流程是一模一样的，唯一改变，就是速度更新公式由公式(1)改为公式(2).所以这里不再赘述该模型的算法流程。

5.源代码

%% =====================================================================%%
%% 粒子群优化：最早版本(doi:10.1109/icnn.1995.488968)
% Github：doFighter
% N: 种群大小
% dim:求解问题的维度
% x_max:解空间上限
% x_min:解空间下限
% iterate_max:最大迭代次数
% fitness:评价函数
%% --------------------------------------------------------------------%%
function[rest]=FunPSO_GBEST(N,dim,x_max,x_min,iterate_max,fitness)
    c = 2*ones(1,2); % 定义加速系数 c1,c2;这里直接放入一个矩阵中
    %由于在文中未能看到上面的 rand 是随机生成亦或是初始化生成，这里使用随机生成
    v_min = x_min*0.2;
    v_max = x_max*0.2;
    x = x_min + (x_max-x_min)*rand(N,dim); % 初始化种群位置；三行两列的矩阵，元素处于-10~10之间。三行寓意为种群，两列为维度.
    v = v_min + (v_max-v_min)*rand(N,dim);   % 初始化种群速度；速度一般在位置的取值范围的10%~20%，这里取10%
    pBest = x;  % 存储粒子的局部最优值,初始为 初始值
    gBest = x(1,:);  % 存储粒子的全局最优值
    for i = 2:N
       if fitness(gBest) > fitness(pBest(i,:))
          gBest = pBest(i,:); 
       end
    end
    iterate = 1;
    res = inf*ones(N,1);
    while iterate < iterate_max+2
        % 下面进行速度与位置的更新，当速度超出 -2~2 时，将其置为边界数，同理，当位置超出 边界 时，将其也置为边界数
        v = v + c(1)*rand(N,dim).*(pBest-x) + c(2)*rand(N,dim).*(gBest-x);
        v(v>v_max) = v_max;
        v(v<v_min) = v_min;
        x = x + v;
        x(x>x_max) = x_max;
        x(x<x_min) = x_min;
        
        for i = 1:N                     % 对每个粒子进行求解
            if fitness(x(i,:)) < fitness(pBest(i,:))   % 判断当前位置是否优于局部最优位置，若优于局部最优，则更新局部最优位置
                pBest(i,:) = x(i,:);
            end
            res(i) = fitness(pBest(i,:));
        end
        if min(res) < fitness(gBest)          % 判断当前所有的局部最优位置是否优于全局最优位置，若是，则更新
            position = find(res == min(res));
            gBest = pBest(position(1),:);
        end
        
        % 迭代计数器加1
        iterate = iterate + 1;
    end
    rest = fitness(gBest);
end



%% =====================================================================%%

%% 二邻域粒子群优化：逻辑相邻(doi:10.1109/MHS.1995.494215)
% Github：doFighter
% Encoding format：utf-8
% N: 种群大小
% dim:求解问题的维度
% x_max:解空间上限
% x_min:解空间下限
% iterate_max:最大迭代次数
% fitness:评价函数
%% --------------------------------------------------------------------%%
function[rest]=FunPSO_LBEST_Logic(N,dim,x_max,x_min,iterate_max,fitness)
    %% ====================================================================
    % LBEST版本，在该版本中，对应粒子会受到邻近粒子的影响而导致搜索方向的改变
    % 在该版本的粒子群算法中，要通过以下两个公式便可(以邻近粒子数等于2为例)
    % 1.v(i,d)=v(i,d)+c1*rand(1,d)*(pBest(i,d)-x(i,d)+c2*rand(2,d)*(x(i-1,d)-x(i,d)+c3*rand(3,d)*(x(i+1,d)-x(i,d));
    % 2.x(i,d)=x(i,d)+v(i,d);
    %% ====================================================================
    % 第一步：初始化变量
    c = 2*ones(1,3); % 定义加速系数 c1,c2;这里直接放入一个矩阵中
    %由于在文中未能看到上面的 rand 是随机生成亦或是初始化生成，这里使用随机生成
    v_min = x_min*0.1;
    v_max = x_max*0.1;
    x = x_min + (x_max-x_min)*rand(N,dim); % 初始化种群位置；三行两列的矩阵，元素处于-10~10之间。三行寓意为种群，两列为维度.
    v = v_min + (v_ma