NSGA-Ⅱ(NSGA2)-多目标优化算法

NSGA-Ⅱ(Non-dominated Sorting Genetic Algorithms-Ⅱ)算法，即带有精英保留策略的快速非支配多目标优化算法，是一种基于Pareto最优解的多目标优化算法。在此感谢 A fast and elitist multiobjective genetic algorithm: NSGA-II 作者及其所著论文。

• NSGA(Non-dominated Sorting Genetic Algorithms)

NSGA-Ⅱ相比，NSGA：

1. 时间复杂度$\mathrm{O}(M{N}^3)$
2. 非精英机制
3. 需要指定共享参数

NSGA-Ⅱ的优势：

1. 提出了快速非支配排序算法 时间复杂度$\mathrm{O}(M{N}^2)$
2. 引入精英策略
3. 在NSGA2中用拥挤比较的方法替换了共享函数方法。要实现这两种方法，首先我们需要定义两个操作：拥挤度和拥挤比较算子

• Pareto最优解

帕累托(Pareto)最优解其实与上述的支配关系类似，即：假定固有的一群人和可分配的资源，从一种分配状态到另一种状态的变化中，在没有使任何人境况变坏的前提下，使得至少一个人变得更好，这就是帕累托改进，这样迭代下去得到最优解的过程被称为帕累托最优化。

• 快速非支配排序

首先想要单独说说什么是支配性策略，以下是百科解释：
        一个博弈当中，无论对手采取什么策略，你若有几个策略，而其中一个策略可以使你得到比采取其他策略更好的结果，那么，这个策略就是你的优势策略。
延伸出来的有一个非常著名的博弈论术语：纳什均衡。简单概括的话，有一个非常常见且现实的例子，就是：价格战。
纳什均衡：
        在一个博弈过程中，无论对方的策略选择如何，当事人一方都会选择某个确定的策略，则该策略被称作支配性策略。如果任意一位参与者在其他所有参与者的策略确定的情况下，其选择的策略是最优的，那么这个组合就被定义为纳什均衡。
        一个策略组合被称为纳什均衡，当每个博弈者的均衡策略都是为了达到自己期望收益的最大值，与此同时，其他所有博弈者也遵循这样的策略。
纳什均衡包括了纯策略纳什均衡 和 混合策略纳什均衡。所谓纯策略是指非0即1的策略概率，即某一策略选择概率为1，则其余策略不被选择，即概率为0。而混合策略是对策略集中的各策略分配了概率，即博弈者根据概率随机选择策略，在某种概率集合下，博弈者$a$可以达到最优效益（这个最优往往是动态的，因为现实情况是每位博弈者的决策相互影响）。

说回快速非支配排序。所谓的支配，是指在优化问题中，某一个解$x_i$在$N$维下的目标函数上均优于$x_j$，这样可以将二者关系称为：$x_i$支配$x_j$。
若$\nexists \forall x_j$优于$x_i$，则标记$x_i$为非支配个体。（即在$N$维目标函数下均不优于$x_i$）
通过上述步骤得到的是非支配个体集是种群的第$n$级非支配层，然后忽略已标记的非支配个体，再次进行该算法，则会得到第$n+1$级非支配层。

• 拥挤度和拥挤比较算子

拥挤度：指的是某解$x_i$周围附近的解的数量，这个数值以最近邻居作为顶点的长方形周长。通常通过拥挤度来保证种群多样性
具体计算步骤：
①初始化每个点的拥挤度：$i_d=0$
②对种群进行非支配排序，令边界的2个个体为无穷大，即：$O_d=I_d=\infty$
③对其余$(i-2)$个个体进行拥挤度计算：$$i_d=\sum _{j=1}^m(|f_j^{i+1}-f_j^{i-1}|)$$
其中：$i_d$表示$i$点的拥挤度，$f_j^{i+1}$表示$i+1$点的第$j$目标函数值，$f_j^{i-1}$表示$i-1$点的第$j$目标函数值。

拥挤比较算子：比较算子的前提是经过了快速非支配排序和拥挤度的计算之后，种群中的每个个体$i$都有两个属性——非支配排序得到的非支配序$i_{rank}$（第几级）和拥挤度$i_d$。
根据这两个属性，我们就可以定义拥挤比较算子：以个体$i$与个体$j$比较，若以下$\forall$条件成立，则个体$i$获胜。
①$i_{rank}<j_{rank}$，个体$i$所处的非支配层优于个体$j$。选择的个体属于较优的非劣等级。
②$i_{rank}=j_{rank}\cap i_d>j_d$，个体$i$和个体$j$拥有相同的层级，但个体$i$有更大的拥挤距离（拥挤度）。即选择较不拥挤区域的个体。

• 精英保留策略

简要来说就是：对遗传算法来说，能否收敛到全局最优解是其首要问题。

接下来展开说，何为精英策略。
将第$t$代的产生的新种群$Q_t$与父代$P_t$结合组成合集$R_t$，种族大小为$2N$。然后对$R_t$进行非支配排序，产生非支配集$Z_i$并计算拥挤度。如图所示$Z_1$为$R_t$中最优的，直接放入新的父代种群$P_{t+1}$中；如果新种群$P_{t+1}$个体数量小于$N$则继续填充下一级非支配集$Z_2$。如图所示，$Z_3$如果加入则超出$N$，此时对$Z_3$中的个体使用拥挤度比较算子，使得$P_{t+1}$个体数量达到N然后通过遗传算子(选择、交叉、变异)产生新的子代种群$Q_{t+1}$。

上流程图~

算法伪码

快速非支配排序算法

$n_p$：解$p$被几个解所支配，是一个数值（左下部分点的个数）
$S_p$：解$p$支配哪些解，是一个解集合（右上部分点的内容）
$F_i$：第$i$层的解集合（显然同一层内的解互相都是非支配解）
$p_{rank}$：解$p$所在层的等级（越小越好）

fast-non-dominated-sort($P$)
for each $p\in P$
    $S_p=\varnothing$
    $n_p=0$
    for each $q\in P$
        if ($p<q$) then                if $p$ dominates $q$
             S p = S p ∪ { q } S_p=S_p\cup\{q\} Sp​=Sp​∪{q}        Add $q$ to the set of solutions dominated by $p$
        else if ($q<p$) then
            $n_p=n_p+1$            Increment the domination counter of $p$
    if $n_p=0$ then                    $p$ belongs to the first front
        $p_{rank}=1$
         F 1 = F 1 ∪ { p } F_1=F_1\cup\{p\} F1​=F1​∪{p}
$i=1$                                      initialize the front counter
while $F_i\ne \varnothing$
    $Q=\varnothing$                                Used to store the members of the next front
    for each $p\in F_i$
        for each $q\in F_i$
            $n_q=n_q-1$
            if $n_q=0$ then            $q$ belongs to the next front
                $q_{rank}=i+1$
                 Q = Q ∪ { q } Q=Q\cup\{q\} Q=Q∪{q}
$i=i+1$
$F_i=Q$

拥挤度距离计算

crowding-distance-assignment($I$)
$l=|I|$                                                    number of solutions in $I$
for each $i$, set $I[i{]}_{distance}=0$              initialize distance
for each objective $m$
    $I=$sort$(I,m)$                                   sort using each objective value
    $I[1{]}_{distance}=I[l{]}_{distance}=\infty$      so that boundary points are always selected
    for $i=2$ to $(l-1)$                           for all other points
        $I[i{]}_{distance}=I[i{]}_{distance}+(I[i+1]\ast m-I[i-1]\ast m)/(f_m^{max}-f_m^{min})$

Python代码实现

"""
优化目标：
    min(f1(x), f2(x))
        f1(x) = -x^2
        f2(X) = -(x-2)^2
    s.t x~[-55, 55]
pop_size = 20
max_gen  =  921
"""
import math
import random
import matplotlib.pyplot as plt


def function1(x):
    """
    First function to optimize
    :param x: 自变量
    :returns: 函数计算结果
    """
    return -x ** 2


def function2(x):
    """
    Second function to optimize
    :param x: 自变量
    :returns: 函数计算结果
    """
    return -(x - 2) ** 2


def index_of(a, list):
    """
    Function to find index of list
    :param a: 待查找的目标值
    :param list: 列表
    :returns: 列表中对于值的索引值，不存在返回-1
    """
    for i in range(0, len(list)):
        if list[i] == a:
            return i
    return -1


def sort_by_values(list1, values):
    """
    Function to sort by values
    :param list1:
    :param values:
    :returns:
    """
    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.inf
    return sorted_list


def fast_non_dominated_sort(values1, values2):
    """
    Function to carry out NSGA-II's fast non-dominated sort
    :param values1: 种群1
    :param values2: 种群2
    :return: 该代的种群集合
    """
    S = [[] for i in range(0, len(values1))]  # 支配解
    front = [[]]  # 集合
    n = [0 for i in range(0, len(values1))]   # 支配个体数
    rank = [0 for i in range(0, len(values1))]  # 支配层级

    for p in range(0, len(values1)):
        S[p] = []
        n[p] = 0
        for 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)  # p支配q, 把q加入p的支配解中
            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] + 1  # 支配解个数+1
        if n[p] == 0:
            rank[p] = 0
            if p not in front[0]:
                front[0].append(p)

    i = 0  # 初始化种群个数
    while front[i]:
        Q = []
        for p in front[i]:
            for q in S[p]:
                n[q] = n[q] - 1  # 支配解个数-1
                if n[q] == 0:
                    rank[q] = i + 1
                    if q not in Q:
                        Q.append(q)
        i = i + 1
        front.append(Q)

    del front[len(front) - 1]
    return front


def crowding_distance(values1, values2, front):
    """
    Function to calculate crowding distance
    :param values1: 种群1
    :param values2: 种群2
    :param front:
    :return:
    """
    distance = [0 for i in range(0, len(front))]
    sorted1 = sort_by_values(front, values1[:])
    sorted2 = sort_by_values(front, values2[:])
    distance[0] = 4444444444444444  # infinite
    distance[len(front) - 1] = 4444444444444444  # infinite
    for k in range(1, len(front) - 1):
        distance[k] = distance[k] + (values1[sorted1[k + 1]] - values2[sorted1[k - 1]]) / (max(values1) - min(values1))
    for k in range(1, len(front) - 1):
        distance[k] = distance[k] + (values1[sorted2[k + 1]] - values2[sorted2[k - 1]]) / (max(values2) - min(values2))
    return distance


def crossover(a, b):
    """
    Function to carry out the crossover
    :param a:
    :param b:
    :return:
    """
    r = random.random()
    if r > 0.5:
        return mutation((a + b) / 2)
    else:
        return mutation((a - b) / 2)


def mutation(solution):
    """
    Function to carry out the mutation operator
    :param solution:
    :return:
    """
    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 = 20  # 种群大小
max_gen = 1821

# Initialization
min_x = -55
max_x = 55
# 随机生成初始种群
solution = [min_x + (max_x - min_x) * random.random()for i in range(0, pop_size)]
gen_no = 0

while gen_no < max_gen:  # 达到最大值停止
    # 自适应度计算
    function1_values = [function1(solution[i]) for i in range(0, pop_size)]  # 根据函数生成 种群数 数量大小的值, 即约束
    function2_values = [function2(solution[i]) for i in range(0, pop_size)]  # 根据函数生成 种群数 数量大小的值, 即约束

    # pareto等级
    non_dominated_sorted_solution = fast_non_dominated_sort(function1_values[:], function2_values[:])
    print("The best front for Generation number ", gen_no, " is")
    for val in non_dominated_sorted_solution[0]:
        print(round(solution[val], 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][:]))
    solution2 = solution[:]  # P+Q
    # Generating offsprings
    while 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]))
    # 计算 P+Q种群的适应度
    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[:])
    # 拥挤度距离计算
    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][:]))
    # 得到下一代种群P1
    new_solution = []  # index
    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]))]
        front22 = sort_by_values(non_dominated_sorted_solution2_1[:], crowding_distance_values2[i][:])
        front = [non_dominated_sorted_solution2[i][front22[j]] for j in
                 range(0, len(non_dominated_sorted_solution2[i]))]
        front.reverse()
        for value in front:
            new_solution.append(value)
            if len(new_solution) == pop_size:
                break
        if len(new_solution) == pop_size:
            break
    solution = [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()