运筹系列85：求解大规模tsp问题的julia代码

1. 大规模tsp问题的挑战

数学模型和精确解法见《运筹系列65：TSP问题的精确求解法概述》和《运筹系列80:使用Julia精确求解tsp问题》：

@variable(m, x[1:n,1:n], Bin,Symmetric) # 0-1约束
@objective(model, Min, sum(x.*distmat)/2)
@constraint(model, [i = 1:n], sum(x[i, :]) == 2)
@constraint(model, [j = 1:n], sum(x[:, j]) == 2)
for subset in subsets # 防subtour约束，需要遍历所有的子集合
	@constraint(model,sum(x[subset[i],subset[i]])<=2*length(subset[i])-2) 
end
optimize!(m)
objective_value(model)

主要的挑战包括：

1. 求解整数规划比较耗时，并且遍历子集难度很大，大规模问题时需要自定义高效的cut和price策略
2. 变量数量过多，需要使用列生成方法来减少求解变量

2. 使用列生成减少求解变量

关于列生成的列子，可以参考《运筹系列8：Set partitioning问题的列生成方法》。对于tsp问题，我们使用列生成的步骤是：

1. 首先将问题变量限制在每个点最邻近的k条边上，求解限制主问题
2. 迭代求解子问题（sub problem）。如果目标函数最优值＜0，就将新生成的列yk+1转入基变量，生成新的限制主问题进行求解。如此往复，直至子问题的目标函数≥0。

我们以u159问题为例，原始模型一共有25281个变量：

using TSPLIB, LinearAlgebra, JuMP,HiGHS,PyPlot
#读取数据
tsp = readTSPLIB(:u159)
n = tsp.dimension
distmat = [round(Int64,norm(tsp.nodes[i,:] - tsp.nodes[j,:])) for i in 1:n, j in 1:n];
for i in 1:n;distmat[i,i] = 10^9;end

#完整模型
model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, 0<=x[i in 1:n,j in js[i]]<=1)
@info "variable size:",length(x)
@objective(model, Min, sum([x[i,j]*distmat[i,j]/2 for i in 1:n for j in js[i]]))
c = @constraint(model, [i = 1:n], sum([x[i, j] for j in js[i]]) >= 2)
e = @constraint(model, [i = 1:n], sum([x[j, i] for j in js[i]]) >= 2)
optimize!(model)
@info "objective:",objective_value(model)

主问题限制在每个点最近的3条边中，此时只有604个变量，第一轮列生成后为686个变量（此时已经得到最优解了），最终经过7轮列生成迭代，变量个数为716：

# js为初始下标合集，我们将变量下标限制在这个集合内。注意需要保持对称性：
idx = Index(2)
add(idx, tsp.nodes)
c = search(idx,tsp.nodes,4)[2];
js = Vector{Vector{Int}}()
for i in 1:n
   push!(js,c[i,2:end])
end
for i in 1:n
   for j in c[i,2:end]
       union!(js[j],i)
   end
end

求解主问题，然后找到检验数<0的列。令对偶变量$p=C_B{B}^{-1}$, 则检验数$\sigma =C_N-p\ast N$ 。我们将所有检验数小于0的列进基，迭代直至an ==0，此部分完整代码如下：

an = 1
while an >0
    model = Model(HiGHS.Optimizer)
    set_silent(model)
    @variable(model, 0<=x[i in 1:n,j in js[i]]<=1)
    @info "variable size:",length(x)
    @objective(model, Min, sum([x[i,j]*distmat[i,j]/2 for i in 1:n for j in js[i]]))
    c = @constraint(model, [i = 1:n], sum([x[i, j] for j in js[i]]) >= 2)
    e = @constraint(model, [i = 1:n], sum([x[j, i] for j in js[i]]) >= 2)
    optimize!(model)
    @info objective_value(model)
    p1 = [dual(c[i]) for i in 1:n]
    p2 = [dual(e[i]) for i in 1:n];
    an = 0
    for i in 1:n
        for j in setdiff(1:n,js[i])
            sigma = distmat[i,j]/2 - p1[i]-p2[j]
            if sigma < 0;push!(js[i],j);push!(js[j],i);an+=2;end
        end
    end
    @info "add size:",an
end

3. 使用行生成添加约束

我们观察att48使用lp求解后的结果：

有很多子环，需要添加去除环的约束。
首先添加subtour约束，去除所有的子环约束：

using Graphs
paths = strongly_connected_components(SimpleDiGraph(round.(x_val,digits = 2)))
if length(paths)>1
    for path in paths
        @constraint(model,sum(x[path,path])<=2*length(path)-2)
    end
    optimize!(model)
    x_val = round.(JuMP.value.(x),digits = 2)
    paths = strongly_connected_components(SimpleDiGraph(x_val))
end

迭代求解，并绘制结果如下：

然后添加comb约束，约束介绍参见《运筹系列65：TSP问题的精确求解法概述》，julia代码参见《运筹系列80: 使用Julia精确求解tsp问题》的4.2节，求解后得到两个comb：

然后添加进约束，迭代求解：

迭代得到如下图：

4. 使用分枝定界求解整数规划问题

我们这里使用queue存储分枝节点，按照最简单的下标顺序，对所有非整数变量进行分枝。完整代码如下：

using TSPLIB,JuMP,DataStructures,GLPK,Graphs
struct Node
    x::Vector{CartesianIndex{2}} 
    b::Vector{Int} # 前面两个用于b&b
    paths::Vector{Vector{Int}} # 用于subtour eliminate
end

# 线性规划求解器
function linPro(node::Node)
    m = Model(GLPK.Optimizer)
    @variable(m, 0<=x[1:n, 1:n]<=1, Symmetric)
    @objective(m, Min, sum(x.*distmat)/2)
    @constraint(m, [i = 1:n], sum(x[i, :]) == 2)
    for i in 1:length(node.x)
        @constraint(m,x[node.x[i]]==node.b[i])
    end
    for path in node.paths
        @constraint(m,sum(x[path,path])<=2*length(path)-2)
    end
    optimize!(m)
    if Int(termination_status(m))==1
        r_x = round.(JuMP.value.(x),digits = 2)
        r_val = objective_value(m)
        return r_x,r_val
    end
    return nothing,nothing
end;

# 简单的分枝定界算法求整数解
function integerPro(paths,t=1.0E-12)
    ub = 10^9
    ub_x = nothing
    PQ = Queue{Node}()
    enqueue!(PQ,Node([],[],paths)) 
    while length(PQ)>0
        v = dequeue!(PQ)
        x,res = linPro(v)
        ind = findfirst(xi->xi>t && xi<1-t ,x)
        if ind == nothing && res < ub
            ub = res
            ub_x = copy(x)
        elseif res < ub
            enqueue!(PQ,Node([v.x;ind],[v.b;0],paths))
            enqueue!(PQ,Node([v.x;ind],[v.b;1],paths))
        end
    end
    return ub_x,ub
end;

tsp = readTSPLIB(:att48) # 使用att48例子
x_val,res = integerPro([])
paths = strongly_connected_components(SimpleDiGraph(round.(x_val,digits = 2)))
paths_list = []
# 行生成算法，不断添加subtour elimination constraints
while length(paths)>1
    @info length(paths)
    paths_list = [paths_list;paths]
    x_val,res = integerPro(paths_list)
    paths = strongly_connected_components(SimpleDiGraph(round.(x_val,digits = 2)))
end

plot_matrix(tsp, x_val)