基于改善解集分布性的多目标优化

王学武 方俊宇 高进 顾幸生

（华东理工大学 能源化工过程智能制造教育部重点实验室，上海 200237）

在实际生产过程中，因客观需求，往往存在2个或者2个以上的优化目标，这类优化问题统称为多目标优化问题（Multi-Objective Problem，MOP）。与单目标优化不同，多目标优化的最优解并不唯一：这些解无法在所有目标函数值均保持最优，且不存在所有目标函数均优于自身的解——非支配解，这些解构成的解集称为帕累托解集（Pareto Set，PS），而解集对应的决策空间目标函数的点则构成帕累托前沿（Pareto Front，PF）。在过去的数十年中，随着计算智能的发展，各种类型的多目标优化算法（Multi-Objective Algorithm，MOA）被开发出来解决不同类型的MOP，均取得了较好的效果，在解集充分逼近真实PF 的同时保证了解的多样性，给决策者提供了多样化的选择。

多目标优化算法中，常见的是基于进化计算的多目标进化优化算法（Multi-Objective Evolutionary Algorithm，MOEA）。按进化机制的不同，MOEA可分为以下3类：①基于分解的MOEA 算法，其代表有MOEA/D［1］、MOEA/D-RWV［2］、R-MOEA/D［3］等；②基于非支配排序的MOEA 算法，其代表有NSGAⅢ［4］、ASDNSGA-Ⅱ［5］等；③基于指标的MOEA 算法，其代表有IBEA、AR-MOEA［6］等。近年来，随着计算机算力的提升，除了各种新的环境选择策略被提出，多名学者还将原有不同类型的策略相结合，甚至将机器学习融入基础策略提出新的算法，在不同类型的MOP 上提升解集的收敛性和分布性：Tian 等［7］利用两个种群分别解决带约束多目标优化问题（Constrained Multi-Objective Problem，CMOP）和无约束多目标优化问题，并设置协同进化框架，有效提升了算法在具有狭小可行域的CMOP上解集的分布性；Panichella［8］设计了一种自适应的、兼顾收敛性和分布性的指标，并将其与非支配排序策略相结合，提出了AGE-MOEA 算法，大大提升了在不规则PF 的MOP 获得解集的分布性；Sonoda 等［9］则采用支持向量机（Support Vector Machine，SVM）生成代理人模型去适应MOEA/D通过分解策略产生的子问题，从而在高维多目标优化问题（Many-Objective Optimization Problem，MaOEA）上获得了分布性较好的解集。

在基于非支配排序策略的多目标优化算法中，NSGAⅡ的基本框架（非支配排序+环境选择策略）成为该类算法的代表。非支配排序用于寻找当前最优前沿解集，环境选择策略则用于提升分布性。SPEA2则通过强壮支配指数［10］搜寻非支配个体，并将其与基于k-近邻聚集密度的分布性指数加权求得个体适应度值，再通过设置适应度阈值进行筛选的方法同时保证收敛性和分布性。SPEA2基于个体适应度选拔父代，这使得通过遗传算法生成的新个体有更大概率同时具有良好的收敛性和分布性，因此其在大多数MOP上的收敛速度要比NSGAⅡ快。然而，SPEA2 存在如下3 个缺陷：①在保证收敛速度和全局搜索能力的同时，局部搜索能力显得有些不足；②受其档案截断策略的影响，种群中非支配个体快速增多，不同个体间的适应度差异也迅速变小，这使得后期选择的父代无法保证基因多样性，通过遗传算法获得更优子代的概率大幅降低，种群进化容易过早陷入停滞；③对于存在不规则PF 的MOP（如IMOP7、ZDT3、VNT1等），SPEA2算法获得的解集分布性有待进一步提高。针对缺陷①，翁理国等［11］提出了一种通过额外设置一个用于存放局部搜索得到的非支配解的归档集的方法，增强了原算法的局部搜索能力。针对缺陷②和③，文中提出一种改进的算法CM-SPEA2（Clustering-Based Modified Strength Pareto Evolutionary Algorithm 2），将层次聚类方法融入到SPEA2的环境选择策略中，并引入改进杂乱度指数参与环境筛选策略。

1 CM-SPEA2的算法原理

1.1 算法总框架

CM-SPEA2 整体沿用SPEA2 的主体框架（即通过遗传算法产生子代、通过强壮支配指数和表征分布性加权得到的适应度值进行环境筛选），只保留适应度“最优”的一半子代个体与原种群合并。

算法的主体框架伪代码如下：

1.2 CM-SPEA2的适应度函数

与多数直接将聚类应用于环境选择的机制不同，CM-SPEA2 将聚类用于计算来改进适应度函数，再与基于适应度大小的环境选择策略结合起来，间接依照杂乱度进行筛选并得到分布性较好的解集。

SPEA2算法的适应度函数如式（1）所示：

式中，R（i）为支配个体i的所有个体强度之和，

其中，S（j）为个体j的支配强度，表示受该个体支配的解的数量，

Pt为第t代的种群，At为第t代产生的归档集。

原始适应度赋值过程需引入密度信息来区分具有相同原始适应度值的个体。SPEA2采用k-近邻方式来计算个体密度值D(i)：

式中，为个体i和第k个相邻个体的距离，其中

N为种群个体数，为归档集个体数。

然而，作为适应度函数中表征种群分布性的算子，D（i）无法区分聚集密度相同、分布均匀程度存在差异的个体；与之相比，杂乱度虽不能表征个体间的聚集密度，但可以衡量聚集密度及相邻个体间距的均匀程度。因此，基于种群的杂乱度分析，CM-SPEA2引入一种新的表征分布性的算子O（i）。

设算法通过聚类得到聚集簇C1，内部包含N个个体；又设第j个个体的杂乱度最低，即Chaos(j)=min{m(i)}，i=1，2，…，N。则聚集簇C1中第i个个体与第j个个体的Manhattan距离D（i）为

式中，C1_O(i)为第i个个体的目标函数向量值。

对D（i）进行升序排序，得Rank（i）。对Rank（i）进行量纲处理，最终得到个体i的分布性指数O(i)：

CM-SPEA2 算法的个体i的适应度如式（8）所示：

1.3 聚类方法

对于不规则PF的MOP，仅凭单一指标（如拥挤度距离、聚集密度等）通常很难反映个体在种群的真实分布情况，需要借助参考点，而参考点的选择很大程度上会影响环境选择策略的优劣。聚类是一种有效寻找参考点的方法，能够适应种群最佳前沿的聚集特性，从而结合相关指标改善分布性。

聚类方法有很多。王学武等［12］采用DBSCAN密度聚类将种群分为不同簇，并挑选领导粒子作为参考，基于拥挤度距离完成环境选择。该方法鲁棒性强，受离群点影响较小；然而，相关参数的设定（邻域半径、最少邻域点数量等）对聚类效果的影响非常大，且不适合形状不规则PF 的多目标优化问题。Hua等［13］在CAMOEA算法中对种群个体进行非支配排序之后，将最劣解个体数作为簇数量上限进行层次聚类（如图1 所示），再计算各个簇的形心，将其作为参考点。该方法充分利用了最劣解集个体多样性最好的特点，能够自适应地根据PF 的形状设置参考点的数量和位置。

图1 层次聚类效果图Fig.1 Effect diagram of hierarchical clustering

受文献［12-13］启发，文中提出一种新的聚类方法，其特点如下：①聚类的对象是种群中的非支配个体（即强壮支配指数为零的个体）；②簇数量最大不超过最劣解的个数；③聚类方法采用不等权重算术平均聚类（Weighted Pair Group Method with Arithmetic Mean，WPGMA）；④个体间的相似度采用带权算术平均Manhattan 距离来代替；⑤聚类完成后，计算个体在所属簇的杂乱度值，选择杂乱度最小的个体作为参考点。算法效果如图1 所示（图中坐标轴f1和f2分别代表目标值1和目标值2）。

1.4 改进种群杂乱度分析

对于非均匀PF 的MOP，传统的种群维护方法容易导致解集分布在PF 的某个集中区域，从而破坏解集分布性。为更好地表征个体在解空间分布的差异性，引入个体杂乱度分析。

对于种群中个体的杂乱度m(i)，李密青等［14］给出了如下定义：

式中，d（i）为个体i在种群生成的一棵欧式最小生成树（Euclidean Minimum Spanning Tree，EMST）中的度数（即与该个体连接的边数），max（l（i））以及min（l（i））分别为EMST 中连接个体i的最长边和最短边。通过该定义，可以轻易找到种群内部杂乱度最大的个体。然而，该计算方法存在以下缺陷：①每次计算都要预先生成EMST，计算消耗大；②只能准确衡量杂乱度最大的个体，不能准确衡量杂乱度最小的个体，尤其是离群边缘个体（如图2所示），此时d(i)=1，而max (l(i))=min (l(i))，从而得到该点的杂乱度m(i)=0，杂乱度达到极小值，这明显与实际不符。

图2 含离群点的EMST示意图Fig.2 Diagrammatic sketch of EMST containing an outlier

此外，随着迭代的进行，个体间距越来越小，使用这种杂乱度计算方法会导致后期解集内不同个体之间的差异越来越小，从而使得环境选择策略难以进一步提升分布性和收敛性。故对式（9）进行如下改进：对于一个特定种群P，用个体i与最近个体的距离和个体间平均最短间距的差异来表征其杂乱度，即

式中，δ（i）为第i个个体与种群P内其他个体的最小Manhattan距离（同样可以通过EMST计算）。对于存在不均匀PF 的MOP，式（10）往往不能很好地表征个体的杂乱度，需要通过聚类方法确定若干个簇，再计算个体在簇内的杂乱度。设个体i通过之前的聚类方法被分在了簇Cj中（Cj中总共有Nj个体），则该个体的杂乱度为

前人的研究已经证实，直接淘汰杂乱度大的个体并不一定能够改善种群分布的均匀特性，尤其当杂乱度大的个体离种群边缘个体较远时，单独删除杂乱度最大的个体会有很大几率增加种群分布的不均匀程度。图3（a）所示为原种群分布图，为方便讨论，除E、F间距设置偏小外，图中其余相邻点（个体）的间距均相同。按照式（9），点F的杂乱度m(F)=LFGLEF-12 ＞0.5，点G的最长边、最短边与点F的相同，故有m(F)=m(G)；而点B、C、D的杂乱度均为0.5，点A、H的杂乱度均为0，故点F、G为杂乱度最大的个体。为此，郑金华等［15］提出删除距杂乱度最大个体最近的个体来改善分布性。图3（b）-3（e）分别为删除点F、G及其各自的近邻个体E、H后的种群分布。

图3 不同筛选策略的对比Fig.3 Comparison of different selecting strategies

可以明显看出，删除个体E或H后的种群分布性优于删除个体F或G后。然而，该方法存在如下缺陷：①一次只能删除1个个体，否则杂乱度会增加；②每次删除操作完成之后需要重新评估杂乱度再重复删除操作，大大增加了计算消耗。为此，笔者提出一种新的策略：由删除杂乱度最大个体的近邻个体转化为保留杂乱度最小个体及其近邻个体。该策略在很大程度上减少了计算消耗，并且由于SPEA2 本身采用的也是设置阈值保留个体的策略，沿用原有环境选择策略框架便可以实现改进策略（当出现2个及2个以上的杂乱度最低的个体时，则选择与杂乱度最大个体间隔最远的个体作为参考点，再保留近邻个体）。同样以图3 为例，由于点A、B、C、D、E的间距设置为相等，故此处将点A作为参考点，在保留近邻7个个体后，得到的种群与删除个体H后的种群分布相同。

新策略以牺牲边界个体为代价，在一定程度上削弱了种群多样性；然而，在迭代早期，非支配个体往往不足以填满整个归档集，此时会接纳一定数量的其他个体，这将有一定几率选择种群或簇的边缘个体，解的多样性在一定程度上得到补偿。

CM-SPEA2 的环境选择策略步骤如下：①初始化种群Pt（此处Pt由原种群和通过遗传算法获得的较优后代合并而来），设置归档集At（大小为N），初始为空集；②计算Pt个体的适应度；③筛选Pt中适应度小于1 的个体，将其纳入归档集；④判断At是否满集，如已满则结束，否则转至⑤；⑤将其他个体按照聚集密度（此处用Manhattan距离代替SPEA2的Euclidean 距离）进行升序排序；⑥依次序将档外个体纳入At直至满集。

1.5 时间复杂度分析

设MOP的目标数为M，种群个数为N。单次迭代过程中，适应度函数的计算可以划分为以下几个步骤：①计算个体的强壮支配指数（记为T1(N)）；②层次聚类（记为T2(N)）；③计算个体杂乱度（记为T3(N)）；④计算分布性指数（记为T4(N)）。

综上所述，单次迭代过程中，计算所有个体的适应度函数的时间复杂度T(N)=T1(N) +T2(N) +T3(N) +T4(N)=O(N3)，可以看出，CM-SPEA2 的时间复杂度与SPEA2的相同。

2 算法仿真及性能对比

实验测试环境如下：Intel（R）Core（TM）i5-9300H CPU @2.40 GHz 2.40 GHz，内存16.0 GB，基于x64 处理器的64 位操作系统Window 10 中文版，软件版本为Matlab R2021a，平台版本为PlatEMO［16］3.4。

为了验证改进后算法的性能（尤其是对不规则PF 的MOP 性能）提升，选择IMOP、ZDT 和VNT 类问题作为测试MOP。IMOP 类问题的a1参数设置为0.4，K参数设置为5。参与性能对比的算法为：原算法SPEA2、基于非支配排序的算法NSGAⅡ、同样采用聚类生成参考点的算法CAMOEA 及在环境选择过程中自适应生成参考向量并致力于改善种群分布性的算法RVEAa［17］。用这些算法求解3类问题（IMOP、ZDT 和VNT）并比较其性能。考虑到RVEAa 算法因基于参考向量的环境选择策略对种群规模有特殊要求，初始种群大小统一设置为105；针对不同类型MOP的收敛速度，将IMOP类问题最大函数评价次数设置为42 000，其余MOP 设置为31 500。由于ZDT5 属于二进制类MOP 问题，需要采用与其他MOP 不同的方式衡量其性能指标，因此其不参与仿真及后续性能的比较。

为避免解集进化过慢或退化为随机搜索，所有算法的遗传算子将交叉概率统一设置为0.5，变异概率p=1/D（D为目标数）。

2.1 性能指标对比

用于比较的指标有：反向世代距离（Inverse Generation Distance，IGD）指标DIGD［18］以及分布均匀程度（Spacing）指标S［19］。

IGD 指标通过计算真实PF 解集中的个体和算法最终得到的最佳前沿的个体的平均Euclidean 距离而获得，其计算公式如式（12）所示：

式中，x*为当前Pareto 解集中的个体，x为问题实际的Pareto 前沿个体，|PS|为Pareto 解集中个体的个数。由式（12）知：IGD 指标越小，解集的综合性能越好。

Spacing 指标通过计算种群个体之间的邻近距离的标准差来表征解集的分布性（均匀程度），其计算公式如式（13）所示：

考虑到MOEA的随机性，算法在各个MOP上单独运行30次，然后计算各个算法获得解集指标的均值和标准差。为更直观对比各个算法在不同MOP上的指标优劣，引入rank-sum 检验（显著水平0.05），检验结果用“+/-/=”表示（三者对应的数值分别表示参与比较的算法与CM-SPEA2 算法相比“更好”“更差”及“相当”的问题的数量）。算法在不同MOP下的IGD和Spacing指标分别见表1和表2。

表1 不同算法在不同MOP下的IGD指标对比1）Table 1 Comparison of IGD of different algorithms in multiple MOPs

由表1 和表2 结果统计可知：在所有测试MOP中，IGD、Spacing指标占优势的MOP分别占总测试MOP的68.75%、87.5%；与SPEA2相比，CM-SPEA2算法的IGD、Spacing 指标提升的MOP 数量则分别占所有MOP 的68.75%、87.5%。对于部分测试问题，SPEA2与NSGAⅡ、CAMOEA、RVEAa相比不占优势的指标经改进后也占优势（如对ZDT6的Spacing指标）。对于规则PF 的MOP（如IMOP1、IMOP2、ZDT1、ZDT2等），在保证和SPEA2的IGD 指标相当的情况下，CM-SPEA2 算法的Spacing 指标基本上大幅提升。为方便讨论，文中将不规则PF的MOP分为两类：不连续PF、连续但形状不规则PF。对于不连续PF 的MOP（如IMOP3、IMOP5、ZDT3 等），CM-SPEA2算法的IGD、Spacing指标占优势的比例分别为75%、100%；对于连续但形状不规则PF的MOP（如IMOP4、IMOP6、VNT1、VNT2 等），CM-SPEA2的IGD、Spacing 指标占优势的比例分别为83.3%、83.3%。与SPEA2 相比，对于不连续PF 的MOP，CM-SPEA2 算法的IGD、Spacing 指标有提升的比例分别为75%、100%；对于连续但形状不规则PF 的MOP，CM-SPEA2 算法的IGD、Spacing 指标有提升的比例则分别为100%、83.3%。这充分表明：文中算法对解集分布均匀程度的改善在一定程度上促进了综合性能的提升，其中，对Spacing 指标的提升最为显著，最高可达16.3%，IGD 指标的提升也普遍达到5%以上。

2.2 最优前沿分布

为了更加直观地比较不同算法获得的最优前沿，此处分别选择CM-SPEA2、SPEA2 和相同采样聚类方法的CAMOEA 算法来测试不连续PF 的MOP（ZDT3和DTLZ7）以及连续但形状不规则PF 的MOP（IMOP7）。除了IMOP5、IMOP7 的最大函数评价次数设置为42 000 外，其余MOP 统一设置为31 500。最终得到的前沿分布如图4-6所示，图中坐标轴f1、f2和f3分别代表目标值1、2和3。

图4 CM-SPEA2、SPEA2 和CAMOEA 在ZDT3 问题上获得的最优前沿的比较Fig.4 Comparison of the optimal fronts among CM-SPEA2，SPEA2 and CAMOEA in ZDT3

图5 CM-SPEA2、SPEA2 和CAMOEA 在IMOP7 问题上获得的最优前沿的比较Fig.5 Comparison of the optimal fronts among CM-SPEA2，SPEA2 and CAMOEA in IMOP7

图6 CM-SPEA2、SPEA2 和CAMOEA 在DTLZ7 问题上获得的最优前沿的比较Fig.6 Comparison of the optimal fronts among CM-SPEA2，SPEA2 and CAMOEA in DTLZ7

通过对比得知：①经过相同次数的迭代后，CM-SPEA2 的最佳前沿能充分逼近问题的实际PF，丝毫不逊色于SPEA2和CAMOEA；②即便CM-SPEA2在环境选择策略中牺牲了种群（簇）的边缘个体，最终获得的解集依然能够充分覆盖整个PF，解集多样性得到了保证。

2.3 策略可行性验证

为进一步探讨接纳参考点附近个体的环境选择策略的性能优劣，分别采用离参考点由近及远（简称策略1）和由远及近（简称策略2）保留个体的环境选择策略来解决多个MOP，得到的IGD 和Spacing指标如表3所示。

表3 不同环境选择策略的CM-SPEA2在多个MOP下的IGD和Spacing指标对比Table 3 Comparison of IGD and Spacing in multiple MOPs for CM-SPEA2 using different environment selection strategies

经对比可知：对不同MOP，采用策略1 时CM-SPEA2 的IGD 指标与采用策略2 时没有显著优劣差别，而Spacing指标占优势的则有9个。这充分表明，删除边远个体策略不仅未明显削弱种群的综合性能，还可以保证种群分布足够均匀，这进一步验证了该策略的可行性。

2.4 演化过程分析

为了体现CM-SPEA2改进策略对算法演化过程的影响，将前述算法在VNT2 问题上进行IGD 性能比较，得到的IGD指标演化曲线如图7所示。

图7 不同MOEA算法在解决VNT2问题时的IGD指标演化曲线Fig.7 IGD evolution curves of different MOEAs when solving VNT2

IGD 指标演化曲线充分展示了不同算法的性能指标提升过程。由图7 可以发现：不改变原SPEA2算法的环境选择策略能够更好地保证算法的收敛速度，且及时在迭代后期进化较为缓慢的阶段；改进后的CM-SPEA2 算法的IGD 指标波动很小，且没有出现类似RVEAa算法指标明显变差的情况。

此外，由图8所示在ZDT3问题上的Spacing指标演化曲线可以看出：迭代早期CM-SPEA2 比SPEA2的Spacing指标下降更快，迭代中后期Spacing指标未出现明显上升趋势，这表明解集分布的均匀程度能快速提升并达到稳定。

图8 不同MOEA算法在解决ZDT3问题时的Spacing指标演化曲线Fig.8 Spacing evolution curves of different MOEAs when solving ZDT3

3 结语

文中提出了改进的SPEA2算法——CM-SPEA2。该算法重新定义了杂乱度计算方法，并将其用于个体适应度计算，再将改进的适应度值应用于父代选择和环境选择策略。环境选择策略通过改进适应度函数被间接调整：当解集分布不均匀时，通过WPGMA 聚类的方法将种群归类为不同的簇，再计算各个簇内个体的杂乱度；找到杂乱度最低的个体，按照就近原则接纳该个体附近的个体，最终填满整个归档集；当解集分布较为均匀时，则直接在种群内找到杂乱度最低的个体，再按照聚集密度从小到大的次序依次纳入归档集。通过Matlab PlatEMO 平台，选择3 种测试MOP（IMOP、ZDT 及VNT）与4 种MOEA（SPEA2、NSGAⅡ、CAMOEA 及RVEAa）进行比较，发现CM-SPEA2在IGD和Spacing指标上均不同程度地占一定优势，且在多数MOP上比SPEA2性能更有提升。后续研究中，将重点探讨动态多目标优化、多目标优化的鲁棒性［20］等，以期进一步改进算法，并将其应用于机器人路径规划。