{{ message }}

Instantly share code, notes, and snippets.

# dirkschumacher/subtour_tsp_r.md

Created Sep 15, 2019
```# 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
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)

to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.