Optimization using Mathematical Programming and the GNU Linear Programming Kit
Michael Hahsler
Thu Dec 1 13:55:43 2016
This work is licensed under the Creative Commons Attribution 4.0 International License. For questions please contact Michael Hahsler.
This example shows how to solve a Traveling Salesman Problem (TSP).
library("Rglpk")## Loading required package: slam## Using the GLPK callable library version 4.57library("seriation") # for pimage
library("TSP") # for TSP supportRead program which creates objectives, constraints and bounds
The program with the data can be found here tsp.mod
x <- Rglpk_read_file("tsp.mod", type = "MathProg")
x## A mixed integer linear program with 480 objective variables,
## 240 are integer and 240 of which are binary variables.
## This problem has 288 constraints with 1440 non-zero values in the constraint matrix.str(x)## List of 5
## $ objective :List of 6
## ..$ i : int [1:240] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ j : int [1:240] 1 1 1 1 1 1 1 1 1 1 ...
## ..$ v : num [1:240] 509 501 312 1019 736 ...
## ..$ nrow : int 480
## ..$ ncol : int 1
## ..$ dimnames: NULL
## ..- attr(*, "class")= chr "simple_triplet_matrix"
## $ constraints:List of 3
## ..$ :List of 6
## .. ..$ i : int [1:1440] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ j : int [1:1440] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ v : num [1:1440] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ nrow : int 288
## .. ..$ ncol : int 480
## .. ..$ dimnames: NULL
## .. ..- attr(*, "class")= chr "simple_triplet_matrix"
## ..$ : chr [1:288] "==" "==" "==" "==" ...
## ..$ : num [1:288] 1 1 1 1 1 1 1 1 1 1 ...
## $ bounds :List of 2
## ..$ lower:List of 2
## .. ..$ ind: int [1:480] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ val: num [1:480] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ upper:List of 2
## .. ..$ ind: int [1:480] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ val: num [1:480] 1 1 1 1 1 1 1 1 1 1 ...
## $ types : chr [1:480] "B" "B" "B" "B" ...
## $ maximum : logi FALSE
## - attr(*, "class")= chr [1:2] "MP_data_from_file" "MILP"
## - attr(*, "n_objective_vars")= int 480
## - attr(*, "n_integer_vars")= int 240
## - attr(*, "n_binary_vars")= int 240
## - attr(*, "n_constraints")= int 288
## - attr(*, "n_nonzeros")= int 1440
## - attr(*, "problem_name")= chr "tsp"
## - attr(*, "objective_name")= chr "total"
## - attr(*, "objective_vars_names")= chr [1:480] "x[1,2]" "x[1,3]" "x[1,4]" "x[1,5]" ...
## - attr(*, "constraint_names")= chr [1:288] "leave[1]" "leave[2]" "leave[3]" "leave[4]" ...
## - attr(*, "file_type")= chr "MathProg"
## - attr(*, "file_name")= chr "/home/hahsler/hahsler_net/public_html/michael/SMU/EMIS7332/R/optimization/tsp.mod"The program contained the distance data in $objective in the same order as the decision variables x.
ns <- attr(x, "objective_vars_names")
ns <- do.call("rbind", strsplit(ns, "\\[|\\]|,"))
take <- ns[,1]=="x"
ns <- ns[take, 2:3]
mode(ns) <- "integer"
n <- max(ns)
d <- matrix(0L, nrow = n, ncol = n)
for(i in 1:nrow(ns)) d[ns[i,1], ns[i,2]] <- x$objective$v[i]
isSymmetric(d)## [1] TRUEd <- as.dist(d)
d## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 2 509
## 3 501 126
## 4 312 474 541
## 5 1019 1526 1516 1157
## 6 736 1226 1184 980 478
## 7 656 1133 1084 919 583 115
## 8 60 532 536 271 996 740 667
## 9 1039 1449 1371 1333 858 470 455 1066
## 10 726 1122 1045 1029 855 379 288 759 328
## 11 2314 2789 2728 2553 1504 1581 1661 2320 1387 1697
## 12 479 958 913 751 677 271 177 493 591 333 1838
## 13 448 941 904 704 651 289 216 454 650 400 1868 68
## 14 479 978 946 720 600 261 207 479 656 427 1841 105 52
## 15 619 1127 1115 783 401 308 343 598 776 622 1789 336 287 237
## 16 150 542 499 455 1033 687 592 206 933 610 2248 417 406 449
## 15
## 2
## 3
## 4
## 5
## 6
## 7
## 8
## 9
## 10
## 11
## 12
## 13
## 14
## 15
## 16 636pimage(d)
The distances in the example are calculated from 2d data using Euclidean distance, so we can reconstruct the relative city positions using MDS.
p <- cmdscale(d)
plot(p)
Solve the MILP
res <- Rglpk_solve_LP(x$objective, x$constraints[[1]], x$constraints[[2]],
x$constraints[[3]], x$bounds, x$types, x$maximum)
res## $optimum
## [1] 6859
##
## $solution
## [1] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [24] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [47] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [70] 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [93] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [116] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## [139] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [162] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## [185] 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [208] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [231] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 15 0 0 0 0 0 0
## [254] 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [277] 0 0 0 0 0 0 0 0 11 0 13 0 0 0 0 0 0 0 0 0 0 0 0
## [300] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 4 0
## [323] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0
## [346] 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [369] 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0
## [392] 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [415] 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## [438] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0
## [461] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 0 0 0 0 0
##
## $status
## [1] 0
##
## $solution_dual
## [1] NA
##
## $auxiliary
## $auxiliary$primal
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [18] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
## [35] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 0 0
## [52] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [69] 0 0 0 0 0 0 0 0 -4 0 -2 0 0 0 0 0 0
## [86] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [103] 0 0 0 -9 0 0 0 0 0 0 -11 0 0 0 0 0 0
## [120] 0 0 0 0 0 0 0 0 0 0 0 0 0 -12 0 0 0
## [137] 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0
## [154] 0 0 0 0 0 0 0 0 -7 0 0 0 0 0 0 0 0
## [171] 0 0 0 0 0 -6 0 0 0 0 0 0 0 0 0 0 -8
## [188] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [205] 0 0 0 0 -13 0 0 0 0 0 0 0 0 0 0 0 0
## [222] 0 0 0 -14 0 0 -15 0 0 0 0 0 0 0 0 0 0
## [239] 0 0 0 0 0 0 0 0 0 -10 0 0 0 0 0 0 0
## [256] 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0
## [273] -15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $auxiliary$dual
## [1] NAInterpret the results
We need to match the decision variables back to a tour. We can use ns from before to convert the solution into an adjacency matrix. We only need the solution for the xs and not the ys.
rs <- res$solution[take]
m <- matrix(0, nrow = n, ncol = n)
for(i in 1:nrow(ns)) m[ns[i,1], ns[i,2]] <- rs[i]
pimage(m)
convert matrix into a tour
o <- integer(n)
o[1] <- 1L
for(i in 2:n) o[i] <- which(m[o[i-1],] >0)
o## [1] 1 8 4 2 3 16 10 9 11 5 15 6 7 12 13 14calculate tour length, look at the reordered distance matrix and the tour.
tour_length(TSP(d), o)## [1] 6859pimage(d, o)
plot(ETSP(p), o)
compare solution with a heuristic
o2 <- solve_TSP(TSP(d))
o2## object of class 'TOUR'
## result of method 'arbitrary_insertion+two_opt' for 16 cities
## tour length: 6941plot(ETSP(p), o2)