NSGA-III 算法实现与测试

1 绪论

1.1 多目标优化问题

• 在实际问题中往往具有多个目标且需要同时满足，即在同一问题模型中同时存在几个非线性目标，而这些目标函数需要同时进行优化处理，并且这些目标又往往是互相冲突的，称这类问题称为多目标优化问题。设有 $m$ 个目标函数， $n$ 维决策变量，多目标优化问题可表示为：

  $$\begin{split} min\;y=&F(x)=(f_1(x),f_2(x),...,f_m(x))\\ subject\;to: &g_i(x)\le 0,i=1,2,...,p\\ &h_j(x)\le 0,i=1,2,...,q\\ &x\in D \end{split}$$

• 其中 $x\in D$ 为 $n$ 维决策变量， $x=(x_1,x_2,…,x_n)$ ， $x_i$ 表示第 $i$ 个决策变量， $D$ 为 $n$ 维决策空间。 $y\in S$ 为 $m$ 维目标变量， $S$ 为 $m$ 维解空间。 $F(x)$ 定义了 $m$ 个由决策空间向解空间映射的目标分量函数，其中 $f_i(x)$ 是目标函数 $F(x)$ 的第 $i$ 个目标分量。 $g_i(x)\le 0(i=1,2,…,p)$ 定义了 $p$ 个不等式约束条件， $h_i(x)\le 0(i=1,2,…,p)$ 定义了 $q$ 个等式约束条件。 $min$ 或 $max$ 表示多目标问题的优化模式，一般情况下优化问题是转化成最小优化模式，即 $min$ 。

• 在对多个目标的处理方法上，主要存在以下两种方式，一种是通过加权的方式将多目标优化问题转化为单目标优化问题，一种是利用 $Pareto$ 支配关系求得一组互不支配的解。第一种方法需要预先制定目标的权重向量，并且求得最终结果受其影响很大，每当权重发生变化，算法都需要重新运行。而在实际应用当中，权重的值大多是无法预知的，并且用户的偏好也会随着环境的变化而发生变化。因此，目前普遍采用第二种方法求得一组互不支配的解集，然后在其基础上根据实际情况进行进一步的决策。

• $Pareto$ 支配，又称 $Pareto$ 占优，设 $x_A,x_B\in X_f$ 是多目标优化问题中的两个可行解，则称与 $x_B$ 相比， $x_A$ 是 $Pareto$ 占优的，当且仅当： $$\forall i=1,2,...,m;f_i(x_A)\le f(x_B)$$ 且至少存在一个 $j\in {1,2,…,m}$ ，使得 $f_i(x_A)<f(x_B)$ 。记为 $x_A\prec x_B$ ， $x_A$ 支配 $x_B$ 。
• 可行解 $x_A\in X_f$ 称为 $Pareto$ 最优解(或非支配解)，当且仅当： $$\neg \exists x\in X_f,x\prec x_A$$ 所有 $Pareto$ 最优解构成的集合称为 $Pareto$ 最优解集。由 $Pareto$ 最优解集中所有 $Pareto$ 最优解所对应的目标函数向量组成的曲面称为 $Pareto$ 前沿。
• 在求解多目标优化问题时，优化目标函数中若存在两个相互冲突的目标分量，那么在优化解空间中是不可能获取一个全局最优的点，而只能获取一组点的集合，即上面所定义的 $Pareto$ 最优解集，然后通过用户的偏好选取一个或一组 $Pareto$ 最优解，即满意解。

1.2 多目标进化算法

• 多目标进化算法 $(MOEA)$ 是一类模拟生物进化机制而形成的全局性概率优化搜索方法，在20世纪90年代中期开始迅速发展，其发展可以分为两个阶段。第一阶段主要有两种方法即不基于 $Pareto$ 优化的方法和基于 $Pareto$ 优化的方法。第二个阶段就是在此基础上提出了外部集这个概念，外部集存放的是当前代的所有非支配个体，从而使解集保持较好的分布度。这个时期提出的多目标进化算法更多地强调算法的效率和有效性。在这两个阶段中，比较典型的多目标进化算法有 $NSGA2$ 、 $PESA2$ 和 $SPEA2$ 。对于这三种算法而言，其优点较多但是其缺点也比较明显的。如 $NSGA2$ 的优点在于运行效率高、解集有良好的分布性，特别对于低维优化问题具有较好的表现，其缺点在于在高维问题中解集过程具有缺陷，解集的多样性不理想。 $PESA2$ 的优点在于其解的收敛性很好，比较容易接近最优面，特别是在高维问题情况下，但其不足之处在于选择操作一次只能选取一个个体，时间消耗很大，而且阶级的多样性不佳。 $SPEA2$ 的优点在于可以取得一个分布度很好的解集，特别是在高维问题的求解上，但是其聚类过程保持多样性耗时较长，运行效率不高。
• 多目标进化算法的基本原理描述如下：多目标进化算法从一组随机生成的种群出发，通过对种群执行选择、交叉和变异等进化操作，经过多代进化，种群中个体的适应度不断提高，从而逐步逼近多目标优化问题的 $Pareto$ 最优解集。与单目标进化算法不同，多目标进化算法具有特殊的适应度评价机制。为了充分发挥进化算法的群体搜索优势，大多数 $(MOEA)$ 均采用基于 $Pareto$ 排序的适应度评价方法。在实际应用中，为使算法更好地收敛到多目标优化问题的 $Pareto$ 最优解，现有的 $(MOEA)$ 通常还采用了精英策略、小生境和设置外部集等关键技术。
• $(MOEA)$ 一般框架所描述的算法思想如下： $(MOEA)$ 通过对种群 $X(t)$ 执行选择、交叉和变异等操作产生下一代种群 $X(t+1)$ 。在每一代进化过程中，首先将种群 $X(t)$ 中的所有非劣解个体都复制到外部集 $A(t)$ 中，然后运用小生境截断算子剔除 $A(t)$ 中的劣解和一些距离较近的非劣解个体，以得到个体分布更为均匀的下一代外部集 $A(t+1)$ ，并且按照概率 $p_e$ 从 $A(t+1)$ 中选择一定数量的优秀个体进入下代种群。在进化结束时，将外部集中的非劣解个体作为最优解输出。目前， $(MOEA)$ 研究取得了大量成果，已被应用于许多领域，如工程领域、工业领域和科学领域。

2 NSGA-III 算法实现

• $NSGA$ (非支配排序遗传算法)系列算法包括 $NSGA$ 、 $NSGA-II$ 和 $NSGA-III$ ，每一代算法都是对上一代算法的优化，提升了运算效率。在此仅对这三个算法进行简单介绍，具体算法细节恕不赘述。

2.1 NSGA 算法原理

• $NSGA$ 与简单的遗传算法的主要区别在于：该算法在选择算子执行之前根据个体之间的支配关系进行了分层。其选择算子、交叉算子和变异算子与简单遗传算法没有区别。
• 算法流程如下，算法首先判断种群是否全部分级，如果已经全部分级，则在分级的基础上，使用基于拥挤策略的小生境技术对虚拟适应度值进行调整，并确定每个种群的虚拟适应度值，然后根据虚拟适应度值的大小，确定优先选择进行处理的种群。

• $NSGA$ 采用非支配分层方法,可以使好的个体有更大的机会遗传到下一代，适应度共享策略则使得准 $Pareto$ 面上的个体均匀分布,保持了群体多样性,克服了超级个体的过度繁殖,防止了早熟收敛。但 $NSGA$ 算法仍存在以下缺陷：
  • 计算复杂度较高，为 $O(mN^3)$ ( $m$ 为目标函数个数， $N$ 为种群大小)，所以当种群较大时，计算相当耗时。
  • 没有精英策略。精英策略可以加速算法的执行速度，而且也能在一定程度上确保己经找到的满意解不被丢失。
  • 需要指定小生境半径，泛化能力较差。

2.2 NSGA-II 算法原理

• $NSGA-II$ 是 $NSGA$ 的改进版本，在 $NSGA$ 的基础上，针对上述缺陷添加了一些策略，使算法效率显著提高。主要改进如下：
  • 提出了快速非支配排序法，降低了算法的计算复杂度。由原来的 $O(mN^3)$ 降到 $O(mN^2)$ ，其中， $m$ 为目标函数个数， $N$ 为种群大小。
  • 引入精英策略，扩大采样空间。将父代种群与其产生的子代种群组合，共同竞争产生下一代种群，有利于保持父代中的优良个体进入下一代，并通过对种群中所有个体的分层存放，使得最佳个体不会丢失，迅速提高种群水平。
  • 提出了拥挤度和拥挤度比较算子，代替了需要指定小生境半径的适应度共享策略，并在快速排序后的同级比较中作为胜出标准，使准 $Pareto$ 域中的个体能扩展到整个 $Pareto$ 域，并均匀分布，保持了种群的多样性。
• 然而，基于拥挤距离的选择机制仅在双目标优化问题中表现较好，当面对三个及其以上目标的多目标优化问题时，如果继续采用拥挤距离的话，算法的收敛性和多样性将受到影响，即得到的解在非支配层上分布不均匀，导致算法陷入局部最优。

2.3 NSGA-III 算法原理

• 为了解决三个及其以上目标的多目标优化问题， $NSGA-III$ 应运而生。拥挤距离度量并不适合求解三个及更多目标的多目标优化问题，因此 $NSGA-III$ 不再采用拥挤距离，而是采用了新的选择机制，该机制会通过所提供的参考点，对种群中的个体进行更加系统地分析，以选择适应度更高的解进入下一代。
• 本文实现的是 $NSGA-III$ 算法，代码见文末。

3 NSGA-III 算法测试

3.1 测试问题

• 考虑到展示效果，仅对双目标和三目标问题进行测试。测试问题选择简单双目标问题、DTLZ系列问题和WFG系列问题。
• 发表在IEEE TEVC 2018，来自 Zapotecas-Martínez, C. A. C. Coello, H.E. Aguirre 和K.Tanaka的文章： $A\;Review\;of\;Features\;and\;Limitations\;of\;Existing\;Scalable\;Multi-Objective\;Test\;Suites$ ，对现有的多目标问题集做了一个综述。该文章不仅介绍了设计多目标演化问题集的方法，还介绍了多目标问题的特征要素，问题特征要素如图3.1所示。另外，文章中还做了一个较为详细的多目标演化问题集相关的综述，并指出不同问题集拥有的不同特征。该文章中详细总结了9个问题集：分别为DTLZ、WFG、LSS、SZDT、LSMOP、MNI、minus-DTLZ \& minus-WFG、ZDT以及用于视觉测试的问题集。本文选取其中的DTLZ和WFG两个系列的测试问题集，计算公式均参考自该文章并直接截取原文中的表格。

3.1.1 简单双目标问题

$$\begin{aligned} &min\;f_1(x1,x2)=x_1^4-10x_1^2+x_1x_2+x_2^4-x_1^2x_2^2\\ &min\;f_2(x1,x2)=x_1^4+x_1x_2+x_2^4-x_1^2x_2^2\\ &s.t.\begin{cases} -5\le x_1\le 5 \\ -5\le x_2\le 5 \end{cases} \end{aligned}$$

3.1.2 DTLZ 系列问题

• Deb等人在2002年首次提出DTLZ测试函数集，并以共同研究者名字首字母命名(Deb，Thiele，Laumanns，Zitzler)，根据不同难度的设置期望，2005年Deb等又在原有七个函数的基础上加入了两个测试函数，共同组成了DTLZ测试函数集。DTLZ测试函数集可以扩展至任意数量的目标 $(m\ge 2)$ 并且可以有任意数目个变量 $(n\ge m)$ 。因为变量数与目标数易于控制，所以DTLZ函数集被广泛应用于多目标优化问题当中作为标准测试函数。本次实验仅测试前七个DTLZ测试函数。

3.1.3 WFG 系列问题

• WFG 测试函数集是Huband等人在2006年提出来的，一共包含9个测试问题，这些问题的目标数目都可以扩展到任意数量。和DTLZ测试函数集比起来，DTLZ问题的变量都是可分离的，因此复杂程度不高，而WFG测试问题的复杂程度更高、处理起来更具有挑战性。WFG测试问题的属性包括可分性或者不可分性、单峰或者多峰、PF形状为凸或者非凸、无偏差参数或有偏差参数。WFG测试函数集可以提供更有效的依据来评估优化算法在各种不同问题上的表现性能。

3.2 测试结果

• 测试结果以散点图形式展示，其中双目标优化问题采用单个二维散点图，三目标优化问题采用两个不同角度的三维散点图。DTLZ和WFG系列问题的 $Pareto$ 前沿参考自Ye Tian、Ran Cheng、Xingyi Zhang和Yaochu Jin发表的文章： $PlatEMO:\;A\;MATLAB\;Platform\;for\;Evolutionary\;Multi-Objective\;Optimization$ ，下面第一张图为原文中截取的图片，作为直观评价测试结果的标准，本节其余图片为测试结果。由于程序运行耗时较长，因此测试中限定双目标优化问题的迭代次数为300代，三目标优化问题的迭代次数为500代，迭代次数均标记在每张图片上方。
• 测试结果反映了算法的拟合性能，同时也反映了不同问题收敛难度的区别。通过与真实 $Pareto$ 前沿的直观对比，可以看出如DTLZ2、DTLZ3等问题的拟合难度比较低，在500次迭代后拟合结果已经十分理想。相反，如DTLZ1、WFG4等问题的你和难度比较高，在500次迭代后拟合结果仍比较粗陋，表明算法的性能有待提高。

4 源码

• 编程环境：python3.6.4 + numpy1.14.3 + scipy1.1.0 + matplotlib2.2.2
• 使用说明：运行前修改代码最后一行 main 函数的参数。其中 Problem 是待测问题，可选择的有 ‘test’ 、 ‘DTLZ1’ ~ ‘DTLZ7’ 、 ‘WFG1’ ~ ‘WFG9’ 。 M 是目标个数，可设为2或3。 Generations 是算法迭代次数，可设为任意正整数。

# -*- coding: utf-8 -*-  
from scipy.special import comb
from itertools import combinations
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def T_weight(H,M):
    N = int(comb(H+M-1,M-1))
    tmp1 = np.array(list(combinations(np.arange(1,H+M),M-1)))
    tmp2 = np.tile(np.arange(0,M-1),(int(comb(H+M-1,M-1)),1))
    Temp = tmp1-tmp2-1
    W = np.zeros((N,M))
    W[:,0] = Temp[:,0]-0
    for i in range(1,M-1):
        W[:,i] = Temp[:,i]-Temp[:,i-1]
    W[:,-1] = H-Temp[:,-1]
    W = W/H
    return N,W

def F_weight(p1,p2,M):
    N,W = T_weight(p1,M)
    if p2 > 0:
        N2,W2 = T_weight(p2,M)
        N = N+N2
        W = np.concatenate([W,W2/2+1/(2*M)])
    return N,W

def F_objective(Operation,p,M,Input):
    if p=='test':
        Output,Boundary = F_test(Operation,p,M,Input)
    elif p[:3]=='WFG':
        Output,Boundary = F_WFG(Operation,p,M,Input)
    elif p[:4]=='DTLZ':
        Output,Boundary = F_DTLZ(Operation,p,M,Input)
    return Output,Boundary

def F_test(Operation,Problem,M,Input):
    Boundary = None
    if Operation=='init':
        D = 10
        MaxValue   = np.ones((1,D))*5
        MinValue   = np.ones((1,D))*(-5)
        Population = np.random.rand(Input,D)
        Population = Population*np.tile(MaxValue,(Input,1))+(1-Population)*np.tile(MinValue,(Input,1))
        Output   = Population
        Boundary = np.array([MaxValue,MinValue])
    elif Operation=='value':
        Population = Input
        FunctionValue = np.zeros((Population.shape[0],M))
        x1 = Population[:,0]
        x2 = Population[:,1]
        FunctionValue[:,0] = x1**4-10*x1**2+x1*x2+x2**4-(x1**2)*(x2**2)
        FunctionValue[:,1] = x2**4-(x1**2)*(x2**2)+x1**4+x1*x2
        Output = FunctionValue
    return Output,Boundary

def F_DTLZ(Operation,Problem,M,Input):
    Boundary = None
    K = [5,10,10,10,10,10,20]
    K = K[int(Problem[-1])-1]
    if Operation=='init':
        D = M+K-1
        MaxValue   = np.ones((1,D))
        MinValue   = np.zeros((1,D))
        Population = np.random.rand(Input,D)
        Population = Population*np.tile(MaxValue,(Input,1))+(1-Population)*np.tile(MinValue,(Input,1))
        Output   = Population
        Boundary = np.array([MaxValue,MinValue])
    elif Operation=='value':
        Population = Input
        FunctionValue = np.zeros((Population.shape[0],M))
        if Problem == 'DTLZ1':
            g = 100*(K+np.sum((Population[:,M-1:]-0.5)**2-np.cos(20*np.pi*(Population[:,M-1:]-0.5)),axis=1))
            for i in range(M):
                FunctionValue[:,i] = 0.5*np.prod(Population[:,0:M-i-1],axis=1)*(1+g)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*(1-Population[:,M-i-1])
        elif Problem == 'DTLZ2':
            g = np.sum((Population[:,M-1:]-0.5)**2,axis=1)
            for i in range(M):
                FunctionValue[:,i] = (1+g)*np.prod(np.cos(0.5*np.pi*Population[:,0:M-i-1]),axis=1)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*np.sin(0.5*np.pi*Population[:,M-i-1])
        elif Problem == 'DTLZ3':
            g = 100*(K+np.sum((Population[:,M-1:]-0.5)**2-np.cos(20*np.pi*(Population[:,M-1:]-0.5)),axis=1))
            for i in range(M):
                FunctionValue[:,i] = (1+g)*np.prod(np.cos(0.5*np.pi*Population[:,0:M-i-1]),axis=1)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*np.sin(0.5*np.pi*Population[:,M-i-1])
        elif Problem == 'DTLZ4':
            Population[:,0:M-2] = Population[:,0:M-2]**100
            g = np.sum((Population[:,M-1:]-0.5)**2,axis=1)
            for i in range(M):
                FunctionValue[:,i] = (1+g)*np.prod(np.cos(0.5*np.pi*Population[:,0:M-i-1]),axis=1)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*np.sin(0.5*np.pi*Population[:,M-i-1])
        elif Problem == 'DTLZ5':
            g = np.sum((Population[:,M-1:]-0.5)**2,axis=1)
            Temp = np.tile(g,(1,M-2))
            if M>2:
                Temp = Temp.reshape(Temp.shape[1],Temp.shape[0])
            Population[:,1:M-1] = (1+2*Temp*Population[:,1:M-1])/(2+2*Temp)
            for i in range(M):
                FunctionValue[:,i] = (1+g)*np.prod(np.cos(0.5*np.pi*Population[:,0:M-i-1]),axis=1)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*np.sin(0.5*np.pi*Population[:,M-i-1])
        elif Problem == 'DTLZ6':
            g = np.sum(Population[:,M-1:]**0.1,axis=1)
            Temp = np.tile(g,(1,M-2))
            if M>2:
                Temp = Temp.reshape(Temp.shape[1],Temp.shape[0])
            Population[:,1:M-1] = (1+2*Temp*Population[:,1:M-1])/(2+2*Temp)
            for i in range(M):
                FunctionValue[:,i] = (1+g)*np.prod(np.cos(0.5*np.pi*Population[:,0:M-i-1]),axis=1)
                if i > 0:
                    FunctionValue[:,i] = FunctionValue[:,i]*np.sin(0.5*np.pi*Population[:,M-i-1])
        elif Problem == 'DTLZ7':
            g = 1+9*np.mean(Population[:,M-1:],axis=1)
            FunctionValue[:,0:M-1] = Population[:,0:M-1]
            if M>2:
                Temp = np.tile(g,(1,M-2))
                Temp = Temp.reshape(Temp.shape[1],Temp.shape[0])
            else:
                Temp = np.tile(g,(1,M-1))
            h = M-np.sum(FunctionValue[:,0:M-1]/(1+Temp)*(1+np.sin(3*np.pi*FunctionValue[:,0:M-1])),axis=1)
            FunctionValue[:,M-1] = (1+g)*h
        Output = FunctionValue
    return Output,Boundary

def F_WFG(Operation,Problem,M,Input):
    K = [4,4,6,8,10,0,7,0,9]
    K = K[M-2]
    L = 10
    Boundary = None
    if Operation == 'init':
        D = K+L
        MaxValue   = np.arange(1,D+1)*2
        MinValue   = np.zeros(D)
        Population = np.random.rand(Input,D)
        Population = Population*np.tile(MaxValue,(Input,1))+(1-Population)*np.tile(MinValue,(Input,1))
        Output   = Population
        Boundary = np.array([MaxValue,MinValue])
    elif Operation == 'value':
        Population = Input
        N = Population.shape[0]
        D = 1
        S = np.tile(np.arange(1,M+1)*2,(N,1))
        if Problem == 'WFG3':
            A = np.concatenate([np.ones((N,1)),np.zeros((N,M-2))],axis=1)
        else:
            A = np.ones((N,M-1))
        z01 = Population/np.tile(np.arange(1,Population.shape[1]+1)*2,(N,1))
        if Problem == 'WFG1':
            t1 = np.zeros((N,K+L))
            t1[:,0:K] = z01[:,0:K]
            t1[:,K:] = s_linear(z01[:,K:],0.35)
            t2 = np.zeros((N,K+L))
            t2[:,0:K] = t1[:,0:K]
            t2[:,K:] = b_flat(t1[:,K:],0.8,0.75,0.85)
            t3 = np.zeros((N,K+L))
            t3 = t2**0.02
            t4 = np.zeros((N,M))
            for i in range(M-1):
                t4[:,i] = r_sum(t3[:,i*K//(M-1):(i+1)*K//(M-1)],np.arange(2*(i*K//(M-1)),2*(i+1)*K//(M-1),2))
            t4[:,M-1] = r_sum(t3[:,K:K+L],np.arange(2*K,2*(K+L),2))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t4[:,M-1],A[:,i]],axis=0)*(t4[:,i]-0.5)+0.5
            x[:,M-1] = t4[:,M-1]
            h = convex(x)
            h[:,M-1] = mixed(x)
        elif Problem == 'WFG2':                    
            t1 = np.zeros((N,K+L))
            t1[:,0:K] = z01[:,0:K]
            t1[:,K:] = s_linear(z01[:,K:],0.35)
            t2 = np.zeros((N,K+L//2))
            t2[:,0:K] = t1[:,0:K]
            for i in range(K,K+L//2):
                t2[:,i] = r_nonsep(t1[:,K+2*(i-K)-2:K+2*(i-K)],2)
            t3 = np.zeros((N,M))
            for i in range(M-1):
                t3[:,i] = r_sum(t2[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t3[:,M-1] = r_sum(t2[:,K:K+L//2],np.ones(L//2))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t3[:,M-1],A[:,i]],axis=0)*(t3[:,i]-0.5)+0.5
            x[:,M-1] = t3[:,M-1]
            h = convex(x)
            h[:,M-1] = disc(x)
        elif Problem == 'WFG3':
            t1 = np.zeros((N,K+L))
            t1[:,0:K] = z01[:,0:K]
            t1[:,K:] = s_linear(z01[:,K:],0.35)
            t2 = np.zeros((N,K+L//2))
            t2[:,0:K] = t1[:,0:K]
            for i in range(K,K+L//2):
                t2[:,i] = r_nonsep(t1[:,K+2*(i-K)-2:K+2*(i-K)],2)
            t3 = np.zeros((N,M))
            for i in range(M-1):
                t3[:,i] = r_sum(t2[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t3[:,M-1] = r_sum(t2[:,K:K+L//2],np.ones(L//2))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t3[:,M-1],A[:,i]],axis=0)*(t3[:,i]-0.5)+0.5
            x[:,M-1] = t3[:,M-1]
            h = linear(x)
        elif Problem == 'WFG4':
            t1 = np.zeros((N,K+L))
            t1 = s_multi(z01,30,10,0.35)
            t2 = np.zeros((N,M))
            for i in range(M-1):
                t2[:,i] = r_sum(t1[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t2[:,M-1] = r_sum(t1[:,K:K+L],np.ones(L))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t2[:,M-1],A[:,i]],axis=0)*(t2[:,i]-0.5)+0.5
            x[:,M-1] = t2[:,M-1]
            h = concave(x)
            h[:,M-1] = disc(x)
        elif Problem == 'WFG5':
            t1 = np.zeros((N,K+L))
            t1 = s_decept(z01,0.35,0.001,0.05)
            t2 = np.zeros((N,M))
            for i in range(M-1):
                t2[:,i] = r_sum(t1[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t2[:,M-1] = r_sum(t1[:,K:K+L],np.ones(L))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t2[:,M-1],A[:,i]],axis=0)*(t2[:,i]-0.5)+0.5
            x[:,M-1] = t2[:,M-1]
            h = concave(x)
            h[:,M-1] = disc(x)
        elif Problem == 'WFG6':
            t1 = np.zeros((N,K+L))
            t1[:,0:K] = z01[:,0:K]
            t1[:,K:] = s_linear(z01[:,K:],0.35)
            t2 = np.zeros((N,M))
            for i in range(M-1):
                t2[:,i] = r_nonsep(t1[:,i*K//(M-1):(i+1)*K//(M-1)],2)
            t2[:,M-1] = r_nonsep(t1[:,K-1:],L)
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t2[:,M-1],A[:,i]],axis=0)*(t2[:,i]-0.5)+0.5
            x[:,M-1] = t2[:,M-1]
            h = concave(x)   
            h[:,M-1] = disc(x)    
        elif Problem == 'WFG7':
            t1 = np.zeros((N,K+L))
            for i in range(K):
                t1[:,i] = b_param(z01[:,i],r_sum(z01[:,i:],np.ones(K+L-i)),0.98/49.98,0.02,50)
            t1[:,K:] = z01[:,K:]
            t2 = np.zeros((N,K+L))
            t2[:,0:K] = t1[:,0:K]
            t2[:,K:] = s_linear(t1[:,K:],0.35)
            t3 = np.zeros((N,M))
            for i in range(M-1):
                t3[:,i] = r_sum(t2[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t3[:,M-1] = r_sum(t2[:,K:K+L],np.ones(L))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t3[:,M-1],A[:,i]],axis=0)*(t3[:,i]-0.5)+0.5
            x[:,M-1] = t3[:,M-1]
            h = concave(x)
            h[:,M-1] = disc(x)
        elif Problem == 'WFG8':
            t1 = np.zeros((N,K+L))
            t1[:,0:K] = z01[:,0:K]
            for i in range(K,K+L):
                t1[:,i] = b_param(z01[:,i],r_sum(z01[:,0:i-1],np.ones(i-1)),0.98/49.98,0.02,50)
            t2 = np.zeros((N,K+L))
            t2[:,0:K] = t1[:,0:K]
            t2[:,K:] = s_linear(t1[:,K:],0.35)
            t3 = np.zeros((N,M))
            for i in range(M-1):
                t3[:,i] = r_sum(t2[:,i*K//(M-1):(i+1)*K//(M-1)],np.ones(K//(M-1)))
            t3[:,M-1] = r_sum(t2[:,K:K+L],np.ones(L))
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t3[:,M-1],A[:,i]],axis=0)*(t3[:,i]-0.5)+0.5
            x[:,M-1] = t3[:,M-1]
            h = concave(x)
            h[:,M-1] = disc(x)
        elif Problem == 'WFG9':
            t1 = np.zeros((N,K+L))
            for i in range(K+L-1):
                t1[:,i] = b_param(z01[:,i],r_sum(z01[:,i:],np.ones(K+L-i)),0.98/49.98,0.02,50)
            t1[:,-1] = z01[:,-1]
            t2 = np.zeros((N,K+L))
            t2[:,0:K] = s_decept(t1[:,0:K],0.35,0.001,0.05)
            t2[:,K:] = s_multi(t1[:,K:],30,95,0.35)
            t3 = np.zeros((N,M))
            for i in range(M-1):
                t3[:,i] = r_nonsep(t2[:,i*K//(M-1):(i+1)*K//(M-1)],K//(M-1))
            t3[:,M-1] = r_nonsep(t2[:,K:],L)
            x = np.zeros((N,M))
            for i in range(M-1):
                x[:,i] = np.max([t3[:,M-1],A[:,i]],axis=0)*(t3[:,i]-0.5)+0.5
            x[:,M-1] = t3[:,M-1]
            h = concave(x)
            h[:,M-1] = disc(x)
        Output = np.tile(D*x[:,M-1].reshape((N,1)),(1,M))+S*h
    return Output,Boundary
        

def b_flat(y,A,B,C):
    Output = A+np.minimum(0,np.floor(y-B))*A*(B-y)/B-np.minimum(0,np.floor(C-y))*(1-A)*(y-C)/(1-C)
    Output = np.round(Output,-6)
    return Output

def b_param(y,Y,A,B,C):
    Output = y**(B+(C-B)*(A-(1-2*Y)*np.abs(np.floor(0.5-Y)+A)))
    return Output

def s_linear(y,A):
    Output = np.abs(y-A)/np.abs(np.floor(A-y)+A)
    return Output

def s_decept(y,A,B,C):
    Output = 1+(np.abs(y-A)-B)*(np.floor(y-A+B)*(1-C+(A-B)/B)/(A-B)+np.floor(A+B-y)*(1-C+(1-A-B)/B)/(1-A-B)+1/B)
    return Output

def s_multi(y,A,B,C):
    Output = (1+np.cos((4*A+2)*np.pi*(0.5-np.abs(y-C)/2./(np.floor(C-y)+C)))+4*B*(np.abs(y-C)/2./(np.floor(C-y)+C))**2)/(B+2)
    return Output

def r_sum(y,w):
    Output = np.sum(y*np.tile(w,(y.shape[0],1)),axis=1)/np.sum(w)
    return Output

def r_nonsep(y,A):
    Output = np.zeros(y.shape[0])
    for j in range(y.shape[1]):
        Temp = np.zeros(y.shape[0])
        for k in range(A-1):
            Temp = Temp+np.abs(y[:,j]-y[:,(j+k)%y.shape[1]])
        Output = Output+y[:,j]+Temp
    Output = Output/(y.shape[1]/A)/np.ceil(A/2)/(1+2*A-2*np.ceil(A/2))
    return Output

def linear(x):
    Output = np.zeros(x.shape)
    for i in range(x.shape[1]):
        Output[:,i] = np.prod(x[:,0:x.shape[1]-i-1],axis=1)
        if i > 0:
            Output[:,i] = Output[:,i]*(1-x[:,x.shape[1]-i-1])
    return Output

def convex(x):
    Output = np.zeros(x.shape)
    for i in range(x.shape[1]):
        Output[:,i] = np.prod(1-np.cos(x[:,0:-1-i]*np.pi/2),axis=1)
        if i > 0:
            Output[:,i] = Output[:,i]*(1-np.sin(x[:,x.shape[1]-i]*np.pi/2))
    return Output

def concave(x):
    Output = np.zeros(x.shape)
    for i in range(x.shape[1]):
        Output[:,i] = np.prod(np.sin(x[:,0:-1-i]*np.pi/2),axis=1)
        if i > 0:
            Output[:,i] = Output[:,i]*(np.cos(x[:,x.shape[1]-i]*np.pi/2))
    return Output

def mixed(x):
    Output = 1-x[:,0]-np.cos(10*np.pi*x[:,0]+np.pi/2)/10/np.pi
    return Output

def disc(x):
    Output = 1-x[:,0]*(np.cos(5*np.pi*x[:,0]))**2
    return Output

def F_mating(Population):
    N,D = Population.shape
    MatingPool = np.zeros((N,D))
    Rank = list(range(N))
    np.random.shuffle(Rank)
    Pointer = 1
    for i in range(0,N,2):
        k = np.zeros((1,2))
        for j in range(2):
            if Pointer >= N:
                np.random.shuffle(Rank)
                Pointer = 1
            k = Rank[Pointer-1:Pointer+1]
            Pointer = Pointer+2
        MatingPool[i,:] = Population[k[0],:]
        if i+1<N:
            MatingPool[i+1,:] = Population[k[1],:]
    return MatingPool

def F_generator(MatingPool,Boundary,MaxOffspring):
    N,D = MatingPool.shape
    ProC = 1    
    ProM = 1/D  
    DisC = 20   
    DisM = 20   
    Offspring = np.zeros((N,D))
    for i in range(0,N,2):
        beta = np.zeros((1,D))
        miu  = np.random.rand(1,D)
        beta[miu<=0.5] = (2*miu[miu<=0.5])**(1/(DisC+1))
        beta[miu>0.5]  = (2-2*miu[miu>0.5])**(-1/(DisC+1))
        beta = beta*(-1)**np.random.randint(2,size=(1,D))
        beta[np.random.rand(1,D)>ProC] = 1
        if i+1<N:
            Offspring[i,:]   = (MatingPool[i,:]+MatingPool[i+1,:])/2+beta*(MatingPool[i,:]-MatingPool[i+1,:])/2
            Offspring[i+1,:] = (MatingPool[i,:]+MatingPool[i+1,:])/2-beta*(MatingPool[i,:]-MatingPool[i+1,:])/2
    Offspring = Offspring[0:MaxOffspring,:]
    if MaxOffspring == 1:
        MaxValue = Boundary[0,:]
        MinValue = Boundary[1,:]
    else:
        MaxValue = np.tile(Boundary[0,:],(MaxOffspring,1))
        MinValue = np.tile(Boundary[1,:],(MaxOffspring,1))
    k = np.random.rand(MaxOffspring,D)
    miu = np.random.rand(MaxOffspring,D)
    Temp = (k<=ProM)*(miu<0.5)
    Offspring[Temp] = Offspring[Temp]+(MaxValue[Temp]-MinValue[Temp])*((2*miu[Temp]+(1-2*miu[Temp])*(1-(Offspring[Temp]-MinValue[Temp])/(MaxValue[Temp]-MinValue[Temp]))**(DisM+1))**(1/(DisM+1))-1)
    Temp = (k<=ProM)*(miu>=0.5)
    Offspring[Temp] = Offspring[Temp]+(MaxValue[Temp]-MinValue[Temp])*(1-(2*(1-miu[Temp])+2*(miu[Temp]-0.5)*(1-(MaxValue[Temp]-Offspring[Temp])/(MaxValue[Temp]-MinValue[Temp]))**(DisM+1))**(1/(DisM+1)))
    Offspring[Offspring>MaxValue] = MaxValue[Offspring>MaxValue]
    Offspring[Offspring<MinValue] = MinValue[Offspring<MinValue]
    return Offspring

def F_sort(FunctionValue):
    Kind = 2
    N,M = FunctionValue.shape
    MaxFront = 0
    cz = np.zeros(N)
    FrontValue = np.zeros(N)+np.inf
    Rank = FunctionValue.argsort(axis=0)[:,0]
    FunctionValue = FunctionValue[np.lexsort(FunctionValue[:,::-1].T)]
    while (Kind==1 and np.sum(cz)<N) or (Kind==2 and np.sum(cz)<N/2) or (Kind==3 and MaxFront<1):
        MaxFront = MaxFront+1
        d = deepcopy(cz)
        for i in range(N):
            if d[i]==0:
                for j in range(i+1,N):
                    if d[j]==0:
                        k = 1
                        for m in range(1,M):
                            if FunctionValue[i,m] > FunctionValue[j,m]:
                                k = 0
                                break
                        if k == 1:
                            d[j] = 1
                FrontValue[Rank[i]] = MaxFront
                cz[i] = 1
    return FrontValue,MaxFront

def F_choose(FunctionValue1,FunctionValue2,K,Z):
    FunctionValue = np.concatenate([FunctionValue1,FunctionValue2])
    N,M = FunctionValue.shape
    N1 = FunctionValue1.shape[0]
    N2 = FunctionValue2.shape[0]
    NoZ = Z.shape[0]
    Zmin = np.min(FunctionValue,axis=0) 
    Extreme = np.zeros(M)                 
    w = np.zeros((M,M))+0.000001+np.eye(M)
    for i in range(M):                  
        Extreme[i] = np.argmin(np.max(FunctionValue/np.tile(w[i,:],(N,1)),axis=1))
    Hyperplane = np.linalg.lstsq(FunctionValue[Extreme.astype('int'),:],np.ones((M,1)))[0]
    a = 1/Hyperplane       
    if np.any(np.isnan(a)):
        a = np.max(FunctionValue,axis=0).T
    FunctionValue = (FunctionValue-np.tile(Zmin,(N,1)))/(np.tile(a.T,(N,1))-np.tile(Zmin,(N,1)))
    Distance = np.zeros((N,NoZ))   
    normZ = np.sum(Z**2,axis=1)**0.5
    normF = np.sum(FunctionValue**2,axis=1)**0.5
    for i in range(N):
        normFZ = np.sum((np.tile(FunctionValue[i,:],(NoZ,1))-Z)**2,axis=1)**0.5
        for j in range(NoZ):
            S1 = normF[i]
            S2 = normZ[j]
            S3 = normFZ[j]
            p = (S1+S2+S3)/2
            Distance[i,j] = 2*np.sqrt(p*(p-S1)*(p-S2)*(p-S3))/S2
    d = np.min(Distance.T,axis=0)
    pii = np.argmin(Distance.T,axis=0) 
    rho = np.zeros(NoZ)
    for i in range(N1):
        rho[pii[i]] = rho[pii[i]]+1
    Choose = np.zeros(N2)      
    Zchoose = np.ones(NoZ)    
    k = 1
    while k <= K:
        Temp = np.argwhere(Zchoose==1).ravel()
        j = np.argmin(rho[Temp])
        j = Temp[j]
        I1 = (np.ravel(Choose)==0).nonzero()
        I2 = (np.ravel(pii[N1:])==j).nonzero()
        I = np.intersect1d(I1,I2)
        if len(I)>0:
            if rho[j] == 0:
                s = np.argmin(d[N1+I])
            else:
                s = np.random.randint(len(I))
            Choose[I[s]] = 1
            rho[j] = rho[j]+1
            k = k+1
        else:
            Zchoose[j] = 0
    return Choose

def main(Problem,M,Generations):
    assert M==2 or M==3
    p1 = [99,13,7,5,4,3,3,2,3]
    p2 = [0,0,0,0,1,2,2,2,2]
    p1 = p1[M-2]
    p2 = p2[M-2]
    N,Z = F_weight(p1,p2,M)
    Z[Z==0] = 0.000001
    Population,Boundary = F_objective('init',Problem,M,N)
    plt.ion()
    F = None
    if M==3:
        fig = plt.figure(figsize=plt.figaspect(0.5))
    for Gene in range(Generations): 
        MatingPool = F_mating(Population)
        Offspring = F_generator(MatingPool,Boundary,N)
        Population = np.concatenate([Population,Offspring])
        FunctionValue = F_objective('value',Problem,M,Population)[0]
        FrontValue,MaxFront = F_sort(deepcopy(FunctionValue))
        Next = np.zeros(N)
        NoN = np.count_nonzero(FrontValue[FrontValue<MaxFront])
        Next[:NoN] = (np.ravel(FrontValue.T)<MaxFront).nonzero()[0]
        Last = (np.ravel(FrontValue.T)==MaxFront).nonzero()[0]
        c = Next[:NoN].astype('int')
        Choose = F_choose(FunctionValue[c,:],FunctionValue[Last,:],N-NoN,Z)
        Next[NoN:N] = Last[Choose.astype(bool)]
        Population = Population[Next.astype('int'),:]
        F = F_objective('value',Problem,M,Population)[0]
        plt.cla()
        if M==3:
            fig.suptitle(Problem+"   Generation="+str(Gene+1))
            ax = fig.add_subplot(1,2,1, projection='3d')
            fig.delaxes(fig.axes[0])
            ax.scatter(F[:, 0], F[:, 1],F[:,2], color='b')
            ax2 = fig.add_subplot(1,2,2, projection='3d')
            ax2.view_init(10,50)
            ax2.scatter(F[:, 0], F[:, 1],F[:,2], color='b')
            if Gene==Generations-1:
                plt.ioff()
                plt.show()
            else:
                plt.pause(0.2)
            
        else:
            plt.title(Problem+"   Generation="+str(Gene+1))
            plt.scatter(F[:,0], F[:,1], color='b')
            if Gene==Generations-1:
                plt.ioff()
                plt.scatter(F[:,0], F[:,1], color='b')
                plt.show()
            else:
                plt.pause(0.2)


if __name__ == "__main__":
    main(Problem="DTLZ7", M=3, Generations=500)