knitr::opts_chunk$set(echo = TRUE)
# Model reference https://www.unc.edu/~pataki/papers/teachtsp.pdf
library(ROI)
library(ROI.plugin.glpk)
library(magrittr)
set.seed(42)

n <- 50
max_x <- 500
max_y <- 500
cities <- data.frame(id = 1:n, x = runif(n, max = max_x), y = runif(n, max = max_y))
distance <- as.matrix(dist(dplyr::select(cities, x, y), diag = TRUE, upper = TRUE))

rmpk <- function(solver = rmpk::ROI_optimizer("glpk")) {
  withr::with_package("rmpk", {
  model <- rmpk::optimization_model(solver)

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

  # a helper variable for the MTZ formulation of the tsp
  u <- model$add_variable("u", i = 1:n, lb = 1, ub = n)

  # 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
  # bounds not yet supported
  model$set_bounds(x[i, i], ub = 0, i = 1:n)

  # leave each city
  model$add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n)

  # visit each city
  model$add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 2:n)

  # ensure no subtours (arc constraints)
  model$add_constraint(u[i] >= 2, i = 2:n)
  model$add_constraint(u[i] - u[j] + 1 <= (n - 1) * (1 - x[i, j]), i = 2:n, j = 2:n)
  })
}

ompr_model <- function() {
  withr::with_package("ompr", {
    MIPModel() %>%
    add_variable(x[i,j], i = 1:n, j = 1:n,
               type = "integer", lb = 0, ub = 1) %>%
    add_variable(u[i], i = 1:n, lb = 1, ub = n) %>%
    set_objective(sum_expr(distance[i, j] * x[i, j], i = 1:n, j = 1:n), "min") %>%
    set_bounds(x[i, i], ub = 0, i = 1:n) %>%
    add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%
    add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 2:n) %>%
    add_constraint(u[i] >= 2, i = 2:n) %>%
    add_constraint(u[i] - u[j] + 1 <= (n - 1) * (1 - x[i, j]), i = 2:n, j = 2:n)
  })
}

bench::mark(rmpk(), ompr_model(), check = FALSE)
system.time(rmpk(solver = rmpk.glpk::GLPK()))


dirkschumacher/rmpk documentation built on Dec. 14, 2021, 5:13 p.m.