基于遗传算法（deap库）的一元函数寻优代码详解

前言

这次做笔记，详解一个简单基于python deap库的一元函数寻优的代码例程。
安装deap库比较简单，在命令终端输入：

pip3 install deap

deap库内集成了遗传算法的主要运行函数，因此在实际使用的时候一般只需要自己编写个体编码方法、评价函数、惩罚函数即可。掌握遗传算法的原理流程与Deap库的程序框架就比较重要，遗传算法的原理这里不加赘述。

问题

问题为一个一元函数在某个取值区间内取极值。
函数为：
取值范围为[-30,30]。
求该函数在以上范围区间内的最大值。

源码

以下为完整代码，可直接运行。

import numpy as np
from deap import base, tools, creator, algorithms
import random
from scipy.stats import bernoulli
GENE_SIZE = 26
#定义问题
creator.create('FitnessMax',base.Fitness,weights=(1.0,))#单变量，求最大值
creator.create('Individual',list,fitness = creator.FitnessMax)#创建individual类
#定义个体
toolbox = base.Toolbox()
toolbox.register('Binary',bernoulli.rvs,0.5)#01随机
toolbox.register('individual',tools.initRepeat,creator.Individual,toolbox.Binary,n=GENE_SIZE)#注册个体def decode(individual):#解码函数num = int(''.join([str(_) for _ in individual]),2)x = -30 + num * 60 /(2**26 - 1)#左边界是-30return x#返回真实值（十进制）def eval(individual):#评价函数，先解码，后评价适应度x = decode(individual)#先解码return ((np.square(x) + x) * np.cos(2*x) + np.square(x) + x),#返回适应度
toolbox.register('evaluate',eval)#注册评价函数
#定义种群
N_POP = 100#种群内个体数量
toolbox.register('population',tools.initRepeat,list,toolbox.individual)#注册种群
pop = toolbox.population(n=N_POP)#得到初始种群fitnesses = map(toolbox.evaluate,pop)#计算初始种群的适应度，先计算得到和种群对应的适应度阵列
for ind,fit in zip(pop,fitnesses):#使用循环将适应度阵列中的元素赋值给种群中个体的适应度变量ind.fitness.values = fitN_GEN = 100#迭代数量
CXPB = 0.5#交叉概率
MUTPB = 0.2#变异概率toolbox.register('tourSel',tools.selTournament,tournsize = 2)#注册锦标赛工具
toolbox.register('crossover',tools.cxUniform)#注册交叉工具
toolbox.register('mutate',tools.mutFlipBit)#注册变异工具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()for gen in range(N_GEN):#利用代数控制循环selectedTour = toolbox.tourSel(pop,N_POP)#对种群进行选择selectedInd = list(map(toolbox.clone,selectedTour))#复制个体for child1,child2 in zip(selectedInd[0::2],selectedInd[1::2]):#将选择好的个体按照单数双数分成两组if random.random() < CXPB:#满足交叉概率toolbox.crossover(child1,child2,0.5)#将这两组进行交叉del child1.fitness.values#交叉完适应度肯定也变了，所以先删掉发生交叉的个体的适应度del child2.fitness.values#同上，毕竟有两组for mutant in selectedInd:#对选择好的种群遍历if random.random() < MUTPB:#满足变异概率toolbox.mutate(mutant,0.3)#变异del mutant.fitness.values#变异了适应度也变了，先把适应度删掉invalid_ind = [ind for ind in selectedInd if not ind.fitness.valid]#选择出删掉适应度的个体fitnesses = map(toolbox.evaluate,invalid_ind)#得到与种群对应的适应度阵列for ind,fit in zip(invalid_ind,fitnesses):#将适应度阵列中的元素赋值给种群个体的适应度变量 ind.fitness.values = fit#适应度赋值pop[:] = selectedInd#重新赋值给原先种群变量，进入下一轮操作record = stats.compile(pop)logbook.record(gen=gen,**record)logbook.header = 'gen','nevals','avg','std','min','max'
print(logbook)
#显示适应度最大的变量
index = np.argmax([ind.fitness for ind in pop])
x = decode(pop[index])
print('best solution:'+str(x)+' function value:'+str(pop[index].fitness))

源码讲解

导入库

在算法中必须导入的库有：numpy（用来实现数学运算操作，算法必备）、deap（遗传算法核心）、random（随机函数库，随机算法必备）、bernouli（这里使用伯努利模块来随机生成二进制码，给个体编码用）。

import numpy as np
from deap import base, tools, creator, algorithms
import random
from scipy.stats import bernoulli

定义问题

在deap中需要使用create方法来定义问题，这个写法是种定式，该具体问题是求解最大值，因此weights=(1.0,),只有一个参数是因为该问题是单目标优化问题。

#定义问题
creator.create('FitnessMax',base.Fitness,weights=(1.0,))#单变量，求最大值
creator.create('Individual',list,fitness = creator.FitnessMax)#创建individual类

求解最小值应该参照以下写法。这里需要注意，求解评价函数的最大值程序与最小值程序只在此处有区别。因为求最小值的时候将适应度换成了负数，因此在数值上还是在计算适应度最大的情况。

#定义问题
creator.create('FitnessMin',base.Fitness,weights=(-1.0,))#单变量，求最小值
creator.create('Individual',list,fitness = creator.FitnessMin)#创建individual类

定义个体

这次使用二进制码来组成个体（染色体）。使用二进制码比使用实数码更为复杂，因为要计算基因长度与编写解码函数。
首先计算基因长度。变量的取值范围是[-30,30]，因此整个范围的长度是30-（-30）=60 。接下来确定答案要精确的位数，如果设定精确度是6位的话，显然在60的范围中总共会有60*1e6个数，二进制的长度必须要能保证包含60*1e6，因此可以算出2的26次方是刚好比这个值更大，因此选择26作为基因长度。

GENE_SIZE = 26
toolbox = base.Toolbox()
toolbox.register('Binary',bernoulli.rvs,0.5)#创建01随机的序列，0,1的概率各是0.5
toolbox.register('individual',tools.initRepeat,creator.Individual,toolbox.Binary,n=GENE_SIZE)#正式注册个体

此外使用bernouli库来对每个个体的26个基因进行随机的预编码供之后使用。如果需要检查自己编写的个体是否正确，可以使用以下代码打印出来一个个体检查：

print(toolbox.individual())

解码函数与评价函数

我们需要求解的值当然是10进制的，但是个体在进行交叉变异等一系列操作的时候是二进制的，显然我们需要一个解码机制，令需要使用10进制的时候将解从二进制转化到10进制。
解码函数如下所示，num为某个个体的绝对长度，这个长度显然需要换算到坐标轴上来得到具体的值。
评价函数在该问题中就是需要求解的一元函数，显然在函数中的值越大，就越接近极值这个答案。

def decode(individual):#解码函数num = int(''.join([str(_) for _ in individual]),2)x = -30 + num * 60 /(2**26 - 1)#左边界是-30return x#返回真实值（十进制）def eval(individual):#评价函数，先解码，后评价适应度x = decode(individual)#先解码return ((np.square(x) + x) * np.cos(2*x) + np.square(x) + x),#返回适应度
toolbox.register('evaluate',eval)#注册评价函数

创建种群与计算适应度

创建种群的代码也是定式的，只需要设定种群个体数量即可。创建后会得到一个变量pop，也就是该程序中的具体种群。
如果想要检查一下种群是不是自己设计的样子，可以使用以下代码打印一整个种群检查：

for ind in pop:print(ind)

得到一个种群后需要立刻为其计算适应度。这样在正式计算前就有了一个拥有适应度的种群pop。

N_POP = 100#种群内个体数量
toolbox.register('population',tools.initRepeat,list,toolbox.individual)#注册种群
pop = toolbox.population(n=N_POP)#得到初始种群fitnesses = map(toolbox.evaluate,pop)#计算初始种群的适应度，先计算得到和种群对应的适应度阵列
for ind,fit in zip(pop,fitnesses):#使用循环将适应度阵列中的元素赋值给种群中个体的适应度变量ind.fitness.values = fit#适应度赋值

设定超参数与注册方法

在遗传算法的迭代计算中，需要执行的最基本的三种操作为选择、交叉、变异。因此需要为这三种操作注册对应的方法，这里使用的是锦标赛选择法，均匀交叉法，位翻转突变法。
此外，在deap中，如果使用deap的内置遗传算法，那么注册的名称必须采用固定的形式：toolbox.evaluate、toolbox.select、toolbox.mate、toolbox.mutate。这里不使用内置遗传算法，因此注册名称可以自定义。
这里也可以看到deap库的工具注册机制是如何实现的：以注册锦标赛工具为例，注册好的方法名为tourSel，它的本体是方法tools.selTournament()，只是换了一个名字添加到了库的工具箱里面。其后的所有内容其实都是方法tools.selTournament()的传参。

N_GEN = 100#迭代数量
CXPB = 0.5#交叉概率
MUTPB = 0.2#变异概率toolbox.register('tourSel',tools.selTournament,tournsize = 2)#注册锦标赛工具
toolbox.register('crossover',tools.cxUniform)#注册交叉工具
toolbox.register('mutate',tools.mutFlipBit)#注册变异工具

迭代循环

使用代数来控制循环，因为在不清楚该问题的适应度情况下，贸然使用适应度作为阈值控制循环并不合适。
在迭代循环中依次实现：选择，交叉，变异。这里需要注意，在交叉，变异之后，该个体的适应度通常会产生变化（自身构造变了），因此需要将原先的适应度删掉，再将这些个体重新赋予适应度。

for gen in range(N_GEN):#利用代数控制循环selectedTour = toolbox.tourSel(pop,N_POP)#对种群进行选择selectedInd = list(map(toolbox.clone,selectedTour))#复制个体for child1,child2 in zip(selectedInd[0::2],selectedInd[1::2]):#将选择好的个体按照单数双数分成两组if random.random() < CXPB:#满足交叉概率toolbox.crossover(child1,child2,0.5)#将这两组进行交叉del child1.fitness.values#交叉完适应度肯定也变了，所以先删掉发生交叉的个体的适应度del child2.fitness.values#同上，毕竟有两组for mutant in selectedInd:#对选择好的种群遍历if random.random() < MUTPB:#满足变异概率toolbox.mutate(mutant,0.3)#变异del mutant.fitness.values#变异了适应度也变了，先把适应度删掉invalid_ind = [ind for ind in selectedInd if not ind.fitness.valid]#选择出删掉适应度的个体fitnesses = map(toolbox.evaluate,invalid_ind)#得到与种群对应的适应度阵列for ind,fit in zip(invalid_ind,fitnesses):#将适应度阵列中的元素赋值给种群个体的适应度变量 ind.fitness.values = fit#适应度赋值pop[:] = selectedInd#重新赋值给原先种群变量，进入下一轮操作

最后计算完成，循环找出种群pop中适应度最大的那个个体，解码，并得到其值，就是这一次运算的解。

#显示适应度最大的变量
index = np.argmax([ind.fitness for ind in pop])
x = decode(pop[index])
print('best solution:'+str(x)+' function value:'+str(pop[index].fitness))