运筹系列65：使用Julia精确求解tsp问题

1. 问题概述

tsp问题的数学模型如下：

1.1 给定upbound的Christofides方法

这是可以给出上界的一个方法，可以证明构造出的路线不超过最优路线的1.5倍。步骤为：
1）构造MST（最小生成树）
2）将里面的奇点连接起来构成欧拉回路称为完美匹配。Edmonds给出了多项式时间内构造最小代价完美匹配的方法，其长度不超过最优解的1.5倍。

证明方法也很直观，奇点最短路径可以拆分成两条完美匹配，其中总有一条的长度$\le$最优奇点路径长度/2$\le$最优完整路径长度/2

3）按欧拉回路顺序逐个扫描点，跳过重复经过的点，即可构造一条完整的路径。

下面是blossom算法的使用示例：

using BlossomV,TravelingSalesmanHeuristics,TSPLIB,Distances,DataStructures
tsp = readTSP("vlsi/bcl380.tsp"); # bays29
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
mst = TravelingSalesmanHeuristics.minspantree(distmat)[1]
x = counter(cat([m[1] for m in mst],[m[2] for m in mst],dims=1))
odds = []
for xi in x
    if xi[2]%2==1
        append!(odds,xi[1])
    end
end
mat = Matching(Float64, length(odds))
for i in 1:length(odds)
    for j in 1:i-1
        add_edge(mat,i-1,j-1,distmat[odds[i],odds[j]])
    end
end
solve(mat)

using PyPlot
tree = mst
PyPlot.scatter(tsp.nodes[:,1],tsp.nodes[:,2],color="black",s=10)
for i in 1:length(tree)
    l = tsp.nodes[[tree[i][1],tree[i][2]],:]
    PyPlot.plot(l[:,1], l[:,2], color="black",linewidth = 1)
end
for i in 1:length(odds)
    j = get_match(mat,i-1)+1
    if j>i
        l = tsp.nodes[[odds[i],odds[j]],:]
        PyPlot.plot(l[:,1],l[:,2],color="r",linewidth = 0.5,linestyle="--")
    end
end

1.2. 给定lowerbound的松弛问题

首先看一个例子：

我们定义$x_{ij}$为边ij是否在路径上。根据TSP问题的要求，每个点只能路过一次，因此每个点连着两条边，有：

1.3 MIP求解方法

可以参考TravelingSalesmanExact 包：

using TravelingSalesmanExact, GLPK
tsp = readTSPLIB(:att48)
cities = [tsp.nodes[i,:] for i in 1:tsp.dimension];
@time tour, cost = get_optimal_tour(cities, GLPK.Optimizer)

我们自己实现的代码如下。如果我们能够调用MIP求解器，那么我们可以使用combination包把所有可能的子集罗列出来求解。当然也可以松弛subtour约束后使用行生成求解：

using TSPLIB,JuMP, GLPK, Distances
tsp = readTSP("25.tsp")
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
m = Model(GLPK.Optimizer)
@variable(m, x[1:N,1:N], Bin)
@objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:N))
for i=1:N 
    @constraint(m, x[i,i] == 0)
    @constraint(m, sum(x[i,1:N]) == 1)
end
for j=1:N
    @constraint(m, sum(x[1:N,j]) == 1)
end
for f=1:N, t=1:N
    @constraint(m, x[f,t]+x[t,f] <= 1)
end

optimize!(m)

结果如下：

接下来不断增加subtour约束重新求解

function is_tsp_solved(m)
    g = JuMP.value.(x)
    N = size(g,1)
    path = Int[]
    push!(path,1)
    while true
        v, idx = findmax(g[path[end],:])
        if idx == path[1]
            break
        else
            push!(path,idx)
        end
    end
    n = length(path)
    if n < N
        @constraint(m,sum(x[path,path])<=n-1)
        return false
    end
    return true
end

while !is_tsp_solved(m)
    optimize!(m)
end

最终结果如下

整数规划还是比较费时间的，att48用了177s才求出最优解。下面进行优化，每次同时把所有的subtour都添加进来。我们修改后的代码如下：

using JuMP, GLPK, Distances

function get_cycles(perm_matrix)
    N = size(perm_matrix, 1)
    remaining_inds = Set(1:N)
    cycles = Vector{Int}[]
    while length(remaining_inds) > 0
        cycle = find_cycle(perm_matrix, first(remaining_inds))
        push!(cycles, cycle)
        setdiff!(remaining_inds, cycle)
    end
    return cycles
end

function find_cycle(perm_matrix, starting_ind = 1)
    cycle = [starting_ind]
    prev_ind = ind = starting_ind
    while true
        next_ind = findfirst(>(0.5), @views(perm_matrix[ind, 1:prev_ind-1]))
        if isnothing(next_ind)
            next_ind = findfirst(>(0.5), @views(perm_matrix[ind, prev_ind+1:end])) +
                       prev_ind
        end
        next_ind == starting_ind && break
        push!(cycle, next_ind)
        prev_ind, ind = ind, next_ind
    end
    return cycle
end



function is_tsp_solved(m,x)
    cycles = get_cycles(JuMP.value.(x))
    if size(cycles,1)>1
        for cycle in cycles
            @constraint(m,sum(x[cycle,cycle])<=2*size(cycle,1)-2)
        end
        return false
    end
    return true
end

function solve_tsp_with_mip(tsp)
    N = tsp.dimension
    distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
    m = Model(GLPK.Optimizer)
    @variable(m, x[1:N,1:N], Symmetric, Bin)
    @objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:i))
    for i=1:N 
        @constraint(m, x[i,i] == 0)
        @constraint(m, sum(x[i,:]) == 2)
    end
    optimize!(m)
    while !is_tsp_solved(m,x)
        optimize!(m)
    end
    return JuMP.value.(x), objective_value(m)
end
@time solve_tsp_with_mip(tsp)

1.4 BFF建模方法

给路径上的边定义一个数，第一条边为1，第二条边为2，……，并且定义方向，则除了depot外每个点的所有边加起来为1，设计模型如下：

using DataFrames, Plots, LinearAlgebra, Random
using TSPLIB,LKH,Distances,PyPlot
tsp=readTSPLIB(:att48)
c = tsp.dimension
distance = 1000;
xpos = tsp.nodes[:,1]; ypos = tsp.nodes[:,2];
Nodes = DataFrames.DataFrame(;k = 1:c,x = xpos,y = ypos,);
Distance = zeros(Int(c*(c-1)/2),3)
global cont = 1
for i = 1:size(Nodes,1)
    for j = i+1:size(Nodes,1)
        global Distance[cont,:] = 
            [i,j,floor(0.5+norm([(Nodes.x[i,1]-Nodes.x[j,1]),(Nodes.y[i,1]-Nodes.y[j,1])]))]
        global cont = cont+1;
    end
end
L = size(Distance,1)
branches = DataFrames.DataFrame(; k = Distance[:,1],m = Distance[:,2], lkm = Distance[:,3],)
# 添加model求解代码
using JuMP, HiGHS
b = size(branches, 1) # 边的个数
A = zeros(c, b) # 点是否为变的起点或者终点
for l = 1:b
    k = Int(branches.k[l, 1]);
    m = Int(branches.m[l, 1])
    A[k, l] = 1; A[m, l] = -1;
end
Branch = 1:b; Cities = 1:c;
model = Model(HiGHS.Optimizer)
set_silent(model)
sd = ones(c, 1); sd[1, 1] = 0; smax = zeros(c, 1); smax[1, 1] = c; flmax = c;
@variable(model, x[l in Branch], Bin) # 是否包含边
@variable(model, f[l in Branch]) # 边上的流量
@variable(model, s[k in Cities] >= 0)
@objective(model, Min, sum(branches.lkm[l, 1] * x[l] for l in Branch)) # 总距离最短
for k in Cities
    @constraint(model, s[k] <= smax[k])
    @constraint(model, sum(abs(A[k, l]) * x[l] for l in Branch) == 2) # 每个点包含2个边
    @constraint(model, s[k] - sd[k] == sum(-A[k, l] * f[l] for l in Branch)) # depot的流为47（因为有一个回来的48），其他点的流为1
end
for l in Branch
    @constraint(model, -flmax * x[l] <= f[l]) # 只有有边才能有流
    @constraint(model, f[l] <= flmax * x[l])
end
@constraint(model, sum(x[i] for i in Branch) == c) # 总共有c条边
@time JuMP.optimize!(model)
@show objective_value(model)
# 添加可视化代码
global sol = [value(x[l]) for l in Branch];
global sou = [value(s[k]) for k in Cities];
global flow = [value(f[l]) for l in Branch];
for l = 1:size(branches, 1)
    if sol[l, 1] >= 1/2
        local index1 = findall(x -> x == branches.k[l, 1], Nodes.k)
        local index2 = findall(x -> x == branches.m[l, 1], Nodes.k)
        local x = [Nodes.x[index1, 1]; Nodes.x[index2, 1]];
        local y = [Nodes.y[index1, 1]; Nodes.y[index2, 1]];
        plot!(x, y, linewidth=3, label="", linecolor=:black)
    end
end
scatter!(Nodes.x, Nodes.y, label="", mc=:blue, ms=5, ma=2)

1.5 BFF加速

这里做一些模型优化，首先松弛调depot的约束，其次每个点只考虑附近的K个点，模型如下：

using DataFrames, Random, TSPLIB,Distances, PyPlot, ProgressBars
plt.figure(figsize=(20, 10))
function plot_tree(tsp,mst)
    s = tsp.nodes;scatter(s[:,1],s[:,2]);for m in mst;plot(s[[m...],1],s[[m...],2],c="b");end #for p in 1:tsp.dimension;text(s[p,1],s[p,2],p);end;
end
function greedyConstruct(dist)
    path = [1];for _ in 1:size(dist)[1];c_dist = dist[path[end],:];c_dist[path].=MaxNum;next_id = findmin(c_dist)[2];push!(path,next_id);end
    return path
end
neighbour_num = 15
MaxNum = 1000000
tsp=readTSPLIB(:d493)
println("Generate branches")
c = tsp.dimension
xpos = tsp.nodes[:,1]; ypos = tsp.nodes[:,2];
Nodes = DataFrames.DataFrame(;k = 1:c,x = xpos,y = ypos,);
branch_list = [Set() for i in 1:c]
for i in tqdm(1:c)
    dist = [floor.(0.5+euclidean(tsp.nodes[i,:],tsp.nodes[j,:])) for j in 1:c];dist[i] = MaxNum
    for j in partialsortperm(dist, 1:neighbour_num)
        if i<j
            push!(branch_list[i],(i,j,dist[j]))
        else
            push!(branch_list[j],(j,i,dist[j]))
        end
    end
end
println("Generate greedy route")
dist_matrix = [floor.(0.5+euclidean(tsp.nodes[i,:],tsp.nodes[j,:])) for i in 1:tsp.dimension,j in 1:tsp.dimension];
path = greedyConstruct(dist_matrix);
for idx in tqdm(1:length(path))
    i = path[idx]
    j = path[idx==length(path) ? 1 : idx+1]
    if i<j
        push!(branch_list[i],(i,j,dist_matrix[i,j]))
    else
        push!(branch_list[j],(j,i,dist_matrix[i,j]))