万有引力算法(GSA)^[1]在2009年首次由Rashedi等人提出，其思想是基于物理学中宇宙间的万有引力定律和运动规律而得出的.基于万有引力算法的优化算法得到广泛关注.受“分裂”思想的启发, S.Sarafrazi等^[2]提出一种具有分裂功能的万有引力算法.徐怀祥等^[3]引入一个变异算子使算法跳出局部最优，给出了一种改进的万有引力算法.谷文祥等^[4]将万有引力算法运用于解决流水线的调度问题.为解决非线性的极大极小值问题，马良等^[5]提出一种新的混沌万有引力算法，对当前的最优位置利用混沌优化进行比较细密的搜索. Nihan Kazak和Alpaslan Duysak针对标准引力算法的收敛速度比较慢的缺点，提出了成员卫星算法^[6]，每个物质都有各自的几个卫星，卫星和成员的同步进化，保证了此算法在当代最优值周围的更为精细的搜索.为提高引力算法的局部收敛能力，一些学者将引力算法与其他算法结合，提出了许多的混合引力算法^[7-10]. Zahir^[11]提出的模糊万有引力算法，给出了数据挖掘的新途径.许多学者通过引入各种变异算子^[12-15]，增加物质种群的多样性，提高算法的搜索效率.

上述算法都存在收敛速度慢, 求解精度低等问题.针对这些问题, 本文提出了一种根据平均粒距大小来判断算法收敛程度, 对物质群中的最优位置和最差位置进行自适应混沌变异的万有引力算法，保持和提高函数求解精度，同时提高算法的收敛速度.该算法的主要思想是在进化的过程中引入平均粒距，如果平均粒距大于给定值，说明算法处于初期，则对最差物质的位置进行混沌变异，加快算法收敛速度.如果平均粒距小于或者等于给定值，即算法处于中后期，陷入局部最优时，则在最优物质的附近进行混沌变异，增加物质种群的多样性，跳出局部收敛.如果变异产生出来的新的物质优于当前的最优物质，则替换原有的最优物质继续进化.测试函数的实验结果表明了该算法的有效性.

在GSA算法模型中，优化问题的每个解与搜索空间中的一个物质相对应，每个物质都有4个特点：位置、施力物、受力物、惯性质量，而物质的位置就是优化问题的解.标准GSA的算法过程如下：

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

首先利用式(1) 和式(2) 计算出每个物质的惯性质量M_i(t).

其中value_i(t)表示物质i在第t代的适应值，x_best(t)(t)，x_worst(t)分别表示第t代种群中最好和最坏的适应值.

万有引力常量定义为

其中，G₀是一个初始值，α为一个常数，t是当前代数，T是最大的迭代次数.对于求最小值问题x_best(t), x_worst(t)的定义为

其次，计算各个物质在每一维间上相互之间的引力.物质j对物质i作用的第d维万有引力定义为

其中，ε是一个极小的值，R_ij(t)表示物质i和物质j之间的欧氏距离，由此可得出物质i在第d维所受到的合力为

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

然后根据牛顿运动规律，计算出在时间t时物质i在第d维的加速度为

最后物质的速度和位置更新为

其中，rand_i为0和1之间的随机数，它可使搜索变得随机化.

由于在标准GSA算法进程中，物质种群是根据追随种群最优位置进行更新的，所以这种算法在运行过程中，如果某物质找到了一个当前最优位置，则其他物质会迅速向这个最优位置聚集.如果这个最优位置是一个局部最优点，万有引力算法就无法在解空间内重新搜索，出现停滞或早熟收敛现象.

试验证明，GSA中无论是全局收敛还是早熟收敛，万有引力算法中的物质都会有“聚集”的现象，此时非常缺乏种群的多样性.种群多样性的缺乏是导致过早收敛以及停滞现象的重要原因.

本文采用种群的平均粒距来描述种群的多样性.设L为搜索空间对角最大长度, n为解空间维数, N为种群规模大小，P_d表示所有物质第d维坐标值均值, x_i^d(t)表示第i个物质的第d维坐标值, 则本文定义第t迭代时种群的平均粒距D(t)为

平均粒距说明了在单位搜索空间上物质的分布情况，描述了种群中各个物质间彼此相互分布的离散程度，D(t)越小，表示种群越集中.即可能陷入局部收敛.

混沌指在一个确定的系统中发生的不规则的运动，行为主要表现为随机性、规律性、遍历性.由于执行简单并且能够有效避免陷入局部极值，因此混沌搜索已成为一种新的局部搜索技术.

本文引入Logistic混沌方程，其表达式为

由于GSA算法本身的局限性使得标准的引力算法比较容易陷入局部最优，而混沌的遍历性、不重复性，可以在一定的搜索空间内跳出局部收敛，因此在GSA算法中引入混沌搜索，获得全局最优解就成为可能.自适应混沌变异GSA算法的主要思想是在标准GSA算法的基础上，引入平均粒矩，当平均粒距小于或等于给定值时，算法陷入早熟收敛状态, 对最优物质进行混沌搜索，引导物质快速跳出局部最优，避免陷入早熟收敛.当物质的平均粒距大于给定值时，对最差物质进行混沌搜索，产生一个优于最差物质的位置, 来替换最差物质，加快算法收敛速度.

当算法的平均粒距大于给定值d₀时，即在算法初期, 对最差物质的位置由下式变异更新.

其中，x_max, x_min为个体的上下界，x_best(t)为第t代时，当前最差位置.α=4.0.x_j(t)为混沌变量.x′_worst(t)为最差位置变异后所得的新的位置.

当算法的平均粒距小于或等于给定值d₀时，此时算法出现早熟收敛，对最优物质的位置由下式变异更新.

其中，x_max, x_min为个体的上下界，x_best(t)为第t代时，当前最优位置.α=4.0, x_j(t)为混沌变量.x′_worst(t)为最优位置变异后所得的新的位置.

考虑到本文中的变异方法可能产生一些越界的物质，对越界物质采用边界变异的方式进行处理.

if x(p) < x_min

 x(p)=x_min

endif

if x(p) > x_max

 x(p)=x_max

endif

其中x_max和x_min为物质在第p维位置的上下界，x(p)为变异后物质第p维的位置.

本文算法主要步骤如下.

步骤1：设定算法参数, 包括种群平均粒距给定值d₀等, 在解空间随机初始化N个物质的位置和速度.

步骤2：计算物质的函数适应值value_i(t).

步骤3：利用式(3)~(4) 计算更新每个物质的x_worst(t)和x_best(t).

步骤4：利用式(9) 计算群体的平均粒距D(t).

步骤5：若D(t)＞d₀，利用式(12) 对最差物质进行混沌变异得到x′_worst(t)，计算的x′_worst(t)适应值如果优于x′_worst(t)，替换最差物质以及计算种群的最优值和最差值，转入步骤7.

步骤6：若D(t)≤d₀，利用式(11) 对最优物质进行混沌变异得到x′_best(t)，计算的x′_best(t)适应值如果优于x_best(t)，x_best(t)替换最优物质以及最优值, 转入步骤7.

步骤7：利用式(1)、(2) 计算每个物质的质量M_i(t).

步骤8：根据式(5)、(6) 计算每个物质在不同维的合力F_i^d(t).

步骤9：根据式(7)、(8) 计算物质的加速度，速度以及更新种群中物质的位置.

步骤10：如果达到了结束的条件, 最大迭代次数或者更新的位置足够好，则迭代结束，否则返回步骤2.

为验证改进后算法的有效性，本文选择了5个标准的测试函数与标准万有引力算法(SGSA)和混沌万有引力算法(CGSA)作比较，如表 1所示.本文改进后的算法简称ACGSA.5个测试函数中f₁~f₃为单模测试函数, f₄~f₅为多模测试函数.在实验中, 种群规模N=30, 最大迭代次数T=50, 平均粒距给定值d₀=0.003, 万有引力常量初值G₀=100, 最大速度v_max=[3, 3, 3], 最小速度v_min=[－3, －3, －3], K₀=100.表 2为3种算法独立运行50次, 得到的最优值(BSF), 最优值的平均值(ABSF)以及最优值的标准方差(STDBSF)的比较结果.

表 1(Table 1)

表 1 测试函数 Table 1 Test function 函数 表达式 最优值 区间 f1 $ {\rm{min}}\left[{\mathit{f}\left( \mathit{x} \right)} \right]{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{n} {\mathit{x}_\mathit{i}^2} $ 0 [-100, 100] f2 $ {\rm{min}}\left[{\mathit{f}\left( \mathit{x} \right)} \right]{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{n} {{{\left( {\sum\limits_{\mathit{j}{\rm{ = 1}}}^\mathit{i} {{\mathit{x}_\mathit{j}}} } \right)}^{\rm{2}}}} $ 0 [65.536, 65.536] f3 $ {\rm{min}}\left[{\mathit{f}\left( \mathit{x} \right)} \right]{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{n} {{{\left( {1 - {\mathit{x}_\mathit{i}}} \right)}^{\rm{2}}}} $ 0 [-2, 2] f4 $ {\rm{min}}\left[{\mathit{f}\left( \mathit{x} \right)} \right]{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{n} {\left( {\begin{array}{*{20}{c}}{\mathit{x}_\mathit{i}^2 - 100{\rm{cos(2 \mathsf{ π} }}{\mathit{x}_\mathit{i}}{\rm{) + 10}}}\end{array}} \right)} $ 0 [-5.12, 5.12] f5 $ \rm{min}\left[ \mathit{f}\left( \mathit{x} \right) \right]\rm{=}\frac{1}{4\ 000}\sum\limits_{\mathit{i}\rm{=1}}^{\mathit{n}}{\mathit{x}_{\mathit{i}}^{2}}-\prod\limits_{\mathit{i}\rm{=1}}^{\mathit{n}}{\rm{cos}\left( \frac{{{\mathit{x}}_{\mathit{i}}}}{\sqrt{\mathit{i}}} \right)+1}$ 0 [-600, 600]

表 1 测试函数 Table 1 Test function

表 2(Table 2)

表 2 计算结果 Table 2 Computational results 函数 SGSA CGSA ACGSA f1 BSF 7.23×10-3 9.16×10-3 3.53×10-5 ABSF 4.64×10-5 1.61×10-5 6.05×10-10 STDBSF 1.53×10-9 1.80×10-4 1.17×10-12 f2 BSF 2.24×10-3 2.94×10-4 1.00×10-7 ABSF 1.06×10-7 2.20×10-6 1.01×10-11 STDBSF 7.17×10-7 2.33×10-7 5.52×10-16 f3 BSF 1.31×10-7 1.78×10-7 3.13×10-13 ABSF 1.40×10-3 2.77×10-4 1.03×10-7 STDBSF 3.89×10-7 8.55×10-7 2.64×1014 f4 BSF 0.48×10-1 0.14×10-1 2.52×10-5 ABSF 1.16×10-5 2.59×10-5 1.98×10-9 STDBSF 1.93×10-5 1.79×10-4 3.13×10-12 f5 BSF 2.59×10-1 1.21×10-1 1.94×10-4 ABSF 3.46×10-3 2.95×10-3 2.77×10-9 STDBSF 2.89×10-3 5.04×10-3 2.25×10-7

表 2 计算结果 Table 2 Computational results

对于单峰函数f₁~f₃, 由图 1可以发现在这3种算法中ACGSA均能准确且快速地找到函数的最优解.同时从表 2的数据结果来看, ACGSA算法搜索到的最优值精度都比SGSA和CGSA高.由标准方差来看, ACGSA算法的稳定性也是最好的.

图 1(Figure 1)

图 1 f1~f5的进化曲线 Figure 1 evolutionary curve of f1~f5

对于多峰函数f₄~f₅, 其他两种算法都快速地陷入局部最优, 出现停滞现象，而ACGSA由于引入混沌搜索, 算法的局部和全局搜索能力的调整更为灵活，提高了搜索过程的效率，又增加了种群的多样性, 使算法继续在解空间搜索，逐步收敛到精度更高的全局最优值.同时表 2中3种算法的标准方差的比较结果也反映出该算法在处理多峰问题时具有较好的稳定性.这说明随着问题复杂程度的增加, ACGSA的优化效果没有减弱.

本文给出一种自适应的混沌变异万有引力优化算法(ACGSA算法).该算法在运行的前期可以经过变异，加快算法收敛速度；在算法运行的后期可以跳出局部收敛，增加种群的多样性，并在全局最优区域进行较为细密的搜索.本文算法搜索速度和收敛精度比标准GSA算法及CGSA算法有较大提高，并且能够以更快的速度找到全局最优解.