1 Intro

Routing optimization generally solves the Vehicle Routing Problem (a simple example being the more widely known Traveling Salesman Problem). A more complex example would be the distribution of goods by a fleet of multiple vehicles to dozens of locations, where each vehicle has certain time windows in which it can operate and each delivery location has certain time windows in which it can be served (e.g. opening times of a supermarket).

In this example we’ll look at a real-world scenario of distributing medical goods during disaster response following one of the worst tropical cyclones ever been recorded in Africa: Cyclone Idai.

Cyclone Idai floods in false color image on 19.03.2019; © Copernicus Sentinel-1 -satellite (modified Copernicus Sentinel data (2019), processed by ESA, CC BY-SA 3.0 IGO)source

In this scenario, a humanitarian organization shipped much needed medical goods to Beira, Mozambique, which were then dispatched to local vehicles to be delivered across the region. The supplies included vaccinations and medications for water-borne diseases such as Malaria and Cholera, so distribution efficiency was critical to contain disastrous epidemics.

We’ll solve this complex problem with the optimization endpoint of openrouteservice.

Before we get to action, we need to make sure the openrouteservice library is installed. If not run the chunk below.

install.packages(c("devtools", "dplyr","sf","geojsonsf","readr","units","googlePolylines",
                   "jsonlite","tmap","ggplot2",,"mapview"))
devtools::install_github("GIScience/openrouteservice-r")

Next we call all libraries required for this scernario.

# main libraries
library(dplyr) # manipulating tabular data
library(sf) # geospatial classes and datastructures - the swiss army knife on GIS in R
library(openrouteservice) # interface to openrouteservice API

# utils
library(geojsonsf) # sf classes to/from geojsons
library(units) # add units to numerics
library(readr) # open csvs, parse timestamps
library(googlePolylines) # decode binary polylines
library(jsonlite)

# # visualization
library(tmap) # visualize spatial data
library(ggplot2) # visualize (non-) spatial data
library(mapview)

In order to use the openrouteservice public api you need to register yourself an account and create an api key.

When you sucessfully created your personal api key, copy it over in the chunk below.

your_api_key <-
  "YOUR_API_KEY"

2 Logistics Setup

In total 20 sites were identified in need of the medical supplies, while 3 vehicles were scheduled for delivery. Let’s assume there was only one type of goods, e.g. standard moving boxes full of one medication. (In reality there were dozens of different good types, which can be modeled with the same workflow, but that would unnecessarily bloat this example).

The vehicles were all located in the port of Beira and had the same following constraints:

  • operation time windows from 8:00 to 20:00
  • loading capacity of 300 [arbitrary unit]

The delivery locations were mostly located in the Beira region, but some extended ~ 200 km to the north of Beira. Their needs range from 10 to 148 units of the arbitrary medication goods (consult the file located at data/idai_health_sites.csv.

In the Chunk below we load the csv file which contains the delivery locations. We also define the depot of our fleet in Beira.

# healthsites to be served
idai_health_sites <- read_csv("data/idai_health_sites.csv",
    col_types = cols(ID = col_integer(),
        Open_From = col_datetime(format = "%Y-%m-%d %H:%M:%S"),
        Open_To = col_datetime(format = "%Y-%m-%d %H:%M:%S"),
        Needed_Amount = col_integer()))


# The vehicles are all located at the port of Beira
depot <- data.frame(lon = 34.835447, lat = -19.818474)

idai_health_sites |> st_as_sf(coords=c("Lon", "Lat"), crs=4326) |>  mapview()

3 The routing problem setup

Now that we have described the setup sufficiently, we can start to set up our actual Vehicle Routing Problem. For this example we’re using the FOSS library of Vroom. VROOM is independently developed from openrouteservice but available through its API.

To properly describe the problem in algorithmic terms, we have to provide the following information:

  • vehicles start/end address: vehicle depot in Beira’s port
  • vehicle capacity: 300
  • vehicle operational times: 08:00 - 20:00
  • service location: delivery location
  • service time windows: individual delivery location’s time window
  • service amount: individual delivery location’s needs

We defineall these parameter down in the code and send a request to openrouteservice optimization endpoint at https://api.openrouteservice.org/optimization.

Let’s have a closer look at vehicles. With the public API we can use up to 3 different vehicles. They can be differentiated by:

  • Routing profile: heavy goods vehicles, regular cars, bikes, pedestrians
  • start and end locations
  • Skills: Does a vehicle offer a specific skill that is demanded by a delivery, like cooling capabilities for instance
  • Capacity: How much goods can be transported
  • Time window: During which time is the vehicle operational
vehicles = vehicles(
  id = 1:3,
  profile = "driving-hgv",
  start = depot,
  end = depot,
  capacity = 300,
  time_window = c("2019-03-22 05:00:00 UTC" |> as.POSIXct(tz = "UTC") |> as.integer(),
                  "2019-03-22 23:00:00 UTC" |> as.POSIXct(tz = "UTC") |> as.integer())
)

And here the delivery locations.

  • service: The amount of time for a delivery to be done at a location
  • amount: The amount of goods need at a delivery location
  • time_window: Periods in which the location can be visited/ a delivery conducted. Can be multiples
open_from <- idai_health_sites$Open_From |> as.integer()
open_to <- idai_health_sites$Open_To |> as.integer()
time_periods <- lapply(seq_along(open_from), function(i) list(c(open_from[i], open_to[i])))

jobs = jobs(
  id = idai_health_sites$ID,
  service = 1200,
  amount = idai_health_sites$Needed_Amount |> lapply(function(x) c(x)),
  location = Map(c,idai_health_sites$Lon,
                   idai_health_sites$Lat),
  time_window = time_periods
  )

Example of how multiple delivery periods are reprsented in the request.

{
  time_window = [[
    # first time_window
    1553238000,
    1553239800],
    [
      # second time window if exist
      1553243400,
      1553248800]],
}

Next we can fire the request against the optimization endpoint of openrouteservice. The list passed to options parameter makes sure that we get the route geometries as well.

res <- ors_optimization(jobs, vehicles, options = list(g=T), api_key=your_api_key)

4 Postprocess

The response is a nested json / list. Therefore we have to employ some looping and lapply to get to our results. See the comments for more context in the Chunk below

# We loop over every route that is planned by VROOM
for (r in 1:length(res$routes)) {

  # Every route bears overview information of the vehicle and its planned trip

  vehicle_steps <-
    res$routes[[r]]$steps |> toJSON(auto_unbox = T) |> fromJSON()

  # Every trip is divided into steps, the travel between delivery locations
  # It contains information like arrival time, load, distance, and the geometry binary encoded

  vehicle_steps$vehicle_id <- res$routes[[r]]$vehicle

  # here we decode it
  vehicle_route <-
    googlePolylines::decode(res$routes[[r]]$geometry)[[1L]] |> st_as_sf(coords = c("lon", "lat"), crs = 4326) |> st_combine() |> st_cast("LINESTRING")
  vehicle_route <- st_sf(vehicle_id = r, geometry = vehicle_route)

  # if first iteration, overwrite, if not append
  if (r == 1) {
    steps <- vehicle_steps
    routes <- vehicle_route
  } else {
    steps <- rbind(steps, vehicle_steps)
    routes <- rbind(routes, vehicle_route)
  }

  routes <- rbind(routes, routes)
}

# not visited locations are found in the unassigned list
unassigned <- do.call(rbind,
                      lapply(res$unassigned, function(x)
                        c(
                          id = x$id,
                          lon = x$location[1],
                          lat = x$location[2]
                        )),) |> data.frame() |> st_as_sf(coords = c("lon", "lat"), crs = 4326)

Now we have extracted the routes and each step into individual dataframes. Much easier to grasp than lists. Next we convert them to sf dataframes to be able to put on a map.

5 Results and evaluation

steps$lon <- lapply(steps$location, function(x) x[[1]])
steps$lat <- lapply(steps$location, function(x) x[[2]])

steps <- steps |> filter(type == "job") |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

steps <- steps |> select(-c(location))

5.1 Map visualization

Perfect, lets have look at our planned trips on a map.

mapview(routes,
        color = c("firebrick3", "darkgreen", "yellow1"),
        zcol = "vehicle_id") +
  mapview(steps, color = "grey75", label = "id") +
  mapview(unassigned, color="orange", label="id")

Looks like vehicle 1 will get far north and vehicle 3 takes care of the locations in the south. However with the current planning, one Job in the north east won’t be visited.

Can you adjust the parameters time_windows or amount to add it?

The following code chunk produces an overview of total quanitity, time and distance traveled by vehicle and serves as overview.

sumry <- do.call(rbind,
                 lapply(res$routes, function(x)
                   c(
                     vehicle = x$vehicle,
                     amount = x$amount,
                     duration = x$duration,
                     distance = x$distance
                   ))) |> data.frame()
sumry
##   vehicle amount duration distance
## 1       1    294    53047   915587
## 2       2    211     3661    42982
## 3       3    284    43053   635175

Vehicle 2 travelled a lot less than the other vehicles.

5.2 Time tables

The other Chunks show the individual timetables for every vehicle.

# timetable vehicle 1

ttable1 <- do.call(rbind,
        res$routes[[1]]$steps |> lapply(
          function(x)
            c(
              station_id = ifelse(x$type %in% c("start","end"),x$type,x$id),
              arrival = ifelse(x$type != "start",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character(),"-"),
              departure = ifelse(x$type == "end","-",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character()),
              load = x$load
            )
        )) |> data.frame()
ttable1
##   station_id             arrival           departure load
## 1      start                   - 2019-03-22 06:30:21  294
## 2         18 2019-03-22 07:00:00 2019-03-22 07:00:00  262
## 3         12 2019-03-22 07:38:32 2019-03-22 07:38:32  114
## 4         11 2019-03-22 08:31:26 2019-03-22 08:31:26   89
## 5         16 2019-03-22 14:31:12 2019-03-22 14:31:12   39
## 6         15 2019-03-22 15:31:02 2019-03-22 15:31:02    0
## 7        end 2019-03-22 22:54:28                   -    0
# timetable vehicle 2

ttable2 <- do.call(rbind,
        res$routes[[2]]$steps |> lapply(
          function(x)
            c(
              station_id = ifelse(x$type %in% c("start","end"),x$type,x$id),
              arrival = ifelse(x$type != "start",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character(),"-"),
              departure = ifelse(x$type == "end","-",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character()),
              load = x$load
            )
        )) |> data.frame()
ttable2
##   station_id             arrival           departure load
## 1      start                   - 2019-03-22 08:46:59  211
## 2          6 2019-03-22 09:04:45 2019-03-22 09:04:45  134
## 3          3 2019-03-22 09:28:59 2019-03-22 09:28:59   86
## 4          7 2019-03-22 10:00:00 2019-03-22 10:00:00   63
## 5          5 2019-03-22 10:34:03 2019-03-22 10:34:03   34
## 6          4 2019-03-22 10:58:50 2019-03-22 10:58:50    0
## 7        end 2019-03-22 11:28:00                   -    0
# timetable vehicle 3

ttable3 <- do.call(rbind,
        res$routes[[3]]$steps |> lapply(
          function(x)
            c(
              station_id = ifelse(x$type %in% c("start","end"),x$type,x$id),
              arrival = ifelse(x$type != "start",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character(),"-"),
              departure = ifelse(x$type == "end","-",x$arrival |>
                                 as.POSIXct(origin="1970-01-01", tz="UTC") |> as.character()),
              load = x$load
            )
        )) |> data.frame()
ttable3
##    station_id             arrival           departure load
## 1       start                   - 2019-03-22 05:00:00  284
## 2           1 2019-03-22 08:19:17 2019-03-22 08:19:17  260
## 3          17 2019-03-22 11:09:41 2019-03-22 11:09:41  212
## 4           8 2019-03-22 13:20:31 2019-03-22 13:20:31  195
## 5           9 2019-03-22 16:00:39 2019-03-22 16:00:39  113
## 6          10 2019-03-22 16:58:49 2019-03-22 16:58:49   69
## 7          13 2019-03-22 17:44:10 2019-03-22 17:44:10   36
## 8          14 2019-03-22 18:21:12 2019-03-22 18:21:12   10
## 9           2 2019-03-22 18:57:04 2019-03-22 18:57:04    0
## 10        end 2019-03-22 19:37:33                   -    0