Implementación, resolución y visualización del problema del viajante con Python | por Carlos J. Uribe

📖 Si no leíste el Artículo anterior, no te preocupes. La formulación matemática también es declarado (pero no derivado)) aquí, con cada componente del modelo junto a su implementación de código. Si no comprende de dónde vienen o significan las cosas, lea el artículo de “carrera 2”, y si desea obtener más contexto, lea el artículo sobre “carrera 1“.

Python se utiliza porque es el lenguaje principal en ciencia de datos y Piomo ya que es una de las mejores bibliotecas (de código abierto) que se ocupa eficazmente de modelos a gran escala.

En esta sección, repasaré cada uno componente del modelo definido en la formulación, y explica cómo se traduce al código Pyomo.. Intenté no dejar espacios en blanco, pero si no crees lo contrario, deja una pregunta en los comentarios.

Descargo de responsabilidad: Se espera que el lector objetivo sea nuevo en Pyomo, e incluso en el modelado, para reducir su carga cognitiva. Se prioriza la implementación concisa y simple sobre las mejores prácticas de programación.. El objetivo ahora es enseñar modelado de optimización, no ingeniería de software. El código se mejora gradualmente en futuras iteraciones a medida que esta prueba de concepto evoluciona hacia un MVP más sofisticado.

1.1. Instalando las dependencias

Para gente con prisa

Instale (o asegúrese de tener ya instaladas) las bibliotecas pyomo, networkx y pandasy el paquete glpk.

📝 El paquete glpk contiene el solucionador de GLPKque es un (código abierto) solucionador externo usaremos para optimizar los modelos que creamos. Pyomo está acostumbrado a crear modelos de problemas y pasarlos a GLPK, que ejecutará los algoritmos que llevan a cabo el proceso de optimización en sí. Luego, GLPK devuelve las soluciones al objeto del modelo de Pyomo, que se almacenan como atributos del modelo, para que podamos usarlos cómodamente sin salir de Python.

La forma recomendada de instalar GLPK es a través de conda para que Pyomo pueda encontrar fácilmente el solucionador GLPK. Para instalar todas las dependencias juntas de una sola vez, ejecute:

conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk

Para gente organizada

Recomiendo crear un separado ambiente virtual en el que instalar todas las bibliotecas necesarias para seguir los artículos de esta serie. Copia este texto

name: ttp  # traveling tourist problem
channels:
- conda-forge
dependencies:
- python=3.9
- pyomo=6.5.0
- pandas
- networkx
- glpk  # external solver used to optimize models
- jupyterlab  # comment this line if you won't use Jupyter Lab as IDE

y guárdelo en un archivo YAML llamado environment.yml. Abra un indicador de Anaconda en la misma ubicación y ejecute el comando

conda env create --file environment.yml

Después de unos minutos, se crea el entorno con todas las dependencias instaladas en su interior. Correr conda activate ttp Para “entrar” en el entorno, inicie Jupyter Lab (ejecutando jupyter lab en la terminal) y ¡comienza a codificar!

1.2. Las matemáticas se convierten en código

Primero, asegurémonos de que Pyomo pueda encontrar el solucionador GLPK

### =====  Code block 3.1  ===== ###
import pandas as pd
import pyomo.environ as pyo
import pyomo.version
import networkx as nx
solver = pyo.SolverFactory("glpk")
solver_available = solver.available(exception_flag=False)
print(f"Solver '{solver.name}' available: {solver_available}")
if solver_available:
print(f"Solver version: {solver.version()}")
print("pyomo version:", pyomo.version.__version__)
print("networkx version:", nx.__version__)

Solver 'glpk' available: True
Solver version: (5, 0, 0, 0)
pyomo version: 6.5
networkx version: 2.8.4

⛔ Si recibiste el mensaje 'glpk' available: False, el solucionador no se instaló correctamente. Pruebe uno de estos para solucionar el problema:

– Vuelva a realizar los pasos de instalación con cuidado.

– correr conda install -y -c conda-forge glpk en el entorno base (predeterminado)

– intente instalar un solucionador diferente que funcione para usted

Luego lea el archivo CSV de datos de distancia

### =====  Code block 3.2  ===== ###
path_data = (
"https://raw.githubusercontent.com/carlosjuribe/"
"traveling-tourist-problem/main/"
"Paris_sites_spherical_distance_matrix.csv"
)
df_distances = pd.read_csv(path_data, index_col="site")
df_distances

Ahora entramos en la “etapa 4” de la [Agile Operations Research workflow]marcado como el bloque verde a continuación:

La tarea consiste en tomar el modelo matemático creado anteriormente e implementarlo en código. exactamente como fue definido matemáticamente.

👁️ Se nos permite crear tantos objetos Python como queramos si eso facilita la implementación del modelo, pero No se nos permite modificar el modelo subyacente de ninguna manera mientras lo codificamos.. Causaría que el modelo matemático y el modelo de computadora no estuvieran sincronizados.lo que dificulta bastante la depuración de modelos posteriores.

Creamos una instancia de un modelo Pyomo vacío, en el que almacenaremos los componentes del modelo como atributos:

model = pyo.ConcreteModel("TSP")

1.2.1. Conjuntos

Para crear el conjunto de sitios 𝕊 = {Louvre, Tour Eiffel,…, hotel}, extraemos sus nombres del índice del marco de datos y lo usamos para crear un Pyomo. Set llamado sites:

### =====  Code block 3.3  ===== ###
list_of_sites = df_distances.index.tolist()
model.sites = pyo.Set(initialize=list_of_sites, 
domain=pyo.Any, 
doc="set of all sites to be visited (𝕊)")

Para crear el conjunto derivado

Guardamos el filtro 𝑖 ≠ 𝑗 dentro de un regla de construcción (la función de Python _rule_domain_arcs), y pasar esta regla al filter palabra clave al inicializar el Set. Tenga en cuenta que este filtro se aplicará al producto cruzado de los sitios (𝕊 × 𝕊) y filtrará aquellos miembros del producto cruzado que no cumplan la regla.

### =====  Code block 3.4  ===== ###
def _rule_domain_arcs(model, i, j):
""" All possible arcs connecting the sites (𝔸) """
# only create pair (i, j) if site i and site j are different
return (i, j) if i != j else None 
rule = _rule_domain_arcs
model.valid_arcs = pyo.Set(
initialize=model.sites * model.sites,  # 𝕊 × 𝕊
filter=rule, doc=rule.__doc__)

1.2.2. Parámetros

El parámetro

𝐷ᵢⱼ ≔ Distancia entre el sitio 𝑖 y el sitio 𝑗

se crea con el constructor pyo.Paramque toma como la primera argumento (posicional) el dominio 𝔸 (model.valid_arcs) que lo indexa, y como el palabra clave argumento initialize otra regla de construcción, _rule_distance_between_sites, que se evalúa para cada par (𝑖, 𝑗) ∈ 𝔸. En cada evaluación, el valor numérico de la distancia se obtiene del marco de datos. df_distancesy vinculado internamente al par (𝑖, 𝑗):

### =====  Code block 3.5  ===== ###
def _rule_distance_between_sites(model, i, j):
""" Distance between site i and site j (𝐷𝑖𝑗) """
return df_distances.at[i, j]  # fetch the distance from dataframe
rule = _rule_distance_between_sites
model.distance_ij = pyo.Param(model.valid_arcs, 
initialize=rule, 
doc=rule.__doc__)

1.2.3. Variables de decisión

Debido a que 𝛿ᵢⱼ tiene el mismo “dominio de índice” que 𝐷ᵢⱼ, la forma de construir este componente es muy similar, con la excepción de que aquí no se necesita ninguna regla de construcción.

model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary, 
doc="Whether to go from site i to site j (𝛿𝑖𝑗)")

El primer argumento posicional de pyo.Var está reservado para su conjunto de índices 𝔸, y el “tipo” de variable se especifica con el argumento de palabra clave within. Con “tipo de variable” me refiero a la rango de valores que la variable puede tomar. Aquí, el rango de 𝛿ᵢⱼ es solo 0 y 1, por lo que es de tipo binario. Matemáticamente, escribiríamos 𝛿ᵢⱼ ∈ {0, 1}, pero en lugar de crear restricciones separadas para indicar esto, podemos indicarlo directamente en Pyomo configurando within=pyo.Binary en el momento que creamos la variable.

1.2.4. Función objetiva

Para construir una función objetivo podemos “almacenar” la expresión en una función que se utilizará como regla de construcción para el objetivo. Esta función solo toma un argumento, el modelo, que se usa para buscar cualquier componente del modelo necesario para construir la expresión.

### =====  Code block 3.6  ===== ###
def _rule_total_distance_traveled(model):
""" total distance traveled """
return pyo.summation(model.distance_ij, model.delta_ij)
rule = _rule_total_distance_traveled
model.obj_total_distance = pyo.Objective(rule=rule, 
sense=pyo.minimize, 
doc=rule.__doc__)

Observe el paralelismo entre la expresión matemática y la declaración de retorno de la función. Especificamos que este es un problema de minimización con el sense palabra clave.

1.2.5. Restricciones

Si recuerda el artículo anterior, una forma conveniente de hacer cumplir que cada sitio se visite solo una vez es hacer que se “ingrese” a cada sitio una vez y se “salga” una vez, simultáneamente.

Cada sitio se ingresa solo una vez.

def _rule_site_is_entered_once(model, j):
""" each site j must be visited from exactly one other site """
return sum(model.delta_ij[i, j] for i in model.sites if i != j) == 1
rule = _rule_site_is_entered_once
model.constr_each_site_is_entered_once = pyo.Constraint(
model.sites, 
rule=rule, 
doc=rule.__doc__)

Se sale de cada sitio solo una vez

def _rule_site_is_exited_once(model, i):
""" each site i must departure to exactly one other site """
return sum(model.delta_ij[i, j] for j in model.sites if j != i) == 1
rule = _rule_site_is_exited_once
model.constr_each_site_is_exited_once = pyo.Constraint(
model.sites, 
rule=rule, 
doc=rule.__doc__)

1.2.6. Inspección final del modelo.

La implementación del modelo está hecha. Para ver cómo se ve en su conjunto, debemos ejecutar model.pprint()y navegue un poco para ver si nos perdimos algunas declaraciones o cometimos algunos errores.

Para tener una idea del tamaño del modelo, imprimimos el número de variables y restricciones que lo componen:

def print_model_info(model):
print(f"Name: {model.name}", 
f"Num variables: {model.nvariables()}",
f"Num constraints: {model.nconstraints()}", sep="\n- ")
print_model_info(model)
#[Out]:
# Name: TSP
#  - Num variables: 72
#  - Num constraints: 18

Al tener menos de 100 restricciones o variables, este es un problema de tamaño pequeño y el solucionador lo optimizará relativamente rápido.