思想往往在最黑暗中的闪光处诞生。

突然发现解 TSP 问题的新方法,即用线性规划的方法。

首先重述下 TSP 问题:

TSP 问题

给定一个带权图,寻找一条最短路径,从结点 1 出发最后仍回到结点 1,且必须且只能经过所有其他结点一次。

也可以简单地说,寻找最短 Hamilton 圈。

我们可以用线性规划的想法来做。令xij=1x_{ij}=1,如果最后的 Hamilton 圈上有从结点ii 到结点jj 的边,而xij=0x_{ij}=0 表示没有这条边。故我们需要最小化函数:

minxi,jwij×xij\min_{x}\sum_{i,j}w_{ij}\times x_{ij}

其中wijw_{ij} 表示iijj 的边权值。

然后为了每个结点只能进出一次,故需要加约束条件:

{j=1nxij=1,i=1,2,...,ni=1nxij=1,j=1,2,...,nxij=0  or  1\begin{cases} \sum_{j=1}^nx_{ij}=1,i=1,2,...,n\\ \sum_{i=1}^nx_{ij}=1,j=1,2,...,n\\ x_{ij}=0\;or\;1 \end{cases}

看起来很好,然而需要考虑个问题就是这个约束不够强,因为它允许子圈的出现。譬如选择路径12,21,34,431\rightarrow 2,2\rightarrow 1,3\rightarrow 4,4\rightarrow 3 也满足上述约束,但 Hamilton 圈不允许分成两部分走。

那么就需要增加约束,而常用的一个消除子圈的约束为 Miller Tucker Zemlin 约束。它引入了nn 个新的变量ui,i=1,...,nu_i,i=1,...,n,要求满足:

uiuj+nxijn1i=1,2,...,nj=2,3,...,nui0u_i-u_j+nx_{ij}\leq n-1\\ i=1,2,...,n\\ j=2,3,...,n\\ u_i\geq 0

有了这个约束有什么用呢?下面证明,若存在满足这样约束,则一定没有子圈。

反设存在子圈,则考虑不包含结点11 的那个子圈i1i2...imi1i_1\rightarrow i_2\rightarrow...\rightarrow i_m\rightarrow i_1。就会有:

ui1ui21ui2ui31...uimui11u_{i_1}-u_{i_2}\leq -1\\ u_{i_2}-u_{i_3}\leq -1\\ ...\\ u_{i_m}-u_{i_1}\leq -1\\

那么把它们全加起来,就会有0m0\leq -m,矛盾,故错误。

反过来,如果存在一个 Hamilton 圈,则把uiu_i 设为ii 在圈上的位置就好(圈上第几个结点),易验证是满足的。

至于ui0u_i\geq 0 的条件为什么要加,首先加上肯定是没有错的。在实验中我发现不加也可以算出正确结果,但是会导致输出时较麻烦(matlab),而且加了也可以缩小空间。

例子代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
% TSP问题
M=1000000; n=19; w = zeros(n);
w(1,[2,3])=[300,600]; w(2,[3,8])=[300,200];
w(3,[5,9,10])=[250,500,500]; w(4,[5,7])=[200,950];
w(5,[6,10])=250; w(6,[7,11])=[500,250]; w(7,19)=1150;
w(8,[9,13])=[200,750]; w(9,[10,16])=[450,600];
w(10,[11,12])=[250,200]; w(11,12)=500; w(12,16)=250;
w(13,14)=200; w(14,15)=200; w(15,17)=600;
w(16,18)=300; w(17,[18,19])=[200,600]; w(18,19)=600;
w=w+w'; w(w==0)=M; % w是图的邻接矩阵
x = optimvar("x",n,n,"Type","integer","LowerBound",0,"UpperBound",1);
u = optimvar("u",n,"LowerBound",0);
prob = optimproblem("Objective",sum(w.*x,"all"));
prob.Constraints.c1 = [sum(x)'==1; sum(x,2)==1;];
prob.Constraints.c3 = repmat(u,1,n-1)-repmat(u(2:n)',n,1)+n*x(:,2:n)<=n-1;
[sol,fval,flag,out] = solve(prob), xx=round(sol.x)
[i,j]=find(xx);
ind = [i'; j'] %输出树的顶点编号