多目标优化算法（四）NSGA3(NSGAIII)论文复现以及matlab和python的代码

前言：最近太忙，这个系列已经很久没有更新了，本次就更新一个Deb大神的NSGA2的“升级版”算法NSGA3。因为multi-objective optimization已经被做烂了，现在学者们都在做many-objective optimization，也就是5个以上的目标函数（悄悄说一句，我觉得这个也要被做烂了）。此次我是用python复现的，这篇文章也主要以python代码讲解为主。在编写代码过程中，一些小技巧借鉴了platEMO，这个平台在我之前的博客里已经介绍过了，里面包含了近乎所有知名的多目标优化算法，想学习的小伙伴可以看我之前的博客介绍：https://blog.csdn.net/qq_40434430/article/details/88366639
为了哪些只想使用NSGA3的同学，我将NSGA3的matlab代码从平台中扣了出来，想要的在我的主页里下载吧，链接如下：（matlab我使用的是2017版的）
https://download.csdn.net/download/qq_40434430/11079440
python代码使用的python3.6，其它版本的也可以，毕竟没有用到什么复杂的库，就只是用到了numpy等常用的库。完整的python代码大家也可以去我CSDN主页下载哈，连接如下：
https://download.csdn.net/download/qq_40434430/11079551

摘要：此次博客主要记录了Kalyanmoy Deb和Himanshu Jain《An Evolutionary Many-Objective Optimization Algorithm Using Reference-point Based Non-dominated Sorting Approach, Part I: Solving Problems with Box Constraints》的论文学习心得[1]，此篇文章提出了处理多个优化目标的进化优化算法NSGAIII。其主要思路是在NSGAII的基础上，引入参考点机制，对于那些非支配并且接近参考点的种群个体进行保留。此次复现处理的优化问题是具有3到15个目标的DTLZ系列[2]，仿真结果反应了NSGAIII良好的搜索帕累托最优解集的能力。

关键字：Many-objective optimization，高维问题，NSGAIII，非支配排序，参考点

I.问题简介

在很多真实的应用中，优化目标往往不只有一两个，而是有4个以上。针对这样的高维目标空间，现有的EMO显得有些力不从心，其主要问题为以下六点：

随着优化目标数量的增加，非支配解在种群中的比例也在增加，因而会导致搜索过程缓慢；
对于高维目标空间，维持多样性的指标计算复杂度过高，解的邻近元素寻找困难；
对于高维目标空间，重组算子的搜索能力过于低效了；
因为目标函数较多，pareto前沿难以表示，决策者无法选择自己需要的解；
性能指标的计算代价过大，算法结果不易评价；
针对高维的目标空间，如何可视化结果也是一个难题。

对于前三个问题，我们可以改造EMO缓解，但是后三个问题目前还没有解决办法。除了以上问题，实际问题中由于目标解集集中在pareto front的一个小区域，因此如何寻找也是算法应用到实际中的一个阻碍。但有研究发现目标函数往往会退化成低维的优化，一般会低2~3维，因此冗余目标的识别也是一个难点。以下提出两种解决1，2，3问题的思路：

使用特殊的支配关系：可以引入自适应离散化pareto front从而解决问题1、2，使用具有较大指数分布的SBX解决问题3。
使用预定义的目标搜索：这种方式又可以分为两种：第一种是预定义一组跨越整个pareto front的搜索方向，这种方法的代表是MOEA/D；另一种是预定义多个参考点，比如算法NSGAIII。

II.NSGAIII算法要点和总结

总体上来说，NSGAIII和NSGAII具有类似的框架，二者区别主要在于选择机制的改变，NSGAII主要靠拥挤度进行排序，其在高维目标空间显然作用不太明显，而NSGAIII对拥挤度排序进行了大刀阔斧的改编，通过引入广泛分布参考点来维持种群的多样性。

以下总结NSGAIII的第t代的步骤。PtP_tPt是第t代的父代，其大小为N，其生成的子代为QtQ_tQt，其大小也为N。第一步将子代和父代结合Rt=Pt∪QtR_t=P_t\cup Q_tRt=Pt∪Qt(大小为2N)并从中选出N个个体。为了实现这个选择过程，首先将RtR_tRt通过非支配排序分为多个非支配层(F1,F2,...F_1,F_2,...F1,F2,...)。然后从F1F_1F1开始构造一个新的种群StS_tSt，直到其大小为N或者第一次超过N。称最后一层为第lll层。第l+1l+1l+1层以上的解将被淘汰出局，在大多数情况下，最后一层被接受层(lll层)仅有部分被接受。在这种情况下，就要用多样性衡量lll层里的解进而进行选择。NSGAII里的这部分使用了拥挤度排序，NSGAIII中我们用以下5步替代。下面先给出这个NSGAIII的第t代的算法步骤如下：

这里需要注意：算法输入里的H是结构化参考点的数目，我在后文中定义为P。
主程序的python代码如下：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # 空间三维画图
from utils import uniformpoint,funfun,cal,GO,envselect,IGD
import copy
import random#参数设置
N_GENERATIONS = 400                                 # 迭代次数
POP_SIZE = 100                                      # 种群大小
name = 'DTLZ1'                                      # 测试函数选择，目前可供选择DTLZ1,DTLZ2,DTLZ3
M = 3                                               # 目标个数
t1 = 20                                             # 交叉参数t1
t2 = 20                                             # 变异参数t2
pc = 1                                              # 交叉概率
pm = 1                                              # 变异概率
#画图部分
if(M<=3):fig = plt.figure()ax = Axes3D(fig)###################################################################################################################################################################
#产生一致性的参考点和随机初始化种群
Z,N = uniformpoint(POP_SIZE,M)#生成一致性的参考解
pop,popfun,PF,D = funfun(M,N,name)#生成初始种群及其适应度值，真实的PF,自变量个数
popfun = cal(pop,name,M,D)#计算适应度函数值
Zmin = np.array(np.min(popfun,0)).reshape(1,M)#求理想点
#ax.scatter(Z[:,0],Z[:,1],Z[:,2],c='r')
#ax.scatter(PF[:,0],PF[:,1],PF[:,2],c='b')#迭代过程
for i in range(N_GENERATIONS):print("第{name}次迭代".format(name=i)) matingpool=random.sample(range(N),N)   off = GO(pop[matingpool,:],t1,t2,pc,pm)#遗传算子,模拟二进制交叉和多项式变异offfun = cal(off,name,M,D)#计算适应度函数mixpop = copy.deepcopy(np.vstack((pop, off)))Zmin = np.array(np.min(np.vstack((Zmin,offfun)),0)).reshape(1,M)#更新理想点pop = envselect(mixpop,N,Z,Zmin,name,M,D)popfun = cal(pop,name,M,D)if(M<=3):ax.cla()type1 = ax.scatter(popfun[:,0],popfun[:,1],popfun[:,2],c='g')plt.pause(0.00001)# 绘制PF
if(M<=3):type2 = ax.scatter(PF[:,0],PF[:,1],PF[:,2],c='r',marker = 'x',s=200)plt.legend((type1, type2), (u'Non-dominated solution', u'PF'))
else:fig1 = plt.figure()plt.xlim([0,M])for i in range(pop.shape[0]):plt.plot(np.array(pop[i,:]))
plt.show()#IGD
score = IGD(popfun,PF)

1.将种群按照非支配层进行划分

将非支配层等级1到lll的种群成员依次放入StS_tSt中，如果∣St∣=N|S_t|=N∣St∣=N，则无需进行下面的操作，直接Pt+1=StP_{t+1}=S_tPt+1=St。如果∣St∣>N|S_t|>N∣St∣>N，那么下一代的一部分解为Pt+1=∪i=1l−1FiP_{t+1}=\cup _{i=1}^{l-1}F_iPt+1=∪i=1l−1Fi，剩余部分(K=N−∣Pt+1∣K=N-|P_{t+1}|K=N−∣Pt+1∣)从FlF_lFl中选择。这个选择过程见步骤2到5。
非支配排序的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
def NDsort(mixpop,N,M):nsort = N#排序个数N,M = mixpop.shape[0],mixpop.shape[1]Loc1=np.lexsort(mixpop[:,::-1].T)#loc1为新矩阵元素在旧矩阵中的位置，从第一列依次进行排序mixpop2=mixpop[Loc1]Loc2=Loc1.argsort()#loc2为旧矩阵元素在新矩阵中的位置frontno=np.ones(N)*(np.inf)#初始化所有等级为np.inf#frontno[0]=1#第一个元素一定是非支配的maxfno=0#最高等级初始化为0while (np.sum(frontno < np.inf) < min(nsort,N)):#被赋予等级的个体数目不超过要排序的个体数目maxfno=maxfno+1for i in range(N):if (frontno[i] == np.inf):dominated = 0for j in range(i):if (frontno[j] == maxfno):m=0flag=0while (m<M and mixpop2[i,m]>=mixpop2[j,m]):if(mixpop2[i,m]==mixpop2[j,m]):#相同的个体不构成支配关系flag=flag+1m=m+1 if (m>=M and flag < M):dominated = 1breakif dominated == 0:frontno[i] = maxfnofrontno=frontno[Loc2]return frontno,maxfno

2.超平面上参考点的确定

NSGAIII使用一组预定义的参考点以确保解的多样性，这一组参考点可以结构化的方式定义，也可以用户根据自己的参考点。以下介绍一种产生结构化参考点的方法叫Das and Dennis’s method，此方法来源于田野老师对产生参考点方法的综述论文[3]。其参考点在一个(M-1)维的超平面上，M是目标空间的维度，即优化目标的个数。如果我们将每个目标划分为HHH份，那么其参考点的数量为
P=(m+H−1H)P=(\begin{matrix} m+H-1\\ H \end{matrix}) P=(m+H−1H)
例如对于一个H=4H=4H=4的三目标问题，其参考点构成了一个三角形，根据公式(1)可知其产生15个参考点，见图1所示。

图1 对于一个H=4的三目标的问题产生的15个参考点

那么它们的坐标该如何计算呢？以下通过田野老师的论文讲解一下如何产生这些参考点。假设由Das and Dennis’s method产生的参考点为s=(s1,s2,...,sM)s=(s_1,s_2,...,s_M)s=(s1,s2,...,sM)，这里
sj∈{0/H,1/H,...,H/H},∑j=1Msj=1s_j\in \left\{ 0/H,1/H,...,H/H \right\},\sum_{j=1}^{M}s_j=1 sj∈{0/H,1/H,...,H/H},j=1∑Msj=1
HHH为每个目标划分的数目。此处举例假设M=3,H=5M=3,H=5M=3,H=5（注意不是上图的例子）。其计算参考点的坐标过程如图2所示。

图2 一个用Das and Dennis’s method产生参考点坐标的过程

这里，我们找到a和b所有使得a,b∈{0,0.2,...,1},a≤ba,b\in \left\{ 0,0.2,...,1 \right\} ,a \leq ba,b∈{0,0.2,...,1},a≤b的组合(0.2是因为1划分为5份)，然后令s1=a−0,s2=b−a,s3=1−bs_1=a-0,s_2=b-a,s_3=1-bs1=a−0,s2=b−a,s3=1−b。这里产生的a,b其实就是为了使用排列组合里的插空法，从而产生所有满足(2)条件的坐标。这个算法流程如下：

这个产生方法其实有个bug，比如当目标函数数目比较大的时候(例如等于10)，那么假设每个函数划分为10段，则会产生92378个参考点，这显然太大了，因此在复现NSGA-III里，我使用以下的参考点产生方法，其叫做Deb and Jain’s Method，此方法也可以在参考文献[3]里找到，其主要想法就是产生内层和外层两个参考点的超平面，这样可以减少参考点的产生，并保证参考点的广泛分布，其流程如下：

例如，一个三目标问题，每个目标函数划分为13份，利用此方法其只产生了210个参考点就保证了广泛分布性。如图3所示。

图3 一个用Deb and Jain’s Method产生的参考点(M=3,H=13)

有人不禁会问这个内外层的P如何确定。在NSGAIII中，我们设定希望产生N个参考点，取最接近但不大于N的个数的H1H_1H1和H2H_2H2参数，来确定内层和外层的参考点个数。其实这种方法里，我们虽希望产生N个参考点，但一般而言，产生的参考点数目PPP要比N少，但是数目接近。产生这样的参考点的主要目的还是为了借助它找到对应的近似帕累托最优解。

参考点生成的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
def uniformpoint(N,M):H1=1while (comb(H1+M-1,M-1)<=N):H1=H1+1H1=H1-1W=np.array(list(combinations(range(H1+M-1),M-1)))-np.tile(np.array(list(range(M-1))),(int(comb(H1+M-1,M-1)),1))W=(np.hstack((W,H1+np.zeros((W.shape[0],1))))-np.hstack((np.zeros((W.shape[0],1)),W)))/H1if H1<M:H2=0while(comb(H1+M-1,M-1)+comb(H2+M-1,M-1) <= N):H2=H2+1H2=H2-1if H2>0:W2=np.array(list(combinations(range(H2+M-1),M-1)))-np.tile(np.array(list(range(M-1))),(int(comb(H2+M-1,M-1)),1))W2=(np.hstack((W2,H2+np.zeros((W2.shape[0],1))))-np.hstack((np.zeros((W2.shape[0],1)),W2)))/H2W2=W2/2+1/(2*M)W=np.vstack((W,W2))#按列合并W[W<1e-6]=1e-6N=W.shape[0]return W,N

3.种群个体的自适应归一化

第一步定义种群StS_tSt的理想点为∪τ∣=0tSτ\cup_{\tau|=0}^tS_{\tau}∪τ∣=0tSτ的最小值(zimin,i=1,2,...,Mz_i^{min},i=1,2,...,Mzimin,i=1,2,...,M)，那么就可以构造出一个理想点zˉ=(z1min,z2min,...,zMmin)\bar{z}=(z^{min}_1,z^{min}_2,...,z^{min}_M)zˉ=(z1min,z2min,...,zMmin)。这个理想点就很厉害了，我们通过这个理想点可以将目标函数转换，其公式如下fi′(x)=fi(x)−ziminf^{'}_i(x)=f_i(x)-z^{min}_ifi′(x)=fi(x)−zimin。第二步就是求每个坐标轴对应的额外点(extreme point)，通过这个额外点我们可以构造出一个超平面。第三步利用超平面与坐标轴的节距，我们就可以自适应归一化目标函数值啦。

产生这个额外点的公式如下：
ASF(x,w)=maxi=1Mfi′(x)/wi,x∈StASF(x,w)=max_{i=1}^{M}f_i^{'}(x)/w_i,x\in S_t ASF(x,w)=maxi=1Mfi′(x)/wi,x∈St

zi,max=s:argmins∈StASF(s,wi),wi=(τ,...,τ),τ=10−6,wji=1z^{i,max}=s:argmin_{s\in S_t}ASF(s,w^i),w^i=(\tau,...,\tau),\tau=10^{-6},w^i_j=1 zi,max=s:argmins∈StASF(s,wi),wi=(τ,...,τ),τ=10−6,wji=1

对于第i个转换目标fi′f_i^{'}fi′，产生一个额外目标向量zi,maxz^{i,max}zi,max。这M个额外目标向量构成一个M维的线性超平面，如图4所示。进而求出截距ai,i=1,...,Ma_i,i=1,...,Mai,i=1,...,M，那么目标函数就可以被归一化为：
fin(x)=fi′(x)/(ai−zimin)=(fi′(x)−zimin)/(ai−zimin),i=1,2...,Mf_i^n(x)=f_i^{'}(x)/(a_i-z^{min}_i)=(f_i^{'}(x)-z^{min}_i)/(a_i-z^{min}_i),i=1,2...,M fin(x)=fi′(x)/(ai−zimin)=(fi′(x)−zimin)/(ai−zimin),i=1,2...,M
这里归一化平面与坐标轴交点的函数值fi′=1f_i^{'}=1fi′=1，此归一化超平面上的点满足∑i=1Mfin=1\sum_{i=1}^{M}f_i^n=1∑i=1Mfin=1。

图4 给出了一个三目标问题的截距计算过程，并从极值点生成超平面(图是未归一化的平面)

以下给出这个过程的伪代码：

在这一步中，主要有三点不好理解及需要说明的地方，以下依次解释：

第一，为什么需要自适应归一化呢？由于后面我们将每一个解和参考点相互联系从而维持多样性，参考点又是均匀分布在目标空间中的，而每个解的各个目标函数值尺度不一样导致解的偏向性不一样，比如目标函数f1的范围为[0,1]，f2的范围为[0,100]，那么在联系解和参考点时，f1和f2起的作用就不“公平”了。

第二，如何理解这个额外点呢？这个额外点其实是指在一个目标值上很大，另外几个目标值上很小的点。以下举一个例子，比如现在有三个个体，他们的目标函数值为(1,2,6),(2,5,2),(4,9,3)。公式(3)(4)其实做的事情就是每一个除以一个w=(1,10−6,10−6)w=(1,10^{-6},10^{-6})w=(1,10−6,10−6)，即固定第一维的坐标。可以理解为沿x轴看去，这三个点y和z值中最大的那个谁更小，也就是谁更接近x平面。这个例子中，(1,2,6)变成了(1,2∗106,6∗106)(1,2*10^6,6*10^6)(1,2∗106,6∗106)，(2,5,2)变成了(2,5∗106,2∗106)(2,5*10^6,2*10^6)(2,5∗106,2∗106)，(4,9,3)变成了(4,9∗106,3∗106)(4,9*10^6,3*10^6)(4,9∗106,3∗106)。三个点对应在y,z轴上的大的那个分别为6∗106,5∗106,9∗1066*10^6,5*10^6,9*10^66∗106,5∗106,9∗106。最小的就是点(2,5,2)。这个作为额外点就很合适啦。其实仔细研究你就会发现，这个产生额外点的算法是有一些问题的，比如图5：

图5 一种新的产生额外点的方法说明

按照此方法产生的额外点是图中的worst point对应的两个边界点，而一个好的额外的点是Nadir point对应的两个边界点，因为额外点是为了自适应归一化，显然worst point对应f1和f2的“缩放”程度与Nadir point不同。这个产生Nadir point的算法见文献[4]。

第三，如果截距不存在怎么办呢？袁源的基于分解的多目标进化算法及其应用[5]里给出里解决方法：如果矩阵E的秩小于m，那么这m个极限点够不成一个平面，有时候就算构造出里超平面，也不一定能得到某些坐标轴的截距，这种情况下我们设置截距为在此目标上的最大值。

说明：原文里算法伪代码中的zjminz_j^{min}zjmin的计算有问题，应该是我文中所说的∪τ∣=0tSτ\cup_{\tau|=0}^tS_{\tau}∪τ∣=0tSτ的最小值，而不仅仅是StS_tSt的最小值。

4.联系个体和参考点

这一步实际上就是找到每个个体距离最近的参考点的参考线，如图6所示。

图6 说明了种群个体与参考点之间的关系

这一步的算法伪代码如下：

5.Niche-Preservation操作

根据步骤4，我们可以发现参考点和种群个体之间的关系可以是一对多，一对一或者没有对应。这一步我们要做的工作就是怎样通过这个关系选出剩余的K个解进入下一个种群。其伪代码如下：

这一步骤中有三点需要解释：

第一，这里的ρj\rho_jρj是指第j个参考点的Niche数，第j个参考点的Niche数是种群Pt+1=St/FlP_{t+1}=S_t/F_lPt+1=St/Fl中和这个参考点相联系的个体数目。

第二，这里挑选出K个解的原则是：对于那些联系少的参考点对应的个体更应该被保留从而维持多样性。因此算法步骤中JminJ_{min}Jmin是联系少的参考点，但其可能不止一个，比如AB两个参考点分别都联系了0个个体，那么j‾\overline{j}j就从JminJ_{min}Jmin中任取一个。Ij‾I_{\overline{j}}Ij就是参考点j‾\overline{j}j所联系的个体，如果Ij‾I_{\overline{j}}Ij为空，那么重新选择参考点；如果Ij‾I_{\overline{j}}Ij不为空，那么看ρj‾\rho_{\overline{j}}ρj是否为0，如果ρj‾=0\rho_{\overline{j}}=0ρj=0，从Ij‾I_{\overline{j}}Ij选择距离j‾\overline{j}j最小的个体进入下一代，如果ρj‾≠0\rho_{\overline{j}}\neq0ρj̸=0，那么从Ij‾I_{\overline{j}}Ij随机选择一个个体进入下一代即可，因为此参考点的多样性已经得到“满足”了。

整个选择过程的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
#求两个向量矩阵的余弦值,x的列数等于y的列数
def pdist(x,y):x0=x.shape[0]y0=y.shape[0]xmy=np.dot(x,y.T)#x乘以yxm=np.array(np.sqrt(np.sum(x**2,1))).reshape(x0,1)ym=np.array(np.sqrt(np.sum(y**2,1))).reshape(1,y0)xmmym=np.dot(xm,ym)cos = xmy/xmmymreturn cosdef lastselection(popfun1,popfun2,K,Z,Zmin):#选择最后一个front的解popfun = copy.deepcopy(np.vstack((popfun1, popfun2)))-np.tile(Zmin,(popfun1.shape[0]+popfun2.shape[0],1))N,M = popfun.shape[0],popfun.shape[1]N1 = popfun1.shape[0]N2 = popfun2.shape[0]NZ = Z.shape[0]#正则化extreme = np.zeros(M)w = np.zeros((M,M))+1e-6+np.eye(M)for i in range(M):extreme[i] = np.argmin(np.max(popfun/(np.tile(w[i,:],(N,1))),1))#计算截距extreme = extreme.astype(int)#python中数据类型转换一定要用astype#temp = np.mat(popfun[extreme,:]).Itemp = np.linalg.pinv(np.mat(popfun[extreme,:]))hyprtplane = np.array(np.dot(temp,np.ones((M,1))))a = 1/hyprtplaneif np.sum(a==math.nan) != 0:a = np.max(popfun,0)np.array(a).reshape(M,1)#一维数组转二维数组#a = a.T - Zmina=a.Tpopfun = popfun/(np.tile(a,(N,1)))##联系每一个解和对应向量#计算每一个解最近的参考线的距离cos = pdist(popfun,Z)distance = np.tile(np.array(np.sqrt(np.sum(popfun**2,1))).reshape(N,1),(1,NZ))*np.sqrt(1-cos**2)#联系每一个解和对应的向量d = np.min(distance.T,0)pi = np.argmin(distance.T,0)#计算z关联的个数rho = np.zeros(NZ)for i in range(NZ):rho[i] = np.sum(pi[:N1] == i)#选出剩余的K个choose = np.zeros(N2)choose = choose.astype(bool)zchoose = np.ones(NZ)zchoose = zchoose.astype(bool)while np.sum(choose) < K:#选择最不拥挤的参考点temp = np.ravel(np.array(np.where(zchoose == True)))jmin = np.ravel(np.array(np.where(rho[temp] == np.min(rho[temp]))))j = temp[jmin[np.random.randint(jmin.shape[0])]]
#        I = np.ravel(np.array(np.where(choose == False)))
#        I = np.ravel(np.array(np.where(pi[(I+N1)] == j)))I = np.ravel(np.array(np.where(pi[N1:] == j)))I = I[choose[I] == False]if (I.shape[0] != 0):if (rho[j] == 0):s = np.argmin(d[N1+I])else:s = np.random.randint(I.shape[0])choose[I[s]] = Truerho[j] = rho[j]+1else:zchoose[j] = Falsereturn choosedef envselect(mixpop,N,Z,Zmin,name,M,D):#非支配排序mixpopfun = cal(mixpop,name,M,D)frontno,maxfno = NDsort(mixpopfun,N,M)Next = frontno < maxfno#选择最后一个front的解Last = np.ravel(np.array(np.where(frontno == maxfno)))choose = lastselection(mixpopfun[Next,:],mixpopfun[Last,:],N-np.sum(Next),Z,Zmin)Next[Last[choose]] = True#生成下一代pop = copy.deepcopy(mixpop[Next,:])return pop

6.遗传算子

为了保证遗传算子的高效性，这里的SBX使用了一个较大的指数分布。
遗传算子的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
def GO(pop,t1,t2,pc,pm):pop1 = copy.deepcopy(pop[0:int(pop.shape[0]/2),:])pop2 = copy.deepcopy(pop[(int(pop.shape[0]/2)):(int(pop.shape[0]/2)*2),:])N,D = pop1.shape[0],pop1.shape[1]#模拟二进制交叉beta=np.zeros((N,D))mu=np.random.random_sample([N,D])beta[mu<=0.5]=(2*mu[mu<=0.5])**(1/(t1+1))beta[mu>0.5]=(2-2*mu[mu>0.5])**(-1/(t1+1))beta=beta*((-1)**(np.random.randint(2, size=(N,D))))beta[np.random.random_sample([N,D])<0.5]=1beta[np.tile(np.random.random_sample([N,1])>pc,(1,D))]=1off = np.vstack(((pop1+pop2)/2+beta*(pop1-pop2)/2,(pop1+pop2)/2-beta*(pop1-pop2)/2))#多项式变异low=np.zeros((2*N,D))up=np.ones((2*N,D))site=np.random.random_sample([2*N,D]) < pm/Dmu = np.random.random_sample([2*N,D])temp = site & (mu<=0.5)off[off<low]=low[off<low]off[off>up]=up[off>up]off[temp]=off[temp]+(up[temp]-low[temp])*((2*mu[temp]+(1-2*mu[temp])*((1-(off[temp]-low[temp])/(up[temp]-low[temp]))**(t2+1)))**(1/(t2+1))-1)temp = site & (mu>0.5)off[temp]=off[temp]+(up[temp]-low[temp])*(1-(2*(1-mu[temp])+2*(mu[temp]-0.5)*((1-(up[temp]-off[temp])/(up[temp]-low[temp]))**(t2+1)))**(1/(t2+1)))return off

7.时间复杂度

假设算法中∣Fl∣=L|F_l|=L∣Fl∣=L，并且L≤2NL\leq2NL≤2N，并且要求N≈P,N>MN\approx P,N>MN≈P,N>M，那么NSGAIII一代的最坏情况下的时间复杂度为O(N2M)O(N^2M)O(N2M)或者O(N2logM−2N)O(N^2log^{M-2}N)O(N2logM−2N)，具体分析可以参考原论文。其中非支配排序里大小为2N目标函数个为M的种群时间复杂度为O(NlogM−2N)O(Nlog^{M-2}N)O(NlogM−2N)。这个分析见文献[6]。

8.NSGAIII参数说明

这个算法的一大优势在于参数较少，只有种群大小，终止参数，交叉变异参数等进化算法的常用参数。有人要问参考点的数目不也是参数吗？其实参考点可以用户自定义也可以结构化生成，因此其数目不算NSGAIII参数。

III.实验

1.DTLZ1~DTLZ4

此处给出DTLZ1~DTLZ4的数学定义式或者说明：

DTLZ1:

DTLZ2:

DTLZ3:

DTLZ3实际上是结合了DTLZ1的g(xM)g(x_M)g(xM)和DTLZ2的优化目标。

DTLZ4:

PF和适应度函数计算的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
def funfun(M,N,name):#种群初始化D=M+4#定义自变量个数为目标个数加4low=np.zeros((1,D))up=np.ones((1,D))pop = np.tile(low,(N,1))+(np.tile(up,(N,1))-np.tile(low,(N,1)))*np.random.rand(N,D)#计算PFif name=='DTLZ1':#g=np.transpose(np.mat(100*(D-M+1+np.sum(((pop[:,(M-1):]-0.5)**2-np.cos(20*np.pi*(pop[:,(M-1):]-0.5))),1))))g=np.array(100*(D-M+1+np.sum(((pop[:,(M-1):]-0.5)**2-np.cos(20*np.pi*(pop[:,(M-1):]-0.5))),1))).reshape(N,1)popfun=np.multiply(0.5*np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),pop[:,:(M-1)]))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.fliplr(pop[:,:(M-1)])))))P,nouse = uniformpoint(N,M)P=P/2elif name=='DTLZ2':#g=np.transpose(np.mat(np.sum((pop[:,(M-1):]-0.5)**2,1)))g=np.array(np.sum((pop[:,(M-1):]-0.5)**2,1)).reshape(N,1)popfun=np.multiply(np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),np.cos(pop[:,:(M-1)]*(np.pi/2))))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.sin(np.fliplr(pop[:,:(M-1)])*(np.pi/2))))))P,nouse = uniformpoint(N,M)        #P = P/np.tile(np.transpose(np.mat(np.sqrt(np.sum(P**2,1)))),(1,M))P = P/np.tile(np.array(np.sqrt(np.sum(P**2,1))).reshape(P.shape[0],1),(1,M))elif name=='DTLZ3':g=np.array(100*(D-M+1+np.sum(((pop[:,(M-1):]-0.5)**2-np.cos(20*np.pi*(pop[:,(M-1):]-0.5))),1))).reshape(N,1)popfun=np.multiply(np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),np.cos(pop[:,:(M-1)]*(np.pi/2))))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.sin(np.fliplr(pop[:,:(M-1)])*(np.pi/2))))))P,nouse = uniformpoint(N,M)        #P = P/np.tile(np.transpose(np.mat(np.sqrt(np.sum(P**2,1)))),(1,M))P = P/np.tile(np.array(np.sqrt(np.sum(P**2,1))).reshape(P.shape[0],1),(1,M))return pop,popfun,P,Ddef cal(pop,name,M,D):N = pop.shape[0]if name=='DTLZ1':g=np.array(100*(D-M+1+np.sum(((pop[:,(M-1):]-0.5)**2-np.cos(20*np.pi*(pop[:,(M-1):]-0.5))),1))).reshape(N,1)popfun=np.multiply(0.5*np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),pop[:,:(M-1)]))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.fliplr(pop[:,:(M-1)])))))elif name=='DTLZ2':g=np.array(np.sum((pop[:,(M-1):]-0.5)**2,1)).reshape(N,1)popfun=np.multiply(np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),np.cos(pop[:,:(M-1)]*(np.pi/2))))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.sin(np.fliplr(pop[:,:(M-1)])*(np.pi/2))))))elif name=='DTLZ3':g=np.array(100*(D-M+1+np.sum(((pop[:,(M-1):]-0.5)**2-np.cos(20*np.pi*(pop[:,(M-1):]-0.5))),1))).reshape(N,1)popfun=np.multiply(np.tile(1+g,(1,M)),(np.fliplr((np.hstack((np.ones((g.shape[0],1)),np.cos(pop[:,:(M-1)]*(np.pi/2))))).cumprod(1))))popfun=np.multiply(popfun,(np.hstack((np.ones((g.shape[0],1)),1-np.sin(np.fliplr(pop[:,:(M-1)])*(np.pi/2))))))return popfun

2.评价指标

这里选用比较IGD指标进行评价，假设A是找到的非支配解集，Z是目标解集，那么IGD定义如下：
IGD(A,Z)=1∣Z∣∑i=1∣z∣minj=1to∣A∣d(zi,aj),d(zi,aj)=∣∣zi−aj∣∣2IGD(A,Z)=\frac{1}{|Z|}\sum^{|z|}_{i=1}min_{j=1to|A|}d(z_i,a_j),d(z_i,a_j)=||z_i-a_j||_2 IGD(A,Z)=∣Z∣1i=1∑∣z∣minj=1to∣A∣d(zi,aj),d(zi,aj)=∣∣zi−aj∣∣2
其值是越小越好，此指标衡量了收敛性和多样性。其更多研究参见文献[7]。

IGD计算的python代码如下：

from scipy.special import comb
from itertools import combinations
import numpy as np
import copy
import math
def EuclideanDistances(A, B):BT = B.transpose()# vecProd = A * BTvecProd = np.dot(A,BT)# print(vecProd)SqA =  A**2# print(SqA)sumSqA = np.matrix(np.sum(SqA, axis=1))sumSqAEx = np.tile(sumSqA.transpose(), (1, vecProd.shape[1]))# print(sumSqAEx)SqB = B**2sumSqB = np.sum(SqB, axis=1)sumSqBEx = np.tile(sumSqB, (vecProd.shape[0], 1))    SqED = sumSqBEx + sumSqAEx - 2*vecProdSqED[SqED<0]=0.0   ED = np.sqrt(SqED)return EDdef IGD(popfun,PF):distance = np.min(EuclideanDistances(PF,popfun),1)score = np.mean(distance)return score

3.实验设置与结果

这里种群大小为100，交叉概率为1，变异概率为1/n1/n1/n，交叉参数为30，变异参数为20。我们分别测试了M=3,5,8,10,15M=3,5,8,10,15M=3,5,8,10,15这5组实验，每组实验分别在DTLZ1~DTLZ4上进行，每次实验运行20次，计算IGD指标的均值，其结果如表1所示。

表1 NSGAIII在DTLZ1~DTLZ3上进行多组实验得到的IGD均值

问题	M	最大迭代次数	IGD均值（标准差）
DTLZ1	3	400	2.0667e-2 (1.16e-4)
DTLZ1	5	600	6.8250e-2 (1.39e-4)
DTLZ1	8	750	1.2004e-1 (3.74e-2)
DTLZ1	10	1000	1.9666e-1 (5.65e-2)
DTLZ1	15	1500	3.2579e-1 (6.39e-2)
DTLZ2	3	250	5.4490e-2 (1.34e-5)
DTLZ2	5	350	2.1231e-1 (4.53e-5)
DTLZ2	8	500	4.3045e-1 (8.08e-2)
DTLZ2	10	750	6.6150e-1 (7.49e-2)
DTLZ2	15	1000	9.3023e-1 (3.22e-2)
DTLZ3	3	1000	5.4733e-2 (2.43e-4)
DTLZ3	5	1000	2.1356e-1 (1.03e-3)
DTLZ3	8	1000	8.1274e-1 (5.96e-1)
DTLZ3	10	1500	2.6353e+0 (2.34e+0)
DTLZ3	15	2000	1.1322e+0 (4.25e-2)

由表可知，算法在DTLZ1~DTLZ3上较好的平衡了收敛性和多样性。

以下画出DTLZ1~DTLZ3在M=3和M=8时得到的非支配集合的分布图。

图7 NSGAIII在M=3的DTLZ1问题上得到的非支配解集分布图

图8 NSGAIII在M=3的DTLZ2问题上得到的非支配解集分布图

图9 NSGAIII在M=3的DTLZ3问题上得到的非支配解集分布图

图10 NSGAIII在M=8的DTLZ1问题上得到的非支配解集分布图

图11 NSGAIII在M=8的DTLZ2问题上得到的非支配解集分布图

图12 NSGAIII在M=8的DTLZ3问题上得到的非支配解集分布图

由上述分布图可知，不管是M=3还是更高的维度M=8，得到的非支配解集都具有良好的收敛性和多样性。

IV 总结

此次博客主要介绍了NSGAIII。对里面的难点进行了解释和总结。对于多个优化目标的问题而言，其确实具有良好的搜索高维空间中解的能力。在算法实现的一些细节中，借鉴了很多其它论文，比如PF的生成[3]，参考点的生成[3]，截距不存在[5]等。

参考文献

[1] Deb K , Jain H . An Evolutionary Many-Objective Optimization Algorithm Using Reference Point-Based Nondominated Sorting Approach, Part I: Solving Problems With Box Constraints[J]. IEEE Transactions on Evolutionary Computation, 2014, 18(4):577-601.

[2] Deb, Kalyanmoy & Thiele, Lothar & Laumanns, Marco & Zitzler, Eckart. (2002). Scalable Multi-Objective Optimization Test Problems.

[3] Tian Y , Xiang X , Zhang X , et al. Sampling Reference Points on the Pareto Fronts of Benchmark Multi-Objective Optimization Problems[C]// 2018 IEEE Congress on Evolutionary Computation (CEC). IEEE, 2018.

[4] Sun Y , Yen G G , Yi Z . IGD Indicator-based Evolutionary Algorithm for Many-objective Optimization Problems[J]. IEEE Transactions on Evolutionary Computation, 2018.

[5] 袁源. 基于分解的多目标进化算法及其应用[D]

[6] H. T. Kung, F. Luccio, and F. P. Preparata, “On finding the maxima of a set of vectors,” Journal of the Association for Computing Machinery, vol. 22, no. 4, pp. 469–476, 1975.

[7] Schütze, Oliver & Esquivel, Xavier & Lara, Adriana & A Coello Coello, Carlos. (2019). Measuring the Averaged Hausdorff Distance to the Pareto Front of a Multi-Objective Optimization Problem.