基于DEAP库的python进化算法-7.多目标遗传算法NSGA-II

文章目录

• • • 一.多目标优化简介
    • • 1.多目标优化问题
      • 2.多目标优化求解思路
    • 二.NSGA-II算法解析
    • • 1.快速非支配排序(Fast non-dominated sort)
      • 2.拥挤距离计算(Crowding distance assignment)
      • 3.精英保存策略(Elitism)
    • 三.NSGA-II算法实现
    • • 1.测试函数

一.多目标优化简介

1.多目标优化问题

在很多实际工程问题中，我们的优化目标不止一个，而是对多个目标函数求一个综合最优解。例如在物流配送问题中，不仅要求配送路径最短，还可能需要参与运输车辆最少等。

多目标优化问题通常有如下特点：

• 各目标函数往往相互冲突，无法同时取到最优解。这种冲突使得解的搜索变得更加复杂；
• 各目标函数的单位不同，不能简单比较。因此多目标优化问题的合理解集通常是Pareto最优解。

2.多目标优化求解思路

对于多目标优化问题，传统方法是将原问题通过加权方式变换为单目标优化问题，进而求得最优解。该方法有两大问题：

• 权重的量化设定困难；
• 多目标优化问题的求解结果是包含一组非支配解的解集，用这种方法每次只能求得一个不同解(运气好的话),要求解具有足够多样性的Pareto最优解集非常困难

遗传算法具有多点多方向搜索的特征，在一次搜索中可以得到多个Pareto最优解，因此更适合求解多目标优化问题。

而当前用于求解多目标优化问题的遗传算法一般有两种思路:

• Pareto排序评价方法：其思想是构造基于解的优劣关系排序的适应值函数；在子代中更多保留支配解，从而在迭代一定次数后获得Pareto前沿。最经典的是Glodberg提出的Pareto ranking genetic algorithm，另外我们这里介绍的NSGA-II也是基于Pareto排序评价的方法。
• 多目标函数加权方法：给多目标函数加以各种权重，转化为单目标优化问题；相对于传统方法的固定权重，这类方法通常会在遗传操作前以一定规律生成权重，指导个体的多方向搜索。代表性的算法有Ishibuchi的random weight GA和Gen-Cheng的adaptive weight GA等。

二.NSGA-II算法解析

NSGA-II(nondominated sorting genetic algorithm II)是2002年Deb教授提出的NSGA的改进型，这个算法主要解决了第一版NSGA的三个痛点：

• 非支配排序的高计算复杂度
• 共享参数难以确定
• 缺少保存精英策略

针对这三个问题，在NSGA-II中，Deb提出了快速非支配排序算子，引入了保存精英策略，并用“拥挤距离”(crowding distance)替代了共享(sharing)。

在介绍NSGA-II的整体流程前，我们需要先了解快速非支配排序与拥挤距离的定义。

1.快速非支配排序(Fast non-dominated sort)

解的支配关系与Pareto最优解




快速非支配排序步骤：
快速非支配排序就是将解集分解为不同次序的Pareto前沿的过程。
它可以描述为：

1. 为每个解p分配两个关键量：支配p的解个数n_p以及被p支配的解集S_p；
2. 设置i=1，将n_p=0的个体归入F_i；
3. 对于F_i中的个体，遍历每个解p的S_p，将其中每个解的n_p减1；
4. i =1，将n_p=0的解归入F_i;
5. 重复3、4，直到解集中所有个体都被归入某一个F_i。
   
   DEAP的实现
   DEAP内置了实现快速非支配排序操作的函数

tools.emo.sortNondominated(individuals, k, first_front_only=False)
参数：
individual:被排序的个体列表
k：需要选择的个体数目，注意这里不一定返回正好k个个体，需要看pareto前沿的个体个数
假设设置k=100，而selectedPop   len(Front[i-1]) <100,返回的个体数是
selectedPop len(Front[i-1])   len(Front[i]),应该是一个大于100的值。
first_front_only:如果设置为True，则对次序为1的前沿排序之后就返回

返回：
Pareto前沿的列表

2.拥挤距离计算(Crowding distance assignment)

拥挤距离的定义
在NSGA-II中，为了衡量在同一个前沿中各个解质量的优劣，作者为每个解分配了一个拥挤距离，其背后的思想是让求得的Pareto最优解在objective space中尽量分散，也就有更大可能让解在Pareto最优前沿上均匀分布。

DEAP实现
DEAP中内置了计算拥挤距离的函数

tools.emo.assignCrowdingDist(individuals)

参数：
individual:用来计算拥挤距离的个体列表
返回：
没有返回值，拥挤距离保存在传入的individuals中的每个个体的individual.fitness.crowding_dist属性中

3.精英保存策略(Elitism)

比较操作
根据快速非支配排序和拥挤距离计算的结果，对族群中的个体进行排序：


在每个迭代步的最后，将父代与子代合为一个族群，依照比较操作对合并后族群中的个体进行排序，然后从中选取数量等同于父代规模的优秀子代，这就是NSGA-II算法中的精英保存策略。
DEAP实现
DEAP内置了实现NSGA-II中基于拥挤度的选择函数用来实现精英保存策略

tools.selNSGA2(individuals, k, nd='standard')

参数：
individuals:被选择的个体列表，在NSGA-II中是父代与子代的并集
k：需要选择的个体个数，通常要比被选择的个体数目小；如果与被选择的个体数量相同，那么效果
等同于基于pareto前沿的排序
nd：选择使用的非支配(non-dominated)算法，可以设置为'standard'或'log'
返回：
被选择的个体列表

三.NSGA-II算法实现

1.测试函数

#!usr/bin/env python
# -*- coding:utf-8 _*-
"""
@author: liujie
@software: PyCharm
@file: NSGA-II实现.py
@time: 2020/11/27 16:17
"""
'''
NSGA-II算法实现
'''
import random
import numpy as np
import matplotlib.pyplot as plt
from deap import creator,tools,base,algorithms


# 定义问题
creator.create('MultiObjMin',base.Fitness,weights=(-1.0,-1.0))
creator.create('Individual',list,fitness=creator.MultiObjMin)

# 个体编码
def uniform(low,up):
    # 均匀分布生成个体
    return [random.uniform(a,b) for a,b in zip(low,up)]

pop_size = 100
NDim = 2
# 变量下界
low = [0] * NDim
# 变量上界
up = [1] * NDim
# 生成个体
toolbox = base.Toolbox()
toolbox.register('Attr_float',uniform,low,up)
toolbox.register('Individual',tools.initIterate,creator.Individual,toolbox.Attr_float)
# 生成种群
toolbox.register('Population',tools.initRepeat,list,toolbox.Individual)
pop = toolbox.Population(n = pop_size)

# ZDT3评价函数,ind长度为2
def ZDT3(ind):
    n = len(ind)
    f1 = ind[0]
    g = 1   9 * np.sum(ind[1:]) / (n-1)
    f2 = g * (1 - np.sqrt(ind[0]/g) - ind[0]/g * np.sin(10*np.pi*ind[0]))
    return f1,f2

# 注册工具
toolbox.register('evaluate',ZDT3)
# 锦标赛选择
toolbox.register('selectGen1',tools.selTournament,tournsize=2)
# selTournamentDCD(individuals, k)
toolbox.register('select',tools.selTournamentDCD)
# tools.cxSimulatedBinaryBounded(ind1, ind2, eta, low, up)
toolbox.register('mate',tools.cxSimulatedBinaryBounded,eta=20.0,low=low,up=up)
# 变异 - 多项式变异
toolbox.register('mutate',tools.mutPolynomialBounded,eta=20.0,low=low,up=up,indpb=1.0/NDim)

# 用数据记录算法迭代过程
# 创建统计信息对象
stats = tools.Statistics(key= lambda ind : ind.fitness.values)
stats.register('avg',np.mean)
stats.register('std',np.std)
stats.register('min',np.min)
stats.register('max',np.max)

# 创建日志对象
logbook = tools.Logbook()


# 遗传算法主程序
# 参数设置
maxGen = 50
cxProb = 0.7
mutateProb = 0.2

# 遗传算法迭代
# 第一代
fitnesses = map(toolbox.evaluate,pop)
for ind,fit in zip(pop,fitnesses):
    ind.fitness.values = fit
record = stats.compile(pop)
logbook.record(gen=0,**record)

# 快速非支配排序操作
fronts = tools.emo.sortNondominated(pop,k=pop_size)
# 将每个个体的适应度设置为pareto前沿的次序
for idx,front in enumerate(fronts):
    for ind in front:
        ind.fitness.values = (idx   1),

# 创建子代
offspring = toolbox.selectGen1(pop,pop_size)    # 锦标赛选择
# algorithms.varAnd进化算法的一部分，仅应用变化部分（交叉和变异）,克隆了个体，因此返回的种群独立于输入种群
offspring = algorithms.varAnd(offspring,toolbox,cxProb,mutateProb)

# 第二代之后的迭代
for gen in range(1,maxGen 1):
    # 合并父代与子代
    combinedPop = pop   offspring
    # 评价族群-更新新族群的适应度
    fitnesses = map(toolbox.evaluate,combinedPop)
    for ind,fit in zip(combinedPop,fitnesses):
        ind.fitness.values = fit

    # 快速非支配排序
    fronts = tools.emo.sortNondominated(combinedPop,k=pop_size,first_front_only=False)
    
    # 拥挤距离计算
    for front in fronts:
        tools.emo.assignCrowdingDist(front)

    # 环境选择--精英保留
    pop = []
    for front in fronts:
        pop  = front

    # 复制
    pop = toolbox.clone(pop)
    # 基于拥挤度的选择函数用来实现精英保存策略
    pop = tools.selNSGA2(pop,k=pop_size,nd='standard')

    # 创建子代
    offspring = toolbox.select(pop,pop_size)
    offspring = toolbox.clone(offspring)
    offspring = algorithms.varAnd(offspring,toolbox,cxProb,mutateProb)

    # 记录数据-将stats的注册功能应用于pop，并作为字典返回
    record = stats.compile(pop)
    logbook.record(gen = gen,**record)

# 输出计算过程
logbook.header = 'gen','avg','std','min','max'
print(logbook)

# 输出最优解
bestInd = tools.selBest(pop,1)[0]
bestFit = bestInd.fitness.values
print('当前最优解:',bestInd)
print('对应的函数最小值为:',bestFit)

front = tools.emo.sortNondominated(pop,len(pop))[0]
gridPop
for ind in front:
    plt.plot(ind.fitness.values[0],ind.fitness.values[1],'r.',ms=2)
plt.xlabel('f_1')
plt.ylabel('f_2')
plt.tight_layout()
plt.show()

运行结果

genavg     std     min      max     
0  2.21829 2.4389  0.03415219.1206  
1  1.34439 2.04583 0.03415219.12091 
2  0.9029851.61735 0.02297189.41342 
3  0.6821121.227   -0.3568729.41342 
4  0.7282421.54161 -0.3568729.41342 
5  0.5973911.24451 -0.6174819.33351 
6  0.5117881.05542 -0.6880999.63604 
7  0.3909280.565392-0.6880992.21917 
8  0.3679410.490075-0.6886862.11735 
9  0.3002650.450415-0.70712 1.08968 
10 0.2986290.459386-0.7081151.08968 
11 0.3044690.459981-0.7081151.36692 
12 0.3091140.467134-0.7346221.3744  
13 0.3253220.470062-0.7451011.3744  
14 0.2749980.445207-0.7638960.98771 
15 0.2695560.447157-0.76485 0.98771 
16 0.2585970.447632-0.7719940.98771 
17 0.2765630.435062-0.7719940.98771 
18 0.2710150.422953-0.7732030.9877  
19 0.2693120.44875 -0.7732250.992029
20 0.2630290.434753-0.7733210.992029
21 0.2610870.431257-0.7733210.992029
22 0.2539560.450363-0.7733210.992029
23 0.2629060.433012-0.7733210.997527
24 0.26855 0.431864-0.7733210.99865 
25 0.2777680.43518 -0.7733210.99865 
26 0.2685710.456225-0.7733490.998648
27 0.2628920.447244-0.7733490.998644
28 0.2677430.449131-0.7733490.998644
29 0.2631410.455621-0.7733490.999438
30 0.2580080.459385-0.7733650.999438
31 0.2721480.441205-0.7733650.999438
32 0.2650670.435414-0.7733680.999438
33 0.2635230.436773-0.7733680.999438
34 0.2673210.427655-0.7733680.998816
35 0.2580970.438233-0.7733680.998401
36 0.2583950.442298-0.7733680.998401
37 0.2643260.439118-0.7733680.998378
38 0.2695270.428407-0.7733680.998378
39 0.2603220.446813-0.7733680.998378
40 0.2720010.425237-0.7733680.998378
41 0.2725340.42689 -0.7733680.998372
42 0.2754430.433953-0.7733681.53123 
43 0.2723290.450508-0.7733681.53123 
44 0.2666520.437602-0.7733681.00659 
45 0.2632180.455956-0.7733681.00659 
46 0.2566010.458183-0.7733680.9984  
47 0.2588210.449391-0.7733680.9984  
48 0.2548440.451403-0.7733680.9984  
49 0.2668760.442811-0.7733680.9984  
50 0.2619860.447793-0.7733680.998367
当前最优解: [2.6721132552406934e-06, 2.0752896061476544e-07]
对应的函数最小值为: (2.6721132552406934e-06, 0.998367206028216)

可视化

可以看出NSGA-II算法得到的pareto最优前沿质量很高：最优解均匀分布在不连续前沿的各个线段上；同时在最优前沿以外没有个体存在。

来源：https://www./content-1-766701.html