¿Cómo calcular una matriz de distancias para modelar un VRP?

Tal como lo hemos abordado ampliamente, una de las aplicaciones más importantes del modelamiento de Cadenas de Suministro, es el diseño de red de abastecimiento, y dentro de esta categoría, el diseño de rutas de transporte (enrutamiento de vehículos).

Los problemas de enrutamiento de vehículos (routing), se encuentran clasificados como problemas de optimización combinatoria, y esto producto de que la cantidad de rutas posibles en un modelo básico, se encuentra determinado por la ecuación (n – 1)!, donde n, es igual al número de ubicaciones que componen el problema de enrutamiento.

Ahora bien, no solo la resolución de los problemas de enrutamiento es considerada como compleja, también lo es la fase preliminar en la que debemos levantar la información de entrada del modelo (inputs), específicamente, la matriz de distancias.

La forma más básica de un modelo VRP, requiere para su resolución la siguiente información mínima de entrada:

La matriz de distancias es un tabulado que registra la longitud entre cada uno de los nodos del modelo, incluido el nodo depósito. Su tamaño será n x n (siendo n el número de nodos del modelo). Asumiendo, por ejemplo, que tenemos un problema que se compone de 1 depósito y 19 ubicaciones, el modelo requerirá de una matriz de distancias con 400 datos.

En artículos anteriores hemos abordado modelos robustos para la resolución de problemas VRP y algunas de sus extensiones (CVRP, por ejemplo). El objetivo de este artículo será la de utilizar funciones espaciales de Python, para obtener una matriz de distancias de acuerdo a ubicaciones reales, con el propósito de resolver un modelo VRP básico.

Caso de aplicación

La Secretaría de Educación de Santiago de Cali, tiene dentro de sus funciones, distribuir el material pedagógico a las diferentes instituciones de formación de la ciudad. La distribución del material se realiza todos los lunes, y para ello, la Secretaría cuenta con dos vehículos con capacidad suficiente para transportar todo el material. El depósito del material (lugar desde donde salen y deben regresar los vehículos), se encuentra en la Alcaldía de Cali; y las instituciones que deben visitarse son 60. Tal como podemos apreciar en la siguiente tabla:

Se desea desarrollar un plan de rutas para cada uno de los vehículos que logre minimizar la distancia total recorrida, al tiempo que asegure la visita de todos los colegios.

¿Qué necesitaremos?

En el desarrollo de este ejercicio emplearemos:

Paso 1: Crear el entorno de trabajo en Colaboratory

Lo primero que vamos a hacer consiste en crear un entorno de trabajo en Google Colaboratory, así que vayamos allá: Abrir cuaderno nuevo.

Verán que tienen un lienzo para programar el modelo, así que en este cuaderno podemos ir generando las líneas de código que explicaremos en los pasos siguientes.

Paso 2: Importar las librerías necesarias

Respecto a las librerías, en la introducción del artículo hicimos una descripción de la funcionalidad de cada una, veamos como importarlas en nuestro entorno:

from sklearn.neighbors import DistanceMetric
from math import radians
import pandas as pd
import numpy as np

De esta manera, tenemos todo lo necesario para empezar a desarrollar nuestro código.

Paso 3: Importar los datos desde Excel

De acuerdo a las necesidades del modelo, podemos desarrollar un código que permita la entrada manual de la información, la captura de los datos desde entornos digitales (Internet, por ejemplo), o podemos, desde luego, alimentar nuestro modelo con información contenida en documentos externos, como es el caso de un archivo de Microsoft Excel.

Esta puede considerarse como una de las ventajas de utilizar Python, su capacidad de integrarse con cualquier fuente de datos. En nuestro caso, toda la información se encuentra contenida en un documento de Excel, el cual presenta el siguiente formato:

Utilizaremos ubicaciones reales, y para eso emplearemos las coordenadas de latitud y longitud.

En Colaboratory, el siguiente fragmento permitirá cargar un archivo al entorno de ejecución:

from google.colab import files

uploaded = files.upload()

for fn in uploaded.keys():
  print('User uploaded file "{name}" with length {length} bytes'.format(
      name=fn, length=len(uploaded[fn])))

Al ejecutar este fragmento de código, se abrirá una ventana emergente del explorador que permitirá cargar nuestra base de datos, en nuestro caso el archivo tienen el nombre de colegios.xlsx.

La siguiente línea de código permitirá almacenar los datos contenidos en el documento en un Dataframe de nuestro entorno, dentro de la variable data.

#Leer el documento de Excel y almacenar los datos en la variable data
data = pd.read_excel('colegios.xlsx')

Podemos en cualquier momento confirmar si la carga de los datos se ha realizado correctamente, para eso imprimiremos las primeras cinco filas del  DataFrame:

data.head()

Al ejecutar esta instrucción tenemos la siguiente salida (Una vista de las 5 primeras filas del marco de datos):

Podemos observar que los datos han sido perfectamente cargados, y que ahora se encuentran almacenados en la variable (DataFrame): data.

Paso 4: Convertir las coordenadas de latitud y longitud en radianes

La mayor parte de las funciones de cálculo de distancias de Sklearn (Scipy) toman las entradas como radianes. Por esta razón, debemos convertir las coordenadas en radianes.

data['Latitud'] = np.radians(data['Latitud'])
data['Longitud'] = np.radians(data['Longitud'])

#Creamos una matriz bidimensional con la latitud y la longitud en radianes
data[['Latitud','Longitud']].to_numpy()

Al ejecutar estas líneas, las coordenadas quedarán convertidas en radianes.

Paso 5: Declarar el tipo de métrica de distancias que se utilizará

En este punto quiero detenerme para mencionar que existen decenas de funciones métricas de distancia rápida. Algunas de las más utilizadas son destinadas a espacios vectoriales de valor real, como: distancias euclidianas, distancias de Manhattan, etc. Prácticas cuando se emplean con coordenadas cartesianas.

En nuestro caso, ya que utilizamos ubicaciones reales y contamos con coordenadas de latitud y de longitud, podemos emplear una función de distancia de vectores bidimensionales que considere la curvatura de la tierra; tal es el caso de las Distancias Haversine (Semiverseno).

Si quieren conocer la fórmula empleada para el cálculo de cada distancia haversine:

Para efectos de nuestro desarrollo, utilizaremos la librería SKLearn para calcular nuestras distancias. Veamos:

dist = DistanceMetric.get_metric('haversine')

dist.pairwise(data[['Latitud','Longitud']].to_numpy())*6373

El fragmento anterior crea un objeto con métricas haversine, y luego utiliza la función pairwise(), para calcular la distancia entre cada uno de los elementos de la matriz (nodos). Cada distancia calculada multiplica el escalar 6373 (radio esférico de la tierra), para calcular las distancias en kilómetros (Para millas multiplicar por 3798).

Paso 6: Crear un marco de datos (tabulado) de matriz de distancias

Una vez que ejecutemos el paso anterior, tendremos nuestras distancias entre nodos calculadas; lo que haremos en este paso será organizar dichos valores en forma de marco de datos o tabulado, obteniendo nuestra matriz de distancias.

distance_matrix = pd.DataFrame(dist.pairwise(data[['Latitud','Longitud']].to_numpy())*6373)

distance_matrix.head()

Al ejecutar esta instrucción tenemos la siguiente salida (Una vista de las 5 primeras filas del marco de datos):

Podemos imprimir la matriz completa:

print(distance_matrix)

Hasta este punto hemos logrado nuestro objetivo principal, que era obtener una matriz de distancias de forma automática tomando como base las coordenadas de longitud y latitud de un conjunto de nodos. Específicamente hemos obtenido 3721 valores de distancia en cuestión de segundos (empleando la métrica haversine).

Lógicamente, estos valores de referencia en la práctica presentan algunas desventajas, como, por ejemplo, la dificultad subyacente de atravesar las ciudades por lugares diferentes que las vías dispuestas para ello. Sin embargo, estos valores pueden ser muy útiles como datos de entrada de un modelo VRP para obtener la secuencia del plan de rutas.

Lo siguiente que haremos, para finalizar nuestro caso de aplicación, será utilizar la matriz de distancias obtenida, en un modelo VRP básico (Si desea profundizar en el modelo, visite el artículo).

Paso 7: Instalar Google Or Tools

Es necesario instalar la librería de Google Or Tools en nuestro entorno de Colaboratory para poder utilizar nuestro modelo VRP.

!pip install ortools

Al ejecutar esta instrucción instalaremos nuestro solucionador del modelo de enrutamiento.

Paso 8: Incorporar el modelo VRP (Previamente formulado)

Como ya lo mencionamos, contamos con un modelo debidamente formulado para resolver problemas VRP básicos. Lo único que modificaremos serán los datos de entrada: Matriz de distancias (distance_matrix), Número de vehículos (2), Depósito (Nodo 0).

"""Problema de enrutamiento de vehículos simple (VRP)

Autor: MSc. Ing. Bryan Salazar López 2021

"""

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def create_data_model():
    """Datos de entrada del modelo"""
    #Llamaremos la matriz de distancia previamente obtenida
    #Emplearemos dos vehículos como lo indica el problema
    #Definiremos el nodo 0 como el depósito (Secretaría)
    data = {}
    data['matriz_distancias'] = distance_matrix
    data['num_vehiculos'] = 2
    data['deposito'] = 0
    return data


def print_solution(data, manager, routing, solution):
    """Imprime la solución sobre la consola"""
    max_route_distance = 0
    for vehicle_id in range(data['num_vehiculos']):
        index = routing.Start(vehicle_id)
        plan_output = 'Ruta para el vehículo {}:\n'.format(vehicle_id)
        route_distance = 0
        while not routing.IsEnd(index):
            plan_output += ' {} -> '.format(manager.IndexToNode(index))
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id)
        plan_output += '{}\n'.format(manager.IndexToNode(index))
        plan_output += 'Distancia de la ruta: {}km\n'.format(route_distance)
        print(plan_output)
        max_route_distance += route_distance
        max_route_distance = max(route_distance, max_route_distance)
    print('Distancia total de todas las rutas: {}km'.format(max_route_distance))



def main():
    """Punto de entrada del programa"""
    # Invocar la data de entrada.
    data = create_data_model()

    # Crea el administrador del índice de rutas.
    manager = pywrapcp.RoutingIndexManager(len(data['matriz_distancias']),
                                           data['num_vehiculos'], data['deposito'])

    # Crea el modelo de enrutamiento.
    routing = pywrapcp.RoutingModel(manager)


    # Crea y registra una devolución de llamada de distancia.
    def distance_callback(from_index, to_index):
        """Retorna la distancia entre dos nodos."""
        # Convierte desde la variable de ruta Index hasta la matriz de distancia NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['matriz_distancias'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define el costo de cada arco.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Adhiere la dimensión de distancia.
    dimension_name = 'Distancia'
    routing.AddDimension(
        transit_callback_index,
        0,  # Sin holgura
        3000,  # Distancia máxima de viaje del vehículo
        True,  # Iniciar el acumulador en cero
        dimension_name)
    distance_dimension = routing.GetDimensionOrDie(dimension_name)
    distance_dimension.SetGlobalSpanCostCoefficient(100)

    # Configurar los parámetros de búsqueda.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Solucionador del problema.
    solution = routing.SolveWithParameters(search_parameters)

    # Imprimir la solución en la consola.
    if solution:
        print_solution(data, manager, routing, solution)
    else:
        print('No se encuentra solución !')


if __name__ == '__main__':
    main()

Al ejecutar el modelo tendremos:

Un plan de rutas para los dos vehículos que parten desde los depósitos (ahí mismo finalizan sus recorridos), y visitan la totalidad de los nodos. La distancia total optimizada es equivalente a 33 km.

Por último, también es posible exportar la matriz de distancias que hemos obtenido, para eso utilizaremos el siguiente fragmento:

distance_matrix.to_csv('distance_matrix.csv')
files.download('distance_matrix.csv')

El código completo de este desarrollo lo puedes encontrar en nuestro cuaderno: Matriz de distancias para modelar un VRP.