1.介绍

2. NSGA-II

2.1 快速非支配排序

2.1.1 NSGA的传统非支配排序

2.1.2 NSGA-II的快速非支配排序

2.2 多样性保护(Diversity Preservation)

2.2.2 NSGA-II的拥挤距离方法(crowded-comparison)

2.3 NSGA-II主循环

3. 代码实现

3.1 第三方库

3.2 自定义

NSGA提出原文 Muiltiobjective Optimization Using Nondominated Sorting in Genetic Algorithms | MIT Press Journals & Magazine | IEEE Xplore

Srinivas, N., & Deb, K. (1994). Muiltiobjective Optimization Using Nondominated Sorting in Genetic Algorithms. Evolutionary Computation, 2(3), 221-248. doi:10.1162/evco.1994.2.3.221

NSGA-II提出原文 A fast and elitist multiobjective genetic algorithm: NSGA-II

Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation, 6(2), 182-197.

NSGA-II作者实验室：Kalyanmoy Deb, Koenig Endowed Chair Professor

pymoo库：实验室创建的python第三方库，实现了各种多目标优化算法pymoo: Multi-objective Optimization in Python

geatpy2库：Geatpy是一个高性能实用型进化算法工具箱，提供许多已实现的进化算法中各项重要操作的库函数，利用“定义问题类 + 调用算法模板”的模式来进行进化优化，可用于求解单目标优化、多目标优化、复杂约束优化、组合优化、混合编码进化优化等 Geatpy

1.介绍

针对多目标优化问题，可以用一些多目标进化算法(multiobjective evolutionary algorithms (MOEAs))找到多个帕累托最优解(Pareto-optimal)，比如非支配排序基因算法(nondominated sorting genetic algorithm (NSGA))。但是NSGA有以下问题

非支配排序时间复杂度太高，为 $O(MN^{3})$ ，其中M为多目标数，N为种群数
缺少精英保留策略(elitism)，研究表明精英策略能提高GA的性能
需要指定共享参数(share parameter)来确保种群多样性

基于上述背景，文章提出NSGA改进版本——NSGA-II，NSGA-II有三个重点的改进：快速非支配排序(fast nondominated sorting)，精英保留策略(elitist-preserving)，无参数的小生境算子(parameterless niching operator)。文章内容如下：

重点介绍NSGA-II的改进和流程
将NSGA-II和其它两种带精英策略的MOEAs比较，Pareto-archived evolution strategy (PAES)和strength- Pareto EA (SPEA)，证明NSGA-II在解质量和收敛性表现更好
强调参数相互作用(parameter interactions)的问题
扩展NSGA-II处理带约束问题(constraint-handling strategy with NSGA-II)

2. NSGA-II

2.1 快速非支配排序

先明白支配解，非支配解，帕累托前沿(Pareto front)的概念，假设任何二解S1 及S2 对所有目标而言，S1均优于S2，则我们称S1 支配(dominate)S2，若S1 的解没有被其他解所支配，则S1 称为非支配解（不受支配解nondominated），也称Pareto解。这些非支配解的集合即所谓的帕累托前沿(Pareto front)。非支配排序就是将种群分成多个不同水平的帕累托前沿。

2.1.1 NSGA的传统非支配排序

在传统方法中，为了在大小为N的种群中获得第一个非支配前沿，对于每一个解，都需要将其和种群中的其它解比较所有目标值，确定该解是否被支配，因此对每个解时间复杂度为 $O(MN)$ ，那么遍历完所有解找到种群中的第一层非支配集合，时间复杂度为 $O(MN^{2})$ 。在找到第一个非支配前沿后，这些解就会移除，在剩下的解中继续重复上述操作，找到接下来所有层的非支配前沿。这种排序方法最坏的情况是，有N个前沿，每个前沿只有一个解，意味着上述操作重复了N词，总时间复杂度为 $O(MN^{3})$ ，空间复杂度为 $O(N)$ 。

2.1.2 NSGA-II的快速非支配排序

NSGA-II中提出的的快速非支配排序将总时间复杂度降到了 $O(MN^{2})$ ，空间复杂度为 $O(N^{2})$ 。

首先，遍历种群中的所有解，对于每个解，计算两个变量值

支配数 $n_p$ ：支配解p的解数量
被支配集合 $S_{p}$ ：解p支配的解集合

如果 $n_p=0$ ，那么解p在第一层非支配前沿，通过上述操作，找到每个解的支配数和被支配集合以及第一层非支配前沿，时间复杂度为 $O(MN^{2})$ 。

接下来，对于每个支配数为0的的解p，遍历其p的支配集合中的每个成员q，q的支配数减少1，如果出现任何q支配数变成了0，则将其放入单独的列表Q中，Q中的成员便是第二层非支配前沿，对于Q中的每个成员重复上述操作找到第三层。重复以上过程直到所有的解已经排序结束。这样排序实际上是两层循环，外循环遍历最新的非支配前沿，内循环遍历非支配前沿中每个成员的支配集合。由于每个解只能在一层非支配前沿里，因此外循环最多N次，而支配集合最多有N-1个成员，内循环最多N-1次，因此时间复杂度为 $O(N^{2})$ 。也可以换一种思路理解，支配数最大为N-1，即对于每个解，到其支配数减少0最多遍历到N-1次，所有解的支配数减到0需要 $O(N^{2})$ 。下面为快速非支配排序的算法流程。

2.2 多样性保护(Diversity Preservation)

在NSGA中，与普通的GA相比，交叉和变异操作相同，区别在于选择操作。首先按照传统非支配排序的方法获得种群的第一非支配前沿(first nondominated front)，给第一非支配前沿中的每个个体分配一个相同的大的虚拟适应值(dummy fitness value)，也就是最初第一非支配前沿的每个个体选择的可能性相同。接下来为了保留每层前沿解的多样性，使用共享函数方法(sharing function)得到每个解的共享适应值(shared fitness value)。共享函数如下，其中 $d_{ij}$ 是当前前沿中任何两个个体的表型距离(phenotypic distance)， $\sigma _{share}$ 为需要设置的共享参数。定义一个变量niche count，niche count为当前前沿所有个体的共享函数值和，最后对于每个解，其共享适应值=虚拟适应值 / niche count。个人通俗理解，同一层的解之间会有相似度，对于每个解，如果与它相似的解很多，那么就降低这个解的适应值，减小它被选择的概率，而同一个生态圈的解数越多，其选择概率越小，但是保留这个生态圈里的解的整体期望还是不变，也维持了解的多样性

2.2.2 NSGA-II的拥挤距离方法(crowded-comparison)

共享函数方法有两个问题

用共享函数维持解的分布的性能很大程度取决于参数 $\sigma _{share}$
每个解都需要和其它所有解比较，共享函数的时间复杂度为 $O(N^{2})$

NSGA-II用拥挤度比较的方法代替共享函数来保持种群的多样性，不仅无参数，还减小了时间复杂度，首先定义一个密度估计度量(density-estimation metric)，然后描述拥挤度比较算法

（1）密度估计(density-estimation)

为了估计种群中某个解被其他解包围的密度，引入拥挤距离(crowding distance) $i_{distance}$ 这个概念。计算拥挤距离，先对每一个等级的非支配解集合，按照每个目标函数值的大小升序排列。因此对每个目标函数，都能找到边界解(目标值最大最小的解)，设置边界解的距离值为无穷大。每个解的拥挤距离为每个目标函数下的距离值之和，求和前先对距离值进行归一化处理，算法流程如下： $\pounds [i].m$ 代表非支配前沿 $\pounds$ 中第i个解的第m个目标函数值， $f_{m}^{max}$ 和 $f_{m}^{min}$ 为第m个目标函数的最大值和最小值。分析时间复杂度：首先排序算法需要logn，外层循环M次，内层循环最多N次，因此总的时间复杂度为 $O(MNlogN)$

从几何角度理解拥挤距离，以两个目标函数为例，下图中黑点和白点分别为两个非支配前沿，对于解i，从与i在同一非支配前沿中选择与解i最相近的两个点i-1和i+1为顶点组成一个长方形(cuboid)，拥挤距离即为长方形周长。

（2）拥挤度比较算子(Crowded-Comparison Operator)

在计算完每个解的拥挤距离后，从一定程度上来说，某个解的拥挤距离越小，这个解被其他解拥挤的程度越高。接下来通过拥挤度比较算子来选择解实现更广的帕累托最优解分布。种群中的每个个体都有两个属性：

非支配等级 $i_{rank}$ ：1是最高等级
拥挤距离 $i_{distance}$

定义一个比较顺序 $\prec _{n}$ ，满足一下判断条件。对于不同非支配等级的两个解，倾向于选择rank值更低的解，如果两个解的等级相同，倾向于选择拥挤距离更大或者说拥挤区域更小的解。

综上，有三个创新点，一个快速非支配排序方法，一个快速估计拥挤距离的方法，一个简单的拥挤度比较算子，接下来描述NSGA-II的主循环。

2.3 NSGA-II主循环

先明确以下概念

elitism：精英保留策略，核心思想是把群体在进化过程中迄今出现的最好个体不进行遗传操作而直接复制到下一代中，理论上证明了具有精英保留的标准遗传算法是全局收敛的
tournament selection：联赛选择算法，每次从种群中取一定数量（n）的个体(放回抽样)，选择其中适应度较好的进入子代种群

随机生成初始种群P(0)，基于非支配排序P(0)，每个解都赋予一个适应值与其前沿等级相同(1是最好的，2其次....），假设问题是要最小化。起初使用联赛选择，重组和变异获得一个子代Q(0)，NSGA-II使用了精英选择策略，因此后面迭代过程会不同。假设现在是第t代种群，下面一步步描述结合了精英选择，快速非支配排序，拥挤距离等方法的NSGA-II算法过程。

首先生成一个基于P(t)和Q(t)得到的组合种群R(t)，R(t)大小为2N，N为种群大小。对R(t)执行快速非支配排序，属于R(t)的第一非支配前沿F1中的解是当前最好的解，因此要重点选择。如果F1长度小于N，那么选择F1添加到新一代种群P(t+1)中，同样操作在F2,F3....，直到P(t+1)剩余解数量不足以再合并一个完整的非支配前沿。假设F(i)是最后一个无法容纳的非支配前沿，按照拥挤距离给F(i)中的元素降序排列，选择排序后靠前的元素添加到P(t+1)中直到P(t+1)大小为N。对P(t+1)执行选择，交叉，变异生成新种群Q(t+1)，此处的选择虽然也是联赛选择方法，但是选择标准是基于拥挤度比较算子。这个算子需要每个解的等级和拥挤距离，因此在组成种群P(t+1)的时候就可以计算这两个属性。P(t+1)是从大小为2N的R(t)中按照前沿等级和拥挤距离生成的大小为N的种群，而R(t)包含了所有之前和现在的种群成员，确保了精英选择策略。算法流程如下图

分析整个过程的时间复杂度如下， NSGA-II的整体时间复杂度为 $O(MN^{2})$ 。实际上在对R(t)做快速非支配排序时，并不需要多整个种群排序，只需找到N个解即可。使用拥挤距离选择下一代也避免了人工参数的设置。

快速非支配排序： $O(MN^{2})$
拥挤距离计算： $O(MNlogN)$
$\prec _{n}$ 的排序： $O(NlogN)$

3. 代码实现

3.1 第三方库

pymoo库了提供最先进的单目标和多目标优化算法，以及与多目标优化相关的更多功能，例如可视化和决策制定。比如用NSGA-II解决一个二元决策问题，示例代码如下

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.problems import get_problem
from pymoo.operators.crossover.pntx import TwoPointCrossover
from pymoo.operators.mutation.bitflip import BitflipMutation
from pymoo.operators.sampling.rnd import BinaryRandomSampling
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

problem = get_problem("zdt5")

algorithm = NSGA2(pop_size=100,
                  sampling=BinaryRandomSampling(),
                  crossover=TwoPointCrossover(),
                  mutation=BitflipMutation(),
                  eliminate_duplicates=True)

res = minimize(problem,
               algorithm,
               ('n_gen', 500),
               seed=1,
               verbose=False)

Scatter().add(res.F).show()

和pymmo相比，geatpy2支持自定义问题，提供模块化的算法框架，示例代码如下

import geatpy as ea
import numpy as np

# 构建问题
r = 1  # 目标函数需要用到的额外数据
@ea.Problem.single
def evalVars(Vars):  # 定义目标函数（含约束）
    f = np.sum((Vars - r) ** 2)  # 计算目标函数值
    x1 = Vars[0]
    x2 = Vars[1]
    CV = np.array([(x1 - 0.5)**2 - 0.25,
                    (x2 - 1)**2 - 1])  # 计算违反约束程度
    return f, CV

problem = ea.Problem(name='soea quick start demo',
                        M=1,  # 目标维数
                        maxormins=[1],  # 目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标
                        Dim=5,  # 决策变量维数
                        varTypes=[0, 0, 1, 1, 1],  # 决策变量的类型列表，0：实数；1：整数
                        lb=[-1, 1, 2, 1, 0],  # 决策变量下界
                        ub=[1, 4, 5, 2, 1],  # 决策变量上界
                        evalVars=evalVars)
# 构建算法
algorithm = ea.soea_SEGA_templet(problem,
                                    ea.Population(Encoding='RI', NIND=20),
                                    MAXGEN=50,  # 最大进化代数。
                                    logTras=1,  # 表示每隔多少代记录一次日志信息，0表示不记录。
                                    trappedValue=1e-6,  # 单目标优化陷入停滞的判断阈值。
                                    maxTrappedCount=10)  # 进化停滞计数器最大上限值。
# 求解
res = ea.optimize(algorithm, seed=1, verbose=True, drawing=1, outputMsg=True, drawLog=False, saveFlag=Tr

关于第三方库的内容可自行查阅官方网站

3.2 自定义

按照前面的算法流程分别实现快速非支配排序模块和拥挤距离计算模块

（1）快速非支配排序模块

"""
两个目标函数为例
输入：种群每个解的两个目标函数值
输入：所有非支配前沿
"""
def fast_non_dominated_sort(values1, values2):
    size = len(values1)  # 种群大小
    s = [[] for _ in range(size)]  # 每个解的被支配集合
    n = [0 for _ in range(size)]  # 每个解的支配数
    rank = [0 for _ in range(size)]  # 每个解的等级
    fronts = [[]]  # 所有非支配前沿

    for p in range(size):  # 遍历所有解
        s[p] = []  # 初始化非支配集合和支配数
        n[p] = 0
        for q in range(size):  # 判断p和q支配情况
            # 如果p支配q，增加q到p的被支配集合
            if values1[p] >= values1[q] and values2[p] >= values2[q] \
                    and ((values1[q] == values1[p]) + (values2[p] == values2[q])) != 2:
                s[p].append(q)
            # 如果q支配p，p的支配数+1
            elif values1[q] >= values1[p] and values2[q] >= values2[p] \
                    and ((values1[q] == values1[p]) + (values2[p] == values2[q])) != 2:
                n[p] += 1
        # n[p]=0的解等级设为0，增加到第一前沿
        if n[p] == 0:
            rank[p] = 0
            fronts[0].append(p)
    # 依次确定其它层非支配前沿
    i = 0
    while fronts[i]:
        Q = []
        for p in fronts[i]:
            for q in s[p]:
                n[q] = n[q] - 1
                if n[q] == 0:
                    rank[q] = i + 1
                    if q not in Q:
                        Q.append(q)
        i = i + 1
        fronts.append(Q)

    del fronts[-1]
    return fronts

（2）拥挤距离计算模块

"""
输入：所有解所有目标的函数值，要计算的前沿，目标数
输出：输入前沿每个解的拥挤距离
"""
def crowed_distance_assignment(total_values, front, objective):
    length = len(front)
    distance = [0 for i in range(length)]
    for i in range(objective):
        value = total_values[i]
        sorted_front = sorted(front, key=lambda x: x in value)
        distance[0] = distance[-1] = np.inf
        for k in range(1, length - 1):
            distance[k] = distance[k] + (value[sorted_front[k+1]] - value[sorted_front[k-1]]) / (max(value) - min(value))
    return distance

完整代码待更新