Multi objective optimization series 1 --- explanation of non dominated sorting function of NSGA2

As a non mathematical 、 Wild research monk programmers who are not computer majors , When learning and practicing multi-objective optimization , Encountered various difficulties , In addition, the communication resources in relevant directions are limited , It also makes the whole process seem slow . Here it is , Thank you very much XD Xiaofeng （https://blog.csdn.net/qq_40434430/article/details/82876572） And other great gods who contributed knowledge from all walks of life 、 communication , His preciseness 、 modesty 、 Enthusiasm enables me not to give up , Into the door .

Why should I explain it specifically NSGA2 Non dominated sort function

The classic of multi-objective optimization is NSGA2,NSGA2 The classic is non-dominated sort. Non dominated sorting is done , To talk about everything behind （ The comprehension of this sentence , Not born to know , But in all kinds of trial and error 、 I learned it in the process of running away ）. Mastered NSGA2 Then we can talk about other kinds of algorithms , Such as ：NSGA3,ε-NSGA2,MOEA,MOEA/D（ The formulation here is not completely correct , Their direct part has no strong coupling relationship , Just for one's own opinion ）.

The basic principle of non dominated sorting function

First, let's look at the mathematical definition ：

Several potential connotations need to be explained ：

• It's all about minimizing （ A conventional rule ）;
• be-all $f_i(X_a)$ Less than or equal to $f_i(X_b)$ , And there is at least one $f_i(X_a)$ Completely less than $f_i(X_b)$
• $f_i(X_a)$ It's exactly the same as $f_i(X_b)$ when ,$X_a$ and $X_b$ There is no dominant relationship （ Focus on ： This is in the actual evolutionary process , There will be a lot of , At present, many codes on the internet ignore this , Will result in Deb The great God is primitive NSGA2 There is no problem running the test , But there will be many problems when we turn to real problems , Such as ： Nonconvergence 、 Instability, etc . Why does this happen ？ because Deb Great God uses one of the standard genetic operators , namely ： Real number coding and polynomial crossing , And defined in the paper PC It's big , There is little possibility of complete replication between generations . However , Application and transformation in reality , Many times, we need to modify the genetic operator , namely ： Nonstandard genetic operator , Thus, all kinds of wonderful problems will appear . The reason why this problem is emphasized here , Because I spent a lot of energy here , All kinds of debugging , No reason found , I finally found out that it was because an equal number pit . Special thanks to the research group Jing Huang Give patience to debug and discuss .）
• 

How to test non dominated sorting function

We need to build test classes specifically for non dominated sorting functions , So as to verify the correctness of the function , The key to building a test class is the input test evaluation value Array The effectiveness of the . I passed a lot of experiments , Refine the basic 、 Correct 、 Three test evaluation values of errors Array.

• Basic test evaluation value Array, You can use it to deduce , This value refers to teacher Zheng Jinhua's book （ Subjective evaluation ： Mr. Zheng's book is a book that has made a systematic introduction to this part in China at present , His scholarly style is rigorous 、 Pragmatic , Instead of directly stacking papers , It is suggested that introductory students can start with this book , Books piled with papers can be learned at a higher level ）
• Correct test evaluation value Array, run ZDT3 From memory Debug Copy Of .
• The test evaluation value of the error report Array, run ZDT6 From memory Debug Copy Of . In the code given below test_fast_non_dominated_sort_error The function will show where he is wrong .

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
        #  Test code 
        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):
        # Modification and adjustment of non dominated sorting written by Huadian xiaozaza 
        # https: // blog.csdn.net / qq_36449201 / article / details / 81046586
        fronts = []  # Pareto The frontier 
        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 control  i,np+1
            set_sp.append(temp)  # i control  j, take  j  Join in  i  In the dominating solution set of 
            if npp[i] == 0:
                fronts[0].append(i)  #  Individual serial number 
                rank[i] = 1  # Pareto The frontier   The first level 
        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  #  The second level 
                        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):
        # Modification and adjustment of non dominated sorting written by Huadian xiaozaza 
        # https: // blog.csdn.net / qq_36449201 / article / details / 81046586
        fronts = []  # Pareto The frontier 
        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 control  i,np+1
            set_sp.append(temp)  # i control  j, take  j  Join in  i  In the dominating solution set of 
            if npp[i] == 0:
                fronts[0].append(i)  #  Individual serial number 
                rank[i] = 1  # Pareto The frontier   The first level 
        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  #  The second level 
                        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):
        # Yes Github Haris Ali Khan Written NSGA2 Adjustment of non dominated sorting function 
        # 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]):
                        #  individual p The dominating set of Sp Calculation 
                        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]):
                        #  Dominated degree Np Calculation 
                        # Np The bigger it is , shows i The worse the individual 
                        npp[i] += 1  # j control  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):
        #  Reference resources MoeaPlat Non dominated sorting 
        # https://blog.csdn.net/qq_40434430/article/details/82876572
        fronts = []  # Pareto The frontier 
        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' The objective function value of is less than the number of individual objective function values 
                    equal = 0  # y' The objective function value of is equal to the number of individual objective function values 
                    greater = 0  # y' The objective function value of is greater than the number of individual objective function values 
                    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 control  i,np+1
                    elif (greater == 0 and equal != self.f_num):
                        temp.append(j)
            set_sp.append(temp)  # i control  j, take  j  Join in  i  In the dominating solution set of 
            if npp[i] == 0:
                fronts[0].append(i)  #  Individual serial number 
                rank[i] = 1  # Pareto The frontier   The first level 
        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  #  The second level 
                        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  Ordinary 
    print(" test 1： Ordinary ")
    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  correct 
    print(" test 2： correct , Focus on the wrong function 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  error 
    print(" test 3：ZDT6  error , Focus on the wrong function 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)

summary

The test function written here , Objectively speaking, it is not good enough , Because it is an exhaustive way to verify that it is less than 、 be equal to 、 Greater than , Two objective functions , It's OK to write like this , But if it is 5 Goals , This is obviously not ok . So it adds test_fast_non_dominated_sort_3() function , The correct approach should be to refer to mooplatm in , Excluded writing , It can be applied to multiple targets , Not just limited to two goals .