满足三角不等式的TSP问题的近似算法

最近因为工作原因，需要用C++实现图论中的部分算法，突然想起在大学时期，已经将所需的算法用C++实现过了，于是开始在许久没有翻牌子的“大学材料”文件夹中翻找起来。无意中，我找到了当年让我头疼不已的《算法设计与分析》这门课中，我唯一实现的一个算法——满足三角不等式的TSP问题的近似算法。

1. 什么是TSP问题

TSP，全称Traveling Salesman Problem，译为“旅行商问题”。那这个TSP问题究竟是什么呢？

TSP问题的文字描述如下：

一个旅行商人必须访问n个城市，这个旅行商人想要进行一次巡回旅行，或经过哈密顿回路，恰好访问每个城市一次，并最终回到出发城市。旅行所需的全部费用是他旅行中从一个城市到另一个城市的费用的总和，作为商人，他希望整个旅行的费用最低，求此次旅行的路径规划。

我们可以将TSP问题中的几个要素抽象一下：将n个城市抽象为一个具有n个顶点的完全图；将从一个城市到另一个城市的费用抽象为完全图中，边的权值。那么，抽象后的TSP问题的描述如下：

在一个具有n个顶点的无向赋权完全图中，找到一条权值总和最小的哈密顿回路。

注意：这里我们所说的TSP问题特指对称TSP问题，两座城市之间来回的距离是相等的，形成的是一个无向图。而在非对称TSP问题中，可能不是双向的路径都存在，或是来回的距离不同，形成的是有向图。

与TSP问题对应的判定问题的形式语言为：

T S P = { < G , c , k > : G = ( V , E ) 是 一 个 完 全 图 , c 是 V × V → Z 上 的 一 个 函 数 , k ∈ Z , G 中 包 含 一 个 最 大 花 费 为 k 的 旅 行 回 路 } TSP=\{<G, c, k>: G = (V, E)是一个完全图, c是V × V \to Z上的一个函数, k \in Z, G中包含一个最大花费为k的旅行回路\} TSP={<G,c,k>:G=(V,E)是一个完全图,c是V×V→Z上的一个函数,k∈Z,G中包含一个最大花费为k的旅行回路}

经常练习算法的小伙伴们可能知道，TSP问题是图论中最著名的问题之一，而这个问题已经被证明是NP完全的，即该问题找不到一个多项式时间算法来计算最优解。啥是NP完全的？？？

2. NP完全问题

在这里我们将要涉及三类问题：P类问题、NP类问题和NPC类问题（也就是NP完全问题）。

P类问题就是在多项式时间内可以解决的问题。或者说P类问题可以在时间$O(n^k)$内解决，其中$k$为某一常量，$n$是此问题输入的规模。

NP类问题是指那些在多项式时间内可以被证明的问题。“可被证明”是啥意思？假如，当前我们已知一个问题解的证书（certificate），那么可以证明在多项式时间内，在该输入规模下此问题可被解决。所有的P类问题也都是NP类问题，因为如果一个问题是P类问题，那么不用任何证书（certificate），就可以在多项式时间内解决此问题。但是NP=P?这是一个开放性的问题。

NPC类问题通俗的讲，就是如果一个NP问题和其他任何NP问题一样“不易解决”，那么就认为这一问题是NPC类问题或NP完全问题。

说了这么多，其实就是想告诉大家，TSP问题是一个不容易解决的问题，那么作为一名工程师，我们当然不能“在一棵树上吊死”，我们可以开发一种近似算法，来解决NPC问题。

3. 满足三角不等式的TSP问题

回过头来，我们再看一下文章开头我所提到的，满足三角不等式的TSP问题的近似算法。

在TSP问题中，输入是一个赋权无向完全图$G=(V,E)$，其中每条边$(u,v)\in E$都有一个非负的整数代价$c(u,v)$，我们希望找出$G$的一条具有最小代价的哈密顿回路。那么什么是满足三角不等式的TSP问题呢？

在很多实际情况中，从一个地方$u$到另一个地方$w$花费的代价总是最小的。如果一条路径经过了某个中转站，则它不可能具有比直接到达更小的代价。将这种情况形式化，即如果对所有的顶点 $u,v,w\in V$，有：

$c(u,w)⩽c(u,v)+c(v,w)$

就称代价函数$c$满足三角不等式。

不过即使代价函数满足三角不等式，也不能改变旅行商问题的NP完全性。因此不能寄希望于找到一个准确解决旅行商问题的多项式时间算法，相反，我们应该寻找一些好的近似算法。

4. 近似算法

对于满足三角不等式的TSP问题，其近似算法也有很多种，例如：最近邻（nearest neighbour, NN）算法、最小生成树（MST）算法、最小生成树-最小对集（MM）算法。这里我用到的是以最小生成树（MST）为基础的算法：

若实例有$m$个城市，则用$O(m^2)$的时间可找出相应赋权完全图的最小生成树。又因为从该实例得出的旅游回路中，删除其中任一条边也可得到一棵生成树，这棵生成树所有边的长度之和总是不大于这条旅游回路的长度。所以，与TSP问题实例相对应的赋权完全图的最小生成树的边长度（权）之和，必不大于该实例的最短旅游回路的长度。那么怎样通过最小生成树，得到最短旅游回路呢？可以考虑将最小生成树变换为欧拉图，再从欧拉图的欧拉回路过渡到最短旅游回路。

上述提到的欧拉图，是指图中的每个顶点的度均为偶数的无向连通图。而欧拉图中的欧拉回路是指通过图的每条边一次且仅一次的回路。

满足三角不等式的TSP问题的最小生成树算法描述如下：

TSP问题的最小生成树（MST）算法
1    对TSP问题实例的赋权完全图G=(V, E)，调用最小生成树算法，得到G的最小生成树T；
2    复制最小生成树T的每条边，形成欧拉图G'；
3    在G'中寻找欧拉回路；
4    采用“抄近路”方法，将欧拉回路变成TSP旅游回路。
1
2
3
4
5

在算法第4步中，所谓“抄近路”，即在欧拉回路的顶点序列中，删除重复顶点，得到每个顶点出现一次且仅一次的序列，及最短旅游回路。这样得到的回路，从一个顶点到达下一个顶点时，不再经过欧拉回路中的重复顶点，即为抄近路。

下面举一个例子来演示一下MST算法的执行过程：

假设有TSP问题的实例图$G$（线太多，懒得画了），如下图所示：

该实例求得的最小生成树如图(b)所示：

将最小生成树的每条边复制一次后，转换为的欧拉图如下图所示，其欧拉回路为：1,3,2,3,4,5,4,6,7,6,8,6,4,3,1

抄近路后得到的最短旅游回路如下图所示，最短旅游回路为：1,3,2,4,5,6,7,8,1

5. 算法分析

采用MST算法解答实例$I$得出的旅游回路的长度不会超过最小生成树所有变得长度（权）之和的2倍，因此也小于该实例最短旅游回路长度的2倍。证明过程如下：

我们假设TSP问题实例为无向完全图$G$的距离$d(●,●)$。由$G$求得的最小生成树为$T$，则$d(T)<\mathrm{OPT}(I)$，其中$d(T)$为树$T$所有边的长度之和。这是因为在一个最短旅游回路中去掉一条边恰为一棵生成树$T_1$，而$d(T)⩽d(T_1)$。所以欧拉回路$E$的总长度$d(E)<2\mathrm{OPT}(I)$，由抄近路得到的旅游回路长度不超过$d(E)$。故$MST(I)<2OPT(I)$。

6. 算法实现

因为这个算法是大学课程中的一个课程设计，要求是用C++来实现算法，后面抽时间用Python实现一下。

6.1 图的邻接矩阵

这里就用最常用的邻接矩阵作为图的存储结构。

typedef string VerTexType;
typedef int ArcType;

typedef struct      //图的邻接矩阵
{
    VerTexType vexs[MVNum];
    ArcType arcs[MVNum][MVNum];
    int vexnum, arcnum;
}AMGraph;
1
2
3
4
5
6
7
8
9

6.2 无向完全图创建

这里为了节省时间，就没有让用户自己一个一个输入点和边，而是用户只需输入总共有多少个点，然后直接创建一个无向完全图，边的权值则为随机值。

void CreateUDN(AMGraph& G)      //创建无向网
{
    srand((unsigned)time(NULL));
    int num;
    cout<<"please input the quantity of city:";
    cin>>num;
    cout<<endl;
    G.vexnum=num;
    G.arcnum=num*(num-1)/2;
    for(int i=0;i<G.vexnum;i++){
        stringstream ss;
        string b = "city";
        ss << b << (i+1);
        G.vexs[i]= ss.str();
    }
    for(int i=0;i<G.vexnum;i++)
        for(int j=0;j<G.vexnum;j++)
            G.arcs[i][j]=MaxInt;
    for(int i=0;i<G.vexnum;i++)
        for(int j=0;j<G.vexnum;j++)
        {
            if(i<j)
                G.arcs[j][i]=G.arcs[i][j]=rand()%100+1;
        }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

6.3 求解最小生成树并转为欧拉图

在这里我将计算最小生成树以及将最小生成树转为欧拉图放到了一起，别问为什么，问就是实现起来更简单。

AMGraph MST_to_EulerGraph(AMGraph& G, VerTexType u)   //生成最小生成树转为欧拉图
{
    AMGraph EulerG;
    int i, j;
    int m, n;
    int k = LocateVex(G, u);
    closedge c;
    EulerG.vexnum = G.vexnum;
    EulerG.arcnum = 0;
    for (m = 0; m < EulerG.vexnum; m++)
    {
        for (n = 0; n < EulerG.vexnum; n++)
            EulerG.arcs[m][n] = 0;
        EulerG.vexs[m] = G.vexs[m];
    }
    for (i = 0; i < MVNum; i++)
	    c[i].lowcost = MaxInt;
    for (j = 0; j < G.vexnum; j++)
    {
        if (j != k)
	    {
            c[j].adjvex = u;
            c[j].lowcost = G.arcs[k][j];
        }
    }
    c[k].lowcost = 0;
    for (i = 1; i < G.vexnum; i++)
    {
        k = Min(c);
        VerTexType u0 = c[k].adjvex;
        VerTexType v0 = G.vexs[k];
        m = LocateVex(G, u0);
        n = LocateVex(G, v0);
        EulerG.arcs[m][n] = EulerG.arcs[n][m] = 1;
        EulerG.arcnum++;
        c[k].lowcost = 0;
        for (j = 0; j < G.vexnum; j++)
        {
            if (G.arcs[k][j] < c[j].lowcost)
            {
                c[j].adjvex = G.vexs[k];
                c[j].lowcost = G.arcs[k][j];
            }
        }
    }
    return EulerG;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

6.4 得到最短旅游回路并输出

这里我将寻找欧拉回路以及“抄近路”两个步骤放在了一个方法中实现，为什么不分开实现呢，我也向我自己问过这个问题，不知道为什么当年的自己如此的“愚蠢”，没有将步骤分开就算了，更过分的是还没有写注释。唉，先把代码放上来吧。

void Get_TSP_circuit(AMGraph& EulerG)
{
    VerTexType Eul_cir[4 * EulerG.arcnum];
    VerTexType TSP_cir[MVNum + 1];
    VerTexType Stack[EulerG.vexnum];
    int top = -1;
    int i = 0;
    int n = 0;
    int x = 0;
    for (int m = 0; m < EulerG.vexnum; m++)
        visited[m] = false;
    Stack[++top] = EulerG.vexs[0];
    Eul_cir[n] = EulerG.vexs[0];
    visited[0] = true;
    while (top != -1)
    {
        bool flag = true;
        i = LocateVex(EulerG, Stack[top]);
        for (int j = 0; j < EulerG.vexnum; j++)
        {
            if (EulerG.arcs[i][j] == 1 && !visited[j])
            {
                Stack[++top] = EulerG.vexs[j];
                Eul_cir[++n] = EulerG.vexs[j];
                visited[j] = true;
                flag = false;
                break;
            }
        }
        if (flag)
        {
            top--;
            if (top != -1)
                Eul_cir[++n] = Stack[top];
        }
    }

    for (int m = 0; m < EulerG.vexnum; m++)
        visited[m] = false;
    for (int k = 0; k <= n; k++)
    {
        i = LocateVex(EulerG, Eul_cir[k]);
        if (!visited[i])
        {
            TSP_cir[x++] = EulerG.vexs[i];
            visited[i] = true;
        }
    }
    TSP_cir[x] = TSP_cir[0];
    cout << endl << endl;
    cout << "Display TSP_circuit:" << endl << endl;
    for (int k = 0; k <= x; k++)
    {
        cout << TSP_cir[k] << " ";
        if(k < x)
            cout << "-> ";
    }
    cout << endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

6.5 程序运行效果展示

7. 总结

工作了快两年了，回看起初学编程时编写的代码，满满的“青涩感”啊。我也借着这个机会，将放下已久的算法重新拾了起来。

项目的源代码我也已经上传至我的Github上了，有需要的同学自行下载：Traveling-Salesman-Problem 。

如果我前面写的有什么错误，真心的希望大家能够指正，或者小伙伴们在看完以后（前提是能够看懂我的写的是啥o(╥﹏╥)o），认为这个算法还有什么可以优化的地方，都可以给我评论留言哦！