Lars Relund NielsenLars Relund Nielsen
Department of Economics and Business Economics
Fuglesangs Allé 4
DK-8210 Aarhus V
Office M330
Phone: +45 871 65145/+45 61 300 299 (mobile)
E-mail: lars@relund.dk

Cand.scient.oecon. (M.Sc.), Ph.D. 2004 in Mathematics and Economics at the University of Aarhus. Currently Associate Professor at Department of Economics and Business Economics, Business and Social Sciences, Aarhus University.

Using LocalSolver to solve VRP problems

R, Teaching 17. May 2017

The last days I have been playing a bit with LocalSolver and tried to model various vehicle routing problems (VRP). LocalSolver is a heuristic solver based on a heuristic search approach combining different optimization techniques. It includes a math modeling language (LSP) and there is a package so that it can be used from R. You need to download and install LocalSolver (I used v7.0) and get an academic license from their webpage.

Below I will formulate models for different VRPs and solve them using LocalSolver. The purpose of these tests is not to compare them with an optimal solution to the problems, but to examine the modelling capabilities of LSP. All the code used can be found in files simulate.R and ReadMe.R at GitHub.

Capacitated VRP (CVRP)

Let us start by consider the classical CVRP. First, we simulate some data:

source("simulate.R")
data <- simulateCVRP(customers = 10, seed = 578)
str(data)
#> List of 6
#>  $ nbTrucks      : int 5
#>  $ truckCapacity : int 10000
#>  $ nbNodes       : int 11
#>  $ nbCustomers   : int 10
#>  $ demands       : int [1:10] 2514 2585 2938 3935 2649 1429 3766 3706 2778 2547
#>  $ distanceMatrix: int [1:11, 1:11] 0 144 5 31 14 97 167 161 68 289 ...

Note that all data are integers since I have chosen to use integers in LocalSolver (input data must be the same data type). Next, we formulate the model using LSP. Here I used the example on LocalSolver’s webpage as a starting point. The learning curve for LSP is quite steep; however, the examples help. I was missing a forum where to ask question though. An important issue is that you must choose where the index in arrays start. I choose C style and let all index start from zero. We need decision variables routeSequences[k] containing a list of numbers in the range zero to customers-1. That is, if routeSequences[k] equals {0,4,6}, then truck number k visit customer 1, 5 and 7 (remember index start from zero).

library(localsolver)
model <- "function model(){
   routeSequences[k in 0..nbTrucks-1] <- list(nbCustomers);      // decision variable
   constraint partition[k in 0..nbTrucks-1](routeSequences[k]);  // visited by exaclty one route

   trucksUsed[k in 0..nbTrucks-1] <- count(routeSequences[k]) > 0;  
   nbTrucksUsed <- sum[k in 0..nbTrucks-1](trucksUsed[k]);

   for [k in 0..nbTrucks-1] {
      local sequence <- routeSequences[k];
      local c <- count(sequence);

      // The quantity needed in each route must not exceed the truck capacity
      routeQuantity[k] <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity[k] <= truckCapacity;

      // Distance travelled by truck k
      routeDistances[k] <- sum(1..c-1, i => distanceMatrix[sequence[i-1]+1][sequence[i]+1]) +
         (c > 0 ? (distanceMatrix[0][sequence[0]+1] + distanceMatrix[sequence[c-1]+1][0]) : 0);
   }

   // Total distance travelled
   totalDistance <- sum[k in 0..nbTrucks-1](routeDistances[k]);

   // Objective: minimize the number of trucks used, then minimize the distance travelled
   minimize nbTrucksUsed;   
   minimize totalDistance;
}"
lsp <- ls.problem(model)
lsp <- set.params(lsp, lsTimeLimit = c(20, 20), indexFromZero = TRUE, lsNbThreads = 4)

The model is added to LocalSolver using the localsolver package. I set the time limit to 20 seconds for each objective. Output returned to R can be added using function add.output.expr.

lsp <- set.temp.dir(lsp, path = getwd())  # add temporay files to this folder
lsp <- add.output.expr(lsp, "routeSequences", c(data$nbTrucks, data$nbCustomers))
lsp <- add.output.expr(lsp, "routeDistances", data$nbTrucks)
lsp <- add.output.expr(lsp, "routeQuantity", data$nbTrucks)
lsp <- add.output.expr(lsp, "nbTrucksUsed")
lsp <- add.output.expr(lsp, "totalDistance")
sol <- ls.solve(lsp, data)

The results are now in the list sol. Three temporary files are used by the solver. input.lsp contains the model (with input and output) passed to the solver, output.txt the solver log and data.txt the input data. Unfortunately, the solver status is not returned! So you may have to check the log to see if a feasible solution was found:

length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] TRUE

We postprocess the solution to get a result in a more readable format:

printSolution <- function(sol) {
    cat("Number of routes:", sol$nbTrucksUsed, "\n")
    cat("Total distance:", sol$totalDistance, "\n")
    cat("Routes:\n")
    for (r in 1:length(sol$routeDistances)) {
        if (max(sol$routeSequences[r, ]) >= 0) {
            cat("Customers:", paste0(sol$routeSequences[r, which(sol$routeSequences[r, ] != -1)], collapse = "-"), 
                "distance:", sol$routeDistances[r], "quantity:", sol$routeQuantity[r], "\n")
        }
    }
}
printSolution(sol)
#> Number of routes: 3 
#> Total distance: 1147 
#> Routes:
#> Customers: 4-1-3 distance: 222 quantity: 9169 
#> Customers: 7-6-0 distance: 559 quantity: 9986 
#> Customers: 2-5-8-9 distance: 366 quantity: 9692

Extension – Workday constraints

Assume that only a single driver is available and that the maximum possible workday is specified. We need workday length, travel times and handling time. Note that we now use maxRoutes instead of trucks which denote the maximum number of routes per driver.

data <- simulateDriver(customers = 10, seed = 578, drivers = 1, maxRoutes = 5)
str(data)
#> List of 11
#>  $ truckCapacity     : int 10000
#>  $ nbNodes           : int 11
#>  $ nbCustomers       : int 10
#>  $ demands           : int [1:10] 2514 2585 2938 3935 2649 1429 3766 3706 2778 2547
#>  $ distanceMatrix    : int [1:11, 1:11] 0 144 5 31 14 97 167 161 68 289 ...
#>  $ maxRoutes         : int 5
#>  $ handlingTime      : int 30
#>  $ drivers           : int 1
#>  $ depotEarliestStart: int 360
#>  $ depotLatestArrival: int 1200
#>  $ travelTimeMatrix  : int [1:11, 1:11] 0 14 28 29 6 26 24 29 13 11 ...

To handle travel time constraints we add variables routeStops[k][i] which is the node visited at stop i (depot = node 0) and similar arrivalTimes[k][i] and departureTimes[k][i].

model <- "function model(){
   routeSequences[k in 0..maxRoutes-1] <- list(nbCustomers);  
   routeStops[k in 0..maxRoutes-1][0] = 0;
   routeStops[k in 0..maxRoutes-1][nbNodes] = 0;
   departureTimes[k in 0..maxRoutes-1][i in 0..nbNodes] <- int(0,depotLatestArrival);
   arrivalTimes[k in 0..maxRoutes-1][i in 0..nbNodes] <- int(0,depotLatestArrival);

   constraint partition[k in 0..maxRoutes-1](routeSequences[k]);  

   for [k in 0..maxRoutes-1] {
      local sequence <- routeSequences[k];
      local c <- count(sequence);
      routesUsed[k] <- c > 0;

      routeQuantity[k] <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity[k] <= truckCapacity;

      routeDistances[k] <- sum(1..c-1, i => distanceMatrix[sequence[i-1]+1][sequence[i]+1]) +
         (c > 0 ? (distanceMatrix[0][sequence[0]+1] + distanceMatrix[sequence[c-1]+1][0]) : 0);

      if (k==0) constraint departureTimes[k][0] >= depotEarliestStart;  // first route 
      for [i in 1..nbNodes] {
         routeStops[k][i] <- routeSequences[k][i-1] + 1;
         arrivalTimes[k][i] <- departureTimes[k][i-1] + travelTimeMatrix[routeStops[k][i-1]][routeStops[k][i]];
         departureTimes[k][i] <- (routeStops[k][i] > 0  ? arrivalTimes[k][i] + handlingTime : 0);
      }
      routeArrivalDepot[k] <- departureTimes[k][0] + (c > 0 ? (travelTimeMatrix[0][sequence[0]+1] +
         travelTimeMatrix[sequence[c-1]+1][0]) : 0) +
         sum(1..c-1, i => travelTimeMatrix[sequence[i-1]+1][sequence[i]+1]) + c*handlingTime;

      if (k>0) {
         constraint departureTimes[k][0] >= routeArrivalDepot[k-1];
         constraint routeArrivalDepot[k] <= depotLatestArrival;
      } else {
         constraint routeArrivalDepot[k] <= depotLatestArrival;
      }
   }

   totalRoutesUsed <- sum[k in 0..maxRoutes-1](routesUsed[k]);
   totalDistance <- sum[k in 0..maxRoutes-1](routeDistances[k]);
   latestArrivalTime <- max[k in 0..maxRoutes-1](routeArrivalDepot[k]);

   minimize totalDistance;
   minimize latestArrivalTime;
}"
lsp <- ls.problem(model)
lsp <- set.params(lsp, lsTimeLimit = 20, indexFromZero = TRUE, lsNbThreads = 4)
lsp <- set.temp.dir(lsp, path = getwd())
lsp <- add.output.expr(lsp, "routeSequences", c(data$maxRoutes, data$nbCustomers))
lsp <- add.output.expr(lsp, "routeDistances", data$maxRoutes)
lsp <- add.output.expr(lsp, "routeQuantity", data$maxRoutes)
lsp <- add.output.expr(lsp, "totalRoutesUsed")
lsp <- add.output.expr(lsp, "totalDistance")
lsp <- add.output.expr(lsp, "latestArrivalTime")
lsp <- add.output.expr(lsp, "routeStops", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "arrivalTimes", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "departureTimes", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "routeArrivalDepot", data$maxRoutes)
sol <- ls.solve(lsp, data)
length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] TRUE

We have minimized the total distance and next the lastest arrival time at the depot and postprocess the solution to get the result in a more readable format:

printSolution <- function(sol) {
    cat("Number of routes:", sol$totalRoutesUsed, "\n")
    cat("Total distance:", sol$totalDistance, "\n")
    cat("Latest depot arrival (min from midnight): ", sol$latestArrivalTime, " (constraint <=", data$depotLatestArrival, 
        ")\n", sep = "")
    cat("Routes:\n")
    for (r in 1:length(sol$routeArrivalDepot)) {
        if (max(sol$routeStops[r, ]) > 0) {
            cat("Customers:", paste0(sol$routeStops[r, 1:(which(sol$routeStops[r, ] == 0)[2])], collapse = "-"), 
                "departure:", sol$departureTimes[r, 1], "finished:", sol$arrivalTimes[r, which(sol$routeStops[r, 
                  ] == 0)[2]], "distance:", sol$routeDistances[r], "quantity:", sol$routeQuantity[r], 
                "\n")
            cat("  Details:", sol$routeStops[r, 1], paste0("(d:", sol$departureTimes[r, 1], ") ->"))
            for (i in 2:(which(sol$routeStops[r, ] == 0)[2] - 1)) {
                cat("\n           ", sol$routeStops[r, i], " (tt:", data$travelTimeMatrix[sol$routeStops[r, 
                  i - 1] + 1, sol$routeStops[r, i] + 1], " a:", sol$arrivalTimes[r, i], " d:", sol$departureTimes[r, 
                  i], ") ", "->", sep = "")
            }
            cat("\n           ", sol$routeStops[r, which(sol$routeStops[r, ] == 0)[2]], " (tt:", data$travelTimeMatrix[sol$routeStops[r, 
                i - 1] + 1, sol$routeStops[r, i] + 1], " a:", sol$arrivalTimes[r, which(sol$routeStops[r, 
                ] == 0)[2]], ")\n", sep = "")
        }
    }
}
printSolution(sol)
#> Number of routes: 4 
#> Total distance: 795 
#> Latest depot arrival (min from midnight): 908 (constraint <=1200)
#> Routes:
#> Customers: 0-8-0 departure: 360 finished: 416 distance: 136 quantity: 3706 
#>   Details: 0 (d:360) ->
#>            8 (tt:13 a:373 d:403) ->
#>            0 (tt:13 a:416)
#> Customers: 0-1-10-7-0 departure: 416 finished: 572 distance: 377 quantity: 8827 
#>   Details: 0 (d:416) ->
#>            1 (tt:14 a:430 d:460) ->
#>            10 (tt:6 a:466 d:496) ->
#>            7 (tt:17 a:513 d:543) ->
#>            0 (tt:17 a:572)
#> Customers: 0-4-2-0 departure: 572 finished: 675 distance: 47 quantity: 6520 
#>   Details: 0 (d:572) ->
#>            4 (tt:6 a:578 d:608) ->
#>            2 (tt:9 a:617 d:647) ->
#>            0 (tt:9 a:675)
#> Customers: 0-5-6-9-3-0 departure: 675 finished: 908 distance: 235 quantity: 9794 
#>   Details: 0 (d:675) ->
#>            5 (tt:26 a:701 d:731) ->
#>            6 (tt:9 a:740 d:770) ->
#>            9 (tt:22 a:792 d:822) ->
#>            3 (tt:27 a:849 d:879) ->
#>            0 (tt:27 a:908)

Note currently there is no time added for filling the truck at the depot. Let us try another instance.

data <- simulateDriver(customers = 20, seed = 578, drivers = 1, maxRoutes = 12)
lsp <- clear.output.exprs(lsp)
lsp <- add.output.expr(lsp, "routeSequences", c(data$maxRoutes, data$nbCustomers))
lsp <- add.output.expr(lsp, "routeDistances", data$maxRoutes)
lsp <- add.output.expr(lsp, "routeQuantity", data$maxRoutes)
lsp <- add.output.expr(lsp, "totalRoutesUsed")
lsp <- add.output.expr(lsp, "totalDistance")
lsp <- add.output.expr(lsp, "latestArrivalTime")
lsp <- add.output.expr(lsp, "routeStops", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "arrivalTimes", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "departureTimes", c(data$maxRoutes, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "routeArrivalDepot", data$maxRoutes)
sol <- ls.solve(lsp, data)
length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] FALSE
printSolution(sol)
#> Number of routes: 7 
#> Total distance: 3682 
#> Latest depot arrival (min from midnight): 1200 (constraint <=1200)
#> Routes:
#> Customers: 0-6-14-16-5-0 departure: 360 finished: 514 distance: 708 quantity: 9133 
#>   Details: 0 (d:360) ->
#>            6 (tt:5 a:365 d:395) ->
#>            14 (tt:5 a:400 d:430) ->
#>            16 (tt:11 a:441 d:471) ->
#>            5 (tt:5 a:476 d:506) ->
#>            0 (tt:5 a:514)
#> Customers: 0-3-19-9-0 departure: 514 finished: 644 distance: 463 quantity: 7914 
#>   Details: 0 (d:514) ->
#>            3 (tt:6 a:520 d:550) ->
#>            19 (tt:7 a:557 d:587) ->
#>            9 (tt:10 a:597 d:627) ->
#>            0 (tt:10 a:644)
#> Customers: 0-1-20-11-0 departure: 644 finished: 780 distance: 545 quantity: 9844 
#>   Details: 0 (d:644) ->
#>            1 (tt:11 a:655 d:685) ->
#>            20 (tt:10 a:695 d:725) ->
#>            11 (tt:12 a:737 d:767) ->
#>            0 (tt:12 a:780)
#> Customers: 0-7-17-0 departure: 780 finished: 869 distance: 437 quantity: 7642 
#>   Details: 0 (d:780) ->
#>            7 (tt:15 a:795 d:825) ->
#>            17 (tt:9 a:834 d:864) ->
#>            0 (tt:9 a:869)
#> Customers: 0-4-12-8-0 departure: 869 finished: 987 distance: 542 quantity: 8861 
#>   Details: 0 (d:869) ->
#>            4 (tt:5 a:874 d:904) ->
#>            12 (tt:11 a:915 d:945) ->
#>            8 (tt:6 a:951 d:981) ->
#>            0 (tt:6 a:987)
#> Customers: 0-2-10-18-0 departure: 973 finished: 1123 distance: 536 quantity: 9062 
#>   Details: 0 (d:973) ->
#>            2 (tt:20 a:993 d:1023) ->
#>            10 (tt:9 a:1032 d:1062) ->
#>            18 (tt:7 a:1069 d:1099) ->
#>            0 (tt:7 a:1123)
#> Customers: 0-13-15-0 departure: 1123 finished: 1200 distance: 451 quantity: 5671 
#>   Details: 0 (d:1123) ->
#>            13 (tt:6 a:1129 d:1159) ->
#>            15 (tt:6 a:1165 d:1195) ->
#>            0 (tt:6 a:1200)

Note no feasible solution can be found, since the driver can not keep the latest arrival at depot limit. That is, the printed solution is the unfeasible solution found when the solver stopped. This is not the same as a proof of that no feasible solution exists since it is a heuristic solver we use.

Extension – Multiple drivers

Assume that multiple drivers are allowed and for simplicity assume that they have the same working time. To handle multiple drivers we need to add a driver index to the model.

data$drivers <- as.integer(3)
model <- "function model(){
   routeSequences[k in 0..maxRoutes-1][d in 0..drivers-1] <- list(nbCustomers);  
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][0] = 0;
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][nbNodes] = 0;
   departureTimes[k in 0..maxRoutes-1][d in 0..drivers-1][i in 0..nbNodes] <- int(-1,depotLatestArrival);
   arrivalTimes[k in 0..maxRoutes-1][d in 0..drivers-1][i in 0..nbNodes] <- int(-1,depotLatestArrival);

   constraint partition[k in 0..maxRoutes-1][d in 0..drivers-1](routeSequences[k][d]);

   for [k in 0..maxRoutes-1][d in 0..drivers-1] {
      local sequence <- routeSequences[k][d];
      local c <- count(sequence);
      routesUsed[k][d] <- c > 0;

      routeQuantity[k][d] <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity[k][d] <= truckCapacity;

      routeDistances[k][d] <- (c > 0 ? (distanceMatrix[0][sequence[0]+1] + 
         distanceMatrix[sequence[c-1]+1][0]) : 0) +
         sum(1..c-1, i => distanceMatrix[sequence[i-1]+1][sequence[i]+1]);

      if (k==0) constraint departureTimes[k][d][0] >= depotEarliestStart;
      for [i in 1..nbNodes] {
         routeStops[k][d][i] <- sequence[i-1] + 1;
         arrivalTimes[k][d][i] <- departureTimes[k][d][i-1] + 
         travelTimeMatrix[routeStops[k][d][i-1]][routeStops[k][d][i]];
         departureTimes[k][d][i] <- (routeStops[k][d][i] > 0  ? arrivalTimes[k][d][i] + handlingTime : -1);
      }
      routeArrivalDepot[k][d] <- departureTimes[k][d][0] + (c > 0 ? (travelTimeMatrix[0][sequence[0]+1] + travelTimeMatrix[sequence[c-1]+1][0]) : 0) +
         sum(1..c-1, i => travelTimeMatrix[sequence[i-1]+1][sequence[i]+1]) + c*handlingTime;

      if (k>0) {
         constraint departureTimes[k][d][0] >= routeArrivalDepot[k-1][d];
         constraint routeArrivalDepot[k][d] <= depotLatestArrival;
      } else {
         constraint routeArrivalDepot[k][d] <= depotLatestArrival;
      }
   }

   totalRoutesUsed <- sum[k in 0..maxRoutes-1][d in 0..drivers-1] (routesUsed[k][d]);
   totalDistance <- sum[k in 0..maxRoutes-1][d in 0..drivers-1](routeDistances[k][d]);
   latestArrivalTime <- max[k in 0..maxRoutes-1][d in 0..drivers-1] (routeArrivalDepot[k][d]);

   minimize totalDistance;
   minimize latestArrivalTime;
}"
lsp <- ls.problem(model)
lsp <- set.params(lsp, lsTimeLimit = c(20, 20), indexFromZero = TRUE, lsNbThreads = 4)
lsp <- set.temp.dir(lsp, path = getwd())
lsp <- add.output.expr(lsp, "routeSequences", c(data$maxRoutes, data$drivers, data$nbCustomers))
lsp <- add.output.expr(lsp, "routeDistances", c(data$maxRoutes, data$drivers))
lsp <- add.output.expr(lsp, "routeQuantity", c(data$maxRoutes, data$drivers))
lsp <- add.output.expr(lsp, "totalRoutesUsed")
lsp <- add.output.expr(lsp, "totalDistance")
lsp <- add.output.expr(lsp, "latestArrivalTime")
lsp <- add.output.expr(lsp, "routeStops", c(data$maxRoutes, data$drivers, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "arrivalTimes", c(data$maxRoutes, data$drivers, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "departureTimes", c(data$maxRoutes, data$drivers, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "routeArrivalDepot", c(data$maxRoutes, data$drivers))
sol <- ls.solve(lsp, data)
length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] TRUE

We considered the same instance as before but now with 3 drivers, minimized the total distance and next the latest arrival time at the depot. The solution are:

printSolution <- function(sol) {
    cat("Number of routes:", sol$totalRoutesUsed, "\n")
    cat("Total distance:", sol$totalDistance, "\n")
    cat("Latest depot arrival (min from midnight):", sol$latestArrivalTime, "\n")
    cat("Routes:\n")
    for (d in 1:length(sol$routeArrivalDepot[1, ])) {
        for (r in 1:length(sol$routeArrivalDepot[, 1])) {
            if (max(sol$routeStops[r, d, ]) > 0) {
                cat("Driver:", d, "Customers:", paste0(sol$routeStops[r, d, 1:(which(sol$routeStops[r, 
                  d, ] == 0)[2])], collapse = "-"), "departure:", sol$departureTimes[r, d, 1], "finished:", 
                  sol$arrivalTimes[r, d, which(sol$routeStops[r, d, ] == 0)[2]], "distance:", sol$routeDistances[r], 
                  "quantity:", sol$routeQuantity[r], "\n")
                cat("  Details:", sol$routeStops[r, d, 1], paste0("(d:", sol$departureTimes[r, d, 1], 
                  ") ->"))
                for (i in 2:(which(sol$routeStops[r, d, ] == 0)[2] - 1)) {
                  cat("\n           ", sol$routeStops[r, d, i], " (tt:", data$travelTimeMatrix[sol$routeStops[r, 
                    d, i - 1] + 1, sol$routeStops[r, d, i] + 1], " a:", sol$arrivalTimes[r, d, i], " d:", 
                    sol$departureTimes[r, d, i], ") ", "->", sep = "")
                }
                cat("\n           ", sol$routeStops[r, d, which(sol$routeStops[r, d, ] == 0)[2]], " (tt:", 
                  data$travelTimeMatrix[sol$routeStops[r, d, i - 1] + 1, sol$routeStops[r, d, i] + 1], 
                  " a:", sol$arrivalTimes[r, d, which(sol$routeStops[r, d, ] == 0)[2]], ")\n", sep = "")
            }
        }
    }
}
printSolution(sol)
#> Number of routes: 7 
#> Total distance: 1476 
#> Latest depot arrival (min from midnight): 728 
#> Routes:
#> Driver: 1 Customers: 0-8-18-6-0 departure: 364 finished: 501 distance: 153 quantity: 9065 
#>   Details: 0 (d:364) ->
#>            8 (tt:6 a:370 d:400) ->
#>            18 (tt:15 a:415 d:445) ->
#>            6 (tt:21 a:466 d:496) ->
#>            0 (tt:21 a:501)
#> Driver: 1 Customers: 0-12-9-14-13-0 departure: 519 finished: 719 distance: 319 quantity: 8426 
#>   Details: 0 (d:519) ->
#>            12 (tt:15 a:534 d:564) ->
#>            9 (tt:30 a:594 d:624) ->
#>            14 (tt:19 a:643 d:673) ->
#>            13 (tt:10 a:683 d:713) ->
#>            0 (tt:10 a:719)
#> Driver: 2 Customers: 0-3-15-10-0 departure: 379 finished: 549 distance: 0 quantity: 0 
#>   Details: 0 (d:379) ->
#>            3 (tt:6 a:385 d:415) ->
#>            15 (tt:20 a:435 d:465) ->
#>            10 (tt:28 a:493 d:523) ->
#>            0 (tt:28 a:549)
#> Driver: 2 Customers: 0-2-5-4-0 departure: 575 finished: 711 distance: 0 quantity: 0 
#>   Details: 0 (d:575) ->
#>            2 (tt:20 a:595 d:625) ->
#>            5 (tt:5 a:630 d:660) ->
#>            4 (tt:16 a:676 d:706) ->
#>            0 (tt:16 a:711)
#> Driver: 3 Customers: 0-7-0 departure: 360 finished: 420 distance: 0 quantity: 0 
#>   Details: 0 (d:360) ->
#>            7 (tt:15 a:375 d:405) ->
#>            0 (tt:15 a:420)
#> Driver: 3 Customers: 0-19-11-20-0 departure: 420 finished: 577 distance: 0 quantity: 0 
#>   Details: 0 (d:420) ->
#>            19 (tt:20 a:440 d:470) ->
#>            11 (tt:7 a:477 d:507) ->
#>            20 (tt:12 a:519 d:549) ->
#>            0 (tt:12 a:577)
#> Driver: 3 Customers: 0-1-17-16-0 departure: 577 finished: 728 distance: 0 quantity: 0 
#>   Details: 0 (d:577) ->
#>            1 (tt:11 a:588 d:618) ->
#>            17 (tt:14 a:632 d:662) ->
#>            16 (tt:17 a:679 d:709) ->
#>            0 (tt:17 a:728)

Note with 3 drivers we can keep the latest arrival at depot limit.

Extension – Delivery over multiple days

Assume that the orders can be delivered at specific days (we ignore the time windows for the moment):

data <- simulateVRPTW(drivers = 2, maxRoutes = 5, customers = 30, days = 7, seed = 789)
str(data)
#> List of 15
#>  $ truckCapacity     : int 10000
#>  $ nbNodes           : int 31
#>  $ nbCustomers       : int 30
#>  $ demands           : int [1:30] 3100 1280 1035 2775 2476 1060 2718 1497 2077 2039 ...
#>  $ distanceMatrix    : int [1:31, 1:31] 0 212 1 214 236 78 226 240 300 86 ...
#>  $ maxRoutes         : int 5
#>  $ handlingTime      : int 30
#>  $ drivers           : int 2
#>  $ depotEarliestStart: int 360
#>  $ depotLatestArrival: int 1200
#>  $ travelTimeMatrix  : int [1:31, 1:31] 0 10 27 19 20 11 9 29 29 11 ...
#>  $ days              : int 7
#>  $ validDays         : int [1:30, 1:7] 0 1 1 1 0 1 1 0 0 0 ...
#>  $ timeWindowsLB     : int [1:31] 360 900 600 900 300 300 900 600 900 600 ...
#>  $ timeWindowsUB     : int [1:31] 1200 1020 660 1020 420 360 1020 720 1020 720 ...
data$validDays
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#>  [1,]    0    1    1    1    1    1    1
#>  [2,]    1    1    1    1    1    1    0
#>  [3,]    1    1    1    1    1    1    1
#>  [4,]    1    1    1    1    1    0    0
#>  [5,]    0    0    0    0    0    1    0
#>  [6,]    1    1    1    1    1    1    1
#>  [7,]    1    1    1    1    1    1    1
#>  [8,]    0    1    1    1    1    1    0
#>  [9,]    0    0    1    1    0    0    0
#> [10,]    0    0    0    0    1    1    0
#> [11,]    1    1    1    1    1    1    1
#> [12,]    0    1    1    1    1    1    1
#> [13,]    1    1    1    1    1    1    0
#> [14,]    1    1    1    1    1    1    1
#> [15,]    0    0    1    1    1    0    0
#> [16,]    0    0    0    1    1    1    1
#> [17,]    0    0    0    0    1    1    0
#> [18,]    0    1    1    1    0    0    0
#> [19,]    0    1    1    0    0    0    0
#> [20,]    0    0    1    1    0    0    0
#> [21,]    1    1    0    0    0    0    0
#> [22,]    0    0    0    1    0    0    0
#> [23,]    1    0    0    0    0    0    0
#> [24,]    1    1    1    1    1    1    1
#> [25,]    1    1    1    1    1    1    1
#> [26,]    0    0    1    1    1    1    1
#> [27,]    1    1    1    1    1    1    1
#> [28,]    0    1    1    1    1    1    0
#> [29,]    0    0    0    1    1    1    1
#> [30,]    0    1    1    1    0    0    0

Valid days for an order corresponds to a one in the matrix validDays. For instance, customer 29 can be visited in days 4, 5, 6, 7. We add a constraint that customer i must be visited at valid days.

model <- "
function model(){
   routeSequences[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] <- list(nbCustomers); 
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][0] = 0;
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][nbNodes] = 0;
   departureTimes[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][0] <- int(0,depotLatestArrival);
   arrivalTimes[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 0..nbNodes] <- int(0,depotLatestArrival);

   constraint partition[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeSequences[k][d][t]);

   for [k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] {
      local sequence <- routeSequences[k][d][t];
      local c <- count(sequence);
      routesUsed[k][d][t] <- c > 0;

      routeQuantity[k][d][t] <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity[k][d][t] <= truckCapacity;

      routeDistances[k][d][t] <-
         (c > 0 ? (distanceMatrix[0][sequence[0]+1] + distanceMatrix[sequence[c-1]+1][0]) : 0) +
         sum(1..c-1, i => distanceMatrix[sequence[i-1]+1][sequence[i]+1]);

      if (k==0) constraint departureTimes[k][d][t][0] >= depotEarliestStart;
      for [i in 1..nbNodes] {
         routeStops[k][d][t][i] <- sequence[i-1] + 1;
         arrivalTimes[k][d][t][i] <- (routeStops[k][d][t][i-1] > 0 || (i==1 && routeStops[k][d][t][i]>0) ?
            departureTimes[k][d][t][i-1] +
            travelTimeMatrix[routeStops[k][d][t][i-1]][routeStops[k][d][t][i]] : 0);
         departureTimes[k][d][t][i] <-
            (routeStops[k][d][t][i] > 0  ? arrivalTimes[k][d][t][i] + handlingTime : 0);
      }
      routeArrivalDepot[k][d][t] <- max[i in 1..nbCustomers](arrivalTimes[k][d][t][i]);

      if (k>0) {
         constraint departureTimes[k][d][t][0] >= routeArrivalDepot[k-1][d][t];
         constraint routeArrivalDepot[k][d][t] <= depotLatestArrival;
      } else {
         constraint routeArrivalDepot[k][d][t] <= depotLatestArrival;
      }
   }

   for [t in 0..days-1][i in 0..nbCustomers-1] {
      visitAtDay[t][i] <- sum[k in 0..maxRoutes-1][d in 0..drivers-1](indexOf(routeSequences[k][d][t],i) > -1);
      constraint visitAtDay[t][i] <= validDays[i][t];
   }

   totalRoutesUsed <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] (routesUsed[k][d][t]);
   totalDistance <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeDistances[k][d][t]);
   latestArrival <- max[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeArrivalDepot[k][d][t]);

   minimize totalDistance;
   minimize latestArrival;

   // modify results so can output to R (cannot output array with more than 3 index)
   for [k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] {
         local r = k + d*maxRoutes + t*maxRoutes*drivers;
         routes[r][i in 0..nbNodes] <- routeStops[k][d][t][i];
         arrivalT[r][i in 0..nbNodes] <- arrivalTimes[k][d][t][i];
         routeDepotArrival[r] <- routeArrivalDepot[k][d][t];
         departureT[r][i in 0..nbNodes] <- departureTimes[k][d][t][i];
         routeDist[r] <- routeDistances[k][d][t];
         routeQuant[r] <- routeQuantity[k][d][t];
         driver[r] <- d+1;
         day[r] <- t+1;
   }
}"
lsp <- ls.problem(model)
lsp <- set.params(lsp, lsTimeLimit = c(10, 10), indexFromZero = TRUE, lsNbThreads = 4)
lsp <- set.temp.dir(lsp, path = getwd())
lsp <- add.output.expr(lsp, "totalRoutesUsed")
lsp <- add.output.expr(lsp, "totalDistance")
lsp <- add.output.expr(lsp, "latestArrival")
maxR <- data$maxRoutes * data$days * data$drivers
lsp <- add.output.expr(lsp, "routes", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "arrivalT", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "departureT", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "driver", maxR)
lsp <- add.output.expr(lsp, "day", maxR)
lsp <- add.output.expr(lsp, "routeDist", maxR)
lsp <- add.output.expr(lsp, "routeDepotArrival", maxR)
lsp <- add.output.expr(lsp, "routeQuant", maxR)
sol <- ls.solve(lsp, data)
length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] TRUE
printSolution <- function(sol) {
    cat("Number of routes:", sol$totalRoutesUsed, "\n")
    cat("Total distance:", sol$totalDistance, "\n")
    cat("Latest depot arrival (min from midnight):", sol$latestArrival, "\n")
    cat("Routes:\n")
    for (r in 1:length(sol$driver)) {
        if (max(sol$routes[r, ]) > 0) {
            cat("Day:", sol$day[r], "driver:", sol$driver[r], "customers:", paste0(sol$routes[r, 1:(which(sol$routes[r, 
                ] == 0)[2])], collapse = "-"), "departure:", sol$departureT[r, 1], "finished:", sol$arrivalT[r, 
                which(sol$routes[r, ] == 0)[2]], "distance:", sol$routeDist[r], "quantity:", sol$routeQuant[r], 
                "\n")
            cat("  Details:", sol$routes[r, 1], paste0("(d:", sol$departureT[r, 1], ") ->"))
            for (i in 2:(which(sol$routes[r, ] == 0)[2] - 1)) {
                cat("\n           ", sol$routes[r, i], " (tt:", data$travelTimeMatrix[sol$routes[r, i - 
                  1] + 1, sol$routes[r, i] + 1], " a:", sol$arrivalT[r, i], " d:", sol$departureT[r, 
                  i], ") ", "->", sep = "")
            }
            cat("\n           ", sol$routes[r, which(sol$routes[r, ] == 0)[2]], " (tt:", data$travelTimeMatrix[sol$routes[r, 
                i - 1] + 1, sol$routes[r, i] + 1], " a:", sol$arrivalT[r, which(sol$routes[r, ] == 0)[2]], 
                ")\n", sep = "")
        }
    }
}
printSolution(sol)
#> Number of routes: 7 
#> Total distance: 1778 
#> Latest depot arrival (min from midnight): 290 
#> Routes:
#> Day: 1 driver: 1 customers: 0-25-4-23-6-21-0 departure: 0 finished: 290 distance: 246 quantity: 9628 
#>   Details: 0 (d:0) ->
#>            25 (tt:19 a:19 d:49) ->
#>            4 (tt:30 a:79 d:109) ->
#>            23 (tt:29 a:138 d:168) ->
#>            6 (tt:28 a:196 d:226) ->
#>            21 (tt:20 a:246 d:276) ->
#>            0 (tt:20 a:290)
#> Day: 3 driver: 1 customers: 0-2-1-28-26-0 departure: 20 finished: 271 distance: 68 quantity: 9421 
#>   Details: 0 (d:20) ->
#>            2 (tt:27 a:47 d:77) ->
#>            1 (tt:26 a:103 d:133) ->
#>            28 (tt:26 a:159 d:189) ->
#>            26 (tt:25 a:214 d:244) ->
#>            0 (tt:25 a:271)
#> Day: 3 driver: 2 customers: 0-9-19-20-18-0 departure: 13 finished: 222 distance: 319 quantity: 9657 
#>   Details: 0 (d:13) ->
#>            9 (tt:11 a:24 d:54) ->
#>            19 (tt:28 a:82 d:112) ->
#>            20 (tt:21 a:133 d:163) ->
#>            18 (tt:13 a:176 d:206) ->
#>            0 (tt:13 a:222)
#> Day: 4 driver: 1 customers: 0-12-8-30-16-0 departure: 89 finished: 290 distance: 237 quantity: 7582 
#>   Details: 0 (d:89) ->
#>            12 (tt:17 a:106 d:136) ->
#>            8 (tt:23 a:159 d:189) ->
#>            30 (tt:30 a:219 d:249) ->
#>            16 (tt:5 a:254 d:284) ->
#>            0 (tt:5 a:290)
#> Day: 4 driver: 2 customers: 0-14-15-13-22-0 departure: 71 finished: 286 distance: 135 quantity: 8790 
#>   Details: 0 (d:71) ->
#>            14 (tt:18 a:89 d:119) ->
#>            15 (tt:24 a:143 d:173) ->
#>            13 (tt:20 a:193 d:223) ->
#>            22 (tt:8 a:231 d:261) ->
#>            0 (tt:8 a:286)
#> Day: 5 driver: 2 customers: 0-27-24-11-7-17-0 departure: 10 finished: 290 distance: 446 quantity: 9975 
#>   Details: 0 (d:10) ->
#>            27 (tt:20 a:30 d:60) ->
#>            24 (tt:27 a:87 d:117) ->
#>            11 (tt:21 a:138 d:168) ->
#>            7 (tt:24 a:192 d:222) ->
#>            17 (tt:30 a:252 d:282) ->
#>            0 (tt:30 a:290)
#> Day: 6 driver: 2 customers: 0-5-3-29-10-0 departure: 1 finished: 213 distance: 327 quantity: 9020 
#>   Details: 0 (d:1) ->
#>            5 (tt:11 a:12 d:42) ->
#>            3 (tt:20 a:62 d:92) ->
#>            29 (tt:9 a:101 d:131) ->
#>            10 (tt:27 a:158 d:188) ->
#>            0 (tt:27 a:213)

Note there may not always be deliveries on a given day.

Time windows

Finally, let us try to add time windows to the model above. I here formulate the time windows as soft constraints and try to minimize the total time not satisfying the time window.

model <- "
function model(){
   routeSequences[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] <- list(nbCustomers); 
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][0] = 0;
   routeStops[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][nbNodes] = 0;
   departureTimes[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][0] <- int(0,depotLatestArrival);
   arrivalTimes[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 0..nbNodes] <- int(0,depotLatestArrival);
   waitingTimes[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 0..nbNodes] <- int(0,600);
   constraint partition[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeSequences[k][d][t]);

   for [k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] {
      local sequence <- routeSequences[k][d][t];
      local c <- count(sequence);
      routesUsed[k][d][t] <- c > 0;

      routeQuantity[k][d][t] <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity[k][d][t] <= truckCapacity;

      routeDistances[k][d][t] <-
         (c > 0 ? (distanceMatrix[0][sequence[0]+1] + distanceMatrix[sequence[c-1]+1][0]) : 0) +
         sum(1..c-1, i => distanceMatrix[sequence[i-1]+1][sequence[i]+1]);

      if (k==0) constraint departureTimes[k][d][t][0] >= depotEarliestStart;
      for [i in 1..nbNodes] {
         routeStops[k][d][t][i] <- sequence[i-1] + 1;
         arrivalTimes[k][d][t][i] <- (routeStops[k][d][t][i-1] > 0 || (i==1 && routeStops[k][d][t][i]>0) ?
            departureTimes[k][d][t][i-1] +
            travelTimeMatrix[routeStops[k][d][t][i-1]][routeStops[k][d][t][i]] +
            waitingTimes[k][d][t][i] : 0);
         // If hard time window constraints
         //constraint arrivalTimes[k][d][t][i] >= (routeStops[k][d][t][i] > 0) * timeWindowsLB[routeStops[k][d][t][i]];
         //constraint arrivalTimes[k][d][t][i] <= timeWindowsUB[routeStops[k][d][t][i]];
         departureTimes[k][d][t][i] <-
            (routeStops[k][d][t][i] > 0  ? arrivalTimes[k][d][t][i] + handlingTime : 0);
      }
      routeArrivalDepot[k][d][t] <- max[i in 1..nbCustomers](arrivalTimes[k][d][t][i]);

      if (k>0) {
         constraint departureTimes[k][d][t][0] >= routeArrivalDepot[k-1][d][t];
         constraint routeArrivalDepot[k][d][t] <= depotLatestArrival;
      } else {
         constraint routeArrivalDepot[k][d][t] <= depotLatestArrival;
      }
   }

   for [t in 0..days-1][i in 0..nbCustomers-1] {
      visitAtDay[t][i] <- sum[k in 0..maxRoutes-1][d in 0..drivers-1](indexOf(routeSequences[k][d][t],i) > -1);
      constraint visitAtDay[t][i] <= validDays[i][t];
   }

   totalRoutesUsed <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] (routesUsed[k][d][t]);
   totalDistance <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeDistances[k][d][t]);
   totalWaiting <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 0..nbNodes]
      (waitingTimes[k][d][t][i]);
   totalViolateTWCtr <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 1..nbCustomers](
      routeStops[k][d][t][i] > 0 && (arrivalTimes[k][d][t][i] < timeWindowsLB[routeStops[k][d][t][i]]  ||
      arrivalTimes[k][d][t][i] > timeWindowsUB[routeStops[k][d][t][i]])
   );
   totalViolateTWMin <- sum[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1][i in 1..nbCustomers](
      (routeStops[k][d][t][i] > 0 && arrivalTimes[k][d][t][i] < timeWindowsLB[routeStops[k][d][t][i]]) *
      (timeWindowsLB[routeStops[k][d][t][i]] - arrivalTimes[k][d][t][i]) +
      (routeStops[k][d][t][i] > 0 && arrivalTimes[k][d][t][i] > timeWindowsUB[routeStops[k][d][t][i]]) *
      (arrivalTimes[k][d][t][i] - timeWindowsUB[routeStops[k][d][t][i]])
   );
   latestArrival <- max[k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1](routeArrivalDepot[k][d][t]);

   minimize totalViolateTWMin;   // total waiting time in minutes
   minimize totalViolateTWCtr;   // times we wait
   minimize totalRoutesUsed;
   minimize totalDistance;
   minimize totalWaiting;
   minimize latestArrival;

   // modify results so can output to R (cannot output array with more than 3 index)
   for [k in 0..maxRoutes-1][d in 0..drivers-1][t in 0..days-1] {
         local r = k + d*maxRoutes + t*maxRoutes*drivers;
         routes[r][i in 0..nbNodes] <- routeStops[k][d][t][i];
         arrivalT[r][i in 0..nbNodes] <- arrivalTimes[k][d][t][i];
         routeDepotArrival[r] <- routeArrivalDepot[k][d][t];
         departureT[r][i in 0..nbNodes] <- departureTimes[k][d][t][i];
         waitingT[r][i in 0..nbNodes] <- waitingTimes[k][d][t][i];
         routeDist[r] <- routeDistances[k][d][t];
         routeQuant[r] <- routeQuantity[k][d][t];
         driver[r] <- d+1;
         day[r] <- t+1;
   }
}"
rm(lsp)
lsp <- ls.problem(model)
lsp <- set.params(lsp, lsTimeLimit = c(20, 20, 20, 20, 20, 20), indexFromZero = TRUE, lsNbThreads = 4)
lsp <- set.temp.dir(lsp, path = getwd())
lsp <- add.output.expr(lsp, "totalRoutesUsed")
lsp <- add.output.expr(lsp, "totalDistance")
lsp <- add.output.expr(lsp, "totalWaiting")
lsp <- add.output.expr(lsp, "totalViolateTWCtr")
lsp <- add.output.expr(lsp, "totalViolateTWMin")
lsp <- add.output.expr(lsp, "latestArrival")
maxR <- data$maxRoutes * data$days * data$drivers
lsp <- add.output.expr(lsp, "routes", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "arrivalT", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "departureT", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "waitingT", c(maxR, data$nbNodes + 1))
lsp <- add.output.expr(lsp, "driver", maxR)
lsp <- add.output.expr(lsp, "day", maxR)
lsp <- add.output.expr(lsp, "routeDist", maxR)
lsp <- add.output.expr(lsp, "routeDepotArrival", maxR)
lsp <- add.output.expr(lsp, "routeQuant", maxR)
sol <- ls.solve(lsp, data)
length(grep("Feasible solution", readLines("output.txt"), value = TRUE, ignore.case = FALSE)) > 0
#> [1] TRUE
printSolution <- function(sol) {
    cat("Number of routes:", sol$totalRoutesUsed, "\n")
    cat("Total distance:", sol$totalDistance, "\n")
    cat("Total waiting time:", sol$totalWaiting, "\n")
    cat("Total time window violations:", sol$totalViolateTWCtr, "\n")
    cat("Total time window violations (min):", sol$totalViolateTWMin, "\n")
    cat("Latest depot arrival (min from midnight):", sol$latestArrival, "\n")
    cat("Routes:\n")
    for (r in 1:length(sol$driver)) {
        if (max(sol$routes[r, ]) > 0) {
            cat("Day:", sol$day[r], "driver:", sol$driver[r], "customers:", paste0(sol$routes[r, 1:(which(sol$routes[r, 
                ] == 0)[2])], collapse = "-"), "departure:", sol$departureT[r, 1], "finished:", sol$arrivalT[r, 
                which(sol$routes[r, ] == 0)[2]], "distance:", sol$routeDist[r], "quantity:", sol$routeQuant[r], 
                "\n")
            cat("  Details:", sol$routes[r, 1], paste0("(d:", sol$departureT[r, 1], ") ->"))
            for (i in 2:(which(sol$routes[r, ] == 0)[2] - 1)) {
                cat("\n           ", sol$routes[r, i], " (tt:", data$travelTimeMatrix[sol$routes[r, i - 
                  1] + 1, sol$routes[r, i] + 1], " w:", sol$waitingT[r, i], " a:", sol$arrivalT[r, i], 
                  " d:", sol$departureT[r, i], ") ", " tw:[", data$timeWindowsLB[sol$routes[r, i] + 1], 
                  ",", data$timeWindowsUB[sol$routes[r, i] + 1], "] ->", sep = "")
            }
            cat("\n           ", sol$routes[r, which(sol$routes[r, ] == 0)[2]], " (tt:", data$travelTimeMatrix[sol$routes[r, 
                i - 1] + 1, sol$routes[r, i] + 1], " a:", sol$arrivalT[r, which(sol$routes[r, ] == 0)[2]], 
                ")\n", sep = "")
        }
    }
}
printSolution(sol)
#> Number of routes: 10 
#> Total distance: 3192 
#> Total waiting time: 3764 
#> Total time window violations: 0 
#> Total time window violations (min): 0 
#> Latest depot arrival (min from midnight): 1015 
#> Routes:
#> Day: 1 driver: 2 customers: 0-23-24-6-21-0 departure: 638 finished: 994 distance: 304 quantity: 6474 
#>   Details: 0 (d:638) ->
#>            23 (tt:22 w:0 a:660 d:690)  tw:[600,660] ->
#>            24 (tt:8 w:0 a:698 d:728)  tw:[600,720] ->
#>            6 (tt:16 w:156 a:900 d:930)  tw:[900,1020] ->
#>            21 (tt:20 w:0 a:950 d:980)  tw:[900,1020] ->
#>            0 (tt:20 a:994)
#> Day: 3 driver: 1 customers: 0-25-19-9-0 departure: 327 finished: 699 distance: 204 quantity: 7676 
#>   Details: 0 (d:327) ->
#>            25 (tt:19 w:74 a:420 d:450)  tw:[300,420] ->
#>            19 (tt:19 w:131 a:600 d:630)  tw:[600,720] ->
#>            9 (tt:28 w:0 a:658 d:688)  tw:[600,720] ->
#>            0 (tt:28 a:699)
#> Day: 3 driver: 1 customers: 0-20-18-12-0 departure: 392 finished: 1000 distance: 299 quantity: 6791 
#>   Details: 0 (d:392) ->
#>            20 (tt:28 w:0 a:420 d:450)  tw:[300,420] ->
#>            18 (tt:13 w:437 a:900 d:930)  tw:[900,1020] ->
#>            12 (tt:23 w:0 a:953 d:983)  tw:[900,960] ->
#>            0 (tt:23 a:1000)
#> Day: 4 driver: 1 customers: 0-15-13-0 departure: 303 finished: 418 distance: 254 quantity: 4050 
#>   Details: 0 (d:303) ->
#>            15 (tt:7 w:0 a:310 d:340)  tw:[300,360] ->
#>            13 (tt:20 w:0 a:360 d:390)  tw:[300,420] ->
#>            0 (tt:20 a:418)
#> Day: 4 driver: 1 customers: 0-30-3-22-0 departure: 420 finished: 992 distance: 345 quantity: 5885 
#>   Details: 0 (d:420) ->
#>            30 (tt:24 w:216 a:660 d:690)  tw:[600,660] ->
#>            3 (tt:17 w:193 a:900 d:930)  tw:[900,1020] ->
#>            22 (tt:7 w:0 a:937 d:967)  tw:[900,960] ->
#>            0 (tt:7 a:992)
#> Day: 4 driver: 2 customers: 0-26-28-1-8-0 departure: 244 finished: 1015 distance: 370 quantity: 9638 
#>   Details: 0 (d:244) ->
#>            26 (tt:27 w:389 a:660 d:690)  tw:[600,660] ->
#>            28 (tt:25 w:0 a:715 d:745)  tw:[600,720] ->
#>            1 (tt:26 w:129 a:900 d:930)  tw:[900,1020] ->
#>            8 (tt:26 w:0 a:956 d:986)  tw:[900,1020] ->
#>            0 (tt:26 a:1015)
#> Day: 5 driver: 1 customers: 0-2-7-29-10-0 departure: 125 finished: 1012 distance: 286 quantity: 9507 
#>   Details: 0 (d:125) ->
#>            2 (tt:27 w:508 a:660 d:690)  tw:[600,660] ->
#>            7 (tt:10 w:13 a:713 d:743)  tw:[600,720] ->
#>            29 (tt:29 w:128 a:900 d:930)  tw:[900,1020] ->
#>            10 (tt:27 w:0 a:957 d:987)  tw:[900,1020] ->
#>            0 (tt:27 a:1012)
#> Day: 5 driver: 2 customers: 0-16-4-17-0 departure: 314 finished: 447 distance: 462 quantity: 5298 
#>   Details: 0 (d:314) ->
#>            16 (tt:6 w:0 a:320 d:350)  tw:[300,360] ->
#>            4 (tt:8 w:0 a:358 d:388)  tw:[300,420] ->
#>            17 (tt:21 w:0 a:409 d:439)  tw:[300,420] ->
#>            0 (tt:21 a:447)
#> Day: 6 driver: 2 customers: 0-5-14-0 departure: 102 finished: 948 distance: 197 quantity: 4874 
#>   Details: 0 (d:102) ->
#>            5 (tt:11 w:247 a:360 d:390)  tw:[300,360] ->
#>            14 (tt:8 w:502 a:900 d:930)  tw:[900,960] ->
#>            0 (tt:8 a:948)
#> Day: 7 driver: 2 customers: 0-11-27-0 departure: 200 finished: 950 distance: 471 quantity: 3880 
#>   Details: 0 (d:200) ->
#>            11 (tt:15 w:145 a:360 d:390)  tw:[300,360] ->
#>            27 (tt:14 w:496 a:900 d:930)  tw:[900,1020] ->
#>            0 (tt:14 a:950)

All time windows has been satisfied.

Pitfalls

To summarize some of the pitfalls I experienced during coding:

  • All input data must have be of the same data type as specified in the model (e.g. integers).
  • Remember when index from zero, this holds for all variables/data.
  • In set.params the parameter lsTimeLimit=c(10,10) must have the same length as the number of objectives. If e.g. lsTimeLimit=10 it corresponds to lsTimeLimit=c(0,10) (only the last objective is optimized) and not lsTimeLimit=c(10,10) which seems more obvious.
  • You must specify the dimensions of output. Hence, when you run a new instance, you have to add output expressions again.
  • The dimension of output arrays cannot be more than 3 for some strange reason. So you may have to transform your output.

Plotting IP and MO-IP models in R

LaTeX, R, Teaching 20. February 2017

I recent released a small package gMOIP which can make 2D plots of linear and integer programming models (LP/IP). With the package you can make plots of the feasible region (or solution space) of an LP, visualize the integer points inside the region and the iso profit curve. Moreover, can also make a plot of a bi-objective criterion space and the non-dominated (Pareto) set. Figures are prepared for LaTeX and can automatically be transformed to TikZ using tikzDevice.

You can install the latest stable release from CRAN:

install.packages("gMOIP")

Alternatively, install the latest development version from GitHub:

install.packages("devtools")
devtools::install_github("relund/gMOIP")

Let us have a look at some examples. First we specify the feasible region and calculate the corner points of the polytope and the integer points inside the feasible region.

# Feasible region Ax<=b, x>=0
A <- matrix(c(9,10,2,4,-3,2), ncol = 2, byrow = TRUE)
b <- c(90,27,3)
# Corner points of the polytope
cPoints<-cornerPoints(A, b)
# Integer points in the polytope
iPoints<-integerPoints(A, b)

It is now easy to make different plots of the polytope using plotPolytope:

# Plot of the polytope
p1 <- plotPolytope(cPoints) + ggtitle("Feasible region only")
p2 <- plotPolytope(cPoints, cPoints, iso = c(7.75, 10), crit = "max") + 
   ggtitle("Solution LP (max)")
p3 <- plotPolytope(cPoints, iPoints, iso = c(7.75, 10), crit = "max") + 
   ggtitle("Solution IP (max)")
p4 <- plotPolytope(cPoints, iPoints, iso = c(3,-3), crit = "min") + 
   ggtitle("Solution IP (min)")
grid.arrange(p1, p2, p3, p4, nrow = 2)

plot of chunk plotPoly

You may also create a tikz file of the plot for LaTeX using:

library(tikzDevice)
tikz(file = "plot_polytope.tex", standAlone=F, width = 7, height = 6)
plotPolytope(cPoints, iPoints, iso = c(3,-3), crit = "min", latex = TRUE)
dev.off()

If you consider a bi-objective IP problem the solution space can be visualized using plotCriterion. First, consider model \(\min\{ \left( (-1,1)x, (1,-1)x\right) | Ax\leq b, x\geq 0, x \text{ integer} \}\). Note that solutions and criterion points are linked using labels.

zPoints<-criterionPoints(iPoints, c1 = c(-1, 1), c2 = c(1, -1), crit = "min")
plotPolytope(cPoints, zPoints, showLbl = TRUE) + ggtitle("Solution space")
plotCriterion(zPoints, addHull = TRUE, addTriangles = TRUE, crit = "min", showLbl = TRUE) + 
   ggtitle("Criterion space")

plot of chunk plotCrit1plot of chunk plotCrit1

Next, consider model \(\max\{ \left( (2,-1)x, (1,4)x\right) | Ax\leq b, x\geq 0, x \text{ integer} \}\). Here solutions corresponding to a non-dominated point is visualized in the solution space using a different shape.

zPoints<-criterionPoints(iPoints, c1 = c(-2, -1), c2 = c(1, 4), crit = "max")
plotPolytope(cPoints, zPoints, showLbl = TRUE, shape = zPoints$nD) + ggtitle("Solution space")
plotCriterion(zPoints, addHull = TRUE, addTriangles = TRUE, crit = "max", showLbl = TRUE) +
   ggtitle("Criterion space")

plot of chunk plotCrit2plot of chunk plotCrit2

Finally, consider model \(\min\{ \left( (2,-1)x, (1,4)x\right) | Ax\leq b, x\geq 0, x \text{ integer} \}\)

zPoints<-criterionPoints(iPoints, c1 = c(-2, -1), c2 = c(1, 4), crit = "min")
plotPolytope(cPoints, zPoints, showLbl = TRUE, shape = zPoints$nD) + ggtitle("Solution space") 
plotCriterion(zPoints, addHull = TRUE, addTriangles = TRUE, crit = "min", showLbl = TRUE) + 
   ggtitle("Criterion space")

plot of chunk plotCrit3plot of chunk plotCrit3

Animation of wind using Shiny

R 19. November 2015

I have played a bit with Shiny the last days. Let us try to create an shiny web application which show a map and add a layer with wind directions and wind speed. First we need weather data for a given latitude and longitude. We will use the API from forecast.io to retrieve data. Your will need an API key for using the service. We define functions to retrieve the data Læs mere »

Distance matrix calculations in R

R 7. July 2015

Many models in Operations Research needs a distance matrix between nodes in a network. If the distance matrix is based on road distances it can be found using R. Let us try to compute the distances between a set for zip codes in Denmark. First we load all zip codes for Jutland in Denmark: Læs mere »

Publishing R markdown to WordPress

R 31. March 2015

In some cases it may be useful to write a WordPress post in R Markdown and afterwards publish it to my blog. This can be done using the RWordPress package. First we setup a link to my blog:

if (!require('RWordPress'))
   install.packages('RWordPress', repos = 'http://www.omegahat.org/R', type = 'source')
library(RWordPress)
options(WordPressLogin = c(<your username> = '<your password>'),
        WordPressURL = 'http://www.research.relund.dk/wp/xmlrpc.php')

Next the post is written in a Rmd file and afterwards published to WordPress:

id<-knit2wp('RWordPress_post.Rmd', title = 'Publishing R markdown to WordPress', 
            categories = c('R'), publish=F )

Now the post with id is a draft on my blog and I can have a look at it before publishing it. If I want to update the post I do:

knit2wp('RWordPress_post.Rmd', title = 'Publishing R markdown to WordPress', postid=id, 
            action='editPost', categories = c('R'), publish = F )

To highlight the above code I use the WP Code Highlight.js plugin. If you use other syntax highlight plugins you may have to change the code a bit.

Migrating from svn (R-forge) to git (GitHub)

R 19. March 2015

I recent migrated my R package mdp from my R-forge svn repo to GitHub. Do the following:

  1. Create an empty folder mdp and open a shell and import the svn to git (I only imported the pkg sub-folder svn+ssh://relund@svn.r-forge.r-project.org/svnroot/mdp/pkg)
    git svn clone svn+ssh://relund@svn.r-forge.r-project.org/svnroot/mdp/pkg .     # import the svn
    git branch -a    # info, should show a git-svn remote branch
    git svn info     # show svn details (also URL)
  2. Now add your local repo to GitHub. First, create an empty GitHub repo and next run from the shell:
    git remote add origin https://github.com/relund/mdp.git
    git push -u origin master

    The first line tells Git that your local repo has a remote version on GitHub, and calls it “origin”. The second line pushes all your current work to that repo.

Now you have a local Git repo and 2 remote repos (one at GitHub and one at R-forge). You want to use GitHub as the main repo. That is, add files, make changes to the master branch and commit to the master branch and GitHub as usual. At some point in time you want to commit the current revision to the remote svn repo at R-forge (so that the package can be checked at different OS etc.):

git svn info               # show svn details (also URL)
git svn dcommit –-dry-run  # show which svn branch you will commit into:
git svn rebase             # pull changes from svn repository
git svn dcommit            # push your local git commits to svn repository

Resources

Using RStudio together with Git and GitHub on Windows

R 11. March 2015

I have started to use Git and GitHub together with RStudio. Git is a distributed version control system which is very useful when doing reproducible research. It is a good way to handle programming/coding. Moreover, Git (via GitHub) allows groups of people to work on the same documents (often code) at the same time, and without stepping on each other’s toes. RStudio is an excellent integrated development environment built specifically for R. Læs mere »

Converting a multidimensional index to a single unique id

R, Research 17. May 2010

I have written a small note about mapping a multidimensional index into a unique id. An unique id might be useful e.g. when having multiple state variables in an
Markov decision process. Here each state at a stage may be identified using the id.

R package testcpp now at R-Forge

R 29. September 2009

I recent wrote a post about testing external pointers with finalization in R using C++. The test package is now available at R-Forge for inspiration.

Accessing an array created in R from C++

R 28. August 2009

Assume that you have created a 3-dim array in R containing vectors of various size:

> a<-array(list(),c(2,2,2))   # array where each element contain a numeric vector
> a[[1,1,1]]<-rnorm(1)
> a[[1,1,2]]<-rnorm(3)
> a[[2,1,1]]<-rnorm(2)
> a[[2,1,2]]<-rnorm(3)
> a[[1,2,1]]<-rnorm(4)
> a[[1,2,2]]<-rnorm(3)
> a[[2,2,1]]<-rnorm(5)
> a[[2,2,2]]<-rnorm(1)
> d = dim(a)
> d
[1] 2 2 2
> i<-1; j<-1; k<-2
> q<-i + j*d[1] + k*(d[1]*d[2]) - (d[1] + d[1]*d[2])
> a[[i,j,k]]      # access an element
[1] -0.2370040 -0.6009635  3.0550405
> a[[q]]          # access the same element using a single index
[1] -0.2370040 -0.6009635  3.0550405

Note you can access the array in two ways. Next you want to copy the array to a vector on the C++ side. For instance in a package you may need to do some operations in C++ to speed up time.

On the C++ side you can create a function: Læs mere »

© Lars Relund Nielsen
Entries RSS Comments RSS Log in