首页 > 编程札记 > 编程

TSP with Miller

阅读：评论：0

TSP with Miller

高级算法课程作业

1. TSP with MTZ模型及解释

模型解释如下：

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

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

2. 求解TSP with MTZ

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

%%
clear
clc
tic%% 读取城市坐标
file = './TSP_';  % 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:nfor j = i:nif i == jc(i, i) = 9999999;elsec(i, j) = sqrt((data(i, 2)-data(j, 2))^2 + (data(i, 3)-data(j, 3))^2);c(j, i) = c(i, j);endend
end%% 设置约束条件、目标函数
% 等式约束
Aeq = zeros(2*n, num_var);
beq = ones(2*n, 1);% 每个点的出度为1
for i = 1:nAeq(i, (i-1)*n+1:i*n) = 1;   
end% 每个点的入度为1
for i = 1:nAeq(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:nfor j = 2:nif i~= jcnt = cnt + 1;A(cnt, n*n+i) = 1;A(cnt, n*n+j) = -1;A(cnt, (i-1)*n+j) = n;endend
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个实例的求解结果如下。