原理

与单目标（遗传算法）最大的不同就是进行选择操作之前进行快速非支配排序，这一步也是为了选择操作而来的，选择哪些、怎么选是通过非快速支配排序来的。这就不像单目标，挑好的选就行了

支配： 在多目标优化问题中，如果个体p至少有一个目标比个体q好，而且个体p中的所有目标都不比个体q差，那么称个体p支配个体q。

序值： 如果p支配q，那么p的序值比q低。如果p和q互不支配，那么p和q有相同的序值。

拥挤距离：用来计算某前端中的某个体与该前端中其他个体之间的距离，用以表征个体间的拥挤程度。希望pareto解出来之后，点与点之间距离是相近的，不要太多的聚集在某个地方。用某个点与前后两个点之间的xy的距离和表示。算法会选择拥挤距离大的去领头。

快速非支配排序：快速非支配排序就是将解集分解为不同次序的Pareto前沿的过程。将一组解分成n个集合：rank1,rank2…rankn,每个集合中所有的解都互不支配，但ranki中的任意解支配rankj中的任意解（i<j）.

代码

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  # 解决保存图像是负号'-'显示为方块的问题

pop_size =400
max_gen = 200

#系数
cost_weight=np.array([0.077,0.077,0.077,0.089,0.089,0.089,0.087,0.087,0.087,0.089,0.101,0.099])
harm_weight=np.array([0.082,0.082,0.082,0.079,0.079,0.079,0.069,0.069,0.069,0.079,0.1,0.1])
fault_fre_weight=np.array([0.085,0.085,0.085,0.08,0.08,0.08,0.087,0.087,0.087,0.079,0.079,0.086])

#==========定义两个目标函数=============
# 定义函数1
def function1(solution):

    value = np.dot(cost_weight,(1/(solution**2)))
    return value
# 定义函数2 可靠性
def function2(solution):
    solution = abs(solution)
    value = np.dot(1/(harm_weight*(1-fault_fre_weight)),solution)
    return value

#误差x的范围
accuracy_min=np.array([-9.7,-9.8,-10.2,-6.2,-6.8,-5.5,-6.6,-5.3,-4.8,-0.0247,-0.026,-0.0254])#12个
accuracy_max=np.array([20.3,20.2,19.8,13.8,13.2,14.5,13.4,14.7,15.2,0.0553,0.0540,0.0546])

# solution=[0 for i in range(12)]
#生成解集   相当于x 单位um
def create_x(pop_size):
    solution = [0 for i in range(12)]
    for i in range(12):
        solution[i]=np.random.uniform(accuracy_min[i],accuracy_max[i],pop_size)
    solution=np.array(solution)
    return solution

#求解集对应的函数值
def creat_y(solution):
    values1=function1(solution)
    values2=function2(solution)
    values = [values1, values2]  # 解集【目标函数1解集，目标函数2解集...】 这是列表[ [values1] [values2] ]
    return values

# #这里的标号1 就是y1这个解 也就是生成的第一列的x带入方程得到的y
# plt.scatter(values1,values2, s=20, marker='o')#绘制散点图 s表示点的大小 marker点的样式
# 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;
  若为空，则表明得到了所有非支配集合，程序结束
'''

#这个values1或者2就是一行的y（共10个）

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))]#每个个体被支配解的个数 。即针对每个解，存放有多少好于这个解的个数 这是n中有len(values11)个0 ，len(values11)得出里面有多少数
    rank = [np.inf for i in range(0, len(values11))]#存放每个个体的级别 INF：Infinity，代表的是无穷大的意思，也是属于浮点类型。

    for p in range(0, len(values11)):#遍历每一个个体
        # ====得到各个个体 的被支配解个数 和支配解集合====
        S[p] = [] #该个体支配解的集合 。即存放差于该解的解
        n[p] = 0  #该个体被支配的解的个数初始化为0  即找到有多少好于该解的 解的个数

        # 将每一个个体与p进行比较
        # 如果在两个函数中都比p好 那么比p好的个体数+1n[p] = n[p] + 1
        # 如果在两个函数中都比p差的个体解序号 S[p].append(q)
        for q in range(0, len(values11)):
            less = 0 #函数值小于p的数目
            equal = 0 #的目标函数值等于p个体的目标函数值数目
            greater = 0 #的目标函数值大于p个体的目标函数值数目
            for k in range(len(values)):  # 遍历每一个目标函数 注意：这里len(values)为2
                #注意：下面这个对比是越小越好  要求的两个目标函数的值p，与其他的值对比
                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 差
            #对对比结果进行处理
            #如果比p好 那么比p好的个体数+1
            # if (less + equal == len(values)) and (equal != len(values)):
            if less == len(values):
                n[p] = n[p] + 1  # q比p,  比p好的个体个数加1

            #如果q比p差，存放比p差的个体解序号
            # elif (greater + equal == len(values)) and (equal != len(values)):
            elif greater  == len(values):
                S[p].append(q)  # q比p差，存放比p差的个体解序号

        #=====找出Pareto 最优解，即n[p]===0 的 个体p序号。=====
        if n[p]==0: #n[p]==0表示没有出现比p好的
            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层各个解 这个解就是x
#         jie.append(solution[j])
#     print('第%d层,解为%s'%(i,jie))


# 拥挤度计算
#1、取单个前沿中个体按照一个目标上的值从小到大排序
#2、将最大目标值作为max,最小目标值保留作为min。并且这两个极值点的拥挤距离都被设置为inf即无穷大。因此注意，一个层中可能有多个具有inf的点，
#   即如果层中有多个点在至少一个目标上相等，并且最大或最小，那么这些点的拥挤距离都是无穷大! !因为目标上呈现垂直的关系也是属于非支配的关系! !
#   如果出现这种情况，说明你算法的多样性很烂! ~或者在某些算法早期可能出现这种情况
# 3、在这个目标上计算每个个体最相邻个体之间的距离，即i-1和i+1的目标值的差。 并使用max和min对次值进行归一化。
# 4、遍历目标，将目标上已经归一化的拥挤距离相加。
# 5、进入下一层front前沿
# 6、拥挤距离越大越好，最后按照拥挤距离重新排序各层，进而排序种群。
def crowding_distance(values,front):
    """
    :param values: 群体[目标函数值1，目标函数值2,...]
    :param front: 群体解的等级，类似[[1], [9], [0, 8], [7, 6], [3, 5], [2, 4]] 不同序号就是装载不同的序值的pareto最优解
    :return: front 对应的 拥挤距离
    """
    distance = np.zeros(shape=(pop_size, ))  # 拥挤距离初始化为0 1行pop_size（群体 10个）个的数组
    for rank in front:  # 遍历每一层Pareto解 这个rank装载的是每一次pareto解对应的序号
        for i in range(len(values)):  # 遍历每一层函数值（先遍历群体函数值1，再遍历群体函数值2...） len(values)为2
            valuesi = [values[i][A] for A in rank]  # 函数值
            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  [-1]就代表最后一个

            #计算rank等级中，除去最优解、最差解外。其余解的拥挤距离
            #算式中另外加上了distance[sort_ranki[j]]，因为这个是循环，距离是两个轴之间的距离和
            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]) #[ [存放的是第一序列的解对应的拥挤度] [] [] ]


#精英选择策略：在得到的父代子代中，选择优秀的群体进行下一代
#1、首先将父代种群C和子代种群D合成一个新种群R
#2、根据以下规则从种群R中，生成一个新的父代种群：首先是将pareto将好的那几层都放进去，直到某一层不能全部放进的话，就按照拥挤度选择

#参数distance就是上面的distanceA
#返回的是 一开始创建的我们需要用到的x的解集的那种结构 就是上面的solution挑选了几个
def elitism(pop_size,front, distance,solution ):
        """
        精英选择策略
        :param front: 父代与子代 组合构成的解的等级
        :param distance:  父代与子代 组合构成的解 拥挤距离
        :param solution:  父代与子代 组合构成的解
        :return:  返回群体解。群体数量=（父代+子代）//2
        """
        X1index = []  # 存储群体编号 他的格式的[1,6,3,5]这样
        pop_size = pop_size // 2  # 保留的群体个数 即（父辈+子辈)//2
        X1=np.array([ [],[],[] ])

        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  [1,2,3,4,5这样的格式]
            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]) #这样的写法，追加后的格式是一维的比如[1,2,3,4,5]这样

            #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]
        X1=np.array([ [],[],[],[],[],[],[],[],[],[] ,[],[]])
        for i in X1index:
            X1=np.c_[X1,solution[:,i]]

        return X1


def creat_new_pop(x1,pop_size):
    solution=create_x(pop_size//2)
    for i in range(0,pop_size//2):
        x1=np.c_[x1,solution[:,i]]
    return x1


print('开始：\n')
b=np.array([[] for i in range(12)])

solution=create_x(pop_size)
for i in range(0,max_gen):
    values=creat_y(solution)
    # plt.scatter(values[0], values[1])
    front=fast_non_dominated_sort(values)
    distance=crowding_distance(values,front)
    x1=elitism(pop_size,front,distance,solution)
    if i==max_gen-1:
        for i in front[0]:
            b=np.c_[b,solution[:,i]]
        values=creat_y(b)
        plt.scatter(values[0], values[1],c='black')
        print(front)
        print(b)

        # for i in front[1]:
        #     b=np.c_[b,solution[:,i]]
        # values=creat_y(b)
        # plt.scatter(values[0], values[1],c='black')
        # for i in front[2]:
        #     b=np.c_[b,solution[:,i]]
        # values=creat_y(b)
        # plt.scatter(values[0], values[1],c='black')
        # for i in front[3]:
        #     b=np.c_[b,solution[:,i]]
        # values=creat_y(b)
        # plt.scatter(values[0], values[1],c='red')

    solution=creat_new_pop(x1,pop_size)

print(front[0])
print(front[1])
print(front[2])
plt.xlim(xmax=300,xmin=50)
plt.xlabel('function1')
plt.ylabel('function2')
plt.title('解的分布示意图')
plt.show()

参考

视频，9.34-45
十分钟了解多目标优化算法
代码实现