目录
遗传算法简介
遗传算法是一种启发式算法,通过模拟自然界中的生物进化过程来解决优化问题。遗传算法的基本思想源于达尔文的进化论,通过模拟自然选择、交叉、变异等操作,逐步优化解空间中的个体,寻找问题的最优解或者较好的解。具体来说,遗传算法包括以下几个步骤:
- 初始化种群:随机生成一定数量的个体作为初始种群。
- 选择:根据适应度函数选出优秀个体,通常选择适应度较高的个体用于后续繁殖。
- 交叉:通过某种交叉方式,将选出的优秀个体进行交叉繁殖,产生下一代个体。
- 变异:对新生代的个体进行一定概率的变异操作,引入新的基因信息,增加种群的多样性。
- 竞争:将经过交叉和变异产生的个体与父代进行竞争,从而保留适应度较高的个体作为下一代。
遗传算法在解决复杂优化问题、搜索空间较大或连续的问题时具有一定优势,并在实际应用中得到了广泛的应用
优化问题
数学优化问题如下,这是一个二维非线性优化问题,目标函数和约束条件中都带有非线性约束。
编程建模
初始化种群
常用的编码方式有二进制编码和实数编码,本文采用实数编码方式进行编码。由于这个问题是一个二维非线性优化问题,因此每个个体包含两个基因片段对应和。初始化种群的主要难点在于如何考虑优化问题的约束条件,使得生成的初始种群尽可能可行,减小可行解搜索的空间。
我们观察到约束1,它的可行域区间实际上是以(0, 0)为圆心,以为半径的圆,而约束2又要求,因此它的可行域区间实际上是第一象限的四分之一圆,并且和还需要满足不等式约束。故可行域为四分之一圆和不等式约束划分的区域内。故生成种群时可以让,这样生成的初始种群虽然不是都可行,但是缩小了可行解空间。代码如下
def initializing_populations(pop_number):
"""
生成初始化种群
:param pop_number: 种群数
"""
temp_pop_array1 = np.random.uniform(0, math.sqrt(5)+1e-7, pop_number).reshape(-1, 1)
temp_pop_array2 = np.random.uniform(0, 2, pop_number).reshape(-1, 1)
pop_array = np.hstack((temp_pop_array1, temp_pop_array2))
return pop_array
适应度计算函数
这个问题是个最大化问题,我们选用各个体目标函数的计算值作为个体的适应度。
def fitness_calculation_function(pop_array):
"""
适应度计算函数,即目标函数,求当前约束下目标函数的最大值
:param pop_array:
:return: 各个体的目标函数值,每一行表示一个个体目标函数值
"""
fit_array = (pop_array[:, 0]-3)**2+(pop_array[:, 1]-2)**2
return fit_array.reshape(-1, 1)
边界约束破坏处理
该优化问题的边界约束都是不等式约束,因此可以采用惩罚函数的方法,对不满足约束条件的个体的适应度进行惩罚,减小其产生下一代的概率。
def penalty_function(pop_array, fit_array):
"""
对于不满足不等式约束的个体的进行惩罚
:param fit_array: 适应度函数
:param pop_array: 种群
:return: 加入惩罚后的适应度
"""
penalty_array = 4-(pop_array[:, 0]+2*pop_array[:, 1])
statis_index = np.where(penalty_array >= 0) # 满足不等式约束的索引
penalty_array[statis_index] = 0 # 将满足约束的惩罚值更改为0
new_fit_array = fit_array+penalty_array.reshape(-1, 1) # 将惩罚值加入适应度数组,得到新的适应度数组
zero_index = np.where(new_fit_array < 0) # 惩罚后的适应度如果小于0,将其更改为0
new_fit_array[zero_index] = 0
return new_fit_array
选择高适应度的个体
遗传算法中常用的选择方法有:
-
轮盘赌选择:根据每个个体的适应度值比例来选择个体,适应度值高的个体有更大的概率被选中。
-
锦标赛选择:随机选取一定数量的个体,根据它们的适应度值选择较优秀的个体作为父代。
-
最优选择:直接选择适应度值最高的个体作为父代。
-
排序选择:将个体按适应度值从高到低进行排序,根据排名选取个体作为父代。
-
随机选择:随机选择个体作为父代,不考虑适应度值
本文选择轮盘赌法,轮盘赌适用于处理适应度为非负数的情况 。根据适应度为每个个体附加相应的杂交概率,适应度越大的个体杂交产生下一代的可能性越大。
def roulette_selection(fit_array):
"""
轮盘赌选择法,根据适应度的大小,对每个个体附加交叉的概率,适应度越大交叉的可能性越大
:param fit_array: 适应度数组
:return: 返回种群内各个个体的交叉概率
"""
crossover_p_array = np.full((fit_array.size, 1), np.nan, dtype=float)
fit_sum = np.sum(fit_array)
if abs(fit_sum) <= 1e-4:
sys.exit("适应度之和不能为0,请更换选择方法")
accumulate_p = 0 # 累积概率
for i in range(fit_array.size):
crossover_p_array[i, 0] = fit_array[i, 0]/fit_sum+accumulate_p
accumulate_p = crossover_p_array[i, 0]
return crossover_p_array
杂交运算
根据输入的杂交概率和种群数量,随机抽取个体进行杂交。
def crossover_operation(pop_number, crossover_p, parent_pop, crossover_p_array):
"""
交叉运算,高适应度的个体交叉的概率大,低适应度的个体交叉的概率小
:param crossover_p_array: 种群的交叉概率数组
:param pop_number: 种群数量
:param parent_pop: 父代
:param crossover_p: 交叉概率
:return: 返回交叉后的种群
"""
new_parent_pop = copy(parent_pop)
for i in range(int(crossover_p*pop_number)): # 总交叉次数
temp_crossover_p1 = random.uniform(0, 1) # 随机选择父代概率1
temp_crossover_p2 = random.uniform(0, 1) # 随机选择父代概率2
parent_index1 = np.where(temp_crossover_p1 <= crossover_p_array)[0][0] # 父代1行索引
parent_index2 = np.where(temp_crossover_p2 <= crossover_p_array)[0][0] # 父代2行索引
random_int = np.random.randint(2) # 生成0-1之间的随机整数确定交叉的位置
# 根据父代进行交叉得到子代
new_parent_pop[parent_index1, random_int] = parent_pop[parent_index2, random_int]
new_parent_pop[parent_index2, random_int] = parent_pop[parent_index2, random_int]
return new_parent_pop
突变运算
种群中所有的个体突变的概率都相同。
def mutation_operation(pop_number, mutation_p, parent_pop):
"""
突变函数,对种群进行突变运算
:param pop_number:
:param mutation_p: 突变概率
:param parent_pop:
:return:
"""
new_parent_pop = parent_pop
for i in range(int(pop_number*mutation_p)):
random_int1 = np.random.randint(pop_number) # 突变个体
random_int2 = np.random.randint(2) # 突变位置
if random_int2 == 0:
new_parent_pop[random_int1, random_int2] = random.uniform(0, math.sqrt(5) + 1e-7)
else:
new_parent_pop[random_int1, random_int2] = random.uniform(0, 2)
return new_parent_pop
主循环
根据输入的最大繁殖代数进行循环迭代。
def main_loop(pop_number, reproduction_number, crossover_p, mutation_p):
"""
遗传算法实现的主循环
:param pop_number: 种群数量
:param reproduction_number: 最大的繁殖代数
:param crossover_p: 杂交概率
:param mutation_p: 突变概率
"""
init_pop_array = initializing_populations(pop_number) # 初始化种群
for i in range(reproduction_number):
if i == 0:
parent_pop = init_pop_array # 父代种群
else:
parent_pop = child_pop_array
fit_array = penalty_function(parent_pop, fitness_calculation_function(parent_pop)) # 获取父本中个体的适应度
crossover_p_array = roulette_selection(fit_array) # 获取各个体的杂交概率
temp_pop_array = crossover_operation(pop_number, crossover_p, parent_pop, crossover_p_array) # 进行交叉运算,生成中间结果的种群
temp_pop_array = mutation_operation(pop_number, mutation_p, temp_pop_array) # 进行突变运算
temp_fit_array = penalty_function(temp_pop_array, fitness_calculation_function(temp_pop_array)) # 临时种群适应度
child_pop_array = competitiveness_operation(pop_number, temp_pop_array,temp_fit_array, parent_pop, fit_array)
结果分析
约束条件为不等式
输入参数:种群数20,最大繁殖代数1000,杂交概率0.9,突变概率0.3
if __name__ == "__main__":
main_loop(20, 1000, 0.9, 0.3)
结果表示适应度值,的值。
12.990864951472929 [0.00074035 0.00117371] # 第一次运行结果
12.983539152035894 [0.00057131 0.00326098] # 第二次运行结果
12.986867969022803 [0.00157115 0.00092712] # 第三次运行结果
检验不等式结果是否正确
首先上述结果都满足约束条件,我们采用开源非线性求解器ipopt来检验结果是否最优。Ipopt(Interior Point OPTimizer)是一种开源软件包,用于求解非线性连续优化问题。它采用内点方法来解决大规模非线性优化问题,可以处理包括等式约束、不等式约束和边界约束等多种类型的约束条件。Ipopt能够高效地求解复杂的非线性优化问题,通常用于在数学建模、工程优化、经济学等领域中的应用。Ipopt在连续优化领域有着广泛的应用,且性能表现良好,是许多优化问题的首选求解器之一。(本文重点不在此,用法不做介绍)
from pyomo.environ import *
model = ConcreteModel()
model.x1 = Var(within=NonNegativeReals)
model.x2 = Var(within=NonNegativeReals)
model.obj_value = Var()
def rule1(m):
return m.x1**2+m.x2**2-5 <= 0
model.constraint1 = Constraint(rule=rule1)
def rule2(m):
return m.x1+2*m.x2-4 <= 0
model.constraint2 = Constraint(rule=rule2)
def get_obj_value(m):
return model.obj_value == (m.x1 - 3) ** 2 + (m.x2 - 2) ** 2
model.get_obj_value_con = Constraint(rule=get_obj_value)
def obj_rule(m):
return model.obj_value
model.objective_constraint = Objective(rule=obj_rule, sense=maximize)
opt = SolverFactory('ipopt')
opt.solve(model, tee=True)
print(value(model.obj_value), value(model.x1), value(model.x2))
所得结果如下:
13.0 0.0 0.0
和遗传算法求解的一致,说明执行的算法是有效的。
约束条件为等式
当我们将约束条件修改为等式之后,采用惩罚方法还能得到稳定的结果吗?
def penalty_function2(pop_array, fit_array):
"""
将约束更改为等式约束
:param fit_array: 适应度函数
:param pop_array: 种群
:return: 加入惩罚后的适应度
"""
penalty_array = 4-(pop_array[:, 0]+2*pop_array[:, 1])
index_lower = np.where(penalty_array >= 0)
index_upper = np.where(penalty_array <= 1e-6)
# 取交集,即既满足大于等于0,又满足小于等于1e-6的索引
statis_index = np.intersect1d(index_lower, index_upper)
penalty_array[statis_index] = 0 # 将满足约束的惩罚值更改为0
new_fit_array = fit_array-np.abs(penalty_array.reshape(-1, 1)) # 将惩罚值加入适应度数组,得到新的适应度数组
zero_index = np.where(new_fit_array < 0)
new_fit_array[zero_index] = 0
return new_fit_array
其余参数和条件不变,运行几次输出结果
8.987463656361355 [1.87469013e-03 1.99841554e+00] # 第一次运行结果
8.998749946120093 [2.44311335e-04 1.42785488e-05] # 第二次运行结果
8.993601963219662 [0.00088544 0.00098629] # 第三次运行结果
从上面可以看出输入条件相同,但是所得结果却相差较大,并且不满足等式约束,无法像不等式约束一样获得较稳定的最优解,那么如何解决这个问题呢?利用等式约束增加一个修正,修正个体的基因编码。
def genetic_modification(pop_array):
"""
对个体的基因型进行修正,根据等式约束
:param pop_array:
:return:
"""
new_pop_array = copy(pop_array)
new_pop_array[:, 1] = 2-0.5*new_pop_array[:, 0]
return new_pop_array
运行三次结果如下:
8.990181777990122 [1.63692857e-03 1.99918154e+00] # 第一次运行结果
8.991384439727147 [1.43635653e-03 1.99928182e+00] # 第二次运行结果
8.998369968364907 [2.71687317e-04 1.99986416e+00] # 第三次运行结果
检验等式结果是否正确
model = ConcreteModel()
model.x1 = Var(within=NonNegativeReals)
model.x2 = Var(within=NonNegativeReals)
model.obj_value = Var()
def rule1(m):
return m.x1 ** 2 + m.x2 ** 2 - 5 <= 0
model.constraint1 = Constraint(rule=rule1)
def rule2(m):
return m.x1 + 2 * m.x2 - 4 == 0
model.constraint2 = Constraint(rule=rule2)
def get_obj_value(m):
return model.obj_value == (m.x1 - 3) ** 2 + (m.x2 - 2) ** 2
model.get_obj_value_con = Constraint(rule=get_obj_value)
def obj_rule(m):
return model.obj_value
model.objective_constraint = Objective(rule=obj_rule, sense=maximize)
opt = SolverFactory('ipopt')
opt.solve(model, tee=True)
print(value(model.obj_value), value(model.x1), value(model.x2))
结果如下:
9.0 0.0 2.0
完整代码:
import sys
import random
from copy import copy
import numpy as np
import math
# 初始化种群
def initializing_populations(pop_number):
"""
生成初始化种群
:param gen_number: 每个个体的基因数量
:param pop_number: 种群数
"""
temp_pop_array1 = np.random.uniform(0, math.sqrt(5)+1e-7, pop_number).reshape(-1, 1)
temp_pop_array2 = np.random.uniform(0, 2, pop_number).reshape(-1, 1)
pop_array = np.hstack((temp_pop_array1, temp_pop_array2))
return pop_array
def fitness_calculation_function(pop_array):
"""
适应度计算函数,即目标函数,求当前约束下目标函数的最大值
:param pop_array:
:return: 各个体的目标函数值,每一行表示一个个体目标函数值
"""
fit_array = (pop_array[:, 0]-3)**2+(pop_array[:, 1]-2)**2
return fit_array.reshape(-1, 1)
def roulette_selection(fit_array):
"""
轮盘赌选择法,根据适应度的大小,对每个个体附加交叉的概率,适应度越大交叉的可能性越大
:param fit_array: 适应度数组
:return: 返回种群内各个个体的交叉概率
"""
crossover_p_array = np.full((fit_array.size, 1), np.nan, dtype=float)
fit_sum = np.sum(fit_array)
if abs(fit_sum) <= 1e-4:
sys.exit("适应度之和不能为0,请更换选择方法")
accumulate_p = 0 # 累积概率
for i in range(fit_array.size):
crossover_p_array[i, 0] = fit_array[i, 0]/fit_sum+accumulate_p
accumulate_p = crossover_p_array[i, 0]
return crossover_p_array
def penalty_function(pop_array, fit_array):
"""
对于不满足等式约束的个体的进行惩罚
:param fit_array: 适应度函数
:param pop_array: 种群
:return: 加入惩罚后的适应度
"""
penalty_array = 4-(pop_array[:, 0]+2*pop_array[:, 1])
statis_index = np.where(penalty_array >= 0) # 满足不等式约束的索引
penalty_array[statis_index] = 0 # 将满足约束的惩罚值更改为0
new_fit_array = fit_array+penalty_array.reshape(-1, 1) # 将惩罚值加入适应度数组,得到新的适应度数组
zero_index = np.where(new_fit_array < 0)
new_fit_array[zero_index] = 0
return new_fit_array
def crossover_operation(pop_number, crossover_p, parent_pop, crossover_p_array):
"""
交叉运算,高适应度的个体交叉的概率大,低适应度的个体交叉的概率小
:param crossover_p_array: 种群的交叉概率数组
:param pop_number: 种群数量
:param parent_pop: 父代
:param crossover_p: 交叉概率
:return: 返回交叉后的种群
"""
new_parent_pop = copy(parent_pop)
for i in range(int(crossover_p*pop_number)): # 总交叉次数
temp_crossover_p1 = random.uniform(0, 1) # 随机选择父代概率1
temp_crossover_p2 = random.uniform(0, 1) # 随机选择父代概率2
parent_index1 = np.where(temp_crossover_p1 <= crossover_p_array)[0][0] # 父代1行索引
parent_index2 = np.where(temp_crossover_p2 <= crossover_p_array)[0][0] # 父代2行索引
random_int = np.random.randint(2) # 生成0-1之间的随机整数确定交叉的位置
# 根据父代进行交叉得到子代
new_parent_pop[parent_index1, random_int] = parent_pop[parent_index2, random_int]
new_parent_pop[parent_index2, random_int] = parent_pop[parent_index2, random_int]
return new_parent_pop
def mutation_operation(pop_number, mutation_p, parent_pop):
"""
突变函数,对种群进行突变运算
:param pop_number:
:param mutation_p: 突变概率
:param parent_pop:
:return:
"""
new_parent_pop = parent_pop
for i in range(int(pop_number*mutation_p)):
random_int1 = np.random.randint(pop_number) # 突变个体
random_int2 = np.random.randint(2) # 突变位置
if random_int2 == 0:
new_parent_pop[random_int1, random_int2] = random.uniform(0, math.sqrt(5) + 1e-7)
else:
new_parent_pop[random_int1, random_int2] = random.uniform(0, 2)
return new_parent_pop
def competitiveness_operation(pop_number, temp_pop_array, temp_fit_array, parent_pop, parent_fit_array):
"""
令生成的temp_pop_array与父代进行竞争, 降序排列,选取前pop_number作为子代
:param parent_fit_array: 父代适应度
:param pop_number: 种群数
:param temp_fit_array: 临时种群适应度
:param temp_pop_array: 临时种群
:param parent_pop: 父代
:return: 子代
"""
fit_array = np.vstack((parent_fit_array, temp_fit_array))
sort_fit_array = sorted(list(fit_array[:, 0]), reverse=True)
child_generation = copy(parent_pop)
for i in range(pop_number):
child_index = np.where(sort_fit_array[i] == fit_array)[0][0]
if child_index > pop_number-1:
child_generation[i, :] = temp_pop_array[child_index-pop_number, :]
else:
child_generation[i, :] = parent_pop[child_index, :]
print(sort_fit_array[0], child_generation[0, :])
return child_generation
def main_loop(pop_number, reproduction_number, crossover_p, mutation_p):
"""
遗传算法实现的主循环,当约束为不等式
:param pop_number: 种群数量
:param reproduction_number: 最大的繁殖代数
:param crossover_p: 杂交概率
:param mutation_p: 突变概率
"""
init_pop_array = initializing_populations(pop_number) # 初始化种群
for i in range(reproduction_number):
if i == 0:
parent_pop = init_pop_array # 父代种群
else:
parent_pop = child_pop_array
fit_array = penalty_function(parent_pop, fitness_calculation_function(parent_pop)) # 获取父本中个体的适应度
crossover_p_array = roulette_selection(fit_array) # 获取各个体的杂交概率
temp_pop_array = crossover_operation(pop_number, crossover_p, parent_pop, crossover_p_array) # 进行交叉运算,生成中间结果的种群
temp_pop_array = mutation_operation(pop_number, mutation_p, temp_pop_array) # 进行突变运算
temp_fit_array = penalty_function(temp_pop_array, fitness_calculation_function(temp_pop_array)) # 临时种群适应度
child_pop_array = competitiveness_operation(pop_number, temp_pop_array,temp_fit_array, parent_pop, fit_array)
def penalty_function2(pop_array, fit_array):
"""
将约束更改为等式约束
:param fit_array: 适应度函数
:param pop_array: 种群
:return: 加入惩罚后的适应度
"""
penalty_array = 4-(pop_array[:, 0]+2*pop_array[:, 1])
index_lower = np.where(penalty_array >= 0)
index_upper = np.where(penalty_array <= 1e-6)
# 取交集,即既满足大于等于0,又满足小于等于1e-6的索引
statis_index = np.intersect1d(index_lower, index_upper)
penalty_array[statis_index] = 0 # 将满足约束的惩罚值更改为0
new_fit_array = fit_array-np.abs(penalty_array.reshape(-1, 1)) # 将惩罚值加入适应度数组,得到新的适应度数组
zero_index = np.where(new_fit_array < 0)
new_fit_array[zero_index] = 0
return new_fit_array
def genetic_modification(pop_array):
"""
对个体的基因型进行修正,根据等式约束
:param pop_array:
:return:
"""
new_pop_array = copy(pop_array)
new_pop_array[:, 1] = 2-0.5*new_pop_array[:, 0]
return new_pop_array
def main_loop2(pop_number, reproduction_number, crossover_p, mutation_p):
"""
遗传算法实现的主循环,当约束为等式
:param pop_number: 种群数量
:param reproduction_number: 最大的繁殖代数
:param crossover_p: 杂交概率
:param mutation_p: 突变概率
"""
init_pop_array = initializing_populations(pop_number) # 初始化种群
for i in range(reproduction_number):
if i == 0:
parent_pop = init_pop_array # 父代种群
else:
parent_pop = child_pop_array
parent_pop = genetic_modification(parent_pop)
fit_array = fitness_calculation_function(parent_pop) # 获取父本中个体的适应度
crossover_p_array = roulette_selection(fit_array) # 获取各个体的杂交概率
temp_pop_array = crossover_operation(pop_number, crossover_p, parent_pop, crossover_p_array) # 进行交叉运算,生成中间结果的种群
temp_pop_array = mutation_operation(pop_number, mutation_p, temp_pop_array) # 进行突变运算
temp_pop_array = genetic_modification(temp_pop_array)
temp_fit_array = fitness_calculation_function(temp_pop_array) # 临时种群适应度
child_pop_array = competitiveness_operation(pop_number, temp_pop_array, temp_fit_array, parent_pop, fit_array)
from pyomo.environ import *
# 检验不等式结果是否成立
def check_result1():
model = ConcreteModel()
model.x1 = Var(within=NonNegativeReals)
model.x2 = Var(within=NonNegativeReals)
model.obj_value = Var()
def rule1(m):
return m.x1**2+m.x2**2-5 <= 0
model.constraint1 = Constraint(rule=rule1)
def rule2(m):
return m.x1+2*m.x2-4 <= 0
model.constraint2 = Constraint(rule=rule2)
def get_obj_value(m):
return model.obj_value == (m.x1 - 3) ** 2 + (m.x2 - 2) ** 2
model.get_obj_value_con = Constraint(rule=get_obj_value)
def obj_rule(m):
return model.obj_value
model.objective_constraint = Objective(rule=obj_rule, sense=maximize)
opt = SolverFactory('ipopt')
opt.solve(model, tee=True)
print(value(model.obj_value), value(model.x1), value(model.x2))
# 检验等式结果是否成立
def check_result2():
model = ConcreteModel()
model.x1 = Var(within=NonNegativeReals)
model.x2 = Var(within=NonNegativeReals)
model.obj_value = Var()
def rule1(m):
return m.x1 ** 2 + m.x2 ** 2 - 5 <= 0
model.constraint1 = Constraint(rule=rule1)
def rule2(m):
return m.x1 + 2 * m.x2 - 4 == 0
model.constraint2 = Constraint(rule=rule2)
def get_obj_value(m):
return model.obj_value == (m.x1 - 3) ** 2 + (m.x2 - 2) ** 2
model.get_obj_value_con = Constraint(rule=get_obj_value)
def obj_rule(m):
return model.obj_value
model.objective_constraint = Objective(rule=obj_rule, sense=maximize)
opt = SolverFactory('ipopt')
opt.solve(model, tee=True)
print(value(model.obj_value), value(model.x1), value(model.x2))
if __name__ == "__main__":
main_loop2(20, 1000, 0.9, 0.3)
check_result2()