旅行商问题（TSP）系列(1)

作者: LeePlusPlus | 来源:发表于2020-07-23 19:19 被阅读0次

尝试使用各种优化算法，对旅行商问题（TSP）及多旅行商问题（MTSP）进行解决。主要使用语言为python。

题目

(1) 在地图上给定30个节点的位置，标号为0--29，每两个节点之间都有一条路径相互连通，求一条回路能遍历每一个节点，并且令其总长度尽可能小。

(2) 在第一问的基础上，求四条回路，它们加起来能遍历每一个节点，并且令这四条回路的总长度尽可能小。

(3) 在第二问的基础上，增加限制条件，要求这四条回路都经过0节点，令这四条回路的总长度尽可能小。

附件：40个节点的位置（按顺序分别为0节点，1节点，……29节点）：

(45, 0), (67, 17), (22, 56), (42, 9), (48, 13),
(4, 66), (58, 9), (5, 0), (81, 41), (97, 21),
(63, 61), (27, 76), (46, 73), (41, 63), (19, 83),
(86, 46), (73, 83), (9, 25), (88, 96), (91, 63),
(4, 76), (51, 56), (62, 71), (76, 63), (90, 40),
(86, 15), (25, 63), (15, 54), (18, 37), (61, 10)

30个点的位置可视化

数据预处理

先对点的位置数据进行预处理，并准备几个常用函数。

# basic.py
# 对点的位置数据进行预处理
# 以及后面会用到的常用函数
from math import sqrt
import numpy as np
import matplotlib.pyplot as plt

# 节点数量
num_pos = 30

# 各节点位置
pos = (
    (45, 0), (67, 17), (22, 56), (42, 9), (48, 13), 
    (4, 66), (58, 9), (5, 0), (81, 41), (97, 21), 
    (63, 61), (27, 76), (46, 73), (41, 63), (19, 83), 
    (86, 46), (73, 83), (9, 25), (88, 96), (91, 63), 
    (4, 76), (51, 56), (62, 71), (76, 63), (90, 40), 
    (86, 15), (25, 63), (15, 54), (18, 37), (61, 10)
    )

# 计算两点间距离
def get_dist(pos_1, pos_2):
    return sqrt((pos_1[0] - pos_2[0])**2 + (pos_1[1] - pos_2[1])**2)

# 各节点之间的距离
dist = [[get_dist(pos_1, pos_2) for pos_2 in pos] for pos_1 in pos]


# 绘制散点图
def show_points(path):
    # 例：{v_0, v_2, v_1, v_3}: path = [0, 2, 1, 3]
    x = [pos[v_i][0] for v_i in path]
    y = [pos[v_i][1] for v_i in path]
    # 绘制散点图
    plt.scatter(x, y)
    # 标号
    for i in range(len(path)):
        plt.annotate(path[i], (x[i] + 1, y[i] - 1))
    # plt.show()


# 绘制回路图
def show_path(path):
    # 例：v_0-v_2-v_1-v_3-v_0: path = [0, 2, 1, 3, 0]
    x = [pos[v_i][0] for v_i in path]
    y = [pos[v_i][1] for v_i in path]
    # 绘制散点图
    plt.scatter(x, y)
    # 绘制连线
    plt.plot(x, y)
    # 标号
    for i in range(len(path)):
        plt.annotate(path[i], (x[i] + 1, y[i] - 1))
    # plt.show()

第一问：状压dp

思路

(1) 在地图上给定30个节点的位置，标号为0--29，每两个节点之间都有一条路径相互连通，求一条回路能遍历每一个节点，并且令其总长度尽可能小。

很显然，第一问便是求最短哈密顿回路，也就是很经典的旅行商问题。旅行商问题的可行解是所有顶点的全排列，随着顶点数的增加，会产生组合爆炸，难以用穷举法解出答案。^[1]但考虑到这道题只有30个点，并不算太多（后面被计算量狠狠打脸），而且也不确定是否能接受非最优解，我一开始决定先用搜索法穷举所有可能，找出全局最优解。

在方法选择上，深度优先搜索+剪枝显然不够给力，于是我选择了状压dp（状态压缩动态规划）。^[2]

用离散数学的知识，对这道题进行抽象化描述：

已知无向完全图 $G=<V,E>$ 中， $V=\{v_0, v_1, v_2, ..., v_{29}\}$ 。设状态 $s$ 为一个长度为30的01串，用来表示 $V$ 的一个子集 $V_s$ 。比如用 $s_{11}=000000000000000000000000001101$ ，表示 $V$ 的一个子集 $V_{11} = \{v_0, v_2, v_3\}$ 。通过这种二进制的表示方法，我们可以把30个节点的所有状态（子集）压缩到一个int里，此时 $s$ 的取值范围是 $[0, 2^{30} -1]$ 。

设 $dist(i, j)$ 为 $v_i$ 和 $v_j$ 间的距离。设 $min\_length(s, i)$ 为： $G$ 的完全无向子图 $G_s=<V_s,E_s>$ 中，以 $v_0$ 为起点， $v_i$ 为终点的哈密顿通路（经过 $s$ 中所有节点一次且仅一次的通路）的最短长度（ $v_0,v_i \in s$ ）。

若 $min\_length(s +v_i, i) \ (v_i \notin s)$ 所指的这条哈密顿通路为 $v_0 ... v_jv_i$ ，则 $d(v_0...v_jv_i) = d(v_0...v_j) + d(v_iv_j)$ ，因而我们可以得到状态转移方程：
$min\_length(s + v_i, i) = \min_{j \in s, \ j \neq 0}{\left(min\_length(s, j) + dist(i, j) \right)} \\ (v_0 \in s, v_i \notin s)$
该状态转移方程的含义是：在到达 $v_i$ 之前，我们先到达的是 $v_j$ ；在到达 $s+v_i$ 这个状态前，我们先到达的是 $s$ 这个状态。因此只要得到所有 $min\_length(s, j)$ 的值，我们就能得知 $min\_length(s + v_i, i)$ 的值。换言之，每次到达一个 $min\_length(s, j)$ 时，我们都可以更新所有 $min\_length(s + v_i, i)$ 的值。 $(j \in s, i \notin s)$

我们发现，压缩到int中时， $s + v_i$ 记作s | (1 << i)，换言之 $int(s + v_i) > int(s)$ 恒成立。这启示我们可以以一个有序（从小到大）的顺序遍历 $s$ ，使得每次抵达一个 $s+v_i$ 时，都可以保证它已经被所有的 $s$ 更新过。这样，我们就可以通过动态规划，得到最终 $G$ 中以 $v_0$ 为起点，每一个以 $v_i$ 为终点的哈密顿通路的最短长度。那么，可以通过 $d(v_0..v_i v_0) = d(v_0...v_i) + d(v_i v_0)$ ，得到最短哈密顿回路的长度 $ans\_length$ ：
$ans\_length = \min_{i \in s, \ i \neq 0}{\left( min\_length(s, i) + dist(i, 0) \right)}$

代码

这种算法用c++跑肯定会快一些，但是python比较方便（可以作图）。

# q_1_1.py
# 状压dp求解最短哈密顿回路
from basic import *

# 先用较少的节点数做测试
num_pos = 10
# 无穷大常量
inf = 987654321
# s中以0为起点，i为终点的最短哈密顿通路长度
min_dist = [[inf for i in range(num_pos)] for s in range(1 << num_pos)]
# s中以0为起点，i为终点的最短哈密顿通路
# 例: v_0-v_3-v_1-v_2表示为[3, 1, 2]
min_path = [[[] for i in range(num_pos)] for s in range(1 << num_pos)]

print('init succeed')


# 状压dp搜索
for s in range(1, 1 << num_pos, 2):
    # 考虑节点集合s必须包括节点0
    if not (s & 1):
        continue
    for j in range(1, num_pos):
        # 终点i需在当前考虑节点集合s内
        if not (s & (1 << j)):
            continue
        if s == int((1 << j) | 1):
            # 若考虑节点集合s仅含节点0和节点j，dp边界，赋予初值
            # print('j:', j)
            min_path[s][j] = [j]
            min_dist[s][j] = dist[0][j]

        # 枚举下一个节点i，更新
        for i in range(1, num_pos):
            # 下一个节点i需在考虑节点集合s外
            if s & (1 << i):
                continue
            # 更新min_dist[s + i][i], min_path[s + i][i]
            if min_dist[s][j] + dist[j][i] < min_dist[s | (1 << i)][i]:
                min_path[s | (1 << i)][i] = min_path[s][j] + [i]
                min_dist[s | (1 << i)][i] = min_dist[s][j] + dist[j][i]


ans_dist = inf
ans_path = [] 
# 求最终最短哈密顿回路
for i in range(1, num_pos):
    if min_dist[(1 << num_pos) - 1][i] + dist[i][0] < ans_dist:
        # 更新，回路化
        ans_path = [0] + min_path[s][i] + [0]
        ans_dist = min_dist[(1 << num_pos) - 1][i] + dist[i][0]

# 输出结果
print()
print('min length:', ans_dist)
print('path:', ans_path)

show_path(ans_path)
plt.show()

运行结果

在第5行注意到，我只使用10个节点做测试，结果如下：

1595563051243.png

Figure_2-1595560139309.png

可以看到10节点时，状压dp跑的很快，解的形状看上去也完全没问题。接下来慢慢放开节点数量，看看跑得如何。

节点数	所用时间	最短长度
10	0.8s	282.9
15	2.1s	325.8
18	14.4s	366.0
20	100.3s	401.3
21	144.0s	405.1
22	1013.3s	409.7
...	...	...

……我佛了。

（这个结论或许可以记一下？22个点以内的TSP都可以用状压dp找最优解）

非常明显，这个算法的时间开销是我们不能接受的（想打死之前觉得30不算多的我）。掐指一算，该算法的时间复杂度是 $O(n^2 \cdot 2^n)$ ，空间复杂度是 $O(n \cdot 2^n)$ ，（光是数组初始化也要很长时间）。这个令人绝望的增长速度也打消了我“用c++重跑一遍”的念头。写到这里，我感觉30个点的最优解求不出来，决定放弃状压dp，改用一些比较随机化的算法求一个较优解。

第一问：模拟退火算法^[3]

简介

模拟退火算法是一种启发式搜索算法，即按照预定的控制策略进行搜索，在搜索过程中获取的中间信息将用来改进控制策略。它的思想借鉴于固体的退火原理，当固体的温度很高的时候，内能比较大，固体的内部粒子处于快速无序运动，当温度慢慢降低的过程中，固体的内能减小，粒子的慢慢趋于有序，最终，当固体处于常温时，内能达到最小，此时，粒子最为稳定。模拟退火算法便是基于这样的原理设计而成。

模拟退火算法从某一高温出发，在高温状态下计算初始解，然后以预设的邻域函数产生一个扰动量，从而得到新的状态，即模拟粒子的无序运动，比较新旧状态下的能量，即目标函数的解。如果新状态的能量小于旧状态，则状态发生转化；如果新状态的能量大于旧状态，则以一定的概率准则发生转化。当状态稳定后，便可以看作达到了当前状态的最优解，便可以开始降温，在下一个温度继续迭代，最终达到低温的稳定状态，便得到了模拟退火算法产生的结果。

在这道题中，将回路 $\Gamma = v_{i_0} v_{i_1} v_{i_2} ... v_{i_{29}} v_{i_0}$ 表示为x = [i_0, i_1, i_2, ..., i_29]，将 $x$ 在附近邻域随机扰动设为令x中多个元素随机交换位置，便可以使 $x$ 的取值范围覆盖整个解空间。

代码

于是我便得到一个（比较粗糙的）模拟退火算法的代码。

# q_1_1.py
# 模拟退火算法求解最短哈密顿回路
from basic import *
from random import randint, uniform, shuffle
from math import exp, log

# 计算回路长度
def get_path_length(x):
    length = 0
    for i in range(len(x) - 1):
        length += dist[x[i]][x[i + 1]]
    length += dist[len(x) - 1][0]
    return length

# 给予x一个随机位移
def get_next_x(x):
    x_next = x.copy()
    # 多次交换
    for i in range(int(T / T_0 * 8 + 2)):
        i_1, i_2 = randint(0, num_pos - 1), randint(0, num_pos - 1)
        x_next[i_1], x_next[i_2] = x_next[i_2], x_next[i_1]
    return x_next

# 采纳概率
# 自己瞎捣鼓出来的
def get_accpectable_probability(y, y_next, T):
    return exp(-T_0 / T * (y_next - y) / y)


# # 温度初值
T_0 = 100000
T = T_0
# 温度终止值
T_min = 30
# x: 回路
# 例: v_0-v_1-v_2-v_0: x = [0, 1, 2]
# x初值随机
x = list(range(num_pos))
shuffle(x)
# y: 回路长度
y = get_path_length(x)
# 内循环次数
k = 100
# 计数器
t = 0           

# 存储求解过程
x_list = [x]
y_list = [y]

x_best = x.copy()
y_best = get_path_length(x)

# 开始模拟退火
while T > T_min:
    # 内循环
    for i in range(k):
        y = get_path_length(x)

        # 给予x随机位移
        x_next = get_next_x(x)

        # 试探新解
        y_next = get_path_length(x_next)
        if y_next < y:
            # 更优解，一定接受新解
            x = x_next
            # 更新最优记录
            if y_next < y_best:
                y_best = y_next
                x_best = x_next.copy()
        else:
            # 更劣解，一定概率接受新解
            p = get_accpectable_probability(y, y_next, T)
            r = uniform(0, 1)
            if r < p:
                x = x_next

    # 做记录
    x_list.append(x)
    y_list.append(y_next)

    t += 1
    if t % 1000 == 0:
        print(T)

    # 降温
    T = T_0 / (1 + t) # 快速降温
    # T = T_0 / log(1 + t) # 经典降温，慢的一批，令人暴躁


x_best += [x_best[0]]
# 输出解
print('min length:', y_best)
print('path:', x_best)

plt.figure(1)
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)

# 绘制连线图
plt.sca(ax1)
show_path(x_best)

# 绘制退火图
plt.sca(ax2)
x = np.linspace(1, len(x_list), len(x_list))
y = np.linspace(1, len(x_list), len(x_list))
for i in range(len(x_list)):
    y[i] = y_list[i]
plt.plot(x, y)

plt.title('模拟退火进程', fontproperties = 'SimHei')
plt.xlabel('遍历次数', fontproperties = 'SimHei')
plt.ylabel('回路长度', fontproperties = 'SimHei')

plt.show()

运行结果

哪怕是在30个节点的情况下，速度依然飞快。多次运行，结果如下：

运行次数	回路长度	运行时间
1	514.5	9.9s
2	486.7	11.3s
3	486.8	10.3s
4	497.4	11.5s
5	463.0	10.5s
6	529.2	11.2s
7	458.5	10.4s
8	455.4	10.8s
...	...	...

一些搜索结果如下，可以看出，有时候找到的解明显不太行，有时候找到的解很舒服。

Figure_3_1.png

Figure_3_2.png

1595566034252.png

总的来说，从右图的模拟退火进程可以看出，模拟退火受初值影响较大，最终的结果主要取决于一开始那一小段能降到多低，在后续大部分时间难以进一步优化。因此，每一次跑的结果差异显著，有时很好，有时很烂，全看运气，不能做到“最终都收敛于相近的答案”。

跑了四五十遍，得到的最短回路长度为453.2，结果如下：

Figure_3_4.png

凭感觉，这个路线已经找不出什么明显的优化空间了，将这个路线作为第一问的解答足以令人满意。

待续

旅行商问题（TSP）系列(2)

参考资料

旅行商问题百度百科 https://baike.baidu.com/item/%E6%97%85%E8%A1%8C%E5%95%86%E9%97%AE%E9%A2%98/7737042 ↩
最短哈密顿路径（位运算+dp） solego https://blog.csdn.net/weixin_43900869/article/details/102934460 ↩
模拟退火算法与其python实现 WFRainn https://blog.csdn.net/wfrainn/article/details/80303138 ↩

旅行商问题（TSP）系列(1)

题目

数据预处理

第一问：状压dp

思路

代码

运行结果

第一问：模拟退火算法^[3]

简介

代码

运行结果

参考资料

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

交通运输工程