Solución de un modelo de transporte mediante un algoritmo de asignación

Es conocido, que los problemas de asignación puros son una variación del problema original del transporte, en cuyos casos las variables de decisión solo pueden tomar valores binarios. ¿Esto qué significa? Pues bien, como algoritmo de red, un problema básico de transporte es un modelo de asignación genérico; modelo en el cual se pretende establecer la asignación entre fuentes y destinos; es decir, cuántas unidades se transportarán desde el origen i hacia el destino j. Se entiende además, que en el caso del transporte, estas variables de asignación no son necesariamente binarias.

El objetivo de este artículo consiste en utilizar las librerías del software Google OR-Tools (Python), para abordar un problema de transporte básico, de acuerdo a un modelo de asignación.

El problema

Con el propósito de evaluar los resultados obtenidos a través del tratamiento de un problema técnicamente formulado y abordado, utilizaremos un caso descrito en el libro Investigación de Operaciones (9na edición), de Hamdy A. Taha (University of Arkansas, Fayetteville), (Ejemplo 5.1-1):

MG Auto cuenta con tres plantas en Los Ángeles, Detroit y Nueva Orleáns, y dos importantes centros de distribución en Denver y Miami. Las capacidades trimestrales de las tres plantas son 1000, 1500 y 1200 automóviles, y las demandas de los dos centros de distribución durante el mismo periodo son de 2300 y 1400 automóviles. La distancia en millas entre las plantas y los centros de distribución aparece en la siguiente tabla: Hamdy A. Taha

Relación de distancias entre plantas y cedis

Distancia dada en millas

El caso también plantea que la compañía de transporte cobra 8 centavos por milla por automóvil. Así entonces, el objetivo del modelo no será minimizar la distancia total, sino minimizar el costo total de transporte. Con base en la distancia dada en millas y el costo de transporte unitario, construimos nuestra tabla de costos (redondeada al dólar más cercano):

Relación de costos de distribución entre plantas y cedis

Costo dado en dólares automóvil

La representación gráfica de la red sería la siguiente:

Resolviendo un problema de transporte mediante Or Tools (Asignación)

De acuerdo a lo mencionado en el artículo de introducción a Google OR-Tools, esta herramienta soporta múltiples lenguajes de programación, así entonces, haremos uso del lenguaje de programación Python.

Utilizaremos un solucionador para programación lineal mixta: SCIP.

Importar librerías

El primer paso consiste en importar la librería de Google Or Tools, esto nos permitirá utilizar todas sus funciones.

from ortools.linear_solver import pywraplp

Declarar el solucionador

Tal como lo mencionamos, vamos a utilizar un solucionador de programación lineal mixta. Esto nos permitirá utilizar tanto variables enteras como variables continuas. Debemos declarar el solucionador SCIP.

solver = pywraplp.Solver.CreateSolver('SCIP')

Crear la data del modelo

La data base d el modelo corresponde a la matriz de costos. Debemos de considerar que el orden en el que se ingresa la información corresponda a la relación precisa entre las variables; veamos:

costos = [
    [ 80, 215],
    [100, 108],
    [102,  68],
]

num_plantas = len(costos)
num_cedis = len(costos[0])

Así mismo, utilizaremos el método len para calcular el número de plantas y el número de cedis. En el caso de las plantas, contará el número de filas de la matriz; en el caso de los cedis, contará en número de columnas (para saber que contará las columnas de la matriz se utiliza costos[0]).

Crear las variables del modelo

Este es un problema básico, y como tal las únicas variables necesarias serán las asignaciones. Es decir, las variables x[i, j], donde i representa el origen y j representa el destino. Este proceso es relativamente sencillo y puede generarse cada variable de manera manual, sin embargo, automatizaremos la creación de las variables por medio de ciclos, con el propósito de preparar nuestro código para problemas más robustos:

x = {}
for i in range(num_plantas):
    for j in range(num_cedis):
        x[i, j] = solver.IntVar(0, solver.infinity(), '')

De esta manera se crean las variables de asignación. Es decir, tendremos tantos i como número de plantas (orígenes), asociados a tantos j como número de cedis (destinos). La naturaleza de las variables es entera (solver.IntVar), ya que hace referencia a la distribución de automóviles. Desde la creación de las variables podemos definir su rango, en este caso diremos que será desde 0 hasta infinito (solver.infinity).

Así entonces, se entiende que, por ejemplo, la variable x[0,0] hace referencia al número de automóviles que se enviarán desde el origen 0 (Los Ángeles) hacia el destino 0 (Denver). El orden de la asignación se basa en la matriz de costos.

El código anterior crea todas las variables de asignación necesarias.

Crear las restricciones

#Restricciones de disponibilidad (oferta en plantas)
solver.Add(solver.Sum([x[0, j] for j in range(num_cedis)]) <= 1000)
solver.Add(solver.Sum([x[1, j] for j in range(num_cedis)]) <= 1500)
solver.Add(solver.Sum([x[2, j] for j in range(num_cedis)]) <= 1200)

#Restricciones de demanda (cedis)
solver.Add(solver.Sum([x[i, 0] for i in range(num_plantas)]) >= 2300)
solver.Add(solver.Sum([x[i, 1] for i in range(num_plantas)]) >= 1400) 

El método Add de solver, nos servirá para agregar las restricciones al modelo. Las restricciones de oferta y demanda necesarias se adicionan tal cual el código anterior. Explicaremos la primera de ellas:

solver.Add(solver.Sum([x[0, j] for j in range(num_cedis)]) <= 80)

La anterior restricción indica que, la sumatoria de todas las variables x[0, j], para los j del conjunto cedis deberá ser menor o igual a 1000. Esto significa que la cantidad enviada desde el nodo (planta) 0 hacia los destinos j (Denver y Miami) no puede ser mayor que su disponibilidad de automóviles.

Crear la función objetivo

Tal como lo planteamos desde la formulación del problema, el criterio de optimización des este modelo corresponde a minimizar el costo total de distribución. El siguiente fragmento calcula el producto entre cada variable x[i, j] y su respectivo costo c[i, j] (obtenida desde la matriz de costos). La sumatoria de productos corresponde a la ecuación objetivo.

for i in range(num_plantas):
    for j in range(num_cedis):
        objective_terms.append(costos[i][j] * x[i, j])
solver.Minimize(solver.Sum(objective_terms))

Luego se declara el criterio de optimización: Minimize (para minimizar) y Maximize (para maximizar). En este caso minimizará la sumatoria de los productos.

Invocar al solucionador

La siguiente línea indica que el modelo debe solucionarse.

status = solver.Solve()

Configurar las salidas del modelo

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Costo total = ', solver.Objective().Value(), '\n')

    for i in range(num_plantas):
        for j in range(num_cedis):
            print('|Planta {} -> Cedi {} - Cantidad: {}|'.format(
            i,
            j,
            x[i, j].solution_value()))

En el caso en el que el modelo encuentre una solución óptima, hemos configurado las salidas de tal manera que nos muestre la relación solución entre fuentes > destinos y cantidades. Al ejecutar el modelo obtendremos la siguiente solución:

Podemos también detallar mucho mejor las salidas del modelo, de acuerdo al nombre de las plantas y los cedis. Para ello, en la data de entrada creamos las listas correspondientes, veamos:

costos = [
    [ 80, 215],
    [100, 108],
    [102,  68],
]

num_plantas = len(costos)
num_cedis = len(costos[0])

plantas = ['Los Ángeles', 'Detroit', 'New Orleans']
cedis = ['Denver', 'Miami']

La información de las listas debe guardar correspondencia con el orden de la matriz de costos.

Ahora, configuramos nuevamente las salidas del modelo:

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Costo total = ', solver.Objective().Value(), '\n')

    for i in range(num_plantas):
        for j in range(num_cedis):
            print('|{:^20} -> {:^20} | Cantidad: {:^20}|'.format(
            plantas[i],
            cedis[j],
            x[i, j].solution_value()))

En este caso hemos remplazado el índice solo de plantas y cedis, por el índice dentro del listado de plantas y cedis. Así entonces, tendremos la siguiente salida del modelo:

Pueden cotejar los resultados obtenidos los cuales corresponden a los mismos resultados presentados por Taha, quien utiliza TORA como solucionador.

¿Cómo ejecutar el modelo?

Alternativa 1, ejecución en nuestro equipo:

Lo primero que debemos considerar, en el caso de que queramos ejecutar este código en nuestro equipo, es que es preciso contar con la instalación de Python 

Instalación de OR-Tools para Python

Ahora, lo recomendable es trabajar con algún editor de código práctico (IDE), por ejemplo: Sublime Text, o Spyder (Una herramienta más completa y por ende más robusta y pesada).

Alternativa 2, ejecución en un entorno virtual (Recomendado):

Podemos utilizar del mismo modo, un entorno virtual. En este caso recomendamos el uso de Colaboratory de Google, un entorno que cuenta con todas las herramientas necesarias para nuestros desarrollos. No tendremos que instalar nada en nuestro equipo, y aprovecharemos la potencia de las máquinas de Google.

El código completo de nuestro desarrollo lo presentamos a continuación.

# Importar la librería de Google OR-Tools
from ortools.linear_solver import pywraplp

# Declarar el solucionador que abordará el modelo
solver = pywraplp.Solver.CreateSolver('SCIP')

# Data del modelo
costos = [
    [ 80, 215],
    [100, 108],
    [102,  68],
]

plantas = ['Los Ángeles', 'Detroit', 'New Orleans']
cedis = ['Denver', 'Miami']

num_plantas = len(costos)
num_cedis = len(costos[0])

# Variables del modelo
x = {}
for i in range(num_plantas):
    for j in range(num_cedis):
        x[i, j] = solver.IntVar(0, solver.infinity(), '')
        
#Restricciones de disponibilidad (oferta en plantas) 
solver.Add(solver.Sum([x[0, j] for j in range(num_cedis)]) <= 1000) 
solver.Add(solver.Sum([x[1, j] for j in range(num_cedis)]) <= 1500) 
solver.Add(solver.Sum([x[2, j] for j in range(num_cedis)]) <= 1200) 

#Restricciones de demanda (cedis) 
solver.Add(solver.Sum([x[i, 0] for i in range(num_plantas)]) >= 2300) 
solver.Add(solver.Sum([x[i, 1] for i in range(num_plantas)]) >= 1400)      



# Función objetivo con criterio de optimización: minimizar
objective_terms = []
for i in range(num_plantas):
    for j in range(num_cedis):
        objective_terms.append(costos[i][j] * x[i, j])

solver.Minimize(solver.Sum(objective_terms))

# Invoca el solucionador
status = solver.Solve()

# Configura los parámetros de impresión, las salidas del modelo
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Costo total = ', solver.Objective().Value(), '\n')

    for i in range(num_plantas):
        for j in range(num_cedis):
            print('|{:^20} -> {:^20} | Cantidad: {:^20}|'.format(
            plantas[i],
            cedis[j],
            x[i, j].solution_value()))

También es posible integrar nuestro desarrollo con data proveniente de fuentes externas, como por ejemplo Excel.

En próximos artículos, resolveremos el problema de transporte básico mediante un algoritmo de flujo mínimo.