多目标优化系列1---NSGA2的非支配排序函数的讲解

作为一名非数学、非计算机专业的野生研究僧程序员，在学习和实践多目标优化时，遇到了各种困难，加之相关方向的交流资源有限，也使得整个过程显得缓慢。在此，非常感谢西电晓风（https://blog.csdn.net/qq_40434430/article/details/82876572）及其他各路贡献知识的大神给予的启发、交流，他的严谨、谦虚、热情让我能不抛弃不放弃，入了门。

为何要专门讲解NSGA2非支配排序函数

多目标优化中的经典是NSGA2，NSGA2的经典是non-dominated sort。非支配排序做好了，才能谈后面的一切（这句话的领悟，非生而知之，而是在各种试错、奔溃的过程中领略到的）。掌握了NSGA2才能谈后面的其他各类算法，如：NSGA3，ε-NSGA2，MOEA，MOEA/D（这里的提法也并非完全正确，他们直接部分并没有强耦合关系，仅为个人一己之见）。

非支配排序函数的基本原理

首先看一下数学定义：

需要说明几点潜在内涵：

• 都是求最小化问题（一个约定俗成的规则）；
• 所有的 $f_i(X_a)$ 小于或者等于$f_i(X_b)$ ，且至少存在一个 $f_i(X_a)$ 完全小于$f_i(X_b)$
• $f_i(X_a)$ 完全等于$f_i(X_b)$时，$X_a$和$X_b$不存在支配关系（划重点：这个在实际的进化过程中，会出现很多，目前网上很多的代码忽略了这一点，会导致在已有Deb大神原始NSGA2测试问题上运行没问题，但是转换到现实问题中会出现很多的问题，如：不收敛、不稳定等。为什么会出现这样的问题呢？因为Deb大神采用的是标准遗传算子中的一种，即：实数编码和多项式交叉，且论文中定义的PC很大，世代之间出现完全复制的可能性小。然而，现实中应用改造问题，很多时候需要改造掉遗传算子，即：非标遗传算子，从而就会展现出各种奇葩问题。之所以在这里强调这问题，是因为我自己在这里耗费了很多的精力，各种调试，找不出原因，最终发现是因为一个等于号坑了我。也特别感谢课题组Jing Huang给予的耐心调试和讨论。）
• 

如何测试非支配排序函数

我们需要构建专门针对非支配排序函数的测试类，从而验证函数的正确性，构建测试类的关键是在于输入的测试评价值Array的有效性。我通过大量的实验，提炼出基础的、正确的、报错的三个测试评价值Array。

• 基础的测试评价值Array，大家可以用来自行推演，该值参考至郑金华老师的书（主观评价：郑老师的书是目前国内对此块做了系统介绍的书，他的治学风格严谨、务实，而非论文的直接堆砌，建议入门学生可以从这本书开始，论文堆砌的书可以到高阶阶段再学习）
• 正确的测试评价值Array，跑ZDT3的时候从内存中Debug Copy的。
• 报错的测试评价值Array，跑ZDT6的时候从内存中Debug Copy的。以下给出的代码中test_fast_non_dominated_sort_error函数会显示出他错在哪。

import random
import numpy as np
from matplotlib.ticker import MultipleLocator
import matplotlib.pyplot as plt


class Test_class():

    def __init__(self,f_num):
        self.f_num=f_num
        # 测试代码
        Y1 = [9, 7, 5, 4, 3, 2, 1, 10, 8, 7, 5, 4, 3, 10, 9, 8, 7, 10, 9, 8]
        Y2 = [1, 2, 4, 5, 6, 7, 9, 1,  5, 6, 7, 8, 9,  5, 6, 7, 9,  6, 7, 9]
        self.objectives_fitness_zjh=np.array([Y1,Y2]).T
        self.objectives_fitness_8 = np.array([[0.6373723585181119, 9.089424920752537],
                                   [0.6307563745957109, 9.484134522661321],
                                   [0.6307726564027054, 9.370453232401315],
                                   [0.9214017573731662, 8.974173267351892],
                                   [0.8106208092655269, 9.00814519794432],
                                   [0.6308299859236132, 9.381311843663337],
                                   [0.9996933421004693, 8.998732212375378],
                                   [0.8106208092655269, 9.087298161567794],
                                   [0.9968206553777186, 9.020037858133483],
                                   [0.9503113437004427, 9.04519027298749],
                                   [0.9214017573731662, 8.974173267351892],
                                   [0.9214017573731662, 9.00187266065025],
                                   [0.6307563745957109, 9.493829448943414],
                                   [0.999999999999489, 9.180616877016963],
                                   [0.8150090640161689, 9.041622216379706],
                                   [0.9990805452551389, 8.980910429010232],
                                   [0.9979468094812165, 9.04561922020985],
                                   [0.7527617476769539, 9.136335451211739],
                                   [0.6364241984468356, 9.20992553291975],
                                   [0.8106208092655269, 9.083868693443087]])

        self.objectives_fitness_20 = np.array([[0.01783689,1.02469787],
                                    [0.04471213,0.93037726],
                                    [0.0236877,0.96322404],
                                    [0.04938809,0.92641409],
                                    [0.07031985,0.83355839],
                                    [0.05199109,0.88398541],
                                    [0.05006115,0.91107188],
                                    [0.,1.18579768],
                                    [0.05358649,0.85849188],
                                    [0.01657041,1.0721905,],
                                    [0.10409921,0.84829613],
                                    [0.06643826,0.84941993],
                                    [0.05358649,0.89876683],
                                    [0.01740193,1.09732744],
                                    [0.04938809,0.94687415],
                                    [0.05199109,0.93536007],
                                    [0.02400905,1.11603965],
                                    [0.05358649,0.89107165],
                                    [0.,1.23996387],
                                    [0.01783689,1.11625582]])

    def test_fast_non_dominated_sort_1(self, objectives_fitness):
        #对华电小扎扎写的非支配排序的修改和调整
        # https: // blog.csdn.net / qq_36449201 / article / details / 81046586
        fronts = []  # Pareto前沿面
        fronts.append([])
        set_sp = []
        npp = np.zeros(2*10)
        rank = np.zeros(2*10)
        for i in range(2 * 10):
            temp = []
            for j in range(2 * 10):
                if j != i:
                    if (objectives_fitness[j][0] >= objectives_fitness[i][0] and objectives_fitness[j][1] > objectives_fitness[i][1]) or \
                        (objectives_fitness[j][0] > objectives_fitness[i][0] and objectives_fitness[j][1] >= objectives_fitness[i][1]) or \
                        (objectives_fitness[j][0] >= objectives_fitness[i][0] and objectives_fitness[j][1] >= objectives_fitness[i][1]):
                        temp.append(j)
                    elif (objectives_fitness[i][0] >= objectives_fitness[j][0] and objectives_fitness[i][1] > objectives_fitness[j][1]) or \
                        (objectives_fitness[i][0] > objectives_fitness[j][0] and objectives_fitness[i][1] >= objectives_fitness[j][1]) or \
                        (objectives_fitness[j][0] > objectives_fitness[i][0] and objectives_fitness[j][1] > objectives_fitness[i][1]):
                        npp[i] += 1  # j支配 i，np+1
            set_sp.append(temp)  # i支配 j，将 j 加入 i 的支配解集里
            if npp[i] == 0:
                fronts[0].append(i)  # 个体序号
                rank[i] = 1  # Pareto前沿面 第一层级
        i = 0
        while len(fronts[i]) > 0:
            temp = []
            for j in range(len(fronts[i])):
                a = 0
                while a < len(set_sp[fronts[i][j]]):
                    npp[set_sp[fronts[i][j]][a]] -= 1
                    if npp[set_sp[fronts[i][j]][a]] == 0:
                        rank[set_sp[fronts[i][j]][a]] = i + 2  # 第二层级
                        temp.append(set_sp[fronts[i][j]][a])
                    a = a + 1
            i = i + 1
            fronts.append(temp)
        del fronts[len(fronts) - 1]
        self.output_fronts(fronts)

    def output_fronts(self, fronts):
        # test code
        sum_coun = 0
        for kk in range(len(fronts)):
            sum_coun += len(fronts[kk])
        print(sum_coun)
        print(fronts)

    def test_fast_non_dominated_sort_error(self, objectives_fitness):
        #对华电小扎扎写的非支配排序的修改和调整
        # https: // blog.csdn.net / qq_36449201 / article / details / 81046586
        fronts = []  # Pareto前沿面
        fronts.append([])
        set_sp = []
        npp = np.zeros(2*10)
        rank = np.zeros(2*10)
        for i in range(2 * 10):
            temp = []
            for j in range(2 * 10):
                if j != i:
                    if (objectives_fitness[j][0] >= objectives_fitness[i][0] and objectives_fitness[j][1] >= objectives_fitness[i][1]) :
                        temp.append(j)
                    elif (objectives_fitness[i][0] > objectives_fitness[j][0] and objectives_fitness[i][1] > objectives_fitness[j][1]):
                        npp[i] += 1  # j支配 i，np+1
            set_sp.append(temp)  # i支配 j，将 j 加入 i 的支配解集里
            if npp[i] == 0:
                fronts[0].append(i)  # 个体序号
                rank[i] = 1  # Pareto前沿面 第一层级
        i = 0
        while len(fronts[i]) > 0:
            temp = []
            for j in range(len(fronts[i])):
                a = 0
                while a < len(set_sp[fronts[i][j]]):
                    npp[set_sp[fronts[i][j]][a]] -= 1
                    if npp[set_sp[fronts[i][j]][a]] == 0:
                        rank[set_sp[fronts[i][j]][a]] = i + 2  # 第二层级
                        temp.append(set_sp[fronts[i][j]][a])
                    a = a + 1
            i = i + 1
            fronts.append(temp)
        del fronts[len(fronts) - 1]
        self.output_fronts(fronts)

    def output_fronts(self, fronts):
        # test code
        sum_coun = 0
        for kk in range(len(fronts)):
            sum_coun += len(fronts[kk])
        print(sum_coun)
        print(fronts)

    def test_fast_non_dominated_sort_2(self, objectives_fitness):
        #对Github Haris Ali Khan写的NSGA2的非支配排序函数的调整
        # coding:utf-8
        # Program Name: NSGA-II.py
        # Description: This is a python implementation of Prof. Kalyanmoy Deb's popular NSGA-II algorithm
        # Author: Haris Ali Khan
        # Supervisor: Prof. Manoj Kumar Tiwari
        set_sp=[[] for i in range(0, np.shape(objectives_fitness)[0])]
        fronts = [[]]
        npp=[0 for i in range(0, np.shape(objectives_fitness)[0])]
        rank = [0 for i in range(0, np.shape(objectives_fitness)[0])]
        for i in range(0, np.shape(objectives_fitness)[0]):
            set_sp[i]=[]
            # npp[i]=0
            for j in range(0, np.shape(objectives_fitness)[0]):
                if i != j:
                    if (objectives_fitness[j][0] >= objectives_fitness[i][0] and objectives_fitness[j][1] > objectives_fitness[i][1]) or \
                        (objectives_fitness[j][0] > objectives_fitness[i][0] and objectives_fitness[j][1] >= objectives_fitness[i][1]) or \
                        (objectives_fitness[j][0] >= objectives_fitness[i][0] and objectives_fitness[j][1] >= objectives_fitness[i][1]):
                        # 个体p的支配集合Sp计算
                        set_sp[i].append(j)
                    elif (objectives_fitness[i][0] >= objectives_fitness[j][0] and objectives_fitness[i][1] > objectives_fitness[j][1]) or \
                        (objectives_fitness[i][0] > objectives_fitness[j][0] and objectives_fitness[i][1] >= objectives_fitness[j][1]) or \
                        (objectives_fitness[j][0] > objectives_fitness[i][0] and objectives_fitness[j][1] > objectives_fitness[i][1]):
                        # 被支配度Np计算
                        # Np越大，则说明i个体越差
                        npp[i] += 1  # j支配 i，np+1
            if npp[i]==0:
                rank[i] = 0
                if i not in fronts[0]:
                    fronts[0].append(i)
        i = 0
        while(fronts[i] != []):
            Q=[]
            for p in fronts[i]:
                for q in set_sp[p]:
                    npp[q] =npp[q] - 1
                    if( npp[q]==0):
                        rank[q]=i+1
                        if q not in Q:
                            Q.append(q)
            i = i+1
            fronts.append(Q)
        del fronts[len(fronts)-1]
        self.output_fronts(fronts)

    def test_fast_non_dominated_sort_3(self, objectives_fitness):
        # 参考MoeaPlat的非支配排序
        # https://blog.csdn.net/qq_40434430/article/details/82876572
        fronts = []  # Pareto前沿面
        fronts.append([])
        set_sp = []
        npp = np.zeros(2 * 10)
        rank = np.zeros(2 * 10)
        for i in range(2 * 10):
            temp = []
            for j in range(2 * 10):
                if j != i:
                    less = 0  # y'的目标函数值小于个体的目标函数值数目
                    equal = 0  # y'的目标函数值等于个体的目标函数值数目
                    greater = 0  # y'的目标函数值大于个体的目标函数值数目
                    for k in range(self.f_num):
                        if (objectives_fitness[i][k] < objectives_fitness[j][k]):
                            less = less + 1
                        elif (objectives_fitness[i][k] == objectives_fitness[j][k]):
                            equal = equal + 1
                        else:
                            greater = greater + 1
                    if (less == 0 and equal != self.f_num):
                        npp[i] += 1  # j支配 i，np+1
                    elif (greater == 0 and equal != self.f_num):
                        temp.append(j)
            set_sp.append(temp)  # i支配 j，将 j 加入 i 的支配解集里
            if npp[i] == 0:
                fronts[0].append(i)  # 个体序号
                rank[i] = 1  # Pareto前沿面 第一层级
        i = 0
        while len(fronts[i]) > 0:
            temp = []
            for j in range(len(fronts[i])):
                a = 0
                while a < len(set_sp[fronts[i][j]]):
                    npp[set_sp[fronts[i][j]][a]] -= 1
                    if npp[set_sp[fronts[i][j]][a]] == 0:
                        rank[set_sp[fronts[i][j]][a]] = i + 2  # 第二层级
                        temp.append(set_sp[fronts[i][j]][a])
                    a = a + 1
            i = i + 1
            fronts.append(temp)
        del fronts[len(fronts) - 1]
        self.output_fronts(fronts)


if __name__ == '__main__':

    test_class=Test_class(f_num=2)

    # test zjh 普通
    print("测试1：普通")
    test_class.test_fast_non_dominated_sort_1(test_class.objectives_fitness_zjh)
    test_class.test_fast_non_dominated_sort_2(test_class.objectives_fitness_zjh)
    test_class.test_fast_non_dominated_sort_3(test_class.objectives_fitness_zjh)
    test_class.test_fast_non_dominated_sort_error(test_class.objectives_fitness_zjh)
    # test ZDT3 正确
    print("测试2：正确，重点看错误的函数test_fast_non_dominated_sort_error")
    test_class.test_fast_non_dominated_sort_1(test_class.objectives_fitness_20)
    test_class.test_fast_non_dominated_sort_2(test_class.objectives_fitness_20)
    test_class.test_fast_non_dominated_sort_3(test_class.objectives_fitness_20)
    test_class.test_fast_non_dominated_sort_error(test_class.objectives_fitness_20)
    # test ZDT6 错误
    print("测试3：ZDT6 错误，重点看错误的函数test_fast_non_dominated_sort_error")
    test_class.test_fast_non_dominated_sort_1(test_class.objectives_fitness_8)
    test_class.test_fast_non_dominated_sort_2(test_class.objectives_fitness_8)
    test_class.test_fast_non_dominated_sort_3(test_class.objectives_fitness_8)
    test_class.test_fast_non_dominated_sort_error(test_class.objectives_fitness_8)

总结

这里写的测试函数，客观来说还是不够好的，因为是采用穷举的方式来验证小于、等于、大于的关系，两个目标函数的时候，这样写还是可以凑合的，但是如果是5个目标，这个显然就不可以了。所以增加了test_fast_non_dominated_sort_3()函数，正确的做法应该是参考mooplatm中，排除的写法，可以适用于多个目标，而不仅仅局限于两个目标。