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

Markov decision theory

Teaching 28. March 2014

During the this semester I am teaching the course “Markov decision theory” at the Department of Mathematics. The course presents the algorithmic aspects of Markov decision theory and illustrates the wide applicability of this theory to a variety of realistic problems. First, the course consider the Poisson process followed by Renewal/Reward Processes. Next we study the theory of Markov chains. Then follows a thorough presentation of the theory of Markov Decision Processes, including some applications. Finally, further applications are discussed based on reports made by the students. A course plan is given here.

Course on supervision

Administration, Teaching 13. May 2013

During the past months I have been participating in a course about supervision. The course, which runs over several days, covers master, PhD and fellow supervision and use results from research, practice and experiences from the participants. It aims at strengthen the participants’ overall supervisory skills in the form of

  • clear and specific criteria-based guidance
  • better feedback when commenting on text from the student
  • establishment and maintenance of a mutually binding supervision agreement
  • initial, ongoing and final project evaluation and feedback
  • using meta-communication of supervision content and processes

The course included video-filming ourselves in doing supervision which was a good way of seeing how one actually perform. The course can be recommended to all university staff doing supervision.

Sets, sums and constraints – A suvival guide for Business School students

Teaching 30. August 2012

Often when I teach students at our Business School they have a hard time understanding compact linear programming (LP) formulations. So here it is,  a short introduction to some of the concepts you need to know for understanding compact LP formulations.

Sets

A set is a group of elements, e.g. $\{1,2,4\}$ is a set with 3 elements, namely, $1,2$ and $4$ and $A=\{(2,3),(4,5),(6,8),(5,6)\}$ is a set called $A$ with 4 elements (pairs), namely,  $(2,3),(4,5),(6,8)$ and $(5,6)$. Note that in the last case each element is a pair $(i,j)$. Sets containing pairs are often used when we formulate LPs based on network problems where the pair $(i,j)$ denote the arc/edge from node $i$ to node $j$. Læs mere »

Using logvrp.com when introducing students to VRP

Teaching 22. December 2011

In my recent course on transportation and distribution systems for undergraduate business students I used a week talking about vehicle routing problems (VRP). Since their math skills aren’t that high the focus was not on algorithms instead I tried to talk about the general structure of the problem. Here the site logvrp.com is very good and provide a way of making a realistic problem case.

Logvrp is a web based vehicle route planner software cable of minimize transportation costs and to minimize fleet mileages. It uses Google maps to display routes and trips. Moreover, it supports many different VRP problems. Current two heuristics are implemented and can be run simultaneously on a problem instance. This is a good thing since the students can see that the heuristics behave differently on the same instance and that there is no general winner.

I created an exercise about a small-sized company located in Aarhus which sell organic fruits and vegetables to companies and private persons. The procedure to create the data to be used at logvrp was as follows: first I created sets of customers using the logvrp site. Each customer group was saved as csv files and modified using R. Next, a sets of orders was generated randomly using R and saved as csv files. The students was now able to answer the different questions using the csv files provided.

routes_day1
A typical solution visualized using crow fly distances.

 

The solution shown using shortest paths.

Transportation and distribution systems

Teaching 3. December 2011

During the fall 2011 I have been teaching the course “Transportation and distribution systems” an optional course for undergraduate students. We have considered various topics within transportation such as its role in the global economy, modal choice, 3PL, VRP problems, flow problems, airline management, ITS, green logistics etc.

Each week the students was given an exercise where they in most cases had to formulate and solve an mathematical programming model. I also have tried to introduce CPLEX OPL studio to them with minor success. Mostly they like to formulate their models in Excel and with OpenSolver it is possible even for quite large models.

Teaching large classes

Teaching 23. February 2011

Wordle: UntitledYesterday I participated in a short course about how to teach large classes. Focus was on how to align teaching for constructive learning, collaborative learning techniques and other active learning tricks. The main point can be summarized in this statement

“If students are to learn desired outcomes in a reasonably effective manner, then the teacher’s fundamental task is to get students to engage in learning activities that are likely to result in their achieving those outcomes… It is helpful to remember that what the student does is actually more important in determining what is learned than what the teacher does.” (Shuell, 1986: 429)

Teaching Management Science

Teaching 15. December 2010

Just finished my last Management Science lecture for this semester. This is the first time I theach this introductory course to Business Students. Topics are linear and integer programming, goal programming, simulation and network models. Læs mere »

Seminar in Business Studies – First Experience

Teaching 4. May 2010

The from February to May this year I have been teaching in “Seminar in Business Studies” a bachelor course (HA – 4. semester).

The course represents a teaching and form of learning where a group of 12 students with one or more teachers discusses selected issues and problems in management science. First the student should write a 10 page report on a specific subject. The report should provide the possibility of in-depth studies of theories, models and methods in management science. While writing the report he can use the assigned advisor for at most 1.5 hours. Second the student, 2 weeks after delivering the report, should give a 10 minutes presentation about their report. Afterwards there is time for 30 minutes of discussion and last 5 minutes for giving a grade.

I was assigned 59 students in 6 groups. In 4 of the groups I was the only teacher. In the last 2 groups we were 3 teachers.

Below you will find some experiences collected during the course. Please post a comment if you have suggestions for improvements. Læs mere »

© Lars Relund Nielsen
Entries RSS Comments RSS Log in