优化的马尔可夫链人工蜂群算法

郭 佳，马朝斌，张绍博，苗萌萌

1.北京交通大学 计算机与信息技术学院，北京100044

2.国家保密科技测评中心，北京100044

1 引言

人工蜂群算法（artificial bee colony algorithm，ABC）是Karaboga[1]于2005 年提出的一种全局优化算法，具有计算简洁、易于实现、控制参数少等优点[2]，但该算法更新迭代过程中随机选择邻域解，蜜源进化率较低。国内外学者对ABC 算法进行了诸多改进，主要有三方面：一是通过算法融合进行改进，如Panniem 等[3]提出Modified ABC，用萤火虫算法产生新的位置来取代侦察蜂阶段未更新的位置。耿璐等[4]利用差分进化具有良好的局部更新能力，将ABC 与DE（differential evolution）差分算子结合提出adaptive modified differential ABC。Tian 等[5]将元胞自动机与ABC 结合，利用元胞间的自动更新保持解的多样性。这些算法将ABC 算法与其他进化算法进行整体融合，提高了搜索精度，但增加了搜索的复杂度。二是利用局部最优信息引导算法，如受粒子群算法启发，Zhu 等[6]提出GABC（gbeat ABC）算法，用全局最优解gbest代替随机生成的邻域解引导搜索。周新宇等[7]在文献[6]的基础上使用正交实验设计来产生新的食物源，确保侦察蜂在不同维度上保存放弃解和全局最优解的有用信息。Jadon 等[8]在文献[3]和文献[6]的基础上增加了选择因子，以避免过快地收敛到局部最优解。宋晓宇等[9]提出自适应搜索策略混合人工蜂群算法，利用目标函数值的反馈作用进行自适应调节，同时利用全局最优解的引导信息来增强开发能力。这些方法均以解空间上相对最优的某个值作为引导，逐步找到函数最优解，但是这类方法容易使算法过分依赖这个值，从而导致陷入局部最优解。三是针对解空间进行改进，如Karaboga[10]于2014年提出快速人工蜂群算法，引入周边食物源的概念，根据周边食物源参数来判别出周边食物源种群大小，从而提高蜂蜜采集效率。Jadon 等[11]利用全局和局部邻域的差分进化，更新食物源信息。Aydin 等[12]提出了增长式蜂群算法，通过迭代让蜂群规模逐渐增长，以解决大规模持续优化问题。Kumar 等[13]提出了一种两步人工蜂群算法，用K-means 算法代替随机初始化食物源。但这些思想仅改进了有限的局部解信息，没有充分利用整个解空间的发展规律，对于多峰函数来说，较容易陷入局部最优解。

马尔可夫链在启发式算法中已有一些工作，比如文献[14]将马尔可夫链与ABC 算法结合提出MABC（Markov ABC）算法，但该算法利用的二维解空间缺少马尔可夫性的理论推导。Wang 等[15]将遗传算法与马尔可夫链结合，采用遗传算法优化初始参数，以马尔可夫模型作为分类器对大型点焊机加工过程的运行状态进行识别。Yang 等[16]将马尔可夫链用于遗传算法运行前的目标分类，并将改进算法应用于基于蓝牙的交通检测。这些算法均将Markov 链用于进化算法运行前或者后的数据处理阶段，提高了进化算法的寻优精度。

本文在文献[14]的基础上进行改进，将原来的二维解空间扩展至三维，从时间维度组成马尔可夫链进行预测，提出IMABC（improved Markov ABC）算法。与文献[15-16]中马尔可夫链仅作用于进化算法部分数据和一段时间不同，IMABC 算法在空间上将马尔可夫链作用于整个解空间，在时间上将马尔可夫链作用于整个算法运行过程。

2 基本算法2.1 马尔可夫链

马尔可夫（Markov）链是描述动态随机现象的数学模型。它基于系统“状态”和“状态转换”的概念[17-18]，根据概率分布，从一种状态转换为另一种状态[19]，该方法已应用于众多实际领域[20-21]，其定义如下：

在有限或可数个值的随机过程{Xn,n=1,2,…}中，把所取可能值的全体记为状态空间S，用{Xn=i}表示过程在时刻n处于状态i，假设每当过程处于状态i，则在下一个时刻将处于状态j的概率是固定的pij，即对任意时刻n，P(Xn+1=jXn=i)=pij，若对任意状态i1,i2,…,in，始终有P(Xn+1=jXn=in,Xn-1=in-1,…,X1=i1)=P(Xn+1=jXn=in)，这样的离散过程被称为马尔可夫链。假设状态空间是有限集，则称为有限状态马尔可夫链。若对于任意的l，有P(Xt+l=jXl=i)=P(Xt=，即其转移概率不依赖于l时,{Xt,t≥0}称为齐次马尔可夫链[16]，转移概率为：

其中，Li表示状态Si出现的总次数，表示状态Si经n步转移到状态Sj的次数，根据转移概率形成状态转移矩阵，并在下一个时刻进行预测：

其中，Δl是状态区间的下限，Δu是状态区间的上限，f(x)是滚动预测值。

2.2 人工蜂群算法

人工蜂群算法将蜜蜂分为雇佣蜂、跟随蜂和侦察蜂。该算法的步骤如下：

（1）初始阶段。确定蜜源大小SN；蜜蜂的数量与蜜源的数量一致，蜜源Xi=(xi1,xi2,…,xiD)代表候选解，其中i∈{1,2,…,SN}，D是蜜源的维度，蜜源范围是[xmin,xmax]，首先按式（3）生成初始种群：

其中，rand是介于（0,1）之间的随机数，j∈{1,2,…,D}，若xij超出了解空间范围，则按式（4）变为边界值：

求解后使用适应度函数来判断蜜源质量：

（2）雇佣蜂阶段。每次采蜜时在蜜源附近按式（6）随机产生候选解：

其中，xij是搜索空间中的第i个蜜源，xkj是不同于xij的另一个蜜源，k∈{1,2,…,SN}，rand是介于（0，1）之间的随机数。雇佣蜂通过贪婪选择算法决定是否由υij替换xij。

（3）跟随蜂阶段。通过轮盘赌算法，按式（7）计算收益率Pi在整个种群中的概率，确定是否选择该蜜源，选择蜜源后，跟随蜂转变为雇佣蜂。

（4）侦察蜂阶段。当一个蜜源达到开采的上限Limit时，将淘汰蜜源，并根据式（3）随机选择一个新的蜜源。

3 算法的改进3.1 齐次马尔可夫链性质的分析

（1）有限性论证。在初始化阶段，解空间的范围为[xmin,xmax]，每个解的维度D有限，蜂群数量为NP，生成时步为1，随着时步的推移，生成时步为t，每个维度的参数都是有限的，故所有蜂群状态空间有限。

（2）齐次马尔可夫性论证。当雇佣蜂更新候选解的时候，t步概率表达式为：

其中，pxij→vij(t)为经t步由状态xij转换到vij的概率，L为通过贪婪算法对解的选择结果，若fit(xij)≤fit(vij)，L=1，否则等于0。根据式（6），为xij至xkj空间中任意一点的概率。由式（8）可以看出，除了考虑空间的上下界外，解从状态xij转换到vij的概率只与xij和随机选择的值xkj有关，也就是说，个体的状态转移只与前一时刻有关。因此，人工蜂群算法的种群状态序列为有限齐次Markov 链。

3.2 优化的改进算法

IMABC 算法将采蜜过程分为两个阶段，分割参数记为ω。第一阶段为原始ABC 算法，运行ω次后进入第二阶段，从初始状态到第一阶段结束，过程中所有历史解记为Markov 解空间，即由xij扩充为xijt，记为Xt=(xij1,xij2,…,xijω)，将Xt按照当前解生成的时步顺序ω进行排序，最小时步的生成解在最前面，依次形成按时步递进的Markov 链解空间，再将Markov 链解空间划分为N个状态子空间，边界值Qn如式（9）所示：

其中，n∈{1,2,…,N+1}，在确定D和SN的前提下，Xmin为Markov 三维解空间在时间维度的最小值，Xmax为Markov 三维解空间在时间维度的最大值，t∈{1,2,…,ω}。将状态子空间记为Sn，按式（10）所示划分：

将Markov解空间内的解Xt分别置于N个子空间中。按式（1）计算m步状态转移概率，m∈{1,2,…,M}，形成M个状态转移概率矩阵P(m)，选取离预测值最近的M个时步的值进行预测，以离预测时步的远近为顺序，时步小在前，在状态转移矩阵中，选取与m时步相对应的起始状态所在的行向量，组成预测矩阵。计算状态子空间Sn转移到Sn′的总概率，如式（11）所示：

Table 1 State transition prediction表1 状态转移预测

根据表1 计算得出：

pbest所对应的状态空间即为预测解所在的最优状态子空间。获得最优状态子空间后，按式（2）计算新解，代替适应度最差的解。随着t的增加，Xt中的元素逐渐增多，优化过程变慢，需对Xt进行处理，过滤掉初始阶段对预测没有意义的解，形成新的Markov 解空间V，如式（12）所示：

其中，t为迭代次数，λ为丢弃参数，介于(1,lnt)之间，用于控制丢弃比例。在早期阶段，基于初始随机解生成的解空间不能充分代表解空间的发展趋势，应被丢弃，故在时步t较小时，前期解丢弃的比例越大对预测进行的干扰越少。随着算法运行时步t的增加，解空间具有两方面的特点：一是范围迅速扩大；二是随着适用度值逐渐增加，越来越趋于某一值。故如果按前期的丢弃数量进行丢弃，解空间范围没有明显缩小，算法运行效率会越来越低；如果仍按前期的丢弃比例进行丢弃，会加大收敛速度，但在局部最优解较多的函数寻优过程中，会减少算法的开拓能力，导致算法寻优结果为局部最优解。故在Markov 解空间的变化过程中应保证：随着时步t的增加，丢弃解的数量增加，但占总时步比例减小。按式（12），取λ=1，计算对应于不同循环时步的新Markov 解空间范围如表2 所示。

Table 2 Range of solution space changing with time表2 随时步增加而变化的解空间范围

从表2 可以看出，在循环次数较小时，早期解丢弃的数量小，但占总时步比例大；随着循环次数的增加，解丢弃的数量大，但占总时步比例小。在算法运行效率和结果准确性之间达到了一定平衡。

IMABC 的伪代码如下，其中SN为蜜源数量，MaxCycles为最大迭代次数，D为解的维度；ub为解空间的上限，lb为解空间的下限，limit为单个蜜源的开采次数上限，ω为分割参数，N为划分的子空间数量。

IMABC 主函数伪代码如下：

3.3 算法性能分析

3.3.1 算法收敛效率分析

IMABC 算法的迭代分为初期、中期和后期三部分。前期为迭代次数小于分割参数阶段，因IMABC算法依赖于已经产生的解空间，故需先运行ABC 算法，产生用于迭代的初始解空间，此时算法的收敛效率和运行时间与原始的ABC 算法一致，按式（6），算法复杂度为TABC=sn×(d+1)，其中搜索的复杂度是sn×d，计算每个食物源适应度值的复杂度是sn。中期为根据初始解空间进行迭代更新寻找最优值的阶段，此时解空间的维度是sn×d×t，其中t从ω开始。在已有初始解空间的基础上，IMABC 算法依时步的增加不断更新解空间，既保留了已有解的有效信息，又避免了依赖某一最优解导致的算法早熟，可以高效准确地找到最优点。但同时，随着迭代次数t的增加，解空间也随之增加，按式（12）进行处理后解空间维度变为，此时算法复杂度为。后期为t较大的阶段，此时IMABC 算法的搜索时间明显变长，故IMABC算法虽然可以显著提高算法的精度，却要牺牲一定的时间成本。但同时也要考虑，在实际寻优过程中，往往会在算法达到一定精度后触发停止条件，因此IMABC 算法会比ABC 提前达到收敛所需精度，减少IMABC 算法的复杂度和运行时间。

3.3.2 算法寻优能力分析

ABC 算法本身具有的特点就是重勘探轻开采，即具有良好的全局寻优能力，但是局部寻优能力差，IMABC 算法在初期采用原始的ABC 算法，且通过分割参数控制原始ABC 算法在整体算法中的比重，可以有效地保留算法的全局寻优能力。

当逐渐趋于最优值的时候，需要尽量多地保留解空间的优质解信息，但在ABC 算法中，雇佣蜂阶段按式（6）找到新蜜源，由于xkj的随机性，蜂群的先验信息丢失，适应度高的蜜源被放弃，造成算法震荡。但是若按以文献[6]为代表的最优值引导收敛算法，会导致过分依赖gbest，在有多个局部极值点的情况下算法易趋向于局部最小值。IMABC 最大限度地保留了解空间的优质解信息，并根据适应度值预测出最佳子解空间，进一步更新迭代，避免了搜寻过程中的盲目和过分引导，提高了局部寻优能力。

4 实验验证

本章包括4 个实验，第1、2 个实验分别从结果的准确性、收敛效率和运行时间上比较了IMABC 算法、GABC 算法和ABC 算法，第3 个实验比较了不同分割参数对IMABC 算法的影响，第4 个实验比较了不同维度对IMABC 算法的影响。

选择9 种经典测试函数进行实验，不同函数的公共参数设置相同，如表3 所示。

这9 个函数[22-24]中f1、f4、f5、f7和f8为单峰函数，仅有一个极值点，其中f5会在给定的间隔上出现不同的阶跃现象，并伴随产生大量局部极值，f2、f3、f6和f9是典型的非线性多模态函数，在解空间中存在大量的局部极小值点，算法对此类函数进行寻优极易陷入局部最优解。

4.1 结果准确性比较

该实验从结果的准确性上比较IMABC 算法、GABC 算法和ABC 算法，设置维度D=30，分割参数ω=1 000，后续实验将扩展比较D和ω。不同文献给出的SN参考值不同，本文采用文献[1]中的取值30。Limit的大小直接决定IMABC 算法的性能，其值越小，食物源被放弃的可能性越大，算法执行率越高。IMABC 算法的Limit参数有3 种典型值100、200和SN×D，本文采用Limit的最小值100，以实现函数的最高执行率。设置最大循环数MaxCycle=10 000。进行30 次实验，记录每个值并计算平均值、标准差、最佳值和最差值。测试结果如表4 所示。

从表4 可以看出，在本节设定的实验参数前提下，ABC 算法、GABC 算法和IMABC 算法均求出了函数f5的最优解，寻优精度一致。在其他函数寻优过程中，改进的IMABC 算法比原始ABC 算法和GABC 算法的结果精度有明显提高。

4.2 收敛速度比较

该实验从算法收敛效率和运行时间上比较了IMABC 算法、GABC 算法和ABC 算法，设置维度D=30，MaxCycle=500，分割参数ω=100，SN=30 。首先，为每个测试功能设置预设精度，如果在最大循环时间内达到预设精度，则认为该操作成功，否则为失败。进行10 次实验，测试结果列于表5。其中，“最小运行次数”为函数重复寻优过程10 次中达到预设精度值时的最小运行次数，“最大运行次数”为函数重复寻优过程10 次中达到预设精度值时的最大运行次数，“算法成功率”为函数重复寻优过程10 次中达到预设精度的比率，“平均运行次数”为10 次函数重复寻优过程所有运行次数的平均值，其中没有达到预设精度的运行次数不算在内，“/”表示函数重复寻优过程中没有一次达到预设精度，“运行时间”为10 次函数寻优过程运行时间的平均数。

Table 3 Standard test function表3 测试函数

Table 4 Comparison of optimization accuracy results of different algorithms表4 不同算法的优化精度结果比较

Table 5 Comparison of algorithm success rate under given accuracy表5 在给定精度下算法成功率的比较

根据表5 可以看出，在指定预设精度的函数寻优过程中，对于函数f2、f3、f6和f8，算法ABC 和GABC 的成功率为0，而IMABC 达到了100%。对于函数f1和f4，算法ABC 的成功率为0，GABC 有部分达到了收敛精度，而IMABC 的收敛成功率达到了100%。对于函数f7，算法IMABC 比ABC 和GABC算法的成功率高。对于函数f5和f9，3 个算法均在500 次循环内找到了函数的最优值，算法成功率达到了100%，但平均运行次数差不多。因为f5为阶跃函数，函数值计算的过程中采用取底运算，f9的函数值计算过程中采用绝对值运算，故这两个函数均存在多个不同的解对应于一个函数值的情况，也就是不同的解对应相同的适应度值，导致IMABC 算法无法有效地判断Markov 解空间内子空间的优劣，故在迭代次数较少时IMABC算法的收敛效率与ABC相差不多。

表5 最后一列为10 次函数重复寻优运行时间的平均数，从函数f9所用时间可以看出，IMABC 算法的运行时间约为ABC 和GABC 算法的4 倍，可见IMABC 算法以时间换取精度。但因在预设精度的寻优过程中，IMABC 算法达到停止条件所用迭代次数少，所用寻优时间明显降低，f3在IMABC 算法下的寻优时间甚至还要少于ABC 及GABC 算法。

图1 为本次实验中达到预设精度时的函数收敛曲线。

从图1 可以看出除函数f5和f9外，其余7 个函数通过IMABC 算法的收敛效率明显高于ABC 算法和GABC 算法。从图1 还可以看出IMABC 算法的收敛曲线相对平滑，ABC 和GABC 算法的收敛曲线呈阶梯状较明显。因为IMABC 算法的寻优过程是基于整个解空间的迭代更新，函数迭代过程中基本每一次均会更新解，而ABC 算法迭代过程中新解的更新依赖于随机选取的邻域解，存在多次迭代仍未更新的现象，故在函数寻优过程中呈阶梯下降状态。同样，GABC 利用最优解引导寻优过程，会出现一次明显的函数值变化后长期处于不更新的状态。

在本次实验中，函数的寻优过程如图2所示，因图形较多，仅选取单峰函数f1、f7和多峰函数f2。

Fig.1 Curve of function convergence result图1 函数收敛结果曲线

根据图2 的（a）至（c）可以看出，因本次实验分割参数设为100，故在迭代次数小于100 时IMABC 与ABC、GABC算法收敛速度基本一致，但是随着迭代次数的增加，根据图2的（d）至（f）可以看出，在150到200阶段，IMABC 算法已能显出较好的收敛速度，随着迭代次数的进一步增加，就达到如图1所示的收敛结果。

4.3 分割参数的扩展比较

分割参数值ω的不同直接影响函数优化结果，本实验验证不同分割参数对IMABC 算法的影响。设维度D=30，MaxCycle=1 000，SN=10。本实验从ω=100 开始，每增加100 进行一次运算。当ω=1 000 时，表示未运行IMABC 算法，等效于ABC 算法。分别选择单峰测试函数f1，多峰测试函数f2和f3，进行30次实验，记录每个值并计算平均值、标准差、最优值和最差值，测试结果如表6 所示。

从表6 可看出，ω值小于1 000 时的实验结果均优于ω等于1 000 时的性能，进一步证明了IMABC算法的性能优于ABC 算法。

Fig.2 Curve of function convergence process图2 函数收敛过程曲线

Table 6 Comparison of optimization results of different segmentation parameters表6 不同分割参数的优化结果比较

同时可以看出，在函数f1的寻优过程中，最优值出现在ω=500 时，而在f2和f3的寻优过程中，最优值出现在ω=300 时。因为f2和f3均是多峰函数，有多个局部极值点，寻找极值点较困难，故越早使用IMABC 算法越有助于函数寻优。

4.4 不同维度的扩展比较

本实验验证不同维度对IMABC 算法性能的影响。设分割参数ω=100，MaxCycle=5 000，SN=10，维度D从10 开始，每增加20 进行一次运算，增至90。分别选择单峰测试函数f1，多峰测试函数f2和f3，进行30 次实验，记录每个值并计算平均值、标准差、最优值和最差值，测试结果如表7 所示。

Table 7 Comparison of optimization results of different dimensions表7 不同维度的优化结果比较

如表7 所示，随着维度的增加，寻优难度增大，精度逐渐降低。

5 结束语

本文将马尔可夫链与人工蜂群算法融合，提出IMABC 算法。IMABC 算法分为两个阶段，第一阶段运行ABC 算法得出初始解空间，第二阶段利用马尔可夫链对第一阶段产生的解空间进行重构，并进一步预测新解。通过实验验证得出IMABC算法具有较高的精度和收敛速度，但算法运行时间较长。在后续研究中，将进一步完善IMABC 算法，研究最优解状态空间的动态变化，提取有用信息，逐渐缩小解空间，减少马尔可夫链预测过程中消耗的运行步骤，从而进一步减小算法的运行时间。