遗传算法
原理
我们知道在自然进化的过程中个体会发生进化，朝着更加适应所处环境的方向发展。简单来说，遗传算法做的就是要模拟这一过程，经过一代又一代的进化使得可行解朝着更优的方向变化。遗传算法是一种全局最优搜索算法,相比于局部优化算法，如梯度下降等，在寻找全局最优解时更加具有鲁棒性。因为局部优化算法容易被初始点的位置所影响，可能会停留在某个局部最优解，而无法找到全局最优解。
相关生物遗传学概念在遗传算法中的含义
- 基因（gene）：染色体中的元素
- 染色体（chromosome）：对应于问题的可行解，是个体的表现形式（i.e.在遗传算法中用染色体表征一个个体）。算法中，染色体编码有很多方式，本文采用便于实现的十进制编码
- 种群（population）：进化发生的单位，由个体组成
- 适应值（fitness）：衡量个体对于环境的适应程度
- 选择（selection）：适应值高的个体存活下来
- 变异（variation）：染色体变异，基因突变
- 交叉互换（crossing）：产生后代时染色体发生交叉互换
注：交叉和变异在遗传算法中是很重要且必要的，可以帮助算法跳出搜索空间增加搜索空间的多样性，避免算法陷入局部最优解
流程
[1] 初始化种群
[2] 交叉，变异产生后代
[3] 自然选择，进化
[4] 重复[2][3]直到算法收敛或到达指定进化次数
在数学建模中，遗传算法主要用于解决函数优化和组合优化问题。下面用遗传算法求解著名的旅行商问题。
遗传算法求解旅行商问题（TSP）
导入相关python库
import numpy as np
import math
from numpy.random import randint,rand
import matplotlib.pyplot as plt # 用于可视化飞行路线
相关参数设置
population_size = 200#种群大小
num_generations = 50 # 迭代次数
mutation_rate = 0.05 # 变异率
crossing_rate = 0.8 # 染色体交叉互换率
length_of_chromosome=102 # 染色体长度
数据读取
读取上表17.1的数据
# 读取城市经纬坐标
file=r"D:\Study\my_notebook\数学建模相关\python数学实验与建模数据\程序及数据\程序及数据\17第17章 智能算法\Pdata17_1.txt"
city_locations=np.loadtxt(fname=file).reshape(-1,2)
start_point=np.array([70,40])
city_locations=np.vstack([start_point,city_locations,start_point])
计算得到各目标间距离的实对称矩阵
这里涉及到已知经纬度求各点间距的问题，使用Haversine公式计算
d
=
2
r
arcsin
sin
2
(
ϕ
2
−
ϕ
1
2
)
+
cos
(
ϕ
1
)
cos
(
ϕ
2
)
sin
2
(
λ
2
−
λ
1
2
)
d = 2r\arcsin\sqrt{\sin^2\left(\frac{\phi_2 - \phi_1}{2}\right) + \cos(\phi_1)\cos(\phi_2)\sin^2\left(\frac{\lambda_2 - \lambda_1}{2}\right)}
d=2rarcsinsin2(2ϕ2−ϕ1)+cos(ϕ1)cos(ϕ2)sin2(2λ2−λ1)
其中，
d
d
d 表示两点之间的距离，
r
r
r 表示地球半径，
ϕ
1
,
ϕ
2
\phi_1,\phi_2
ϕ1,ϕ2 和
λ
1
,
λ
2
\lambda_1,\lambda_2
λ1,λ2 分别表示两点的纬度和经度。
# 计算巡航点之间的距离
def calculate_distance(city_locations):
r= 6371*1e3 # 地球半径m
distances=np.zeros((len(city_locations),len(city_locations))) # 101x101的实对称数组,存储i到j的距离
longitude,latitude=city_locations[:,0],city_locations[:,1]
for i in range(len(city_locations)):
for j in range(len(city_locations)):
phi1 = math.radians(latitude[i])
phi2 = math.radians(latitude[j])
delta_phi = math.radians(latitude[j] - latitude[i])
delta_lambda = math.radians(longitude[j] - longitude[i])
a = math.sin(delta_phi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2) ** 2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
distance=r*c # 通过经纬度计算地球上两点距离的公式
distances[i,j]=distance
return distances
distances = calculate_distance(city_locations)
个体适应度
本题适应度就取目标值（巡航总距离），由于是极小化目标函数值，所以在本题中适应度（觉得别扭可以把目标函数转为极大化，那么适应值越大就越能活下来了）越低越能活下来。
# 计算个体所表示的解对应的目标函数值（适应度）,染色体的形式是0->...->i->j->k...->length_of _chromosome-1
def fitness(chromosome):
fitness=0
for i in range(length_of_chromosome-1):
fitness+=distances[chromosome[i],chromosome[i+1]] # distances[chromosome[i],chromosome[i+1]]表示染色体给出的解相邻两个巡逻点的距离
return fitness
[1]初始化种群
为了是算法尽快收敛或接近收敛，我们需要一个好的初始解，考虑使用改良圈算法生成初始种群。
注：其实随机生成初始解也可以，由于遗传算法是全局搜索的优化算法，并不必担心会陷入局部最优，只不过要靠近最优解需要更多次迭代罢了
# 初始化种群（使用改良圈算法）
def initial_population(population_size):
population=[] # 种群
for _ in range(population_size):
chromosome=np.hstack([0,np.random.permutation(np.arange(1,101)),101]) # 随机生成的个体（染色体/可行解）shape(102,)
improved=True
while improved:
improved=False
for gene1_position in range(1,length_of_chromosome-3):
for gene2_position in range(gene1_position+1,length_of_chromosome-2):
# 反转序列
new_chromosome=chromosome.copy()
new_chromosome[gene1_position:gene2_position+1]=new_chromosome[gene2_position:gene1_position-1:-1]
if fitness(new_chromosome) < fitness(chromosome):
improved=True
chromosome=new_chromosome
population.append(chromosome)
return population
[2] 交叉与变异
# 染色体变异,如果变异发生则随机选择基因突变为其他基因
def chromosome_variation(chromosome):
# 生成均匀分布（0,1）区间随机数，如果该随机数小于mutamutation_rate对应的概率恰好为mutamutation_rate于是基因变异发生
if np.any(np.random.rand(length_of_chromosome)<= mutation_rate) :
chromosome[np.where(np.random.rand(length_of_chromosome)<0.05)]=np.random.randint(1,100,(len(np.where(np.random.rand(length_of_chromosome)<0.05))))
return chromosome
# 染色体交叉互换
def chromosome_crossing(chromosome1,chromosome2):
child=chromosome1
if np.any(np.random.rand(length_of_chromosome)<= crossing_rate) :
crossing_position=np.random.randint(1,100)
child[crossing_position:]=chromosome2[crossing_position:]
return child
[3] 自然选择
# 自然选择与进化
def selection_and_evolution(new_population): # 新种群由原来的种群和子代种群构成
new_population_size=len(new_population)
fitnesses=[]
# 计算新种群个体适应度
for i in range(new_population_size):
fitnesses.append(fitness(new_population[i]))
pos=np.argsort(fitnesses)[:population_size]
return new_population[pos]
执行遗传算法
# 遗传算法
def genetic_algorithm():
population=np.array(initial_population(population_size)) #初始化种群
for _ in range(num_generations):
# 交叉互换+变异生成新种群
new_population=[]
child=[]
for i in range(population_size-1):
chromosome1,chromosome2=population[i],population[i+1]
child.append(chromosome_variation(chromosome_crossing(chromosome1,chromosome2)))
new_population=selection_and_evolution(np.vstack([population,child])) # 新种群由父代和子代自然选择后得到
# 适应度计算,选择适应度最小的染色体
fitnesses=[]
for i in range(population_size):
fitnesses.append(fitness(new_population[i]))
pos=np.argsort(fitnesses)[0]
return new_population[pos],fitnesses[pos]
代码汇总
import numpy as np
import math
from numpy.random import randint,rand
import matplotlib.pyplot as plt # 用于可视化飞行路线
population_size = 200#种群大小
num_generations = 50 # 迭代次数
mutation_rate = 0.05 # 变异率
crossing_rate = 0.8 # 染色体交叉互换率
length_of_chromosome=102 # 染色体长度
# 读取城市经纬坐标
file=r"D:\Study\my_notebook\数学建模相关\python数学实验与建模数据\程序及数据\程序及数据\17第17章 智能算法\Pdata17_1.txt"
city_locations=np.loadtxt(fname=file).reshape(-1,2)
start_point=np.array([70,40])
city_locations=np.vstack([start_point,city_locations,start_point])
# 计算巡航点之间的距离
def calculate_distance(city_locations):
r= 6371*1e3 # 地球半径m
distances=np.zeros((len(city_locations),len(city_locations))) # 101x101的实对称数组,存储i到j的距离
longitude,latitude=city_locations[:,0],city_locations[:,1]
for i in range(len(city_locations)):
for j in range(len(city_locations)):
phi1 = math.radians(latitude[i])
phi2 = math.radians(latitude[j])
delta_phi = math.radians(latitude[j] - latitude[i])
delta_lambda = math.radians(longitude[j] - longitude[i])
a = math.sin(delta_phi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2) ** 2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
distance=r*c # 通过经纬度计算地球上两点距离的公式
distances[i,j]=distance
return distances
distances = calculate_distance(city_locations)
# 计算个体所表示的解对应的目标函数值（适应度）,染色体的形式是0->...->i->j->k...->length_of _chromosome-1
def fitness(chromosome):
fitness=0
for i in range(length_of_chromosome-1):
fitness+=distances[chromosome[i],chromosome[i+1]] # distances[chromosome[i],chromosome[i+1]]表示染色体给出的解相邻两个巡逻点的距离
return fitness
# 初始化种群（使用改良圈算法）
def initial_population(population_size):
population=[] # 种群
for _ in range(population_size):
chromosome=np.hstack([0,np.random.permutation(np.arange(1,101)),101]) # 随机生成的个体（染色体/可行解）shape(102,)
improved=True
while improved:
improved=False
for gene1_position in range(1,length_of_chromosome-3):
for gene2_position in range(gene1_position+1,length_of_chromosome-2):
# 反转序列
new_chromosome=chromosome.copy()
new_chromosome[gene1_position:gene2_position+1]=new_chromosome[gene2_position:gene1_position-1:-1]
if fitness(new_chromosome) < fitness(chromosome):
improved=True
chromosome=new_chromosome
population.append(chromosome)
return population
# 染色体变异,如果变异发生则随机选择基因突变为其他基因
def chromosome_variation(chromosome):
# 生成均匀分布（0,1）区间随机数，如果该随机数小于mutamutation_rate对应的概率恰好为mutamutation_rate于是基因变异发生
if np.any(np.random.rand(length_of_chromosome)<= mutation_rate) :
chromosome[np.where(np.random.rand(length_of_chromosome)<0.05)]=np.random.randint(1,100,(len(np.where(np.random.rand(length_of_chromosome)<0.05))))
return chromosome
# 染色体交叉互换
def chromosome_crossing(chromosome1,chromosome2):
child=chromosome1
if np.any(np.random.rand(length_of_chromosome)<= crossing_rate) :
crossing_position=np.random.randint(1,100)
child[crossing_position:]=chromosome2[crossing_position:]
return child
# 自然选择与进化
def selection_and_evolution(new_population): # 新种群由原来的种群和子代种群构成
new_population_size=len(new_population)
fitnesses=[]
# 计算新种群个体适应度
for i in range(new_population_size):
fitnesses.append(fitness(new_population[i]))
pos=np.argsort(fitnesses)[:population_size]
return new_population[pos]
# 遗传算法
def genetic_algorithm():
population=np.array(initial_population(population_size)) #初始化种群
for _ in range(num_generations):
# 交叉互换+变异生成新种群
new_population=[]
child=[]
for i in range(population_size-1):
chromosome1,chromosome2=population[i],population[i+1]
child.append(chromosome_variation(chromosome_crossing(chromosome1,chromosome2)))
new_population=selection_and_evolution(np.vstack([population,child])) # 新种群由父代和子代自然选择后得到
# 适应度计算,选择适应度最小的染色体
fitnesses=[]
for i in range(population_size):
fitnesses.append(fitness(new_population[i]))
pos=np.argsort(fitnesses)[0]
return new_population[pos],fitnesses[pos]
if __name__ == "__main__":
chromosome,fitness=genetic_algorithm()
求解结果与可视化
# %matplotlib widget 命令是 Jupyter Notebook 中的一个魔法命令，用于在 notebook 中启用交互式图形界面，即使用 ipympl 后端
%matplotlib widget
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
velocity=1000 #km/h
print('最短飞行时间:{0}h'.format(fitness/(velocity*1000)),'飞行路线图:',sep='\n')
plt.scatter(city_locations[:,0],city_locations[:,1],c='c',label='Cruise point',alpha=0.7,s=15,marker='o')
plt.annotate('$BASE$',start_point,[start_point[0]-10,start_point[1]-2],arrowprops=dict(facecolor='red', arrowstyle='->'))
for i in range(len(chromosome)-1):
point1=city_locations[chromosome[i]]
point2=city_locations[chromosome[i+1]]
if i==0:
plt.plot([point1[0],point2[0]],[point1[1],point2[1]],c='r',alpha=1,linestyle=':',label='Departure direction')
elif i==(len(chromosome)-2):
plt.plot([point1[0],point2[0]],[point1[1],point2[1]],c='k',alpha=1,linestyle=':',label='Return direction')
else:
plt.plot([point1[0],point2[0]],[point1[1],point2[1]],c='b',alpha=0.4,linestyle='--')
plt.legend(loc='best')
plt.grid()
plt.title('Flight route map')
plt.xlabel('latitude')
plt.ylabel('longitude')
plt.savefig('./flight cruise.png')
plt.show()
6634
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
