0 引言
近些年，用多目标优化技术来解决约束问题是个热点[1-2].首先要找到全局内的可行解集，其次才能考虑解集的分布性和收敛性，其中设计有效的约束处理技术是关键所在.因此，将进化算法引入到约束多目标优化问题中是必须的.除了要解决多目标优化过程中面临的提高进化算法搜索能力、避免陷入局部最优、合理设置参数等多个问题外，还必须保证算法获得可靠的求解性能和全局的优化效果.因此，约束多目标优化算法的研究内容较多，不仅需要重点研究约束处理技术，还必须研究与之相适应的进化策略、多样性维护策略以及精英选择策略等多项关键技术.
孟红云等[3]提出一种用于约束多目标优化问题的双群体差分进化算法.该算法同时使用两个群体，一个群体用于保存搜索过程中找到的可行解； 另一个用于记录在搜索过程中得到的部分具有某些优良特性的不可行解, 避免了构造罚函数和直接删除不可行解.但是，该方法在更新可行解集时会存在丢失边界解和不利于保持种群的整体多样性的缺陷；而在更新不可行解集时会出现目标函数较差的个体，从而影响种群的收敛.Deb提出了一种可行性法则[4]来比较个体, 但Deb准则认为可行解优于不可行解，没有利用较优不可行解所携带的信息,不利于保持种群的多样性.最近，一些学者基于精英保留机制提出了双档案集方法，将进化过程产生的可行解和不可行解个体分别存入可行解档案和不可行档案.例如，文献[5]提出了基于遗传算法的双档案集机制；文献[6-9]提出了基于差分算法的双档案方法.然而这些方法仅粗略地以约束违反度值作为依据，优先选择靠近可行边界的不可行解，没有区分出最优不可行解，将直接影响算法的寻优能力.
基于此，笔者提出一种两阶段三存档集约束优化算法，主要创新点包括：①提出新的三存档集处理约束条件，避免了处于Pareto边缘的解与函数值较差的解采用同样的进化策略；②对非支配解存档集和支配解存档集采用侧重点不同的搜索方式.非支配解存档集采用混合策略进化，既进行局部搜索又进行全局搜索，保证在Pareto前沿附近的解快速收敛，又不会陷入局部最优；支配解存档集采用自适应变异，将变异率和个体信息联系起来，引导种群进化.
1 约束优化问题及相关定义
不失一般性，一个约束的多目标优化问题可定义如下[1]：
(1)
式中:x=[x1, x2,…,xn]∈S是决策变量，S表示n维搜索空间，它满足如下边界约束条件：li≤xi≤ui，1≤i≤n;F(x)是目标函数；q表示不等式约束条件的个数；m表示总约束条件的个数；gj(x)为第j个不等式约束条件;hj(x)为第j-q个等式约束条件；ui和li分别为分量xi的上界和下界.
满足所有约束条件的解是可行解，而所有的可行解组成可行域Ω∈S；不满足任意一个约束条件的解被称为不可行解，所有的不可行解组成不可行域.如果gj(x)=0(j=1,2,…,q)在一个可行解处成立，则称该不等式约束条件在该可行解处是活跃的.
个体x在第j个约束条件上的约束违反度表示为：
(2)
式中:δ为等式约束的容忍参数，该参数可以定义为一个特定的值(一般设为0.000 1)或采用只适应变化的方式定义.所以个体x总的约束违反度表示为：
2 两阶段三存档集约束优化算法(TSDA)
通过对三存档集进行搜索机制不同的进化方式，来引导种群向真实Pareto前沿进行进化.利用非支配解存档集双重寻优，一方面使得靠近Pareto前沿的解加速向前沿靠近；另一方面使存档集避免陷入局部收敛.支配解存档集利用自适应的变异概率，将目标函数值与种群个体的信息联系起来，提高了算法的搜索效率.最后通过非支配可行解存档集保留种群当中的精英个体，提高了算法约束优化性能.
2.1 最优不可行解
在非支配解存档集中让最优不可行解进入存档集参与进化，不仅能加快算法的收敛速度，同时也增强了可行域边界附近的寻优能力.判定一个不可行解是否位于可行域附近，具体如下：若一个不可行解的约束违反度小于等于当代的ε(t)值(可行阈值)，则其为靠近可行域边界的不可行解.
定义1 ε值(可行阈值)：
(3)
式中：t是当前代数；T为进化总代数；ε(0)是种群中个体根据约束违反度升序排序后第0.05·N个个体，N为种群的规模.
定义2 最优不可行解：靠近可行域边界，且在可行域中找不到能支配该不可行解的可行解，则称该不可行解为最优不可行解，如图1所示.
图1 最优不可行解
Fig.1 The optimal infeasible solution
图1中点i及k皆为靠近可行域边界(图中的阴影部分为可行域)的不可行解，不可行解i虽然靠近可行域边界，但其远离Pareto前沿，对种群进化没有引导作用，因此，在选择最优不可行解时，需要跳过这种不可行解.而不可行解k既满足靠近可行域边界又满足在可行域中找不到能支配k的可行解，则k为最优不可行解，让其参与进化.
2.2 三存档集模型
笔者基于NSGA-Ⅱ的基础上加入了两阶段三存档模型，该算法将进化过程划分为两个阶段.首先，将种群利用快速非支配排序法进行分层.根据Pareto等级、约束违反度及ε(t)值划分为非支配解(包括可行解及不可行解)、支配解及非支配可行解，分别存入3个存档集：非支配解存档集、支配解存档集以及非支配可行解存档集.然后,种群进行两阶段进化操作.第一阶段：非支配解存档集及支配解存档集采取不同的策略分别进行进化；第二阶段：判断非支配可行解存档集中解的个数是否超过种群数量，如果超过则将其进行快速非支配排序，计算出非支配可行解存档集中个体的rank值及拥挤度距离值，选出种群规模大小的优秀个体.
2.3 非支配解存档集进化
定义3 非支配解存档集(NDA)：根据2.1节所计算出当代的ε(t)值区分出种群中约束违反度小于等于ε(t)值的解，其中既包括可行解，也包括不可行解.若该解rank值等级为1，则为非支配解，非支配解的集合称为非支配解存档集.
非支配解存档集(NDA)既进行局部搜索又进行全局搜索.NDA存档集中的解虽然皆为Pareto等级为1的解，但这些解有靠近真实Pareto前沿的、也有远离Pareto前沿的.若只进行局部搜索，很可能陷入局部最优当中，因此，笔者采用非支配解存档集局部搜索和全局搜索同时进行.因为非支配解存档集中的不可行解为靠近可行域且在可行域当中没有能支配它的可行解，所以非支配不可行解只进行局部搜索可以更好地引导种群向Pareto前沿进化.非支配解存档集当中的可行解分为两种情况：一种情况为该可行解位于Pareto前沿附近，采用局部搜索策略可以使其尽快收敛；另外一种情况为该可行解远离Pareto前沿，采用全局搜索策略可以扩大搜索的空间，避免陷入局部最优.鉴于以上两种情况，非支配解存档集当中的可行解为不确定因素，故笔者采用两种不同交叉变异概率对其进行进化.首先将可行解复制一份放入集合CNDFA中，然后非支配解存档集整体进行局部搜索，而复制的可行解集CNDFA进行全局搜索.非支配解存档集进化过程中，交叉算子取0.5，变异算子取0.05，在进化过程中生成一个随机数，若随机数u≤0.5，u为[0,1]上均匀分布的随机数，则进行SBX交叉操作；若随机数u≥0.95，则进行多项式变异操作；0.5<若随机数u<0.95，则按照约束违反度值从小到大，从非支配解存档集合中依次选择两个个体进入下一代.CNDFA集合仍采用原NSGA-Ⅱ的SBX交叉和多项式变异操作.
2.4 支配解存档集进化
定义4 支配解存档集(DA)：若该解不为非支配解，则其为支配解，支配解的集合称为支配解存档集.
支配解存档集进行全局搜索，具体操作如下.
(1)交叉操作：为了跳出局部最优以及加强空间搜索能力，笔者将正态分布引入到SBX算子中.由于正态分布的引入，算法可搜索到的空间更为广阔，因此更容易跳出局部最优，从而使Pareto前沿更具延展性，均匀地分布在Pareto域上，保证种群的多样性.NDX算子的产生方式如下：
如果u≤0.5，则
(4)
如果u>0.5，则
(5)
式中:c1/2,k为子代染色体上对应的第k个变量；p1,k、p2,k分别为父代两个染色体上对应的第k个变量；|N(0,1)|为正态分布随机变量.
(2)变异操作：经典NSGA-Ⅱ算法中多项式变异算子与种群个体的信息并没有联系，无论种群进化到何种程度，变异概率都为固定值，种群没有特定的方向进行进化.因此，将种群个体当前信息与变异算子结合起来时，才能使种群向着真实的Pareto前沿进化.自适应变异算子[10]产生方式如下：
pm=
(6)
式中：fmax、fmin分别为当前种群中目标函数的最大值、最小值；favg为当前种群所有个体平均目标函数值.
在进化的初期，个体的最大目标函数值fmax与最小目标函数值fmin差异较大，目标函数值较大和较小的个体数目大致相同，此时favg近似等于fmax与fmin的平均值.由式(6)知，(fmax-favg)/(fmax-fmin)的值大致为0.5，此时求得的变异概率较大，有助于算法初期进行全局搜索，寻找最优解集.算法进化到后期，大部分种群个体求得的目标函数值大致相同，此时favg是一个略大于fmin的值.当所有个体都进化到最优解时，达到一种极限条件，即fmin=favg，这时公式中(fmax-favg)/(fmax-fmin)的值接近1，交叉变异概率相对调整为较小的数，增进局部寻优的能力，与种群个体寻优的方向相一致，有利于最优解集的搜索.
2.5 非支配可行解存档集
笔者使用非支配可行解存档集来保存进化过程中产生的所有可行非劣解.首先，在种群每次迭代过程中，将非支配存档集中约束违反度为0的解提取出来，加入到非支配可行解集当中；然后，通过快速非支配排序，选出的N个优秀个体进行输出，此时输出的解全部为最优的可行解.
2.6 算法流程
两阶段三存档集约束优化算法流程：
Step 1开始gen=1.
Step 2初始化种群.从决策空间S随机产生N个个体的初代种群P(t)={x1,t,…,xN,t}，并根据目标函数计算每个个体的目标函数值.
Step 3非支配排序.对Step 2产生的种群进行非支配排序，按照非支配等级由小到大排序，并计算每个个体的约束违反度.
Step 4存档集划分.将种群中每个个体的约束违反度值与可行阈值ε(t)值相比，再根据rank值大小，分类到非支配解存档集(NDA)、支配解存档集(DA)及非支配可行解存档集(NDFA)当中.
Step 5基因操作.对非支配解存档集和支配解存档集采取不同的进化策略.
Step 5.1非支配解存档集采用混合搜索策略.根据2.3节中的内容分别对NDA和CNDFA两个集合进行局部搜索和全局搜索，合并进化后的两个集合，记作U(t).
Step 5.2支配解存档集采用全局搜索.利用正态分布算子和自适应变异算子对DA进行交叉变异操作，从而生成子种群Q(t).
Step 6合并子父代种群.将U(t)，Q(t)及父代种群P(t)进行合并，记作Y(t).
Step 7选择操作.将Y(t)进行快速非支配排序，根据rank值及拥挤度距离值大小从Y(t)中选择出N个个体作为下一代种群P(t+1).
Step 8判断当前代数是否大于设定代数，如满足则跳转至Step 9，进行第二阶段的操作；否则，gen=gen+1，并跳转至Step 3.
Step 9对非支配可行解存档集(NDFA)进行快速非支配排序，根据rank值及拥挤度距离选择出N个个体，输出最终的Pareto最优解，算法终止.
2.7 计算复杂度分析
假设非支配解存档集规模为N1，支配解存档集规模为N2，非支配可行解存档集规模为N3.则计算非支配解存档集、支配解存档集和非支配可行解存档集目标函数和约束违反度的时间复杂度为O(m(N1+N2+N3)).其中，m为目标函数个数；N为种群规模；非支配解存档集变异和交叉操作的时间复杂度为O(mN1)；支配解存档集变异和交叉操作的时间复杂度为O(mN2)；选择操作的时间复杂度为O(m(N1+N2+N3)2).综上，TSDA算法迭代一次的最坏时间复杂度为O(m(N1+N2+N3))+O(mN1)+O(mN2)+O(m(N1+N2+N3)2)，所以本文算法时间复杂度为O(mN2).
3 TSDA算法实验结果与分析
3.1 测试函数与参数设置
所有实验在硬件配置为Intel Pentium、G2030 CPU、4G内存、主频3.0 GHz、Win10系统的计算机上进行，程序采用Matlab R 2010 编写.
为了验证TSDA算法的性能，选择了SRN、TNK、OSY作为测试函数，其函数表达式见表1.并与运用可行性规则的NSGA-Ⅱ算法(FRNSGA-Ⅱ)进行比较.
TSDA算法参数设置：随机生成的种群规模为N=200，进化代数为2 000代，交叉算子为NDX(正态分布)交叉算子，采用自适应变异算子，交叉变异概率详见2.4节.
FRNSGA-Ⅱ算法参数设置：随机生成的种群规模为N=200，进化代数为2 000代，交叉算子为SBX算子，交叉概率为0.95，变异算子为多项式变异，变异概率为0.05.
测试函数特点：SRN函数的Pareto前沿连续；TNK函数为非线性约束函数，且其Pareto前沿不连续；OSY函数共有6个不等式约束条件，既有线性约束又有非线性约束;Pareto前沿由三段折线组成.约束测试函数的真实Pareto前沿源自参考文献[11].
3.2 评价指标
3.2.1 GD指标
世代距离(generational distance，GD)[12]，用于评价所求得的近似Pareto前沿Pknow相对于真实Pareto前沿Ptrue的逼近程度，该指标的定义如下：
(7)
式中：NPF为Pknow中个体的数量；di为Pknow中第i个个体的目标向量到Ptrue中最近个体的欧式距离;GD越小，表明Pknow越逼近Ptrue，收敛性越好.
表1 测试函数表达式
Tab.1 Description of test functions
测试函数目标函数min F(x)=(f1(x),f2(x))变量约束决策变量范围SRNmin f1(x)=2+(x1-2)2+(x2-1)2min f2(x)=9x1-(x2-1)2 x21+x22≤225x1-3x2+10≤0 -20≤xi≤20i=1,2TNKmin f1(x)=x2min f2(x)=x1 -x21-x22+1+0.1cos(16arctan(x1/x2))≤0(x1-0.5)2+(x2-0.5)2≤0.5 0≤xi≤πi=1,2OSYmin f1(x)=4x21+4x22min f2(x)=(x1-5)2+(x2-5)2 x1+x2-2≥06-x1-x2≥02+x1-x2≥02-x1+3x2≥04-(x3-3)2-x4≥0(x5-3)2+x6-4≥0 0≤xi≤10i=1,2,61≤x3,x5≤50≤x4≤6
3.2.2 SP指标
Schott在1995年提出了Spacing指标[13],用来计算所求得的解集在目标空间上的分布均匀性，其计算公式如下所示：
(8)
式中:
表示离第i个个体最小的目标距离和；d′为所有di的均值；NPF为所求得解集的个体数量；S越小表示所求得的解集分布的越均匀，当S为0时，表示所求解集中的所有解都是等距离分布的.
表2为2种算法对上述3种测试函数的SP、GD的统计结果.可以看出，本文算法TSDA在测试函数TNK、OSY上的GD值明显小于对比算法——FRNSGA-Ⅱ算法，说明笔者所提算法采用的策略避免了种群陷入局部最优，对种群收敛到Pareto前沿起着重要的作用.在SP值方面，算法也有着明显的优势，这说明所提算法的分布性更好，保证了种群的多样性.
3.3 仿真结果
图2、图4、图6为本文算法TSDA所求得的Pareto前沿，图3、图5、图7为基于可行性规则的NSGA-Ⅱ算法所求得的Pareto前沿.通过以上的对比可以看出:对于SRN函数，两种算法虽然都找到了真实的Pareto前沿，但FRNSGA-Ⅱ算法的分布性和收敛性较差；对于TNK测试函数，本文算法TSDA有更好的收敛性，并保持良好的分布性；对于OSY测试函数，FRNSGA-Ⅱ算法所得到的解集明显偏离真实前沿.由此可以得出，笔者提出的TSDA算法在3个约束测试问题上，其最优Pareto前沿在逼近性和分布性上都明显优于FRNSGA-Ⅱ.
4 结论
针对约束优化算法处理位于Pareto边缘的解与函数值较差的解采用相同的进化策略使得寻优结果较差的问题，提出两阶段三存档集约束优化算法.利用两个存档集不同的寻优方式，很好地兼顾了探索开发和收敛的平衡.并且通过与其他算法在3种测试函数上的对比表明，本文算法在处理约束多目标优化问题时，分布性及收敛性均具有一定优势.
表2 2种进化算法求解标准测试函数所得的GD、SP值
Tab.2 GD and SP values obtained by two evolutionary algorithms for solving standard test functions
测试函数GDSPTSDAFRNSGA-ⅡTSDAFRNSGA-Ⅱ平均值方差平均值方差平均值方差平均值方差SRN2.8346E-042.5521E-030.022 30.528 60.093 20.001 41.176 20.972 3TNK1.137E-048.1149E-037.499 57.042 70.134 87.865E-042.435 73.071 5OSY0.612 50.304 02.178 50.718 63.246 81.559 78.523 62.461 8
图2 TSDA算法SRN函数仿真结果
Fig.2 Simulation result of TSDA algorithm on SRN function
图3 FRNSGA-Ⅱ算法SRN函数仿真结果
Fig.3 Simulation result of FRNSGA-Ⅱ on SRN function
图4 TSDA算法TNK函数仿真结果
Fig.4 Simulation result of TSDA algorithm on TNK function
图5 FRNSGA-Ⅱ算法TNK函数仿真结果
Fig.5 Simulation result of FRNSGA-Ⅱ on TNK function
图6 TSDA算法OSY函数仿真结果
Fig.6 Simulation result of TSDA algorithm on OSY function
图7 FRNSGA-Ⅱ算法OSY函数仿真结果
Fig.7 Simulation result of FRNSGA-Ⅱ on OSY function
然而，必须承认，本文算法TSDA在计算量上要高于FRNSGA-Ⅱ算法.接下来的研究将致力于如何降低算法的时间复杂度及本文算法的实际应用.
[1] 李智勇,黄滔,陈少淼,等.约束优化进化算法综述[J]. 软件学报,2017,28(6):1529-1546.
[2] 梁静,刘睿,瞿博阳,等.进化算法在大规模优化问题中的应用综述[J].郑州大学学报(工学版),2018,39(3):15-21.
[3] 孟红云,张小华,刘三阳.用于约束多目标优化问题的双群体差分进化算法[J]. 计算机学报, 2008, 31(2):228-235.
[4] DEB K. An efficient constraint handling method for genetic algorithms[J]. Computer methods in applied mechanics & engineering, 2000, 186(2):311-338.
[5] HAMMACHE A, BENALI M, AUBÉ F. Multi-objective self-adaptive algorithm for highly constrained problems: novel method and applications[J]. Applied energy, 2010, 87(8):2467-2478.
[6] 俞国燕,李鹏,何真,等.一种用于多目标约束优化的改进进化算法[J].计算机集成制造系统,2009, 15(6):1172-1178.
[7] 毕晓君,张磊,肖婧.基于双种群的约束多目标优化算法[J].计算机研究与发展,2015,52(12):2813-2823.
[8] 暴励,曾建潮.一种双种群差分蜂群算法[J].控制理论与应用,2011,28(2):266-272.
[9] 毕晓君,张磊.基于混合策略的双种群约束优化算法[J].控制与决策,2015,30(4):715-720.
[10] 肖易寒,李明逵,陈立伟.基于改进NSGA-Ⅱ算法的多光谱测温数据处理[J].应用科技, 2017, 44(1):33-39.
[11] DEB K, PRATAP A, MEYARIVAN T. Constrained test problems for multi-objective evolutionary optimization[C]//Proceedings of First International Conference on Evolutionary Multi-Criterion Optimization, EMO 2001. Zurich, Switzerland:EMO, 2001:284-298.
[12] 韩红艳. 基于Pareto支配的高维多目标进化算法研究[D]. 大连：大连理工大学计算机科学与技术研究, 2016.
[13] JIANG S, YANG S. A steady-state and generational evolutionary algorithm for dynamic multi-objective optimization [J]. IEEE transactions on evolutionary computation, 2017, 21(1): 65-82.