TSP with Miller-Tucker-Zemlin (MTZ) model

高级算法课程作业

1. TSP with MTZ模型及解释

模型解释如下：

1. 目标函数z，对应TSP问题哈密顿回路各有向路径的距离之和；
2. 约束1的数学含义是任意节点 (城市) $i$ 的出度为1，在TSP问题中对应着每个城市只有一个出边；
3. 约束2的数学含义是任意节点 (城市) $j$ 的入度为1，在TSP问题中对应着每个城市只有一个入边。约束1和约束2限制了每个城市最多访问一次，且必须被访问一次。
4. 约束3为了消除子环约束。约束3分为两部分：
   $u_i-u_j+n{x}_{ij}\le n-1$：若一条可行的访问路径中经过$<i,j>$，即$x_{ij}=1$，则不等式等价变换为$u_i\le u_j-1$，即$u_i<u_j$。当一条路径存在重复节点（形成环）时，例如$i\to j\to k\to i$，则有$u_i<u_j<u_k<u_i$，出现矛盾$u_i<u_i$。故$u_i-u_j+n{x}_{ij}\le n-1$能限制路径不包含任意环；
   由于TSP问题本身是一个经过所有节点的最大的环，我们的目的是消除子环而不是消除 (最大的) 环，故使用$1<i\ne j\le n$，即剔除第一个节点，只针对剩余节点，使其不能成环 (是子环)，从而实现消除子环 。

MTZ方法消除子环，只增加了n个决策变量和 (n-1) ²个约束，方法巧妙，方便实现且性能高效。

2. 求解TSP with MTZ

使用MATLAB的混合整数线性规划(MILP)求解器intlinprog来求解。
完整代码如下。

%%
clear
clc
tic

%% 读取城市坐标
file = './TSP_data/burma14.txt';  % burma14  bays29  eil51  eil76  eil101
data = importdata(file);
n = data(1);
num_var = n*n + n;  % 总的决策变量数目
data(1) = [];
data = reshape(data, 3, n);
data = data';

%% 计算两两城市间的欧式距离
c = zeros(n, n);
for i = 1:n
    for j = i:n
        if i == j
            c(i, i) = 9999999;
        else
            c(i, j) = sqrt((data(i, 2)-data(j, 2))^2 + (data(i, 3)-data(j, 3))^2);
            c(j, i) = c(i, j);
        end
    end
end

%% 设置约束条件、目标函数
% 等式约束
Aeq = zeros(2*n, num_var);
beq = ones(2*n, 1);

% 每个点的出度为1
for i = 1:n
    Aeq(i, (i-1)*n+1:i*n) = 1;   
end

% 每个点的入度为1
for i = 1:n
    Aeq(i+n, i:n:n*n) = 1;
end

% 不等式约束 MTZ
A = zeros((n-1)*(n-1) - (n-1), num_var);
b = ones((n-1)*(n-1) - (n-1), 1) .* (n-1);
cnt = 0;
for i = 2:n
    for j = 2:n
        if i~= j
            cnt = cnt + 1;
            A(cnt, n*n+i) = 1;
            A(cnt, n*n+j) = -1;
            A(cnt, (i-1)*n+j) = n;
        end
    end
end

% 01整数约束
intcon = 1:n*n;
bound_lower = zeros(num_var, 1);
bound_lower(n*n+1:num_var) = -inf;
bound_upper = ones(num_var, 1);
bound_upper(n*n+1:num_var) = inf;

% 目标函数
f = zeros(num_var, 1);
f(1:n*n) = reshape(c, n*n, 1);

%% 求解器求解
options = optimoptions('intlinprog','AbsoluteGapTolerance',1e-3,...
                       'MaxTime', 1000);
[x, y]=intlinprog(f,intcon, A, b, Aeq, beq, bound_lower, bound_upper);
toc

由于部分实例 (后三个) 求解时间较长，对其设置最大求解时间，所有5个实例的求解结果如下。

burma14	bays29	eil51	eil76	eil101
30.8785	9074.1480	436.8716	599.9859	681.2858

所有的运行截图在result文件下，为避免杂乱，这里仅给出burma14的运行信息。

3. 实验总结与心得

1. 通过对经典规划问题TSP的建模与求解，了解了TSP问题的建模方式；
2. 学习并理解了TSP问题建模的核心难点，即如何消除子环，并对MTZ方法的原理进行了剖析；
3. 用求解器求解TSP问题，进一步熟悉了MATLAB混合整数规划求解器的编程。