遗传算法
任务描述
本关任务:利用遗传算法编写一个解决 TSP 问题的程序。
相关知识
为了完成本关任务,你需要掌握遗传算法及其步骤。
遗传算法简介
能够生存下来的往往不是最强大的物种,也不是最聪明的物种,而是最能适应环境的物种。 你也许在想:这句话和遗传算法有什么关系?其实遗传算法的整个概念就基于这句话。让我们用一个简单的例子来解释 :
我们先假设一个情景,现在你是一国之王,为了让你的国家免于灾祸,你实施了一套法案:
-
你选出所有的好人,要求其通过生育来扩大国民数量。
-
这个过程持续进行了几代。
-
你将发现,你已经有了一整群的好人。
这个例子虽然不太可能,但是我用它是想帮助你理解概念。也就是说,我们改变了输入值(比如:人口),就可以获得更好的输出值(比如:更好的国家)。现在,我假定你已经对这个概念有了大致理解,认为遗传算法的含义应该和生物学有关系。那么我们就快速地看一些小概念,这样便可以将其联系起来理解。 接着我们就可以理解遗传算法的概念: 让我们回到前面讨论的那个例子,并总结一下我们做过的事情。
-
首先,我们设定好了国民的初始人群大小。
-
然后,我们定义了一个函数,用它来区分好人和坏人。
-
再次,我们选择出好人,并让他们繁殖自己的后代。
-
最后,这些后代们从原来的国民中替代了部分坏人,并不断重复这一过程。
遗传算法实际上就是这样工作的,也就是说,它基本上尽力地在某种程度上模拟进化的过程。 总的来说,遗传算法是模拟自然界物种进化的自然选择和遗传学机理的生物进化的计算模型,它通过模拟生物进化来搜索最优解的方法,在机器学习、组合优化等方面有广泛的用途。其基本流程如下:
图1 遗传算法流程图
接着我们需要了解一些杂乱的概念: 在生物学上,种群是在一定空间范围内同时生活着的同种生物的全部个体。而在遗传算法中,个体通常为某个问题的一个解,而把一组一组的可能解的个体的集合就叫做种群。 生物的 DNA 有四种碱基对,分别是 ACGT , DNA 的编码可以看作是 DNA 上碱基对的不同排列,不同的排列使得基因的表现出来的性状也不同(如单眼皮双眼皮)。在计算机中,我们可以模仿这种编码,但是碱基对的种类只有两种,分别是 0 , 1 。只要我们能够将不同的实数表示成不同的 0 ,1二进制串表示,就完成了编码。也就是说其实我们并不需要去了解一个实数对应的二进制具体是多少,我们只需要保证有一个映射能够将十进制的数编码为二进制即可,至于这个映射是什么,其实可以不必关心。将个体(可能解)编码后的二进制串叫做染色体,染色体(或者有人叫 DNA)就是个体(可能解)的二进制编码表示。而解码就是一个逆向的过程。 接着我们详细的讨论一下这些步骤。
种群初始化
在初期选择时,我们首先要选择一些比较优秀的个体,然后在利用一些近似算法比如改良圈算法得到一个较好的种群,在 TSP 问题中,其算法思想是先任意交换两个城市的顺序位置,如果总距离减少,则更新改变染色体,如此重复直至不能再修改。 本实训中初始化进行如下操作:设置进化计数器 t=0,设置最大进化数 T ,随机生成 M 个个体作为初始种群 P(0) ,然后再进行一个初次改良。
适应度和选择
我们已经得到了一个种群,现在要根据适者生存规则把优秀的个体保存下来,同时淘汰掉那些不适应环境的个体。现在摆在我们面前的问题是如何评价一个个体对环境的适应度?在我们的求 TSP 问题中,可以直接用可能解(个体)对应的距离的大小来评估。 接着就是选择操作,在这里,我们首先选择 20% 适应度强的个体直接生存(适者生存);其次使用轮盘赌法选择一定数量的个体( K )。根据轮盘赌的过程可知,染色体的适应度越大,则其被选中的概率就越大。 对于轮盘选择算法,又称比例选择方法。其基本思想是,各个个体被选中的概率与其适应度大小成正比,基本操作如下: (1)计算出群体中每个个体的适应度f(i=1,2,…,M),M 为群体大小; (2)计算出每个个体被遗传到下一代群体中的概率; (3)计算出每个个体的累积概率; (q[i]称为染色体x[i] (i=1, 2, …, n)的积累概率) (4)在[0,1]区间内,产生一个均匀分布的伪随机数 r; (5)若 r<q[1],则选择个体 1,否则,选择个体 k,使得:q[k-1]<r≤q[k] 成立; (6)重复(4)、(5)共 M 次。 在本实训中,我们使用 numpy 中的 np.random.choice( )可以达到一样的效果。
交叉
交叉是指每一个个体是由父亲和母亲两个个体繁殖产生,子代个体的 DNA (二进制串)获得了一半父亲的 DNA ,一半母亲的 DNA ,但是这里的一半并不是真正的一半,这个位置叫做交配点,是随机产生的,可以是染色体的任意位置。通过交叉子代获得了一半来自父亲一半来自母亲的 DNA ,但是子代自身可能发生变异,使得其 DNA 既不来自父亲,也不来自母亲,在某个位置上发生随机改变,通常就是改变 DNA 的一个二进制位( 0 变到 1 ,或者 1 变到 0 )。交叉和变异是遗传算法的核心步骤,也是遗传算法的难点。他们是实现群体多样性的一种手段,同时也是全局寻优的保证。 交叉是一个种群中,不同个体的部分信息交换。传统的交叉有单点交叉、双点交叉以及部分匹配交叉。为了避免城市重复和城市遗漏等问题,本文采用次序交叉法。 次序杂交算法首先随机地在双亲中选择两个杂交点,再交换杂交段,其他位置根据双亲城市的相对位置确定。 例如:本例子中分为两个阶段。 第一阶段:随机选择两条染色体,再选择两个杂交点。 随机交叉点索引:4,7 p1:0,8,5,4,7,3,1,6,9,2 p2:1,5,4,6,9,0,3,2,8,7 最终子代: P1:&,&,&,&,9,0,3,2,&,& P2:&,&,&,&,7,3,1,6,&,& 第二阶段:转换父代染色体 p1 和 p2 顺序得到亚父代。 亚父代: p3:9,2,0,8,5,4,7,3,1,6 p4:8,7,1,5,4,6,9,0,3,2 去掉交叉片段(p1 应该去掉 p2 的片段 9-0-3-2,p2 同理)后得到: p5:8,5,4,7,1,6 p6:8,5,4,9,0,2 最后,第二个杂交点开始填入,得到子代: P1:4,7,1,6,9,0,3,2,8,5 P2:4,9,0,2,7,3,1,6,8,5
变异
变异是对种群中选定个体的操作,与交叉不同。 具体设计如下,按照给定的变异率,对选定变异的个体,在染色体随机地取三个整数,满足 1<u<v<w<32 ,把 v、u 之间(包括 u 和 v)的基因段插到 w 后面。 例如,对于 p1,随机选择三个整数:1,3,5。 p1:0,8,5,4,7,3,1,6,9,2 变异后: P1:0,7,3,8,5,4,1,6,9,2
编程要求
根据提示,本关的编程任务是补全右侧代码片段 Begin 至 End 中间的代码,计算并分别输出 TSP 问题的迭代次数,最优值以及最优路径( TSP 算法各城市的距离在代码中已给出)。 参数设置如下:
种群规模:200
强者生存率:0.2
变异概率:0.1
迭代次数:200
改良次数:500
测试说明
平台会对你编写的代码进行测试, 最优值低于 170 即可通过测试。
预期输出1:
迭代次数: 200
最优值是: 157.925959920418
最优路径: [5, 23, 30, 29, 21, 15, 14, 16, 22, 17, 18, 19, 20, 12,
8, 7, 6, 4, 10, 9, 13, 25, 3, 24, 11, 2, 1, 28, 26, 0, 27]
路径长度符合要求
预期输出2:
迭代次数: 200
最优值是: 158.90947572684246
最优路径: [1, 2, 3, 25, 13, 12, 9, 10, 4, 6, 7, 8, 14, 15, 16, 22,
17, 18, 20, 19, 21, 29, 30, 23, 5, 28, 27, 0, 26, 24, 11]
路径长度符合要求
import numpy as np
import math
import random
import time
start = time.time()
# 31个城市的坐标
city_condition=[[106.54,29.59],[ 91.11,29.97],[ 87.68,43.77],[106.27,38.47],[111.65,40.82],
[108.33,22.84],[126.63,45.75],[125.35,43.88],[123.38,41.8 ],[114.48,38.03],[112.53,37.87],
[101.74,36.56],[117.0,36.65],[113.6,34.76],[118.78,32.04],[117.27,31.86],
[120.19,30.26],[119.3,26.08],[115.89,28.68],[113.0,28.21],[114.31,30.52],
[113.23,23.16],[121.5,25.05],[110.35,20.02],[103.73,36.03],[108.95,34.27],
[104.06,30.67],[106.71,26.57],[102.73,25.04],[114.1,22.2 ],[113.33,22.13]]
# 距离矩阵
city_count = 31
Distance = np.zeros([city_count, city_count])
for i in range(city_count):
for j in range(city_count):
Distance[i][j] = math.sqrt(
(city_condition[i][0] - city_condition[j][0]) ** 2 + (city_condition[i][1] - city_condition[j][1]) ** 2)
# 种群数
count = 200
# 改良次数
improve_count = 500
# 进化次数
iteration = 200
# 设置强者的定义概率,即种群前20%为强者
retain_rate = 0.2
# 变异率
mutation_rate = 0.1
# 设置起点
index = [i for i in range(city_count)]
#总距离
def get_total_distance(path_new):
distance = 0
for i in range(city_count - 1):
# count为30,意味着回到了开始的点,此时的值应该为0.
distance += Distance[int(path_new[i])][int(path_new[i + 1])]
distance += Distance[int(path_new[-1])][int(path_new[0])]
return distance
# 改良
#思想:随机生成两个城市,任意交换两个城市的位置,如果总距离减少,就改变染色体。
def improve(x):
i = 0
distance = get_total_distance(x)
while i < improve_count:
u = random.randint(0, len(x) - 1)
v = random.randint(0, len(x) - 1)
if u != v:
new_x = x.copy()
## 随机交叉两个点,t为中间数
t = new_x[u]
new_x[u] = new_x[v]
new_x[v] = t
new_distance = get_total_distance(new_x)
if new_distance < distance:
distance = new_distance
x = new_x.copy()
else:
continue
i += 1
# 适应度评估,选择,迭代一次选择一次
def selection(population):
# 对总距离从小到大进行排序
graded = [[get_total_distance(x), x] for x in population]
graded = [x[1] for x in sorted(graded)]
# 选出适应性强的染色体
# ********** Begin **********#
retain_length = int(len(graded) * retain_rate)
# ********** End **********#
#适应度强的集合,直接加入选择中
# ********** Begin **********#
parents = graded[:retain_length]
# ********** End **********#
## 轮盘赌算法选出K个适应性不强的个体,保证种群的多样性
s = graded[retain_length:]
# 挑选的不强的个数
k = count * 0.2
# 存储适应度
a = []
for i in range(0, len(s)):
a.append(get_total_distance(s[i]))
sum = np.sum(a)
b = np.cumsum(a / sum)
while k > 0: # 迭代一次选择k条染色体
t = random.random()
for h in range(1, len(b)):
if b[h - 1] < t <= b[h]:
parents.append(s[h])
k -= 1
break
return parents
# 交叉繁殖
def crossover(parents):
# 生成子代的个数,以此保证种群稳定
target_count = count - len(parents)
# 孩子列表
children = []
while len(children) < target_count:
male_index = random.randint(0, len(parents) - 1)
female_index = random.randint(0, len(parents) - 1)
#在适应度强的中间选择父母染色体
if male_index != female_index:
male = parents[male_index]
female = parents[female_index]
left = random.randint(0, len(male) - 2)
right = random.randint(left + 1, len(male) - 1)
# 交叉片段
gene1 = male[left:right]
gene2 = female[left:right]
#得到原序列通过改变序列的染色体,并复制出来备用。
child1_c = male[right:] + male[:right]
child2_c = female[right:] + female[:right]
child1 = child1_c.copy()
child2 = child2_c.copy()
#已经改变的序列=>去掉交叉片段后的序列
for o in gene2:
child1_c.remove(o)
for o in gene1:
child2_c.remove(o)
#交换交叉片段
# ********** Begin **********#
child1[left:right] = gene2
child2[left:right] = gene1
child1[right:] = child1_c[0:len(child1) - right]
child1[:left] = child1_c[len(child1) - right:]
child2[right:] = child2_c[0:len(child1) - right]
child2[:left] = child2_c[len(child1) - right:]
children.append(child1)
children.append(child2)
# ********** End **********#
return children
# 变异
def mutation(children):
#children现在包括交叉和优质的染色体
for i in range(len(children)):
if random.random() < mutation_rate:
child = children[i]
#产生随机数
u = random.randint(0, len(child) - 4)
v = random.randint(u + 1, len(child) - 3)
w = random.randint(v + 1, len(child) - 2)
child = child[0:u] + child[v:w] + child[u:v] + child[w:]
children[i] = child
return children
# 得到最佳纯输出结果
def get_result(population):
graded = [[get_total_distance(x), x] for x in population]
graded = sorted(graded)
return graded[0][0], graded[0][1]
模拟退火算法
任务描述
本关任务:利用模拟退火算法解决 TSP 问题。 TSP 问题是在给定的一些城市和这些城市之间的距离中,找出访问每一座城市一次并最终回到起始城市的最短回路。
相关知识
为了完成本关任务,你需要掌握模拟退火算法。
退火
在冶金学中,金属通常要进行原子重排,这是在退火过程中实现的。金属中的原子按照局部能量最小化进行排列。为了以较低的能量重新排列这些原子,首先需将金属加热,直到液化,然后将熔融的金属缓慢冷却,直到凝固。退火后的金属表现出许多人们期望的性能,例如它的韧性和硬度都得到了提升。
图1 金属退火过程
模拟退火
如图 2 所示:假设这是一个求最小值的函数图像,如果采用贪心策略,那么从 A 点开始试探,如果函数值继续减少,那么试探过程就会继续。而当到达点 B 时,显然我们的探求过程就结束了(因为无论朝哪个方向努力,结果只会越来越大)。最终我们只能找到一个局部最优解 B 。
图2 求最小值的函数图像
模拟退火其实也是一种贪心算法,但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解,因此有可能会跳出这个局部的最优解,达到全局的最优解。以上图为例,模拟退火算法在搜索到局部最优解 B 后,会以一定的概率接受向右继续移动。也许经过几次这样的不是局部最优的移动后,会到达 B 和 C 之间的峰点,于是就跳出了局部最小值 B 。
模拟退火算法——概率取值
刚刚我们说了,模拟退火算法是以一定的概率来接受比当前解要差的值,那么这个概率是多少呢?由于我们之前说的模拟退火算法是源自于金属退火的原理,所以这个概率的计算也是来源于它。根据 Metropolis 准则,粒子在温度 T 时,趋于平衡的概率为 exp(−ΔE/(kT)),其中 E 为温度 T 时的内能,ΔE 为其改变数,k 为 Boltzmann 常数。Metropolis 准则常表示为:
Metropolis 准则表明,在温度为 T 时,出现能量差为 dE 的降温的概率为 P(dE),表示为:P(dE)=exp(dE/(kT))。其中 k 是一个常数,exp 表示自然指数,且 dE<0 ,所以 P 和 T 正相关。这条公式就表示:温度越高,出现一次能量差为 dE 的降温的概率就越大;温度越低,则出现降温的概率就越小。又由于 dE 总是小于0(因为退火的过程是温度逐渐下降的过程),因此 dE/kT<0,所以 P(dE) 的函数取值范围是 (0,1)。随着温度 T 的降低, P(dE) 会逐渐降低。应用到实际优化问题中,我们将内能 E 看作目标函数值 F,温度 T 代表控制参数 t。控制参数 t 控制函数的迭代次数,迭代期间不断地进行“产生新解-计算差值-Metropolis 准则判断是否接受新解 -t衰减”的循环。
那么 Metropolis 准则可以转化如下:
简单地来说,就是判断当前的解是否比原来的解好,如果新解更好,那么以 1 的概率接受它;如果比它差,那么以 p 的概率接受它,p 的取值与当前解、新解与控制参数有关。
模拟退火算法步骤
伪代码如下:
/*
* J(y):在状态y时的评价函数值
* Y(i):表示当前状态
* Y(i+1):表示新的状态
* r: 用于控制降温的快慢
* T: 系统的温度,系统初始应该要处于一个高温的状态
* T_min :温度的下限,若温度T达到T_min,则停止搜索
*/
while( T > T_min )
{
dE = J( Y(i+1) ) - J( Y(i) ) ;
if ( dE >=0 ) //表达移动后得到更优解,则总是接受移动
Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
else
{
// 函数exp( dE/T )的取值范围是(0,1) ,dE/T越大,则exp( dE/T )也
if ( exp( -dE/T ) > random( 0 , 1 ) )
Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
}
T = r * T ; //降温退火 ,0<r<1 。r越大,降温越慢;r越小,降温越快
/*
* 若r过大,则搜索到全局最优解的可能会较高,但搜索的过程也就较长。若r过小,则搜索的过程会很快,但最终可能会达到一个局部最优值
*/
i ++ ;
}
编程要求
根据提示,本关的编程任务是补全右侧编辑器 Begin 至 End 中间的代码,saa()函数的任务是利用模拟退火算法求出 TSP 问题的最短路径。 注意:saa()函数返回三个参数,分别为 ans 、 result 、 cnt,其中 ans 为最短路径,result 为路径长度,cnt 为降温次数( TSP 各城市的距离在代码中已给出)。
测试说明
平台将自动编译补全后的代码,并生成若干组测试数据,接着根据程序的输出判断程序是否正确,其中路径长度只需低于 50000 并且降温次数为 1448 即可通过实训。
以下是平台的测试样例:
预期输出1:
路径如下:[14, 13, 19, 18, 15, 12, 9, 28, 7, 17, 27, 11, 26, 20, 6,
16, 30, 2, 21, 25, 22, 8, 0, 10, 29, 4, 5, 23, 3, 1, 24]
路径长度:46860.89447968355
降温次数:1448
路径长度符合要求
预期输出2:
路径如下:[6, 27, 29, 17, 3, 4, 24, 14, 15, 18, 30, 11, 12, 7, 9,
2, 19, 21, 25, 22, 28, 23, 26, 8, 16, 10, 5, 1, 0, 13, 20]
路径长度:44583.241163416074
降温次数:1448
路径长度符合要求
import math
import time
import random
# 初始温度
T = 50000
# 最低温度
T_end = 1e-8
# 在每个温度下的迭代次数
L = 100
# 退火系数
delta = 0.98
# 31个城市的坐标
citys = [[1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535], [3326, 1556], [3238, 1229], [4196, 1004],
[4312, 790], [4386, 570], [3007, 1970], [2562, 1756], [2788, 1491], [2381, 1676], [1332, 695], [3715, 1678],
[3918, 2179], [4061, 2370], [3780, 2212], [3676, 2578], [4029, 2838], [4263, 2931], [3429, 1908], [3507, 2367],
[3394, 2643], [3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826], [2370, 2975]]
# 存储两个城市之间的距离
d = [[0 for i in range(31)] for j in range(31)]
# 存储一条路径
ans = []
# 计算降温次数
cnt = 0
# 计算两个城市之间的距离
def get_city_distance():
for i in range(len(citys)):
for j in range(i, len(citys)):
d[i][j] = d[j][i] = math.sqrt((citys[i][0] - citys[j][0]) ** 2 + (citys[i][1] - citys[j][1]) ** 2)
# 使用随机交换路径中两个城市的位置来产生一条新路径
def create_new(a):
#a=list(range(len(a)))
i = random.randint(0, len(a) - 1)
j = random.randint(0, len(a) - 1)
a = list(a)
a[i], a[j] = a[j], a[i]
return a
# 获取路径的长度
def get_route_distance(a):
dist = 0
for i in range(len(a) - 1):
dist += d[a[i]][a[i + 1]]
return dist
#模拟退火
def saa():
# ********** Begin **********#
get_city_distance()
cnt=0
ans=range(0,len(citys))
t=T
result=0
while t>=T_end:
for i in range(0,L):
ans_new = create_new(ans)
d1,d2=get_route_distance(ans),get_route_distance(ans_new)
de=d2-d1
result=d1
if de<0:
ans=ans_new
result=d2
else:
if (math.e ** (-de/T)>random.random()):
ans=ans_new
result=d2
t=t * delta
cnt+=1
# ********** End **********#
#函数返回三个参数
#ans:路径,例如:[14, 13, 19, 18, 15, 12, 9, 28, 7, 17, 27, 11, 26, 20, 6, 16, 30, 2, 21, 25, 22, 8, 0, 10, 29, 4, 5, 23, 3, 1, 24]
#result:路径长度
#cnt:降温次数
#路径长度在50000以内即可,降温次数最终返回值应为1448
return ans,result,cnt