干货|十分钟教你用动态规划算法解Travelling Salesman Problem（TSP）问题，附代码……

乍一看标题，大家是不是觉得“动态规划”这四个字组合在一起有点眼熟？似乎哪会儿学过来着……但是吧，细细一琢磨，又忘了它具体是什么、怎么用、用来解决哪些问题了。

莫方，小编出现就是为了解决大家一切在学（zhuang）习（bi）上的需求的。动态规划忘了是吧，那今天小编就陪你好好回忆一下。

什么是TSP和动态规划

简单来说，Travelling Salesman Problem (TSP) 是最基本的路线问题。它寻求的是旅行者由起点出发，通过所有给定的需求点后，再次返回起点所花费的最小路径成本，也叫旅行商问题、旅行推销员问题、货郎担问题……

当然，如果你非要把TSP理解成“内容服务提供者”（Telematics Service Provider）小编也不会打你……计算机网络学得不错啊，四级过了吗？

说完TSP问题，咱们再来聊聊什么是动态规划。

动态规划算法（Dynamic Programming，简称DP）通常用于求解具有某种最优性质的问题，其基本思想是将待求解问题分解成若干个子问题，先求解子问题，然后由这些子问题的解再得到原问题的解。

看到这里想必你已经明白了，动态规划恰是一种求解TSP问题的好方法，具体如何求解，我们可以举例实操一下。

实例操作

假设现在有四个城市，它们分别是0，1，2，3，他们之间往来的代价如下图所示：

为了方便起见，我们把它化成二维表的形式：

好了这里要敲黑板划横线了！现在，我们要从城市0出发，期间1，2，3每个城市都必须经过并且只能经过一次，最后回到0，使得路上花费的代价最小。请问你要怎么走？

事实上，这是一个最基本的TSP问题，且构成最优子结构性质，所以可以使用动态规划求解，下面来验证一下此方法求解的可行性。

设 s，s1，s2…s为满足题意的最短回路。假设从s到s1的路径已经确定，则问题转化为从s1到s的最短路径问题。而很显然，s1，s2…s一定可以构成一条最短路径，所以构成最优子结构性质，可以用动态规划求解。

明确问题可解，那下一步就是列方程求解了。

简单推导一下动态规划方程

用 V’ 表示一个点的集合，假设从顶点 s 出发， d ( i , V’ ) 表示当前到达顶点 i，经过 V’ 集合中所有顶点一次的最小花费。

① 当 V’ 为仅包含起点的集合，也就是：

d ( s , { s } ) = 0 ;

② 其他情况，则对子问题求最优解。需在 V’ 这个城市集合中，尝试每一个城市结点，并求出最优解。

③ 最后的求解方式为:

其中 S 为包含所有点的集合。

把公式一套，题就解了。是不是很简单？但是，小编还有更简单的方法。

其实，绝大部分TSP问题都比例子中复杂许多，用程序求解是更好的选择。在这里小编给大家提供一种较为简单的方法，只要把动态规划算法原理掌握好了，代码自然就不难理解了。

用代码前，你需要做哪些准备？

理解状态压缩DP

所谓状态压缩，就是利用二进制以及位运算来实现对于本来应该很大的数组的操作。而求解动态规划问题，很重要的一环就是状态的表示，一般来说，一个数组即可保存状态。但是有这样的一些题目，它们具有DP问题的特性，但是状态中所包含的信息过多，如果要用数组来保存状态的话需要四维以上的数组。于是，我们就需要通过状态压缩来保存状态，而使用状态压缩来保存状态的DP就叫做状态压缩DP。

例题TSP的动态规划方程中，V’ 是一个集合，而对于集合的状态表示最简单的办法就是利用C++中STL里的set，但是这个时候就要考虑一个问题，在代码实现的时候，我们不能用一个集合去做一个数组的下标。自然而然，我们想到可以利用集合的特征值，但这个方法很复杂，而且不容易实现。

小编在这里给大家普及一下位运算的知识。最简单的与（and），或 ( or )，非 ( not )， 大家都很熟悉，和逻辑电路是相通的。而对于异或 ( xor ), 则是很有趣的一种位运算，它的运算规则是相同为 0，不同为 1。例如：

1 xor 1 = 0，0 xor 0 = 0，

1 xor 0 = 1，0 xor 1 = 1；

它的运算满足交换律以及结合律。除此之外，xor 还有很多神奇的操作，有兴趣的同学可以自己去查阅。

再复杂一点的有左移 ( shl )， 右移 ( shr )，相当于对于二进制数的位置移动。例如10001(2) shl 1，就是10001(2)左移一位，变成了100010(2)，换算成十进制，相当于扩大了 2 倍，同理右移则是缩减两倍。那么对于任意的一个二进制数，左移 k 位就是乘 2k, 右移就是整除 2k 。

|||| 小贴士：

在C++中，位运算操作符分别是：

与 &，或 |，非 ~，异或 ^，

左移 <<，右移 >>

推到动态规划方程时，我们注意到 V’ 是一个数的集合，而且解决的问题规模比较小，于是可以用一个二进制数来存储这个集合。简单来说就是——如果城市 k 在集合 V’ 中，那么存储集合的变量 i 的第 k 位就为 1，否则为 0。由于有 n 个城市，所有的状态总数我们用 M 来表示，那么很明显：M = 2^n，而 0 到 2^n -1 的所有整数则构成了 V’ 的所有状态。这样，结合位运算，动归方程的状态表示就很容易了。

准备所需工具

还有两样你需要准备的东西，那就是城市数据文件和编译软件。代码中使用的城市数据文件可以有两种保存格式：一种是上例提到的矩阵式，也可以是 “城市名 城市X坐标 城市Y坐标” 式。大家可以根据实际情况自行调整。

至于编译软件，小编在这里给大家提供的是C++代码，用你用得最顺手的编译器就可以了。小编在这里强烈推荐DEV-CPP！体积小，编译方便，代码还很美观。

好了不啰嗦了，上代码~

代码示例（C++）

#include<bits/stdc++.h>
using namespace std;
// 定义常量
const int INF = 0x3f3f3f3f;
#define sqr(x) ((x)*(x))
// 定义变量
string file_name;
int type; // type == 1 满秩矩阵格式, type == 2 二维坐标式
int s;
int N;// 城市结点数量
int init_point;
double **dp; // 动态规划状态数组dp[i][j]，i表示集合V’，j表示当前到达的城市结点
double **dis; // 两个城市结点之间的距离
double ans;
// 定义结构体
struct vertex{
  double x, y; // 城市结点的坐标
  int id; // 城市结点的id
  int input(FILE *fp){
    return fscanf(fp, "%d %lf %lf", &id, &x, &y);
  }
}*node;
double EUC_2D(const vertex &a, const vertex &b){
  return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
void io(){ // 数据读入
  printf("input file_name and data type\n");
  cin >> file_name >> type;
  FILE *fp = fopen(file_name.c_str(), "r");
  fscanf(fp, "%d", &N);
  node = new vertex[N + 5];
  dis = new double*[N + 5];
  if (type == 1){
    for (int i = 0; i < N; i ++){
      dis[i] = new double[N];
      for (int j = 0; j < N; j ++)
      fscanf(fp, "%lf", &dis[i][j]);
    }
  }
  else{
    for (int i = 0; i < N; i ++)
    node[i].input(fp);
    for (int i = 0; i < N; i ++){
      dis[i] = new double[N];
      for (int j = 0; j < N; j ++)
      dis[i][j] = EUC_2D(node[i], node[j]);// 计算城市之间的距离
    }
  }
  fclose(fp);
  return;
}
void init(){ // 数据初始化
  dp = new double*[(1 << N) + 5];
  for(int i = 0; i < (1 << N); i++){
    dp[i] = new double[N + 5];
    for(int j = 0; j < N; j++)
    dp[i][j] = INF; 
  } // 初始化，除了dp[1][0]，其余值都为INF
  ans = INF;
  return;
}
double slove(){
  int M = (1 << N); 
  // M就是第四部分所说的V’状态总数，1<<N表示2^N，总共有2^N种状态
  dp[1][0] = 0; 
  // 假设固定出发点为0，从0出发回到0的花费为0。TSP只要求是一个环路，所以出发点可以任选
  for (int i = 1; i < M; i ++){ 
  // 枚举V’的所有状态
    for (int j = 1; j < N; j ++){ 
    // 选择下一个加入集合的城市
      if (i & (1 << j)) continue; 
      // 城市已经存在于V’之中
      if (!(i & 1)) continue; 
      // 出发城市固定为0号城市
      for (int k = 0; k < N; k ++){ 
      // 在V’这个城市集合中尝试每一个结点，并求出最优解
        if (i & (1 << k)){ 
        // 确保k已经在集合之中并且是上一步转移过来的结点
           dp[(1 << j) | i][j] = min(dp[(1 << j) | i][j], dp[i][k] + dis[k][j]); // 转移方程
          } // 将j点加入到i集合中
      }
    }
  }
  for (int i = 0; i < N; i ++)
  ans = min(dp[M - 1][i] + dis[i][0], ans);
  // 因为固定了出发点，所以要加上到城市0的距离。另外要从所有的完成整个环路的集合V’中选择，完成最后的转移
  return ans;
}
int main(){
  io();
  init();
  string tmp = file_name + ".sol";
  FILE *fp = fopen(tmp.c_str(), "w");
  fprintf(fp, "%.2lf\n", slove());
  delete[] dp;
  delete[] node;
  delete[] dis;
  fclose(fp);
  return 0;
}

算例运行

例1：满秩矩阵式（type==1）

就拿前文的例子，其文件存储的格式应如下：

4

0 3 6 7

5 0 2 3

6 4 0 2

3 7 5 0

运行结果

11

例2：二维坐标式（type==2）

若城市数据文件如下所示：

16 1 38.24 20.42 2 39.57 26.15 3 40.56 25.32 4 36.26 23.12 5 33.48 10.54 6 37.56 12.19 7 38.42 13.11 8 37.52 20.44 9 41.23 9.10 10 41.17 13.05 11 36.08 -5.21 12 38.47 15.13 13 38.15 15.35 14 37.51 15.17 15 35.49 14.32 16 39.36 19.56

运行结果

73.99

总结

动态规划通过迭代方式寻找每一个子问题的最优解法，因此该解法可以得出TSP的最优解。但算法时间效率较差，因此在问题规模逐渐变大的过程中计算量会急剧膨胀。所以，本算法只适用于小规模求精确解的TSP问题，但对于你平常遇到的大多数TSP问题，这也足够了。

所以，你学会了吗？

-end-