多物态拟态物理学优化算法

安泽晔，谢丽萍

(太原科技大学计算机科学与技术学院，太原 030014)

求解最优化问题的启发式算法[1-3]中，包含受生物群体行为启发的微粒群算法[4]，模拟自然进化规律的遗传算法[5]，还有受物理学原理启发的类电磁机制算法[6]，中心力算法[7]，模拟退火算法[8]和拟态物理学算法[9]等。其中，拟态物理学算法(APO)是从自然物理规律中得到启发，模拟两个物体之间虚拟的作用力，然后仿照牛顿第二力学定律进行算法寻优的过程，该算法作用于种群个体之间，建立起算法中的寻优个体与物理学中的个体之间的对应关系。

近年来，对拟态物理学优化算法的研究有很多。有对其算法框架的研究[10-11]，对算法性能的研究[12]，对收敛性分析及对参数G的设置的研究[13]，对其在实际中的应用的研究[14-18]。但 APO算法在迭代过程中只遵循一种运动规则，使所有寻优个体体现相同的运动特征，不能有效的平衡全局搜索和局部搜索能力，而单一的搜索策略往往不能得到更好的解，也不能更好的贴合实际，解决问题。由此为了增加种群的多样性，引入多物态的概念，并对物态划分的方法进行了研究。

物理学中的宏观物质状态包含了三种基本的状态，即固态，液态和气态，而这三种基本状态在微观世界中由于分子间的布朗运动又可以相互转化，不同物态分子的热运动随分子之间的结构不同而运动方式不同。将个体映射到物理学中具有物态特性的个体，不同“物态”的个体采用不同的搜索策略，并且由于物理环境中温度变化而导致的状态转化行为，算法同样映射了不同物态个体的行为方式。个体的运动规则融合了不同的物理学运动规律，包括万有引力中的引斥力机制，模拟退火算法的单点运动规则以及宇宙大爆炸算法中的质心爆炸过程。这三种运动规律近似地模拟了固态，液态，气态的运动特性，在本文中，固态集合的运动规则采用APO算法的运动规则，液态个体的运动仿照模拟退火算法的单点运动规则，气态个体运动仿照宇宙大爆炸算法每代从质心位置向外爆炸的过程，使个体尽可能遍历更多的点。本文利用这三种物态的粒子运动为解决算法中单一的搜索策略易于早熟的状况。

1 标准APO算法

APO算法[9]主要包括三个部分，首先初始化种群。设种群规模为N，维数为D，个体位置xi和运动速度vi在可行域内均匀取值。其次计算个体质量，质量公式为

(1)

在t时刻个体i受到的作用力合力表示为：

(2)

最后更新个体运动的速度和位置。

∀i≠best

(3)

xi,k(t+1)=xi,k(t)+vi,k(t)

(4)

其算法流程如下：

Step1：初始化种群的个体位置和速度，并计算适应值，保留算法的最优适应值及最优个体；

Step2：计算个体质量及个体所受合力；

Step3：更新个体速度和位置；

Step4：计算个体适应值，更新最优个体及其适应值；

Step5：判断是否满足程序结束条件，若满足，则输出最优结果并结束程序，否则返回Step2，继续迭代。

2 多物态拟态物理学算法的融合

多物态物理学算法的融合主要包括三个部分：初始化种群，物态划分，计算不同物态的个体运动。初始化种群后，按一定的划分标准划分固态，液态和气态个体，首先分离出气态个体，然后根据其他个体与选出的中心点个体间的距离设定固态和液态的邻域集合，集合内含两个或两个以上的个体归为固态集合，只有一个个体的归为液态集合。物态划分后，个体按不同的运动规则更新速度和位置，计算每个个体的适应值，然后根据个体的适应值及环境，判断某些个体的物态是否需要发生变化，若需要，则重新划分物态，并按相应的运动规则进行运动。

2.1 运动规则

本文根据不同物态的个体运动规律不同，选用了三种不同的运动规则。

规则一：标准APO算法的运动规则

规则二：仿照模拟退火算法的单点运动规则。

速度和位置计算公式如下

∀i≠best

xi,k(t+1)=xi,k(t)+(2*rand-1)L

(5)

vi,k(t+1)=vi,k(t)+xi,k(t+1)

(6)

规则三：仿照宇宙大爆炸算法，每代从质心位置向外爆炸，个体运动方向和运动步长随机，其中质量m计算如公式(1)。

(7)

(8)

vi,k(t+1)=xi,k(t+1)-xi,k(t)

(9)

2.2 物态划分方法

APO算法只遵循一种作用力规则，个体运动方式较为单一，算法的种群多样性相对较差。在物理学中，固态和液态物质的运动比较局限，而气态物质扩散性好，映射到算法中，适应值较好且被划分为固态和液态的个体适于做局部搜索，而适应值较差的个体则对应于气体的扩散性质，适于做全局搜索。这种划分方法增加了种群的多样性，简单易操作。但若是直接通过将适应值大小进行排序，挑选前K个较好个体作为固态或液态个体，那么K的选择至关重要，有可能出现第k个和k+1个个体适应值相同或相似的情况，从而导致物态的划分不合理。当k值过大，可能导致全局搜索能力较弱，k值过小，又可能导致算法收敛过慢的情况。由此，这里引入了一种物态(包括固态、液态和气态等)的划分标准Cls，如公式(10).该公式将个体的适应值映射到(0,1)的区间内，避免具有相同适应值的不同个体被划分为不同的物态，增加了划分的合理性。

(10)

当xi适应值较大，clsi较小，将该个体设为气态个体，用于全局搜索；当xi适应值较小，clsi较大，将该个体设为固态或液态个体，用于局部搜索。

该算法中取cls_standard为划分的标准参数，当个体的clsi>cls_standard时，划分为固态或液态个体，相反划分为气态个体。该标准基于每一代适应值的最大值和最小值，用较为相对的标准划分物态，增加了物态划分的合理性。

算法描述

该方案在APO的基础上引入物态划分的概念，以公式(10)的cls作为物态划分的标准，在初始化种群后，便进行物态的划分，以期提高种群的多样性。划分之后，各物态集合的个体按照不同的运动规则运动，之后保留最优值并判断是否达到结束条件，达到则结束程序，未达到则重新划分物态，进行运动。在此算法中，划分为固态集合的个体按照APO算法中的引斥力规则运动，划分为液态集合的个体按照模拟退火运动，划分为气态集合的个体按照宇宙大爆炸的方法运动。该算法的算法步骤如下：

Step1：初始化。在可行域内均匀产生个体的位置与速度，并计算相应的适应值，保留最优个体及其适应值。

Step2：物态划分。按照公式(10)计算每个个体的cls值，clsir的个体作为液态个体。

Step3：各集合按相应的运动规则运动一次。每个固态集合都按照规则一运动，液态个体按照规则二运动，气态个体按照规则三运动。

Step4：计算适应值，保留最优个体及适应值。判断是否达到结束条件，是则结束运行并输出结果，否则返回步骤2.

3 实例仿真

为了测试上述方案的可行性，选用CEC2013 Benchmark 测试集中的Sphere Function, Rotated Rosenbrock’s Function, Rotated Schaffers F7 Function, Rotated Ackley’s Function, Rotated Griewank’s Function, Composition Function 1(n=5,Rotated)六个测试函数。这些测试函数的参数如表1所示，其中问题维数(n)、函数的搜索空间(Search range)、函数的已知最有值(Known optimum)、算法使用的种群规模(Npop)以及算法迭代的最大代数(MAX-ITER)已列出。

表1 测试函数参数列表
Tab.1 The parameter list of test function

FunNNpopSearch rangeKnown optimumMAX-ITERSphere Function3030[-100,100]n-140010000Rotated Rosenbrock’s Function3030[-100,100]n-90010000RotatedSchaffers F7 Function3030[-100,100]n-80010000Rotated Ackley’s Function3030[-100,100]n-70010000Rotated Griewank’s Function3030[-100,100]n-50010000Composition Function 1(n=5,Rotated)3030[-100,100]n70010000

3.1 MSAPO与APO算法性能分析

在仿真实验中，为了方便算法的比较和分析，将引入cls参数的多物态物理学算法命名为MSAPO.MSAPO中模拟退火运动规则的温度初始值T=100，衰减系数k=0.5，并按照经验值，邻域半径r=30，引斥力运动规则中的G值为0.1，物态划分标准cls_standard=0.5.原APO算法中G值取0.1.在相同的实验条件下，表1的测试函数在APO算法和MSAPO算法下的对比实验结果如表2所示。每个实验重复执行10次，统计最好适应值(Best)，平均最好适应值(Mean)，标准方差(STD)和平均误差(Mean_error).

表2 APO、MSAPO算法性能比较结果
Tab.2The performance comparison results of algorithm APO and MSAPO

FuncDimAlgorithmBestMeanSTDMean_errorSphere Function30APO-1.4000E+03-1.4000E+039.7019E-086.3109E-08MSAPO-1.4000E+03-1.4000E+038.3226E-114.7066E-11Rotated Rosenbrock’sFunction30APO-8.8355E+02-8.5805E+022.6578E+014.1950E+01MSAPO-8.8492E+02-8.5874E+022.5519E+014.1264E+01Rotated Schaffers F7 Function30APO-6.4900E+02-6.1893E+022.1656E+011.8107E+02MSAPO-7.8956E+02-7.5823E+023.5565E+014.1766E+01Rotated Ackley’s Function30APO-6.7909E+02-6.7905E+023.0563E-022.0949E+01MSAPO-6.7910E+02-6.7905E+023.0341E-022.0950E+01Rotated Griewank’s Function30APO-4.9985E+02-4.9900E+027.0865E-011.0045E+00MSAPO-4.9995E+02-4.9984E+025.9681E-021.5627E-01Composition Function 1(n=5,Rotated)30APO9.0000E+021.7743E+036.7707E+021.0743E+03MSAPO1.0000E+031.0287E+036.0524E+013.2871E+02

表2给出了相同代数下所测函数在APO与MSAPO算法下的对比实验结果。MSAPO算法在以上六个测试函数中找到的最优适应值均比APO算法好或者相等，其性能整体上也优于APO算法。其中，复杂函数Composition在MSAPO算法下的最优值比APO算法要差一些，但平均值和标准方差均优于APO算法。其他基本多峰函数在MSAPO算法下的性能都相对较好，除了Rotated Schaffers F7 函数，其在MSAPO算法下的标准方差和平均误差值较APO差。

Rotated Schaffers F7函数是复杂非线性多模态函数，拥有广泛的搜索空间和大量的局部极小点，MSAPO算法在该函数上平均值和最优值都得到明显的提升，这是因为APO算法具有较好的全局搜索能力，但算法的搜索方式较为单一，算法种群多样性差，局部搜索能力较弱，而MSAPO算法加入了物态划分和多种运动规则，使算法的搜索方式多样化，可以根据个体当前的状态选择相应的运动规则，增加了APO算法的种群多样性，相对于运动规则单一的APO算法而言，更容易找到测试函数的最优值。但SchaffersF7与compodition1函数表现不佳之处说明了该算法种群多样性虽然好，但是在局部极小值较多的复杂函数上的精细搜索不够，较不稳定。

3.2 物态划分标准参数对算法的影响

在引入物态划分的APO算法中，物态划分的标准参数为种群中各物态的个体个数分布比例情况提供了重要信息，其选择的好坏能够影响APO算法的性能。实验通过研究该物态划分标准参数的选择策略判断对算法性能的影响。

本算法取cls_standard为划分物态的标准参数，当个体的Clsi达到标准参数cls_standard时，划分为气态个体，相反则划分为固态或液态个体。该标准参数基于公式(10)给出，Clsi的取值范围为(0,1).当种群中划分为固态或液态的个体个数较多时，由于算法中固态或液态个体的运动规则偏向于精细搜索，可能使得整个算法的局部搜索能力随之增强，而全局搜索能力也就相应减弱，这样不利于遍历搜索空间中尽可能多的点。当种群中划分为气态的个体个数较多时，算法中气态个体的运动偏向于全局搜索，可能使得算法的全局搜索能力增强，但是容易丢失潜在较好解，且在较好解附近的精细搜索能力不足。本算法中Clsi基于当前最差和最好的适应值，cls_standard的选择对算法整体性能有一定的影响，实验选取几个不同的值，测试表1中函数的最优值，平均值，平均方差和平均误差值。

表3 MSAPO算法在不同参数下性能比较结果
Tab.3 The performance comparison results of MSAPO under different parameters

FuncCls_standardBestMeanSTDMean_errorSphere Function0.3-1.4000e+03-1.4000e+034.2298e-101.9463e-100.5-1.4000E+03-1.4000E+038.3226E-114.7066E-110.7-1.4000E+03-1.3907E+032.2787E+019.2907E+00Rotated Rosenbrock’s Function0.3-8.8469e+02-8.5381e+022.4726e+014.6187e+010.5-8.8492E+02-8.5874E+022.5519E+014.1264E+010.7-8.8510E+02-8.7796E+024.3128E+002.2035E+01Rotated Schaffers F7Function0.3-7.8512e+02-7.5826e+022.5490e+014.1737e+010.5-7.8956E+02-7.5823E+023.5565E+014.1766E+010.7-7.8872E+02-7.7323E+021.3612E+012.6773E+01Rotated Ackley’sFunction0.3-6.7916e+02-6.7906e+025.2205e-022.0936e+010.5-6.7910E+02-6.7905E+023.0341E-022.0950E+010.7-6.7911E+02-6.7903E+024.8177E-022.0970E+01Rotated Griewank’sFunction0.3-4.9994e+02-4.9984e+025.1473e-021.5896e-010.5-4.9995E+02-4.9984E+025.9681E-021.5627E-010.7-4.9995E+02-4.9988E+027.6477E-021.1633E-01Composition Function 1(n=5,Rotated)0.31.0000e+031.0431e+036.9338e+013.4306e+020.51.0000E+031.0287E+036.0524E+013.2871E+020.79.0000E+021.0474E+038.8151E+013.4742E+02

由表3的实验结果可得，对于简单函数和复杂函数而言，随着cls_standard的取值越大，算法性能变化不同。Cls_standard的取值关系到算法中不同运动规则的个体个数，当cls_standard的取值从0.3到0.5，简单函数sphere的平均值有所下降，而对于基本多峰函数和复杂函数而言，除了Ackley函数，其他函数的平均值均有变差趋势，但最优适应值都变好。这是因为当cls_standard的值越大，算法中被划分为固态和液态的个体随之减少，气态个体随之增加，算法越趋向于按照气态个体的运动规则进行全局搜索，对多峰函数和复杂函数的寻优影响较好，但可能会导致在潜在较好解附近的精细搜索不足，由此产生平均值结果变差的问题。当cls_standard的取值从0.5到0.7变化时，尤其对于Rosenbrock和Schaffers F7这两个典型的多峰函数，平均方差和平均误差值均变小，说明较大的cls_standard值对于多峰函数的影响较好。

实验结果表明该算法在选取cls_standard的值时，需要选取适当才能有效提高算法性能。这是由于MSAPO算法采用了多种运动规则，对不同状态的粒子采用不同的运动规则，而cls_standard值的选取能有效平衡算法的全局搜索能力和局部搜索能力。当cls_standard很小时，算法中做全局搜索的粒子较少，相反，当cls_standard很大时，算法中做局部搜索的粒子较少。对于简单函数而言，局部极小值较少，cls_standard值较小时容易找到最优值且收敛较快，对于多峰函数局部极小值较多，过小或过大的cls_standard值都不利于提高算法的性能。在实际应用中可以根据问题模型的复杂程度选择cls_standard的值。

4 结论

本文在原拟态物理学优化算法的基础上引入多物态的概念，对不同的个体划分不同的物态，并采用符合各物态运动规律的运动规则进行寻优。仿真实验表明该算法的有效性，后续工作将对MSAPO算法在求解最优化问题时的平均最优值精度进行进一步研究，通过引入环境因素，更合理地模拟物理学中分子的热运动现象，实现部分物态个体因环境变化而实现的状态转化，提高寻优结果的精度，改善算法的时间复杂度。