The basic problem is that we are given a set of co-ordinates $(a_i, b_i)_{i=1}^n$ corresponding to a set of $n$ nodes in a graph, and we need to decide which $n$ of the $n^2$ arcs in the graph we will use to tour the graph, wherein we visit each node exactly once, at minimum total cost. Let's formulate this as an integer program!
Let the cities be indexed from 1 to N. Let $d_{ij}$ be the distance between city $i$ and city $j$.
Decision variables: $x_{ij}=\begin{cases} 1,\quad \text{if city $i$ and city $j$ are adjacent in the shortest tour}\\ 0,\quad \text{otherwise.}\end{cases}$
N.B. $x_{ij}$ and $x_{ji}$ are redundant ($x_{ij}=x_{ji}$), so we only define the variable $x_{ij}$ for $i < j$. Then we can formulate the following integer program.
$$ \underset{x}{\min}\ \sum_{i=1}^{N-1}\sum_{j=i+1}^N d_{ij}x_{ij} \\ \text{s.t.}\quad \sum_{j=i+1}^N x_{ij} + \sum_{j=1}^{i-1}x_{ji} = 2 \quad\forall i, 1\le i \le N \\ x_{ij}\in\{0,1\}\quad\forall i,j \text{ s.t. } 1\le i < j \le N $$However, if we attempt to solve this problem, we will find that we get an infeasible solution, with multiple disjoint subtours.
Yikes! Our formulation is missing something! What are some potential ways to fix it?
One common way is subtour elimination constraints, to prevent the final solution from having any small cycles, i.e. cycles that do not include all the nodes.
Given a subtour $S\subset \{1,\ldots,N\}$, a subtour elimination constraint looks like: $$\sum_{i\in S} \left(\sum_{j\notin S, j > i}x_{ij}+\sum_{j\notin S, j < i}x_{ji}\right) \ge 2.$$
As $N$ grows larger, the number of subtour elimination constraints grows exponentially. It is therefore impractical to add all of these constraints into the model.
Instead, we generate these constraints lazily. Every time Gurobi has an incumbent solution, we find the shortest subtour in the solution, and add a lazy constraint eliminating this particular subtour.
We begin by loading the required modules.
using JuMP, Gurobi, Test, LinearAlgebra, DelimitedFiles, Random
Next, we define a function to extract a n+1 dimensional vector representing a tour from an n x n symmetric matrix representing a solution provided by a solver.
function extractTour(n, sol)
tour = [1] # Start at city 1 always
cur_city = 1
while true
# Look for first arc out of current city
for j = 1:n
if sol[cur_city,j] >= 0.5-1e-6
# Found next city
push!(tour, j)
# Don't ever use this arc again
sol[cur_city, j] = 0.0
sol[j, cur_city] = 0.0
# Move to next city
cur_city = j
break
end
end
# If we have come back to 1, stop
if cur_city == 1
break
end
end # end while
return tour
end
extractTour (generic function with 1 method)
Next, we define a function which acts as a seperation oracle, which is a fancy way of saying that it either identifies a subtour which should be banned from the set of all possible solutions, or decides that the current solution is optimal.
# Input:
# n Number of cities
# sol n-by-n 0-1 symmetric matrix representing solution
# Outputs:
# subtour n length vector of booleans, true iff in a particular subtour
# subtour_length Number of cities in subtour (if n, no subtour found)
function findSubtour(n, sol)
# Initialize to no subtour
subtour = fill(false,n)
#=
# Always start looking at city 1
cur_city = 1
=#
# Start looking at a random city: much faster because we explore different subtours
cur_city=rand(1:n)
subtour[cur_city] = true
subtour_length = 1
while true
# Find next node that we haven't yet visited
found_city = false
indices = shuffle(1:n)
for j = 1:n
if !subtour[indices[j]]
if sol[cur_city, indices[j]] >= 1 - 1e-6
# Arc to unvisited city, follow it
cur_city = indices[j]
subtour[indices[j]] = true
found_city = true
subtour_length += 1
break # Move on to next city
end
end
end
if !found_city
# We are done
break
end
end
return subtour, subtour_length
end
findSubtour (generic function with 1 method)
Next, we define a function which solves TSP, given a matrix of city locations, using an optimization solver.
# Inputs:
# cities n-by-2 matrix of (x,y) city locations
# Output:
# path Vector with order to cities are visited in
function solveTSP(cities; time_limit=30.0)
n = size(cities)[1]
# Calculate pairwise distance matrix
dist = zeros(n, n)
for i = 1:n
for j = i:n
d = norm(cities[i,1:2] - cities[j,1:2])
dist[i,j] = d
dist[j,i] = d
end
end
m = Model(Gurobi.Optimizer)
set_optimizer_attribute(m, "TimeLimit", time_limit)
# x[i,j] is 1 iff we travel between i and j, 0 otherwise.
# Although we define all n^2 variables, we will only use the (strict) upper triangle.
@variable(m, x[1:n,1:n], Bin)
# Minimize total length of tour
@objective(m, Min, dot(dist, x))
# Make x_ij and x_ji be the same thing (undirectional TSP)
@constraint(m, x.==x')
# Don't allow self-arcs, by ensuring diagonal is vector of 0s
@constraint(m, diag(x).==zeros(n))
# We must enter and leave every city once and only once
for i = 1:n
@constraint(m, sum(x[i,j] for j=1:n) == 2)
end
# Lazy constraint
lazy_called = false
function subtour(cb)
lazy_called = true
# Find any set of cities in a subtour
x_val = callback_value.(Ref(cb), x) # In previous versions, you'd simply use getvalue(x)
# println(x_val)
subtour, subtour_length = findSubtour(n, x_val)
if subtour_length == n
# This "subtour" is actually all cities, so we are done with this node of the branch and bound tree
return
end
# Subtour found - add lazy constraint
arcs_from_subtour = zero(AffExpr)
for i = 1:n
if subtour[i]
# If this city isn't in subtour, skip it
for j = 1:n
# Want to include all arcs from this city, which is in the subtour,
# to all cities not in the subtour
if (i !=j) && !(subtour[j])
# j isn't in subtour
arcs_from_subtour += x[i,j]
end
end
end
end
# Add the subtour elimination constraint
con = @build_constraint(arcs_from_subtour >= 2)
# Submit built constraint to model via MOI
MOI.submit(m, MOI.LazyConstraint(cb), con)
# Here's how you'd do this in previous JuMP versions:
# @lazyconstraint(cb, arcs_from_subtour >= 2)
end
MOI.set(m, MOI.LazyConstraintCallback(), subtour)
optimize!(m)
return extractTour(n, value.(x))
end
solveTSP (generic function with 1 method)
Next, we define a function to plot the solution
plot_instance(pts) = plot(x = pts[1,:], y = pts[2,:], Geom.point, Guide.xlabel(nothing), Guide.ylabel(nothing))
function plot_solution(pts, path, extras = [])
ptspath = pts[:,path]
plot(x = ptspath[1,:], y = ptspath[2,:], Geom.point, Geom.path, Guide.xlabel(nothing), Guide.ylabel(nothing), extras...)
end
plot_solution (generic function with 2 methods)
Next, we solve a small 6 city example (which you might recognize from the diagram above) to verify correctness of our code.
n = 6
cities = [50 200;
100 100;
100 300;
500 100;
500 300;
550 200]
tour = solveTSP(cities)
#plot_solution(cities', tour)
Set parameter TokenServer to value "flexlm" Set parameter TimeLimit to value 30 Set parameter TimeLimit to value 30 Set parameter LazyConstraints to value 1 Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64) Thread count: 48 physical cores, 96 logical processors, using up to 32 threads Optimize a model with 48 rows, 36 columns and 102 nonzeros Model fingerprint: 0x7b9aa57c Variable types: 0 continuous, 36 integer (36 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+02, 5e+02] Bounds range [0e+00, 0e+00] RHS range [2e+00, 2e+00] Found heuristic solution: objective 4189.4701349 Presolve removed 42 rows and 21 columns Presolve time: 0.00s Presolved: 6 rows, 15 columns, 30 nonzeros Variable types: 0 continuous, 15 integer (15 binary) Root relaxation: objective 1.694427e+03, 4 iterations, 0.00 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time * 0 0 0 2494.4271910 2494.42719 0.00% - 0s Cutting planes: Lazy constraints: 1 Explored 1 nodes (8 simplex iterations) in 0.04 seconds (0.00 work units) Thread count was 32 (of 96 available processors) Solution count 2: 2494.43 4189.47 Optimal solution found (tolerance 1.00e-04) Best objective 2.494427191000e+03, best bound 2.494427191000e+03, gap 0.0000% User-callback calls 117, time in user-callback 0.04 sec
7-element Vector{Int64}: 1 2 4 6 5 3 1
What's the quickest tour around the 48 US state capitals in the mainland US?
# Source: https://people.sc.fsu.edu/~jburkardt/datasets/tsp/att48.tsp
n=48
citiesdata=[6734 1453;2233 10;5530 1424;401 841;3082 1644;7608 4458;7573 3716;7265 1268;6898 1885;1112 2049;5468 2606;
5989 2873;4706 2674;4612 2035;6347 2683;6107 669;7611 5184;7462 3590;7732 4723;5900 3561;4483 3369;6101 1110;5199 2182;
1633 2809;4307 2322;675 1006;7555 4819;7541 3981;3177 756;7352 4506;7545 2801;3245 3305;6426 3173;4608 1198;23 2216;
7248 3779;7762 4595;7392 2244;3484 2829;6271 2135;4985 140;1916 1569;7280 4899;7509 3239;10 2676;6807 2993;5185 3258;3023 1942]
tour = solveTSP(citiesdata);
# plot_solution(citiesdata', tour)
Set parameter TokenServer to value "flexlm" Set parameter TimeLimit to value 30 Set parameter TimeLimit to value 30 Set parameter LazyConstraints to value 1 Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64) Thread count: 48 physical cores, 96 logical processors, using up to 32 threads Optimize a model with 2400 rows, 2304 columns and 6864 nonzeros Model fingerprint: 0xc697e34d Variable types: 0 continuous, 2304 integer (2304 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+02, 8e+03] Bounds range [0e+00, 0e+00] RHS range [2e+00, 2e+00] Presolve removed 2352 rows and 1176 columns Presolve time: 0.00s Presolved: 48 rows, 1128 columns, 2256 nonzeros Variable types: 0 continuous, 1128 integer (1128 binary) Root relaxation: objective 6.333851e+04, 78 iterations, 0.00 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 63338.5149 0 14 - 63338.5149 - - 0s H 0 0 75451.043760 63338.5149 16.1% - 0s H 0 0 75399.577115 63338.5149 16.0% - 0s 0 0 64641.6215 0 6 75399.5771 64641.6215 14.3% - 0s 0 0 65133.1956 0 12 75399.5771 65133.1956 13.6% - 0s 0 0 65472.9195 0 13 75399.5771 65472.9195 13.2% - 0s 0 0 66807.1344 0 10 75399.5771 66807.1344 11.4% - 0s H 0 0 73613.524473 66807.1344 9.25% - 0s 0 0 66849.6523 0 14 73613.5245 66849.6523 9.19% - 0s 0 0 66849.6523 0 14 73613.5245 66849.6523 9.19% - 0s 0 0 66849.6523 0 10 73613.5245 66849.6523 9.19% - 0s 0 0 66899.0515 0 18 73613.5245 66899.0515 9.12% - 0s H 0 0 69831.516009 66899.0515 4.20% - 0s H 0 0 67413.691982 66899.0515 0.76% - 0s 0 0 66939.3648 0 10 67413.6920 66939.3648 0.70% - 0s 0 0 66939.3648 0 17 67413.6920 66939.3648 0.70% - 0s H 0 0 67047.417015 66939.3648 0.16% - 0s Cutting planes: MIR: 1 Zero half: 3 Explored 1 nodes (309 simplex iterations) in 0.31 seconds (0.07 work units) Thread count was 32 (of 96 available processors) Solution count 6: 67047.4 67413.7 69831.5 ... 75451 Optimal solution found (tolerance 1.00e-04) Best objective 6.704741701487e+04, best bound 6.704741701487e+04, gap 0.0000% User-callback calls 314, time in user-callback 0.17 sec
Let's try to solve a TSP with 200 cities using the above code.
ludata=readdlm("bcl380tsp.txt")
n = 200 #38
ludata = ludata[1:n, 2:3] #drop index
tour_mip = solveTSP(ludata, time_limit=120.)
Set parameter TokenServer to value "flexlm" Set parameter TimeLimit to value 120 Set parameter TimeLimit to value 120 Set parameter LazyConstraints to value 1 Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64) Thread count: 48 physical cores, 96 logical processors, using up to 32 threads Optimize a model with 40400 rows, 40000 columns and 119800 nonzeros Model fingerprint: 0x2cfe5ede Variable types: 0 continuous, 40000 integer (40000 binary) Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+00, 2e+02] Bounds range [0e+00, 0e+00] RHS range [2e+00, 2e+00] Presolve removed 40200 rows and 20100 columns Presolve time: 0.04s Presolved: 200 rows, 19900 columns, 39800 nonzeros Variable types: 0 continuous, 19900 integer (19900 binary) Root relaxation: objective 1.619982e+03, 281 iterations, 0.01 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 1619.98221 0 24 - 1619.98221 - - 0s 0 0 1644.98988 0 37 - 1644.98988 - - 1s 0 0 1647.49528 0 59 - 1647.49528 - - 1s 0 0 1649.16532 0 6 - 1649.16532 - - 1s 0 0 1649.16532 0 6 - 1649.16532 - - 1s 0 2 1649.16532 0 6 - 1649.16532 - - 1s 31 64 1677.95395 5 26 - 1676.46135 - 12.4 5s 95 128 1692.03883 6 14 - 1678.76686 - 13.6 11s 159 192 1700.30701 7 8 - 1691.68859 - 19.6 22s 191 253 1700.30701 8 6 - 1692.03883 - 20.6 27s 252 310 1703.15592 9 6 - 1692.03883 - 17.8 38s 309 368 1700.92951 10 16 - 1692.03883 - 15.4 46s 367 400 1700.98772 11 12 - 1692.03883 - 14.3 54s 370 400 1719.20875 12 8 - 1692.03883 - 14.2 55s 392 400 1724.71078 8 24 - 1692.03883 - 13.9 60s 456 527 1703.60298 13 6 - 1692.03883 - 13.0 72s 526 581 1705.96270 14 15 - 1692.03883 - 12.2 85s 568 581 1722.49980 12 25 - 1692.03883 - 11.7 92s 582 630 1713.16866 16 6 - 1692.03883 - 11.5 95s 611 630 1736.58883 11 22 - 1692.03883 - 11.3 100s 670 714 1704.63399 16 22 - 1692.03883 - 10.9 109s 713 744 1706.69201 16 12 - 1692.03883 - 10.4 117s * 716 744 17 1723.3132254 1692.03883 1.81% 10.4 117s Cutting planes: Zero half: 21 Lazy constraints: 689 Explored 745 nodes (7969 simplex iterations) in 121.39 seconds (2.00 work units) Thread count was 32 (of 96 available processors) Solution count 1: 1723.31 Time limit reached Best objective 1.723313225417e+03, best bound 1.692038832216e+03, gap 1.8148% User-callback calls 2143, time in user-callback 114.00 sec
201-element Vector{Int64}: 1 14 13 12 35 34 33 32 31 52 50 49 11 ⋮ 17 45 44 56 58 55 43 42 53 41 40 1
Let's test our new solver! Notice that the solution reached is provably optimal (up to solver tolerances)! This also verifies that the heuristic solution obtained earlier is indeed of high-quality.
This material is adapted from previous versions of this course, which have been designed by numerous ORC students.
Some of the sources used to create this year's version include: