NSGA-II 遗传多目标算法（python示例）

一、前言

最近在准备毕业论文，研究了一下主流的多目标算法，对于NSGA-II，网上大部分代码是全部是面向过程来实现的，本人更喜欢采用面向对象的方式，故采用python面向对象实现了一个示例，实现了对于二元多目标问题的求解。

二、算法基本流程

三、核心思想

1、非支配排序

这个简单的例子说明了帕累托最优的概念。上面我们有4个成员A, B, C和D，有两个特征:身高和工资。现在，如果我们同时比较他们的身高和薪水，我们会发现这不是很直观，因为他们有多个目标。

既然这两个目标越大越好，我们可以简单地对它们进行比较。首先，我们观察到A和B都比C和D多，所以我们说A和B在身高和薪水上“支配”C和D。同理，C支配D，D可被A,B,C支配。

A和B呢?A比B高，但是工资低。相反，B面临着同样的情况。我们称这种情况为“非支配”。如果我们能找到一组解它们不互相支配，也不受其他解支配，我们称之为"帕累托最优"解。在上面的例子中，A和B都在帕累托最优前沿。

几个概念：

非支配解：假设任何二解S1 及S2 对所有目标而言，S1均优于S2，则我们称S1 支配S2，若S1 的解没有被其他解所支配，则S1 称为非支配解（不受支配解），也称Pareto解（帕雷托解）

支配解：若解S2的所有目标均劣于S1，则称S1优于S2，也称S1支配S2，S2为受支配解。

Pareto前沿面：找到所有Pareto解之后，这些解组成的平面叫做Pareto前沿面（Non-dominated front）。在目标函数较多时，前沿面通常为超曲面。

2、拥挤度

通俗的来讲，当需要舍弃某一rank平面的部分节点时，由于同一平面中的所有节点rank相同，不能通过rank来舍弃，而是要通过拥挤度来舍弃，以上就是拥挤度的作用。

算法更倾向于稀疏的点，也就是让节点更可能的分散，可以有效地方式早熟和过拟合现象

3、精英选择策略

每次都将父类与子类想结合，依次采用非支配排序、计算拥挤度来选择父代。

四、python实现

github: 源代码地址

"""
author: pym
time: 2023.10.29
ide: pycharm2
"""from collections import defaultdict
import numpy as np
import random
import matplotlib.pyplot as plt
import mathclass Individual(object):def __init__(self):self.solution = None    # 实际为nparray类型，方便四则运算self.objective = defaultdict()self.n = 0              # 解p被几个解支配self.rank = 0           # 解p所在层数self.S = []             # 解p支配解的集合self.distance = 0       # 拥挤度距离def bound_process(self, bound_min, bound_max):"""对解向量 solution 中的每个分量进行定义域判断；超过最大值赋为最大值:param bound_min: 定义域下限:param bound_max: 定义域上限:return:"""for i, item in enumerate(self.solution):if item > bound_max:self.solution[i] = bound_maxelif item < bound_min:self.solution[i] = bound_mindef calculate_objective(self, objective_fun):"""计算目标值:param objective_fun: 目标函数:return:"""self.objective = objective_fun(self.solution)def __lt__(self, other):"""重载小于号，只有当solution中全部小于对方，才判断小于:param other: 比较的个体:return: 1：小于 0：大于"""v1 = list(self.objective.values())v2 = list(other.objective.values())for i in range(len(v1)):if v1[i] > v2[i]:return 0return 1def fast_non_dominated_sort(P):"""非支配排序:param P: 种群P:return: F：分层结果，返回值类型为dict，键为层号，值为list（该层中的个体）"""F = defaultdict(list)for p in P:p.S = []p.n = 0for q in P:if p < q:       # p支配qp.S.append(q)elif q < p:     # q支配pp.n += 1if p.n == 0:p.rank = 1F[1].append(p)i = 1while F[i]:Q = []for p in F[i]:for q in p.S:q.n -= 1if q.n == 0:q.rank = i + 1Q.append(q)i += 1F[i] = Qreturn Fdef crowding_distance_assignment(L):"""计算拥挤度:param L: F[i]，是个list，为第i层的节点集合:return:"""l = len(L)# 初始化距离for i in range(l):L[i].distance = 0# 遍历每个目标方向（有几个优化目标，就有几个目标方向）for m in L[0].objective.keys():L.sort(key=lambda x: x.objective[m])    # 使用objective值排序L[0].distance = float('inf')L[l - 1].distance = float('inf')f_max = L[l - 1].objective[m]f_min = L[0].objective[m]# 当某一个目标方向上的最大值和最小值相同时，会出现除0错误try:for i in range(1, l - 1):L[i].distance = L[i].distance + (L[i + 1].objective[m] - L[i - 1].objective[m]) / (f_max - f_min)except Exception:print(str(m) + "目标方向上，最大值为：" + str(f_max) + " 最小值为：" + str(f_min))def binary_tornament(ind1, ind2):"""二元锦标赛：先选非支配排序靠前的，再选拥挤度低（即距离远）；如果都不行，则随机:param ind1: 个体1:param ind2: 个体1:return: 返回较优的个体"""if ind1.rank != ind2.rank:return ind1 if ind1.rank < ind2.rank else ind2elif ind1.distance != ind2.distance:return ind1 if ind1.distance > ind2.distance else ind2else:return ind1def crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun):"""交叉：二进制交叉算子（SBX），变异：多项式变异（PM）:param parent1: 父代1:param parent2: 父代2:param eta: 变异参数，越大则后代个体越逼近父代:return:"""poplength = len(parent1.solution)   # 解向量维数# 初始化两个后代个体offspring1 = Individual()offspring2 = Individual()offspring1.solution = np.empty(poplength)offspring2.solution = np.empty(poplength)# 二进制交叉for i in range(poplength):rand = random.random()if rand < 0.5:beta = (rand * 2) ** (1 / (eta + 1))else:beta = (1 / (2 * (1 - rand)))**(1 / (eta + 1))offspring1.solution[i] = 0.5 * ((1 + beta) * parent1.solution[i] + (1 - beta) * parent2.solution[i])offspring2.solution[i] = 0.5 * ((1 - beta) * parent1.solution[i] + (1 + beta) * parent2.solution[i])# 多项式变异for i in range(poplength):mu = random.random()if mu < 0.5:delta = 2 * mu ** (1 / (eta + 1))else:delta = (1 - (2 * (1 - mu)) ** (1 / (eta + 1)))# 只变异一个offspring1.solution[i] = offspring1.solution[i] + deltaoffspring1.bound_process(bound_min, bound_max)offspring2.bound_process(bound_min, bound_max)offspring1.calculate_objective(objective_fun)offspring2.calculate_objective(objective_fun)return [offspring1, offspring2]def make_new_pop(P, eta, bound_min, bound_max, objective_fun):"""选择交叉变异获得新后代:param P: 父代种群:param eta: 变异参数，越大则后代个体越逼近父代:param bound_min: 定义域下限:param bound_max: 定义域上限:param objective_fun: 目标函数:return: 子代种群"""popnum = len(P)     # 种群个数Q = []# 二元锦标赛选择for i in range(int(popnum / 2)):# 从种群中随机选择两个个体，进行二元锦标赛，选择一个parenti = random.randint(0, popnum - 1)j = random.randint(0, popnum - 1)parent1 = binary_tornament(P[i], P[j])parent2 = parent1while (parent1.solution == parent2.solution).all():     # 小细节alli = random.randint(0, popnum - 1)j = random.randint(0, popnum - 1)parent2 = binary_tornament(P[i], P[j])Two_offspring = crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun)Q.append(Two_offspring[0])Q.append(Two_offspring[1])return Qdef KUR(x):"""计算各个目标方向上的目标值:param x: 解向量:return: 字典：各个方向上的目标值（key：目标方向；value：目标值）"""f = defaultdict(float)poplength = len(x)f[1] = 0f[2] = 0for i in range(poplength - 1):f[1] = f[1] + (-10) * math.exp((-0.2) * (x[i] ** 2 + x[i + 1] ** 2) ** 0.5)for i in range(poplength):f[2] = f[2] + abs(x[i]) ** 0.8 + 5 * math.sin(x[i] ** 3)return fdef plot_P(P):"""给种群绘图:param P: 种群集合:return:"""X = []Y = []for ind in P:X.append(ind.objective[1])Y.append(ind.objective[2])plt.xlabel('F1')plt.ylabel('F2')plt.scatter(X, Y)def main():# 初始化参数generations = 250   # 迭代次数popnum = 100        # 种群大小eta = 1             # 变异分布参数poplength = 3       # 单个个体解向量的维数bound_min = -5bound_max = 5objective_fun = KUR# 生成第一代种群P = []for i in range(popnum):P.append(Individual())P[i].solution = np.random.rand(poplength) * (bound_max - bound_min) + bound_minP[i].bound_process(bound_min, bound_max)    # 越界处理P[i].calculate_objective(objective_fun)     # 计算目标值# 快速非支配排序fast_non_dominated_sort(P)Q = make_new_pop(P, eta, bound_min, bound_max, objective_fun)P_t = P     # 当前这一代的父代种群Q_t = Q     # 当前这一代的子代种群for gen_cur in range(generations):R_t = P_t + Q_tF = fast_non_dominated_sort(R_t)P_n = []    # 即为P_t+1，表示下一代的父代i = 1# 依次将最高级别的支配平面中的节点放入到P_n中，之后更新非支配，直到达到要求的规模while len(P_n) + len(F[i]) < popnum:crowding_distance_assignment(F[i])P_n += F[i]i += 1# 按照支配排序选完之后，再按照拥挤度来选择F[i].sort(key=lambda x: x.distance)P_n = P_n + F[i][:popnum - len(P_n)]Q_n = make_new_pop(P_n, eta, bound_min, bound_max, objective_fun)# 将下一届的父代和子代成为当前的父代和子代P_t = P_nQ_t = Q_n# 可视化plt.clf()plt.title("current generation: " + str(gen_cur + 1))plot_P(P_t)plt.pause(0.1)plt.show()return 0if __name__ == "__main__":main()