source("../R/scm_functions.R") # nearest_neighbor(), two_opt(), tour_length(), savings_cw()
library(ompr); library(ompr.roi); library(ROI.plugin.glpk); library(dplyr)
# --- TSPLIB fri26
D <- as.matrix(read.table("https://people.math.sc.edu/Burkardt/datasets/tsp/fri26_d.txt"))
dimnames(D) <- NULL
opt <- 937
nn <- sapply(1:26, function(s) tour_length(nearest_neighbor(D, start = s), D))
summary((nn - opt) / opt)
# --- Exakte Lösung (MTZ) für die ersten n Städte
tsp_mtz <- function(D) {
n <- nrow(D)
MIPModel() |>
add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") |>
add_variable(u[i], i = 1:n, lb = 1, ub = n) |>
set_objective(sum_expr(D[i, j] * x[i, j], i = 1:n, j = 1:n), "min") |>
add_constraint(x[i, i] == 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 = 1:n) |>
add_constraint(u[i] - u[j] + 1 <= (n - 1) * (1 - x[i, j]),
i = 2:n, j = 2:n, i != j) |>
solve_model(with_ROI(solver = "glpk"))
}
system.time(sol <- tsp_mtz(D[1:8, 1:8])); objective_value(sol)
# --- CVRPLIB A-n32-k5 (Datei vorher herunterladen)
lies_vrp <- function(datei) {
z <- trimws(readLines(datei))
cap <- as.numeric(sub(".*:", "", grep("^CAPACITY", z, value = TRUE)))
a <- grep("NODE_COORD_SECTION", z); b <- grep("DEMAND_SECTION", z)
e <- grep("DEPOT_SECTION", z)
xy <- read.table(text = z[(a + 1):(b - 1)])
dm <- read.table(text = z[(b + 1):(e - 1)])
list(cap = cap, xy = xy[, 2:3], q = dm[, 2])
}
inst <- lies_vrp("A-n32-k5.vrp")
Dv <- round(as.matrix(dist(inst$xy)))
ids <- c("Depot", paste0("K", 2:nrow(Dv)))
dimnames(Dv) <- list(ids, ids)
q <- setNames(inst$q[-1], ids[-1])
vrp <- savings_cw(Dv, q, cap = inst$cap, depot = "Depot")
length(vrp$touren)
sum(sapply(vrp$touren, function(t) tour_length(c("Depot", t), Dv)))