旅行商问题（TSP）状压DP Python代码

本文主要是介绍旅行商问题（TSP）状压DP Python代码，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

来自Wikipedia的定义

The travelling salesman problem (also called the travelling salesperson problem or TSP) asks the following question: "Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?"

 最近为了学习Kool大佬的论文，用神经网络指导动态规划解决TSP，于是乎打算自己做一遍这个思路。反正先不管神经网络用在哪，先把dp写好，然后美丽的马正气又开始全网搜罗代码整理整理就用了呢。。。为了节省我的时间，下面只解释几个重点部分，dp的基本流程可以参考Reference中的链接。

我用0, 1, ..., N - 1来给N个城市编号，比较绕的是以下这一组对应关系：

 城市编号从1到N当然也是可以的吧，但是像正气姐姐这么懒惰的当然是没去研究啦。就按0到N-1这么编号，解释一下这几组对应关系

• 前两行是十进制和对应的二进制位之间的关系，比如二进制数1对应的是十进制数，二进制数10对应十进制数，二进制数11对应十进制数。
• 它们与移位运算的关系是二进制数1对应十进制数，使用移位运算1<<0可以得到相同的数值，十进制数与移位的对应关系是：幂次即需要移位的个数，比如1<<(k-1)即十进制。
• 以上三行对应关系与城市编号无关，而最后一行的城市编号是自己规定的，也可以定义成其他的，我这里用二进制的最末位的0或1表示城市1是否存在在集合中，而城市0永远不会出现在集合中，所以这里有1到N-1共N-1个城市需要每个城市用二进制的一位进行标记。

 上图是从Reference链接中复制过来的，表格画出了dp[i][j]的情况，其中“索引”一行是j的十进制数值，下面一行集合即该数值所表示的集合，比如j为0时其二进制中每一位均为0，即所有城市均不在其所表示的集合中，所以为空集。j为1时其二进制只有最末为为1，即集合中只有城市1，而j为3时其二进制为11，所以城市1和城市2在集合中，j为4时二进制为100，此时只有城市3在集合中。

  依此类推，j所表示的集合最大时即城市1到N-1均在集合内，此时j的二进制中N-1位均为1，对应的十进制数值为。从对应关系的图中可以看到城市N-1对应的位权是，也即其十进制数值，而在对应关系中特别列出了一项，其意义在于我们计算时经常用到的一个技巧：

•  =   对应的二进制即是111 = 1000 - 1这种关系。当二进制所有位均为1时将一位一位按照位权加起来比较麻烦，通常的做法先计算比它大1的数的十进制数值再减去1。再根据十进制与移位的对应关系，这个数值我们可以用(1<<(N-1)) - 1来计算。

 状压dp需要思路在二进制与十进制之间来回切换，其中又涉及到二进制对应的城市编号及用移位运算计算十进制数值，应该先把它们之间的关系理清楚再继续。

 它们之间的对应关系确定之后我们再看其中的核心部分：dp[i][j]这个数组的计算。

dp[i][j]：从城市i出发，经过j表示的集合中的各个城市仅一次然后回到城市0的最短路径长度

  这里我们默认从城市0出发然后回到城市0，所以有了以上dp[i][j]的定义。其中dp[0][]的含义是从城市0出发，经过城市1到N-1然后回到城市0的最短路径长度，即为所求，而当城市i包含在j表示的集合中时，dp[i][j]无意义，因为TSP规定每个城市仅访问一次，从i出发再经过i是访问了两次所以无意义。在写代码时我们跳过数组中这样的位置。dp[i][0]表示从城市i出发，经过空集到城市0的最短路径长度，即边i->0的长度。在dp开始前我们先把dp[i][0]的数值填进去，后续dp[i][j]需要根据它们计算。

  下面3个for循环即核心部分，前面两侧遍历i和j，dp[i][j]应该按列来填，因为后面的列需要根据前面列来计算，所以在里层遍历i，i负责遍历0到N-1共N个城市，而j负责遍历以j表示的集合，该集合空集在初始化时已经填好，这里只需要遍历1到j的最大值。j的最大值应为(1<<(N-1)) - 1，我们将变量M定义为(1<<(N-1))，因为range是左闭右开的，M比需要遍历的最大值大1比较方便。另外，在这两层循环中，我们跳过i包含在j表示的集合中的情况，可以用(j >> (i - 1)) & 1 == 1来表示，一种特殊情况是i为0，此时i必然不在j表示的集合里，是有意义的。我们跳过那些无意义的dp[i][j]。

  最里层的for循环是用来计算当前dp[i][j]的数值的，运用递归的思想，计算从城市i出发，经过j表示的集合中的各个城市仅一次然后回到城市0的最短路径长度，我们可以把j表示的集合中的每个城市都拿出来试一次，用k遍历j中各个城市，代码上即k从1到N-1遍历一遍，再用一个if判断在不在j里，不在就跳过。对于j中的各个城市k，我们都计算一个从i先到k，再从k出发经过其他所有城市仅一次然后回到城市0的最短长度，这个最短长度可以通过dis[i][k] + dp[k][j ^ (1 << (k - 1))]来表示，这个式子里前一半dis[i][k]表示先从i走到k，后一半表示从k出发经过j中除k外其他城市仅一次然后回到城市0所用的最短路径。我们将j中所有的城市k都计算一遍，其中最短的那个就是此时从i出发的最优选择，也即dp[i][j]的数值。这里j ^ (1 << (k - 1)) 是从集合j中删去城市k，将j中城市k对应的二进制位置0即可，删去城市k后j的数值一定比删去之前的小，因为我们j是从小到大去遍历的，所以dp[i][j]中j更小的位置一定被计算过了，此时就可以直接拿来用，这就是为什么dp[i][j]的表格一定要一列一列去填。

  最后我们根据dp[i][j]去反向确定路径，计算得到dp[0][]即为最短的路径长度，此时已知的路径可以表示为 ? ? ... ? ? ? 0。旅行商从城市0出发访问过所有城市后又回到城市0，我们倒序确定路径，先将最后回到的城市0作为后继城市（suc），查找倒数第二个城市，即最后一个0前面的?是1->N-1中的哪个城市。该城市就是计算dp[0][]时最小的城市k，我们重新把dis[i][k] + dp[k][j ^ (1 << (k - 1))]都计算一遍，这里i是后继城市(suc)即城市0，而j从开始。确定k之后我们就确定了 ? ? ... ? ? k 0。此时将找到的城市k设为新的suc，并将k从j表示的集合中删除。反向查找知道所有位置均被确定即为最优路径。

import numpy as np
from scipy.spatial import distance_matrix
import time

N = 20
train_graphs = np.genfromtxt('./data/tsp{}_train.txt'.format(N)).reshape(-1, N, 2)

for graph in train_graphs:
    total_st = dis_st = time.time()
    dis = distance_matrix(graph, graph)
    dis_ed = time.time()
    
    
    # DP
    dp_st = time.time()
    M = 1 << (N - 1) # j的二进制每一位代表该编号的城市存在与否，
                     # 集合j最大为，除城市0外，城市1至N-1均在j表示的集合中,
                     # 此时j的值为2^(N-1)+...+2^1 = 2^N - 1 即 1 << (N - 1) - 1
    
    # dp[i][j]代表从城市i出发，经过j表示的集合中的所有城市仅一次然后回到城市0所用的最短距离
    dp = np.ones((N, M)) * np.inf  
    dp[:,0] = dis[:,0]  # dp[i][0]为从城市i出发，经过空集然后回到城市0所用的最短距离，即边i -> 0的长度
    
    
    for j in range(1, M):  # 左闭右开，j遍历[1, M)即[1，M - 1]
        for i in range(N):
            if i != 0 and (j >> (i - 1)) & 1 == 1:  # j表示的集合中包含城市i则dp[i][j]无意义，跳过
                continue
            for k in range(1, N):
                if (j >> (k - 1)) & 1 != 1: # 遍历在j表示的集合中的各个城市k，不在j中的城市跳过
                    continue
                if dp[i][j] > dis[i][k] + dp[k][j ^ (1 << (k - 1))]: #从i出发先到城市k，再到j中除k外的各个城市
                    dp[i][j] = dis[i][k] + dp[k][j ^ (1 << (k - 1))]      
    dp_ed = time.time()                
    
    
    # 倒序输出路径
    tour_st = time.time()
    suc = 0  # 在计算dp[0][M-1]过程中当前城市的后一个城市编号
    j = M - 1 # j表示的集合从包含除城市0外其他所有城市的集合开始，逐步删除城市直到空集
    tour = [0]
    for i in range(N - 1): # 除城市0外其他城市均需比较长度得到，共需N - 1次
        min_len = np.inf
        for k in range(1, N):
            if (j >> (k - 1)) & 1 != 1: # 遍历在j表示的集合中的所有城市，不在集合内的跳过
                continue
            if min_len > dis[suc][k] + dp[k][j ^ (1 << (k - 1))]: # 找到使dp[suc][j]最小的k，即为当前寻找的城市
                min_len = dis[suc][k] + dp[k][j ^ (1 << (k - 1))]
                tmp = k
        suc = tmp
        tour.append(suc)
        j = j ^ (1 << (suc - 1))
    total_ed = tour_ed = time.time()
    
    
    print("tour: {} tour len: {}".format(tour, dp[0][M - 1]))
    print("Total duration: {}s \nDistances duration: {}s \nDP duration: {}s \nTour duration: {}s \n" \
          .format(total_ed - total_st, dis_ed - dis_st, dp_ed - dp_st, tour_ed - tour_st))
    #break

给一组测试样例，是随机生成的坐标然后用Gurobi求的最优解

0.976064721802 0.245184117574 0.838277998552 0.949267213072 0.0676232302886 0.86665492153 0.954081065089 0.0962588315714 0.634296398269 0.404612925754 0.317284691223 0.540314315697 0.925689168873 0.0987753221956 0.761323280947 0.246835271084 0.374806546136 0.0575070527685 0.617327136243 0.678832372194 0.597783701296 0.489010648008 0.032855462271 0.883801127308 0.268434390934 0.675948261943 0.58618911836 0.755495660388 0.0376380290029 0.596520826603 0.834560198828 0.207960071941 0.879295374647 0.384123896197 0.300851987203 0.384237016803 0.643051333936 0.322582832126 0.000905072843543 0.685928559335 

0 3 6 15 7 8 17 5 14 19 11 2 12 13 1 9 10 4 18 16

 最后，贴一个4个城市的dp数组全貌给大家看看，里面inf的位置是无意义的dp[i][j]，到了最后一列出了城市0，其他城市都在j表示的集合里了，所以最后一列只有一个位置有值哦。

[[0.         1.43487726 2.20135586 2.5931866  0.30107821 1.72881085 2.42566057 2.81749131]
 [0.71743863        inf 1.87574797        inf 1.01137222        inf 2.10005268        inf]
 [1.10067793 1.49250867        inf        inf 1.32498264 1.78644226        inf        inf]
 [0.1505391  1.57827174 2.27512147 2.66695221        inf        inf        inf        inf]]

Reference

旅行商问题（动态规划方法，超级详细的）_仁者乐山智者乐水的博客-CSDN博客_旅行商问kais

这篇关于旅行商问题（TSP）状压DP Python代码的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！