具有万有引力加速机理的布谷鸟搜索算法

一、理论基础

1、布谷鸟算法(CS)及其存在的局限性

布谷鸟采用一种特殊的寄生宿主巢穴的方式孕育繁殖，它将孵育的蛋置入寄生宿主的巢穴，让寄生宿主孵化布谷鸟蛋。由于布谷鸟幼雏能发出比寄生宿主幼雏更闪亮的叫声，因此获得更多的食物，具有更高的存活率。在某些情形下，寄生宿主发现巢穴内不是宿主幼雏，则会遗弃该巢穴，并重新选择新的巢穴孵育繁殖。Yang X. S.和Deb等人基于上述孵育寄生机理，提出布谷鸟算法。它遵循以下3个理想规则：
规则1. 每只布谷鸟1次有且仅有孵化1个蛋，并随机选择1个寄生宿主巢穴存放。
规则2. 在随机选择的1个寄生宿主巢穴中，最优的寄生宿主巢穴将被保留至下一代。
规则3. 可利用的寄生宿主巢穴数量是固定的，1个寄生宿主巢穴的宿主发现寄生布谷鸟蛋的概率为$P_a$。
基于上述3个理想规则，可以得到1个宿主巢穴代表1个候选解。因此，布谷鸟算法的基本算法流程包含 3部分：首先，在当前候选解的基础上，以Levy飞行随机游动生成新的候选解，评价并保留较好的候选解；其次，按照发现概率$P_a$舍弃部分候选解；最后，按照偏好随机游动方式产生新的候选解并替代舍弃的解，保留较好解后进入下一次进化，直至满足算法的收敛条件。
设布谷鸟算法进化至第$r$代时第$m$个候选解为${\mathbf{X}}_{r,m}$，${\mathbf{X}}_{r,m}=(x_{r,m}^{(1)},x_{r,m}^{(2)},\cdots ,x_{r,m}^{(j)},\cdots ,x_{r,m}^{(D)}),j\in [1,D]$。采用Levy飞行随机游动产生新的个体(候选解)${\mathbf{X}}_{r+1,m}$更新律的表达式为$${\mathbf{X}}_{r+1,m}={\mathbf{X}}_{r,m}+({\mathbf{X}}_{r,m}-{\mathbf{X}}_{r,gb})\cdot (\gamma \oplus L(\beta ))\text{\hspace{1em}(1)}$$其中，${\mathbf{X}}_{r,gb}$为当前搜索到的全局最优解，$L(\beta )$表示Levy飞行随机游动路径，${\gamma }_0$表示初始搜索步长，$\oplus$表示点对点乘法，$t$为飞行时间：$$L(\beta )\sim \vartheta =t^{-1-\beta },0<\beta \le 2\text{\hspace{1em}(2)}$$通过数学代换，式(2)等价于：$$L(\beta )\sim \frac{\vartheta }{|v{|}^{1/\beta }}\cdot {\left(\frac{\Gamma (1+\beta )\sin (\pi \cdot \beta /2)}{\Gamma \left(\frac{1+\beta }{2}\right)\cdot \beta \cdot 2^{(\beta -1)/2}}\right)}^{1/\beta }\text{\hspace{1em}(3)}$$其中，$\Gamma (\cdot )$为伽玛函数，取$\beta =1.5$；$\vartheta ,v$为标准高斯分布随机数。
在偏好随机游动方式中，按照发现概率$P_a$舍弃部分候选解后，生成相同数量的新解：$${\mathbf{X}}_{r+1,m}={\mathbf{X}}_{r,m}+\varphi \cdot ({\mathbf{X}}_{r,k}-{\mathbf{X}}_{r,s})\text{\hspace{1em}(4)}$$其中，$\varphi$为算法的控制缩放系数，满足$\varphi \in U(0,1)$；${\mathbf{X}}_{r,k}$和${\mathbf{X}}_{r,s}$分别为第$r$代时第$k$个和第$s$个随机候选解。
基本CS算法虽然具有一定的收敛能力，但是存在以下局限性：
（1）采用Levy飞行随机游动产生新候选解的搜索步长设置困难：如果搜索步长较大，则寻优速度增大，但由于Levy飞行具有无限跳跃特性，使得算法执行过程中极易跳过全局最优解；如果搜索步长较小，虽然能以较大概率获得全局最优解，但是寻优速度会降低。考虑到Levy飞行具有无限大方差并且其增量服从二阶分布，因此，设计合理的搜索步长变得尤为困难。
（2）CS算法采用发现概率$P_a$舍弃偏好随机游走方式产生新的候选解。这种方法虽然加强了算法的全局搜索性能，但随机处理会舍弃较好的种群个体，而保留适应度较差的种群个体，进而降低算法搜索速度和收敛精度。
（3）CS算法仅通过两个随机个体的位置向量差获得下次进化的方向，未能充分利用当前宿主巢穴的足够信息，因此在搜索过程中具有盲目性，不能对种群个体进行局部精细化搜索，导致算法具有较差的局部搜索能力。纵然算法具有一定的全局搜索能力，但是仍然存在全局搜索和局部搜索不平衡的问题。

2、万有引力加速机理的布谷鸟算法(GASCS)

2.1 万有引力算法

受万有引力定律启发，伊朗学者Esmat Rashedi等人于2009年提出了万有引力搜索算法(Gravitational search algorithm, GSA)。GSA算法的原理是：物质个体之间由于引力作用相互吸引，并且引力沿质量较大的个体方向移动。物质个体的质量决定其万有引力的大小，质量越大的个体所受的引力越大。在外有引力的作用下，物质个体位置不断变化，最终，整个群体都聚集在物质个体质量最大的周围，即逼近优化问题的全局最优解。
在万有引力算法中，借鉴理论物理学中质量定义的3种属性描述函数适应度值，即主动引力质量$M_a$、被动引力质量$M_p$和惯性质量$M_i$。每个物质个体的位置对应优化问题解，其万有引力和惯性质量$M_i$共同决定相应的函数适应度数值。事实上，GSA算法通过调节$M_a,M_p$和$M_i$来控制算法进化，促使群体中的物质个体聚集在万有引力最大的物质个体附近，进而搜索到最优解。
GSA算法模型包含物质个体的质量计算、引力计算、加速度计算和速度更新与位置更新。算法首先对解空间内的物质个体位置和速度初始化。设$D$维空间维度中进化至第$r$代时第$m$个物质个体的位置和速度分别为$${\mathbf{X}}_{r,m}=(x_{r,m}^{(1)},x_{r,m}^{(2)},\cdots ,x_{r,m}^{(j)},\cdots ,x_{r,m}^{(D)}),{\mathbf{V}}_m=(v_{r,m}^{(1)},v_{r,m}^{(2)},\cdots ,v_{r,m}^{(j)},\cdots ,v_{r,m}^{(D)}),$$其中，$j=1,2,\cdots ,D,x_{r,m}^{(j)}$和$v_{r,m}^{(j)}$分别表示第$m$个物质个体在第$j$维的位置分量和速度分量。通过物质个体位置对应的适应度数值，确定物质个体的质量$M_i$和受到的万有引力。根据牛顿第二定律计算物质个体的加速度，并更新物质个体位置${\mathbf{X}}_{r,m}$和速度${\mathbf{V}}_{r,m}$。GSA算法主要包括以下4个步骤。
（1）物质个体质量计算
物质个体$m$的质量定义为$$q_{r,m}=\frac{f_{r,m}-f_{r,worst}}{f_{r,best}-f_{r,worst}}\text{\hspace{1em}(5)}$$$$M_{r,m}=\frac{q_{r,m}}{\sum _{m=1}^N{q}_{r,m}}\text{\hspace{1em}(6)}$$其中，$f_{r,m}$和$M_{r,m}$分别表示GSA算法第$r$次进化时，物质个体$m$的函数适应度值和相应的质量。假设所求解的为最小化问题，$f_{r,best}$和$f_{r,worst}$分别表示第$r$次进化时物质个体最优的适应度数值和最差适应度数值，其数学表达式如下： f r , b e s t = min ⁡ m ∈ { 1 , 2 , ⋯   , N } f r , m (7) f_{r,best}=\min_{m\in\{1,2,\cdots,N\}}f_{r,m}\tag{7} fr,best​=m∈{1,2,⋯,N}min​fr,m​(7) f r , w o r s t = max ⁡ m ∈ { 1 , 2 , ⋯   , N } f r , m (8) f_{r,worst}=\max_{m\in\{1,2,\cdots,N\}}f_{r,m}\tag{8} fr,worst​=m∈{1,2,⋯,N}max​fr,m​(8)（2）物质个体万有引力计算
根据万有引力定律，在维度$j$上，物质个体$m$对物质个体$k$的万有引力定义如下：$$F_{r,mk}^{(j)}=G_r\cdot \frac{M_{r,pm}\cdot M_{r,ak}}{R_{r,mk}+\epsilon }\cdot (x_{r,m}^{(j)}-x_{r,k}^{(j)})\text{\hspace{1em}(9)}$$其中，$\epsilon$为一个无穷小的常数；$M_{r,pm}$表示被作用于物质个体$m$第$r$次进化时的惯性质量，$M_{r,ak}$表示为作用于物质个体$m$第$r$次进化时的惯性质量，并且$M_{r,pm}=M_{r,ak}=M_{r,m}$；$R_{r,mk}$为物质个体$m$和物质个体$k$的欧氏空间距离，即$R_{r,mk}=||{\mathbf{X}}_{r,m}-{\mathbf{X}}_{r,k}|{|}_2$；$G_r$表示第$r$次进化时物质个体的万有引力常数，其表达式如式(10)所示：$$G_r=G(G_0,r)=G_0\cdot e^{-\theta \cdot r/W}\text{\hspace{1em}(10)}$$其中，$G_0$为进化初始时物质个体的万有引力系数；$\theta$为算法控制参数，一般取$\theta =20$；$W$为算法进化代数的最大值。
在维度$j$上，物质个体$m$所受的万有引力合力为$F_{r,m}^{(j)}$：$$F_{r,m}^{(j)}=\sum _{k\in b,k\ne m}^N{d}_j\cdot F_{r,mk}^{(j)}\text{\hspace{1em}(11)}$$其中，$d_j$为区间$(0,1)$内均匀分布的一个随机数，$b$为物质个体数目。
（3）物质个体加速度计算
由牛顿第二定律，物质个体$m$在维度$j$上第$r$次进化时的加速度定义如下：$$a_{r,m}^{(j)}=\frac{F_{r,m}^{(j)}}{M_{r,mm}}\text{\hspace{1em}(12)}$$其中，$M_{r,mm}=M_{r,pm}=M_{r,ak}=M_{r,m}$。
（4）物质个体更新速度和位置
算法在每次进化过程中，物质个体按照式(13)和式(14)分别更新速度$v_{r,m}^{(j)}$和位置$x_{r,m}^{(j)}$，其具体的数学表达式分别为：$$v_{r+1,m}^{(j)}=d_j\cdot v_{r,m}^{(j)}+a_{r,m}^{(j)}\text{\hspace{1em}(13)}$$$$x_{r+1,m}^{(j)}=x_{r,m}^{(j)}+v_{r+1,m}^{(j)}\text{\hspace{1em}(14)}$$

2.2 GASCS算法

传统的CS算法是基于Levy飞行随机游动方式和偏好随机游动方式两种搜索机理，虽然具有一定的收敛性能，但仍然存在第1节所阐述的3点局限性。为此，文献[1]提出具有万有引力加速机理的布谷鸟算法。它基于万有引力搜索无需学习外部环境因素的变化亦能感知全局最优信息的特点，将布谷鸟巢穴赋予不同的个体质量，其在优化过程中不仅遵循Levy飞行规律，而且遵循万有引力定律。
根据牛顿第二定律，宿主巢穴在进化过程中所受的加速度由巢穴所受的合外力及自身质量共同决定。基于此，将布谷鸟算法的个体更新律(1)和更新律(4)改进为式(15)和式(16)：$${\mathbf{X}}_{r+1,m}={\mathbf{a}}_r-({\gamma }_0\oplus L(\beta ))\cdot ({\mathbf{a}}_r-{\mathbf{X}}_{r,gb})\text{\hspace{1em}(15)}$$$${\mathbf{X}}_{r+1,m}={\mathbf{X}}_{r,gb}+p\cdot ({\mathbf{X}}_{r,k}-{\mathbf{X}}_{r,s}-{\mathbf{a}}_r)\text{\hspace{1em}(16)}$$其中，式(15)和式(16)分别表示Levy飞行随机游动和偏好随机游动方式的个体位置，$p\in (0,+\infty )$，${\mathbf{a}}_r$为${\mathbf{F}}_{r,t}$产生的合加速度。
针对第3点局限性中存在全局搜索和局部搜索不平衡的问题，GASCS算法提出一种概率变异方法增大种群进化的多样性，避免算法陷入局部搜索以及算法执行后期出现的迟滞现象。
假设${\mathbf{X}}_{r,m}=(x_{r,m}^{(1)},x_{r,m}^{(2)},\cdots ,x_{r,m}^{(j)},\cdots ,x_{r,m}^{(D)})$，并且$x_{r,m}^{(j)}\in [l^{(j)},u^{(j)}]$。执行概率变异操作时，首先从${\mathbf{X}}_{r,m}$中以概率$1/D$选取种群个体$x_{r,m}^{(j)},j\in [1,D]$；其次，利用式(17)的${\tilde{x}}_{r,m}^{(j)}$取代$x_{r,m}^{(j)}$，其中，$k$为$[1,D]$内的随机整数，$\xi$为$(0,1)$内的均匀分布随机数，$l^{(j)}$和$u^{(j)}$分别为$x_{r,m}^{(j)}$的上界和下界；最后，通过概率变异将候选解${\mathbf{X}}_{r,m}$更替为${\tilde{\mathbf{X}}}_{r,m}$，并进入算法执行过程，其中，${\tilde{\mathbf{X}}}_{r,m}=({\tilde{x}}_{r,m}^{(1)},{\tilde{x}}_{r,m}^{(2)},\cdots ,{\tilde{x}}_{r,m}^{(j)},\cdots ,{\tilde{x}}_{r,m}^{(D)})$。经过概率变异后的种群进化获得更大的多样性，避免算法寻优过程单一化而造成的全局搜索效率较低的问题。通过两个随机个体的位置向量差、当前搜索到的最优解及群体中宿主巢穴所受到的万有引力加速度，获得下一次进化的方向，充分利用宿主巢穴及群体内部进化过程出现的有效信息：$${\tilde{x}}_{r,m}^{(j)}=\{\begin{array}{l}(u^{(j)}+l^{(j)})/2+(1-2\cdot \xi )\cdot (u^{(j)}-l^{(j)})/2,j=k\\ x_{r,m}^{(j)},j\ne k\end{array}\text{\hspace{1em}(17)}$$GASCS算法在进化过程中不仅利用了Levy飞行随机游动方式的无限跳跃特性易于跳出局部最优解的特点，而且采用概率变异增大算法的多样性，确保算法能够跳离局部最优解。因此，GASCS算法较好地解决了CS算法存在全局搜索和局部搜索不平衡的问题。

2.3 GASCS算法实现

基于CS算法存在的3点局限性，文献[1]提出了具有万有引力加速机理的GASCS算法。该算法通过当前搜索到的最优解、Levy飞行随机游动和偏好随机游动方式产生的个体更新位置及宿主巢穴所受到的万有引力加速度共同作用，获得下一次进化的方向及最佳巢穴位置，并采用概率变异方法，引导宿主巢穴沿全局最优解方向移动，进而搜索到全局最优解。
GASCS算法的执行步骤如下所示。
步骤1. 算法初始化，种群为$N$个宿主巢穴，优化问题的空间维度$D$，最大进化代数$W$，发现概率 $P_a$，初始时刻万有引力系数$G_0$，算法控制参数$\theta$，宿主巢穴的初始移动速度${\mathbf{V}}_{r,m}$，其中，进化代数$r=1$。
步骤2. 计算随机化宿主巢穴位置${\mathbf{X}}_{r,m}(m=1,2,\cdots ,N)$对应的适应度函数值$f({\mathbf{X}}_{r,m})$。
步骤3. 由式(10)计算$G_r$，同时，由式(7)和式(8)计算当前最优值$f_{r,best}$和最差值$f_{r,worst}$，以及对应的最优解${\mathbf{X}}_{r,gb}$。
步骤4. 根据式(5)和式(6)计算宿主巢穴的质量$q_{r,m},M_{r,m}$。
步骤5. 分别根据式(11)和式(12)计算当前进化代数的宿主巢穴所受到的万有引力合力$F_{r,m}^{(j)}$和加速度$a_{r,m}^{(j)}$。
步骤6. 采用式(15)的Levy飞行随机游动方式产生新的宿主巢穴，按照发现概率$P_a$舍弃候选解${\mathbf{X}}_{r+1,m}$。
步骤7. 采用式(16)的偏好随机游动方式产生新的宿主巢穴，替换步骤6中被舍弃的候选解。
步骤8. 采用概率变异方法产生候选解${\tilde{\mathbf{X}}}_{r+1,m}$对应的函数适应度值$f({\tilde{\mathbf{X}}}_{r+1,m})$。
步骤9. 计算步骤8中种群产生的候选解对应的函数适应度值$f({\tilde{\mathbf{X}}}_{r+1,m})$，同时更新当前最优值$f_{r+1,best}$和最差值$f_{r+1,worst}$以及对应的最优解${\mathbf{X}}_{r+1,gb}$。
步骤10. 若满足算法终止条件，则输出当前进化的最优值及最优解，并停止算法；否则，转向步骤3继续执行算法。

二、实验仿真及分析

将GASCS与GSA、CS、PSO和ABC进行对比，以常用23个测试函数中的F1、F5(单峰函数/30维)、F10、F13(多峰函数/30维)、F14、F16(固定维度多峰函数/2维、2维)为例，实验设置种群规模为20，最大迭代次数为400，每种算法独立运算30次，结果显示如下：

• 函数适应度值的收敛曲线对比及收敛结果

函数：F1
GSA：最优值: 0.23077, 最差值: 588.5136, 平均值: 81.2518, 标准差: 122.5703
CS：最优值: 15.7073, 最差值: 128.0775, 平均值: 57.6246, 标准差: 23.681
PSO：最优值: 162.3483, 最差值: 940.4387, 平均值: 456.4764, 标准差: 184.8215
ABC：最优值: 0.067679, 最差值: 1184.4898, 平均值: 159.9748, 标准差: 274.4239
GASCS：最优值: 1.5303e-149, 最差值: 1.4195e-125, 平均值: 4.9421e-127, 标准差: 2.5901e-126
函数：F5
GSA：最优值: 29.3923, 最差值: 1088.5277, 平均值: 197.5399, 标准差: 183.1622
CS：最优值: 968.9382, 最差值: 107929.0128, 平均值: 8600.9133, 标准差: 19473.3501
PSO：最优值: 4539.3427, 最差值: 135094.6679, 平均值: 35655.2356, 标准差: 28678.7236
ABC：最优值: 12.9133, 最差值: 11851.841, 平均值: 1325.5531, 标准差: 2925.8172
GASCS：最优值: 27.7906, 最差值: 28.8253, 平均值: 28.6015, 标准差: 0.18926
函数：F9
GSA：最优值: 36.8135, 最差值: 123.5968, 平均值: 76.1909, 标准差: 22.7362
CS：最优值: 71.9859, 最差值: 136.2665, 平均值: 107.9521, 标准差: 15.6606
PSO：最优值: 66.685, 最差值: 165.3538, 平均值: 101.6716, 标准差: 23.317
ABC：最优值: 278.7571, 最差值: 340.2926, 平均值: 317.6587, 标准差: 15.2433
GASCS：最优值: 0, 最差值: 0, 平均值: 0, 标准差: 0
函数：F13
GSA：最优值: 9.6037, 最差值: 41.6958, 平均值: 27.8702, 标准差: 9.3281
CS：最优值: 10.2221, 最差值: 390.2833, 平均值: 55.3644, 标准差: 71.777
PSO：最优值: 35.7413, 最差值: 30597.4275, 平均值: 5411.6015, 标准差: 9015.082
ABC：最优值: 0.039083, 最差值: 34.6121, 平均值: 6.3655, 标准差: 7.976
GASCS：最优值: 0.26526, 最差值: 1.2775, 平均值: 0.65661, 标准差: 0.28374
函数：F14
GSA：最优值: 0.998, 最差值: 13.0697, 平均值: 5.2528, 标准差: 3.6007
CS：最优值: 0.998, 最差值: 0.998, 平均值: 0.998, 标准差: 1.2555e-12
PSO：最优值: 0.998, 最差值: 13.6186, 平均值: 4.7318, 标准差: 3.3877
ABC：最优值: 0.998, 最差值: 0.998, 平均值: 0.998, 标准差: 2.0623e-15
GASCS：最优值: 0.998, 最差值: 0.998, 平均值: 0.998, 标准差: 1.2333e-12
函数：F16
GSA：最优值: -1.0316, 最差值: -1.0316, 平均值: -1.0316, 标准差: 4.7012e-16
CS：最优值: -1.0316, 最差值: -1.0316, 平均值: -1.0316, 标准差: 1.5863e-15
PSO：最优值: -1.0316, 最差值: -1.0316, 平均值: -1.0316, 标准差: 6.4539e-16
ABC：最优值: -1.0316, 最差值: -1.0316, 平均值: -1.0316, 标准差: 1.0006e-05
GASCS：最优值: -1.0316, 最差值: -1.0316, 平均值: -1.0316, 标准差: 1.6951e-10

• 独立运行30次的函数最优适应度值比较特性图

实验结果表明：提出的GASCS算法与其他算法相比具有更优秀的全局寻优性能、更好的鲁棒性、更快的搜索速度和更高的收敛精度。

三、参考文献

[1] 傅文渊. 具有万有引力加速机理的布谷鸟搜索算法[J]. 软件学报, 2021, 32(5): 1480-1494.
[2] Yang X S, Deb S. Engineering Optimisation by Cuckoo Search[J]. International Journal of Mathematical Modelling & Numerical Optimisation, 2010, 1(4): 330-343.
[3] Esmat Rashedi, Hossein Nezamabadi-pour, Saeid Saryazdi. GSA: A Gravitational Search Algorithm[J]. Information Sciences, 2009, 179(13): 2232-2248.