Problema de enrutamiento de vehículos

El problema de enrutamiento de vehículos (VRP por sus siglas en inglés) es un problema de optimización logística que consiste en hayar el conjunto óptimo de rutas para que una flota de vehículos entregue los pedidos de un conjunto de clientes. Es lo que se conoce como un problema NP-duro de optimización combinatoria y tiene múltiples aplicaciones reales en la industria.

A continuación demuestro cómo se puede solucionar un VRP a pequeña escala usando R. Aunque es cierto que existen herramientas más desarrolladas en otros lenguajes (principalmente el paquete OT tools de Google) existe una librería de R que sirve para solucionar algunas versiones de este problema.

Esta librería es vrpoptima, y permite solucionar problemas de enrutamiento. En concreto permite solucionar las variantes de enrutamiento con capacidad (es decir con límites de capacidad de entrega, por ejemplo considerar que los vehículos pueden hacer un número máximo de entregas), con múltiples depósitos (más de un lugar de origen de los vehículos), y con restricción de distancia (considerar que los vehículos no deben superar un kilometraje máximo).

Manos a la obra

Para el ejemplo que quiero mostrar, he muestreado un conjunto de direcciones geolocalizadas en la provincia de Madrid usando los datos del servicio Cartociudad del Instituto Geográfico Nacional que representarán los puntos de entrega de pedidos de los clientes. Para el origen de los vehículos he geolocalizado manualmente cuatro centros logísticos de Amazon en la provincia usando Google Maps.

La siguiente imagen muestra los puntos de entrega (en negro) y los centros logísticos (rojo).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

# Load the package
library(dplyr)
library(ggplot2)
library(sf)
library(tidyterra)

orders <- read_sf("data/orders.gpkg")
vehicles <- read_sf("data/vehicles.gpkg")

tile <-
  mapSpain::esp_getTiles(
    orders,
    type = "IGNBase.Todo",
    zoommin = 2
  )

orders %>% 
  ggplot() +
  geom_spatraster_rgb(data = tile, maxcell = 10e6) +
  geom_sf() +
  geom_sf(data = vehicles, color = "red", size = 3) +
  theme_void()

Con este planteamiento, lo siguiente es calcular la matriz de distancias de los puntos de entrega. Se trata de una matriz de tantas filas y columnas como puntos de entrega en la que tenemos calculada la distancia de cada punto de entrega respecto al resto.

Lo ideal, por supuesto, es que esta matriz de distancia no mida tanto la distancia en linea recta, sino el tiempo de viaje. Esto es costoso porque requiere usar servicios de pago como Google o Mapbox. Ambos servicios tienen la posibilidad de calcular matrices de forma gratuita con ciertos límites.

Para el ejemplo que nos atañe utilizo Mapbox para calcular la matriz de tiempo de viaje. Utilizo la librería mapboxapi), que facilita hacer llamadas a la API de Mapbox.

1
2
3
4
5
6
7
8
9

# Matriz de tiempos de viaje usando la API de Mapbox
# Se necesita una clave API de Mapbox: https://account.mapbox.com/
dist_mat <-
  mapboxapi::mb_matrix(
    origins = orders,
    destinations = orders,
    fallback_speed = 30, # Velocidad cuando no se encuentra ruta entre dos puntos
    allow_large_matrix= T
  )

Una vez tenemos esta matriz podemos usar la función VehicleRouting del paquete vrpoptima. Esta función es muy flexible y permite introducir una serie de parámetros. Si queréis profundizar en cada uno de ellos os recomiendo leer la documentación, pero aquí explico algunos:

• cost_type: Controla el tipo de coste que se pretende optimizar. La opción 1 intenta minimizar la distancia acumulada de todos los vehículos. La opción 2 intenta minimizar la ruta más larga.
• max_tour_distance: La distancia máxima deseable para las rutas. Es una restricción de distancia para controlar que las rutas no sean demasiado largas. El algoritmo intentará que las rutas no sobrepasen este valor, pero podrán hacerlo con una penalización que se puede controlar con el parámetro penalty_factor.
• min_tour/max_tour_visits: El número de paradas mínimo máximo que pueden hacer los vehículos. Con esto se controla la restricción de capacidad del VRP.
• num_generations: El número de evoluciones de la solución. Valores más altos producen mejores soluciones, pero son más costosos computacionalmente.
• population_size: El número de soluciones que el algoritmo puede observar en cada paso. Valores más altos exploran un espacio de soluciones mayor, pero pueden ralentizar el proceso.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# La función VehicleRouting necesita los puntos como una matriz de coordenadas 
# y no como un objeto sf
orders_mat <- st_coordinates(orders)
vehicles_mat <- st_coordinates(vehicles)

colnames(orders_mat) <- c("lat", "lon")
colnames(vehicles_mat) <- c("lat", "lon")

library(vrpoptima)
solution <-  VehicleRouting(
  visit_points = orders_mat,
  num_agents = nrow(vehicles_mat),
  agent_points = vehicles_mat,
  cost_type = 2,
  max_tour_distance = Inf,
  max_tour_visits = 100,
  distance_metric = 'Geodesic',
  distance_matrix = dist_mat,
  min_tour = 2,
  population_size = 96,
  num_generations = 15000, 
  distance_truncation = TRUE, 
  seed = 42
)

Una vez que el algoritmo termina de evaluar las opciones posibles (dentro de los límites especificados con los parámetros num_generations y population_size) nos devuelve un objeto con la información sobre la mejor solución. En este caso, la suma total del tiempo estimado para las rutas de los 4 vehículos es de 3524 minutos (recordemos que la matriz empleada es una de tiempos de viaje en minutos, no de distancia). Esto sugiere que los vehículos tardarían una media de 15 horas en repartir los aproximadamente 64.5 pedidos asignados a cada uno. Esto es poco realista, pero tengamos en cuenta que el ejemplo tampoco lo es: hemos considerado que el reparto de 250 pedidos lo realizarían 4 vehículos.

1
2
3

print(solution$tour_lengths)

print( sum(solution$tour_lengths) )

Pero el objeto que nos ha devuelto la función VehicleRouting también contiene la ruta que deberá realizar cada vehículo. A continuación procesamos los datos para visualizar el resultado en un mapa.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

routes <- solution$routes
rownames(routes) <- 1:nrow(routes)
routes_list = RoutesDataPrep(routes = solution$routes, 
                             visit_points = orders_mat, 
                             agent_points = vehicles_mat)

cast_route <- function(route_list) {
  route_list %>% 
    st_as_sf(coords = c("lat", "lon"), na.fail = FALSE, crs = "epsg:4258") %>% 
    summarize(do_union=FALSE) %>%
    st_cast("LINESTRING") 
}

routes <- purrr::map_df(routes_list, cast_route, .id = "route")

destinations <- 
  routes %>% 
  st_cast("POINT") %>% 
  group_by(route) %>% 
  mutate(size = if_else(geometry == first(geometry), 3.5, 2.3))

Ahora visualizamos el resultado.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

routes %>% 
  ggplot() +
  geom_spatraster_rgb(data = tile, maxcell = 10e6) +
  geom_sf(aes(color = route), linewidth = 1.2) +
  geom_sf(aes(shape = as.factor(size), size = size+.8), data = destinations, color = "white") +
  geom_sf(aes(color = route, shape = as.factor(size), size = size), data = destinations) +
  facet_wrap(~route) +
  scale_size_identity() +
  theme_void() +
  theme(legend.position = "none")

Esto ha sido un ejemplo sencillo de cómo solucionar un problema de enrutamiento de vehículos a muy pequeña escala usando R. Sin embargo, para implementar estos algoritmos en casos reales en los que se trabaja con miles de pedidos sería necesario buscar formas de optimizar el número de paradas y el cálculo de la matriz de distancias. En esta publicación los autores proponen usar una versión recursiva del algoritmo de clustering dbscan para optimizar el cálculo de una matriz de distancia.