什么是多目标优化算法，见链接：
十分钟了解完多目标优化算法
支配：假设小明9岁，50斤，小红8岁，45斤，小明无论是岁数还是体重都比小红大，所以小明支配小红。
互不支配：假设小明7岁，50斤，小红8岁，45斤，小明岁数比小红小，但体重比小红大，所以小明和小红互不支配。
帕累托集：在这个集合中，任意两个解互不支配。 如果有2个目标函数，帕累托集应该分布成一条曲线；如果有3个目标函数，帕累托集应该分布成一个超平面。 常规的2个目标函数，解法是目标加权得到是一个点，一个解；帕累托集，目标函数之间没有加权关系，所以得到是一条曲线。
非支配排序：将一组解分成n个集合：rank1,rank2…rankn,每个集合中所有的解都互不支配，但ranki中的任意解支配rankj中的任意解（i<j）.
多目标遗传算法是用来分析和解决多目标优化问题的一种进化算法，其核心就是协调各个目标函数之间的关系，找出使得各个目标函数都尽可能达到比较大的（或者比较小的）函数值的最优解集。在众多目标优化的遗传算法中，NSGA2算法是影响最大和应用范围最广的一种多目标遗传算法。在其出现后，由于它简单有效以及比较明显的优越性，使得该算法已经成为多目标优化问题中的基本算法之一。
介绍NSGA2，首先来介绍以下NSGA算法。
NSGA通过基于非支配排序的方法保留了种群中的优良个体,并且利用适应度共享函数保持了群体的多样性，取得了非常良好的效果。但实际工程领域中发现NSGA算法存在明显不足，这主要体现在如下3个方面：
为了克服非支配排序遗传算法(NSGA)的上述不足，印度科学家Deb于2002年在NSGA算法的基础上进行了改进，提出了带精英策略的非支配排序遗传算法(Elitist Non-Dominated Sorting Genetic Algorithm，NSGA-II)，NSGA-II 算法针对NSGA的缺陷通过以下三个方面进行了改进[16]：
解的支配关系与Pareto最优解
快速非支配排序步骤：
快速非支配排序就是将解集分解为不同次序的Pareto前沿的过程。
它可以描述为：
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: yudengwu(余登武）
# @Date : 2021/10/26
#@email:1344732766@qq.com
#============导入相关包========
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#==========定义两个目标函数=============
# 定义函数1
def function1(x):
value = -(x+2) ** 2+2*x
return value
# 定义函数2
def function2(x):
value = -(x - 2) ** 2
return value
#=========定义群体，并绘制初始解的分布图=====
pop_size = 10
max_gen = 100
# 迭代次数
#Initialization
min_x=-10
max_x=10
np.random.seed(10)#固定随机数种子，使每次生成的初始解集一样
solution=np.random.uniform(min_x,max_x,pop_size) #生成的初始解集
#函数1 对应的初始目标函数值
values1=map(function1,solution) #python中的map用法，可以对遍历一个列表如solution，然后将列表元素传入到函数function1。得到解。得到的格式为对象
values1=[i for i in values1] #因为上一步得到是对象格式，需要转换成列表格式
#函数2 对应的初始目标函数值
values2=map(function2,solution) #python中的map用法，可以对遍历一个列表如solution，然后将列表元素传入到函数function2。得到解。得到的格式为对象
values2=[i for i in values2] #因为上一步得到是对象格式，需要转换成列表格式
plt.scatter(values1,values2, s=20, marker='o')
for i in range(pop_size):
plt.annotate(i, xy=(values1[i], values2[i]), xytext=(values1[i] - 0.05, values2[i] - 0.05),fontsize=18)
plt.xlabel('function1')
plt.ylabel('function2')
plt.title('解的分布示意图')
plt.show()
#=================快速非支配排序==============
'''
1.n[p]=0 s[p]=[]
2.对所有个体进行非支配判断，若p支配q，则将q加入到S[p]中，并将q的层级提升一级。
若q支配p，n[p]+1.
3.找出种群中np=0的个体，即最优解，并找到最优解的支配解集合。存放到front[0]中
4 i==0
5.判断front是否为空，若不为空，将front中所有的个体sp中对应的被支配个体数减去1，(存放np==0的解序号进front[i+1]);i=i+1,跳到2;
若为空，则表明得到了所有非支配集合，程序结束
'''
values=[values1,values2] #解集【目标函数1解集，目标函数2解集...】
def fast_non_dominated_sort(values):
"""
优化问题一般是求最小值
:param values: 解集【目标函数1解集，目标函数2解集...】
:return:返回解的各层分布集合序号。类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]] 其中[1]表示Pareto 最优解对应的序号
"""
values11=values[0]#函数1解集
S = [[] for i in range(0, len(values11))]#存放 每个个体支配解的集合。
front = [[]] #存放群体的级别集合，一个级别对应一个[]
n = [0 for i in range(0, len(values11))]#每个个体被支配解的个数 。即针对每个解，存放有多少好于这个解的个数
rank = [np.inf for i in range(0, len(values11))]#存放每个个体的级别
for p in range(0, len(values11)):#遍历每一个个体
# ====得到各个个体 的被支配解个数 和支配解集合====
S[p] = [] #该个体支配解的集合 。即存放差于该解的解
n[p] = 0 #该个体被支配的解的个数初始化为0 即找到有多少好于该解的 解的个数
for q in range(0, len(values11)):#遍历每一个个体
less = 0 #的目标函数值小于p个体的目标函数值数目
equal = 0 #的目标函数值等于p个体的目标函数值数目
greater = 0 #的目标函数值大于p个体的目标函数值数目
for k in range(len(values)): # 遍历每一个目标函数
if values[k][p] > values[k][q]: # 目标函数k时，q个体值 小于p个体
less = less + 1 # q比p 好
if values[k][p] == values[k][q]: # 目标函数k时，p个体值 等于于q个体
equal = equal + 1
if values[k][p] < values[k][q]: # 目标函数k时，q个体值 大于p个体
greater = greater + 1 # q比p 差
if (less + equal == len(values)) and (equal != len(values)):
n[p] = n[p] + 1 # q比p, 比p好的个体个数加1
elif (greater + equal == len(values)) and (equal != len(values)):
S[p].append(q) # q比p差，存放比p差的个体解序号
#=====找出Pareto 最优解，即n[p]===0 的 个体p序号。=====
if n[p]==0:
rank[p] = 0 #序号为p的个体，等级为0即最优
if p not in front[0]:
# 如果p不在第0层中
# 将其追加到第0层中
front[0].append(p) #存放Pareto 最优解序号
# =======划分各层解========
"""
#示例，假设解的分布情况如下，由上面程序得到 front[0] 存放的是序号1
个体序号 被支配个数 支配解序号 front
1 0 2,3,4,5 0
2, 1， 3,4,5
3， 1， 4,5
4， 3， 5
5 4， 0
#首先 遍历序号1的支配解，将对应支配解[2,3,4,5] ，的被支配个数-1（1-1,1-1,3-1,4-1）
得到
表
个体序号 被支配个数 支配解序号 front
1 0 2,3,4,5 0
2, 0， 3,4,5
3， 0， 4,5
4， 2， 5
5 2， 0
#再令 被支配个数==0 的序号 对应的front 等级+1
得到新表...
"""
i = 0
while (front[i] != []): # 如果分层集合为不为空
Q = []
for p in front[i]: # 遍历当前分层集合的各个个体p
for q in S[p]: # 遍历p 个体 的每个支配解q
n[q] = n[q] - 1 # 则将fk中所有给对应的个体np-1
if (n[q] == 0):
# 如果nq==0
rank[q] = i + 1
if q not in Q:
Q.append(q) # 存放front=i+1 的个体序号
i = i + 1 # front 等级+1
front.append(Q)
del front[len(front) - 1] # 删除循环退出 时 i+1产生的[]
return front #返回各层 的解序号集合 # 类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]]
front=fast_non_dominated_sort(values)
#=================打印结果=======================
#遍历各层
for i in range(len(front)):
print('第%d层,解的序号为%s'%(i,front[i]))
jie=[]
for j in front[i]:#遍历第i层各个解
jie.append(solution[j])
print('第%d层,解为%s'%(i,jie))
结果如图，发现成功找到最优解
拥挤距离的定义
在NSGA-II中，为了衡量在同一个前沿中各个解质量的优劣，作者为每个解分配了一个拥挤距离，其背后的思想是让求得的Pareto最优解在objective space中尽量分散，也就有更大可能让解在Pareto最优前沿上均匀分布。
拥挤距离主要是维持种群中个体的多样性。具体而言，一般来说 是指种群按照支配关系进行非支配排后单个Rank层中个体的密集程度。常用于支配关系的多目标算法中，例如NSGA-1I.
主要步骤如下:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: yudengwu(余登武）
# @Date : 2021/10/27
#@email:1344732766@qq.com
#============导入相关包========
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#==========定义两个目标函数=============
# 定义函数1
def function1(x):
y = (x - 1) ** 2
return y
# 定义函数2
def function2(x):
y = np.cos(x)
return y
#=========定义群体，并绘制初始解的分布图=====
pop_size = 50
max_gen = 100
# 迭代次数
#Initialization
min_x=-10
max_x=10
np.random.seed(10)#固定随机数种子，使每次生成的初始解集一样
solution=np.random.uniform(min_x,max_x,pop_size) #生成的初始解集
#函数1 对应的初始目标函数值
values1=map(function1,solution) #python中的map用法，可以对遍历一个列表如solution，然后将列表元素传入到函数function1。得到解。得到的格式为对象
values1=[i for i in values1] #因为上一步得到是对象格式，需要转换成列表格式
#函数2 对应的初始目标函数值
values2=map(function2,solution) #python中的map用法，可以对遍历一个列表如solution，然后将列表元素传入到函数function2。得到解。得到的格式为对象
values2=[i for i in values2] #因为上一步得到是对象格式，需要转换成列表格式
plt.scatter(values1,values2, s=20, marker='o')
for i in range(pop_size):
plt.annotate(i, xy=(values1[i], values2[i]), xytext=(values1[i] - 0.05, values2[i] - 0.05),fontsize=18)
plt.xlabel('function1')
plt.ylabel('function2')
plt.title('解的分布示意图')
plt.show()
#=================快速非支配排序==============
values=[values1,values2] #解集【目标函数1解集，目标函数2解集...】
def fast_non_dominated_sort(values):
"""
优化问题一般是求最小值
:param values: 解集【目标函数1解集，目标函数2解集...】
:return:返回解的各层分布集合序号。类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]] 其中[1]表示Pareto 最优解对应的序号
"""
values11=values[0]#函数1解集
S = [[] for i in range(0, len(values11))]#存放 每个个体支配解的集合。
front = [[]] #存放群体的级别集合，一个级别对应一个[]
n = [0 for i in range(0, len(values11))]#每个个体被支配解的个数 。即针对每个解，存放有多少好于这个解的个数
rank = [np.inf for i in range(0, len(values11))]#存放每个个体的级别
for p in range(0, len(values11)):#遍历每一个个体
# ====得到各个个体 的被支配解个数 和支配解集合====
S[p] = [] #该个体支配解的集合 。即存放差于该解的解
n[p] = 0 #该个体被支配的解的个数初始化为0 即找到有多少好于该解的 解的个数
for q in range(0, len(values11)):#遍历每一个个体
less = 0 #的目标函数值小于p个体的目标函数值数目
equal = 0 #的目标函数值等于p个体的目标函数值数目
greater = 0 #的目标函数值大于p个体的目标函数值数目
for k in range(len(values)): # 遍历每一个目标函数
if values[k][p] > values[k][q]: # 目标函数k时，q个体值 小于p个体
less = less + 1 # q比p 好
if values[k][p] == values[k][q]: # 目标函数k时，p个体值 等于于q个体
equal = equal + 1
if values[k][p] < values[k][q]: # 目标函数k时，q个体值 大于p个体
greater = greater + 1 # q比p 差
if (less + equal == len(values)) and (equal != len(values)):
n[p] = n[p] + 1 # q比p, 比p好的个体个数加1
elif (greater + equal == len(values)) and (equal != len(values)):
S[p].append(q) # q比p差，存放比p差的个体解序号
#=====找出Pareto 最优解，即n[p]===0 的 个体p序号。=====
if n[p]==0:
rank[p] = 0 #序号为p的个体，等级为0即最优
if p not in front[0]:
# 如果p不在第0层中
# 将其追加到第0层中
front[0].append(p) #存放Pareto 最优解序号
# =======划分各层解========
i = 0
while (front[i] != []): # 如果分层集合为不为空
Q = []
for p in front[i]: # 遍历当前分层集合的各个个体p
for q in S[p]: # 遍历p 个体 的每个支配解q
n[q] = n[q] - 1 # 则将fk中所有给对应的个体np-1
if (n[q] == 0):
# 如果nq==0
rank[q] = i + 1
if q not in Q:
Q.append(q) # 存放front=i+1 的个体序号
i = i + 1 # front 等级+1
front.append(Q)
del front[len(front) - 1] # 删除循环退出 时 i+1产生的[]
return front #返回各层 的解序号集合 # 类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]]
front=fast_non_dominated_sort(values)
#=============拥挤距离================
def crowding_distance(values,front):
"""
:param values: 群体[目标函数值1，目标函数值2,...]
:param front: 群体解的等级，类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]]
:return: front 对应的 拥挤距离
"""
distance = np.zeros(shape=(pop_size, )) # 拥挤距离初始化为0
for rank in front: # 遍历每一层Pareto 解 rank为当前等级
for i in range(len(values)): # 遍历每一层函数值（先遍历群体函数值1，再遍历群体函数值2...）
valuesi = [values[i][A] for A in rank] # 取出rank等级 对应的 目标函数值i 集合
rank_valuesi = zip(rank, valuesi) # 将rank,群体函数值i集合在一起
sort_rank_valuesi = sorted(rank_valuesi, key=lambda x: (x[1],x[0])) # 先按函数值大小排序，再按序号大小排序
sort_ranki = [j[0] for j in sort_rank_valuesi] # 排序后当前等级rank
sort_valuesi = [j[1] for j in sort_rank_valuesi] # 排序后当前等级对应的 群体函数值i
#print(sort_ranki[0],sort_ranki[-1])
distance[sort_ranki[0]] = np.inf # rank 等级 中 的最优解 距离为inf
distance[sort_ranki[-1]] = np.inf # rank 等级 中 的最差解 距离为inf
#计算rank等级中，除去最优解、最差解外。其余解的拥挤距离
for j in range(1, len(rank) - 2):
distance[sort_ranki[j]] = distance[sort_ranki[j]] + (sort_valuesi[j + 1] - sort_valuesi[j - 1]) / (
max(sort_valuesi) - min(sort_valuesi)) # 计算距离
# 按照格式存放distances
distanceA = [[] for i in range(len(front))] #
for j in range(len(front)): # 遍历每一层Pareto 解 rank为当前等级
for i in range(len(front[j])): # 遍历给rank 等级中每个解的序号
distanceA[j].append(distance[front[j][i]])
return distanceA
distanceA=crowding_distance(values,front)
#打印拥挤距离结果
for i in range(len(front)):
print('当前等级 解的序号为:',front[i])
print('当前等级 解的拥挤距离为:',distanceA[i])
多目标优化过程中，得到父辈和子代解后，我们需要选择优秀的群体进行下一代。
这里采用精英选择策略
#这里的front，distance，solution 数据格式同前文
def elitism(self,front, distance, solution):
"""
精英选择策略
:param front: 父代与子代 组合构成的解的等级
:param distance: 父代与子代 组合构成的解 拥挤距离
:param solution: 父代与子代 组合构成的解
:return: 返回群体解。群体数量=（父代+子代）//2
"""
X1index = [] # 存储群体编号
pop_size = len(solution) // 2 # 保留的群体个数 即（父辈+子辈)//2
for i in range(len(front)): # 遍历各层
rank_distancei = zip(front[i], distance[i]) # 当前等级 与当前拥挤距离的集合
sort_rank_distancei = sorted(rank_distancei, key=lambda x: (x[1], x[0]),
reverse=True) # 先按拥挤距离大小排序，再按序号大小排序,逆序
sort_ranki = [j[0] for j in sort_rank_distancei] # 排序后当前等级rank
sort_distancei = [j[1] for j in sort_rank_distancei] # 排序后当前等级对应的 拥挤距离i
if (pop_size - len(X1index)) >=len(sort_ranki): # 如果X1index还有空间可以存放当前等级i 全部解
X1index.extend([A for A in sort_ranki])
#print('已存放len(X1index)', len(X1index))
#print('当前等级长度', len(sort_ranki))
#print('需要存放的总长度,popsize)
#num = pop_size-len(X1index)# X1index 还能存放的个数
elif len(sort_ranki) > (pop_size-len(X1index)): # 如果X1空间不可以存放当前等级i 全部解
num = pop_size - len(X1index)
X1index.extend([A for A in sort_ranki[0:num]])
X1 = [solution[i] for i in X1index]
return X1
本文通定义两个目标函数，然后生成初始解。利用初始解，解释了核心知识点：快速非支配排序、拥挤距离、精英选择策略。本文的文字部分有借鉴他人博客。代码为自己纯手打。下一步，将快速非支配排序、拥挤距离、精英选择策略写进遗传算法,解决常规约束、带复杂约束问题（可能需要修改下快速非支配排序代码，因为复杂约束不仅仅考虑适应度，还有惩罚项在里面）构成多目标优化问题代码模板。
多目标遗传优化算法nsga2[python源码实现]
多目标遗传优化算法nsga2求解复杂约束问题【python源码实现】
经实验：快速非支配排序、拥挤距离、精英选择策略 无法应用于粒子群算法。因为无法定义pbest,gbest。
作者：电气-余登武。原创不易，禁止抄袭。
