subtour_tsp_r.md

# an example of the TSP solved through solver callbacks
# follows the formulation of the Gurobi example
# http://examples.gurobi.com/traveling-salesman-problem/
# and from the TSP vignette for the MTZ formulation

# all experimental

library(ggplot2)
suppressPackageStartupMessages(library(dplyr))
library(rmpk)
library(rmpk.glpk)

n <- 25

# from 0 to ...
max_x <- 500
max_y <- 500
set.seed(123456)
cities <- data.frame(id = 1:n, x = runif(n, max = max_x), y = runif(n, max = max_y))


ggplot(cities, aes(x, y)) + 
  geom_point()

distance <- as.matrix(dist(select(cities, x, y), diag = TRUE, upper = TRUE))


solver <- GLPK()
model <- MIPModel(solver)

# we create a variable that is 1 iff we travel from city i to j
model$add_variable(x[i, j], i = 1:n, j = 1:n, 
                   type = "binary", lb = 0, ub = 1)

# minimize travel distance
model$set_objective(sum_expr(distance[i, j] * x[i, j], i = 1:n, j = 1:n), "min")

# you cannot go to the same city
model$set_bounds(x[i, i], ub = 0, i = 1:n)

# the model is undirected
model$add_constraint(x[i, j] == x[j, i], i = 1:n, j = 1:n)

# if you enter a node you have to leave it as well twice
model$add_constraint(sum_expr(x[i, j], j = 1:n) == 2, i = 1:n)

solver$set_irowgen_callback(function() {
  # first we construct an adjacency matrix
  # we also stop if we find fractional values
  adjacency_matrix <- matrix(0L, nrow = n, ncol = n)
  all_integral <- TRUE
  for (i in 1:n) {
    for (j in 1:n) {
      adjacency_matrix[i, j] <- solver$glpk_get_col_prim(x[i, j])
      if (adjacency_matrix[i, j] %% 1 != 0) {
        all_integral <- FALSE
        break
      }
    }
    if (!all_integral) {
      break
    }
  }
  
  if (!all_integral) {
    return()
  }
  
  # find the shortest circle and check if that tour does not contain all the nodes
  # if yes, we add a subtour elimination constraint
  graph <- igraph::graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected")
  shortest_circle <- igraph::girth(graph, TRUE)
  circle_size <- shortest_circle$girth
  if (circle_size < n && circle_size > 1) {
    first_v <- shortest_circle$circle[[1L]]
    prev_v <- first_v
    expr <- 0
    for (v in shortest_circle$circle[2:circle_size]) {
      expr <- expr + x[prev_v, v]
      prev_v <- v
    }
    expr <- expr + x[prev_v, first_v]
    model$add_constraint(expr <= circle_size - 1)
  }
})

model$optimize()

tour <- model$get_variable_value(x[i,j])
obj <- model$objective_value() / 2

tour <- tour[tour$value == 1, ]
paths <- select(tour, i, j) %>% 
  rename(from = i, to = j) %>% 
  mutate(trip_id = row_number()) %>% 
  tidyr::gather(property, idx_val, from:to) %>% 
  mutate(idx_val = as.integer(idx_val)) %>% 
  inner_join(cities, by = c("idx_val" = "id"))

ggplot(cities, aes(x, y)) + 
  geom_point() + 
  geom_line(data = paths, aes(group = trip_id)) + 
  ggtitle(paste0("Optimal route with cost: ", round(obj, 2)))

^{Created on 2019-09-15 by the reprex package (v0.3.0)}