Optimizando con Pyomo: una guía completa y paso a paso.

Configuración de Pyomo

Antes de sumergirnos en el modelado, necesitamos configurar nuestro entorno instalando Pyomo y eligiendo un solucionador adecuado.

1. Prerequisitos

Para usar pyomo, debes tener Python 3.6 o superior. Pyomo se puede instalar a través de pip.

pip install pyomo

Este tutorial ha sido creado utilizando la versión 6.8.0 de Pyomo.

import pyomo
print(pyomo.__version__)

Salida:

>>> 6.8.0

2. Elección y instalación del solucionador correcto

En optimización, los solucionadores son esenciales ya que son los algoritmos que encuentran la solución óptima al problema que has definido. Diferentes solucionadores se adecuan mejor dependiendo del tipo de problema (por ejemplo, lineal, no lineal, entero). Pyomo es una herramienta de modelado que depende de solucionadores externos para realizar la computación real.

Vamos a revisar algunos de los solucionadores más comunes.

Solucionadores de código abierto

1. GLPK (GNU Linear Programming Kit)

GLPK es una herramienta popular para resolver problemas de programación lineal (LP) y de programación entera mixta (MIP).

Es una excelente opción para tareas básicas de optimización lineal y es ampliamente utilizada en aplicaciones académicas y industriales.

Instalación

• Windows: Descargar e instalar GLPK.
• macOS: 

brew install glpk

• Linux: 

sudo apt-get install glpk-utils

2. CBC (Coin-or Branch and Cut)

CBC es un solucionador de código abierto para problemas de programación lineal (LP) y de programación entera mixta (MIP).

Ofrece características avanzadas y mejores prestaciones en algunos casos en comparación con GLPK, lo que lo hace una opción fuerte para tareas de optimización más complejas.

CBC se puede instalar a través del gestor de paquetes conda.

conda install -c conda-forge coincbc

3. IPOPT (Optimizador de Punto Interior)

IPOPT es un poderoso solucionador diseñado para problemas deprogramación no lineal (NLP) a gran escala.

Este es particularmente adecuado para manejar modelos complejos y no lineales, lo que lo hace una excelente opción para problemas que superan la optimización lineal.

IPOPT se puede instalar mediante el gestor de paquetes conda.

!conda install -c conda-forge ipopt

Solucionadores comerciales

1. CPLEX

CPLEX es un solucionador de optimización de punta de arte que maneja eficientemente problemas de programación lineal (LP), programación entera mixta (MIP) y problemas cuadráticos (QP). 

Requiere una licencia de IBM, pero está disponible gratuitamente para usuarios académicos, lo que lo convierte en una excelente opción para fines de investigación y educativos.

2. Gurobi

Gurobi es un solucionador comercial líder conocido por su velocidad y eficiencia en la resolución de problemas de LP, MIP, QP y programación no lineal (NLP).

Al igual que CPLEX, requiere una licencia pero ofrece acceso gratuito a usuarios académicos. Así, es una herramienta estándar de la industria para la optimización avanzada.

Solucionadores de código abierto vs. solucionadores comerciales.

Solucionadores de código abierto como GLPK y CBC son gratuitos y suficientes para la mayoría de las necesidades básicas de optimización. Son excelentes opciones para proyectos de pequeña escala y fines educativos.

Sin embargo, los solucionadores comerciales como CPLEX y Gurobi ofrecen generalmente un mejor rendimiento, especialmente para problemas más grandes y complejos. Estos solucionadores tienen características avanzadas, que incluyen soporte mejorado para programación cuadrática y no lineal, y están optimizados para aplicaciones industriales a gran escala.

Mientras que los solucionadores de código abierto pueden manejar muchas tareas de optimización rutinaria, los solucionadores comerciales son a menudo una mejor opción cuando se trata de requisitos más complejos y de alto rendimiento.

Ten en cuenta que los solucionadores comerciales requieren una licencia, aunque están disponibles gratuitamente para usuarios académicos.

Ahora, veamos cómo configurar un solucionador en Pyomo. En este caso usaremos GLPK.

3. Configurar un solucionador en Pyomo

Primero, asegúrese de que el ejecutable del solucionador está en su PATH del sistema después de la instalación.

A continuación, cree un script de Python y agregue lo siguiente:

from pyomo.environ import SolverFactory
solver = SolverFactory('glpk')

Para confirmar que tanto Pyomo como su solucionador estén instalados correctamente, vamos a resolver un problema de prueba simple.

Problema de prueba: programa lineal simple.

Objetivo: Minimizar Z=x+y

Subject to:

• x + 2y ≥ 4
• x – y ≤ 1
• x ≥ 0
• \y ≥ 0\

Este problema se refiere a encontrar el valor más pequeño posible de Z, que es la suma de dos variables, x y y. Sin embargo, x y y deben cumplir ciertas condiciones.

Primero, cuando añades x y dos veces y, el resultado debe ser de al menos 4. Segundo, x menos y debe ser menor o igual a 1. Finalmente, ambos x y y deben ser cero o números positivos (no pueden ser negativos). 

El objetivo es encontrar valores de x y y que satisfagan estas condiciones mientras se hace Z lo más pequeño posible.

Implementando utilizando Pyomo:

import pyomo.environ as pyo
# Crear un modelo
model = pyo.ConcreteModel()
# Definir variables
model.x = pyo.Var(within=pyo.NonNegativeReals)
model.y = pyo.Var(within=pyo.NonNegativeReals)
# Definir objetivo
model.obj = pyo.Objective(expr=model.x + model.y, sense=pyo.minimize)
# Definir restricciones
model.con1 = pyo.Constraint(expr=model.x + 2 * model.y >= 4)
model.con2 = pyo.Constraint(expr=model.x - model.y <= 1)
# Seleccionar resolutor
solver = pyo.SolverFactory('glpk')
# Resolver el problema
result = solver.solve(model)
# Mostrar resultados
print('Status:', result.solver.status)
print('Termination Condition:', result.solver.termination_condition)
print('Optimal x:', pyo.value(model.x))
print('Optimal y:', pyo.value(model.y))
print('Optimal Objective:', pyo.value(model.obj))

Si todo funciona correctamente, la salida esperada será:

Status: ok
Termination Condition: optimal
Optimal x: 0.0
Optimal y: 2.0
Optimal Objective: 2.0

Vamos a recorrer el código anterior: primero, define dos variables, x y y, que solo pueden tomar valores no negativos. El objetivo del modelo es minimizar la suma de x y y (x + y). El código define el solucionador como glpk para encontrar los valores óptimos de x y y que satisfacen estas restricciones mientras minimizan el objetivo.

Después de ejecutar el código, observamos que los valores óptimos para las variables son x = 0.0 y y = 2.0, los cuales minimizan la función objetivo Z = x + y. Por tanto, el valor mínimo de la función objetivo es 2.0, lo que satisface las restricciones dadas.

Conocimientos básicos de modelado con Pyomo

Es necesario comprender cómo definir los componentes básicos de un modelo de optimización en Pyomo para configurar y resolver problemas de optimización efectivamente.

1. Definición de variables

Las variables representan las decisiones que se deben tomar en un problema de optimización. En Pyomo, las variables son las cantidades que el solucionador ajustará para optimizar la función objetivo mientras se satisfacen todas las restricciones.

Variables escalares

Una variable escalar es una sola variable que no se indexa sobre ningún conjunto. Para definir una variable escalar en Pyomo, se utiliza la clase Var del módulo pyomo.environ.

from pyomo.environ import Var
model.x = Var()

Primero importamos Var y creamos una variable x usando Var(). Esta variable no tiene límites especificados, lo que significa que puede tomar cualquier valor real, a menos que se restrinja de otra manera en el modelo.

Agregar límites

Puede restringir los valores que puede tomar una variable especificando límites. Los límites se definen como una tupla (lower_bound, upper_bound):

from pyomo.environ import Var
model.x = Var(bounds=(0, None))

Especificar dominios

Pyomo proporciona dominios predefinidos que puede utilizar para especificar el tipo de valores que puede tomar una variable, como NonNegativeReals, Integers o Binary: 

from pyomo.environ import Var, NonNegativeReals
model.x = Var(domain=NonNegativeReals)

Variables indexadas

Cuando se trata de varias variables similares en naturaleza, como por ejemplo variables que representan diferentes periodos de tiempo o elementos, es eficiente utilizar variables indexadas. Las variables indexadas son variables que se definen sobre un conjunto.

import pyomo.environ as pyo
model.I = pyo.Set(initialize=[1, 2, 3])
model.y = pyo.Var(model.I, domain=pyo.NonNegativeReals)

Supongamos que estás modelando las cantidades de producción para tres productos. Puedes definir:

model.Products = pyo.Set(initialize=['A', 'B', 'C'])
model.production = pyo.Var(model.Products, domain=pyo.NonNegativeReals)

Ahora, model.production['A'], model.production['B'], y model.production['C'] representan las cantidades de producción para los productos A, B y C respectivamente.

2. Definición de objetivos

La función objetivo es lo que estamos tratando de optimizar (minimizar o maximizar). Define el objetivo del modelo, como minimizar costos o maximizar beneficios, y generalmente se expresa como una ecuación matemática que involucra a las variables de decisión.

Estos se definen usando la clase Objective:

from pyomo.environ import ConcreteModel, Var, Objective, minimize, maximize, NonNegativeReals
# Crear un modelo
model = ConcreteModel()
# Definir variables
model.x = Var(within=NonNegativeReals)
model.y = Var(within=NonNegativeReals)
# Minimización (costo)
model.cost = Objective(expr=2 * model.x + 3 * model.y, sense=minimize)
# Cuando la maximización de beneficio - (solo puede haber un objetivo a la vez)
# model.profit = Objective(expr=5 * model.x + 4 * model.y, sense=maximize)

3. Agregar restricciones

Las restricciones definen las limitaciones o requisitos del problema:

from pyomo.environ import Constraint
model.con1 = Constraint(expr=model.x + model.y >= 10)

El ejemplo anterior define una restricción en el modelo de Pyomo utilizando la clase Constraint. La restricción model.con1 especifica que la suma de las variables x y y debe ser mayor o igual que 10.

4. Parametrización de modelos

Los parámetros son valores fijos utilizados en el modelo para representar cantidades conocidas o constantes que no cambian durante el proceso de optimización.

Ayudan a definir las relaciones entre variables y restricciones, proporcionando estructura al modelo al incorporar datos o asumidos de la realidad:

from pyomo.environ import Param
model.p = Param(initialize=5)

El código de arriba define un parámetro p en el modelo de Pyomo utilizando la clase Param y lo inicializa con un valor fijo de 5. El parámetro p ahora puede ser utilizado en el modelo para representar un valor constante que no cambia durante el proceso de optimización.

Ahora, trabajemos en un problema de optimización punta a punta!

Ejemplo de Punta a Punta de Pyomo

Veamos un ejemplo de punta a punta de cómo resolver un problema de optimización utilizando Pyomo. Modelaremos una situación real donde una fábrica produce dos productos, y el objetivo es maximizar el beneficio considerando las restricciones de tiempo de la máquina.

1. Declaración del problema

Una fábrica produce dos productos, P1 y P2. La ganancia por unidad es:

• P1: $40
• P2: $50

Tiempo de máquina disponible:

• Máquina A: 100 horas
• Máquina B: 80 horas
• Máquina C: 90 horas

Tiempo requerido por unidad:

Objetivo: Maximizar beneficio.

Variables de decisión:

• x₁: Unidades de P1 para producir.
• x₂: Unidades de P2 para producir.

2. Formulación matemática

Función objetivo:

Maximizar Z = 40x₁ + 50x₂

Restricciones:

1. Capacidad de la máquina A: 1x₁ + 2x₂ ≤ 100
2. Capacidad de la máquina B: 2x₁ + 1x₂ ≤ 80
3. Capacidad de la máquina C: 3x₂ ≤ 90
4. No negatividad: x₁, x₂ ≥ 0

3. Implementación

Basándonos en el objetivo y las restricciones del problema, aquí tenemos el código de Python para modelarlo, nuevamente, utilizando GLPK.

Paso 1: Importar bibliotecas
import pyomo.environ as pyo
Paso 2: Crear un modelo concreto
model = pyo.ConcreteModel()
Paso 3: Definir variables de decisión (unidades de P1 y P2 para producir)
model.x1 = pyo.Var(within=pyo.NonNegativeReals)
model.x2 = pyo.Var(within=pyo.NonNegativeReals)
Paso 4: Definir la función objetivo (maximizar beneficio)
model.profit = pyo.Objective(expr=40 * model.x1 + 50 * model.x2, sense=pyo.maximize)
Paso 5: Definir restricciones
Restricción de capacidad de la máquina A: 1x1 + 2x2 <= 100
model.machine_a = pyo.Constraint(expr=1 * model.x1 + 2 * model.x2 <= 100)
Restricción de capacidad de la máquina B: 2x1 + 1x2 <= 80
model.machine_b = pyo.Constraint(expr=2 * model.x1 + 1 * model.x2 <= 80)
Restricción de capacidad de la máquina C: 3x2 <= 90
model.machine_c = pyo.Constraint(expr=3 * model.x2 <= 90)
Paso 6: Resolver el modelo usando el solucionador GLPK
solver = pyo.SolverFactory('glpk')
result = solver.solve(model)
Paso 7: Analizar resultados
Mostrar el estado del solucionador y las condiciones de finalización
print('Solver Status:', result.solver.status)
print('Termination Condition:', result.solver.termination_condition)
Obtener y mostrar los valores óptimos para x1, x2, y el máximo beneficio
x1_opt = pyo.value(model.x1)
x2_opt = pyo.value(model.x2)
profit_opt = pyo.value(model.profit)
print(f'Optimal production of P1 (x1): {x1_opt}')
print(f'Optimal production of P2 (x2): {x2_opt}')
print(f'Maximum Profit: ${profit_opt}')

Salida:

>>> Solver Status: ok
>>> Termination Condition: optimal
>>> Optimal production of P1 (x1): 25.0
>>> Optimal production of P2 (x2): 30.0
>>> Maximum Profit: $2500.0

En el código de arriba, definimos un modelo de optimización lineal para maximizar el beneficio de la producción de dos productos (P1 y P2). La función objetivo está configurada para maximizar el beneficio, con cada unidad de P1 aportando $40 y cada unidad de P2 aportando $50.

Imponemos tres restricciones que representan los límites de tiempo de las máquinas para las Máquinas A, B, y C.

Finalmente, utilizamos el solucionador GLPK para resolver el problema.

La respuesta final es producir 25 unidades de P1 y 30 unidades de P2 donde nuestro máximo beneficio será de $2,500.

Funciones avanzadas en Pyomo

En la sección anterior, vimos lo fácil que resulta implementar un problema de optimización completo con Pyomo. Sin embargo, la mayoría de los problemas de la vida real no son fáciles de resolver.

En esta sección, presento algunas funciones avanzadas que puede utilizar para resolver escenarios más complejos.

1. Optimización no lineal

La optimización no lineal minimiza o maximiza una función objetivo no lineal sometida a restricciones no lineales. Veamos un ejemplo en el que minimizamos la suma de las diferencias cuadradas sometidas a una restricción circular.

Estatuto del problema

Minimizar el objetivo: Z = (x – 1)² + (y – 2)²

Subject to:

• x² + y² ≤ 4
• x, y ≥ 0

En Pyomo, podemos definir las variables de decisión x y y con límites de 0 para asegurar que no sean negativas. La función objetivo se escribe como la suma de las diferencias cuadradas con respecto a puntos específicos, y la restricción garantiza que la solución se encuentre dentro de un círculo de radio 2.

En este caso, el solucionador IPOPT es adecuado por su capacidad para resolver problemas de optimización no lineal:

import pyomo.environ as pyo
model = pyo.ConcreteModel()
# Defina variables con límites inferiores
model.x = pyo.Var(bounds=(0, None))
model.y = pyo.Var(bounds=(0, None))
# Función objetivo: minimizar (x - 1)² + (y - 2)²
model.obj = pyo.Objective(expr=(model.x - 1)**2 + (model.y - 2)**2, sense=pyo.minimize)
# Restricción: x² + y² ≤ 4 (circulo de radio 2)
model.circle = pyo.Constraint(expr=model.x**2 + model.y**2 <= 4)
solver = pyo.SolverFactory('ipopt')
result = solver.solve(model)
print('Optimal x:', pyo.value(model.x))
print('Optimal y:', pyo.value(model.y))
print('Minimum Z:', pyo.value(model.obj))

2. Programación entera mixta (MIP)

La programación entera mixta se utiliza cuando algunas variables de decisión son enteras (a menudo binarias) y otras son continuas. Es valioso para problemas de toma de decisiones como la ubicación de instalaciones y la planificación de producción.

Ecuación del problema

Una empresa debe decidir si abrir almacenes en ubicaciones A, B, y C. El objetivo es minimizar el costo total, que incluye los costos fijos de apertura de almacenes y los costos de transporte.

Comencemos inicializando los datos, que incluyen los costos fijos de apertura de almacenes, los costos de transporte, los límites de capacidad y la demanda total:

locations = ['A', 'B', 'C']
FixedCost = {'A': 1000, 'B': 1200, 'C': 1500}
TransportCost = {'A': 5, 'B': 4, 'C': 6}
Capacity = {'A': 100, 'B': 80, 'C': 90}
Demand = 150
model = pyo.ConcreteModel()
# Variable binaria: 1 si el almacén está abierto, 0 en caso contrario
model.y = pyo.Var(locations, domain=pyo.Binary)
# Variable continua: cantidad de bienes transportados
model.x = pyo.Var(locations, domain=pyo.NonNegativeReals)
model.cost = pyo.Objective(
    expr=sum(FixedCost[i] * model.y[i] + TransportCost[i] * model.x[i] for i in locations),
    sense=pyo.minimize
)
# Restricción de demanda
model.demand = pyo.Constraint(expr=sum(model.x[i] for i in locations) >= Demand)
# Restricciones de capacidad
def capacity_rule(model, i):
    return model.x[i] <= Capacity[i] * model.y[i]
model.capacity = pyo.Constraint(locations, rule=capacity_rule)
solver = pyo.SolverFactory('cbc')
result = solver.solve(model)
for i in locations:
    print(f"Warehouse {i}: Open={pyo.value(model.y[i])}, Transported={pyo.value(model.x[i])}")
print('Minimum Total Cost:', pyo.value(model.cost))

El modelo incluye dos tipos de variables de decisión: una variable binaria y que representa si un almacén está abierto (1 si está abierto, 0 en caso contrario), y una variable continua x que representa la cantidad de bienes transportados desde cada almacén.

La función objetivo suma los costos fijos y de transporte de cada almacén y minimiza el total. Las restricciones aseguran que la cantidad total de mercancías transportadas cumple con la demanda y que no se excede la capacidad de cada almacén si está abierto.

3. Manejar múltiples objetivos

A veces, los problemas de optimización implican múltiples objetivos que pueden entrar en conflicto, como maximizar el beneficio mientras se minimiza el impacto ambiental. Una aproximación común es el método de suma ponderada, donde cada objetivo se asigna una ponderación para equilibrar su importancia.

Ecuación del problema

Nos proponemos maximizar el beneficio mientras minimizamos el impacto ambiental:

• Beneficio:Z₁ = 3x + 5y
• Impacto ambiental:Z₂ = 2x + y

Podemos combinar estos objetivos utilizando ponderaciones w1=0.6, w2=0.4, donde el objetivo total se convierte en una suma ponderada:

w1 = 0.6
w2 = 0.4
model.obj = pyo.Objective(
    expr=w1 * (3 * model.x + 5 * model.y) - w2 * (2 * model.x + model.y),
    sense=pyo.maximize
)

En este objetivo combinado, maximizamos el beneficio mientras minimizamos el impacto ambiental mediante la ajustación de las ponderaciones.

4. Usando fuentes de datos externas

Cuando se manejan grandes conjuntos de datos, importar datos de fuentes externas, como archivos CSV, es a menudo útil. Pyomo funciona bien con Pandas para leer y utilizar datos externos.

Podemos leer un archivo CSV utilizando Pandas y usar los datos para inicializar conjuntos y parámetros en nuestro modelo:

import pandas as pd
data = pd.read_csv('parameters.csv')
# Definir conjunto a partir de datos CSV
model.I = pyo.Set(initialize=data['index'].unique())
# Definir parámetro inicializado a partir de datos CSV
param_dict = data.set_index('index')['value'].to_dict()
model.param = pyo.Param(model.I, initialize=param_dict)

Consejos y Mejores Prácticas para Usar Pyomo

Al trabajar con Pyomo, es importante mantener tus modelos eficientes, bien documentados y fáciles de depurar.

1. Depuración y resolución de problemas

Mientras se construyen modelos de optimización en Pyomo, es común encontrar problemas como soluciones infeasibles, fallas del resolutorio o resultados incorrectos. Aquí tienes algunas mejores prácticas para la depuración:

• Constraintes de verificación: Revise sus restricciones si su modelo no está generando una solución factible. Restricciones demasiado estrictas pueden hacer que un problema no sea factible. Use el método .display() de Pyomo para imprimir los valores de las variables y restricciones para verificar que están funcionando como se espera.
• Salida del solucionador: Active los registros detallados del solucionador pasando tee=True cuando llame al método solve(). Esto puede proporcionar una idea de dónde podría tener dificultades el solucionador, como variables sin límites o infeasibilidad.
• Pruebe primero modelos sencillos: Cuando se trata de modelos complejos, pruebe una versión simplificada. Esto puede ayudar a aislar problemas potenciales sin la sobrecarga de un modelo completamente especificado.

La depuración es mucho más fácil si se aborda de manera sistemática, analizando restricciones, la función objetivo y la retroalimentación del solucionador.

2. Eficiencia en modelado

Los problemas de optimización pueden resultar computacionalmente costosos a medida que aumenta el tamaño del modelo. Para garantizar un modelado eficiente, considere las siguientes recomendaciones:

• Usar esparcidad: Evite recorrer índices innecesarios al definir restricciones o objetivos. Aprovechar la esparcidad en su problema reduce el tiempo de cómputo.
• Variables binarias versus continuas: En cuanto sea posible, reduzca el número de variables binarias o enteras. Las variables continuas son más fáciles de manejar para los solucionadores, lo que conduce a soluciones más rápidas.
• Fórmulación de restricciones: Mantenga las restricciones lo más simples posibles, tanto en forma matemática como en su implementación. Evite las no linealidades innecesarias y descomponga las restricciones complejas en pequeñas y manejables.

Los modelos eficientes resuelven más rápido y son más fáciles de depurar y mantener.

3. Documentación y mantenimiento

El mantenimiento de modelos de Pyomo bien documentados es una buena práctica para un uso a largo plazo y colaboración. Una buena documentación también facilita revisar y actualizar los modelos con el tiempo:

• Usar comentarios en línea: Siempre añade comentarios para explicar el propósito de las variables, las restricciones y la función objetivo. Esto es especialmente importante en modelos de optimización donde la lógica podría no ser inmediatamente obvia.
• Modularizar su código: Divide su modelo en secciones lógicas o incluso funciones separadas. Este enfoque modular puede mejorar la lectura y hacer que el depurado y la modificación de partes específicas del modelo sean más fáciles.
• Controlar cambios del modelo: Mainte una historia de versiones de su modelo, especialmente si evoluciona. Utilice herramientas de control de versiones como Git para rastrear cambios y asegurar que cualquier actualización o mejora pueda ser trazada atrás.

Un documentación adecuada y un código estructurado harán que sus modelos de Pyomo sean más accesibles para futuros colaboradores y más fáciles de escalar o modificar a medida que evolucionen sus requisitos.

Conclusión

Pyomo es una herramienta poderosa y flexible para construir y resolver modelos de optimización en Python. A lo largo de este tutorial, exploramos cómo Pyomo permite a los usuarios modelar diversos problemas de optimización, desde la programación lineal hasta la no lineal y la programación de enteros mixtos.

Con su sintaxis amigable y su integración con solucionadores, Pyomo hace que la formulación y resolución de problemas de optimización reales sea accesible tanto para principiantes como para usuarios avanzados.

Si estás interesado en aprender más sobre la resolución de problemas reales con optimización, echa un vistazo al curso gratuito Introducción a la Optimización en Python en DataCamp!