FICO Xpress Optimization Examples Repository
 FICO Optimization Community FICO Xpress Optimization Home

Solve Traveling Salesperson Problems using callbacks or delayed rows

Description

In this example, we show how to use the FICO Xpress R-interface to model and solve the famous Travelling Salesperson Problem (TSP). The TSP is the classic combinatorial optimization problem to find the shortest so-called tour through a complete graph.

Further explanation of this example: Xpress R Reference Manual

 tsp_r.zip [download all files]

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.
 tsp.R [download]

tsp.R

#####################################
# This file is part of the          #
# Xpress-R interface examples       #
#                                   #
#   (c) 2021 Fair Isaac Corporation #
#####################################
#' ---
#' title: "Finding a Shortest Tour"
#' subtitle: "The Travelling Salesperson Problem"
#' author: Gregor Hendel
#' date: Dec. 2020
#' ---
#'
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(results = "hold")

#'
## ----Load libraries-----------------------------------------------------------
suppressMessages(library(xpress))

# we use the subsets function from CombMSC
suppressMessages(library(CombMSC))

# load necessary libraries for plotting
suppressMessages(library(reshape2))
suppressMessages(library(magrittr))
suppressMessages(library(ggplot2))
suppressMessages(library(dplyr))
suppressMessages(library(tibble))

#'
#' In this example, we show how to use the FICO Xpress R-interface to model and
#' solve the famous **Travelling Salesperson Problem (TSP)**. The TSP is
#' the classic combinatorial optimization problem to find the shortest so-called
#' tour through a complete graph.
#'
#' This longer example discusses different approaches for solving TSPs with Xpress.
#' First we solve the problem without any subtour elimination constraints, which
#' will give an invalid solution. Then we solve the problem initialized with all
#' subtour elimination constraints, which does not perform well. Finally we show
#' how to use callbacks for adding subtour elimination constraints
#' dynamically during the search.
#'
#' We dive right into the mathematical formulation.
#'
#' # Mathematical formulation
#'
#' There are different variants how to formulate the TSP. We use the classical
#' formulation using binary variables and exponentially many constraints to
#' eliminate all so-called subtours.
#'
#' ## Graphs and Tours
#'
#' We assume we are given a complete Graph $K_n$ on $n$ vertices $V$. A graph is
#' called complete if it contains all $\binom{n}{2}$ possible edges, one between
#' each pair of nodes $v \neq w \in V \times V$. Denote the set of edges by $E$. We
#' further have a function that assigns a cost to each edge,
#'
#' $$#' d: E \rightarrow \mathbb{Q} #'$$
#'
#' A path in $K_n$ is a sequence of adjacent edges $(e_1,e_2,\dots,e_k)$ such that
#' consecutive edges have a common node, and no node is visited twice (except,
#' possibly, for the endpoints). A tour is a path of length $n$, which starts and
#' ends at the same node and in between visits each node exactly once. The cost of
#' a tour is the sum of the costs of its edges. The TSP demands to find a tour of
#' minimum cost.
#'
#' For the mathematical formulation, we introduce $n(n-1)/2$ binary variables
#' $x_{e}$ for each edge $e \in E$. Of course, $x_{e} = 1$ corresponds to $e$ being
#' part of the tour represented by a solution to our program.
#'
#' #' \begin{align} #' && \min \sum\limits_{e \in E} d(e) x_{e}\\ #' & \text{s.t.} &\sum\limits_{e \in \delta(v)} x_{e} &= 2 & \forall v\in V\\ #' & &\sum\limits_{e \in \delta(S)} x_e & \geq 2 & \forall \emptyset \neq S \subsetneq V\\ #' & & x_e &\in \{0,1\} & \forall e \in E #' \end{align} #'
#'
#' The first set of $n$ constraints requires that each node should have exactly two
#' adjacent edges in a solution. This set of constraints is often referred to as
#' "node-degree" constraints. The second set of constraints formulate that each
#' nonempty strict subset of nodes be entered at least once and left at least once.
#' We refer to this type of constraints as "subtour-elimination" constraints.
#'
#' It is easy to verify that a tour, encoded as a characteristic vector of its
#' edges, is indeed a feasible solution for the above program. The reverse is also
#' true.
#'
#' It is obvious that the above formulation is computationally critical for larger
#' graphs due to its sheer size. The number of subtour elimination constraints
#' grows exponentially with the number $n$ of nodes in the graph, which becomes
#' intractable even for moderate graph sizes.
#'
#' # Creating the Formulation Explicitly
#'
#' We generate $n = 10$  points $(x,y)$. The locations are drawn uniformly at
#' random from [-3,3]. We create a distance data frame that holds the Euclidean
#' distance for each pair of nodes Vi, Vj, $i < j$
#' Each row of this data frame will correspond to 1 binary variable of our TSP.
#'
#' ## Data Generation
## ----Data Generation----------------------------------------------------------
# the number of nodes
n <- 10

# random seed for reproducibility
set.seed(3701)
loc.matrix <- matrix(runif(n * 2, -3, 3), ncol = 2)

#
# each row of dist.df will correspond to 1 binary variable of our TSP
dist.df <- dist(loc.matrix) %>%
as.matrix() %>%
melt(varnames=c("Vi","Vj"), value.name="Dist") %>%
dplyr::filter(Vi < Vj)

tibble::glimpse(dist.df)

# plot the n nodes in the plane and label them
ggplot(data.frame(loc.matrix),
aes(x=loc.matrix[,1], y=loc.matrix[,2], label=1:n)
) +
geom_point() +
geom_label(nudge_y=0.3) +
labs(x="X", y="Y")

#'
#' ## Loading Node Degree Constraints
#'
#' We declare a function load_tsp that receives the above distance data frame as
#' input and creates a new Xpress Problem with only the node degree constraints
#' loaded. We align the columns we create with the dist.df data frame.
#'
#'
## ----Load TSP with Node Degree Constraints------------------------------------

# create a TSP problem with node degree constraints
load_tsp <- function(dist.df) {
p <- createprob()
nedges <-  nrow(dist.df)
problemdata <- list()

# column information
problemdata$lb <- rep(0,nedges) problemdata$ub <- rep(1,nedges)

# column names
problemdata$colname <- sprintf("x(%d,%d)",dist.df$Vi, dist.df$Vj) # objective coefficients problemdata$objcoef <- dist.df$Dist problemdata$columntypes <- rep("B", nedges)

# right hand side and row type
problemdata$rowtype <- rep("E", n) problemdata$rhs <- rep(2, n)

# we input the node degree constraints as sparse matrix in column major format
# directly using the corresponding problemdata names,
# see the documentation of
# xprs_loadproblemdata for more information

# each edge appears in exactly two node degree rows, once for each endpoint
# the number of elements in each column
problemdata$collen <- rep(2,nedges) # 1 additional marker for the end # the start index of the coefficients of each column. we need to pass a 0-based index into Xpress! problemdata$start <- 2 * (0:(nedges - 1)) # 0-based index!

# matrix coefficients, all 1's in the node degree constraints
problemdata$rowcoef <- rep(1, 2 * nedges) # row indices of each coefficient. We first create an array with 0's problemdata$rowind <- rep(0, 2 * nedges)

# The first endpoint Vi corresponds to the first row index of each edge
problemdata$rowind[2 * (1:nedges) - 1] <- dist.df$Vi - 1 # 0-based index!

# The second endpoint Vj corresponds to the second row index of each edge
problemdata$rowind[2 * (1:nedges)] <- dist.df$Vj - 1 # 0-based index!

# assign row names to identify the constraints easily in MPS or LP format
problemdata$probname <- "TSP" problemdata$rowname <- sprintf("NodeDeg_%d", 1:n)

# this invisibly returns the XPRSprob p from the function
return(xprs_loadproblemdata(p, problemdata))
}

#'
#' ## Solving the Node Degree TSP Without Subtour Elimination
#'
#' We solve this initial formulation using Xpress and plot the corresponding
#' solution.
#'
## ----Load and Solve with Node Degree Constraints------------------------------
p <- load_tsp(dist.df)
xprs_optimize(p)

#'
#' Let's plot the TSP solution. For this, we translate
#' the vertex locations $(x,y)$ into
#' edges (segments in ggplot2 lingo) with start and endpoints.
#'
## ----Plot the Solution--------------------------------------------------------
plot_tsp_solution <- function(p) {

# convert (x,y) locations of the nodes into a data frame.
loc.df <- as.data.frame(loc.matrix)
colnames(loc.df) <- c("X", "Y")
loc.df$V <- 1:n solution <- xprs_getsolution(p) # copy dist.df. For plotting, we must store the (x,y)-coordinates # of both endpoints plot.df <- dist.df[solution == 1,] plot.df <- plot.df %>% dplyr::mutate( Xi = loc.df$X[Vi],
Yi = loc.df$Y[Vi], Xj = loc.df$X[Vj],
Yj = loc.df$Y[Vj] ) plot.df %>% ggplot(aes(x=Xi, xend=Xj, y=Yi, yend=Yj, label=Vi)) + geom_segment() + geom_point() + geom_label(nudge_y = 0.3) + geom_label(aes(x=Xj, y=Yj, label=Vj), nudge_y = 0.3) + geom_point(aes(x=Xj, y=Yj)) } summary(p) plot_tsp_solution(p) #' #' ## Adding Subtour Elimination Constraints #' #' In a first shot, we add all subtour elimination constraints #' explicitly to the model. We can do this because the #' example is relatively small. #' A row always enters the problem as last row (at position xpress:::ROWS - 1). #' After adding a row, we declare it immediately as a delayed row. #' #' ## ----Adding Subtour Elimnation Constraints------------------------------------ add_subtour_elim <- function(p, dist.df) { nedges <- nrow(dist.df) # we enumerate all subsets that contain the first node, # i.e, all subsets of v_2, ... ,v_n with at most n - 3 elements. for (r in 1:(n - 3)) { allsubsets <- subsets(n-1, r, 2:n) # each row in allsubsets represents 1 subset for (s in 1:nrow(allsubsets)) { viins <- (dist.df$Vi == 1) | (dist.df$Vi %in% allsubsets[s,]) vjins <- (dist.df$Vj == 1) | (dist.df$Vj %in% allsubsets[s,]) edgeins <- (viins != vjins) # exactly one endpoint in S # use addrows to add each row explicitly addrows(p, rowtype = "G", rhs = 2, start = c(0,(r+1) * (n - r - 1)), colind = (1:nedges)[edgeins] - 1, rowcoef = rep(1, (r+1) * (n - r - 1)) ) # declare this row as delayed row. loaddelayedrows(p, rowind=getintattrib(p, xpress:::ROWS) - 1) } } # return p again to the caller p } #' #' ## Solving the Full TSP #' #' When we solve the full TSP with the added subtour elimination constraints, we #' obtain a correct tour, albeit with higher cost than the subtour solution from #' the previous run. #' ## ----------------------------------------------------------------------------- p <- load_tsp(dist.df) %>% add_subtour_elim(dist.df) %>% xprs_optimize() print(p) summary(p) plot_tsp_solution(p) #' #' # Separate Violated Subtour Elimination Constraints. #' #' It is much more efficient to separate subtour elimination constraints #' dynamically during the search. #' We initially omit all subtour elimination constraints. #' If we encounter an integer feasible LP solution, we verify if #' the solution is indeed a tour with full length. If not, #' we determine the subtour in which the$V1$lies. #' This subtour cannot reach all other nodes in the graph. #' #' We define a function get_subtour. This will be used #' for the separation of violated subtour elimination constraints. #' It will also be used to check the feasibility of any #' integer feasible solution. #' #' ## Callback Definition #' ## ----------------------------------------------------------------------------- get_subtour <- function(prob, sol) { # subset the distance data frame to the currently used edges dist.df.sol <- dist.df[sol >= 1 - 1e-6,] subtourindices <- NULL # initialize a logical to indicate which nodes we have visited. visited <- logical(n) # edge rep maps each node to its 2 adjacent edges in the solution edgerep <- matrix(0, ncol = 2, nrow=n) for (i in 1:nrow(dist.df.sol)) { vi <- dist.df.sol$Vi[i]
vj <- dist.df.sol$Vj[i] edgerep[vi, 1 + (edgerep[vi, 1] != 0)] <- i edgerep[vj, 1 + (edgerep[vj, 1] != 0)] <- i } # iterate through the edges of the solution. curredge <- 1 pathlength <- 0 currnode <- dist.df.sol$Vi[curredge]
while (!visited[currnode] && pathlength <= n) {
# label node as visited
visited[currnode] <- TRUE
# go to the other node of the current edge
currnode <- ifelse(dist.df.sol$Vi[curredge] == currnode, dist.df.sol$Vj[curredge], dist.df.sol$Vi[curredge]) #find the other edge to walk to edges <- edgerep[currnode,] curredge <- ifelse(edges[1] == curredge, edges[2], edges[1]) pathlength <- pathlength + 1 } # if the path length was smaller than n, we found a subtour that is not connected # to the rest of the graph. The subset is part of the visited array if (pathlength < n) { subtourindices <- which(visited) # indices of the corresponding nodes. } return(subtourindices) } #' ## ----Optnode Callback Definition---------------------------------------------- optnodecallback <- function(prob) { # return if there are fractional variables in the LP solution if (getintattrib(prob, xpress:::MIPINFEAS) > 0) return() currlpsol <- getlpsol(prob, x=TRUE)$x

# query the indices of a subtour
subtourindices <- get_subtour(prob, currlpsol)

if (!is.null(subtourindices)) {
viins <- (dist.df$Vi %in% subtourindices) vjins <- (dist.df$Vj %in% subtourindices)
edgeins <- (viins != vjins) # exactly one endpoint in S
subtourlength <- length(subtourindices)

presrow <- presolverow(p, rowtype = "G",
origcolind = which(edgeins) - 1, #
origrowcoef = rep(1, subtourlength * (n - subtourlength)),
origrhs = 2)

if (presrow$status == 0) { addcuts(p, cuttype=1, rowtype = "G", rhs = presrow$rhs,
start = c(0, presrow$ncoefs), colind = presrow$colind,
cutcoef = presrow$rowcoef) cat(paste("Added cut for subtour", paste(subtourindices, collapse = ", "), "\n")) } } # we always return infeasible = 0, # we cannot decide infeasibility at this point. return(list(infeasible = 0)) } preintsolcallback <- function(prob, soltype, cutoff) { # use getlpsol to access the current solution # pass in FALSE for optional arrays that need not be returned for our purpose sol <- getlpsol(prob, slack = F, duals = F, djs = F)$x

subtourindices <- get_subtour(prob, sol)
reject <- !is.null(subtourindices)

return (list(reject=reject))
}

#'
#'
#' ## Callback registration
#'
#' We register the two callbacks above.
#' We load the initial TSP only with node degree constraints.
#' We display the optimal solution to verify that
#' the callbacks worked correctly.
#'
## -----------------------------------------------------------------------------
p <- load_tsp(dist.df) # do not preload the subtour elimination constraints!
p <- p %>%
setintcontrol(xpress:::MIPDUALREDUCTIONS, 0) %>%
addcboptnode(optnodecallback) %>%
addcbpreintsol(preintsolcallback) %>%
setoutput()

xprs_optimize(p)
summary(p)

#'
#' The resulting solution represents a tour!
#'
## -----------------------------------------------------------------------------
plot_tsp_solution(p)

#'
#'
#'
#'
#'
#'
#'
#'