/
tsp.jl
209 lines (189 loc) · 6.21 KB
/
tsp.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# tsp.jl
#
# Solves the travelling salesman problem using integer programming and
# lazy generation of the subtour elimination constraints.
#############################################################################
using JuMP
using Gurobi
using Base.Test
# extractTour
# Given a n-by-n matrix representing the solution to an undirected TSP,
# extract the tour as a vector
# Input:
# n Number of cities
# sol n-by-n 0-1 symmetric matrix representing solution
# Output:
# tour n+1 length vector of tour, starting and ending at 1
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] >= 1-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
# findSubtour
# Given a n-by-n matrix representing solution to the relaxed
# undirected TSP problem, find a set of nodes belonging to a subtour
# 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
subtour[cur_city] = true
subtour_length = 1
while true
# Find next node that we haven't yet visited
found_city = false
for j = 1:n
if !subtour[j]
if sol[cur_city, j] >= 1 - 1e-6
# Arc to unvisited city, follow it
cur_city = j
subtour[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
# solveTSP
# Given a matrix of city locations, solve the TSP
# Inputs:
# n Number of cities
# cities n-by-2 matrix of (x,y) city locations
# Output:
# path Vector with order to cities are visited in
function solveTSP(n, cities)
# 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
# Create a model that will use Gurobi to solve
m = Model(solver=GurobiSolver())
# 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 upper triangle
@variable(m, x[1:n,1:n], Bin)
# Minimize length of tour
@objective(m, Min, sum(dist[i,j]*x[i,j] for i=1:n for j=i:n))
# Make x_ij and x_ji be the same thing (undirectional)
# Don't allow self-arcs
for i = 1:n
@constraint(m, x[i,i] == 0)
for j = (i+1):n
@constraint(m, x[i,j] == x[j,i])
end
end
# 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
function subtour(cb)
# Optional: display tour starting at city 1
println("----\nInside subtour callback")
println("Current tour starting at city 1:")
print(extractTour(n, getvalue(x)))
# Find any set of cities in a subtour
subtour, subtour_length = findSubtour(n, getvalue(x))
if subtour_length == n
# This "subtour" is actually all cities, so we are done
println("Solution visits all cities")
println("----")
return
end
# Subtour found - add lazy constraint
# We will build it up piece-by-piece
arcs_from_subtour = zero(AffExpr)
for i = 1:n
if !subtour[i]
# If this city isn't in subtour, skip it
continue
end
# Want to include all arcs from this city, which is in
# the subtour, to all cities not in the subtour
for j = 1:n
if i == j
# Self-arc
continue
elseif subtour[j]
# Both ends in same subtour
continue
else
# j isn't in subtour
arcs_from_subtour += x[i,j]
end
end
end
# Add the new subtour elimination constraint we built
println("Adding subtour elimination cut")
println("----")
@lazyconstraint(cb, arcs_from_subtour >= 2)
end # End function subtour
# Solve the problem with our cut generator
addlazycallback(m, subtour)
solve(m)
# Return best tour
return extractTour(n, getvalue(x))
end # end solveTSP
# Create a simple instance that looks like
# + +
# + +
# + +
# The optimal tour is obvious, but the initial solution will be
# /--+ +--\
# + | +
# \--+ +--/
n = 6
cities = [ 50 200;
100 100;
100 300;
500 100;
500 300;
550 200]
tour = solveTSP(n, cities)
println("Solution: ")
println(tour)