【MATLAB第10期】基于MATLAB的粒子群及若干改进的 ...

DomDomm

一、经典粒子群PSO算法

1 思想来源

粒子群优化(Particle Swarm Optimization,PSO) 作为进化计算的一个分支，是由 Eberhart 和 Kennedy 于 1995 提出的一种搜索算法。该算法在提出时是为了优化非线性函数，亦即求解非线性函数在某个求解范围内的最优解。
粒子群优化是一种模拟自然界生物的活动以及群体智能的随机搜索算法，该算法吸取了人工生命(Artificial Life)、鸟群觅食(Birds Flocking)、鱼群学习(Fish Schooling) 和 群理论(Swarm Theory) 的思想，通过对以上群体智能行为的抽象模仿，从而缔造出粒子群优化 。从这我们也可以看出，粒子群优化，其本质就是一种 群体智能优化 ，该优化算法就是通过模拟生物群体在觅食或者躲避天敌时展现的群体智能行为而产生的。
2 基本原理

2.1 群体智能(Swarm Intelligence)
在自然界中，鸟类是通过什么养的机制发现猎物的呢？或许，我们很少见到鸟群协同寻找食物，但是在动物世界中，我们经常可以看到庞大的鱼群聚集在一起躲避被猎杀。以及狼群和猎狗等等，它们是如何发现猎物的呢？
当我们不断剖析其中原因，会发现，狼群以及猎狗在发现猎物之后会通过声音向同伴传递信息，使得同伴不断向猎物所在地不断聚集。而对于其他的亦是如此，虽然我们无法理解鸟群以及鱼群是如何进行协同配合来躲避天敌或是寻找食物，但是在它们之间一定是存在某种可供交流的机制，使得各个个体之间的信息、经验能够互享，随着群体中经验的累加，就会展现出一种超越个体的 群体智能 ，这也是我们为什么说人多力量大，其实人类也是具有群体智能的，并且我们老祖宗很早之前就已经得出这个结论。
2.2 算法原理
在上文我们已经分析了群体智能的基本原理，现在最关键的，便是如何将其要素进行抽象，最后变成能够进行编程的规则要求。
在之前的场景中已经提出，一群分散寻找食物的狼或者猎狗，我们以狼群为例，一开始并不知道猎物的所在，但是会有一个间接的机制会让狼群知道猎物的所在，比如猎物留下的足迹以及气味等。于是狼就会随着这些痕迹不断分散搜索，寻找气味最浓的位置方向，并通过声音或其他信息交流的方式通知其他狼，以此不断调整整个群体的位置，保证不断向猎物逼近，以至于最后发现猎物并捕食。这样，狼群中的每个个体都有了一个指导方向，总体的方向是猎物的方向，在这个过程中，所有个体不断地配合，根据自身的经验以及整个群体的经验调整自己的速度以及跟踪方向。

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

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

狼群捕食	粒子群优化算法
狼群	在搜索空间进行搜索的所有有效解(也就是种群规模N)
捕食空间	问题搜索空间(即解的维数，也就是未知数个数)
捕食范围	求解问题的范围(即解的取值范围)
追赶速度	解的速度向量$v_i=[v_i^1,v_i^2,\dotsb,v_i^D]$
自身位置	有效解的位置$x_i=[x_i^1,x_i^2,\dotsb,x_i^D]$
个体认知	粒子根据自身历史最优位置更新自身速度和位置
群体认知	粒子根据全局最优位置更新自身速度和位置
捕获猎物	算法结束，输出当前全局最优解

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

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

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

• 自身历史最优($pBest$) 在该算法中，每个粒子还需要维护自身的历史最优位置 pBest
  pBest
  。也就是每次粒子到达一个位置，都会通过指定的评价函数以及评价方式对该位置的解进行评价，如若此时得到的解的质量要优于之前的历史最优解，则更新该粒子的历史最优位置。所以在该算法不断优化的过程中，该粒子的最优位置也就在不断地改变，$pBest$ 向量也就在不断地更新
  
• 全局最优($gBest$) 除了需要进行局部最优的维护，另一个就是全局最优位置的维护，使用$gBest$来表示。有了局部最优就相当于狼群中的每只狼都有自己搜索得到最好的结果，然后传递给同伴，也就是整个狼群。狼群通过分析所有的局部最优结果，得到一个全局的最优结果。所以全局最优就是在局部最优内再挑选出其内最优的结果。这个全局最优将会对所有粒子都会有引导作用，使得整个种群往全局最优位置收敛。
  
• PSO 的算法步骤如下：
  

• 种群初始化：初始化种群中各粒子的速度和位置，分别按照公式(1)和公式(2)。并将个体历史最优$pBest$设为当前值，将种群最优的个体作为当前$gBest$值。
  
• 迭代搜索，计算每次迭代后的粒子位置的适应度函数值(也就是评判粒子的质量)。
  
• 判断当前粒子的质量是否优于历史最优，如是，则将历史最优值设置为当前值；如不是，则不进行操作。
  
• 判断当前全局最优值是否优于历史的全局最优值，如是，则更新全局最优值；否则，不进行操作。
  
• 对粒子进行速度和位置的更新,更新方式如公式$(3)$，公式$(4)$。
  
• 如果没有达到结束条件，转到步骤(2)，否则输出 gBest并结束。
  

一般在开始时，就会需要确定求解函数解的范围，所以$x_{min}以及x_{max}$就能在此时进行确定。而$v_{min}以及v_{max}$则一般取解的范围的10%~20%。


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


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



PSO流程图

在讲述了关于粒子群优化的原理，以及算法的整体流程结构，在这一部分，我们将会介绍在早期经典的PSO算法，并使用 MATLAB 以及 python 编写代码。
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+c1\times rand_i^d\times (pBest_i^d-x_i^d)+c2\times rand_i^d\times (gBest^d-x_i^d)$
  
• Step6：判断粒子速度是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。
  
• Step7：按照公式更新粒子位置，更新粒子位置公式为：$x_i^d=x_i^d+v_i^d$
  
• Step8：判断粒子位置是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。
  
• 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+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)$;在这里需要注意，我们采取首尾相连的措施，也就是第一个粒子的前一个粒子是最后一个粒子，同理，最后一个粒子的后一个粒子就是第一个粒子。
  
• Step5：判断粒子速度是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。
  
• Step6：按照公式更新粒子位置，更新粒子位置公式为：$x_i^d=x_i^d+v_i^d$
  
• Step7：判断粒子位置是否超过取值范围，将超过取值范围的重新赋值（可以赋值为边界值，或者重新随机生成）。
  
• Step8：判断迭代次数是否大于最大迭代次数，如不大于，则跳至Step3。否则输出最优解
  

4.2.2 非常重要
在邻域模型中，存在着很多邻域结构：例如二邻域、三邻域等多邻域结构，但是这些类似结构的编程方法都类似，无需多言。
但是在邻域模型中有一个非常重要的点，那就是对邻域的理解，或者说定义。一般说来可以分为两种：

• 第一种：逻辑上的邻域。在粒子群进行搜索的时候，每个粒子都会有一个位置下标，就如同我们给每个粒子进行编号一样，并且该编号不会随着后续程序的运行而发生改变。在该基础上实现多邻域粒子群优化代码会更加简单。但是会存在问题，那就是所谓的邻域粒子除了逻辑上相邻之外，在物理上并不相邻，因此在物理上或许并不能称为邻域。
  
• 第二种：物理上的邻域。从上面就可以很快明白，物理上的邻域就是在搜索过程中在粒子 ii 身边的粒子，因此在这种邻域结构中，就需要额外计算个粒子的相邻粒子，再根据计算的结果进行邻域更新。因此第二种邻域结构更符合正常的思维，但是要付出更大的空间以及时间。
  

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


但是当没有第一部分时，粒子的飞行速度将只会受到局部最优位置以及全局最优位置的影响。当该是全局最优的时候，那么该粒子就会停在原地，直到其他粒子找到一个更优的解。
所以从上面的描述我们可以看出，当粒子速度更新公式中缺少第一部分的时候，那么该粒子群展示的是局部搜索能力：“因为此时所有粒子都会围绕在全局最优处”。而当保留第一部分的时候，它其实是扩大了粒子群搜索空间，使得他们更有能力搜索到新的区域。
因此，通过添加第一部分，更有可能具有全局搜索能力。局部搜索和全局搜索都有助于解决某些类型的问题。 全局搜索和局部搜索之间存在权衡对于不同的问题，局部搜索能力和全局搜索能力之间应该有不同的平衡。 考虑到这一点，将惯性权重w引入等式（1）中，如等式（2）所示。 这个 w 起到平衡全局搜索和局部搜索的作用，它可以是正常数，甚至是时间的正线性或非线性函数。
$$ 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} $$ 同时，作者在论文也做了相关测试，从测试的数据可以发现：“范围[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_max-v_min)*rand(N,dim);   % 初始化种群速度；速度一般在位置的取值范围的10%~20%，这里取10%
    pBest = x;  % 存储粒子的局部最优值
    iterate = 1;
    results = inf*ones(iterate_max,1);
    res = inf*ones(N,1);

    % 第二步：迭代求最优
    while iterate < iterate_max+1
        % 下面进行速度与位置的更新，当速度超出 -2~2 时，将其置为边界数，同理，当位置超出 -10~10 时，将其也置为边界数
        for i=1:N
            pre = i - 1;
            nex = i + 1;
            if pre < 1
                pre = N;
            end
            if nex > N
                nex = 1;
            end
            v(i,:) = v(i,:) + c(1)*rand(1,dim).*(pBest(i,:)-x(i,:)) + c(2)*rand(1,dim).*(pBest(pre,:)-x(i,:)) + c(3)*rand(1,dim).*(pBest(nex,:)-x(i,:));
        end
        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
        results(iterate) = min(res);    % 保存当前最优解
        
        % 迭代计数器加1
        iterate = iterate + 1;
    end
    rest = min(results);
end
%% =====================================================================%%
%% 惯性权重线性递减粒子群优化(doi:10.1109/ICEC.1998.699146)
% Github：doFighter
% Encoding format：utf-8
% N: 种群大小
% dim:求解问题的维度
% x_max:解空间上限
% x_min:解空间下限
% iterate_max:最大迭代次数
% fitness:评价函数
%% --------------------------------------------------------------------%%
function[rest]=FunPSO_Omega(N,dim,x_max,x_min,iterate_max,fitness)
    c = 2*ones(1,2); % 定义加速系数 c1,c2;这里直接放入一个矩阵中
    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+1
        % 下面进行速度与位置的更新，当速度超出 -2~2 时，将其置为边界数，同理，当位置超出 -10~10 时，将其也置为边界数
        omega = 0.9 - 0.5*((iterate-1)/iterate_max);
        v = omega*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
参考文献

[1]Eberhart, R. and J. Kennedy (1995). A new optimizer using particle swarm theory. Mhs95 Sixth International Symposium on Micro Machine & Human Science. [2]Shi, Y. (1998). A Modified Particle Swarm Optimizer. Proc of IEEE Icec Conference.
[2]Kennedy, J. and R. Eberhart (1995). Particle Swarm Optimization. Icnn95-international Conference on Neural Networks. IEEE, 1995, 4: 1942-1948. [2]张军,詹志辉等.计算智能[M].北京:清华大学出版社,2009.
二、一种基于 Sigmoid 函数的自适应加权粒子群优化器(AWPSO)

1 引言

在该文中，提出了一种新的基于Sigmoid函数的自适应加权粒子群优化 (PSO) 算法。

• 算法主要创新：开发了一种基于 sigmoid 函数的加权策略来自适应地调整加速度系数。
  

新提出的自适应加权策略同时考虑了粒子到全局最佳位置的距离和粒子到其个人最佳位置的距离，从而具有提高收敛速度的显着特点。受神经网络激活函数的启发，采用新策略通过 sigmoid 函数更新加速系数。

论文DOI:10.1109/TCYB.2019.2925015
2 算法介绍

2.1 改进起源
在PSO算法中，加速度系数被用来激励粒子移动到pbest和gbest。每个粒子的位置到它的pbest和gbest的距离在确定粒子的运动时起决定性作用。另一方面，控制参数的适应是寻求具有令人信服的效率和准确性的最佳解决方案的重要因素。因此，为了有效地控制PSO算法，在论文中，作者努力提出一种新的自适应加权机制，随着迭代的进行，加速度系数可以自适应地调整。
2.2 权重自适应策略
在经典的PSO算法中，单个粒子的速度会根据粒子到其pbest和gbest的距离而加速。因此，选择合适的加速度系数对于通过问题空间找到全局最优解至关重要。在这种情况下，根据上述距离自适应地更新加速系数迭代，以有效提高 PSO 算法的搜索能力，具有理论和实践意义。
在文中，作者提到，在很多的PSO变体仅以实时迭代次数的方式调整加速度系数(例如公式(1)所示)，而没有考虑种群演化的信息。


作者提出，在进化过程的早期，鼓励所有个体尽可能多地探索整个搜索空间。然后，在优化过程的后期，个体被激励收敛到全局最优并尽快找到优化解决方案。从经典粒子速度更新公式中可以看出，粒子更新的速度取决于粒子到它们自己的pbest和gbest的距离。在这种情况下，根据每个粒子到其pbest和gbest的距离来调整加速度系数是合理的。
考虑到所有上述问题，作者提出了一种自适应加权策略来自适应地控制加速度系数。

• 主要作用：加快粒子搜索以尽可能快地找到最优解，从而提高收敛速度。
  
• 区别：与实时更新策略不同，加速度系数根据粒子与其gbest和pbest的距离而改变。如果粒子远离其pbest和gbest，则使用相对较大的加速系数来加速粒子。
  

同时，加速度系数的值被限制在一个适当的范围内，以避免过早收敛，这意味着速度应该是有界的，以保证算法的搜索能力。受前面讨论的启发，作者认为自适应加权更新函数适合描述加速度系数和距离（从粒子到它的pbest和gbest）之间的关系。换句话说，前一个加速度系数的更新应该适应后一个距离，从而充分证明粒子运动的速度朝向全局最优。从数学的角度来看，所提出的自适应加权更新规则可以描述如下：


其中，函数$F(·)$表示自适应加权更新函数，后面将会介绍；$g_{pi}(k)$和$g_{gi}(k)$定义为: $$ \begin{cases} g_{pi}(k)=p_i(k)-x_i(k)\ g_{gi}(k)=g_i(k)-x_i(k) \end{cases}\tag{3} $$ 分别表示第$k$次迭代时粒子$i$到其$pbest$和$gbest$的距离。
2.3 自适应加权更新函数的选择
直观地说，自适应加权更新函数应该具有以下两个性质：

• 更新函数是单调递增的；
  
• 更新函数是有界的。
  

第一个属性主要是由于加速度系数的特性。众所周知，加速度系数是将粒子拉到pbest和gbest的加权项。远离其pbest和gbest的粒子需要快速移动到其pbest和gbest。因此，需要一个单调递增的函数。
第二个属性是由约束优化问题的搜索空间是正常有界的事实证明的。一旦一个粒子接近它的pbest和gbest，运动应该放慢，以避免错过它的pbest和gbest。因此，加速度系数应该有界以控制粒子的速度。
在寻找既单调递增又统一有界的适当更新函数时，神经网络中使用的激活函数似乎是理想的候选者。神经网络有一些流行的激活函数，例如阶跃函数和sigmoid函数，其中我们决定选择sigmoid函数作为自适应加权更新函数，原因有以下三个：

• sigmoid函数是单调且有界的；
  
• sigmoid函数的曲线是S形的，这样可以避免控制参数发生不希望的突变；
  
• sigmoid函数是平滑可微的，从而反映了权重更新迭代的自适应/动态特性。
  

因此，该文采用sigmoid函数来调整加速度系数，如公式(4)所示。


其中e是自然对数，a表示曲线的陡度，是一个常数值，b表示曲线的峰值，c表示曲线中心点的横坐标值，d是一个正常数值， D是由(3)确定的函数的输入。具体来说，D是粒子与其认知加速系数的pbest之间的距离。对于社会加速度系数，D表示粒子与gbest之间的距离。
在上面的描述中，例如b表示的不是公式(4)的峰值，公式(4)的峰值为b+d。但是总的来说，这并不会影响我们对论文的理解和复现。
总之，所提出的基于 sigmoid 函数的自适应加权策略的三个主要优点总结如下。

• 加速度系数自适应控制在合理范围内，自适应加权策略保证了速度更新过程的效率。
  
• 自适应加权更新函数，选择sigmoid函数，用于反映加速度系数的单调但相对平滑的变化，其中较大的距离将导致较大的加速度系数值。
  
• 激发粒子尽可能快地寻求最优解，从而提高准确性和收敛性。
  

2.4 AWPSO 算法的框架
在该文中，速度及位置更新方式如下(5)，(6)式所示。


其中$\omega$是惯性权重；$g_{pi}(k)$和$g_{gi}(k)$分别表示第$k$次迭代时粒子$i$到其$pbest$和$gbest$的距离；$c_{g_{pi}}(k)$表示由$g_{pi}(k)$确定的加速度常数；$c_{g_{gi}}(k)$表示由$g_{gi}(k)$确定的加速度常数。 其中需要稍微指出的是，$r_1,r_2$是随机数。
流程图如下所示：


3 参数设置

在原文中，参数设置分别如下：

• a=0.000035m
  
• b=0.5
  
• c=0
  
• d=1.5
  

其中m表示优化问题的搜索范围。

4 代码复现

在该文中，对于距离的计算并没有较为清晰的说明。首先我们通过分析文中的描述以及资料可知，加速系数$c_1,c_2$是一个常数，因此距离计算变成了此时最重要的问题，由于文中并未描述因此在实现时，我对以下计算公式进行计算。

上面三个公式中$dim$代表问题的求解维度。

公式(7)和公式(9)可以看成是同一个计算方法，因此我们只需要验证公式(7)和公式(8)。
但是从根本上说，距离计算起始应该是公式(9)中描述的，这也是我们常用的，因此在这里我还是倾向于公式(9)的距离计算方法，当然，上面已经说了，公式(9)可以化简为公式(7)的模样，因此为了方便代码编写，可以直接使用公式(7)。
5 参考文献

[1]W. Liu, Z. Wang, Y. Yuan, N. Zeng, K. Hone and X. Liu, "A Novel Sigmoid-Function-Based Adaptive Weighted Particle Swarm Optimizer," in IEEE Transactions on Cybernetics, vol. 51, no. 2, pp. 1085-1093, Feb. 2021, doi: 10.1109/TCYB.2019.2925015.
三、双中心粒子群优化（DCPSO）

1 引言

粒子群优化（$PSO$）算法是一种新兴的群体智能优化技术，由于其原理简单、参数少、效果好等 优点已经广泛应用于求解各类复杂优化问题．而影响该算法收敛速度和精度的２个主要因素是粒子个 体极值与全局极值的更新方式．通过分析粒子的飞行轨迹和引入广义中心粒子和狭义中心粒子，提出 双中心粒子群优化（$double; center; particle; swarm; optimization,DCPSO$）算法，在不增加算法复杂度条 件下对粒子的个体极值和全局极值更新方式进行更新，从而改善了算法的收敛速度和精度。
在双中心粒子群优化中，主要有以下改变：

• 改变轨迹探索方式
  
• 引入广义中心粒子
  
• 引入狭义中心粒子
  

2 算法精讲

首先需要明确的是，该算法基本沿用 A;Modified;Particle;Swarm;Optimizer（这里称之为OPSO）中的算法结构，在该结构的基础上进行了三点改变：

• 飞行轨迹搜索改变
  
• 广义中心粒子
  
• 狭义中心粒子
  

因此，只要理解更改的部分便能理解该算法。而这三点改变分别体现在两处，分别是：

• 个体最优位置更新
  
• 全局最优位置更新
  

2.1 个体最优位置更新
在 OPSO中，粒子的更新方式为：


因此，最终位置是通过上面的三项进行合成得到。在双中心粒子群优化中，个体最优位置更新是将上面三项进行拆分，在矢量合成的过程中，会经过三个位置，分别为：


如下图所示：


在上图中粒子飞行路径为 xoldABxnew，就单个粒子而言，折线运动经历的区域比直线方式复 杂，其捕食范围能够被最大程度地覆盖在整个搜索 空间内，体现在$PSO$算法中，每个粒子能够在更为 广阔的搜索空间内寻找最优解，这样粒子的个体极 值就能够得到显著的改善；从粒子群整体来说，每个粒子个体极值的改善将直接影响着种群全局极值的增强，这一影响使得种群能更快地收敛于最优解．
因此在更新时分别会对这三个位置进行比较更新，更新公式如下：

注意事项：在个体最优位置更新时需要注意，在更新规则方面论文没有详细描述，通过有限测试，在这几个阶段速度不进行约束，但是位置需要约束，不能超过边界。并且最后要对速度进行归正，即执行公式（1），并且需要对最后的位置进行约束。该方式效果会更好！

• 探讨
  

其实如果深究，就会发现如果在算法中对速度进行约束，那么肯定就无法得到理想中的三个位置，因为在 OPSO中更新是整体更新，当速度超出范围就会被限制，但是我们无法知道在哪一处被限制。当然，倘若只对位置进行约束，那的确是不会出现上述现象。但是在参考论文中，速度是有限制的。因此我们无法获得准确的轨迹位置，而是依据上面的思想进行编写！
2.2 全局最优位置更新
在全局位置更新中，双中心粒子群优化引入了广义中心粒子($general center particle,GCP$)及狭义中心粒子($special center particle,SCP$)，两者计算公式为：


最终全局最优位置更新从所有的历史最优位置以及广义中心粒子、狭义中心粒子中竞选，公式如(4)式：


3 参考文献

[1].汤可宗,柳炳祥,杨静宇,孙廷凯.双中心粒子群优化算法[J].计算机研究与发展,2012,49(05):1086-1094.