一、问题描述
利用遗传算法求解一些典型的二元单目标函数优化问题,对五个二元最优化函数(函数表达式、决策变量取值范围)进行求解,结果要尽可能精确。
五个函数分别为:
二维球形函数: f 1 ( x , y ) = x 2 + y 2 f _ { 1 } ( x , y ) = x ^ { 2 } + y ^ { 2 } f1(x,y)=x2+y2, x , y ∈ [ − 5.12 , 5.12 ] x , y \in [ - 5.12,5.12 ] x,y∈[−5.12,5.12]
De-Jong函数: f 2 ( x , y ) = 100 ( y − x 2 ) 2 + ( x − 1 ) 2 f _ { 2 } ( x , y ) = 100 \left( y - x ^ { 2 } \right) ^ { 2 } + ( x - 1 ) ^ { 2 } f2(x,y)=100(y−x2)2+(x−1)2, x , y ∈ [ − 2.048 , 2.048 ] x , y \in [ - 2.048,2.048 ] x,y∈[−2.048,2.048]
Himmelblau函数: f 3 ( x , y ) = ( x 2 + y − 11 ) 2 + ( x + y 2 − 7 ) 2 f _ { 3 } ( x , y ) = \left( x ^ { 2 } + y - 11 \right) ^ { 2 } + \left( x + y ^ { 2 } - 7 \right) ^ { 2 } f3(x,y)=(x2+y−11)2+(x+y2−7)2, x , y ∈ [ − 6 , 6 ] x , y \in [ - 6,6 ] x,y∈[−6,6]
SIX-HUMP CAMEL函数: f 4 ( x , y ) = 4 x 2 − 2.1 x 4 + x 6 / 3 + x y + 4 y 4 − 4 y 2 f _ { 4 } ( x , y ) =4 x ^ 2 - 2.1 x ^ 4 + x ^6 / 3 + xy + 4 y ^4 - 4 y^ 2 f4(x,y)=4x2−2.1x4+x6/3+xy+4y4−4y2, x , y ∈ [ − 5 , 5 ] x , y \in [ - 5,5 ] x,y∈[−5,5]
BOHACHEVSKY函数: f 5 ( x , y ) = x 2 + 2 y 2 − 0.3 cos 3 π x cos 4 π y + 0.3 f _ { 5 } ( x , y ) = x ^ { 2 } + 2 y ^ { 2 } - 0.3 \cos 3 \pi x \cos 4 \pi y + 0.3 f5(x,y)=x2+2y2−0.3cos3πxcos4πy+0.3, x , y ∈ [ − 1 , 1 ] x , y \in [ - 1,1 ] x,y∈[−1,1]
二、解决思路和方案
遗传算法是一种借鉴生物界自然选择和遗传机制的随机化搜索最优解的算法。它模拟自然选择和遗传过程中发生的环境选择、有性繁殖和基因突变现象,在每次迭代都从种群中选取较优的个体,利用交叉算子和变异算子使个体发生交叉(产生新个体的主要方式)和变异(产生新个体的辅助方式,有利于跳出函数的局部最优解),产生适应性更好(使函数更加最优解)的新一代种群,重复此过程,直到满足某种收敛指标为止。应用遗传算法需确定的几个主要内容是:初始种群的产生、适应度函数的确定、遗传算子(选择、交叉、变异)的确定、终止条件。
(1)初始种群的产生
产生初始种群的编码采取二进制编码方式,在产生初始种群之前,要确定染色体的基因数目,这是由算法想达到的精度决定的,如果要求最终求得的决策变量精确为5(精确到小数点后五位),根据公式 2 l j − 1 < ( b j − a j ) × 1 0 5 ≤ 2 l j − 1 2 ^ { l _ { j } - 1 } < \left( b _ { j } - a _ { j } \right) \times 10 ^ { 5 } \leq 2 ^ { l _ { j } } - 1 2lj−1<(bj−aj)×105≤2lj−1,即可确定该决策变量需要的基因数 l j l _ { j } lj,其中 b j b_ { j } bj 和 a j a_ { j } aj 为该决策变量取值的上下界,一个个体所需的基因数目为所有决策变量对应的 l j l _ { j } lj之和。
在基因数目确定了之后,即可产生初始种群,采取二进制编码随机生成的方式,
(2)适应度函数
适应度评价函数反映了个体对环境的适应能力,在以上几个最小化的函数中,f函数值越小说明其适应能力越强。为了符合常规习惯,将求 f 最小值的问题转化为求 -f 的最大值,那么适应度评价函数即为 -f, 其值越大,适应能力越强,越应该被保留。
(3)遗传算子的确定
对于选择步骤,采取轮盘赌选择,其思想是适应度更大的个体更大概率会被保留。另外加入了精英保留策略,即每次选择时都将适应度最高的个体保留下来不淘汰,这样可以优化种群质量,提升收敛速度。
对于交叉步骤,采用单点交叉方式,交叉算子设置为0.75,交叉是针对个体而言的,如果随机产生的概率小于0.75,则该个体可以参与交叉,随机与另外一个个体交叉,生成新的子代。
对于变异步骤,采用单点变异,变异算子设置为0.01,变异是针对染色体的每个基因而言的,如果随机产生的概率小于0.01,则该基因位发生变异,0和1互换。
(4)终止条件
首先设置了一个进化代数的上下限,最多迭代1000代,最小100代。另外,如果持续很多代最优值的改善幅度很小,便会提前终止,在算法执行时间和精确性之间取得了一个平衡。这里,如果持续200代最优值的改善幅度不超过1e-6(10^-6)便提前终止。
使用python语言进行仿真实现,求解每个最优化函数的最优解和对应的决策变量值,并可视化进化代数和最优值的关系。决策变量的精度为5,目标变量精度为6。
三、仿真结果及分析
(1)函数 f1 仿真结果
函数
f
1
(
x
,
y
)
=
x
2
+
y
2
f _ { 1 } ( x , y ) = x ^ { 2 } + y ^ { 2 }
f1(x,y)=x2+y2,
x
,
y
∈
[
−
5.12
,
5.12
]
x , y \in [ - 5.12,5.12 ]
x,y∈[−5.12,5.12]在迭代了307代后收敛,图3.1为迭代结束后的结果展示,在x=-0.00049,y=0.00285处取得了最优值8e-06,基本收敛到精确解0。其他函数同,其他函数便不进行解释,只进行结果展示。
图3.1 f1 求解结果
图3.2可视化了函数 f1 和x, y之间的关系,红点表示最后函数收敛到了最小。其他函数同。
图3.2 f1 可视化
图3.3表示了随着进化代数的增加,“当代最优值”(绿色曲线)和“目前为止最优值”(红色曲线)的变化情况,可以看出,红色曲线始终在绿色的下方,说明最优值不会变差,当代最优值虽有波动,但整体来看,也在变好。并且,图中每隔50代标注了目前为止的最优值(具体数值见表3.1),可以看出,随着代数的增加,最优值的改善不会很明显,同时函数也基本收敛。其他函数同。
图3.3 f1 进化代数-最优值变化图
表3.1 f1 进化代数-最优值变化表
进化代数 | 1 | 51 | 101 | 151 | 201 | 251 | 301 |
---|---|---|---|---|---|---|---|
最优值 | 0.247035 | 0.001725 | 1.2e-05 | 8e-06 | 8e-06 | 8e-06 | 8e-06 |
(2)函数 f2 仿真结果
函数
f
2
(
x
,
y
)
=
100
(
y
−
x
2
)
2
+
(
x
−
1
)
2
f _ { 2 } ( x , y ) = 100 \left( y - x ^ { 2 } \right) ^ { 2 } + ( x - 1 ) ^ { 2 }
f2(x,y)=100(y−x2)2+(x−1)2,
x
,
y
∈
[
−
2.048
,
2.048
]
x , y \in [ - 2.048,2.048 ]
x,y∈[−2.048,2.048]
图3.4 f2 求解结果
图3.5 f2 可视化
图3.6 f2 进化代数-最优值变化图
表3.2 f2 进化代数-最优值变化表
进化代数 | 1 | 51 | 101 | 151 | 201 |
---|---|---|---|---|---|
最优值 | 1.418274 | 0.00648 | 0.00648 | 0.00518 | 0.002211 |
进化代数 | 251 | 301 | 351 | 401 | 451 |
最优值 | 0.000868 | 0.00085 | 0.00085 | 0.00085 | 0.00085 |
(3)函数 f3 仿真结果
函数
f
3
(
x
,
y
)
=
(
x
2
+
y
−
11
)
2
+
(
x
+
y
2
−
7
)
2
f _ { 3 } ( x , y ) = \left( x ^ { 2 } + y - 11 \right) ^ { 2 } + \left( x + y ^ { 2 } - 7 \right) ^ { 2 }
f3(x,y)=(x2+y−11)2+(x+y2−7)2,
x
,
y
∈
[
−
6
,
6
]
x , y \in [ - 6,6 ]
x,y∈[−6,6]
图3.7 f3 求解结果
图3.8 f3 可视化
图3.9 f3 进化代数-最优值变化图
表3.3 f3 进化代数-最优值变化表
进化代数 | 1 | 51 | 101 | 151 | 201 | 251 | 301 |
---|---|---|---|---|---|---|---|
最优值 | 37.36007 | 0.177718 | 0.177718 | 0.06862 | 0.021468 | 0.001839 | 0.001839 |
(4)函数 f4 仿真结果
函数
f
4
(
x
,
y
)
=
4
x
2
−
2.1
x
4
+
x
6
/
3
+
x
y
+
4
y
4
−
4
y
2
f _ { 4 } ( x , y ) =4 x ^ 2 - 2.1 x ^ 4 + x ^6 / 3 + xy + 4 y ^4 - 4 y^ 2
f4(x,y)=4x2−2.1x4+x6/3+xy+4y4−4y2,
x
,
y
∈
[
−
5
,
5
]
x , y \in [ - 5,5 ]
x,y∈[−5,5]
图3.10 f4 求解结果
图3.11 f4 可视化
图3.12 f4 进化代数-最优值变化图
表3.4 f4 进化代数-最优值变化表
进化代数 | 1 | 51 | 101 | 151 | 201 | 251 |
---|---|---|---|---|---|---|
最优值 | -0.836399 | -1.029728 | -1.030549 | -1.030549 | -1.030549 | -1.031549 |
进化代数 | 301 | 351 | 401 | 451 | 501 | 551 |
最优值 | -1.031576 | -1.031617 | -1.031625 | -1.031625 | -1.031625 | -1.031625 |
(5)函数 f5 仿真结果
函数
f
5
(
x
,
y
)
=
x
2
+
2
y
2
−
0.3
cos
3
π
x
cos
4
π
y
+
0.3
f _ { 5 } ( x , y ) = x ^ { 2 } + 2 y ^ { 2 } - 0.3 \cos 3 \pi x \cos 4 \pi y + 0.3
f5(x,y)=x2+2y2−0.3cos3πxcos4πy+0.3,
x
,
y
∈
[
−
1
,
1
]
x , y \in [ - 1,1 ]
x,y∈[−1,1]
图3.13 f5 求解结果
图3.14 f5 可视化
图3.15 f5 进化代数-最优值变化图
表3.5 f5 进化代数-最优值变化表
进化代数 | 1 | 51 | 101 | 151 | 201 | 251 |
---|---|---|---|---|---|---|
最优值 | 0.261395 | 0.023016 | 0.005106 | 0.000206 | 0.0 | 0.0 |
代码:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 函数f1 x和y的取值范围
x_bound = [-5.12, 5.12]
y_bound = [-5.12, 5.12]
# 函数f2 x和y的取值范围
# x_bound = [-2.048, 2.048]
# y_bound = [-2.048, 2.048]
# 函数f3 x和y的取值范围
# x_bound = [-6, 6]
# y_bound = [-6, 6]
# 函数f4 x和y的取值范围
# x_bound = [-5, 5]
# y_bound = [-5, 5]
# 函数f5 x和y的取值范围
# x_bound = [-1, 1]
# y_bound = [-1, 1]
X_precision = 5 # 决策变量的精度
y_precision = 6 # 目标变量的精度
pop_size = 30 # 种群个体数
crossover_rate = 0.75 # 交叉算子
mutation_rate = 0.01 # 变异算子
max_generation = 1000 # 最多迭代的进化次数
min_generation = 100 # 最少需迭代的进化次数
without_optim_tolerate = 200 # 如果持续200代最优值改善很小(1e-6)的话,提前终止迭代
# 利用决策变量的精度和取值区间计算基因数
def cal_DNA_size(x_bound, y_bound, X_precision):
x_size, y_size = 1, 1
discriminant_X = (x_bound[1] - x_bound[0]) * 10 ** X_precision
discriminant_Y = (y_bound[1] - y_bound[0]) * 10 ** X_precision
while 2 ** (x_size - 1) < discriminant_X:
if discriminant_X < 2 ** x_size - 1:
break
x_size = x_size + 1
while 2 ** (y_size - 1) < discriminant_Y:
if discriminant_Y < 2 ** y_size - 1:
break
y_size = y_size + 1
return x_size, y_size
x_size, y_size = cal_DNA_size(x_bound, y_bound, X_precision)
DNA_size = x_size + y_size
# 适应度计算,即函数值,函数表达式乘-1
def get_fitness(x, y):
return -1 * (x ** 2 + y ** 2)
# return -1 * (100 * (y - x ** 2) ** 2 + (x - 1) ** 2)
# return -1 * ((x ** 2 + y - 11) ** 2 + (x + y ** 2 - 7) ** 2)
# return -1 * (4 * x ** 2 - 2.1 * x ** 4 + x ** 6 / 3 + x * y + 4 * y ** 4 - 4 * y ** 2)
# return -1 * (x ** 2 + 2 * y ** 2 - 0.3 * np.cos(3 * np.pi * x) * np.cos(4 * np.pi * y) + 0.3)
# 二进制转十进制,变换到决策变量取值区间
def binary_to_decimal(pop, bound, size):
return np.around(bound[0] + pop.dot(2 ** np.arange(size)[::-1]) * (bound[1] - bound[0]) / (2 ** x_size-1), decimals=X_precision)
# 模拟自然选择,适应度越大,越大概率被保留
def select(pop, fitness):
idx = np.random.choice(np.arange(pop_size), size = pop_size - 1, replace = True,
p = (fitness + (abs(min(fitness)) if min(fitness) < 0 else 0)) / (fitness.sum() + abs(min(fitness) if min(fitness) < 0 else 0) * pop_size))
return pop[idx]
# 模拟交叉,生成新的子代
def crossover(parent, pop):
if np.random.rand() < crossover_rate:
i_ = np.random.randint(0, pop_size, size=1)
cross_points = np.random.randint(0, DNA_size, size=1)
cross_points = cross_points[0]
return np.append(parent[0:cross_points], pop[i_, cross_points:])
return parent
# 模拟变异
def mutate(child):
for point in range(DNA_size):
if np.random.rand() < mutation_rate:
child[point] = 1 if child[point] == 0 else 0
return child
pop = np.random.randint(2, size=(pop_size, DNA_size))
plt.ion() # 打开plt交互模式
fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(x_bound[0], x_bound[1], (x_bound[1] - x_bound[0])/50)
Y = np.arange(y_bound[0], y_bound[1], (y_bound[1] - y_bound[0])/50)
X, Y = np.meshgrid(X, Y)
Z = -1 * get_fitness(X, Y)
# 绘制3D曲面
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='coolwarm')
best_F = float('-inf')
worse_F = float('inf')
best_x, best_y = 0, 0
actual_generation = 0
cur_best_F_list = []
best_F_list = []
# 核心代码,进行迭代进化
for i in range(max_generation):
x = binary_to_decimal(pop[:, 0:x_size], x_bound, x_size)
y = binary_to_decimal(pop[:, x_size:], y_bound, y_size)
fitness = np.around(get_fitness(x, y), decimals = y_precision)
cur_best_x = binary_to_decimal(pop[np.argmax(fitness), 0:x_size], x_bound, x_size)
cur_best_y = binary_to_decimal(pop[np.argmax(fitness), x_size:], y_bound, y_size)
print("generation:", i + 1)
print("most fitted DNA: ", pop[np.argmax(fitness), :])
print("var corresponding to most fitted DNA: ", cur_best_x, cur_best_y)
print("F_values corresponding to DNA: ", -1 * fitness[np.argmax(fitness)])
if fitness[np.argmax(fitness)] > best_F:
best_F = fitness[np.argmax(fitness)]
best_x = cur_best_x
best_y = cur_best_y
if fitness[np.argmax(fitness)] < worse_F:
worse_F = fitness[np.argmax(fitness)]
# 判断是否需要提前终止迭代
if i + 1 > min_generation and (best_F - best_F_list[i - (without_optim_tolerate if i > without_optim_tolerate else i)]) < 10 ** (-y_precision):
actual_generation = i + 1
break
# 逐代绘制,绘制前先清除上一代的
if 'sca' in globals():
sca.remove()
sca = ax.scatter(best_x, best_y, fitness[np.argmax(fitness)], s=200, lw=0, c='red', alpha=0.5)
plt.pause(0.001)
# 精英保留
pop = np.vstack((select(pop, fitness), pop[np.argmax(fitness), :]))
pop_copy = pop.copy() # parent会被child替换,所以先copy一份pop
for parent in pop:
child = crossover(parent, pop_copy)
child = mutate(child)
parent[:] = child
cur_best_F_list.append(fitness[np.argmax(fitness)])
best_F_list.append(best_F)
# 迭代完输出最优值和对应的决策变量
print("best F_value is", -1 * best_F)
print("var corresponding to best F_value is", best_x, best_y)
plt.ioff()
plt.show()
# 绘制进化代数和最优解的关系图
best_F_list = -1 * np.array(best_F_list)
cur_best_F_list = -1 * np.array(cur_best_F_list)
plt.plot(range(actual_generation - 1), best_F_list, label='best solution so far', color='red')
plt.plot(range(actual_generation - 1), cur_best_F_list, label='best solution of current generation', color='green')
# 每隔50代标注一下最优值
l = [i for i in range(actual_generation) if i % 50 == 0]
for x, y in zip(l, best_F_list[l]):
print(x + 1, y)
plt.text(x, y + 0.001, f'%.{y_precision}f' % y, ha='center', va= 'bottom',fontsize=12)
# 参数标注
plt.text(actual_generation * 0.5, (best_F + worse_F) / -2, "crossover_rate=%.2f, mutation_rate=%.2f, actual_generation=%d" % (crossover_rate, mutation_rate, actual_generation), fontdict={'size': 15, 'color': 'black'})
plt.legend(loc='best', fontsize=20)
plt.title('generation vs. F_value of f1', fontsize=25, color='black')
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
fig = plt.gcf()
fig.set_size_inches(18, 10)
plt.savefig('f1_generation_F_value.png')
plt.show()