NSGA-Ⅱ算法Python实现详解：完整代码与超详细笔记

参考文章：

多目标进化算法——NSGA-II（python实现）_nsga python-CSDN博客

声明：我是学习这个作者的代码，代码都是他原创的，我只是在学习过程中写了一些笔记帮助理解

✅1. 使用OOP定义解类型

✨为什么 NSGA-II 用面向对象编程（Object-Oriented Programming） 好？

1. 清晰建模：用类描述“个体”、“种群”、“问题”等，结构清晰。

2. 模块复用：个体类可以在遗传操作（变异、交叉）、选择、评估中反复使用。

3. 封装性好：Individual 的内部实现可以变动，而算法主体不用改。

首先定义一个“类”（Class）叫 Individual，表示 NSGA-II 算法中的一个“个体”。它是现实问题中“一个解”的抽象。

class Individual(object):
    def __init__(self):

 ✨类和对象的区别

class Individual 是一个 类：是“模板”或“蓝图”。

individual = Individual() 是一个 对象（实例）：是“具体的个体”。

✨（1）接下来为这个类添加属性，使得使用这个类生成的所有实例都具备这些属性。

# 面向对象编程：定义一个类，叫做Individual，后续就可以使用Ind1=Individual()来生成一个实例
class Individual(object): 
    def __init__(self):
        self.solution = None
        # 给self一个属性solution,并且赋值为none
        # solution代表个体的决策变量向量（通常是一个 numpy 数组），即一个解的具体值。
        # 后续将赋值为形如 np.array([x1, x2, ..., xn]) 的形式。
        self.objective = defaultdict()
        # 给self一个属性为objective，表示个体的目标函数值，可以是多个目标
        # defaultdict()表示objective属性可以是任何类型目标的
        self.n = 0  # 解p被几个解所支配，是一个数值
        self.rank = 0  # 解所在第几层。如果rank = 0 表示第 1 层（最优非支配解）
        self.S = []  # 解p支配哪些解，是一个解集合
        self.distance = 0  # 拥挤度距离

注意这里写class Individual(object)是python2中遗留的习惯：

在 Python 2 中：

class A: 会创建一个“旧式类”

class A(object): 才是“新式类”——支持 super(), property, __slots__ 等高级特性

在 Python 3 中：

所有类都是新式类，即使写 class Individual:，它也隐式继承自 object

所以 class Individual: 和 class Individual(object): 是等价的

def __init__(self)是构造函数，当使用ind1 = Individual() 生成实例时，会自动调用这个函数（方法）

之后我们就可以创建很多个这样的对象，例如

ind1 = Individual()
ind2 = Individual()

就可以得到两个完全独立的解，互不影响

✨（2）上述部分定义了这个类的属性，接下来在这个class中继续完善，定义这个类的行为：

    # 对解的每个分量进行边界处理，确保它们在 [bound_min, bound_max] 区间内
    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_max
            elif item < bound_min:
                self.solution[i] = bound_min
    # 输入一个目标函数 objective_fun，用当前的 solution 计算目标值
    def calculate_objective(self, objective_fun):
        self.objective = objective_fun(self.solution)

    # 重载小于号“<”，定义解之间的排序规则（比较两个解谁更好）
    def __lt__(self, other):
        v1 = list(self.objective.values())
        v2 = list(other.objective.values())
        for i in range(len(v1)):
            if v1[i] > v2[i]:
                return 0  # 但凡有一个位置是 v1大于v2的 直接返回0,如果相等的话比较下一个目标值
        return 1
        # 如果所有目标都小于等于对方（即该解在所有目标上不劣于对方），返回 1。

✨ 举个类比（类 = 模板，人 = 对象）：

假设你设计一个 Person 类：

属性（状态）：

self.name = "Alice"

self.age = 25

方法（行为）：

def birthday(self): —— 让这个人过生日 self.age += 1

def speak(self): —— 让这个人说话

所以方法是对象自己可以“做的事情”。

✨ 回到 Individual 类，这些方法各自有非常清晰的目的：

方法名	功能	为什么要定义在类里？
bound_process()	对当前个体的 solution 做边界修正	这属于“个体自己的行为”，不同个体的解不同，自然该行为由个体来执行
calculate_objective()	用当前个体的变量去计算目标值	每个个体都要单独评估它的目标函数值
__lt__()	定义个体间的比较规则（用于排序）	比较两个个体，必须访问它们各自的 objective 值

✅2.主函数部分

NSGA-II 主流程框架分为 4 个阶段：

1. 初始化种群
2. 非支配排序 + 拥挤度计算
3. 选择 + 交叉 + 变异产生子代
4. 种群更新 & 可视化迭代

（1）第一阶段：初始化种群

def main():
    # ----------- 第一阶段：参数设置 + 初始化种群 ----------- 
    generations = 250  # 迭代次数
    popnum = 100  # 种群大小
    eta = 1  # 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1

    # poplength = 30  # 单个个体解向量的维数
    # bound_min = 0  # 定义域
    # bound_max = 1
    # objective_fun = ZDT1

    poplength = 3  # 单个个体解向量的维数
    bound_min = -5  # 定义域
    bound_max = 5
    objective_fun = KUR  # 指定目标函数（比如 KUR 是一个2目标测试函数）
    
    # 生成第一代种群
    P = []
    for i in range(popnum):
        P.append(Individual())
        P[i].solution = np.random.rand(poplength) * (bound_max - bound_min) + bound_min  # 随机生成个体可行解
        # 生成 popnum 个体，初始化其 solution 为在 [bound_min, bound_max] 区间内的随机值。
        P[i].bound_process(bound_min, bound_max)  # 修正越界变量。
        P[i].calculate_objective(objective_fun)  # 计算目标函数值

❶ eta = 1 是什么？为什么设置成 1？

这是 NSGA-II 中模拟二进制交叉（SBX，Simulated Binary Crossover） 和 多项式变异（Polynomial Mutation） 的分布指数参数，通常记作：η_c 或 η_m。

它决定了交叉或变异产生的子代“接近父代”的程度”。

在交叉时：

子代 = 父代 ± 某个扰动

eta 越大 → 产生的小变动越多 → 子代接近父代

eta 越小 → 产生的大变动越多 → 子代更分散、跳跃性更强

为什么 Deb 建议设为 1？

这是经验值，来自 Deb 博士（NSGA-II 作者）和大量实验：

"Setting η = 1 gives a good balance between exploration (跳出局部) and exploitation (局部精细搜索)."

❷ objective_fun = KUR 是什么？

这是在设置 目标函数（多目标测试函数），也叫 benchmark function。

KUR 是 Kurasawe function

是经典的 2 目标测试函数之一

通常定义在 [-5, 5] 的定义域中

有多个局部最优，适合测试算法的探索能力和 Pareto 前沿的逼近能力

KUR 的形式如下（3个变量）

def KUR(x):
    f = defaultdict(float)
    poplength = len(x)

    f[1] = 0
    f[2] = 0
    for 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)

所以代码中也可以看到poplength = 3 ，表明我们解的是一个 3维变量

✅注意这里代码的写法：

# 生成第一代种群
    P = []
    for i in range(popnum):
        P.append(Individual())
        P[i].solution = np.random.rand(poplength) * (bound_max - bound_min) + bound_min  # 随机生成个体可行解
        # 生成 popnum 个体，初始化其 solution 为在 [bound_min, bound_max] 区间内的随机值。
        P[i].bound_process(bound_min, bound_max)  # 修正越界变量。
        P[i].calculate_objective(objective_fun)  # 计算目标函数值

P[i].calculate_objective(objective_fun) 执行之后，目标函数值保存到哪里了？

✔️ 没错，它把目标函数值 计算出来并保存在了这个个体自身的属性 self.objective 中

那么objective_fun这个目标函数又在哪里定义的？

objective_fun = KUR  # 指定目标函数（比如 KUR 是一个2目标测试函数）

因此这部分的流程是：

把 KUR 函数作为一个变量（函数本身）传给 calculate_objective()

在 calculate_objective() 中，用 objective_fun(self.solution) 调用它

把返回的值保存到 self.objective

（2）第二阶段：非支配排序计算

该函数将种群中的个体根据支配关系划分为多个“非支配层”，也称为“Pareto 前沿”，并标记每个个体的 rank 值（表示它在哪一层）。

✅ 非支配排序的核心概念

我们说：

个体 p 支配 个体 q（记作 p < q），如果：

在所有目标上 p 都不比 q 差，并且

在至少一个目标上 p 比 q 更好

所以：

第一层（F₁）：不被任何人支配的个体

第二层（F₂）：只被 F₁ 中的个体支配的个体

第三层（F₃）：只被前两层支配的个体

依此类推…

注意这里的 "<" 是在class中重新定义过的

for i in range(len(v1)):
            if v1[i] > v2[i]:
                return 0  # 但凡有一个位置是 v1大于v2的 直接返回0,如果相等的话比较下一个目标值
        return 1
        # 如果所有目标都小于等于对方（即该解在所有目标上不劣于对方），返回 1。

来看代码：

def fast_non_dominated_sort(P):
    """
    非支配排序
    :param P: 种群 P
    :return F: F=(F_1, F_2, ...) 将种群 P 分为了不同的层， 返回值类型是dict，键为层号，值为 List 类型，存放着该层的个体
    """
    F = defaultdict(list) # 用一个字典 F 来存储所有 Pareto 层

    # 先来看哪些可以排进rank=1
    for p in P: # 遍历P中每个p
        p.S = []  # 存放被p支配的解
        p.n = 0   # 被多少个解支配
        for q in P: # 查看其他P中的其他个体q
            if p < q:  # p越小越好，if p dominate q
            # 注意这里的 < 是在clss中重新定义过的，要求是p中的所有值都小于等于q，才视为p<q
                p.S.append(q)  # 说明q被p支配，放入p.S中
            elif q < p:
                p.n += 1  # Increment the domination counter of p
        if p.n == 0: # 查看完所有的q，发现都没有能够支配p的
            p.rank = 1 # 那么p的rank=1
            F[1].append(p)
    
    i = 1
    while F[i]: # 只要F[i]不为空，这个循环就会继续
        Q = []
        for p in F[i]:
            for q in p.S:  # q属于p.S（被p支配的解）
                q.n = q.n - 1  #q.n 表示q被多少个解支配，现在-1表示解除p对他的支配
                if q.n == 0: # 如果解除之后发现，他就不被任何其他解支配了
                    q.rank = i + 1 # 那么这个解就属于下一个rank了
                    Q.append(q) # 此时这些新一轮rank的q被放入Q中
        i = i + 1
        F[i] = Q  # Q成为新的F[i]

    return F

（3）第二阶段：生成新个体（选择（非支配排序+拥挤度）-交叉-变异）

Q = make_new_pop(P, eta, bound_min, bound_max, objective_fun)
    # 基于选择、交叉、变异机制生成一代新的“子代种群” Q

# 生成新种群：选择-交叉-变异
def make_new_pop(P, eta, bound_min, bound_max, objective_fun):
    """
    use select,crossover and mutation to create a new population Q
    :param P: 父代种群
    :param eta: 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1
    :param bound_min: 定义域下限
    :param bound_max: 定义域上限
    :param objective_fun: 目标函数
    :return Q : 子代种群
    """
    popnum = len(P)
    Q = [] # 空的子代种群列表
    # binary tournament selection 二元锦标赛选择
    for k in range(int(popnum / 2)): # 循环(popnum / 2)次即可，因为每次生成2个个体
        # 从种群中随机选择两个个体，进行二元锦标赛，选择出一个 parent1
        i = random.randint(0, popnum - 1)
        j = random.randint(0, popnum - 1)
        parent1 = binary_tournament(P[i], P[j])

        # 从种群中随机选择两个个体，进行二元锦标赛，选择出一个 parent2
        i = random.randint(0, popnum - 1)
        j = random.randint(0, popnum - 1)
        parent2 = binary_tournament(P[i], P[j])

        while (parent1.solution == parent2.solution).all():  # 如果选择到的两个父代完全一样，则重选另一个
            i = random.randint(0, popnum - 1)
            j = random.randint(0, popnum - 1)
            parent2 = binary_tournament(P[i], P[j])

        # parent1 和 parent1 进行交叉，变异 产生 2 个子代
        Two_offspring = crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun)

        # 产生的子代进入子代种群
        Q.append(Two_offspring[0])
        Q.append(Two_offspring[1])
    return Q

🎯其中，二元锦标赛怎么个事儿呢？

def binary_tournament(ind1, ind2):
    """
    二元锦标赛
    :param ind1:个体1号
    :param ind2: 个体2号
    :return:返回较优的个体
    """
    if ind1.rank != ind2.rank:  # 如果两个个体有支配关系，即在两个不同的rank中，选择rank小的
        return ind1 if ind1.rank < ind2.rank else ind2
    elif ind1.distance != ind2.distance:  # 如果两个个体rank相同，比较拥挤度距离，选择拥挤读距离大的
        return ind1 if ind1.distance > ind2.distance else ind2
    else:  # 如果rank和拥挤度都相同，返回任意一个都可以
        return ind1

🎯拥挤度又是什么？

拥挤度距离 是用来衡量一个解在“同一 Pareto 层”中是否“孤立”或“稀疏”的指标，
它越大（注意是越大），说明这个解周围越稀疏，越有代表性，算法更倾向保留它。

NSGA-II 中，我们不仅要追求最优解，还要保证解的多样性。

所以我们：

用 Pareto rank 保证解的质量（越靠前越好）

用 拥挤度 distance 保证解的分布广泛、均匀

可以发现一个小细节：对于第一代种群，其实并没有计算distance,而是直接使用了在生成实例时赋予的默认值0。这是种常见但非必须的实践策略。

在 Deb 的原始论文（2002）中，NSGA-II 的流程如下：

1. 初始种群 P_0 随机生成

2. 对 P_0：

3. 做非支配排序 → 得到 rank

4. 对每一层计算 拥挤度 distance

5. 从 P_0 选择父母 → 生成 Q_0

所以在理论上，第一代也应该计算拥挤度。

那这个实现为什么没做？

这其实是编码实践中常见的一种简化处理，原因如下：

 原因 1：第一代选择时 rank 就够用了

在 binary_tournament() 中，如果两者 rank 不同，就只看 rank

第一代的选择主要目的是快速构建初始基因多样性

.distance 对于前期探索意义不大

 原因 2：节省一点点计算（尤其在目标函数计算很耗时时）

拥挤度计算涉及多维目标的排序和归一化，虽然不是复杂操作，但如果只是第一代，意义不大

 原因 3：实现逻辑简单

省去了第一代就引入多样性计算的复杂结构

避免了在未分层（或只一层）情况下拥挤度计算的“边界处理”

✅ 接下来看 交叉-变异

给两个父代（parent1, parent2）做 SBX 交叉，生成两个子代，再对其中一个子代做 PM 变异，最后修复越界并计算目标函数值。

beta 是 SBX 的核心参数，用于控制子代在父代之间的偏移量

eta 控制这个偏移范围（越大越靠近父代）

# 交叉&变异
def crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun):
    """
    交叉方式使用二进制交叉算子（SBX），变异方式采用多项式变异（PM）
    :param parent1: 父代1
    :param parent2: 父代2
    :param eta: 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1
    :param bound_min: 定义域下限
    :param bound_max: 定义域上限
    :param objective_fun: 目标函数
    :return: 2 个子代
    """
    poplength = len(parent1.solution)

    offspring1 = Individual() # 建两个新 Individual 实例
    offspring2 = Individual()
    offspring1.solution = np.empty(poplength) # 初始化其 solution 向量，准备存储交叉结果
    offspring2.solution = np.empty(poplength)

    # 交叉：SBX 模拟二进制交叉
    for i in range(poplength):
        rand = random.random()
        beta = (rand * 2) ** (1 / (eta + 1)) if rand < 0.5 else (1 / (2 * (1 - rand))) ** (1.0 / (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])

    # 多项式变异
    # TODO 变异的时候只变异一个，不要两个都变，不然要么出现早熟现象，要么收敛速度巨慢 why？
    # 如果两个子代都变异，可能破坏了交叉产生的优势结构（尤其早期）
    
    for i in range(poplength):
        mu = random.random()
        delta = 2 * mu ** (1 / (eta + 1)) if mu < 0.5 else (1 - (2 * (1 - mu)) ** (1 / (eta + 1)))
        offspring1.solution[i] = offspring1.solution[i] + delta

    # 定义域越界处理
    offspring1.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]

🎯 举个完整三维 SBX 的例子（手动计算）

▶️ 假设两个父代如下：

parent1.solution = [2.0, 5.0, 7.0]
parent2.solution = [4.0, 3.0, 9.0]
eta = 1

👉 我们将为每一维单独生成一个 rand 并计算出相应的 beta 和子代值。

第一维：

x1 = 2.0, x2 = 4.0
rand = 0.3 → beta = (2*0.3)^(1/(1+1)) = 0.7746

offspring1_x = 0.5 * [(1+0.7746)*2 + (1-0.7746)*4] ≈ 2.225
offspring2_x = 0.5 * [(1-0.7746)*2 + (1+0.7746)*4] ≈ 3.775

第二维：

x1 = 5.0, x2 = 3.0
rand = 0.8 → beta = (1 / (2 * (1 – 0.8)))^(1/2) = (1 / 0.4)^0.5 = 1.58

offspring1_y = 0.5 * [(1+1.58)*5 + (1-1.58)*3] = 0.5 * [12.9 + (-1.74)] ≈ 5.58
offspring2_y = 0.5 * [(1-1.58)*5 + (1+1.58)*3] ≈ 2.42

第三维：

x1 = 7.0, x2 = 9.0
rand = 0.45 → beta = (2*0.45)^(1/2) = 0.9487

offspring1_z = 0.5 * [(1+0.9487)*7 + (1-0.9487)*9] ≈ 7.05
offspring2_z = 0.5 * [(1-0.9487)*7 + (1+0.9487)*9] ≈ 8.95

最终结果

维度	父代1	父代2	子代1（offspring1）	子代2（offspring2）
x	2.0	4.0	≈ 2.225	≈ 3.775
y	5.0	3.0	≈ 5.58	≈ 2.42
z	7.0	9.0	≈ 7.05	≈ 8.95

📌变异的时候只变异一个，不要两个都变，不然要么出现早熟现象，要么收敛速度巨慢 why？

如果两个子代都变异，可能破坏了交叉产生的优势结构（尤其早期）

全部变异容易出现：

早熟：局部陷入极值

过慢：全局探索太强，收敛速度下降

保持“一个稳定，一个探索”是更合理的折中

（4）第四阶段：更新迭代

#  ----------- 第四阶段：种群更新 & 可视化迭代  ----------- 
    P_t = P  # 当前这一届的父代种群
    Q_t = Q  # 当前这一届的子代种群

    for gen_cur in range(generations):  # 开始迭代
        R_t = P_t + Q_t  # 父代 + 子代 合并为总种群
        F = fast_non_dominated_sort(R_t)  # 非支配排序，得到多个 Pareto 前沿
        # fast_non_dominated_sort函数return的是 F，也就是多个 Pareto 层 F[1], F[2], ...

        # 精华部分：选择下一代的父代
        P_n = []  # 即为P_t+1,表示下一届的父代
        i = 1
        while len(P_n) + len(F[i]) < popnum:  # 从F[1]开始，逐层加入到P_n中，直到到达规定个体数popnum
            crowding_distance_assignment(F[i])  # calculate crowding-distance in F_i
            P_n = P_n + F[i]  # include ith non dominated front in the parent pop
            i = i + 1  # check the next front for inclusion
        F[i].sort(key=lambda x: x.distance)  # sort in descending order using <n，因为本身就在同一层，所以相当于直接比拥挤距离
        P_n = P_n + F[i][:popnum - len(P_n)]
        # 保证种群规模固定为 popnum，在最后一层内按拥挤度挑选
        Q_n = make_new_pop(P_n, eta, bound_min, bound_max,
                           objective_fun)  # use selection,crossover and mutation to create a new population Q_n

        # 求得下一届的父代和子代成为当前届的父代和子代，，进入下一次迭代 《=》 t = t + 1
        P_t = P_n
        Q_t = Q_n

        # 绘图
        plt.clf()
        plt.title('current generation:' + str(gen_cur + 1))
        plot_P(P_t)
        plt.pause(0.1)

    plt.show()

注意这里有个细节：

F[i].sort(key=lambda x: x.distance)  # sort in descending order using <n，因为本身就在同一层，所以相当于直接比拥挤距离
        P_n = P_n + F[i][:popnum - len(P_n)]

当最后一层可以加入P_n的F[i]加入之后，P_n中的个体数量完全有可能超出了规定的popnum

此时应该做一个筛选

F[i] 是一个个体列表（Pareto 第 i 层的所有个体）

.sort() 是 Python 列表的排序方法

key=... 是告诉 .sort() 你按照什么字段来排？

lambda x: x.distance 是一个匿名函数，意思是：

def key_function(x):
    return x.distance

所以也可以等价地写成：

F[i].sort(key=lambda individual: individual.distance)

然后再取到满足popnum个即可

完整代码

'''
本程序用于学习NSGA-II算法
'''

from collections import defaultdict
import numpy as np
import random
import math
import matplotlib.pyplot as plt

# 面向对象编程：定义一个类，叫做Individual，后续就可以使用Ind1=Individual()来生成一个实例
class Individual(object): 
    def __init__(self):
        self.solution = None
        # 给self一个属性solution,并且赋值为none
        # solution代表个体的决策变量向量（通常是一个 numpy 数组），即一个解的具体值。
        # 后续将赋值为形如 np.array([x1, x2, ..., xn]) 的形式。
        self.objective = defaultdict()
        # 给self一个属性为objective，表示个体的目标函数值，可以是多个目标
        # defaultdict()表示objective属性可以是任何类型目标的
        self.n = 0  # 解p被几个解所支配，是一个数值
        self.rank = 0  # 解所在第几层。如果rank = 0 表示第 1 层（最优非支配解）
        self.S = []  # 解p支配哪些解，是一个解集合
        self.distance = 0  # 拥挤度距离

    # 对解的每个分量进行边界处理，确保它们在 [bound_min, bound_max] 区间内
    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_max
            elif item < bound_min:
                self.solution[i] = bound_min
    # 输入一个目标函数 objective_fun，用当前的 solution 计算目标值
    def calculate_objective(self, objective_fun):
        self.objective = objective_fun(self.solution)

    # 重载小于号“<”，定义解之间的排序规则（比较两个解谁更好）
    def __lt__(self, other):
        v1 = list(self.objective.values())
        v2 = list(other.objective.values())
        for i in range(len(v1)):
            if v1[i] > v2[i]:
                return 0  # 但凡有一个位置是 v1大于v2的 直接返回0,如果相等的话比较下一个目标值
        return 1
        # 如果所有目标都小于等于对方（即该解在所有目标上不劣于对方），返回 1。

# 双目标函数KUR
def KUR(x):
    f = defaultdict(float)
    poplength = len(x)

    f[1] = 0
    f[2] = 0
    for 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 f

def fast_non_dominated_sort(P):
    """
    非支配排序
    :param P: 种群 P
    :return F: F=(F_1, F_2, ...) 将种群 P 分为了不同的层， 返回值类型是dict，键为层号，值为 List 类型，存放着该层的个体
    """
    F = defaultdict(list) # 用一个字典 F 来存储所有 Pareto 层

    # 先来看哪些可以排进rank=1
    for p in P: # 遍历P中每个p
        p.S = []  # 存放被p支配的解
        p.n = 0   # 被多少个解支配
        for q in P: # 查看其他P中的其他个体q
            if p < q:  # p越小越好，if p dominate q
            # 注意这里的 < 是在clss中重新定义过的，要求是p中的所有值都小于等于q，才视为p<q
                p.S.append(q)  # 说明q被p支配，放入p.S中
            elif q < p:
                p.n += 1  # Increment the domination counter of p
        if p.n == 0: # 查看完所有的q，发现都没有能够支配p的
            p.rank = 1 # 那么p的rank=1
            F[1].append(p)
    
    i = 1
    while F[i]: # 只要F[i]不为空，这个循环就会继续
        Q = []
        for p in F[i]:
            for q in p.S:  # q属于p.S（被p支配的解）
                q.n = q.n - 1  #q.n 表示q被多少个解支配，现在-1表示解除p对他的支配
                if q.n == 0: # 如果解除之后发现，他就不被任何其他解支配了
                    q.rank = i + 1 # 那么这个解就属于下一个rank了
                    Q.append(q) # 此时这些新一轮rank的q被放入Q中
        i = i + 1
        F[i] = Q  # Q成为新的F[i]

    return F

def crowding_distance_assignment(L):
    """ 传进来的参数应该是L = F(i)，类型是List"""
    l = len(L)  # number of solution in F

    for i in range(l):
        L[i].distance = 0  # initialize distance

    for m in L[0].objective.keys():
        L.sort(key=lambda x: x.objective[m])  # sort using each objective value
        L[0].distance = float('inf')
        L[l - 1].distance = float('inf')  # so that boundary points are always selected

        # 排序是由小到大的，所以最大值和最小值分别是 L[l-1] 和 L[0]
        f_max = L[l - 1].objective[m]
        f_min = L[0].objective[m]


        # 当某一个目标方向上的最大值和最小值相同时，此时会发生除零错，这里采用异常处理机制来解决
        try:
            for i in range(1, l - 1):  # for all other points
                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_tournament(ind1, ind2):
    """
    二元锦标赛
    :param ind1:个体1号
    :param ind2: 个体2号
    :return:返回较优的个体
    """
    if ind1.rank != ind2.rank:  # 如果两个个体有支配关系，即在两个不同的rank中，选择rank小的
        return ind1 if ind1.rank < ind2.rank else ind2
    elif ind1.distance != ind2.distance:  
        # 如果两个个体rank相同，比较拥挤度距离，选择拥挤读距离大的
        # 注意在class中已经定义了属性self.distance = 0
        #   
        return ind1 if ind1.distance > ind2.distance else ind2
    else:  # 如果rank和拥挤度都相同，返回任意一个都可以
        return ind1


# 交叉&变异
def crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun):
    """
    交叉方式使用二进制交叉算子（SBX），变异方式采用多项式变异（PM）
    :param parent1: 父代1
    :param parent2: 父代2
    :param eta: 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1
    :param bound_min: 定义域下限
    :param bound_max: 定义域上限
    :param objective_fun: 目标函数
    :return: 2 个子代
    """
    poplength = len(parent1.solution)

    offspring1 = Individual() # 建两个新 Individual 实例
    offspring2 = Individual()
    offspring1.solution = np.empty(poplength) # 初始化其 solution 向量，准备存储交叉结果
    offspring2.solution = np.empty(poplength)

    # 交叉：SBX 模拟二进制交叉
    for i in range(poplength):
        rand = random.random()
        beta = (rand * 2) ** (1 / (eta + 1)) if rand < 0.5 else (1 / (2 * (1 - rand))) ** (1.0 / (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])

    # 多项式变异
    # TODO 变异的时候只变异一个，不要两个都变，不然要么出现早熟现象，要么收敛速度巨慢 why？
    # 如果两个子代都变异，可能破坏了交叉产生的优势结构（尤其早期）
    
    for i in range(poplength):
        mu = random.random()
        delta = 2 * mu ** (1 / (eta + 1)) if mu < 0.5 else (1 - (2 * (1 - mu)) ** (1 / (eta + 1)))
        offspring1.solution[i] = offspring1.solution[i] + delta

    # 定义域越界处理
    offspring1.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):
    """
    use select,crossover and mutation to create a new population Q
    :param P: 父代种群
    :param eta: 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1
    :param bound_min: 定义域下限
    :param bound_max: 定义域上限
    :param objective_fun: 目标函数
    :return Q : 子代种群
    """
    popnum = len(P)
    Q = [] # 空的子代种群列表
    # binary tournament selection 二元锦标赛选择
    for k in range(int(popnum / 2)): # 循环(popnum / 2)次即可，因为每次生成2个个体
        # 从种群中随机选择两个个体，进行二元锦标赛，选择出一个 parent1
        i = random.randint(0, popnum - 1)
        j = random.randint(0, popnum - 1)
        parent1 = binary_tournament(P[i], P[j])

        # 从种群中随机选择两个个体，进行二元锦标赛，选择出一个 parent2
        i = random.randint(0, popnum - 1)
        j = random.randint(0, popnum - 1)
        parent2 = binary_tournament(P[i], P[j])

        while (parent1.solution == parent2.solution).all():  # 如果选择到的两个父代完全一样，则重选另一个
            i = random.randint(0, popnum - 1)
            j = random.randint(0, popnum - 1)
            parent2 = binary_tournament(P[i], P[j])

        # parent1 和 parent1 进行交叉，变异 产生 2 个子代
        Two_offspring = crossover_mutation(parent1, parent2, eta, bound_min, bound_max, objective_fun)

        # 产生的子代进入子代种群
        Q.append(Two_offspring[0])
        Q.append(Two_offspring[1])
    return Q

def 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 show_some_ind(P):
    # 测试使用
    for i in P:
        print(i.solution)


"""
NSGA-II 主流程框架分为 4 个阶段：
1. 初始化种群
2. 非支配排序 + 拥挤度计算
3. 选择 + 交叉 + 变异产生子代
4. 种群更新 & 可视化迭代
"""

def main():
    # ----------- 第一阶段：参数设置 + 初始化种群 ----------- 
    generations = 250  # 迭代次数
    popnum = 100  # 种群大小
    eta = 1  # 变异分布参数，该值越大则产生的后代个体逼近父代的概率越大。Deb建议设为 1

    # poplength = 30  # 单个个体解向量的维数
    # bound_min = 0  # 定义域
    # bound_max = 1
    # objective_fun = ZDT1

    poplength = 3  # 单个个体解向量的维数
    bound_min = -5  # 定义域
    bound_max = 5
    objective_fun = KUR  # 指定目标函数（比如 KUR 是一个2目标测试函数）
    
    # 生成第一代种群
    P = []
    for i in range(popnum):
        P.append(Individual())
        P[i].solution = np.random.rand(poplength) * (bound_max - bound_min) + bound_min  # 随机生成个体可行解
        # 生成 popnum 个体，初始化其 solution 为在 [bound_min, bound_max] 区间内的随机值。
        P[i].bound_process(bound_min, bound_max)  # 修正越界变量。
        P[i].calculate_objective(objective_fun)  # 计算目标函数值
        # 注意这里的P仍然是种群，并且每个个体P[i]已经调用了class中的函数，拥有了solution和objective属性值

    # -----------  第二阶段：非支配排序 + 生成初代子代 ----------- 
    fast_non_dominated_sort(P)
    # 对种群进行快速非支配排序，标记每个个体的 rank

    #  -----------  第三阶段：选择 + 交叉 + 变异产生子代  ----------
    Q = make_new_pop(P, eta, bound_min, bound_max, objective_fun)
    # 注意，对于第一代的种群来说，此时还没有计算拥挤度，大家的拥挤度都是初始赋予的0.在后面的迭代中才计算了拥挤度
    # 这是一种常见但非必须的做法
    # 基于选择、交叉、变异机制生成一代新的“子代种群” Q


    #  ----------- 第四阶段：种群更新 & 可视化迭代  ----------- 
    P_t = P  # 当前这一届的父代种群
    Q_t = Q  # 当前这一届的子代种群

    for gen_cur in range(generations):  # 开始迭代
        R_t = P_t + Q_t  # 父代 + 子代 合并为总种群
        F = fast_non_dominated_sort(R_t)  # 非支配排序，得到多个 Pareto 前沿
        # fast_non_dominated_sort函数return的是 F，也就是多个 Pareto 层 F[1], F[2], ...

        # 精华部分：选择下一代的父代
        P_n = []  # 即为P_t+1,表示下一届的父代
        i = 1
        while len(P_n) + len(F[i]) < popnum:  # 从F[1]开始，逐层加入到P_n中，直到到达规定个体数popnum
            crowding_distance_assignment(F[i])  # calculate crowding-distance in F_i
            P_n = P_n + F[i]  # include ith non dominated front in the parent pop
            i = i + 1  # check the next front for inclusion
        F[i].sort(key=lambda x: x.distance)  # sort in descending order using <n，因为本身就在同一层，所以相当于直接比拥挤距离
        P_n = P_n + F[i][:popnum - len(P_n)]
        # 保证种群规模固定为 popnum，在最后一层内按拥挤度挑选
        Q_n = make_new_pop(P_n, eta, bound_min, bound_max,
                           objective_fun)  # use selection,crossover and mutation to create a new population Q_n

        # 求得下一届的父代和子代成为当前届的父代和子代，，进入下一次迭代 《=》 t = t + 1
        P_t = P_n
        Q_t = Q_n

        # 绘图
        plt.clf()
        plt.title('current generation:' + str(gen_cur + 1))
        plot_P(P_t)
        plt.pause(0.1)

    plt.show()

    return 0



if __name__ == '__main__':
    main()

作者：阿斑ABam