1 引言
在实际问题中，绝大多数优化问题由多个子目标所组成，并且这些子目标之间相互矛盾，相互冲突，一个子目标性能的改善可能会引起其余子目标性能的降低，这种在优化过程中要求多个子目标同时在指定区域内达到最优值的问题即为多目标优化问题[1]。针对多目标优化问题及其研究，越来越多的群体智能进化算法在自然行为的启发下得到发展，如多目标粒子群优化算法(MOPSO)[2]、非支配排序遗传算法(NSGAII)[3]、多目标分布估计算法(MOEDA)[4]等。
APO算法是一种新的基于群体智能的随机优化算法。通过牛顿第二定律的启发，建立个体之间的虚拟力规则，更新个体的速度和位置，使个体移动到最优目标，并围绕全局最优收敛[5-6]。APO算法已成功应用于无约束多目标优化问题。Wang Y等[7]针对APO算法自身的性能，从多目标优化问题出发，以虚拟力的角度分析多目标算法中个体-优劣，提出一种基于虚拟力排序的约束多目标拟态物理学优化算法。Sun B等[8]在MOAPO算法的基础上结合约束违反度的判断准则，提出一种基于序值和拥挤度的拟态物理学多目标算法(improved constrained rank multi-objective artificial physics optimization,ICRMOAPO)，但在解的筛选过程中将新的个体与archive解集中的个体进行比较，容易将优质解丢失，使其无法得到一组均匀且广泛的最优解。
近年来，对于多目标优化算法的选解机制，国内外学者做了大量的研究。Sedighizadeh D等[9]提出了一种连续空间粒子群优化算法，通过在速度更新方程中引入新的项，使得在搜索空间中更好地进行搜索。Cai D等[10]将目标空间分解为多个子空间，并在每个子空间上保留最优解以保持解集的分布性。Kohler M等[11]将多样性插入到群体中，扩大搜索空间的覆盖率，提高了算法的效率和收敛速度。Garcia I C等[12]提出了一种基于超体积指标的多目标粒子群算法，利用外部存档来存储进化过程中发现的全局非支配解，为求解多目标问题提出了一种新的思考方向。Li F团队[13]针对HV指标难以有效权衡计算时间复杂度和计算精度之间的矛盾进行研究，发现R2指标可以替代HV指标求解多目标问题，将R2指标和PBI分解策略结合求解多目标优化问题,提出了R2-MOPSO算法。
针对ICRMOAPO算法的不足，为了得到收敛于Pareto前沿且均匀分布的解，利用R2指标综合评价解集的特质，提出一种基于R2指标的拟态物理学约束多目标优化算法(R2-ICRMOAPO)。算法将非支配排序和R2指标结合，并根据R2指标的贡献值对外部存储集进行更新和维护，保留分布均匀的精英解。
2 约束多目标优化问题
一般地,约束多目标优化问题的数学模型可以被描述为：
min f(x)=[f1(x), f2(x),…, fk(x)]
s.t. g(x)=(g1(x),g2(x).…,gm(x))<0
h(x)=(h1(x),h2(x).…,hm(x))=0
(1)
式中，x=[x1,x2,…,xn]∈X为决策变量，决策空间X∈Rn，
满足边界约束条件：li≤x≤ui,1≤i≤n，li,ui分别为xi的上界和下界； f(x)为目标函数；g(x)<0为不等式约束条件，h(x)=0为等式约束条件，约束条件决定可行域的范围。
定义1 pareto支配
x Pareto支配y，记为x<y，其中x为支配解，y为非支配解，当且仅当∀i∈{1,2,…,m}, fi(x)≤fi(y)，∃j∈{1,2,…,m},s.t. fj(x)<fj(y)。
定义2 Pareto最优解集
对于一个约束多目标优化问题，它的可行决策空间中所有最优解的集合称为Pareto最优解集，令其为P。
P={x∈Xf|
∃x′∈Xf, x′>x}
(2)
定义3 Pareto前沿
Pareto最优解集中每个解对应的目标值向量组成的集合称为Pareto前沿。
约束多目标优化问题就是要得到一组解集，使之逼近Pareto前沿。在目标空间中则表示为，非支配解集覆盖的目标空间区域尽可能大，且与真实前沿的距离尽可能小[8]。
3 拟态物理学算法
拟态物理学是由美国怀俄名州立大学的Spear等人以自然物理力量为动力所提出来的，受启发于物理力学定律，故称之为“拟态”，最初被应用于机器人群体的分布式控制研究。
文献[5]中利用拟态物理学的思想处理全局优化问题，提出一种新的优化算法APO。该方法本质上是通过对牛顿第二定律的模拟，将问题可行域中采样的个体粒子看作是物理中的质点，在粒子之间虚构一种相互作用力，给出该作用力的计算规则，结合个体的速度和位移更新粒子的运动，实现对整个种群的控制。
在APO算法中，种群规模为N，个体i在k维上的速度表示为vi,k，个体i在k维上位置表示为xi,k。第i个个体的质量为mi(i=1,2,…,n)，且个体的质量随着适应值的改变而改变，适应值小的个体质量大，适应值大的个体质量反而小。此外个体质量还受到种群中适应度值最优和最差个体的影响，进而影响个体所受的力。所以，计算个体质量函数定义为：
(3)
式中：xbest、xworst分别表示最优、最差适应值的位置。
由每个个体质量，仿照万有引力公式，计算出各个个体所受作用力的大小，故个体在第k维上所受的第j个(i≠j)个体的虚拟力为Fij,k：
(4)
个体在第k维上所受作用力合力的计算公式为：
(5)
式中，G为引力因子。
个体i的速度和位置更新公式为：
vi,k(t+1)=wvi,k(t)+λFi,k/mi, ∀i≠best
(6)
xi,k(t+1)=xi,k(t)+vi,k(t+1), ∀i≠best
(7)
式中：w∈(0,1)为惯性权重；λ为分布在区间[0，1]上的一个随机数；xi∈[xmin,xmax]； vi∈[vmin,vmax]。
4 R2-ICRMOAPO算法
随着约束多目标问题的广泛应用，相应的约束多目标优化算法也层出不穷。非支配解的选取是解决约束多目标问题的重要技术之一，直接影响算法的性能。ICRMOAPO 算法通过比较新的非支配个体和archive集中的个体来获得新的外部存储集，在一定程度上对问题进行了求解，但由于没有严格地对解集进行筛选，可能会使一些优质个体缺失。R2指标可以很好地解决这一问题，利用个体的R2指标贡献值，对个体进行筛选，可以得到一组性能较好的解集。
4.1 R2指标
R2指标是基于效用函数，区分候选解的优劣，从而选择效用较大的候选解，最初被用于评价2组个体的相对质量。利用具有特定参考点的标准加权切比雪夫函数，可以用该指标相对于参考点来评估单个个体集合的质量。因R2指标可以综合评价种群的收敛性和多样性,并可以通过快速计算获得均匀分布的解集,已被应用于求解多目标优化问题[14-15]。
在此利用加权切比雪夫函数作为R2指标的效用函数，通过一组个体(集合A)和参考点z*对个体的质量进行评价：
(8)
式中：W为一组均匀的权重向量，每个权重向量μ=(μ1， μ2，…， μm)∈W，均匀的分布在目标空间中，1/|W|代表每个权重向量被选择的概率。
如果效用函数集中每个函数由于不同的权重向量而确定,则个体的R2指标的贡献值被定义为：
CR2(x,A,W,z*)=R2(A,W,z*)-R2(A\{x},W,z*)
(9)
通过R2指标的贡献值可以对种群中的个体进行全排序,R2指标的贡献值越大说明解集的分布性能越好，解的质量越高，将其保留到下一代种群的概率就越大。
4.2 R2-ICRMOAPO算法基本思想
外部存储集非支配个体的选取对算法的收敛性和分布性具有直接的影响，在算法的迭代进化过程中，非支配解的数量也在同步增加，并且可能会超过预先设定的解集规模，因此外部存储集的更新和维护对于解集的质量起着非常重要的作用。本研究在ICRMOAPO算法的基础上，将非支配排序和R2指标结合，首先利用Pareto支配关系得到一组非支配解集，计算出非支配解集中个体的R2指标贡献值，然后将R2指标贡献值作为外部存储集的更新机制，通过删减R2指标贡献值较小的个体，实现对整个种群的控制，从而选择出效用更好的候选解。针对上述思想构造出一种R2-ICRMOAPO算法，为求解约束多目标优化问题提供一种新的思路。算法的关键要素如下：
1) 个体序值
设在第t代生成的种群P(t)由N个个体所组成，ni(t)表示种群P(t)中支配个体i的所有个体数量，则个体i在第t代的序值定义为：
ri(t)=1+ni(t)
(10)
依据个体序值的定义，为种群中所有的支配个体，分配相对应的序值。但若在某一代中出现具有相同序值的个体，则需要从其他角度对此问题进行解决。
2) 邻域半径
针对序值相同的个体，下面引入邻域半径的定义。
定义个体的ε邻域为该个体半径ε内所包含的区域。为了避免算法中邻域半径过小，故将邻域半径定义为：
(11)
式中，max|Ei|为初始化时第i维上个体间的最大欧式距离。
3) 邻域拥挤度
邻域拥挤度是指在个体间支配与被支配的作用下所计算出某个个体在对应邻域半径区域内所包含个体数量的大小。
当出现序值大小相同，邻域拥挤度也相同的个体时，将邻域半径重新调整。若调整后，具有相同序值的个体在相应邻域半径内还包含相同的个体数量，将这些个体随机进行排列，分配编号即可。
4) 质量函数
为了考虑个体间的支配关系和个体所处邻域半径内的拥挤程度，并充分体现多目标优化问题本身的特点，将质量函数定义为[8]：
(12)
5) 引力因子
引力因子是影响算法性能的重要因素之一，在一定程度上会决定种群做收敛运动的个体数量。在标准APO算法中，将引力因子G规定为10。R2-ICRMOAPO算法为了使引力因子G随种群的进化而产生动态的变化，将其进行改进，令G为：
(13)
式中：Gi和Gf为引力因子G的初始值和终值，由经验可知，算法中Gi=70，迭代到最后，G的终取值为Gf=1；iter为迭代次数；Maxit为最大迭代次数。
6) 惯性权重
惯性权重的选取影响算法的性能和效率，迭代初期较大的惯性权重使算法保持较强的全局搜索能力，迭代后期较小的惯性权重有利于算法进行更精确的局部搜索。
在R2-ICRMOAPO算法中采用线性递减的惯性权重，使惯性权重随着迭代的进行而动态减小：
(14)
式中：wi为初始值，算法中取0.9；wf为终值，算法中取0.4。
7) 不可行度
不可行度可以作为一种准则来判断更新个体是否落入可行域中。
R2-ICRMOAPO算法中定义不可行度为：
(15)
式中: gi、hp分别为多目标优化问题的不等式约束和等式约束。当xi为可行解时，不可行度为0。
4.3 R2-ICRMOAPO算法流程
R2-ICRMOAPO算法流程如图1所示。
图1 R2-ICRMOAPO算法流程框图
Fig.1 R2-ICRMOAPO algorithm flowchart
R2-ICRMOAPO算法的主要步骤如下：
Step 1：初始化种群：种群规模N，外部存储集规模archivesize，权重向量W，最大迭代次数Maxit，随机产生所有个体的初始速度和位置，初始化参考点。
Step 2：计算所有个体在每个目标下的适应度值，确定非支配个体，存储于archive集。对非支配个体的序值定义为1，其余个体的序值定义为ni+1。
Step 3：对archive集中的非支配个体计算R2指标贡献值，按降序原则进行排序。
Step 4：根据个体序值的定义，对种群中所有个体按照升序原则进行排序。
Step 5：利用欧式距离计算邻域半径，若出现序值相等的个体，重新确定个体邻域拥挤度，按升序原则对全部个体重新进行排序。
Step 6：出现邻域半径内包含相同个体数量的同序值个体时，将该个体标记为flag(i)=1。
Step 7：flag(i)=1时，对邻域半径进行调整并重新排序，对于调整后序值大小和拥挤度相等的个体，进行随机排序。对所有个体赋予1～N的自然数，作为每个个体序号，记作ranki(t)。
Step 8：计算每个个体的质量mi，个体在第k维上所受的第j个(i≠j)个体的虚拟力Fij,k，在第k维上所受作用力的合力Fi,k。
Step 9：根据式(6)、式(7)更新个体的速度和位置。
Step 10：利用不可行度的判断准则来确定所更新的个体是否落在可行域内，将未落入可行域内的个体用最优迭代更新的方法重新更新其位置。
Step 11：确定新的非支配个体，更新archive集，当archive集中都是非支配个体时，对所有个体计算CR2值，按照降序原则进行排序，若个体数量超过archive集规模时，将archive集中R2指标贡献值较小的个体删除，直到个体数量满足设定的规模，获得新的archive集。
Step 12：若达到最大迭代次数则退出，否则返回步骤5。
5 实验结果与分析
为了评估R2-ICRMOAPO算法在解决约束多目标问题时的性能，通过R2-ICRMOAPO算法和其他四种优化算法：基于序值与拥挤度的拟态物理学多目标算法(improved constrained rank multi-objective artificial physics optimization，ICRMOAPO)、带精英策略的非支配排序的遗传算法 (non-dominated sorting genetic algorithm，NSGAII)、基于R2指标的多目标粒子群算法(multi-objective particle swarm optimization on R2 indicator,R2MOPSO)、ε多目标进化算法(epsilon multi-objective evolutionary algorithm,ε-MOEA)在标准约束多目标测试函数上进行实验对比，并将超体积(Hypervolume,HV)[16]和反转世代距离(inverted generational distance，IGD)[17]作为算法的性能评价指标，从收敛性和分布性对算法进行综合评价，全面说明算法的有效性。
本文选定的测试函数如表1所示。
表1 约束多目标标准测试函数
Table 1 Constrained multi-objective standard test functions
测试函数变量范围目标函数约束域Binh(2)x1∈[0,5]x2∈[0,3]min f1(x1,x2)=4x21+4x22min f2(x1,x2)=(x1-5)2+(x2-5)2(x1-5)2+x22-25≤0-(x1-8)2-(x2+3)2+7.7≤0SRNxi∈[-20,20]i=1,2min f1(x1,x2)=(x1-2)2+(x2-1)2+2min f2(x1,x2)=9x1-(x2-1)2x21+x22-225≤0x1-3x2+10≤0DEBx1∈[0.1,1.0]x2∈[0,5]min f1(X)=x1min f2(X)=(1+x2)/x1-9x1-x2+6≤0-9x1+x2+1≤0
实验中，设定种群规模和外部存储集规模均为100，进化代数为50。图2～图4表示算法R2-ICRMOAPO、ICRMOAPO、R2MOPSO、NSGAII、和ε-MOEA对测试函数生成的Pareto前沿和真实Pareto前沿曲线。
图2 Binh(2)的Pareto前沿曲线
Fig.2 The Pareto Frontier of Binh (2)
图3 SRN的Pareto前沿曲线
Fig.3 The Pareto Frontier of SRN
图4 DEB的Pareto前沿曲线
Fig.4 The Pareto Frontier of DEB
图5～图7显示了R2-ICRMOAPO算法所求得Pareto 解集中更新个体的位置。
图5 Binh(2)的Pareto解集中更新个体位置曲线
Fig.5 The updated position of Binh(2)
图6 SRN的Pareto解集中更新个体位置曲线
Fig.6 The updated position of SRN
图7 DEB的Pareto解集中更新个体位置曲线
Fig.7 The updated position of DEB
大多数研究使用逼近实际解决方案的程度和设置解决方案的分布情况来测量多目标优化算法的性能。本文采用HV和IGD算法性能指标对5种算法进行对比。
HV指标是一种能够同时衡量算法的收敛性和分布性的综合性能指标，表示非支配解集覆盖的目标空间区域大小。HV函数的定义为：
(16)
式中： δ为Lebesgue测度,用来测量体积；|S|为非支配解的数目；vi为参照点与解集中第i个解所构成的超体积。算法求得的HV值越大则意味着解集收敛性和分布性越好。
IGD指标是一种综合评价解集质量的指标，衡量的是真实前沿的个体到所求得的近似解集之间最小距离的平均值。IGD的函数定义为：
(17)
式中：
为目标空间上每一个个体与真实Pareto前沿中最近个体之间的欧氏距离；
为目标空间中个体的数量。IGD值越小则意味着找到的解集收敛性和分布性能整体效果越好，也说明Pareto解集更加近似于真实解。
对测试函数进行了100次独立实验后，表2～表3分别展示了5种算法关于HV和IGD性能指标的平均值、最优值和方差。此外，还采用统计学Wilcoxon秩和检验来表明结果的显著性水平，显著性差异水平为0.05，如果P值大于0.05，则表示2种算法之间没有明显差异[3]。
表2 测试函数的HV指标值
Table 2 HV performance indicators on test functions
测试函数R2-ICRMOAPOICRMOAPONSGAIIε-MOEAR2MOPSOBinh(2)均值3.824E-013.475E-013.507E-013.028E-012.580E-01方差5.237E-025.942E-023.515E-025.543E-017.928E-02P2.518E-063.712E-062.368E-091.722E-07SRN均值6.974E-016.860E-014.943E-011.021E0-23.424E-01方差8.023E-028.738E-023.529E-029.763E-021.471E-01P2.748E-036.285E-054.064E-053.002E-07DEB均值3.855E-031.721E-057.198E-032.365E-039.842E-04方差5.747E-022.648E-041.065E-048.846E-042.314E-03P3.049E-032.781E-014.001E-053.616E-08
表3 测试函数的IGD指标值
Table 3 IGD performance indicators on test functions
测试函数R2-ICRMOAPOICRMOAPONSGAIIε-MOEAR2MOPSOBinh(2)均值2.937E-022.963E-022.962E-023.452E-025.031E-02方差4.872E-035.341E-035.303E-032.278E-027.835E-03P2.940E-043.031E-036.342E-084.232E-08SRN均值5.236E-025.264E-026.052E-021.419E-017.314E-02方差7.424E-036.943E-032.379E-021.908E-022.389E-02P7.908E-023.028E-052.973E-036.326E-09DEB均值7.322E-031.486E-025.818E-037.494E-032.771E-02方差3.262E-052.937E-032.685E-052.619E-045.159E-05P5.322E-072.915E-015.367E-024.685E-06
从表2和表3可以看出，R2-ICRMOAPO算法的HV值高于ICRMOAPO算法，且IGD值低于ICRMOAPO算法。数据表明，利用R2指标贡献值对个体进行筛选，提高了算法的收敛性和分布性。
从图2～图4和表2～表3可知，R2-ICRMOAPO算法在测试函数Binh(2)和SRN上对Pareto前沿产生了更好的近似，得到的解大部分够收敛于真实的Pareto前沿上，在DEB函数上略差于NSGAII算法。以上数据和图表表明，R2-ICRMOAPO算法在大多数测试问题上优于对比算法，说明R2-ICRMOAPO算法在解决约束多目标问题方面具有很大的潜力，可进一步推广应用于实际问题中。
6 结论
在APO算法的基础上，针对约束多目标优化算法，构造基于R2指标的多目标约束优化算法，将非支配排序和R2指标结合，利用R2指标贡献值区分候选解的优劣，更新最优个体，维护外部存储集；将其与4种算法在HV和IGD指标下进行对比实验。主要得到以下结论：
1) 算法利用个体间的支配关系，考虑个体的拥挤度，对搜索空间中个体适应值合理排序，增强了Pareto 解集的多样性。
2) 设计R2指标贡献值选择个体，对外部存储集进行维护和更新，增强了种群的分布，在处理约束多目标问题中取得较好的效果。
3) 该算法在求解约束多目标的问题中收敛性和分布性均得到提高，但对于提高收敛速度还需进行研究。
[1] Tian Y,Cheng R,Zhang X,et al.An indicator based multi-objective evolutionary algorithm with reference point adaptation for better versatility[J].IEEE Transactions on Evolutionary Computation,2017,22(4):609-622.
[2] Coello C,Lechuga M S.MOPSO:A proposal for multiple objective particle swarm optimization[C]//Congress on Evolutionary Computation.IEEE Service Center,2002:1051-1056.
[3] Deb K,Pratap A,Agarwal S,et al.A fast and elitist multiobjective genetic algorithm:NSGA-II[J].IEEE Transactions on Evolutionary Computation,2002,6(02):182-197.
[4] Wang W X,Li K S,Tao X Z,et al.An improved MOEA/D algorithm with an adaptive evolutionary strategy[J].Information Sciences,2020,539:1-15.
[5] Xie L P,Zeng J C,Formato R A.Convergence analysis and performance of the extended artificial physics optimization algorithm[J].Applied Mathematics and Computation,2011,218(8):4000-4011.
[6] 王艳.多目标拟态物理学优化算法及其应用研究[D].兰州:兰州理工大学,2011.
Wang Y.Research on Multi-objective Artificial physics optimization algorithm and its application[D].Lanzhou:Lanzhou University of Technology,2011.
[7] Yan W,Zeng J,Ying T.An artificial physics optimization algorithm for multi-objective problems based on virtual force sorting proceedings[C]//International Conference on Swarm.Springer Berlin Heidelberg,2010.
[8] 孙宝,孙大刚,李占龙,等.基于序值与拥挤度的拟态物理学多目标算法[J].系统工程与电子技术,2014,36(12):2442-2448.
Sun B,Sun D G,Li Z L,et al.Artificial physics multi-objective algorithm based on sequence value and crowding degree[J].Systems Engineering and Electronics,2014,36(12):2442-2448.
[9] Sedighizadeh D,Masehian E,Sedighizadeh M,et al.GEPSO:A new generalized particle swarm optimization algorithm[J].Mathematics and Computers in Simulation,2021,179:194-212.
[10] Cai D,Wang P Y,Miao Y.A new multi-objective particle swarm optimization algorithm based on decomposition[J].Information Sciences,2015,325(1):541-557.
[11] Kohler M,Vellasco M,Tanscheit R.PSO+:A new particle swarm optimizetion algorithm for constrained problems[J].Applied Soft Computing,2019,85:105865.
[12] Garcia I C,CoelloCoello C A,Arias-Montano A.MOPSOhv:A new hypervolume-based multi-objective particle swarm optimizer[C]//Evolutionary Computation.IEEE,2014.
[13] Li F,Liu J C,Tan S,et al.R2-MOPSO:A multi-objective particle swarm optimizer based on R2-indicator and decomposition[C]//Evolutionary Computation.IEEE,2015.
[14] Wei L X,Li X,Fan R,et al.A hybrid multi-objective particle swarm optimization algorithm based on R2 indicator[J].IEEE Access,2018,6(99):14710-14721.
[15] Chabane B,Basseur M,Hao J K.R2-IBMOLS applied to a practical case of the multiobjective knapsack problem[J].Expert Systems with Applications,2016,71(APR.):457-468.
[16] Zhou H,Wang L P,Ling M,et al.New research of updating neighborhood strategy in MOEA/D[J].Journal of Chinese Computer Systems,2017,38(4):852-856.
[17] Cheng R,Jin Y,Olhofer M,et al.A reference vector guided evolutionary algorithm for many-objective optimization[J].IEEE Transactions on Evolutionary Computation,2016,20(05):773-791.