使用deap实现遗传算法

发表于 |

最近在做毕业论文，其中涉及到一项是检验某个为有限元程序优化参数的遗传算法的正确性。
deap是Python上的一个遗传算法库，主要封装了和遗传算法相关的生成population、select、mutant这些相关的操作。比较方便的是deap框架允许自定义个体类型、种群生成方法、评价函数等，所以灵活性比较高

任务简述

可以把有限元程序想象成一个有六个参数的黑箱函数f(X | args)，现在对于数据X计算得结果Y = f(X | args)和实际值T有偏差，现在希望通过参数优化（参数值必须在一定范围内）使得计算值尽可能接近实际值。不考虑泛化能力等问题，

deap框架

前期工作

安装deap，pip install一如既往地失败了，幸亏还可以通过setup.py安装。

creator和register

这两个模块起到语法糖的作用

creator基于给定的类创建一个派生类，生成的类通过creator.ClassName访问

不同于creator，register创建为一个函数创建别名，并且可以绑定其中的一些参数，生成的函数通过tools.func_name访问

1
toolbox.register("select", tools.selTournament, tournsize=3)

表示注册一个tools.select作为tools.selTournament的别名，并绑定tools.selTournament的参数tournsize为3

评价函数

定义评价值类

我们希望用一个scalar或者一组scalar来量化对个体的评价，它的值由评价函数tools.evaluate计算得到。在使用deap时，评价值应当继承自base.Fitness，这个类包含一个values字段，它是一个tuple，实际上是评价值。
下面创建一个评价值类

1
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))

这里的weights表示每个scalar的权重，当weights=(-1.0,)时说明评价值是一个scalar，并且值越小，评价越高。
所以我们看到不直接使用元组，而选择继承封装元组的的Fitness类是为了更方便地比较评价值

评价函数

使用均方误差作为评价函数。整个评价函数的实现流程应当是：

  1. 使用当前参数值调用有限元程序（fortran77）
  2. 获得有限元程序中返回结果
  3. 计算均方误差
    由于整个源程序是5000+行的比较乱的fortran77代码，还要外部link一个obj文件，所以无论是强行转换到C++还是在上面继续实现都比较麻烦，所以使用Python通过管道来调用fortran77程序中的有限元函数，并获得返回结果。

定义初始化个体

一个个体应该包含若干数量的基因，也就是我们要优化的参数，通常可以用一个list来存储，类似于Fitness类的方法，我们不直接使用这个list，而是通过creator来创建一个继承自list的类
creator.create(“Individual”, list, fitness=creator.FitnessMin)
这个类包含有fitness字段，也就是它的评价值
下面我们使用这个类生成若干个体，这可以分解为两步：

  1. individual函数为某个个体的基因提供随机的初始值
    可以使用tools.initIterate函数，这个函数在/tools/init.py中定义如下

    1
    2
    def initIterate(container, generator):
    	return container(generator())

    其中container就是你定义个体的类Individual
    generator是一个生成器来提供初始值。假设一个个体有6个基因，可以参照以下代码

    1
    2
    3
    def f():
    	for gene_index in xrange(6):
    		yield random_value_for_gene_index

  2. population通过调用若干次individual函数生成种群
    方便起见可以使用/tools/init.py钟提供的另一个tools.initRepeat函数：

    1
    2
    def initRepeat(container, func, n):
    	return container(func() for _ in xrange(n))

    这个函数将func执行n次

运行遗传算法

在deap的文档中，可以下载到deap_onemax的源码
这份代码使用了定义的四个函数，分别是evaluate、mate、mutant、select，分别对应计算评价值、交配、突变、选择等四个原子操作
下面给出的是遗传算法的主程序，在修改时需要注意其中的突变函数toolbox.mutate(mutant)，如果不希望突变，应该把这里注释掉，否则应该手动实现一个自己的突变函数，否则容易随机出无效值。特别地，在deap_onemax的源码中，这个函数时是直接01取反，所以新的程序中出现除0错误很可能是这个原因。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
pop = toolbox.population(n=20)
CXPB, MUTPB, NGEN = 0.5, 0.2, 4
# Evaluate the entire population
# fitnesses = list(map(toolbox.evaluate, pop))
fitnesses = []
for (i, ind) in zip(count(0), pop):
	mse = toolbox.evaluate(ind)
	fitnesses.append(mse)
for ind, fit in zip(pop, fitnesses):
	ind.fitness.values = fit	
# print("  Evaluated %i individuals" % len(pop))
print "Start GA Loop"	
for g in range(NGEN):
	print("-- Generation %i --" % g)
	offspring = toolbox.select(pop, len(pop))
	offspring = list(map(toolbox.clone, offspring))
	
	# offspring是个体组成的`list`
	for child1, child2 in zip(offspring[::2], offspring[1::2]):

		if random.random() < CXPB:
			toolbox.mate(child1, child2)
			del child1.fitness.values
			del child2.fitness.values

	for mutant in offspring:
		if random.random() < MUTPB:
			toolbox.mutate(mutant)
			del mutant.fitness.values

	invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
	fitnesses = map(toolbox.evaluate, invalid_ind)
	for ind, fit in zip(invalid_ind, fitnesses):
		ind.fitness.values = fit
	
	# print("  Evaluated %i individuals" % len(invalid_ind))
	
	pop[:] = offspring
	fits = [ind.fitness.values[0] for ind in pop]
	
	length = len(pop)
	mean = sum(fits) / length
	sum2 = sum(x*x for x in fits)
	std = abs(sum2 / length - mean**2)**0.5
	
	print("  Min %s" % min(fits))
	print("  Max %s" % max(fits))
	print("  Avg %s" % mean)
	print("  Std %s" % std)
	best_ind = tools.selBest(pop, 1)[0]

best_ind = tools.selBest(pop, 1)[0]
print("Best individual is %s, %s" % (best_ind, best_ind.fitness.values))
return best_ind