万有引力算法(GSA)在2009年首先被Rashedi等人提出来^[1]，这种算法是基于宇宙间的万有引力定律和运动定律而得出的.对于许多复杂的优化问题如目标控制、模式识别、图像处理以及目标函数优化等万有引力算法都能很好地加以解决.

在GSA的进化过程中，物质向着当前最优值的附近移动，出现了停滞和收敛速度慢等现象. Nihan Kazak和Alpaslan Duysak^[2]在研究的过程中，针对标准引力算法收敛速度慢的缺点，给出了一种成员卫星算法，物质群每个物质称为一个成员个体，每个成员都有各自的若干个卫星，成员和卫星一起进化，保证了算法在当前最优值附近进行更加精细的搜索.刘勇和马良^[3]提出了一种解决非线性极大极小问题的混沌万有引力算法，利用混沌优化对当前的最优位置进行精细搜索, 为解决早熟收敛现象提供了一个方向.而针对粒子群算法(PSO)的早熟现象以及收敛速度慢的缺点，Natsuki Higashi和Hitoshi Iba^[4]提出了一种基于高斯变异的粒子群算法，将应用在遗传算法(GA)的高斯变异引入PSO，提高了PSO算法的性能.在此基础上，Tang Jun和Zhao Xiao-juan^[5]提出了一种带有自适应局部搜索算子的PSO，在进化过程中围绕着全局最优个体附近进行局部搜索，提高了最优值附近的寻优能力.Sarafrazi, Nezamabadi-pour以及Saryazdi^[6]受天文学中“分裂”思想的启发，提出了一种具有分裂功能的万有引力算法，将满足一定条件的物质分裂，避免了算法在进化过程中出现停滞的状况.

变异是一种跳出局部最优解的有效算子，并且对最有物质变异操作对万有引力算法的效率有重大影响，为此本文提出了一种围绕最优物质变异的万有引力算法，在保持和提高问题解的质量的同时提高引力算法的搜索能力和收敛速度.该算法的主要思想是在进化的过程中在最优物质的附近进行变异.如果变异产生出来的新物质优于当前的最优物质，则替换原有的最优物质继续进化.对8个常用的测试函数的实验结果表明，与标准的引力算法(SGSA)和混沌引力算法(CGSA)的比较，新算法有较强的搜索能力和较高的稳定性.

1 基于变异的万有引力算法

优化问题的模型为minf(x), 其中，x∈Rⁿ, R为实数域, x =(x₁, x₂, …, x_n)^T, x_i∈[a_i, b_i](i=1, 2, …，n).

1.1 万有引力算法简介

万有引力算法是根据宇宙间的引力法则和物质间的相互作用的启发而提出来的新的优化算法.根据牛顿定律任何两个物质之间都是相互吸引的，正是由于引力的存在使得物质向着最优解附近运动，在这种算法里，质量比较大的物质比质量小的物质运动得慢，从而相互作用的物质之间都朝着质量大的运动.物质的位置就是问题的解, 当算法满足结束的条件时，优化问题也就达到了最优解.

假设每个物质被定义在一个n维的搜索空间里，由N个物质组成的一个物质群为X ={x₁, x₂, …, x_n}，其中第i个物质的位置是X_i=(x_i¹, x_i², …, x_i^d, …, x_iⁿ)^T, x_i^d为物质i在空间第d维的位置，i=1, 2, …, N.

万有引力常量定义如下：

$ G\left( t \right) = {G_0}{{\rm{e}}^{\left( { - \alpha \frac{t}{T}} \right)}}, $

(1)

其中，G₀为一个初始值，α为一个常数，t为当前代数，T为最大迭代次数.

当前物质的质量是由适应值计算出来的，对于极小值问题，有

$ \begin{array}{l} {\rm{best}}\left( t \right) = {\min _{j = \left\{ {1, \cdots ,N} \right\}}}{\rm{fi}}{{\rm{t}}_j}\left( t \right),\\ {\rm{worst}}\left( t \right) = {\max _{j = \left\{ {1, \cdots ,N} \right\}}}{\rm{fi}}{{\rm{t}}_j}\left( t \right). \end{array} $

(2)

其中best(t)，worst(t)分别表示在第t代物质群中最好和最坏的适应值. fit_j(t)表示物质j在第t代的适应值，每一个物质在第t代的质量为

$ \begin{array}{l} {M_{ai}} = {M_{pi}} = {M_{ii}} = {M_i},i = 1,2, \cdots ,N.\\ {m_i}\left( t \right) = \frac{{{\rm{fi}}{{\rm{t}}_i}\left( t \right) - {\rm{worst}}\left( t \right)}}{{{\rm{best}}\left( t \right) - {\rm{worst}}\left( t \right)}}, \end{array} $

(3)

$ {M_i}\left( t \right) = \frac{{{m_i}\left( t \right)}}{{\sum\nolimits_{j = 1}^N {{m_j}\left( t \right)} }}. $

(4)

在第t代，第d维定义物质j对物质i作用的万有引力是

$ F_{ij}^d\left( t \right) = G\left( t \right)\frac{{{M_{pi}}\left( t \right) \times {M_{aj}}\left( t \right)}}{{{R_{ij}}\left( t \right) + \varepsilon }}\left( {X_j^d\left( t \right) - X_i^d\left( t \right)} \right), $

(5)

其中，R_ij(t)为物质i和物质j之间的欧氏距离，ε为一个很小的值.物质i在第d维所受到的合力为

$ F_i^d\left( t \right) = \sum\limits_{j \in k{\rm{best}},j \ne i}^N {{\rm{ran}}{{\rm{d}}_j}F_{ij}^d\left( t \right)} , $

(6)

rand_j为介于0和1之间的随机数，kbest从k₀到1随着代数t线性地减小，使得算法到最后在搜索空间中只有最好的解的那个物质去作用于其他的物质.

由牛顿运动法则，在第t代物质i在第d维的加速度为

$ a_i^d\left( t \right) = \frac{{F_i^d\left( t \right)}}{{{M_{ii}}\left( t \right)}}, $

(7)

物质的速度更新为

$ \left. \begin{array}{l} V_i^d\left( {t + 1} \right) = {\rm{ran}}{{\rm{d}}_i}V_i^d\left( t \right) + a_i^d\left( t \right)\\ X_i^d\left( {t + 1} \right) = X_i^d\left( t \right) + V_i^d\left( {t + 1} \right) \end{array} \right\}, $

(8)

其中，rand_j为0和1之间的随机数.当终止条件达到以后，算法结束.

1.2 基于最优变异的引力算法

考虑到标准的引力算法在进化的过程中，系统中的物质都向着当前最优物质的位置移动，如果算法达到循环结束条件后输出的解并不是所期望的最优解，则算法陷入了局部最优，一旦陷入了局部最优就很难使其跳出当前的状态^[7-9].针对该缺点，本文给出了一种变异算子，该算子根据当前的最优物质进行变异，赋予了算法更强的搜索能力^[10-15].

1.2.1 变异算子

$ {\rm{gbes}}{{\rm{t}}_k}\left( {t + 1} \right) = {\rm{gbes}}{{\rm{t}}_k}\left( t \right) + 0.5rn\;{\rm{gbes}}{{\rm{t}}_k}\left( t \right), $

(9)

其中，gbest_k(t)为第t代全局最优物质在第k维的位置，gbest_k(t+1)为经过上式变异后在第k维新生成的位置.高斯分布的概率密度函数为：$ f\left( x \right) = \frac{1}{{\sqrt {2\pi } \sigma }} \times {{\rm{e}}^{{\rm{-}}\frac{{\left( {x-\mu } \right)2}}{{2{\sigma ^2}}}}}$, 其中σ和μ为密度函数的参数，rn为均值为0方差为1的高斯分布的随机数. 图 1为正态分布的一个图像.

根据该图像的意义，遵从正态分布的随机变量的概率规律为取0附近的值的概率大，而取距离0越远的值的概率越小.由式(9)可以看出，变异后的取值在当前最优位置上下浮动，并且兼顾了算法所要考虑的探索能力和开发能力.

1.2.2 约束处理

考虑到这种变异算子可能产生一些越界的物质，本文将采用下面的方式进行处理，假设变异后的物质i的位置是X_i′=(x_i¹, x_i², …, x_i^d, …, x_iⁿ)^T, x_i^d代表物质i在空间第d维的位置，i=1, 2, …, N.则有

$ \begin{array}{l} {\rm{if}}\left( {x_i^d < {x^d}\_\min } \right){\rm{or}}\left( {x_i^d > {x^d}\_\max } \right),\\ x_i^d = {x^d}\_\min + {\rm{rand}} * \left( {{x^d}\_\max - {x^d}\_\min } \right),\\ {\rm{end}}. \end{array} $

其中x^d_max和x^d_min为最优物质在第d维的上下界，rand为介于0和1之间的随机数.

本文算法的主要步骤的流程图如图 2所示.

2 数值实验

在实验过程中，为了验证本文算法的效果，这里用了8个最小值优化函数来比较3种算法的效果，这8个测试函数中f₁～f₄为单模测试函数，f₅～f₈为多模测试函数，维数都取30维.测试函数如表 1所示.

3种算法分别是标准引力算法(SGSA)^[1]、混沌引力算法(CGSA)^[3]以及本文设计的改进的引力算法(OMGSA).实验中，种群规模N取100，最大迭代次数为2 000，G₀取100，k₀取100.每个算法独立运行30次，得出了循环最后的最优值(BSF)，最优值的平均值(ABSF)以及最优值的标准差(STDBSF)，如表 2所示.

从表 2中可以看出，在函数f₁～f₈中，OMGSA的最优值以及最优值的平均值都比SGSA和CGSA好，对于4个单模测试函数，OMGSA明显优于CGSA和SGSA，对于f₁, CGSA只是略优于SGSA，而OMGSA将标准遗传算法的最优值提高了多个数量级，无论是对于BSF、ABSF还是STDBSF都是如此.相似的情况也可以从其他的3个单模函数的表现得到体现.f₅作为一个容易发生搜索越界的多模函数，算法也会很容易陷入停滞状态，而OMGSA很好地解决了这个问题，得到的最优解明显优于其他两个算法.对于f₆本文的算法在搜索后输出的最优解可以达到已知的全局最优值0.另外两个多模函数从表中的数据也可以看出本文算法无论在求解精度和稳定性方面都表现出了很好的性能.表 2的数值结果表明，OMGSA明显优于另外的两个算法.

图 3是各算法的进化曲线.通过收敛曲线图可以发现，OMGSA不仅具有良好的求解精度, 也具有很好的收敛速度.对于f₁、f₂和f₄的进化过程中发现OMGSA在进化的前期就能很快地收敛到最优解，由于单峰函数只有一个极值点，物质之间相互吸引，加上一个变异算子使得算法可以在当前最优解附近探索到更好的解，从而引导物质群向更好的位置进化.对于f₅～f₈多峰函数，从f₅和f₈可以看出，没有这种变异算子的标准引力算法在进化的中后期陷入了停滞状态, 无法找到周围更好的解，而加上变异算子后算法在进化过程中总能够探索出更好的解, 使得算法不断地进化.f₆和f₇的收敛速度比较其他两个都有很大的提升.总体上可以看出改进的算法突破了局部约束探索到周围更好的解.相比于其他的两个算法，本文的算法能够以更短的代数达到精确的解.

由f₄、f₆～f₈的进化曲线可以看出，本文的算法在进化的前期就以很快的速度达到最优值.由f₃和f₅的进化曲线可以看出，本文算法在进化过程中可以突破出局部极值，向最优值进化.

3 结论

本文提出了一种围绕最优个体附近进行变异的万有引力算法，通过一系列的优化函数对其性能进行了测试.测试的结果表明改进后的OMGSA算法求解精度更高，收敛速度更快，是一种稳健的算法.