本文共 7088 字，大约阅读时间需要 23 分钟。

前言

  这次做笔记，详解一个简单基于python deap库的一元函数寻优的代码例程。

  安装deap库比较简单，在命令终端输入：

pip3 install deap

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

问题

  问题为一个一元函数在某个取值区间内取极值。

  函数为：   取值范围为[-30,30]。   求该函数在以上范围区间内的最大值。

源码

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

import numpy as npfrom deap import base, tools, creator, algorithmsimport randomfrom scipy.stats import bernoulliGENE_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)#左边界是-30    return 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 npfrom deap import base, tools, creator, algorithmsimport randomfrom 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 = 26toolbox = base.Toolbox()toolbox.register('Binary',bernoulli.rvs,0.5)#创建01随机的序列，0,1的概率各是0.5toolbox.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)#左边界是-30    return 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))