运筹系列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)

主要的挑战包括：

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

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

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

首先将问题变量限制在每个点最邻近的k条边上，求解限制主问题
迭代求解子问题（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:npush!(js,c[i,2:end])
end
for i in 1:nfor j in c[i,2:end]union!(js[j],i)end
end

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

an = 1
while an >0model = 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 = 0for i in 1:nfor 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;endendend@info "add size:",an
end

3. 使用行生成添加约束

我们观察att48使用lp求解后的结果：
在这里插入图片描述

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

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

迭代求解，并绘制结果如下：
在这里插入图片描述