基于种群个体数自适应的多尺度量子谐振子优化算法

焦育威 王 鹏

伴随人工智能的发展,大量的启发式算法相继提出,如遗传算法、粒子群优化算法等,这些算法在优化领域中有着出色的表现.随着优化问题复杂程度的加深,这些算法的弊端开始显现,如算法“早熟”或种群多样性丧失导致精度下降、计算成本上升等问题,为解决这些问题近年来大量的改进策略被提出,这些改进策略可以划分为两大类.

第1 类改进方案属于经典力学范畴,如对算法的参数进行调节,文献[1-3] 对粒子群优化算法(Particle swarm optimization,PSO) 的惯性系数进行研究,分析惯性系数对全局搜索和局部搜索的影响,以此提高 PSO 的性能.或是融合其他算法的某种机制以提高自身性能,如文献[4]将粒子群优化算法的局部快速收敛性与模拟退火算法的全局收敛性结合,提出一种协同进化的算法,新算法在全局收敛能力和收敛速度上都有所提高;或是以算法迭代过程中的种群为入手点,改变种群的更新方式或将单一种群进行分组或划分,使其多种群化,如基于生物学上物种划分的小生境遗传算法[5].2010 年Luo 等[6]基于人工智能体的原始结构,将人工智能体划分为几个独立的子群体,不同的子群体中进行通信,这一方法被证实有效.同时也有学者基于拓扑结构或层次结构对 PSO 中粒子进行种群划分[7-9],以此平衡算法的勘探能力和开采能力[10].

第2 类改进方案是在原算法基础上引入量子机制,从量子物理与优化问题在概率意义上的相似性出发,通过引入量子机制提高算法的性能,在优化算法中应用较为广泛的量子机制可分为如下几种:

1) 量子波动的引入.1994 年 Finnila 等[11]在模拟退火算法(Simulated annealing,SA)中引入量子波动,利用系统中逐级递减的量子波动寻找经典系统的基态,提出量子退火算法(Quantum annealing,QA).量子波动的引入将函数优化问题转化为求解系统基态的问题,目标函数f(x) 则对应量子系统中的势阱,伴随量子系统的演化,粒子有一定的概率发生隧道效应,在量子系统中粒子能通过势垒贯穿出现在经典禁区(可映射为最优解区间),这称为隧道效应.隧道效应增强了算法跳出局部最优的能力,但量子退火算法并未对隧道效应的发生概率进行主观控制.

2) 势阱等效.通过薛定谔方程建立采样函数与目标函数f(x) 之间的关系,文献[12-13]以 Delta势阱下的波函数模方引导粒子位置更新,提出量子粒子群优化算法(Quantum PSO,QPSO),将粒子寻优空间等效为势阱,粒子的搜索行为受势阱约束,全局最优位置映射为势能最低点位置.王鹏等[14]通过模拟量子谐振子从高能态向基态收敛过程,将目标函数以谐振子势进行等效,提出用于解决高维函数优化问题的多尺度量子谐振子优化算法(Multiscale quantum harmonic oscillator algorithm,MQHOA),并于 2019 年通过模拟量子系统下自由粒子的运动过程解决函数优化问题,提出了基于自由粒子模型的优化算法框架[15].通过势阱等效,采用某一势阱下的波函数对粒子搜索过程进行概率描述,较经典力学的种群确定性演化过程而言该体系下算法的随机性更强.

3) 叠加态原理的应用.首先是 Narayanan 等[16]在 1996 年通过多种群的方式将量子多宇宙的概念引入遗传算法;随后,Han 等[17]将量子编码和量子旋转门应用在种群的编码和更新过程中,提出量子遗传算法,叠加态原理在遗传算法中的应用提升了算法的寻优能力,采用量子编码的染色体可表达多个态的叠加,通过量子手段增加了种群的多样性,类似的改进工作有量子蚁群算法[18];同时文献[19]采用量子编码的方式对差分进化算法进行改进,并在电力系统的优化问题中取得了不错的效果,将量子叠加态的概率特性以量子编码的形式应用在已有的算法改进工作中已被证明是有效的.

无论是量子退火算法中以递减的量子波动寻找量子系统基态,还是通过势阱等效的方法使粒子在量子系统下寻找势能最低点,在算法的迭代过程中,粒子都有一定的概率发生量子隧道效应,发生隧道效应的概率在函数优化问题中可衡量采样点从局部最优跳到全局最优的能力.隧道效应的发生是量子系统特有的现象,目前应用量子机制改进算法的方案中虽然涉及隧道效应的发生,但未对隧道效应发生概率有效控制以及隧道效应概率大小对算法性能影响进行研究.势垒与粒子簇对量子隧道效应的发生概率具有较大的影响,因此在算法迭代过程中对这两个变量信息的获取与利用变得尤为重要.本文以经典力学中多种群搜索方式为基础,将子种群的个体数映射为某一粒子簇中粒子的数量,将子种群间分布的相对位置用于评估当前粒子簇(采样点)到函数最优解间的势垒宽度,将经典力学下的多种群搜索策略赋予量子意义,通过动态调节迭代过程中子种群规模而增加种群整体量子隧道效应的发生概率,达到增强算法性能的目的,从控制隧道效应发生概率的角度为算法改进提供一种新思路.

本文结构安排如下.第1 节阐述优化算法的隧道效应,并建立多种群搜索策略与隧道效应发生概率的关系;第2 节和第3 节以 MQHOA 为例,说明其优化过程背后的物理背景,同时从发生隧道效应的角度分析该算法的不足之处,证明通过多种群手段增加隧道效应发生概率的必要性,然后基于蒙特卡洛评估的思想提出子种群规模自适应的策略;第4 节给出改进后具体的算法MQHOA-d (Multiscale quantum harmonic oscillator algorithm based on dynamic subpopulation),并以双阱函数为例设计实验跟踪算法在寻优过程中发生量子隧道效应点的数量变化情况,同时说明隧道效应对算法寻优的影响;第5 节设计实验并对改进后的算法(MQHOA-d) 从寻优误差、计算资源消耗等方面进行分析.

1 多种群策略与隧道效应

在经典力学体系下,将优化算法的采样方式由单一种群采样改进为多种群采样,从某种程度上增加了算法迭代过程中解的多样性,由此提高原算法寻优能力.但该方法仍存在弊端,对于将原算法单一种群采样变为多种群采样这一类改进方案,普遍存在一个共性,即当算法对子种群的划分完成后,其子种群的数目是固定的,随着迭代进行,子种群的个体数并不会动态变化.

对于类似多种群化的改进策略,除种群间的协作或学习等方式对算法性能的提高外,采样过程中子种群规模的变化对算法性能也具有影响,大量的算法改进专注于种群间的行为方式,而忽略子种群规模自身对性能的影响.如果从量子力学的角度看待多种群策略,可将以多种群的形式代表的整体采样点中子种群的数量映射为粒子簇的数量,某一子种群规模(个体数) 映射为粒子簇中粒子的数量.当前子种群到最优解区间的距离用于势垒的评估,目标函数f(x) 用于势能评估,详细对应关系见表1.由此可将函数优化过程看作种群寻找势能最低点位置的过程,在量子体系下粒子向势能最低点的搜索过程中,粒子存在穿越势垒到达势能最低点的可能性,即采样点从局部最优跳到全局最优区间,这种现象称为隧道效应,也是量子系统特有的现象.而隧道效应的发生概率受子种群个体数与当前子种群到最优解区间的距离(势垒) 所影响,通过势垒信息自适应地调节子种群规模,可改变总体隧道效应的发生概率,达到增强算法寻优能力的目的,同时弥补了经典力学下多种群改进方案忽略子种群规模自身影响算法性能的弊端.

表1 量子理论与优化算法的对应关系Table 1 The relationship between quantum theory and optimization algorithm

通过定义优化问题的波函数,可将最优区域[x0-ε,x0+ε] (x0为最优解,ε为求解误差) 发生隧道效应的概率通过波函数进行计算,此处记为Pi.则通过量子隧道效应在下一次迭代中出现在最优区域的个体数为×Pi,其中k为种群数,dyni为第i个子种群的规模(个体数).在寻优过程中,发生隧道效应的粒子将出现在最优解区域,由于多种群采样方式具有一定的信息交流,发生量子隧道效应的个体可引导所属的子种群向最优解方向迁移,以此达到由各子种群组成的种群整体向最优区域收敛.

这里不妨假设一种情况,某一子种群i的规模(个体数) 为dyni,且当前区域内发生隧道效应的概率为Pi,则该子种群在下一次迭代过程中通过隧道效应达到全局最优区间的个体数期望为Pi×dyni.无限制的增大子种群规模dyni可以增强该种群向最优解区间的迁移能力,但该操作会改变种群采样率.文献[20]指出增加采样点(种群规模)的数量会提高采样率,而采样率越小算法效率越高,出于上述考虑,基于蒙特卡洛估计的思想[21],提出可以根据种群适应度动态调节子种群个体数的策略.该策略在不改变算法采样率的基础上改变了迭代过程中种群发生隧道效应的概率,由此提高了算法性能,同时这一策略可以广义地应用在各种优化算法的改进工作中.

第2 节将以 MQHOA 为例,从控制隧道效应的角度对该算法进行改进.选择 MQHOA 作为改进案例的原因在于该算法是基于量子理论提出的优化算法,若选择经典力学体系下的优化算法作为目标算法,量子机制的引入容易从多个维度增强算法随机性(如种群在量子体系下的叠加性),虽然能增强算法性能达到改进的目的,但无法准确地验证通过子种群规模自适应改变隧道效应的发生对提高算法寻优能力的有效性.而 MQHOA 以薛定谔方程为基础,建立了优化问题与量子系统之间的关系,将优化问题转化为求解粒子在约束下的基态波函数.MQHOA在迭代过程中会自发地产生隧道效应,以 MQHOA为目标算法,通过动态调节子种群个体数的方式改变其迭代过程中发生隧道效应的粒子数,可以更好地验证通过隧道效应改进优化算法的有效性和可行性.

2 标准 MQHOA

2.1 MQHOA 物理模型

量子力学中波函数可以描述粒子在某位置出现的概率,通过文献[22]可知,薛定谔方程可以描述量子系统,不随时间变化的薛定谔方程表示为

通过薛定谔方程的变换,可将求解目标函数最小值问题转化为求解谐振子势阱下的基态波函数的问题.经过近似操作后,波函数在基态的概率密度函数为

由式(3) 可知,目标函数的最优解的概率分布为高斯分布.MQHOA 算法以高斯采样为基本操作,k个高斯采样点作为粒子,利用k个高斯采样叠加作为波函数描述当前最优解的概率分布.在 MQHOA 的收敛过程中,谐振子波函数从高能态逐步变向基态,从高能态多个正态分布的概率叠加,逐步收敛到基态单一正态分布的稳定状态,此时算法获得最优解.

2.2 MQHOA 基本流程

MQHOA 算法流程如算法1 所示.

算法 1.MQHOA 算法

算法1 中,参数k为采样中心数,S为第1 次采样区间,σ0和σmin分别为初始尺度和最小尺度,其中σmin对算法解的精度具有重要影响.f(x) 为目标函数.伪代码中从内到外的3 层循环分别对应MQHOA 的3 种物理过程,分别是能级稳定、能级下降、尺度下降.第1 层循环即对应能级稳定过程,在此过程中以k个采样点为中心分别按正态分布N(xi,)在每一个维度进行采样,以此生成新的采样点.计算当前采样点的目标函数值f() 是否小于当前采样中心的目标函数值f(xi),满足则进行替换,不满足则暂不替换当前采样中心.此外,如采样过程中受优化空间限制,则使用映射操作将越界点映射到优化空间内,本文不讨论上述映射机制.

当k个位置的采样中心全部完成采样后,计算x0～xk的方差σk,然后与迭代前的方差做差,计算出迭代前后方程的变化量 Δσk.如果Δσk ＜σs则判定算法能级稳定过程结束,进行能级下降操作,以k个采样点的均值替换当前目标函数值最差的采样点,即xworstxmean.当算法基态判据满足σk ＜σs时,可认为算法达到了尺度σs下的能量基态,将进行尺度下降过程,这个过程通过尺度减半的操作,即σs/2 降低当前尺度,这一过程决定了搜索精度,随着尺度下降这一步骤的进行,算法从全局搜索逐步过渡到局部搜索.

3 标准 MQHOA 存在的问题及解决方案

3.1 标准 MQHOA 存在的问题

原算法(MQHOA)忽视了不同采样中心在整个迭代过程中的适应度值是变化的,不同采样中心发生隧道效应概率不同,因此在每个中心位置采取相同的采样行为是不合理的.在解空间中,当部分采样点陷入局部最优区域时应增强该区域种群通过隧道效应跳出该区域的能力,由于原算法存在上述缺陷,由此对部分函数寻优时容易发生算法“早熟”现象.

MQHOA 搜索过程为k个正态概率分布不断向最优解位置迭代的过程,即算法的波函数为k个正态分布的叠加,可表示为

以一维函数为例,xi(i1,2,3,···,k) 为当前各正态分布期望,ε为误差上限,存在k个正态分布Ni(xi,).对k个正态分布约束下的任一采样点xi,在当前采样中能找到最优解的概率为

此时以该正态分布进行采样的种群发生量子隧道效应采样点的数量计算式为

其中,x0为最优解位置,dyni为以当前正态分布采样的子种群个体数.MQHOA 为单一种群采样,可看作k个规模为1 (采样点) 的子种群集合所构成的种群整体(单一种群),因此对 MQHOA 而言,dyni1,i1,2,3,···,k.由于采样中心xi的不同,进行 Ni(xi,) 采样时,采样点落在区间[x0-ε,x0+ε]的概率具有差异性(在高维函数解空间中任意两采样中心x1,x2相等的概率极低,因此不考虑这种情况).图1 所示的2 条高斯曲线分别代表采样中心x12,x28 的概率密度分布情况,X轴上的黑点代表最优解x0的位置,阴影部分为采样点落在区间 [x0-ε,x0+ε] 的概率P(xi).以x2 为中心的高斯曲线下阴影部分面积大于以x8 为中心的高斯曲线下的阴影部分面积,即在这次采样过程中以x2 为中心的采样区域能获得最优解位置的概率大于以x8 为采样中心的区域.因此,采样成功的概率会随采样中心的变化而产生差异.

图1 不同采样中心概率密度曲线Fig.1 Probability density curves of different sampling centers

本文将采样中心在当前迭代中能找到最优解的概率称作适应度值,适应度的大小以采样点与最优解间的势垒宽度评估.在采样点所构成的集合中,采样点xi通过隧道效应跳到最优解区域的概率最大,本文仅需比较k个采样中心适应度的排名,不需要进行精确计算,且当前适应度值最大的采样中心xbest应满

式中,xr为除具有最大适应度采样点xbest外的任一采样点.原算法忽略了对适应度不同的采样中心进行区别采样,导致种群在适应度较小的区域采样时,不能准确获取最优解位置信息,即能级稳定过程缓慢且易陷入局部最优解,同时影响收敛速度.实验表明,当目标函数的维度较高时,在规定的迭代次数中,MQHOA 的求解成功率会大幅度降低.当算法随着迭代次数的增加,从全局搜索过渡到局部搜索后,这种采样策略在搜索能力上的不足开始体现,寻优误差较大.

3.2 基于蒙特卡洛估计的种群数分配

函数的优化问题从某方面而言是具有黑盒特性的,在优化问题中目标函数的表达式f(x) 是不确定的,在迭代过程中,每一个子种群发生mi次测量,mi为当前子种群的个体数(规模),实际上可以将采样映射为对函数f(x) 积分的评估.通过第2 节可知,增加以当前中心xi的采样次数mi可以增大找到最优解概率,本节内容具体分析对不同mi之间是否应该具有差异.基于蒙特卡洛思想[23],可知采样越多,越近似于最优解.对于函数优化,在最优解xo所在误差范围 [x0-ε,x0+ε] 内,应满足

由于MQHOA 采样过程中尺度σs是不断递减的,对于k个子种群搜索过程而言,实际上是评估k个子种群在当前尺度下积分最小的问题,表达为

由于波函数概率密度函数为正态分布,函数采样区间有限,这里做近似处理

式(9)是以xi(i1,···,k) 为中心的子种群对当前区域的评估蒙特卡洛算子,mi(i1,···,k)为该子种群的规模(个体数),P(xi)(i1,···,k)为当前第i个子种群的概率密度函数,其中,N为采样次数,它与mi相对应,f(x)为目标函数,P(x)为x的概率密度函数,记f(xi)/P(xi)g(xi),据大数定律可知

由此可以推出区间内积分在N次重复采样下的表现形式为

对于不同中心Xi而言,定义域内的采样概率密度函数具有差异,而固定区间内对函数的积分为定值,此时应分析概率密度函数对积分评估的影响.我们知道所有采样中心点服从正态分布的方差是相同的,不同的是期望,而期望值为采样中心的坐标信息.此时计算方差为

记g(xi)方差为σ2,此时g(xi) 相互独立,因此蒙特卡洛方差与σ2的关系为

4 基于多种群的MQHOA-d

4.1 MQHOA-d 流程

k个种群以当前种群中心的波函数为基准进行采样,每个子种群的个体数mi(i1,···,k) 的更新是根据当前中心到适应度最高种群中心的欧氏距离决定的,更新式为

其中,n为维度大小,l为当前维度,Xi为当前种群中心,Xbest为适应度最高的采样中心,依据欧氏距离的种群规模更新式为

式 (14)通过各子种群中心Xi到当前最优子种群中心Xbest的欧氏距离评估各子种群所处解空间的适应度,距当前最优子种群中心距离越远,则判断该子种群适应度越差,由此通过式 (15) 更新该子种群规模(个体数).该操作在量子力学上的意义在于通过评估粒子簇与最优解区间的势垒宽度(式(14)对子种群中心欧氏距离的计算) 调节粒子数量(式(15)对子种群规模的调节),从而改变该区域发生量子隧道效应的概率.

能级稳定判据对算法在当前尺度下种群进行充分搜索有重要影响,当k个种群迁移速度过快,达到能级稳定判据标准时,会进行能级下降,种群分布聚集,由此容易导致算法早熟,全局搜索不充分.因此,需要在每轮迭代过程中对当前尺度下的种群空间进行限定,即实现种群迁移步长控制.MQHOAd 在迭代过程中,首先初始化k个种群中心,以当前尺度σs2生成n维协方差矩阵,具体计算为

其中,cov(g)表示第g代的协方差矩阵,在当前种群中心Xi处按照正态分布 N(Xi,cov(g))生成个采样点,完成当前种群个体初始化.当种群方差满足能级稳定判据时 (Δσk ≤σs),MQHOA-d 进行能级下降操作,需要将适应度最差的种群所有个体淘汰.取剩下k个种群中心位置的均值坐标,初始化一个新的种群中心m(g),具体为

以Xi为中心通过式 (14)和式(15) 生成种群个体,用新种群替换当前最差种群,完成替换过程,实现能级下降.当所有的个体满足能级下降判据,此时在当前尺度的搜索过程完成,下一步在更小的尺度进行搜索,进行尺度下降操作.MQHOAd 的尺度下降步骤与MQHOA 相同,仍为尺度减半操作.当满足基态判据时,输出当前最优解,寻优结束.

图2 表示k为30、λ为1 000 时,MQHOA-d所生成的种群在 Ellipsoidal 函数上种群分布示意图.函数最优解位置为(1,2,3,···,n),粗点代表当前种群中心Xi,细点为采样点(当前子种群个体).本文中子种群是指以当前位置Xi为均值,根据正态分布 N(Xi,cov(g))生成的个采样点,即属于同一个子种群中的个体 (采样点) 变更满足N(Xi,cov(g)).由图2 可见,依据式 (14)和式(15)生成的种群特点如下: 适应度差的区域分布子种群规模大,即个体数多;适应度优的区域子种群分布规模小.这一操作的目的是降低第2.1 节中对适应度差的子种群空间的评估误差,以此保证整体评估的稳定性.

图2 MQHOA-d 所生成种群在Ellipsoidal 函数二维空间分布示意图Fig.2 Schematic diagram of spatial distribution of subpopulations generated by MQHOA-d in Ellipsoidal function of 2D

MQHOA-d 伪代码见算法2,k为种群的数量,S为第1 次采样区间,σs和σmin分别为当前尺度和最小尺度,其中,σs初始化为σs=UB-LB,初始化的k个采样中心.(g1,2,···,k) 代表每一个种群的规模 (个体数).

算法 2.MQHOA-d 算法

4.2 量子隧道效应跟踪实验

多种群化可以增加种群发生量子隧穿的概率,为了更好地分析该策略与粒子发生量子隧道效应的关系,用 MQHOA/MQHOA-d 做对照实验.实验组1 中 MQHOA 采样方式为单一种群采样,种群数为1,个体数通常为30,这里为控制参数一致,将个体数设置为200;实验组2 中 MQHOA-d 采样方式为多种群采样,设置子种群数为30,个体总数为200.在迭代过程中可以根据当前波函数的概率密度计算发生量子隧道效应点的数量,即

图3(a) 中实线为实验所用函数的一维示意图,其中加粗部分为最优区间,图3(b)为子种群中部分采样个体发生量子隧道效应,贯穿势垒跳出局部最优示意图.由式(18)可知在下一次迭代过程中,该子种群发生隧道效应的个体数量与当前子种群规模正相关,式(15) 通过适应度的差异性调节各子种群的规模,由此增加适应度差的子种群个体数,以此提高当前能级下该子种群量子隧穿到最优解区域的概率,增强该子种群向最优解区域的迁移能力.隧道效应的发生是粒子数量与势垒(本文映射为适应度) 共同影响的结果,固定子种群规模的搜索策略忽略了子种群个体数量对隧道效应的影响,因此,基于子种群自适应策略改进的 MQHOA-d 在相同的环境下较固定种群规模搜索策略的算法发生隧道效应的概率更大,由此可以提升算法搜索能力.

图3 一维双阱函数图与隧道效应示意图Fig.3 One dimensional double well function diagram and tunnel effect diagram

实验的目的之一是计算种群中通过量子隧道效应出现在最优区间的点的数量,虚线为随迭代次数增加采样概率密度函数尺度递减示意图,当迭代次数增加到一定数量时,概率密度函数中σs收敛于0,此时发生量子隧道效应的概率也趋近于0,因此会发生图4 中随迭代次数增加,量子隧道效应的点的数量先增加后递减甚至趋于0.图4 中三角形折线和圆形折线分别代表 MQHOA-d和MQHOA种群中发生量子隧道效应的点数量与迭代次数的关系图,MQHOA-d 与 MQHOA 的迭代过程中,在相同的迭代次数下,MQHOA-d 发生量子隧道效应点的数量要大于 MQHOA,实验设置的输出条件为迭代次数达到10 000.通过图4 可知,MQHOA-d 在迭代次数达到约2 500 时尺度σs收敛于0,即种群已完成搜索,而 MQHOA 需要迭代约9 500 次才完成搜索.从发生量子隧道效应点数量与迭代次数的关系图可知,快速地达到隧道效应点数量最大化有利于缩减算法搜索时间,提高种群收敛速度,这也是对第3.1 节中通过欧氏距离评估适应度而动态调节子种群规模使整体发生量子隧道效应最大化的有效性验证.

图4 MQHOA-d 与MQHOA 发生隧道效应点数量与迭代次数关系图Fig.4 The relationship between the number of tunneling points and the number of iterations between MQHOA-d and MQHOA

发生量子隧道效应的点数量越多可映射为当前种群分布寻找到最优解区域的概率越大.MQHOA-d在迭代次数较小时即可达到整体种群发生隧道效应的最大状态,此时尺度的缩减程度低,可以说明MQHOA-d 在全局搜索时即可快速定位全局最优区域,面对复杂函数或多峰函数时更具有优势;而MQHOA 在整体种群发生隧道效应达到最大状态时需要的迭代次数较大,此时尺度σs缩减程度高,在大尺度下对全局的评估不如 MQHOA-d,由此可以推测当测试函数为结构复杂的函数时,MQHOA容易陷入局部最优区域,且 MQHOA 的寻优能力不如 MQHOA-d.从蒙特卡洛估计积分的角度上看,合理降低适应度差的子种群采样方差可以提高对当前种群整体分布的评估准确度,有利于解空间的全局评估.

5 实验及分析

为了全面分析 MQHOA-d 的性能,选取烟花算法 (Fireworks algorithm,FWA)[24-25]的衍生EFWA (Enhanced FWA)[26]、AFWA (Adaptive FWA)[27]以及 MQHOA、SPSO2011[28]、QPSO[29-30]作为对照组.

所有的算法在 CEC2013 标准函数测试集[31]上进行,每组实验重复51 次,单次运算的最大迭代次数为10 000 ×d,其中d为维度,实验维度设置为30.测试集一共包括28 个测试函数,根据函数结构不同可划分为单峰函数 (f1～f5)、多峰函数(f6～f20)和复合函数 (f21～f28).所有函数的定义域为[-100,100],后续对实验结果的误差进行分析.实验环境为Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz,2.21 GHz,16.0 GB,程序采用 MATLAB2018a 编写.其中,MQHOA、MQHOA-d 的种群数设置为30,MQHOA-d 的个体总数λ设置为200;QPSO和SPSO2011 粒子总数设置为30.具体实验数据见表1(EFWA、AFWA 实验数据来自于文献[32]).表1中 MeanErr 表示误差均值,Std 代表解的标准差,最优的实验结果在表2 中用粗体标出.

5.1 实验数据分析

为更直观地分析数据,将6 组实验数据的误差均值进行等级划分(1 级～ 6 级),均值误差最小的为1 级,最大的为6 级,同时计算对照组的平均等级并绘制成图5.

图5 CEC2013 测试集上各算法误差等级图(维度为30,重复51 次)Fig.5 Rank sum test results of all algorithms on CEC2013 test set (repeated 51 times in the 30-dimensional case)

与原 MQHOA 算法相比较,通过表2和图6可以看出,MQHOA-d 在单峰函数上 (除f3) 表现较差,这是由于 MQHOA-d 每一个高斯采样空间中产生的个采样点并不能实现向最优解空间的快速迁移,一次采样所需要的计算次数为而这一操作对应算法的能级稳定过程: MQHOA 仅在每一个高斯采样空间中生成一个采样点,即一次采样仅需要计算k次.而结构简单的单峰函数并无过多局部最优解,多种群搜索的策略虽可以提高跳出局部最优解的概率,但针对这类型的函数,需要的是整个种群快速向最优解区域收敛,即搜索空间从大尺度向小尺度的过渡,这一过程对应 MQHOA/MQHOA-d 的能级下降过程,算法需要尽快从全局过渡到局部.而 MQHOA-d 的多种群生成决定了它将大量的计算资源消耗在能级稳定的过程中,由于实验计算次数存在上限,因此在有限的计算中,MQHOA-d 进行能级下降的次数比 MQHOA 少,即尺度缩减程度不如 MQHOA.

从多峰函数和复合函数上看,MQHOA-d 的求解误差要远小于原 MQHOA,这也验证了第4.2 节中根据量子隧道效应实验数据对算法性能的推测是可靠的,MQHOA-d 在大尺度的状态下尽可能多地进行搜索,全局搜索更加全面.此策略在适应度低的区域生成更多的采样个体,对整个种群的计算能力进行合理分配,提高了种群多样性,增加跳出局部最优解概率.因此,MQHOA-d 对结构复杂的多峰函数和复合函数更具有针对性.相比于 SPSO2011、QPSO、EFWA、AFWA 等算法,通过图6 可知,MQHOA-d 在单峰函数下处于劣势,但多峰函数和复合函数的平均等级最小,且总体平均误差等级最小.综合28 个测试函数分析,MQHOA-d 的平均误差等级值最小,仅为1.93.从表2 中可知,MQHOAd 的标准差在多峰函数和复合函数中较小,具有较强的鲁棒性.通过每组实验重复51 次的数据绘制箱体图(具体见图6).这些箱体图清晰地显示出MQHOA-d 在多峰函数和复合函数中的准确性和稳定性占据主导地位.

MQHOA-d 性能的增强得益于搜索随机性的增强,在量子模型下该算法的随机性可由量子隧道效应衡量.种群个体数自适应过程扩大适应度差的子种群规模,以此增强搜索过程的随机性,同时降低种群对目标函数的评估误差.在量子模型下解释为当势垒宽度较宽时,增加粒子数量可以加大粒子簇发生隧道效应的概率,非量子模型下的优化算法通常以增加扰动项或设置变异操作的方式增强算法的随机性.如文献[30]中通过评估种群所处的阶段实现变异策略的反馈调节,由此提出基于状态估计反馈的策略自适应差分进化算法 SEFDE,且 SEFDE 在与 CEC2013 中性能最好的 DE 改进算法SHADE 的对比实验中表现出较强竞争力.将表2 中MQHOA-d 的 CEC2013 测试数据与文献[33]中SEFDE 在同等条件下的 CEC2013 测试数据进行对比,结果显示,MQHOA-d 在16 个函数优化上优于SEFDE,该16 个函数除f2外,都为多峰函数或复合函数,与此同时,MQHOA-d 在函数f8的优化上与SEFDE 具有相似的优化效果.总体而言,MQHOAd 在多峰函数和复合函数的优化上优于 SEFDE.这也验证了基于蒙特卡洛评估提出的动态种群策略的合理性,增大适应度差的种群规模可以提高算法搜索随机性,有利于降低整体评估误差.

统计51 次重复实验的平均迭代时间,绘制成图7.在图7 中,方框中数字代表算法在对应测试函数中所消耗的平均时间,单位为s.将时间归一化后进行涂色,颜色越深的区域代表消耗时间越多.

本次实验中,算法MQHOA 在测试函数f1上的平均迭代次数为 3.49×104,小于上限D×10 000,其他实验组的平均迭代次数均达到上限D×10 000.针对测试函数f1而言,MQHOA 所消耗时间较其他3 种算法少一个数量级,且图6 中箱体图显示MQHOA 的51 次求解数据分布不稳定,可见 MQHOA 在f1测试函数上存在算法“早熟” 现象,过早地收敛在某一局部最优无法跳出,也验证了增大种群初期发生量子隧道效应的概率更不易陷入局部最优.在测试函数f9上,MQHOA-d 所消耗时间较其他3 种算法具有极大的优势,平均误差仅次于 QPSO,且误差范围与 QPSO 在同一个数量级,当实际工业生产面对结构类似于f9函数时,可在节约计算成本的情况下使用 MQHOA-d.总体而言,MQHOA 算法消耗时间最少,但寻优精度远不如其他3种算法;MQHOA-d 消耗时间仅次于 MQHOA,但寻优精度和鲁棒性强于 MQHOA、QPSO和SPSO2011.

6 结束语

本文基于蒙特卡洛估计中不同子种群对总体评估误差的差异性,提出动态调节子种群规模的寻优策略,将此策略应用于 MQHOA 的改进,改进后的算法(MQHOA-d)在未增加随机化特征,或者添加扰动等多余步骤的前提下解决了原算法的“早熟”现象[34].核心在于动态调节种群规模方式增大了种群迭代初期发生量子隧道效应的概率,使其对解空间全局评估更加精准.在迭代次数相同的情况下MQHOA-d 对多峰函数和复合函数寻优误差更小,根据适应度评估的多种群采样方式动态地调节了种群规模,在解空间中更加合理地进行勘探和开发,同时种群利用率得到了提高.改进后的算法虽然结构上较之前相比变得复杂,但时间复杂度上仍优于几大主流算法(SPSO2011、QPSO).目前,MQHOAd 仍有局限性,如: 1) MQHOA 系列算法的随机性和隐含并行性尚不明确[35-36];2) MQHOA-d 种群个体行为所抽象的随机游走模型[37-38]尚未研究.