运筹系列80: 使用Julia精确求解tsp问题

Posted 2023-03-10 IE06

1. 基础模型

使用flow约束构建基础模型

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)

结果如下：

2. subtour约束

接下来不断增加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才求出最优解

3. 调用包求解

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)

用时12.7s，主要的优化是每次同时把所有的subtour都添加进来。我们修改后的代码如下：

using JuMP, GLPK, Distances

function get_cycles(perm_matrix)
    N = size(perm_matrix, 1)
    remaining_inds = Set(1:N)
    cycles = VectorInt[]
    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)

4. 仅使用lp求解

4.1 使用subtour约束

第一步，将所有的subtour消除：

using TravelingSalesmanExact, GLPK,TSPLIB,Distances,JuMP
function get_tours(perm_matrix, thresh = 0.1)
    N = size(perm_matrix, 1)
    remaining_inds = Set(1:N)
    tours = VectorInt[]
    while length(remaining_inds) > 0
        tour = find_tour(perm_matrix, first(remaining_inds), thresh)
        push!(tours, tour)
        setdiff!(remaining_inds, tour)
    end
    return tours
end

function find_tour(perm_matrix, starting_ind, thresh)
    tour = []
    search = [starting_ind]
    n = size(perm_matrix,1)
    while length(search)>0
        v=pop!(search)
        if !(v in tour)
            push!(tour,v)
            for i in 1:n
                if perm_matrix[v,i]>=thresh && !(i in tour)
                    push!(search, i)
                end
            end
        end
    end
    return tour
end

tsp = readTSPLIB(:bays29); # bays29
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, 0<=x[1:N,1:N]<=1, Symmetric)
@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
@time optimize!(m)
thresh = 0.1
tours = get_tours(JuMP.value.(x), thresh)
while length(tours) > 1
    for tour in tours
        @constraint(m,sum(x[tour,tour])<=2*size(tour,1)-2)
    end
    @time optimize!(m)
    tours = get_tours(JuMP.value.(x), thresh)
end
# println(JuMP.objective_value(m))
plot_matrix(tsp,round.(JuMP.value.(x); digits=1))

最终结果如下：

4.2 使用comb约束

下面仅作例子：

function find_comb(sol, odds, starting_ind = 1)
    comb = VectorVectorInt()
    teeths = VectorVectorInt()
    tour = VectorInt()
    search = [starting_ind]
    n = size(odds,1)
    N = size(sol,1)
    while length(search)>0
        v = odds[pop!(search)]
        if !(v in tour)
            push!(tour,v)
            for i in 1:n
                if sol[v,odds[i]]>0.01 && sol[v,odds[i]]<0.99  && !(odds[i] in tour)
                    push!(search, i)
                end
            end
            teeth = [v]
            for i in 1:N
                if sol[v,i]==1
                    push!(teeth, i)
                end
            end
            if length(teeth)>1
                push!(teeths,teeth)  
            end
        end
    end
    for teeth in teeths
        if !(teeth[2] in tour)
            push!(comb,teeth)
        end
    end
    push!(comb,tour) 
    return comb
end

function get_combs(sol)
    N = size(sol, 1)
    odds = VectorInt()
    for i in 1:N
        for j in 1:N
            if 0< sol[i,j] && sol[i,j]<1
                if !(i in odds);push!(odds,i);end
                if !(j in odds);push!(odds,j);end
            end
        end
    end
    combs = VectorVectorVectorInt()
    while length(odds) > 0
        comb = find_comb(sol, odds)
        push!(combs, comb)
        setdiff!(odds, comb[end])
    end
    return combs
end

combs = get_combs(JuMP.value.(x))
while length(combs)>0
	for comb in combs
	    k = length(comb)-1
	    @constraint(m,sum([sum(x[teeth,setdiff(1:N,teeth)]) for teeth in comb])>=3*k+1)
	end
	optimize!(m)
	combs = get_combs(JuMP.value.(x))
end

plot_matrix(tsp,round.(JuMP.value.(x); digits=1))