遗传算法 python 简书_基于DEAP库的Python进化算法从入门到入土--（三）遗传算法求解TSP问题...

前言

在用基于DEAP设计的遗传算法求解函数极值后，我们想要进一步解决一些更加困难点的问题。TSP问题就是很好的实验算法的舞台。本文将会介绍以下内容：

TSP问题的定义与运用遗传算法求解

对比遗传算法与贪心算法及其改进型在用时和精确度两方面的表现

探讨可能的改进方向

TSP问题的定义与GA求解

TSP问题简介

旅行商问题(Traveling Salesman Problem(TSP))描述：

给出

个不同的城市

以及城市两两之间的距离

，求一个商人从某个城市出发，走遍所有城市并回到起点所需要经过的路径总距离为多少？

即求解一个路径(也即城市的排列)

，使得

取最小值

如果对两个城市

有

，那该问题被称为对称旅行商问题(Symmetric Traveling Salesman Problem(STSP))。否则就是非对称旅行商问题(ATSP)。

尽管描述很简单，TSP问题是一个非常经典而重要的问题：

自从1930年代 Merrill M.Flood在学界引起了对该问题的重视以来，学者们已经从图论、组合优化、线性规划等不同角度对该问题进行了深入广泛的研究。

TSP问题的一个特殊性在于它必然存在可以找到的最优解(通过穷举法一定能在

时间内找到该最优解)，但是在计算复杂度方面它是一个经典的NP-hard问题。因而快速有效寻找TSP的最优解仍然是个有吸引力的问题。

在实际应用方面，它可以被视为物流调度，交通规划，芯片布线优化等等问题的抽象简化。

TSP问题的遗传算法求解

通过之前对DEAP的解读，我们很容易就能基于该框架，写出一个求解TSP问题的GA算法。在TSPLib上测试该算法之前，我们先在随机生成的数据集上进行测试与改进：

## 环境设置

import matplotlib.pyplot as plt

import numpy as np

%matplotlib inline

from scipy.spatial import distance

from deap import creator, base, tools, algorithms

import random

params = {

'font.family': 'serif',

'figure.figsize': [4.0, 3.0],

'figure.dpi': 300,

'savefig.dpi': 300,

'font.size': 12,

'legend.fontsize': 'small'

}

plt.rcParams.update(params)

#--------------------------------

## 定义TSP中的基本元素

# 用[n,2]的np.array存储城市坐标；每行存储一个城市

def genCity(n, Lb = 100 ,Ub = 999):

# 生成城市坐标

# 输入：n -- 需要生成的城市数量

# 输出: nx2 np array 每行是一个城市的[X,Y]坐标

np.random.seed(42) # 保证结果的可复现性

return np.random.randint(low = Lb, high = Ub, size=(n,2))

# 计算并存储城市距离矩阵

def cityDistance(cities):

# 生成城市距离矩阵 distMat[A,B] = distMat[B,A]表示城市A，B之间距离

# 输入：cities -- [n,2] np array， 表示城市坐标

# 输出：nxn np array， 存储城市两两之间的距离

return distance.cdist(cities, cities, 'euclidean')

def completeRoute(individual):

# 序列编码时，缺少最后一段回到原点的线段

return individual + [individual[0]] # 不要用append

# 计算给定路线的长度

def routeDistance(route):

# 输入：

# route -- 一条路线，一个sequence

# 输出：routeDist -- scalar，路线的长度

if route[0] != route[-1]:

route = completeRoute(route)

routeDist = 0

for i,j in zip(route[0::],route[1::]):

routeDist += cityDist[i,j] # 这里直接从cityDist变量中取值了，其实并不是很安全的写法，单纯偷懒了

return (routeDist), # 注意DEAP要求评价函数返回一个元组

# 路径可视化

def plotTour(tour, cities, style = 'bo-'):

if len(tour)>1000: plt.figure(figsize = (15,10))

start = tour[0:1]

for i,j in zip(tour[0::], tour[1::]):

plt.plot([cities[i,0],cities[j,0]], [cities[i,1],cities[j,1]], style)

plt.plot(cities[start,0],cities[start,1],'rD')

plt.axis('scaled')

plt.axis('off')

#--------------------------------

## 设计GA算法

nCities = 30

cities = genCity(nCities) # 随机生成nCities个城市坐标

# 问题定义

creator.create('FitnessMin', base.Fitness, weights=(-1.0,)) # 最小化问题

creator.create('Individual', list, fitness=creator.FitnessMin)

# 定义个体编码

toolbox = base.Toolbox()

toolbox.register('indices', random.sample, range(nCities), nCities) # 创建序列

toolbox.register('individual', tools.initIterate, creator.Individual, toolbox.indices)

# 生成族群

N_POP = 100

toolbox.register('population', tools.initRepeat, list, toolbox.individual)

pop = toolbox.population(N_POP)

# 注册所需工具

cityDist = cityDistance(cities)

toolbox.register('evaluate', routeDistance)

toolbox.register('select', tools.selTournament, tournsize = 2)

toolbox.register('mate', tools.cxOrdered)

toolbox.register('mutate', tools.mutShuffleIndexes, indpb = 0.2)

# 数据记录

stats = tools.Statistics(key=lambda ind: ind.fitness.values)

stats.register('avg', np.mean)

stats.register('min', np.min)

# 调用内置的进化算法

resultPop, logbook = algorithms.eaSimple(pop, toolbox, cxpb = 0.5, mutpb = 0.2, ngen = 200, stats = stats, verbose = True)

本次运行的结果为：

gen nevals avg min

0 100 14201 11314.5

1 58 13675.4 11505.4

2 58 13410.6 11300.1

3 60 13174.9 11064.3

4 57 12962.4 10907.4

5 63 12796.9 10470.6

6 66 12585.8 10644.2

7 64 12245.2 10232

8 68 12067 9592.99

9 62 12059.6 9592.99

10 70 12206.6 9592.99

11 57 11908.6 9496.14

12 42 11468.3 9496.14

13 60 11393.1 9062.07

14 47 11109.7 9062.07

15 63 11123.8 9062.07

16 58 10923.4 8926.45

17 61 10995.1 8545.11

18 58 11174.7 9062.07

19 76 11091.9 9062.07

20 53 10993.1 8888.39

21 58 11016.3 8842.55

22 52 10914.2 8842.55

23 47 10745.8 8659.72

24 54 10526 8659.72

25 60 10870.4 8717.75

26 56 11078.5 8717.75

27 66 10904.5 8717.75

28 55 10796.9 8717.75

29 58 10840.5 8568.54

30 70 10783.7 8714.81

31 75 10945.5 8714.81

32 62 10906.3 8714.81

33 55 10883.5 8714.81

34 61 10768.4 8717.75

35 64 10861.7 9045.8

36 72 10839 9097.6

37 73 10819.7 8860.66

38 70 10528.5 8736.59

39 55 10250.1 8736.59

40 67 10563.3 8513.88

41 50 10602.1 8513.88

42 59 10518 8513.88

43 61 10596.9 8931.09

44 57 10780 9109.96

45 67 10757.6 8771.09

46 60 10504.8 8771.09

47 63 10421.2 8771.09

48 69 10526.5 8771.09

49 58 10442.3 8691.01

50 65 10535.4 8771.09

51 62 10683 8771.09

52 68 10532.6 8804.67

53 54 10488.4 8804.67

54 62 10583.1 8408.69

55 53 10514.8 8833.92

56 57 10524.8 9052.01

57 65 10599.3 8732.31

58 58 10227.1 8732.31

59 62 10358.2 8732.31

60 66 10532 8730.2

61 51 10348.5 8671.23

62 67 10542 8671.23

63 60 10584.6 8663.02

64 57 10417.8 8522.64

65 58 10552.6 8522.64

66 61 10413.8 8367.42

67 66 10422.5 8671.23

68 57 10559.6 8671.23

69 58 10263.6 8285.6

70 60 10463.1 8285.6

71 65 10499 8285.6

72 55 10330.4 8529.22

73 67 10350.1 8576.73

74 73 10271 7570.84

75 51 10252.9 8289.29

76 60 10343.2 8112.17

77 61 9767.09 7565.35

78 57 9751.91 7565.35

79 50 9684.03 7565.35

80 53 9808.93 7565.35

81 59 9509.4 7565.35

82 58 9758.84 7565.35

83 65 9852.58 6437.84

84 56 9866.73 6437.84

85 55 9572.01 6298.22

86 61 9217.73 6298.22

87 63 9311.42 6298.22

88 61 9057.95 6298.22

89 63 9103.11 6298.22

90 66 9116.68 6298.22

91 51 9024.19 6454.46

92 45 8683.05 6454.46

93 64 8584.42 6709.33

94 59 8942 6709.33

95 64 9082.05 6709.33

96 71 9037.46 6709.33

97 60 9154.49 6933.71

98 67 9384.12 6933.71

99 75 9682.95 7115.49

100 61 9300.16 7266.03

101 50 9035.16 7266.03

102 52 8740.91 6960.39

103 62 9063.54 7145.63

104 62 9220.89 7252.07

105 57 9210.75 6956.42

106 52 9228 6956.42

107 62 9095.18 6956.42

108 55 8764.35 6956.42

109 65 8961.6 6956.42

110 62 9029.89 6911.78

111 56 8757.81 6558.41

112 75 8436.12 6261.76

113 48 8489.56 6261.76

114 56 8379.92 6261.76

115 50 8567.57 6348.63

116 58 8505.19 6405.25

117 56 8485.4 5930.98

118 58 8689.05 5930.98

119 68 8504.73 6167.82

120 70 8123.5 6167.82

121 58 7992.72 6128.46

122 60 8441.54 6128.46

123 59 8270.66 6128.46

124 58 7944.21 5924.28

125 61 8259.92 5924.28

126 68 7954.45 5924.28

127 63 7913.86 6063.51

128 62 8507.67 6063.51

129 52 8079.8 6063.51

130 62 8180.93 5693.15

131 79 8299 6063.51

132 60 8070.93 6255.88

133 70 8149.17 6048.78

134 54 8125.64 6268.16

135 65 8161.05 6080.82

136 62 8093.43 6221.13

137 54 8046.15 6221.13

138 67 8064.43 5715.42

139 57 8108.32 5447.17

140 63 7910.74 5362.81

141 56 8040.57 5362.81

142 61 7936.12 5362.81

143 59 7951.14 5362.81

144 60 8222.17 5362.81

145 70 7815.39 5362.81

146 57 7660.78 5362.81

147 61 7545.06 5362.81

148 62 7619.84 5362.81

149 51 7387.11 5362.81

150 66 7396.11 5362.81

151 57 7224.2 5362.81

152 62 7062.38 5280.67

153 60 7450.55 5280.67

154 57 7058.69 5280.67

155 51 7452.15 5218.92

156 64 7332.41 5280.67

157 55 7261.36 5280.67

158 53 6949.59 5280.67

159 51 7034.86 5280.67

160 57 6985.01 5280.67

161 67 6554.21 5280.67

162 62 6671.5 5280.67

163 61 6938.98 5121.92

164 57 7247.23 5121.92

165 69 7018.4 5039.78

166 67 6764.25 5039.78

167 61 7172.14 5039.78

168 58 7172.76 5039.78

169 48 6760.36 5039.78

170 66 6591.99 5039.78

171 59 7039.59 5039.78

172 51 6917.34 5039.78

173 60 6631.73 5039.78

174 58 6683.84 4997.91

175 65 6559.69 4997.91

176 69 6908.49 4910.62

177 61 7016.93 4910.62

178 60 6772.34 4910.62

179 58 6492.51 4757.57

180 67 6389.65 4757.57

181 49 6018.27 4757.57

182 55 6364.89 4757.57

183 68 6781.2 4757.57

184 63 6500.91 4757.57

185 60 6451.73 4757.57

186 73 7277.07 4757.57

187 57 6232.07 4757.57

188 74 6360.29 4757.57

189 60 6608.44 4757.57

190 52 6222.14 4757.57

191 51 6482.11 4757.57

192 52 6118.93 4732.05

193 57 6170.32 4732.05

194 65 6303.12 4732.05

195 52 6060.32 4732.05

196 57 6133.72 4757.57

197 60 6120.28 4757.57

198 66 5602.95 4757.57

199 68 6006.72 4757.57

200 60 6057.53 4757.57

对结果进行可视化：

tour = tools.selBest(resultPop, k=1)[0]

tourDist = tour.fitness

tour = completeRoute(tour)

print('最优路径为:'+str(tour))

print('最优路径距离为：'+str(tourDist))

plotTour(tour, cities)

可视化结果：

最优路径为:[7, 15, 0, 14, 12, 6, 17, 25, 8, 27, 21, 20, 23, 24, 18, 13, 11, 1, 10, 29, 22, 9, 3, 28, 4, 5, 19, 26, 16, 2, 7]

最优路径距离为：(4757.574098243668,)

从肉眼可见的角度上，这个解距离最优解还是有差距的。

30cities_SGA

对训练过程可视化如下：

fitnessVSiterations_SGA

算法对比与分析

nnTSP + 2OPT

最近邻算法：nn TSP

穷举法能够保证求出TSP问题的最优解，但是其耗时会随着问题规模的扩大而剧增。假设求解一个规模为10的TSP问题需要1秒，那么可以估计对于更大规模问题，穷举法所需要的时间：

问题规模(城市数量)

求解时间

10

1 秒

11

11 秒

12

2.2 分

13

28.6 分

14

6.7 小时

15

100 小时

16

1601 小时

而在应用于实际时，TSP问题的规模往往成千上万，用穷举法求解显然不现实，因此我们需要更加有效的求解方法。

在考虑问题时，一般算法设计的策略有以下几种:

暴力破解：列出所有可能性，一一尝试并最终给出答案；

近似算法：如果问题的最优解或准确解很难给出，那么可以尝试找出问题的次优解或近似解；

贪心算法：对于多步骤算法，可以在每一步都取局部最优，并重复直到问题解决；

迭代算法：用一个算法求解，另一个算法来提升解的有效性；

集合算法：在算法集中准备一系列算法，依次尝试并取最好结果；

分而治之：讲一个问题划分为一系列子问题，依次求解后组合成最后结果；

站在巨人肩膀上：参考已有算法(just google it!)。

最近邻算法就是同时包含了近似算法和贪心算法思想的一种算法。用于求解TSP问题时，nnTSP可以描述为：

从某个城市出发；

移动到离现在所处城市最近的且没有访问过的邻居去；

其实现还是非常简单的：

# 在未访问列表中寻找当前位置的最近邻

def nearestNeighbor(cityDist, currentPos, unvisited):

# 输入：cityDist -- [n,n] np array,记录城市间距离

# currentPos -- 一个数字，指示当前位于的城市

# unvisited -- 一个列表，指示可选邻居列表

# 输出：

# nextToVisit -- 一个数字，给出最近邻的index

# dist -- 一个数字，给出到最近邻的距离

neighborDist = cityDist[currentPos,unvisited]

neighborIdx = np.argmin(neighborDist)

nextToVisit = unvisited[neighborIdx]

dist = neighborDist[neighborIdx]

return nextToVisit, dist

# 用贪婪算法求解TSP

def nnTSP(cities, start = 0):

cityList = list(range(cities.shape[0]))

tour = [start]

unvisited = cityList.copy()

unvisited.remove(start)

currentPos = start

tourDist = 0

cityDist = cityDistance(cities)

while unvisited:

# 当unvisited集合不为空时，重复循环

# 找到当前位置在unvisited中的最近邻

nextToVisit, dist = nearestNeighbor(cityDist, currentPos, unvisited)

tourDist += dist

currentPos = nextToVisit

tour.append(nextToVisit)

unvisited.remove(nextToVisit)

# 重新回到起点

tour.append(start)

tourDist += cityDist[currentPos, start]

return tour,tourDist

在上面相同的问题上应用nnTSP求解：

tour,tourDist = nnTSP(cities)

print('nnTSP寻找到的最优路径为：' + str(tour))

print('nnTSP寻找到最优路径长度为：' + str(tourDist))

plotTour(tour, cities)

结果为：

nnTSP寻找到的最优路径为：[0, 15, 7, 14, 2, 16, 19, 26, 5, 4, 9, 28, 3, 22, 29, 1, 11, 13, 10, 12, 6, 21, 20, 27, 25, 17, 8, 23, 18, 24, 0]

nnTSP寻找到最优路径长度为：5156.612343174932

30cities_nnTSP

repNNTSP

从nnTSP的算法设计上可以看到，该算法求出的最优路径是和起点相关的，对于不同的起点，会产生不同的最优路径。因而一个简单的改进思路是我们以每个城市为起点进行一次nnTSP搜索，然后在所有得到的路线中选取最优的。

代码实现如下：

# 重复nnTSP但是每次用不同的起点，挑选最佳的结果

def repNNTSP(cities):

optimizedRoute = [] # 最优路径

minDistance = np.Inf # 最优路径长度

for i in range(len(cities)):

tour,tourDist = nnTSP(cities, start = i)

if tourDist < minDistance:

optimizedRoute = tour

minDistance = tourDist

return optimizedRoute, minDistance

结果为：

repNNTSP寻找到的最优路径为：[20, 21, 27, 25, 17, 8, 15, 0, 7, 14, 2, 16, 19, 26, 5, 4, 9, 28, 3, 22, 29, 1, 11, 13, 10, 12, 6, 23, 18, 24, 20]

repNNTSP寻找到的最优路径长度为：4529.115576833992

30cities_repnnTSPtour

如预期一般，这个结果比nnTSP更优秀，而且感觉上距离最优解已经很接近了。

2-OPT

在nnTSP给出的结果中，包含了一些交叉。这些交叉会导致路径长度的增加：

cross

可以看到当出现X交叉时，从三角不等式可知，两根红线表示的长度一定小于两根蓝线，故而解开这样的交叉一定能减少路径长度。

我们可以考虑运用迭代思想，在使用nnTSP得到一条可行的路径之后，对路径进行优化，从而得到比nnTSP更优，同时比穷举法更快的解法。

这里我们可以在nnTSP建立了可行路径之后，选用2-Opt作为优化算法，其步骤为：

设路径

长度为

，选择起点

和终点

；

对位于

之间的子路径

进行反转；

如果翻转后的路径较优，将其设置为最优路径

# 2OPT

def opt(cityDist, route, k=2):

# 用2-opt算法优化路径

# 输入：cityDist -- [n,n]矩阵，记录了城市间的距离

# route -- sequence，记录路径

# 输出： 优化后的路径optimizedRoute及其路径长度

nCities = len(route) # 城市数

optimizedRoute = route # 最优路径

minDistance = routeDistance(route) # 最优路径长度

for i in range(1,nCities-2):

for j in range(i+k, nCities):

if j-i == 1:

continue

reversedRoute = route[:i]+route[i:j][::-1]+route[j:]# 翻转后的路径

reversedRouteDist = routeDistance(reversedRoute)

# 如果翻转后路径更优，则更新最优解

if reversedRouteDist < minDistance:

minDistance = reversedRouteDist

optimizedRoute = reversedRoute

return optimizedRoute, minDistance

用2-OPT迭代优化nnTSP:

optimizedRoute, minDistance = opt(cityDist, tour)

print('nnTSP + 2OPT优化后的最优路径为：' + str(optimizedRoute))

print('nnTSP + 2OPT优化后的最优路径长度为：' + str(minDistance))

plotTour(optimizedRoute, cities)

结果：

nnTSP + 2OPT优化后的最优路径为：[0, 15, 7, 14, 2, 16, 19, 26, 5, 4, 9, 28, 3, 22, 29, 1, 11, 13, 10, 12, 24, 18, 23, 8, 17, 25, 27, 20, 21, 6, 0]

nnTSP + 2OPT优化后的最优路径长度为：(4847.935437722887,)

30cities_nnTSP_2OPT

可以看到在应用2-OPT之后，相比优化之前，路径长度减少了308.677；但是另一方面，通过简单的对sub-tour进行翻转并没有能完全uncross所有的交叉。

这个优化方法同样可以用于repNNTSP，来检查是否能找到更优解。用2-OPT迭代优化repNNTSP的代码:

repnnTSPtourOptimized, repnnTSPtourDistOptimized = opt(cityDist, repnnTSPtour, 2)

print('repnnTSPtour + 2OPT优化后的最优路径为：' + str(repnnTSPtourOptimized))

print('repnnTSPtour + 2OPT优化后的最优路径长度为：' + str(repnnTSPtourDistOptimized))

plotTour(repnnTSPtourOptimized, cities)

结果：

repnnTSPtour + 2OPT优化后的最优路径为：[20, 21, 27, 25, 8, 17, 15, 0, 7, 14, 2, 16, 19, 26, 5, 4, 9, 28, 3, 22, 29, 1, 11, 13, 10, 12, 6, 23, 18, 24, 20]

repnnTSPtour + 2OPT优化后的最优路径长度为：(4479.933406593657,)

30cities_repnnTSPtour_2OPT

这个结果看起来已经非常像最优解了：路径之间没有交叉，干净利落的沿着各点形成的convex hull游走。

算法对比

用时

在笔者的工作站上不那么精确地分别测试几种算法在求解规模为30的TSP时的耗时：

在SGA种群规模为100，进化200代时，通常的用时在700ms左右：

707 ms ± 6.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

而nnTSP的用时为：

142 µs ± 5.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

repNNTSP则不出所料是nnTSP的大约30倍：

4.21 ms ± 61.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

而作为迭代改进算法的2-Opt在一条路线上的用时为：

2.64 ms ± 18.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

总结算法的耗时：

问题规模

SGA

nnTSP

repNNTSP

2-OPT

30

707 ms

0.142 ms

4.21 ms

2.64 ms

100

1.76 s

0.727 ms

68.3 ms

2.85 ms

300

4.71 s

3.95 ms

1.19 s

2.87 ms

500

7.63 s

9.6 ms

4.86 s

2.71 ms

1000

35.9 s

34.7 ms

15.2 s

2.64 ms

time_consumption

精度

由于遗传算法并非每次都能稳定收敛于一个值，此处列出的解是5次计算所得的均值

问题规模

SGA

nnTSP

repNNTSP

nnTSP+2-OPT

repNNTSP+2-OPT

30

5156.61234

4529.11558

4847.93543

4479.93341

100

9597.18163

8664.28408

9388.79648

8532.54589

300

14958.54994

14753.59721

14033.02485

13914.14821

为了在速度上稍有竞争力，只有200代迭代的情况下，SGA往往还没有收敛到最佳，波动还是相对较大的。但是规模问题越大，SGA收敛需要的迭代步数越多，强行在200步结束迭代的结果就是获取的解毫无竞争力。

小结

在随机生成的TSP问题中，我们的SGA基本上被其他算法吊打：

SGA需要的时间远大于nnTSP和repNNTSP；而且它对问题的规模更加敏感，在大规模问题上想要收敛到满意的解需要更多迭代步与更大的种群，这又会造成计算时间进一步增加；

可能的改进方向之一是加入环境选择而不是用生成的子代完整替代父代，保证留下精英族群，来加速收敛；

可能的改进方向之二是对种群规模、迭代步数、交叉和突变概率等超参数进行调参，以使算法在TSP问题上更快收敛，并且获得更好的结果；

可能的改进方向之三在于利用问题特性，在SGA的基础上加入局部搜索算法，例如用2-OPT改善精英解；

下一期将会试着按照这三个方向进行SGA的改进。