前言
在上一篇中，我们已经介绍了如何在DEAP中实现进化算法的基本操作，这一篇中我们试图将各个操作组装起来，用进化算法解决一个简单的一元函数寻优问题
一.问题描述与分析
给定一个函数，求解其最大值
该函数的图像为：
该函数的最大值应该出现在28.309042309042308处，值为1657.4235763265594。
可以看到该函数有很多局部极值作为干扰项，如果进化算法过早收敛，很容易陷入某个局部最优。
二.问题的编码与解码
对于该问题，可以选取很多不同的编码方式。本文计划采用二进制编码，精确度要求到6位，那么首先应当确定编码的长度与解码方式。由于有：
所以我们需要的二进制编码长度为26
三.利用DEAP自带算法求解
1.DEAP自带的进化算法介绍
DEAP自带的算法都比较基础，通常可以用来测试问题描述、编码方式和交叉突变操作组合的有效性。需要比较复杂的进化算法时，可以通过在已有的算子上进行扩充
简单进化算法
deap.algorithms.eaSimple
DEAP中预置的简单进化算法流程描述如下：
1.根据工具箱中注册的toolbox.evaluate评价族群
2.根据工具箱中注册的toolbox.select选择与父代相同个数的育种个体
3.在族群中进行第一次循环，用工具箱中注册的toolbox.mate进行配种，并用生成的两个子代替换对应父代
4.在族群中进行第二次循环，用工具箱中注册的toolbox.mutate进行变异，用变异后的子代替换对应父代
5.从1开始重复循环，直到达到设定的迭代次数
需要注意的是在这个过程中，生成子代有四种情况：受到配种影响；受到变异影响；既受到配种也受到变异影响；既不受配种影响也不受变异影响
对应的伪代码可以表述为：
evaluate(population)
for g in range(ngen):
population = select(population, len(population))
offspring = varAnd(population, toolbox, cxpb, mutpb)
evaluate(offspring)
population = offspring
(μ + λ)进化算法
deap.algorithms.eaMuPlusLambda
该算法的流程如下：
1.根据工具箱中注册的toolbox.evaluate评价族群
2.在族群中进行循环，在每次循环中，随机选择crossover，mutation和reproduction三者之一：
如果选择到crossover，那么随机选择2个个体，用工具箱中注册的toolbox.mate进行配种，<将生成的第一个子代加入到后代列表中，第二个子代丢弃；
如果选择到mutation，用工具箱中注册的toolbox.mutate进行变异，将变异后的子代加入到后代列表中；
如果选择到reproduction，随机选择一个个体，将其复制加入到后代列表中
3.根据工具箱中注册的toolbox.select，在父代+子代中选择给定数量的个体作为子代
4.从1开始重复循环，直到达到设定的迭代次数
注意在这个子代生成的过程中，子代不会同时受到变异和配种影响。
对应的伪代码可以表述为：
evaluate(population)
for g in range(ngen):
offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
evaluate(offspring)
population = select(population + offspring, mu)
(μ，λ)进化算法
deap.algorithms.eaMuCommaLambda
与(μ + λ)进化算法基本相同，唯一的区别在于生成子代族群时，只在产生的子代中选择，而丢弃所有父代。
对应伪代码可以表述为
evaluate(population)
for g in range(ngen):
offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
evaluate(offspring)
population = select(offspring, mu)
2.调用DEAP自带的进化算法
在调用DEAP自带的算法时需要注意的是，由于内置算法调用的alias已经提前给定，因此我们在register的时候，需要按照给定名称注册。
例如toolbox.register('crossover', tools.cxUniform)就不能被内置算法识别，而应当按照
要求，命名为mate，并且显示给出交叉概率：
toolbox.register('mate', tools.cxUniform, indpb = 0.5)
按照要求，使用预置的算法需要注册的工具有：
toolbox.evaluate：评价函数
toolbox.select：育种选择
toolbox.mate：交叉操作
toolbox.mutate：突变操作
from deap import algorithms, creator, tools, base
from scipy.stats import bernoulli
import numpy as np
import random
random.seed(42) # 保证结果可以复现
# 定义问题
creator.create('FitnessMax', base.Fitness, weights=(1.0,)) # 单目标优化，最大值问题
creator.create('Individual', list, fitness=creator.FitnessMax)
# 生成个体
gene_size = 26 # 26位编码
toolbox = base.Toolbox()
toolbox.register('Binary', bernoulli.rvs, 0.5)
toolbox.register('Individual', tools.initRepeat, creator.Individual, toolbox.Binary, n=gene_size)
# 解码 - 二进制转换为十进制
def decode(individual):
# 解码到10进制
num = int(''.join([str(_) for _ in individual]), 2)
# 映射到-30到30区间
x = -30 + num * 60 / (2 ** 26 - 1)
return x
# 评价函数-适应度
def eval(individual):
# 转化到区间内实数
x = decode(individual)
return ((np.square(x) + x) * np.cos(2 * x) + np.square(x) + x),
# 生成初始族群
pop_size = 100 # 族群中的个体数量
toolbox.register('Population', tools.initRepeat, list, toolbox.Individual)
pop = toolbox.Population(n=pop_size)
# 在工具箱中注册遗传算法需要的工具
toolbox.register('evaluate', eval)
toolbox.register('select', tools.selTournament, tournsize=2) # tournsize=2的锦标赛选择,参数少k
toolbox.register('mate', tools.cxUniform, indpb=0.5) # 均匀交叉，参数少ind1,ind2
toolbox.register('mutate', tools.mutFlipBit, indpb=0.5) # 位翻转变异
# 注册计算过程中需要记录的数据
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)
# 调用DEAP内置算法
resultPop, logbook = algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=50, stats=stats, verbose=False)
# 输出计算过程
# nevals代表该迭代中调用evalute函数的次数
logbook.header = 'gen', 'nevals', 'avg', 'std', 'min', 'max'
print(logbook)
计算过程输出
gen nevals avg std min max
0 100 279.668 376.539 -0.403678 1517.68
1 53 375.377 388.779 0.0422919 1641.3
2 59 335.908 371.105 0.0446177 1516.72
3 52 353.956 324.161 -0.280488 1505.81
4 59 445.318 406.245 0.841933 1528.44
5 62 519.898 415.008 1.39083 1606.41
6 71 499.191 430.979 0.184617 1417.72
7 62 511.865 399.538 0.0560637 1586.01
8 67 491.792 421.77 0.224034 1578.36
9 55 482.469 398.13 1.03128 1426.2
10 59 565.195 445.841 0.0425501 1427.1
11 50 655.547 444.508 0.127139 1547.92
12 49 670.214 466.569 0.0357719 1541.85
13 60 616.596 498.718 0.0158924 1544.06
14 56 573.002 484.627 -0.36737 1544.06
15 61 588.575 473.83 0.00965977 1533.09
16 57 671.817 476.531 0.0386515 1533.09
17 72 692.467 557.545 0.0437906 1639.39
18 66 700.388 527.192 0.0205807 1485.58
19 67 750.856 526.877 0.00740906 1542.36
20 64 719.56 535.471 -0.0687704 1387.21
21 58 691.961 546.773 0.00459343 1469.56
22 60 635.223 525.815 0.0344323 1602.17
23 49 759.276 519.709 0.0329074 1602.31
24 47 928.791 519.078 -0.395298 1628.81
25 57 948.308 530.236 0.0115315 1628.42
26 53 1032.62 496.975 0.0556756 1653.9
27 71 920.054 552.981 -0.24814 1656.97
28 52 1000.62 527.772 0.975663 1656.97
29 61 993.985 556.678 0.0756073 1653.9
30 56 988.46 553.916 2.33859 1653.91
31 68 938.748 598.801 -0.343211 1653.91
32 58 968.058 571.579 -0.374508 1653.91
33 61 1052.52 584.106 -0.395098 1657.17
34 60 1018.82 593.899 0.0344873 1657.17
35 50 1094.81 558.391 0.180522 1657.29
36 40 1149.37 577.618 -0.0701393 1657.29
37 51 1192.23 620.757 -0.0860435 1657.29
38 48 1284.57 593.595 0.434814 1657.29
39 55 1314.9 557.095 5.45569 1657.29
40 58 1279.23 632.164 0.00597406 1657.3
41 54 1339.26 578.89 -0.236449 1657.31
42 53 1329.61 591.209 0.0531222 1657.31
43 60 1319.97 607.591 0.834223 1657.31
44 65 1237 645.968 0.000313783 1657.31
45 71 1284.9 626.286 -0.364857 1657.31
46 51 1392.22 517.275 0.16591 1657.31
47 66 1310.41 592.914 0.304726 1657.31
48 66 1161.6 691.608 0.0934323 1657.31
49 50 1277.25 638.279 0.000197178 1657.31
50 61 1288.16 618.87 0.168657 1657.31
输出结果
# 输出最优解
index = np.argmax([ind.fitness for ind in resultPop])
x = decode(resultPop[index])
print('当前最优解:'+ str(x) + '\t对应的函数值为:'+ str(resultPop[index].fitness))
# 当前最优解:28.294433329916494 对应的函数值为:(1657.069165527869,)
可以看出遗传算法成功避开了局部最优，给出的结果非常接近全局最优解28.309
3.自行编写算法求解
自行编写通常来说会需要比内置函数更长的篇幅，但是也能获得更大的自由度。下面是一个用锦标赛交叉、位翻转突变求解同样问题的例子
import random
import numpy as np
from scipy.stats import bernoulli
from deap import creator,base,tools
random.seed(42)
# 定义问题
creator.create('FitnessMax',base.Fitness,weights=(1.0,))
creator.create('Individual',list,fitness=creator.FitnessMax)
# 生成个体-二进制编码
gene_size = 26
toolbox = base.Toolbox()
toolbox.register('Binary',bernoulli.rvs,0.5)
toolbox.register('Individual',tools.initRepeat,creator.Individual,toolbox.Binary,n=gene_size)
# 生成初始种群
pop_size = 100 # 族群中个体的数量
toolbox.register('Population',tools.initRepeat,list,toolbox.Individual)
pop = toolbox.Population(n=pop_size)
# 评价函数-二进制解码到十进制并计算适应度
def eval(individual):
# 二进制解码到十进制
num = int(''.join([str(_) for _ in individual]),2)
# 转换为区间内的对应实数
x = -30 + 60 * num /(2**26 - 1)
return ((np.square(x) + x) * np.cos(2*x) + np.square(x) + x),
# 注册遗传算法工具-评价函数
toolbox.register('evaluate',eval)
# 评价初始族群
fitnesses = map(toolbox.evaluate,pop)
for ind,fit in zip(pop,fitnesses):
ind.fitness.values = fit
# 进化迭代
N_GEN = 50 # 迭代次数
CXPB = 0.5 # 交叉概率
MUTPB = 0.2 # 突变概率
# 注册进化过程中需要的工具：配种选择，交叉，变异
toolbox.register('select',tools.selTournament,tournsize = 2) # tournsize = 2的锦标赛选择
toolbox.register('mate',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.select(pop,pop_size) # 选择pop_size个个体
# 复制个体，供交叉变异用
selectedInd = list(map(toolbox.clone,selectedTour))
# 对选出的育种族群两两进行交叉，对于被改变的个体，删除其适应度值
# list[start:end:step]
for child1,child2 in zip(selectedInd[::2],selectedInd[1::2]):
if random.random() < CXPB:
toolbox.mate(child1,child2,0.5) # 均匀交叉
del child1.fitness.values
del child2.fitness.values
# 对选出的育种族群进行变异，对于被改变个体，删除适应度值
for mutant in selectedInd:
if random.random() < MUTPB:
toolbox.mutate(mutant,0.5) # 位翻转变异
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','avg','std','min','max'
print(logbook)
计算过程输出
gen avg std min max
0 466.229 468.083 -0.361262 1639.96
1 523.428 500.73 0.000387387 1639.96
2 591.327 499.476 0.761878 1639.96
3 511.309 484.249 -0.404202 1639.96
4 542.566 455.348 -0.367558 1549.18
5 511.607 477.01 0.0601817 1630.98
6 520.779 496.691 0.00812873 1630.98
7 588.768 520.461 -0.298492 1631.01
8 541.561 499.122 -0.202241 1631.01
9 730.499 538.137 0.0242079 1631.01
10 839.125 559.088 0.43618 1631.01
11 829.414 574.283 0.233167 1631.01
12 797.474 641.566 0.0437757 1654.98
13 813.062 648.566 0.999945 1657.27
14 890.102 685.853 -0.154962 1657.28
15 937.362 677.619 0.394182 1657.28
16 1107 656.154 -0.0508697 1657.28
17 1216.31 648.938 1.39068 1657.28
18 1272.64 625.1 -0.270614 1657.28
19 1218.76 612.227 0.0275138 1657.28
20 1082.98 678.12 0.0401458 1657.28
21 1100.7 630.675 0.0618022 1657.28
22 1244.65 595.943 7.8221e-05 1657.28
23 1235.14 641.172 -0.180385 1657.3
24 1219.11 660.604 0.051994 1657.3
25 1400.85 531.75 0.0735181 1657.3
26 1161.86 692.244 -0.22549 1657.3
27 1248.31 642.029 0.150307 1657.3
28 1303.98 619.694 1.40528 1657.3
29 1244.15 652.702 0.0203885 1657.3
30 1199.76 676.044 0.0233112 1657.3
31 1202.59 672.031 -0.136965 1657.3
32 1204.52 667.777 -0.0454565 1657.31
33 1248.02 630.396 0.283428 1657.31
34 1340.9 582.728 14.4446 1657.31
35 1359.91 597.02 -0.378353 1657.31
36 1286.13 640.359 -0.109894 1657.31
37 1310.39 618.029 0.898644 1657.31
38 1308.29 568.212 10.3651 1657.31
39 1283 635.098 0.0838442 1657.31
40 1338.32 584.495 -0.241494 1657.31
41 1334.45 591.224 0.0600735 1657.31
42 1327.41 606.735 1.83592 1657.31
43 1236.79 648.824 0.000148321 1657.31
44 1283.51 627.575 -0.384135 1657.36
45 1393.7 521.303 0.159144 1657.36
46 1310.79 590.335 0.529636 1657.36
47 1161.16 692.163 0.0592846 1657.36
48 1271.36 641.495 0.000130297 1657.36
49 1284.52 625.65 0.213304 1657.36
输出最优解
# 输出最优解
index = np.argmax([ind.fitness for ind in pop])
x = decode(pop[index])
print('当前最优解:'+str(x)+'\t对应的函数值为:'+str(pop[index].fitness))
# 当前最优解:28.308091138423848 对应的函数值为:(1657.4220750846391,)
结果可视化
import matplotlib.pyplot as plt
# 用select方法从logbook中提取迭代次数,最大值、均值
gen = logbook.select('gen')
fit_maxs = logbook.select('max')
fit_average = logbook.select('avg')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(gen,fit_maxs,'b-',linewidth=2.0,label='Max Fitness')
ax.plot(gen,fit_average,'r-',linewidth=2.0,label='Average Fitness')
ax.legend(loc='best')
ax.set_xlabel('Generation')
ax.set_ylabel('Fitness')
plt.tight_layout()
plt.savefig('Generation_Fitness.png')
