python遗传算法（应用篇1）--求解一元函数极值

下面我们使用遗传算法尝试求解一元函数的最值
$$y=sin(x^2-1)+2cos(2x)，x\in [0,10]$$

遗传算法求解过程

算法参数

# 种群数量
m_population = 50
# 迭代次数（种群进化代数）
epoch = 99
# 基因个数
L = 40
# 交叉概率
mp = 0.5
# 变异概率
mm = 0.01

# 最大函数值
mf = 0
# 最大函数值对应的自变量
mx = 0
# 变量约束
bound = [0, 10]
# 存储每一代的最优个体
f_list = []
x_list = []

构建初始化种群

生成二进制数组，形状为（种群个体数，个体基因个数），即(m_population, L)

population = np.random.randint(0, 2, (m_population, L))

运行结果

解码（二进制>十进制）

1. 将生成的二进制数组先转换为十进制数字，便于计算个体适应度
   二进制转十进制有多种方法，这里笔者提供一种：由于个体基因是以二进制数组的形式存在，我们先将二进制数组转换为二进制字符串，然后使用int()函数实现转换过程

for i in range(m_population):
    transform_middle = int(''.join(map(lambda x: str(x), population[i])), 2)

2. 然后进行标准化。由于二进制字符串的长度为40，转化得到的十进制数字将会非常大，我们的自变量约束区域为$bound=[0,10]$，我们将得到的十进制数字都转化到$0-10$之间，

    transform_[i] = bound[0] + (bound[1] - bound[0]) * transform_middle / pow(2, L)	

3. 求不同个体的的适应度

f = fun(transform_)

自然选择

这是事关重要的一步，我们知道适应度大的个体有更大的概率存活下来，如何使用代码实现这一过程？
我们使用频数来代替概率，在原来种群的基础之上，我们通过概率选择一些个体出来，数量等于原来的种群个体数量相等，这样选择出来的个体在组成一个新的种群，这样我们便完成了自然选择这一过程

bingo = np.random.choice(np.arange(m_population), m_population, True, f / sum(f))
select_populaton = population[bingo]

# 生成目前种群的两个复制版本，实际上不生成也行，但是这样有利于模拟两个体的交配过程
last_gen = select_populaton
new_gen = select_populaton.copy()

交叉

交叉相当于两个个体（染色体）交配。遗传算法中有多种方式实现交叉，笔者在这里采用一种叫做单点交叉的方式

match = np.random.choice(np.arange(m_population), m_population, False)
for i in range(m_population // 2):# 种群中的个体两两结合
    if np.random.rand() < mp:# 以交叉概率
        location = np.random.randint(0, L)# 交叉点
        last_gen[match[i]][:location], new_gen[match[L - i - 1]][:location] = \
            new_gen[match[L - i - 1]][:location], last_gen[match[i]][:location]

变异

for i in range(m_population):
    for j in range(L):
        if np.random.rand() < mm:
            new_gen[i][j] = 1 if new_gen[i][j]==0 else 0

此时上一代繁衍基本结束，new_gen即为新一代种群

将population更新为新种群

population = new_gen

解码（新种群>十进制）

transform_1 = np.zeros(m_population)
for i in range(m_population):
    transform_middle = int(''.join(map(lambda x: str(x), new_gen[i])), 2)
    transform_1[i] = bound[0] + (bound[1] - bound[0]) * transform_middle / pow(2, L)

计算新种群的适应度

# 计算适应度
new_f = fun(transform_1) - 4
# 选择出最优的个体
mx = transform_1[np.argmax(new_f)]
# 最优适应度
mf = new_f.max()
# 将每一代最优个体及其适应度添加到数组中便于统计
f_list.append(mf)
x_list.append(mx)

运行结果

完整代码及其可视化版本

可视化，即使用matplotlib库动态绘制遗传算法求解该问题的过程，主要添加的代码如下

1. 绘制该函数图像
2. 以散点图的形式绘制每一代个体的最优个体及其适应度

import numpy as np
import matplotlib.pyplot as plt

# 定义适应度函数
def fun(x):
    return np.sin(x ** 2 - 1) + 2 * np.cos(2 * x) + 4

plt.ion()
inputs=np.arange(0,10,0.1)
output=[fun(i)-4 for i in inputs]
plt.plot(inputs,output)
# 种群数量
m_population = 50
# 迭代次数
epoch = 99
# 基因个数
L = 40
# 交叉概率
mp = 0.5
# 变异概率
mm = 0.01

# 最大函数值
mf = 0
# 最大函数值对应的自变量
mx = 0

bound = [0, 10]

f_list = []
x_list = []
# 生成二进制种群
population = np.random.randint(0, 2, (m_population, L))
transform_ = np.zeros(m_population)
# 二进制>十进制，并标准化
for i in range(m_population):
    transform_middle = int(''.join(map(lambda x: str(x), population[i])), 2)
    transform_[i] = bound[0] + (bound[1] - bound[0]) * transform_middle / pow(2, L)
# 计算当前种群适应度
f = fun(transform_)
for num in range(epoch):
	# 自然选择
    bingo = np.random.choice(np.arange(m_population), m_population, True, f / sum(f))
    select_populaton = population[bingo]
	# 生成目前种群的两个复制版本，实际上不生成也行，但是这样有利于模拟两个体的交配过程
    last_gen = select_populaton
    new_gen = select_populaton.copy()
	# 交叉
    match = np.random.choice(np.arange(m_population), m_population, False)
    for i in range(m_population // 2):
        if np.random.rand() < mp:
            location = np.random.randint(0, L)
            last_gen[match[i]][:location], new_gen[match[L - i - 1]][:location] = \
                new_gen[match[L - i - 1]][:location], last_gen[match[i]][:location]
	# 变异
    for i in range(m_population):
        for j in range(L):
            if np.random.rand() < mm:
                new_gen[i][j] = 1 - new_gen[i][j]
    # 现在new_gen就是新生成的种群，我们需要将其转化为十进制
    # 我们先生成一个空数组，然后使用它储存转化后的种群十进制
    transform_1 = np.zeros(m_population)
    for i in range(m_population):
        transform_middle = int(''.join(map(lambda x: str(x), new_gen[i])), 2)
        transform_1[i] = bound[0] + (bound[1] - bound[0]) * transform_middle / pow(2, L)
    # 计算适应度
    new_f = fun(transform_1) - 4
    # 选择出最优的个体
    mx = transform_1[np.argmax(new_f)]
    # 将当前的种群视为最优的种群
    population = new_gen
    # 最优适应度
    mf = new_f.max()
    f_list.append(mf)
    x_list.append(mx)
    plt.scatter(mx,mf)
    plt.title(f'epoch:{num+1}')
    plt.pause(0.1)
plt.ioff()
plt.show()

print(max(x_list), max(f_list))
>>> 9.93222096441059 2.9685353661257627

运行结果

其他

numpy中的随机数

1. np.random.randint

numpy.random.randint(low, high=None, size=None, dtype=int)

返回$(low,high]$之间的随机整数或形状为size的随机整数数组

2. np.random.choice

numpy.random.choice(a, size=None, replace=True, p=None)

从一维数组a中抽取随机样本

参数		描述
size	int or tuple of ints, optional	随机样本的数组形状
replace	boolean, optional	表示样本是否可以重复出现
p	1-D array-like, optional	每个样本被抽取到的概率