using JuMP
# Adapt the following two lines to the solver of your choice.
using Gurobi
function getSolver()
solver = GurobiSolver(OutputFlag=false)
end
n_cities = 7;
cost_vector = vec([2 4 1 3 7 1 7]);
distance_matrix = [
[0 4 3 3 5 2 5];
[4 0 5 3 3 4 7];
[3 5 0 4 6 0 4];
[3 3 4 0 4 4 6];
[5 3 6 4 0 5 8];
[2 4 0 4 5 0 3];
[5 7 4 6 8 3 0]
] # Check that it is symmetric: distance_matrix' == distance_matrix.
function showEdges(n_cities, sol)
println("This solution has the following edges: ")
for i = 1:n_cities
for j = (i+1):n_cities
if sol[i, j] > .99
println(" ", i, " to ", j)
end
end
end
end
function getPCTSP(n_cities::Int, cost_vector::Array{Int64,1}, matrix::Array{Int64,2})
m = Model(solver=getSolver())
@defVar(m, 0 <= x[i=1:n_cities, j=(i+1):n_cities] <= 1) # Not integer: linear relaxation. Upper triangular, zero diagonal!
@defVar(m, 0 <= y[1:n_cities] <= 1) # Not integer: linear relaxation.
@setObjective(m, Max, dot(cost_vector, y) - sum{distance_matrix[i, j] * x[i, j], i=1:n_cities, j=(i+1):n_cities})
for i in 1:n_cities
@addConstraint(m, sum{x[i, j], j=(i+1):n_cities} + sum{x[j, i], j=1:(i-1)} == 2 * y[i]) # Don't forget the matrix is triangular: from i to all the others, then from all the others to i.
for j in (i+1):n_cities
# Trivial GSECs.
@addConstraint(m, x[i, j] <= y[i])
@addConstraint(m, x[i, j] <= y[j])
end
end
@addConstraint(m, y[1] == 1)
return m, x, y # Solution ensured to be integer! Except if more constraints...
end
# First step.
m1, m1_x, m1_y = getPCTSP(n_cities, cost_vector, distance_matrix)
solve(m1)
showEdges(n_cities, getValue(m1_x)) # See the two tours: (1, 3, 6, 7, 1) and (2, 4, 5, 2).