遗传算法简介
遗传算法（ Genetic Algorithm ）是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。
遗传算法是从代表问题可能潜在的解集的一个种群（ population ）开始的，而一个种群则由经过基因（ gene ）编码的一定数目的个体（ individual ）组成。每个个体实际上是染色体（ chromosome ）带有特征的实体。染色体作为遗传物质的主要载体，即多个基因的集合，其内部表现（即基因型）是某种基因组合，它决定了个体的形状的外部表现，如黑头发的特征是由染色体中控制这一特征的某种基因组合决定的。
因此，在一开始需要实现从表现型到基因型的映射即编码工作。由于仿照基因编码的工作很复杂，我们往往进行简化，如二进制编码，初代种群产生之后，按照适者生存和优胜劣汰的原理，逐代（ generation ）演化产生出越来越好的近似解，在每一代，根据问题域中个体的适应度（ fitness ）大小选择（ selection ）个体，并借助于自然遗传学的遗传算子进行组合交叉（ crossover ）和变异（ mutation ），产生出代表新的解集的种群。这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境，末代种群中的最优个体经过解码（ decoding ），可以作为问题近似最优解。
遗传算法的流程图如下图所示：
编码实现
初始化
首先导入依赖：
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython import display
算法的初始化部分需要完成以下工作：
- 检验参数初始化设定是否合理：每个变量的取值范围设定是否合理（下限<上限）、每个变量 v 的初始值是否在取值范围之内 (下限 < v < 上限)；
- 种群初始化，并对种群进行编码。编码的具体过程我们留在 encode 函数中实现，encode的实现细节在下一小节介绍。
class GA():
# nums: m * n, m is the size of population, n is the num of variables
# bound: n * 2, [[min, max], [min, max], [min, max],...]
# func: target function
# DNA_SIZE: length of DNA encoding
# cross_rate: probability of cross
# mutation_rate: probability of mutation
# iter_time: times of iteration
def __init__(self, nums, bound, func, DNA_SIZE=None, cross_rate=0.8, mutation=0.05, iter_time=100):
nums = np.array(nums)
self.bound = bound = np.array(bound)
# check if each variable has a constraint interval.
if nums.shape[1] != bound.shape[0]:
raise Exception('number of bounds is inconsistent with variables, you have {} variables, {} bounds.'.format(nums.shape[1], bound.shape[0]))
# check if each constraint interval is legal.
for min_bound, max_bound in bound:
if max_bound < min_bound:
raise Exception('{} - {} is not a reasonable range of variables.'.format(min_bound, max_bound))
# check if the value of each variable is legal, i.e. v_min < v < v_max.
for var in nums:
for index, var_curr in enumerate(var):
if var_curr < bound[index][0] or var_curr > bound[index][1]:
raise Exception('variable is not in the setting range.')
# calculate DNA_SIZE
min_nums, max_nums = np.array(list(zip(*bound)))
self.var_len = var_len = max_nums - min_nums
bits = np.ceil(np.log2(var_len + 1))
# calculate the length of DNA code
if DNA_SIZE == None:
DNA_SIZE = int(np.max(bits))
self.DNA_SIZE = DNA_SIZE
# encode population
self.pop = pop = self.encode(nums)
# backup for reset
self.copy_pop = pop.copy()
# cross rate, mutation rate, iteration times and target function
self.cross_rate = cross_rate
self.mutation = mutation
self.iter_time = iter_time
self.func = func
Python笔记
zip() 函数
| zip() 函数用于将可迭代的对象作为参数，将对象中对应的元素打包成一个个元组，然后返回由这些元组组成的对象。我们可以使用 list() 转换来输出列表。如果各个迭代器的元素个数不一致，则返回列表长度与最短的对象相同，利用 * 号操作符，可以将元组解压为列表。 |
| a = [1, 2, 3] |
| b = [4, 5, 6] |
| c = [7, 8, 9, 10] |
| d = [ [1, 2], [10, 20], [100, 200] ] |
| zip(a, b) |
| < zip at 0x224685c1708 > |
| list(zip(a, b)) |
| [(1, 4), (2, 5), (3, 6)] |
| list(zip(a, b, c)) |
| [(1, 4, 7), (2, 5, 8), (3, 6, 9)] |
| list(zip(*d)) |
| [(1, 10, 100), (2, 20, 200)] |
编码与解码
编码是指将实际值转换为二进制编码值，解码过程与之相反。
由上图很容易推算出二进制编码算法的编码公式和解码公式：
其中 DNA_SIZE 为 DNA 编码长度，这个参数开放出来给用户设定，如果用户未设定 DNA_SIZE，则计算出能表示所有变量取值范围的最小编码长度。需要说明的是， 编码长度设定越长，所表示的数值的精度越高。
举个栗子说明：
假设变量 v 的取值范围是[0, 10]，DNA_SIZE 设定为 4 就可以表示 0—10 的所有整数了（10 二进制编码为 1010）。
按照我们上面推算的编码公式，如果 v = 7.5 则 encode_num = (24 - 1) × 7.5 / (10 - 0) = 11.25 ，11.25 舍去小数位后为 11，对应的二进制编码为 1011。现在我们尝试将1011解码回去：decode_num = 11 × （10 - 0）/ (24 - 1) = 7.33 ，可见，原来的 7.5 经过 4 位 DNA 编码再解码后变成了 7.33，精度损失 0.17。
如果我们用 SIZE 为 6 的DNA进行编码会如何呢？同样，当v=7.5则 encode_num = (26 - 1) × 7.5 / (10 - 0) = 47.25 ，47.25 舍去小数位后为 47，对应的二进制编码为 101111。现在再将 101111 解码回去： decode_num = 47 × （10 - 0）/ (26 - 1) = 7.46 ，可见，原来的 7.5 经过 6 位DNA编码再解码后变成了 7.46 ，精度损失减小为 0.04 。
由此可见，DNA_SIZE 设定长度越大，数值求解精度越高；不过因为 DNA_SIZE 变大了，解空间也相应增大，搜索空间的增大导致求解过程的耗时也会随之增加。
class GA():
def encode(self, nums):
pop_encode = np.zeros((*nums.shape, self.DNA_SIZE))
for i in range(nums.shape[0]):
for j in range(nums.shape[1]):
num = int(round((nums[i, j] - self.bound[j][0]) * ((2**self.DNA_SIZE - 1) / self.var_len[j])))
pop_encode[i, j] = [int(k) for k in ('{0:0' + str(self.DNA_SIZE) + 'b}').format(num)]
return pop_encode
def decode(self):
w_vector = np.array([2**i for i in range(self.DNA_SIZE)]).reshape(self.DNA_SIZE, 1)[::-1]
pop_decode = self.pop.dot(w_vector).reshape(self.pop.shape[0:2])
for i in range(pop_decode.shape[0]):
for j in range(pop_decode.shape[1]):
pop_decode[i, j] = pop_decode[i, j] / ((2**self.DNA_SIZE - 1) / self.var_len[j]) + self.bound[j][0]
return pop_decode
python笔记
numpy.zeros 函数创建零数组
| np.zeros((2, 10)) |
| array( [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] ] ) |
| np.zeros((2, 3, 10)) |
| array( [ [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] ], [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] ] ] ) |
format 函数实现进制转换
| ('{0:0' + str(5) + 'b}').format(10) |
| '01010' |
| ('{0:0' + str(10) + 'b}').format(12) |
| '0000001100' |
| [int(k) for k in ('{0:0' + str(10) + 'b}').format(12)] |
| [ 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 ] |
numpy的 点积运算与矩阵相乘
| a = np.array( [ [ 1, 1 ], [ 0, 1 ] ] ) |
| b = np.array( [ [ 2, 0 ], [ 3, 4 ] ] ) |
| a * b # a * b = [ [ 1 * 2, 1 * 0 ], # [ 0 * 3, 1 * 4 ] ] |
| array( [ [ 2, 0 ], [ 0, 4 ] ] ) |
| a.dot(b) # a * b = [ [ 1 * 2 + 1 * 3, 1 * 0 + 1 * 4 ], # [ 0 * 2 + 1 * 3, 0 * 0 + 1 * 4 ] ] |
| array( [ [ 5, 4 ], [ 3, 4 ] ] ) |
| np.dot(a, b) # a * b = [ [ 1 * 2 + 1 * 3, 1 * 0 + 1 * 4 ], # [ 0 * 2 + 1 * 3, 0 * 0 + 1 * 4 ] ] |
| array( [ [ 5, 4 ], [ 3, 4 ] ] ) |
reshape 函数实现array形状转换
| np.array( [ [ 1, 2, 3 ], [ 4, 5, 6 ] ] ) |
| array( [ [ 1, 2, 3 ], [ 4, 5, 6 ] ] ) |
| np.array( [ [ 1, 2, 3 ], [ 4, 5, 6 ] ] ).reshape(3, 2) |
| array( [ [ 1, 2 ], [ 3, 4 ], [ 5, 6 ] ] ) |
| np.array( [ 2**i for i in range(10) ] ) |
| array( [ 1, 2, 4, 8, 16, 32, 64, 128, 256, 512 ] ) |
| np.array( [ 2**i for i in range(10) ] ).reshape(10, 1) |
| array( [ [ 1 ], [ 2 ], [ 4 ], [ 8 ], [ 16 ], [ 32 ], [ 64 ], [ 128 ], [ 256 ], [ 512 ] ] ) |
[::-1] 实现 list 逆序
| np.array( [ 2**i for i in range(10) ] ).reshape(10, 1)[::-1] |
| array( [ [ 512 ], [ 256 ], [128 ], [ 64 ], [ 32 ], [ 16 ], [ 8 ], [ 4 ], [ 2 ], [ 1 ] ] ) |
种群适应度
适应度指的是当代种群对环境的适应能力，对于具体的目标优化函数而言，指的是当前解代入目标函数的映射值。
class GA():
def get_fitness(self):
fitness = self.func(*np.array(list(zip(*self.decode()))))
return fitness
此处为了一次求解出整个种群中所有个体的适应度，需要对解码后的结果稍作转换，因为解码后的结果是 个体 × 变量 形式的，而self.func([1, 2, 3], [10, 20, 30]) 才能一次求解出func(1, 10), func(2, 20), func(3, 30)，即要求输入是 变量 × 个体 形式，这就是代码中 zip(*self.decode()) 所做的工作：
选择
计算完了适应度，就可以对种群进行自然选择了，优胜劣汰适者生存，适应度越高的个体被选择的概率越大。在解实际问题的时候，如果求解目标函数的最小值怎么办呢？那就给目标函数取反，这样求出的最优解就是原目标函数取最大值时的最优解。
class GA():
def select(self):
fitness = self.get_fitness()
fitness_norm = (fitness-np.min(fitness))/(np.max(fitness)-np.min(fitness))
self.pop = self.pop[np.random.choice(np.arange(self.pop.shape[0]), size=self.pop.shape[0], replace=True, p=fitness_norm/np.sum(fitness_norm))]
python笔记
抽样函数 np.random.choice(a, size, replace, p)
| a 为原数组， size 为抽样大小， replace 设置是否为有放回抽样， p 为 a 中各个元素被抽取的概率。需要注意的是 p 应该与 a 等长，且 p 中所有概率之和应该等于1。 |
| np.random.choice( a = np.arange(5), size = 5, replace = True, p = [ 0.1, 0.1, 0.2, 0.3, 0.3 ] ) |
| array( [ 0, 4, 0, 1, 1 ] ) |
交叉
交叉过程模拟的是自然界种群中个体之间的交配行为。对于种群中的每个个体，首先在种群中为其随机选择一个配偶个体，然后再随机选择一些基因点位进行等位交叉。
class GA():
def crossover(self):
for i in range(len(self.pop)):
if np.random.rand() < self.cross_rate:
spouse=np.random.randint(0, self.pop.shape[0], size = 1)
cross_genes=np.random.randint(0, 2, size=(len(self.var_len), self.DNA_SIZE)).astype(bool)
temp = self.pop[i, cross_genes]
self.pop[i, cross_genes] = self.pop[spouse, cross_genes]
self.pop[spouse, cross_genes] = temp
python笔记
np.random.randint(low, high, size) 生成随机整数
| low 代表最小值（包含）， high 代表最大值（不包含），即生成的随机数区间是 [ low, high ) ， size 代表生成随机数的个数。 |
| np.random.randint(0, 2, size = 10) |
| array( [ 1, 0, 0, 0, 1, 0, 1, 0, 0, 0 ] ) |
astype 强制类型转换
| np.array( [ 1, 9, 0 ] ).astype(bool) |
| array([ True, True, False]) |
np.random.rand() 生成 0 - 1 之间的随机小数。
变异
变异模拟的是自然界中的基因突变，基因突变会使种群的多样性进一步增加。基因突变并不一定会造成坏的影响，相反还可能带来意想不到的益处。艾滋病治愈第一人就是拥有 CCR5 基因突变的人。
class GA():
def mutate(self):
for individual in self.pop:
for chromosome in individual:
for gene in range(self.DNA_SIZE):
if np.random.rand() < self.mutation:
chromosome[gene] = 1 if chromosome[gene] == 0 else 1
进化
历经自然选择、交叉、变异的过程即称之为一次进化：
class GA():
def evolution(self):
self.select()
self.crossover()
self.mutate()
获取最优解
经历数代进化，输出最新一代种群中的最优个体作为最优解。当然这不是唯一的方案，在其它参考文献中，有将进化历史中的最优个体保留作为最优解的做法。
class GA():
def search(self):
for _ in range(self.iter_time):
self.evolution()
# self.log()
pop = self.decode()
fitness = self.get_fitness()
optimal_fitness = np.max(fitness)
optimal_individual = pop[(list)(fitness).index(optimal_fitness)]
return optimal_individual
其它
还需要实现算法的 reset 以及 log 打印功能：
class GA():
def reset(self):
self.pop = self.copy_pop.copy()
def log(self):
df_log = pd.DataFrame(np.hstack((self.decode(), self.get_fitness().reshape(len(self.pop), 1))),
columns= [f'x{i}' for i in range(len(self.var_len))] + ['F'])
print(df_log)
python笔记
np.hstack() 函数实现array的拼接。
| a = np.array( [ [ 1, 2, 3 ], [ 10, 20, 30 ] ] ) b = np.array( [ [ 6 ], [ 60 ] ] ) np.hstack((a, b)) |
| array( [ [ 1, 2, 3, 6 ], [ 10, 20, 30, 60 ] ] ) |
由于 self.getfitness() 函数的 shape 是一维的 array（1 × m），需要将其转换成 m × 1 才可以完成与 decode() 结果的拼接操作。
实验结果演示
二维空间实验
目标优化函数：
x ⋅ s i n ( 10 x ) + x ⋅ c o s ( 2 x ) x ·sin(10x) + x·cos(2x) x⋅sin(10x)+x⋅cos(2x) x ∈ ( 0 , 10 ) x \in (0, 10) x∈(0,10)
def example_2D():
func = lambda x : np.sin(10 * x) * x + x * np.cos(2 * x)
nums = [[np.random.rand() * 10] for _ in range(50)]
bound = [(0, 10)]
DNA_SIZE = 10
ga = GA(nums = nums, bound = bound, DNA_SIZE = DNA_SIZE, func = func)
plt.ion()
for _ in range(100):
plt.cla()
x = np.linspace(*bound[0], ga.var_len[0] * 50)
plt.plot(x, func(x))
ga.evolution()
x = ga.decode().reshape(len(ga.pop))
plt.scatter(x, func(x), s = 100, lw = 0, c = 'red', alpha = 0.8)
display.clear_output(wait = True)
display.display(plt.gcf())
plt.close()
三维空间实验
目标优化函数：
e − x 2 − y 2 − e − ( x − 1 ) 2 − ( y − 1 ) 2 e^{-x^2 - y^2} - e^{-(x - 1)^2 - (y - 1)^2} e−x2−y2−e−(x−1)2−(y−1)2 x ∈ ( − 3 , 3 ) , y ∈ ( − 2 , 2 ) x \in (-3, 3), y \in (-2, 2) x∈(−3,3),y∈(−2,2)
def example_3D():
func = lambda x, y: np.exp(-x**2 - y**2) - np.exp(-(x-1)**2 - (y-1)**2)
nums = list(zip(np.arange(-3, 3, 0.125), np.arange(-2, 2, 0.125)))
bound = [(-3, 3), (-2, 2)]
DNA_SIZE = 20
ga = GA(nums = nums, bound = bound, DNA_SIZE = DNA_SIZE, func = func)
plt.ion()
for _ in range(100):
x = np.arange(-3, 3, 0.125)
y = np.arange(-2, 2, 0.125)
x, y = np.meshgrid(x, y)
z = func(x, y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, rstride = 1, cmap = plt.get_cmap('rainbow'))
ax.set_xlabel('x0', color='g')
ax.set_ylabel('x1', color='g')
ax.set_zlabel('F', color='b')
ga.evolution()
ax.scatter(*list(zip(*ga.decode())), ga.get_fitness(), s = 100, lw = 0, c = 'blue', alpha = 0.9)
display.clear_output(wait = True)
display.display(plt.gcf())
plt.close(fig)
参考文献
本文中遗传算法的实现部分参考了crossous的博客，原文已经介绍的非常详尽。我在原文的基础上，补充实现了实验部分，修复了几个 bug，并添加了一些注解，算作我的学习笔记，以备后用。再次向原作者crossous表示感谢。
扫一扫
365
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
