【动态规划】二进制思维（状压DP）巧妙求解TSP（旅行商问题），输出最短路径及其长度（C++）

思路参考：https://blog.csdn.net/joekwok/article/details/4749713

以旅行商问题为关键词搜索，找了一圈没有发现比较优雅的解决方法，大多是固定大小与数据，所以打算自己写一个。

这种算法思想叫做状态压缩动态回归（状压DP），感兴趣的话还可以以此为关键词搜索了解更多有关内容。本文只介绍该算法思想在旅行商问题（Travelling salesman problem，TSP）上的应用。

旅行商问题的讲解可以参考这个链接，我对其代码进行了改进，使得能针对不同的带权图，修改宏表示的城市个数N、全局变量表示的带权图距离矩阵distances[N][N]就能解决不同的TSP（旅行商问题）。本代码中distances[N][N]采用随机数生成，可根据自己需要更改。

由于内存的限制，本代码在城市个数N $\le$ 15时运行正常；当城市个数N$>$ 15时，建议将数据存入文件，通过读取文件实现。

代码的关键是用二进制01来代表城市集合，例如对于城市1，2，3组成的集合 { 1 , 2 , 3 } \{1,2,3\} {1,2,3}，用二进制$111$表示，即十进制的$7$，二维数组dp的第7列；对于城市1，3组成的集合 { 1 , 3 } \{1,3\} {1,3}，用二进制$101$表示，即十进制的$5$，二维数组dp的第5列。

所以，我们的目标是$dp[0][111_2]$，即$dp[0][7_{10}]$，代表从$0$号城市出发，将要走过城市群 { 1 , 2 , 3 } \{1,2,3\} {1,2,3}，动态规划的过程参考上面提到的链接。

下面代码的理解上，重点是看懂里面很多的对二进制数的操作，其他就不难理解了。

填表法的动态规划结果打印都是递归的方法，path[i][j]存储：从i号城市出发，将要经过j的二进制形式表示的城市群，i号城市的下一个城市的编号。

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <iomanip>
using namespace std;

#define N 15
#define MAX 0x3f3f3f3f

// int distances[N][N] = {{MAX, 3, 6, 7}, {5, MAX, 2, 3}, {6, 4, MAX, 2}, {3, 7, 5, MAX}};

int removeCity(int j, int k) { // 从j二进制表示的城市集合中去除k号城市（k位设为0）
    return j - (1 << (k-1));
}

int TSP(int dp[][1 << (N-1)], int path[][1 << (N-1)], int distances[][N]) {
    if (N == 1) return 0; // 只有一个结点

    for (int i = 1; i < N; i++) {
        dp[i][0] = distances[i][0]; // 初始化所有城市到0号城市的距离
    }

    for (int j = 1; j < (1 << (N-1)); j++) {
        for (int i = 1; i < N; i++) {
            if ((j >> (i-1)) % 2 == 0) { // i号城市不在j二进制形式表示的城市集合里
                int min = MAX;
                int next_city = i;
                for (int k = 1; k < N; k++) {
                    if ((j >> (k-1)) % 2 != 0) { // k号城市在j二进制形式表示的城市集合里
                        int temp = distances[i][k] + dp[k][removeCity(j, k)]; // 表示把k号城市从j城市集合中去除
                        if (temp < min) {
                            min = temp;
                            next_city = k;
                        }
                    }
                }
                dp[i][j] = min;
                path[i][j] = next_city;
            }
        }
    }
    
    // 填左上角元素
    int min = MAX;
    int next_city = 0;
    int j = (1 << (N-1)) - 1; // 代表除0号城市外的所有城市组成的集合
    for (int k = 1; k < N; k++) {
        int temp = distances[0][k] + dp[k][removeCity(j, k)];
        if (temp < min) {
            min = temp;
            next_city = k;
        }
    }
    dp[0][j] = min;
    path[0][j] = next_city;

    
    /*for (int i = 0; i < N; i++) { // 打印动态规划表
        for (int j = 0; j < 1 << (N-1); j++) {
            if (dp[i][j] == 0)
                cout << "\\" << " ";
            else
                cout << dp[i][j] << " ";
        }
        cout << endl;
    }*/

    return dp[0][j];
}

void printPath(int path[][1 << (N-1)], int i, int j) { // 打印路线 i为出发城市编号 j为剩下城市组成的集合
    if (j != 0) {
        cout << i << " -> ";
        int next_city = path[i][j];
        printPath(path, next_city, removeCity(j, next_city));
    }
    else {
        cout << i << " -> " << 0;
    }
}

int main() {
    int distances[N][N] = {0};
    int path[N][1 << (N-1)] = {0};
    int dp[N][1 << (N-1)] = {0};
    // 1 << (N-1) == pow(2, N-1)
    // 纵轴为0，1，2...表示城市编号
    // 横轴为000，001，010，...第x位为1，代表该集合里面有x号城市

    // 随机初始化distances[N][N]
    srand((unsigned int) time(NULL));
    for (int i = 0; i < N; i++) {
        for (int j = i; j < N; j++) {
            if (i == j) distances[i][j] = MAX;
            else {
                int temp = rand();
                while (temp == 0) {
                    temp = rand();
                }
                distances[i][j] = distances[j][i] = temp;
            }
        }
    }

    // 打印代价矩阵
    cout << "结点数为：[" << N << "]" << endl;
    cout << "代价矩阵为：" << endl;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (distances[i][j] == MAX) cout << setw(5) << "INF" << " ";
            else cout << setw(5) << distances[i][j] << " ";
        }
        cout << endl;
    }

    // 对动态规划过程计时(毫秒 ms)
    clock_t start = clock();
    int d = TSP(dp, path, distances);
    clock_t spent = (double) (clock() - start) / CLOCKS_PER_SEC * 1000;
    cout << "花费时间：" << spent << "ms" << endl;
    cout << "最短距离为：" << d << endl;
    cout << "旅行路线为：";
    printPath(path, 0, (1 << (N-1))-1);

    return 0;
}

代码运行结果：

结点数为：[15]
代价矩阵为：
  INF  3729 13549 16032  5214 16144 28132 31901 18181 27523 31894 28842 30136 32028 13307
 3729   INF 13870 20822  6156 17167 11654 12498  3912 27011 11758  2954 17805 25778  3749
13549 13870   INF 28116  8097 10456 19276  7130 27613 14974 18319 22185 25320 15417 30216
16032 20822 28116   INF 30596  9669  1172 26208 21886 23824 20555 32440 10568 23940 25786
 5214  6156  8097 30596   INF 25671 25415 30145 13972 22340 10217 18764 26840  1542 26170
16144 17167 10456  9669 25671   INF 23768 14337  7018 15538 14276  9598  3189  8273 16100
28132 11654 19276  1172 25415 23768   INF 30561 26938 31215 11396 24667 13566 11725 18858
31901 12498  7130 26208 30145 14337 30561   INF 25576  2459 16751 19612 19143  2594  5038
18181  3912 27613 21886 13972  7018 26938 25576   INF   162 15277 26820 23305  6925 19514
27523 27011 14974 23824 22340 15538 31215  2459   162   INF  1011 24124 24952  7724 24260
31894 11758 18319 20555 10217 14276 11396 16751 15277  1011   INF  6932  1358 23763 30872
28842  2954 22185 32440 18764  9598 24667 19612 26820 24124  6932   INF 20353 22669 11973
30136 17805 25320 10568 26840  3189 13566 19143 23305 24952  1358 20353   INF 18427  6381
32028 25778 15417 23940  1542  8273 11725  2594  6925  7724 23763 22669 18427   INF 16329
13307  3749 30216 25786 26170 16100 18858  5038 19514 24260 30872 11973  6381 16329   INF
花费时间：13ms
最短距离为：79328
旅行路线为：0 -> 1 -> 11 -> 10 -> 9 -> 8 -> 13 -> 6 -> 3 -> 5 -> 12 -> 14 -> 7 -> 2 -> 4 -> 0