# Copyright (c) 2022 Daniel Schermer and contributors #src
# #src
# Permission is hereby granted, free of charge, to any person obtaining a copy #src
# of this software and associated documentation files (the "Software"), to deal #src
# in the Software without restriction, including without limitation the rights #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell #src
# copies of the Software, and to permit persons to whom the Software is #src
# furnished to do so, subject to the following conditions: #src
# #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software. #src
# #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE #src
# SOFTWARE. #src
# # [Traveling Salesperson Problem](@id tsp_lazy)
# **This tutorial was originally contributed by Daniel Schermer.**
# This tutorial describes how to implement the
# [Traveling Salesperson Problem](https://en.wikipedia.org/wiki/Travelling_salesman_problem)
# in JuMP using solver-independent lazy constraints that dynamically separate
# subtours. To be more precise, we use lazy constraints to cut off infeasible
# subtours only when necessary and not before needed.
# It uses the following packages:
using JuMP
import GLPK
import Random
import Plots
# ## [Mathematical Formulation](@id tsp_model)
# Assume that we are given a complete graph $\mathcal{G}(V,E)$ where $V$ is the
# set of vertices (or cities) and $E$ is the set of edges (or roads). For each
# pair of vertices $i, j \in V, i \neq j$ the edge $(i,j) \in E$ is associated
# with a weight (or distance) $d_{ij} \in \mathbb{R}^+$.
# For this tutorial, we assume the problem to be symmetric, that is,
# $d_{ij} = d_{ji} \, \forall i,j \in V$.
# In the Traveling Salesperson Problem, we are tasked with finding a tour with
# minimal length that visits every vertex exactly once and then returns to the
# point of origin, that is, a *Hamiltonian cycle* with minimal weight.
#
# To model the problem, we introduce a binary variable,
# $x_{ij} \in \{0,1\} \; \forall i, j \in V$, that indicates if edge $(i,j)$ is
# part of the tour or not. Using these variables, the Traveling Salesperson
# Problem can be modeled as the following integer linear program.
# ### [Objective Function](@id tsp_objective)
# The objective is to minimize the length of the tour (due to the assumed
# symmetry, the second sum only contains $i 1
pop!(unvisited, current)
end
neighbors =
[j for (i, j) in edges if i == current && j in unvisited]
end
if length(this_cycle) < length(shortest_subtour)
shortest_subtour = this_cycle
end
end
return shortest_subtour
end
# Let us declare a helper function `selected_edges()` that will be repeatedly
# used in what follows.
function selected_edges(x::Matrix{Float64}, n)
return Tuple{Int,Int}[(i, j) for i in 1:n, j in 1:n if x[i, j] > 0.5]
end
# Other helper functions for computing subtours:
subtour(x::Matrix{Float64}) = subtour(selected_edges(x, size(x, 1)), size(x, 1))
subtour(x::AbstractMatrix{VariableRef}) = subtour(value.(x))
# ### Iterative method
# An iterative way of eliminating subtours is the following.
# However, it is reasonable to assume that this is not the most efficient way:
# whenever a new subtour elimination constraint is added to the model, the
# optimization has to start from the very beginning.
# That way, the solver will repeatedly discard useful information encountered
# during previous solves (for example, all cuts, the incumbent solution, or lower
# bounds).
# !!! info
# Note that, in principle, it would also be feasible to add all subtours
# (instead of just the shortest one) to the model. However, preventing just
# the shortest cycle is often sufficient for breaking other subtours and
# will keep the model size smaller.
iterative_model = build_tsp_model(d, n)
optimize!(iterative_model)
time_iterated = solve_time(iterative_model)
cycle = subtour(iterative_model[:x])
while 1 < length(cycle) < n
println("Found cycle of length $(length(cycle))")
S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
@constraint(
iterative_model,
sum(iterative_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
)
optimize!(iterative_model)
global time_iterated += solve_time(iterative_model)
global cycle = subtour(iterative_model[:x])
end
objective_value(iterative_model)
#-
time_iterated
# As a quick sanity check, we visualize the optimal tour to verify that no
# subtour is present:
function plot_tour(X, Y, x)
plot = Plots.plot()
for (i, j) in selected_edges(x, size(x, 1))
Plots.plot!([X[i], X[j]], [Y[i], Y[j]]; legend = false)
end
return plot
end
plot_tour(X, Y, value.(iterative_model[:x]))
# ### Lazy constraint method
# A more sophisticated approach makes use of _lazy constraints_. To be more
# precise, we do this through the `subtour_elimination_callback()` below, which
# is only run whenever we encounter a new integer-feasible solution.
lazy_model = build_tsp_model(d, n)
function subtour_elimination_callback(cb_data)
status = callback_node_status(cb_data, lazy_model)
if status != MOI.CALLBACK_NODE_STATUS_INTEGER
return # Only run at integer solutions
end
cycle = subtour(callback_value.(cb_data, lazy_model[:x]))
if !(1 < length(cycle) < n)
return # Only add a constraint if there is a cycle
end
println("Found cycle of length $(length(cycle))")
S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
con = @build_constraint(
sum(lazy_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
)
MOI.submit(lazy_model, MOI.LazyConstraint(cb_data), con)
return
end
set_attribute(
lazy_model,
MOI.LazyConstraintCallback(),
subtour_elimination_callback,
)
optimize!(lazy_model)
objective_value(lazy_model)
# This finds the same optimal tour:
plot_tour(X, Y, value.(lazy_model[:x]))
# Surprisingly, for this particular model with GLPK, the solution time is worse
# than the iterative method:
time_lazy = solve_time(lazy_model)
# In most other cases and solvers, however, the lazy time should be faster than
# the iterative method.