0 引言

萤火虫算法(Firefly Algorithm, FA)是受萤火虫发光相互吸引和移动的启发而设计的一类新型的群智能优化算法^[1]。已有的仿真结果表明该算法具有精度较高、收敛速度较快、可调参数少、操作简单、易于实现等特点，因而该算法近年来已引起许多学者的广泛关注，并应用到函数优化^[2]、云计算调度^[3]、流水车间调度^[4]、疾病预测^[5]、蛋白质网络^[6]、旅行商问题(Travelling Salesman Problem, TSP)^[7]等诸多领域。

自2008年Yang^[1]首次提出萤火虫算法以来，不少学者从相对吸引力计算、位置更新、自适应步长、算法混合等角度对标准萤火虫算法进行深入研究。例如欧阳喆等^[8]通过引入萤光因子自适应调整萤火虫的步长来平衡算法的全局探索和局部寻优的能力。针对初始种群产生困难、相对吸引力与亮度无关、惯性权重固定、位置更新固定等不足，王吉权等^[9]提出了改进萤火虫算法，该算法使用遗传算法快速生成初始种群，给出了与绝对亮度相关的相对吸引力计算公式，并利用自适应惯性权重动态调整和基于压缩因子的位置更新等方法提高收敛速度，仿真实验结果表明与标准萤火虫算法相比，改进算法具有运算速度快、收敛速度快、进化次数少等优势。针对萤火虫算法存在种群早熟和容易陷入局部最优等不足，王铭波等^[10]提出了一种多种群萤火虫算法。该算法将种群均分为若干子种群，引入动态权重以调整萤火虫的视野范围，并按一定概率对子种群进行局部搜索增强算法跳出局部最优的能力，仿真结果表明该算法在4个标准测试函数上均能找到最优解。针对种群的多样性保持和全局搜索能力等问题，冯艳红等^[11]提出了基于混沌搜索的萤火虫算法，该算法采用立方映射生成混沌序列初始化萤火虫位置，进化过程中使用动态方式执行混沌搜索操作，以一定概率对全局最优解高斯变异以增强算法跳出局部最优的能力，仿真结果表明该算法在6个标准测试函数上较好地平衡了全局搜索能力、收敛速度和精度等。张军丽等^[12]将Powell方法作为局部搜索算子嵌入到萤火虫算法中，提出一种基于Powell方法局部优化的人工萤火虫算法；结果表明，改进后的萤火虫算法在收敛速度、精度和稳定性方面都优于标准萤火虫算法。

尽管萤火虫算法及其改进算法具有不少优点，但仍然存在种群过早收敛、容易陷入局部最优等不足。鉴于此，首先，本文从数学理论上分析并揭示萤火虫陷入局部最优的原因；进而，借鉴Memetic算法的思想，提出一种结合模拟退火算子和萤火虫算法的改进算法。该算法包括两个阶段：1)利用标准萤火虫算法对上一代种群进行全局搜索得到当代种群，进化过程中保持种群多样性和算法全局探索能力；2)对当代种群中的部分最优个体使用模拟退火进行局部搜索，以较大概率接受较好的解，小概率接受较差的解，增强算法的局部开发能力，避免算法陷入局部最优。

值得一提的是，至目前为止，尚未发现与萤火虫算法容易陷入局部最优相关的理论分析。此外，与已有的改进算法^[8-10]不同，本文算法采取吸引过程和退温过程同步进行的算法策略，有效地降低了算法的计算规模和对局部最优解过度优化的风险。

1 优化问题和标准萤火虫算法概述 1.1 优化问题

定义1  不失一般性，优化问题定义为：

$ \begin{array}{l} \min f(\mathit{\boldsymbol{x}})\\ s.t.\;\;\;{x_i} \in [L{b_i}, U{b_i}], i = 1, 2, ..., n \end{array} $

(1)

其中x=(x₁, x₂, …, x_n)为n维决策变量。

1.2 标准萤火虫算法

萤火虫算法是通过模拟萤火虫因发光而相互吸引和移动的现象设计的新兴的群智能优化算法。该算法利用萤火虫发出荧光信号吸引同伴，并朝向区域内发光亮度高的位置移动的机制来实现区域内种群位置的进化寻优^[13]。萤火虫的吸引和移动取决于萤火虫自身的发光亮度I(x)和相对吸引度β：前者由其自身位置的亮度大小所决定，亮度越高表示位置越优，从而越能吸引同伴向它移动；后者决定了萤火虫的移动方向，吸引度越大移动距离越大。

定义2  萤火虫的发光亮度I(x)：

$ I(\mathit{\boldsymbol{x}}) \propto f(\mathit{\boldsymbol{x}}) $

(2)

萤火虫的发光亮度正比于其目标函数值，目标函数值越优，萤火虫的发光亮度越亮。

定义3  萤火虫和萤火虫的相对吸引度β_ij:

$ {\beta _{ij}}{\rm{ = }}{\beta _0}{{\rm{e}}^{ - \gamma {r^2}}}, r = \left\| {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right\| $

(3)

β₀是r=0处的相对吸引度, 通常取1, γ为传播介质的光强吸收系数。式(3)表明了相对吸引度在传播介质随距离的增大而减小。

定义4  萤火虫被萤火虫吸引的位置更新：

$ x_i^{t + 1} = x_i^t + {\beta _0}{{\rm{e}}^{ - \gamma {r_{ij}}^2}}(x_j^t - x_i^t) + {\alpha _t}(rand - 0.5) $

(4)

其中：第一项是萤火虫i当前位置，第二项是被吸引而移动的距离，第三项是随机扰动项，α_t通常为[0, 1]步长因子，rand服从[0, 1]上的均匀分布。

标准萤火虫算法的流程^[2]如下：

1)设定算法基本参数，包括种群规模PopSize，最大进化次数MaxIter，相对吸引度β₀，吸收系数γ，步长因子α_t，目标函数f(x);

2)随机初始化萤火虫种群x_i(i=1, 2, …, n)，计算发光亮度I(x_i)∝f(x_i)，iter=0；

3)开始进入萤火虫种群吸引进化：

①i:=1, j:=1；

②如果f(x_j) < f(x_i)，按照式(3)更新萤火虫i的位置；

③计算新的目标函数值；

④i:=i+1, j:=j+1；

⑤如果i:=i+1, j:=j+1，执行4)，否则重复执行②~⑤；

4) iter:=iter+1；

5)如果iter>MaxIter，算法终止，否则重复执行3)~4)。

2 混合萤火虫Memetic算法 2.1 研究动机

根据式(2)、式(3)可知，萤火虫个体之间的吸引程度随着距离的增加而呈指数衰减，特别地：

$ \left\{ \begin{array}{l} {r_{ij}} \to 0 \Rightarrow {\kern 1pt} {\beta _{ij}} = {\beta _0}{e^{ - \gamma r_{ij}^2}} \to {\beta _0}\\ {r_{ij}} \to \infty \Rightarrow {\beta _{ij}} = {\beta _0}{e^{ - \gamma r_{ij}^2}} \to 0 \end{array} \right. $

(5)

特别的，当β₀=0, γ=1，简单计算可知：

$ \begin{array}{l} {r_{ij}} \in [0, {\rm{2}}] \Rightarrow {\beta _{ij}} \in [1, {\rm{0}}{\rm{.018}}]\\ {r_{ij}} \in [2, 5] \Rightarrow {\beta _{ij}} \in [{\rm{0}}{\rm{.018}}, {\rm{1}}{\rm{.39}} \times {\rm{1}}{{\rm{0}}^{{\rm{-}}11}}] \end{array} $

根据式(5)可知，近距离内的萤火虫个体因相对吸引度较大而靠近，距离较远的个体因吸引系数较小，从而使得移动距离β₀re^-γr_ij2几乎为0。此时：

$ x_i^{t + 1} \approx x_i^t + {\alpha _t}(rand - 0.5) $

(6)

根据式(6)可以看出，x_i^t+1与x_i^t相差随机扰动项α_t(rand-0.5)，这相当于随机搜索解空间。

进一步，对式(6)两边取期望，可知：

$ \begin{array}{l} E(x_i^{t + 1}) \approx E(x_i^t) + E({\alpha _t}(rand - 0.5)) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} E(x_i^t) + {\alpha _t}E(rand - 0.5){\kern 1pt} {\kern 1pt} = E(x_i^t){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \end{array} $

(7)

由式(7)可知，当萤火虫j与萤火虫i的距离达到一定程度时，第t+1代萤火虫的期望位置E(x_i^t+1)和第t代萤火虫的位置E(x_i^t)基本一致，此时算法陷入局部最优。换言之，标准萤火虫算法存在种群早熟、容易陷入局部最优等不足。

因而，针对上述不足，需引入合理的策略，使得算法早期具有一定的全局探索能力，随着进化代数的增加，增强算法的局部搜索能力，同时能以一定概率能够跳出局部最优。尽管模拟退火算法已被证明在具有较强局部搜索能力的同时，能够以一定概率跳出局部最优，但是因其单点搜索特点而使得它的全局优化能力较弱。鉴于此，本文结合萤火虫算法的全局搜索和模拟退火的局部搜索能力，提出模拟退火增强的混合萤火虫Memetic算法。

2.2 模拟退火

模拟退火算法来源于金属固体退火过程，即先将固体加热至高温，然后逐渐冷却，高温时，固体内分子杂乱，处于各个状态的概率均匀分布，低温时，固体内分子逐渐趋于稳定，处于平衡状态的概率最大，从而使内能达到最小值^[15]。模拟退火的核心思想是按照Metropolis准则以一定概率接受较差的解，从而使算法跳出局部最优。

记能量函数为f(S)，当前解S_old的能量为f(S_old)，新解S_new的能量为f(S_new)，能量增量为df=f(S_new)-f(S_old)，则Metropolis准则下接受新解S_new的概率为p：

$ p = \left\{ \begin{array}{l} 1, \;\;\;\;\;\;\;\;\;\;df < 0\\ {{\rm{e}}^{ - df/T}}, \;\;\;\;df \ge 0 \end{array} \right. $

(8)

也就是在温度为时，如果增量df < 0，则以概率1接受新解S_new，否则以概率1-e^-df/T接受当前解S_old。

2.3 FA_SA

Memetic算法是近年来提出的结合全局搜索和局部搜索的一种进化算法设计框架，通过组合不同的全局搜索算子和局部搜索算子可以构成不同的Memetic算法^[14]。其基本框架为：1)确定全局搜索算子和局部搜索算子；2)初始化种群；3)利用全局搜索算子对种群执行全局搜索；4)利用局部搜索算子对个体执行局部搜索；5)根据3)和4)更新种群；6)判断条件，如果满足终止，否则返回2)。

本文提出基于模拟退火增强的混合萤火虫Memetic算法，该算法首先使用萤火虫算法进行全局搜索，其次对排名靠前的部分个体采用模拟退火执行局部搜索，增强算法跳出局部最优的能力。

混合萤火虫Memetic算法流程如下：

1)设定算法参数，包括种群规模PopSize，最大进化次数MaxIter，相对吸引度β₀，吸收系数γ，步长因子α_t，目标函数f(x)，初温T=T₀，链长L，退火因子α，执行模拟退火的萤火虫个数k。

2)随机初始化萤火虫种群x_i(i=1, 2, …, n)，计算发光亮度I(x_i)∝f(x_i)，iter=0。

3)开始进入萤火虫种群吸引进化：

①i:=1, j:=1；

②如果f(x_j) < f(x_i)，按照式(3)更新萤火虫i的位置；

③计算新的目标函数值；

④i:=i+1, j:=j+1;

⑤如果i>PopSize, j>PopSize，执行4)，否则重复执行②~④。

4)对排名并找到前k个萤火虫个体，使用模拟退火算子执行局部搜索：

①产生每个个体X_k的新解X_k^new;

②计算增量df_k=f(X_k^new) -f(X_k)；

③按照式(8)以一定概率接受新解X_k^new；

④重复①~③共L次，结束模拟退火局部搜索。

5) iter:=iter+1，T:=αT；

6)如果iter>MaxIter，算法终止；否则重复执行3)~5)。

值得一提的是，与文献[10, 12]将萤火虫算法与局部搜索算法简单结合不同，本文提出的混合算法同步进行萤火虫的吸引过程与模拟退火的退温过程，两者共用一个外层循环，这样有助于降低算法的计算复杂度和避免对局部最优解的过度优化。

3 仿真实验与分析 3.1 测试函数

为验证FA_SA的寻优性能，本文对10个标准测试函数^{[10, 15-16]}进行测试，与萤火虫算法(FA)、粒子群优化(Particle Swarm Optimization, PSO)、蝙蝠算法(Bat Algorithm, BA)、花朵授粉算法(Flower Pollination Algorithm, FPA)、改进花朵授粉算法(Simulated Flower Pollination Algorithm, SFPA)、人工蜂群算法(Artificial Bee Colony, ABC)、多种群萤火虫算法(Simulated-based Multi-group Firefly Algorithm, MFA_SA)进行对比。

F1~F2为单峰函数，具有一个全局极值点；F3~F6是多峰函数，具有多个局部极值点；F7~F10是F1~F3和F5经过平移和旋转后的多峰函数。

标准测试函数的信息和形态参见表 1和图 2。

3.2 实验环境

实验环境：Window 7，CPU奔腾T4400，主频2.2 GHz, 内存为2 GB，编程语言为Matlab 2014b。

3.3 参数设置

为保证结果的公平性及客观性，标准萤火虫算法(FA)和混合萤火虫Memetic算法(FA_SA)采取相同参数，其算法参数设置来源于文献[2]：

1) FA和FA_SA参数：种群规模，最大迭代次数，相对吸引度，吸收系数分别为20、300、1、1；

2)模拟退火参数：初温T₀=100, 链长L=40, 退火因子α=0.5，执行模拟退火的萤火虫个体数k=1。

3) FA和FA_SA分别独立运行20次，函数评价次数达到6000次终止进化；

4) PSO、BA、FPA、SFPA等4种算法在F1、F3~F6上的结果来源于文献[15]，种群规模为20，最大迭代次数为500，最大函数评价次数为20×500=10000；

5) ABC、PSO、MFA_SA等3种算法在F2上的结果来源于文献[2], 种群规模为100，进化次数为250，最大函数评价次数为100×250=25000。

3.4 结果分析

为了全面地评价FA_SA的性能，本文分别比较了PSO、BA、FPA、SFPA、FA、FA_SA等6种算法在F1, F3~F6上的寻优性能，ABC、PSO、MFA_SA、FA、FA_SA等5种算法在F2上的寻优性能以及FA、FA_SA在F7~F10上的寻优性能。

3.4.1 单峰函数寻优性能分析

表 2汇总了不同算法在单峰函数F1~F2上的最优值、最差值、平均值、标准差等寻优统计结果。

从表 2中可以看出，对于Sphere球形函数，相对于表现最好的SFPA，FA_SA使用60%的计算代价(6000/10000=60%)，找到了排名第2的最优值，并且方差也比其他3种算法要小，表现出更强的鲁棒性和适应性。

Rosenbrock的全局最优点位于一个平滑、狭长的抛物线型山谷地带，由于其方向较难确定，难以达到全局最优值，但相对于ABC、PSO、MFA_SA、FA等4种算法，FA_SA使用24%的计算代价(6000/25000=24%)，找到了排名第2的较优解，其最优值、平均值、最差值和标准差均位于第二位。

特别地，相对于MFA_SA，FA_SA结果无论是从数量级还是鲁棒性上都与MFA_SA一致，可以认为增加计算资源FA_SA将找到更优化的解。

3.4.2 多峰函数寻优性能分析

表 3汇总了不同算法在多峰函数F3~F6上的最优值、最差值、平均值、标准差等寻优统计结果。从表 3可知，对于Rastrigin、Ackley、Shaffer三个函数问题，与PSO、BA、FPA、FA等4种算法相比，FA_SA表现了更强的寻优性能：使用60%的计算代价(6000/10000=60%)，找到了更优的最优值，并且方差也比其他算法要小，表现出更强的鲁棒性和适应性。特别地，FA_SA在Ackley、Shaffer问题上3000次左右的评价时就达到了较好的寻优精度(参见图 3(d)与图 3(f))。

因此，与其他算法相比，尤其是标准FA，FA_SA在寻优精度、收敛速度和鲁棒性方面有较优的性能，同时求解多峰极值优化问题时，FA_SA因为引入模拟退火算子执行局部搜索时具备跳出局部最优的能力，利用Memetic机制可有效平衡搜索过程中的全局探索和局部开发能力。

3.4.3 进化收敛速度分析

为了更加清晰、直观地对比FA和FA_SA的寻优精度和收敛速度，图 3是FA和FA_SA在6个标准测试函数下的进化收敛曲线。

从图 3可以看出：相对于FA，FA_SA因综合了全局搜索和局部搜索的优点，表现出更优的寻优精度；FA_SA在3000~4000次评价后基本能找到最优解，收敛速度优于FA；FA在进化一段时间后收敛曲线基本不变，FA在搜索过程中陷入了局部最优，而FA_SA因引入了模拟退火操作使算法能够有效地跳出局部最优从而继续寻找更优解，这在Shaffer函数上的表现尤为明显。

3.4.4 复合多峰函数寻优性能分析

为进一步考察FA_SA对复合多峰函数优化的适用性，使用FA和FA_SA对文献[16]中的经平移和旋转变换后的多峰函数进行对比测试，表 4汇总了FA与FA_SA在复合多峰函数中的寻优结果。

表 4的结果表明，对于经过平移和旋转变换后的复合多峰函数F7~F10，FA_SA的最优值、平均值、最差值和标准差等寻优效果均优于FA。FS_SA在Shifted Sphere和Shifted Rastrigin十分接近最优值，在Shifted Griewank达到最优值。

由于Rosenbrock的全局最优点位于一个平滑、狭长的抛物线型山谷地带，优化经平移和旋转变换后Shifted Rosenbrock更加困难。尽管FA_SA的寻优结果离最有结果有一定偏差，但是其最优值、平均值、最差值和标准差等均优于FA。

3.4.5 种群规模对FA_SA性能的影响

为探讨种群大小对FA_SA寻优性能的影响，以F2、F3、F7和F9作为三类函数的典型代表，进一步研究不同种群规模对FA_SA的寻优性能规律。该组实验种群大小依次为30、40、50，函数最大评价次数为6000，每组实验独立运行200次。

从表 5的结果可知，随着种群规模的增大，FA和FA_SA的寻优性能均有不同程度下降，但考虑到不同种群下的计算资源为6000次，而这种寻优性能差异均没有数量级的下降，因而可以认为种群规模对FA_SA的寻优性能影响较小。

从表 5可以看出，在相同的种群规模和同样评价次数时，除了少数几种情况，FA_SA寻优的平均值和标准差均优于FA。这些结果表明了FA_SA寻优性能的鲁棒性和最优性均优于标准的萤火虫算法FA。

3.4.6 模拟退火的萤火虫个数K对FA_SA性能的影响

为探讨执行模拟退火的萤火虫个数K对FA_SA寻优性能的影响，以F2、F3、F7和F9作为三类函数的典型代表，进一步研究参数K对FA_SA的寻优性能规律。该组实验种群大小依次为0、1、5、10、15、20，函数最大评价次数为6000，每组实验独立运行200次。

从表 6的汇总结果可知，随着执行模拟退火萤火虫个数K的增大，除Shifted Sphere以外，FA_SA的寻优性能整体上有一定的程度的改进，但改进程度不是很大，因而可以看出使用模拟退火对过多的个体进行局部搜索不完全能使得搜索性能得到很大改善，甚至可能会下降，这是因为模拟退火算子以一定概率接收劣解，太多的模拟退火操作可能增加搜索非最优解的可能性。因而，实际操作中推荐执行模拟退火操作的萤火虫个数K设置为1~5。

4 结语

针对萤火虫算法存在过早收敛、容易陷入局部最优等不足，本文提出了基于模拟退火算子的混合萤火虫Memetic算法。该算法使用萤火虫算法执行全局搜索，使用模拟退火对当前种群中的部分个体执行局部搜索，增强算法的跳出局部最优的能力。对10个标准测试函数的对比实验表明，与其他算法尤其是FA相比，FA_SA在3个标准测试函数中找到最优解，在其他3个测试函数中找到了与最优解数量级和鲁棒性一致的次优解，在寻优精度、收敛速度、鲁棒性方面表现出较优的性能，且算法能够有效地跳出局部最优，在4个经平移和旋转变换后的复合多峰函数上效果均优于萤火虫算法。此外，不同种群规模对FA_SA性能的影响的结果表明了FA_SA寻优性能的鲁棒性和较优性。

进一步的研究工作包括：1)从FA的位置更新方程出发，深入研究FA陷入局部最优的机制，据此改进并设计更有效的吸引算子；2)研究算法参数对FA_SA的寻优性能影响；3)研究FA与其他局部搜索算子混合的Memetic算法^[17]；4)研究FA_SA在工程优化、车间调度^[18]等问题中的应用。