基于DEAP库的NSGA2代码详解

源代码

完整代码如下所示，评价函数使用的是经典的ZDT3测试函数。

import numpy as np
from deap import base, tools, creator, algorithms
import random
import matplotlib.pyplot as plt
#定义问题
creator.create('MultiObjMin',base.Fitness,weights=(-1.0,-1.0))#两个目标，都求最小值
creator.create('Individual',list,fitness = creator.MultiObjMin)#创建individual类def genInd(low,up):return [random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]#实数编码，一个个体里两个基因toolbox  = base.Toolbox()
Ndim = 2#两个变量
low = [0,0]#两个变量的下界
up = [1,1]#两个变量的上界toolbox.register('genInd',genInd,low,up)
toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.genInd)#注册个体生成工具
#print(toolbox.individual())#打印一个个体检查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)#注册评价函数N_POP = 200#种群内个体数量，参数过小，搜索速度过慢
toolbox.register('population',tools.initRepeat,list,toolbox.individual)#注册种群生成工具
pop = toolbox.population(n=N_POP)#建立种群pop
#for ind in pop:#打印一个种群检查
#   print(ind)toolbox.register('selectGen1', tools.selTournament, tournsize=2)
toolbox.register('select', tools.emo.selTournamentDCD) # 该函数是binary tournament，不需要tournsize
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)NGEN = 200#迭代步数，参数过小，在收敛之前就结束搜索
CXPB = 0.8#交叉概率，参数过小，族群不能有效更新
MUTPB = 0.2#突变概率，参数过小，容易陷入局部最优#开始迭代
#第一代与第二代之后的代数操作不同
fitnesses = toolbox.map(toolbox.evaluate,pop)#为初始种群进行评价，得到适应度
for ind,fitness in zip(pop,fitnesses):ind.fitness.values = fitness#这里先有一个适应度才能进行快速非支配排序fronts = tools.emo.sortNondominated(pop,k=N_POP,first_front_only=False)#快速非支配排序，得到不同前沿的pareto层集合fronts
#print(fronts)#打印一个pareto层集合fronts检查，每层前沿都是一个列表for idx,front in enumerate(fronts):#使用枚举循环得到各层的标号与pareto解#print(idx,front)#打印检查前沿标号与每层的pareto解for ind in front:ind.fitness.values = (idx+1),#将个体的适应度设定为pareto解的前沿次序offspring = toolbox.selectGen1(pop,N_POP)#进行选择
offspring = algorithms.varAnd(offspring,toolbox,CXPB,MUTPB)#只做一次交叉与变异操作#从第二代开始循环
for gen in range(1,NGEN):combinedPop = pop + offspring#将子代与父代结合成一个大种群fitnesses = toolbox.map(toolbox.evaluate,combinedPop)#对该大种群进行适应度计算for ind,fitness in zip(combinedPop,fitnesses):ind.fitness.values = fitnessfronts = tools.emo.sortNondominated(combinedPop,k=N_POP,first_front_only=False)#对该大种群进行快速非支配排序for front in fronts:tools.emo.assignCrowdingDist(front)#计算拥挤距离pop = []for front in fronts:pop += frontpop = toolbox.clone(pop)pop = tools.selNSGA2(pop,k=N_POP,nd='standard')#基于拥挤度实现精英保存策略offspring = toolbox.select(pop,N_POP)#选择offspring = toolbox.clone(offspring)offspring = algorithms.varAnd(offspring,toolbox,CXPB,MUTPB)#交叉变异bestInd = tools.selBest(pop,1)[0]#选择出种群中最优个体
bestFit = bestInd.fitness.values
print('best solution:',bestInd)
print('best fitness:',bestFit)front = tools.emo.sortNondominated(pop,len(pop))[0]#返回的不同前沿的pareto层集合fronts中第一个front为当前最优解集
#图形化显示
for ind in front:plt.plot(ind.fitness.values[0],ind.fitness.values[1],'r.',ms=2)
plt.xlabel('f1')
plt.ylabel('f2')
plt.tight_layout()
plt.show()

问题定义

测试函数为ZTD3，取n=2，也就是有两个目标。两个目标都取最小值，因此代码中目标参数要写成weights=(-1.0,-1.0)。

creator.create('MultiObjMin',base.Fitness,weights=(-1.0,-1.0))#两个目标，都求最小值
creator.create('Individual',list,fitness = creator.MultiObjMin)#创建individual类

个体编码与生成

测试函数ZTD3要求自变量的取值范围是[0,1]，我们有两个变量，说明一个个体中含有两个基因。考虑到基因为复数并且对精度没有特殊需求，因此使用实数编码法生成个体。
以下函数的编写是常见的实数编码法应用，在一个列表中生成两个均匀随机分布的实数。

def genInd(low,up):return [random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]

编码完毕之后设定变量的上界与下界，再注册到工具箱中即可。

def genInd(low,up):return [random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]#实数编码，一个个体里两个基因toolbox  = base.Toolbox()
Ndim = 2#两个变量
low = [0,0]#两个变量的下界
up = [1,1]#两个变量的上界toolbox.register('genInd',genInd,low,up)
toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.genInd)#注册个体生成工具
#print(toolbox.individual())#打印一个个体检查

评价函数（ZTD3）

ZTD3是一个常用来测试代码性能的函数，ZTD3的内容如下图所示：

这里取m=2（无视上图中的m=30，那样目标太多难以使用图像直观观察）。从函数中可以看到，双目标的求解其实是计算出一组解，令其在f1,f2上求得最小值。

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)#注册评价函数

此外，和单目标程序最明显的区别为，评价函数返回了两个值。这样的表现也可以理解为该程序有两个评价函数，分别用于评价两个目标上的适应度。

生成种群

种群的生成方法与单目标程序无异：

N_POP = 200#种群内个体数量，参数过小，搜索速度过慢
toolbox.register('population',tools.initRepeat,list,toolbox.individual)#注册种群生成工具
pop = toolbox.population(n=N_POP)#建立种群pop
#for ind in pop:#打印一个种群检查
#   print(ind)

注册工具与超参数设置

注册工具中注册了两个选择函数，一个是selectGen1，一个是select。selectGen1函数是第一次进行选择时使用的，而第二代之后的选择就需要考虑到拥挤度比较了，因此需要另一个选择函数tools.emo.selTournamentDCD。
此外，在测试了其他的交叉，变异函数后发现运行结果并不理想，说明该多目标算法对交叉，变异的手段有特定的需求，这一点将在今后继续研究。
填坑，后续的项目测试证明，个体采用实数编码，有取值范围的情况，变异与交叉应当采用含边界方式的变异与交叉方法。

toolbox.register('selectGen1', tools.selTournament, tournsize=2)
toolbox.register('select', tools.emo.selTournamentDCD) # 该函数是binary tournament，不需要tournsize
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)NGEN = 200#迭代步数，参数过小，在收敛之前就结束搜索
CXPB = 0.8#交叉概率，参数过小，族群不能有效更新
MUTPB = 0.2#突变概率，参数过小，容易陷入局部最优

第一次迭代

在NSGA2代码中，第一次迭代与后续的迭代操作不同，主要区别有：
（1）第一次迭代不计算拥挤距离。
（2）第一次迭代中选择的依据——适应度，用的是个体对应的pareto解的前沿次序。
如果熟悉之前deap程序的话，可以看到多了些陌生但在NSGA2程序中常见的变量：fronts、front。
首先解释一下fronts，fronts是一个列表，列表内记录了一个种群中个体按照非支配排序后形成的不同层的前沿pareto解。由下图所示，完整的整个列表就是fronts。红框处框选了一个列表，该列表正是第一层pareto解，也就是pareto最优解，其中又包含了两个小列表，每个小列表显然是一个解。
也就是说，该fronts里面的第一层pareto解里面有两个解，同理可以看出，第二层里有4个解，第三层里有4个解。

front的概念就更好理解了，front就是fronts中某一层pareto解的集合，以上图举例，红框框选的那个列表就是一个front，该front里面有两个个体（ind）。
在程序中按照枚举循环的方式打印fronts可以看到内容如下图所示（这里换了一个种群）：

前方箭头指向的数字0~5就是pareto解的层数，该种群中总共有6层pareto解。数字之后的列表就是每个front的内容，可以看到，第一层front有3个解（个体），第二层front有5个解（个体）等等……

#开始迭代
#第一代与第二代之后的代数操作不同
fitnesses = toolbox.map(toolbox.evaluate,pop)#为初始种群进行评价，得到适应度
for ind,fitness in zip(pop,fitnesses):ind.fitness.values = fitness#这里先有一个适应度才能进行快速非支配排序fronts = tools.emo.sortNondominated(pop,k=N_POP,first_front_only=False)#快速非支配排序，得到不同前沿的pareto层集合fronts
#print(fronts)#打印一个pareto层集合fronts检查，每层前沿都是一个列表for idx,front in enumerate(fronts):#使用枚举循环得到各层的标号与pareto解#print(idx,front)#打印检查前沿标号与每层的pareto解for ind in front:ind.fitness.values = (idx+1),#将个体的适应度设定为pareto解的前沿次序offspring = toolbox.selectGen1(pop,N_POP)#进行选择
offspring = algorithms.varAnd(offspring,toolbox,CXPB,MUTPB)#只做一次交叉与变异操作

第二次及之后的迭代

在NSGA2代码中，后续的迭代与第一次迭代操作不同，主要区别有：
（1）引入了父代子代合并后选择的精英选择机制。
（2）计算拥挤距离，利用拥挤距离进行种群内选择。

#从第二代开始循环
for gen in range(1,NGEN):combinedPop = pop + offspring#将子代与父代结合成一个大种群fitnesses = toolbox.map(toolbox.evaluate,combinedPop)#对该大种群进行适应度计算for ind,fitness in zip(combinedPop,fitnesses):ind.fitness.values = fitnessfronts = tools.emo.sortNondominated(combinedPop,k=N_POP,first_front_only=False)#对该大种群进行快速非支配排序for front in fronts:tools.emo.assignCrowdingDist(front)#计算拥挤距离pop = []for front in fronts:pop += frontpop = toolbox.clone(pop)pop = tools.selNSGA2(pop,k=N_POP,nd='standard')#基于拥挤度实现精英保存策略offspring = toolbox.select(pop,N_POP)#选择offspring = toolbox.clone(offspring)offspring = algorithms.varAnd(offspring,toolbox,CXPB,MUTPB)#交叉变异

图形化显示

front = tools.emo.sortNondominated(pop,len(pop))[0]#返回的不同前沿的pareto层集合fronts中第一个front为当前最优解集
#图形化显示
for ind in front:plt.plot(ind.fitness.values[0],ind.fitness.values[1],'r.',ms=2)
plt.xlabel('f1')
plt.ylabel('f2')
plt.tight_layout()
plt.show()

可以看到pareto最优解集将ztd3函数的部分轮廓较为完整地勾勒出来。

后记

最近项目中可用的优化影响因素较少，因此暂时使用单目标优化的程序，更多的约束打算放在惩罚函数中处理，多目标优化的学习要暂时耽搁。