????????在接觸學(xué)習(xí)多目標優(yōu)化的問題上，經(jīng)常會被提及到多目標遺傳算法NSGA-II，網(wǎng)上也看到了很多人對該算法的總結(jié)，但真正講解明白的以及配套用算法實現(xiàn)的文章很少，這里也對該算法進行一次詳解與總結(jié)。會有側(cè)重點的闡述，不會針對NSGA以及算法復(fù)雜度等等進行，有相關(guān)需求可自行學(xué)習(xí)了解。文末附有Python與Matlab兩種的源代碼案例。

---

在串講算法原理之前，先說明幾個關(guān)鍵詞。

1、什么是非支配遺傳算法

????????非支配排序遺傳算法NSGA (Non-dominated Sorting Genetic Algorithms)是由Srinivas和Deb于1995年提出的。這是一種基于Pareto最優(yōu)概念的遺傳算法，它是眾多的多目標優(yōu)化遺傳算法中體現(xiàn)Goldberg思想最直接的方法。該算法就是在基本遺傳算法的基礎(chǔ)上，對選擇再生方法進行改進:將每個個體按照它們的支配與非支配關(guān)系進行分層，再做選擇操作，從而使得該算法在多目標優(yōu)化方面得到非常滿意的結(jié)果。

2、支配與非支配關(guān)系

????????在最小化問題中，
理解起來就是，任意一個目標分量，的值都不大于，并且在這些目標分量中至少存在一個目標分量的值嚴格比的值小，記作.

2.1、種群分層

????????當出現(xiàn)大量滿足支配與非支配的個體時，會出現(xiàn)一種如下的分層現(xiàn)象。因為優(yōu)化問題，通常存在一個解集，這些解之間就全體目標函數(shù)而言是無法比較優(yōu)劣的，其特點是：無法在改進任何目標函數(shù)的同時不削弱至少一個其他目標函數(shù)。
????????這種解稱作非支配解（nondominated solutions）或Pareto最優(yōu)解（Pareto optimal solutions），形象展示如下：

2.2、如何進行分支配排序

（1）設(shè)i=1；

（2）所有的且，按照以上定義比較個體和個體之間的支配與非支配關(guān)系；

（3）如果不存在任何一個個體優(yōu)于，，則標記為非支配個體

（4）令i=i+1，轉(zhuǎn)到步驟（2），直到找到所有的非支配個體。

???? 通過上述步驟得到的非支配個體集是種群的第一級非支配層，然后忽略這些已經(jīng)標記的非支配個體（即這些個體不再進行下一輪比較），再遵循步驟（1）-（4），就會得到第二級非支配層。依此類推，直到整個種群被分層。

3、擁擠度（密度評估）

????????為了獲得對種群中某一特定解附近的密度估計， 我們沿著對應(yīng)目標計算這個解兩側(cè)的兩點間的平均距離。 這個值是對用最近的鄰居作為頂點形成的長方體的周長的估計（稱之為擁擠距離）。在圖 1 中， 第i 個解在其前沿的擁擠距離（用實心圓標記）是長方體的平均邊長（用虛線框表示。（用實心圓標記的點是同一非支配前沿的解）

????????擁擠距離的計算需要根據(jù)每個目標函數(shù)值的大小對種群進行升序排序。此后，對于每個目標函數(shù)，邊界解（對應(yīng)函數(shù)最大最小值的解）被分配一個無窮距離值。所有中間的其他解被分配一個與兩個相鄰解的函數(shù)值的歸一化差的絕對值相同的距離值。 對其他目標函數(shù)繼續(xù)進行這種計算?？偟膿頂D距離值為對應(yīng)于每個目標的各個距離之和。在計算擁擠距離之前， 歸一化每個目標函數(shù)。 下面所示的算法概述了非支配性解集 中所有解的擁擠距離計算過程。對非支配性解集進行擁擠距離分配

?????????這里， 是指集合中第個個體的第個目標函數(shù)值，參數(shù)? 和是第m 個目標函數(shù)的最大和最小值.

4、快速非支配性排序方法

????????對于每個個體 都設(shè)有以下兩個參數(shù) 和 ，： 在種群中支配個體 的解個體的數(shù)量。： 被個體 所支配的解個體的集合。

?
1）找到種群中所有 n(i)=0 的個體， 將它們存入當前集合F(1)；
2）對于當前集合 F(1) 中的每個個體 ， 考察它所支配的個體集 ， 將集合 中的每個個體 k 的 n(k) 減去1， 即支配個體 k 的解個體數(shù)減1；（ 對其他解除去被第一層支配的數(shù)量， 即減
1） 如 n(k)-1=0則將個體 k 存入另一個集H。
（ 3） 將 F(1) 作為第一級非支配個體集合， 并賦予該集合內(nèi)個體一個相同的非支配序 i(rank)， 然后繼續(xù)對 H 作上述分級操作并賦予相應(yīng)的非支配序， 直到所有的個體都被分級。 其計算復(fù)雜度為，m為目標函數(shù)個數(shù)， N為種群大小。按照1 )、 2） 的方法完成所有分級,示意圖如下。

?

5、偏序關(guān)系

????????由于經(jīng)過了排序和擁擠度的計算，群體中每個個體都得到兩個屬性：非支配序和擁擠度，則定義偏序關(guān)系為：當滿足條件，或滿足且時，定義。也就是說：如果兩個個體的非支配排序不同，取排序號較小的個體（分層排序時，先被分離出來的個體）；如果兩個個體在同一級，取周圍較不擁擠的個體。

6、 什么是精英策略

????????NSGA-II 算法采用精英策略防止優(yōu)秀個體的流失， 通過將父代和子代所有個體混合后進行
非支配排序的方法。 精英策略的執(zhí)行步驟：

7、NSGA-II算法的基本思想

1） 隨機產(chǎn)生規(guī)模為N的初始種群Pt，經(jīng)過非支配排序、 選擇、 交叉和變異， 產(chǎn)生子
代種群Qt， 并將兩個種群聯(lián)合在一起形成大小為2N的種群Rt,；
2）進行快速非支配排序， 同時對每個非支配層中的個體進行擁擠度計算， 根據(jù)非支
配關(guān)系以及個體的擁擠度選取合適的個體組成新的父代種群Pt+1﹔
3） 通過遺傳算法的基本操作產(chǎn)生新的子代種群Qt+1， 將Pt+1與Qt+1合并形成新的種
群Rt， 重復(fù)以上操作， 直到滿足程序結(jié)束的條件

??Matlab版本的數(shù)據(jù)案例（重點以擁擠度函數(shù)與非支配排序為主），網(wǎng)上案例很多都是使用matlab，請自行查閱學(xué)習(xí)。

function [FrontValue,MaxFront] = NonDominateSort(FunctionValue,Operation)
% 進行非支配排序
% 輸入: FunctionValue, 待排序的種群(目標空間)，的目標函數(shù)
%       Operation,     可指定僅排序第一個面,排序前一半個體,或是排序所有的個體, 默認為排序所有的個體
% 輸出: FrontValue, 排序后的每個個體所在的前沿面編號, 未排序的個體前沿面編號為inf
%       MaxFront,   排序的最大前面編號if Operation == 1Kind = 2; elseKind = 1;  %√end[N,M] = size(FunctionValue);MaxFront = 0;cz = zeros(1,N);  %%記錄個體是否已被分配編號FrontValue = zeros(1,N)+inf; %每個個體的前沿面編號[FunctionValue,Rank] = sortrows(FunctionValue);%sortrows：由小到大以行的方式進行排序，返回排序結(jié)果和檢索到的數(shù)據(jù)(按相關(guān)度排序)在原始數(shù)據(jù)中的索引%開始迭代判斷每個個體的前沿面,采用改進的deductive sort，Deb非支配排序算法while (Kind==1 && sum(cz)<N) || (Kind==2 && sum(cz)<N/2) || (Kind==3 && MaxFront<1)MaxFront = MaxFront+1;d = cz;for i = 1 : Nif ~d(i)for j = i+1 : Nif ~d(j)k = 1;for m = 2 : Mif FunctionValue(i,m) > FunctionValue(j,m)  %比較函數(shù)值，判斷個體的支配關(guān)系k = 0;break;endendif k == 1d(j) = 1;endendendFrontValue(Rank(i)) = MaxFront;cz(i) = 1;endendend
end%% 非支配排序偽代碼
% def fast_nondominated_sort( P ):
% F = [ ]
% for p in P:
% Sp = [ ] #記錄被個體p支配的個體
% np = 0  #支配個體p的個數(shù)
% for q in P:
% if p > q:                #如果p支配q，把q添加到Sp列表中
% Sp.append( q )  #被個體p支配的個體
% else if p < q:        #如果p被q支配，則把np加1
% np += 1  #支配個體p的個數(shù)
% if np == 0:
% p_rank = 1        #如果該個體的np為0，則該個體為Pareto第一級
% F1.append( p )
% F.append( F1 )
% i = 0
% while F[i]:
% Q = [ ]
% for p in F[i]:
% for q in Sp:        #對所有在Sp集合中的個體進行排序
% nq -= 1
% if nq == 0:     #如果該個體的支配個數(shù)為0，則該個體是非支配個體
% q_rank = i+2    #該個體Pareto級別為當前最高級別加1。此時i初始值為0，所以要加2
% Q.append( q )
% F.append( Q )
% i += 1

function CrowdDistance = CrowdDistances(FunctionValue,FrontValue)
%分前沿面計算種群中每個個體的擁擠距離[N,M] = size(FunctionValue);CrowdDistance = zeros(1,N);Fronts = setdiff(unique(FrontValue),inf);for f = 1 : length(Fronts)Front = find(FrontValue==Fronts(f));Fmax = max(FunctionValue(Front,:),[],1);Fmin = min(FunctionValue(Front,:),[],1);for i = 1 : M[~,Rank] = sortrows(FunctionValue(Front,i));CrowdDistance(Front(Rank(1))) = inf;CrowdDistance(Front(Rank(end))) = inf;for j = 2 : length(Front)-1CrowdDistance(Front(Rank(j))) = CrowdDistance(Front(Rank(j)))+(FunctionValue(Front(Rank(j+1)),i)-FunctionValue(Front(Rank(j-1)),i))/(Fmax(i)-Fmin(i));endendend
end

?下面將給出Python版本的案例實現(xiàn)。

# Program Name: NSGA-II.py
# author；lxy
# Time:2023/03/11#Importing required modules
import math
import random
import matplotlib.pyplot as plt#First function to optimize
def function1(x):value = -x**2return value#Second function to optimize
def function2(x):value = -(x-2)**2return value#Function to find index of list,且是找到的第一個索引
def index_of(a,list):for i in range(0,len(list)):if list[i] == a:return ireturn -1#Function to sort by values 找出front中對應(yīng)值的索引序列
def sort_by_values(list1, values):sorted_list = []while(len(sorted_list)!=len(list1)):if index_of(min(values),values) in list1:sorted_list.append(index_of(min(values),values))values[index_of(min(values),values)] = math.infreturn sorted_list#Function to carry out NSGA-II's fast non dominated sort
def fast_non_dominated_sort(values1, values2):S=[[] for i in range(0,len(values1))]   #len(values1)個空列表front = [[]]n=[0 for i in range(0,len(values1))]rank = [0 for i in range(0, len(values1))]#將front0全部整理出來了，并未對front1-n等進行整理for p in range(0,len(values1)):S[p]=[]n[p]=0for q in range(0, len(values1)):if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (values1[p] > values1[q] and values2[p] >= values2[q]):if q not in S[p]:S[p].append(q)elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (values1[q] > values1[p] and values2[q] >= values2[p]):n[p] = n[p] + 1if n[p]==0:rank[p] = 0if p not in front[0]:front[0].append(p)i = 0#該循環(huán)能將所有的個體全部進行分類，顯然最后一層的個體中，沒有可以支配的個體了while(front[i] != []):Q=[]for p in front[i]:for q in S[p]:n[q] =n[q] - 1if( n[q]==0):rank[q]=i+1if q not in Q:Q.append(q)i = i+1front.append(Q)del front[len(front)-1]    #刪除了最后一層無支配個體的front層,最后一層是空集return front#Function to calculate crowding distance  同層之間的一個計算
def crowding_distance(values1, values2, front):# distance = [0 for i in range(len(front))]lenth= len(front)for i in range(lenth):distance = [0 for i in range(lenth)]sorted1 = sort_by_values(front, values1[:])  #找到front中的個體索引序列sorted2 = sort_by_values(front, values2[:])  #找到front中的個體索引序列distance[0] = 4444distance[lenth-1] = 4444for k in range(2,lenth-1):distance[k] = distance[k]+ (values1[sorted1[k+1]] - values1[sorted1[k-1]])/(max(values1)-min(values1))# print("/n")print("k:",k)print("distance[{}]".format(k),distance[k])for k in range(2,lenth-1):distance[k] = distance[k]+ (values2[sorted2[k+1]] - values2[sorted2[k-1]])/(max(values2)-min(values2))return distance# #Function to carry out the crossover
def crossover(a,b):r=random.random()if r>0.5:return mutation((a+b)/2)else:return mutation((a-b)/2)# #Function to carry out the mutation operator
def mutation(solution):mutation_prob = random.random()if mutation_prob <1:solution = min_x+(max_x-min_x)*random.random()return solution#Main program starts here
pop_size = 10
max_gen = 100#Initialization
min_x=-55
max_x=55
solution=[min_x+(max_x-min_x)*random.random() for i in range(0,pop_size)]
print('solution',solution)
gen_no=0
while(gen_no<max_gen):print('\n')print('gen_no:迭代次數(shù)',gen_no)function1_values = [function1(solution[i])for i in range(0,pop_size)]function2_values = [function2(solution[i])for i in range(0,pop_size)]print('function1_values:',function1_values)print('function2_values:', function2_values)non_dominated_sorted_solution = fast_non_dominated_sort(function1_values[:],function2_values[:])print('front',non_dominated_sorted_solution)# print("The best front for Generation number ",gen_no, " is")# for valuez in non_dominated_sorted_solution[0]:#     print("solution[valuez]",round(solution[valuez],3),end=" ")
#     print("\n")crowding_distance_values=[]for i in range(0,len(non_dominated_sorted_solution)):crowding_distance_values.append(crowding_distance(function1_values[:],function2_values[:],non_dominated_sorted_solution[i][:]))print("crowding_distance_values",crowding_distance_values)solution2 = solution[:]#Generating offspringswhile(len(solution2)!=2*pop_size):a1 = random.randint(0,pop_size-1)b1 = random.randint(0,pop_size-1)solution2.append(crossover(solution[a1],solution[b1]))print('solution2',solution2)function1_values2 = [function1(solution2[i])for i in range(0,2*pop_size)]function2_values2 = [function2(solution2[i])for i in range(0,2*pop_size)]non_dominated_sorted_solution2 = fast_non_dominated_sort(function1_values2[:],function2_values2[:])  #2*pop_sizeprint('non_dominated_sorted_solution2', non_dominated_sorted_solution2)# print("\n")crowding_distance_values2=[]for i in range(0,len(non_dominated_sorted_solution2)):crowding_distance_values2.append(crowding_distance(function1_values2[:],function2_values2[:],non_dominated_sorted_solution2[i][:]))print('crowding_distance_values2',crowding_distance_values2)new_solution= []for i in range(0,len(non_dominated_sorted_solution2)):non_dominated_sorted_solution2_1 = [index_of(non_dominated_sorted_solution2[i][j],non_dominated_sorted_solution2[i] ) for j in range(0,len(non_dominated_sorted_solution2[i]))]print('non_dominated_sorted_solution2_1:',non_dominated_sorted_solution2_1)front22 = sort_by_values(non_dominated_sorted_solution2_1[:], crowding_distance_values2[i][:])print("front22",front22)front = [non_dominated_sorted_solution2[i][front22[j]] for j in range(0,len(non_dominated_sorted_solution2[i]))]print('front',front)front.reverse()for value in front:new_solution.append(value)if(len(new_solution)==pop_size):breakif (len(new_solution) == pop_size):breaksolution = [solution2[i] for i in new_solution]gen_no = gen_no + 1
#
# Lets plot the final front now
function1 = [i * -1 for i in function1_values]
function2 = [j * -1 for j in function2_values]
plt.xlabel('Function 1', fontsize=15)
plt.ylabel('Function 2', fontsize=15)
plt.scatter(function1, function2)
plt.show()

?最終得出的Pareto前沿解的圖像為：

????????如還有其他問題，歡迎留言溝通。

????????后續(xù)應(yīng)該會針對具體的現(xiàn)實問題如生產(chǎn)調(diào)度問題，車輛路徑問題等進行算法實現(xiàn)。