Tutorial Pyomo (II) — Parámetros, Sets y variables indexadas

Modela problemas de optimización complejos con Pyomo. Parámetros, Sets y variables indexadas en Pyomo, útiles en modelos temporales

En el Tutorial Pyomo (I) — Introducción a la programación lineal en Python vimos una breve introducción a Pyomo. A través del artículo modelamos y resolvimos nuestro primer modelo de programación lineal. En este artículo, profundizaremos cómo enfocar el modelado de problemas más complejos.

Todo el código del proyecto puede encontrarse en:

RikiSot/pyomo-tutorial: Tutorial series to understand linear programming and Pyomo basics (github.com)

Una de las mayores ventajas de Pyomo es la facilidad con la que se pueden modelar problemas a priori difíciles gracias al uso de los componentes de su librería AML. Sin estos componentes, sería muy complicado ya no solo modelar cierto tipo de problemas, sino depurar y verificar cada variable durante el proceso de validación del modelo.

Qué mejor que un ejemplo para comprender todo esto. Hasta ahora, los problemas a los que nos hemos enfrentado se componían de variables sencillas, variables que estaban definidas por un valor escalar a la hora de resolver el problema.

from pyomo.environ import *

model = ConcreteModel()

# Definición de variables
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

Esto es útil cuando cada variable representa un único resultado y el problema es sencillo. Vamos a usar este ejemplo para formular un problema sencillo y a continuación iremos complicando el modelo para comprender mejor la nueva problemática.

Planteamiento del problema

Imaginemos que queremos dimensionar un rebaño de ovejas. Tenemos ovejas de dos tipos: churras y merinas. Como bien dice el refranero español:

“No hay que mezclar churras con merinas”

En algún momento de la historia alguien decidió dedicarle un refrán a estas ovejas

Por lo que seguimos su consejo y asignamos a X la cantidad de churras y a Y la cantidad de merinas.

Disclaimer: no tengo ni idea de ganadería así que seguramente el problema que voy a plantear carezca de sentido real.

Vamos a plantear el siguiente problema de optimización lineal:

Función objetivo

Minimizar el coste total del rebaño:

Parámetros

• C_ch es el coste de cada churra y Cₘ el de cada merina
• Cada churra produce Prod_ch kg de lana
• Cada merina produce Prodₘ kg de carne
• Capacidad máxima de animales en nuestras instalaciones: CAP
• Tamaño mínimo del rebaño: Oₘᵢₙ
• Necesidades de demanda de lana: Dₗ
• Necesidades de demanda de carne: D_c

Restricciones

• El número total de ovejas en el rebaño tiene que ser mayor de Oₘᵢₙ (si no no compensa hacerse pastor!)
• Las merinas tienen una carne de mayor calidad mientras que las churras ofrecen una mejor lana. La demanda de nuestra zona es mayor para la carne con lo que tenemos que asegurarnos al menos el doble de merinas que de churras
• Tenemos capacidad para alojar como máximo CAP animales
• Tenemos que cubrir las necesidades de demanda tanto de lana como de carne

Con esta información, podemos modelar el problema en Pyomo:

from pyomo.environ import *

model = ConcreteModel()

# Parámetros
model.C_ch = Param(within=NonNegativeReals, default=10)
model.C_m = Param(domain=NonNegativeReals, default=15)
model.Prod_ch = Param(domain=NonNegativeReals, default=2)
model.Prod_m = Param(domain=NonNegativeReals, default=3)
model.CAP = Param(domain=NonNegativeIntegers, default=100)
model.O_min = Param(domain=NonNegativeIntegers, default=10)
model.D_l = Param(domain=NonNegativeReals, default=20)
model.D_c = Param(domain=NonNegativeReals, default=30)

# Variables
model.X = Var(domain=NonNegativeIntegers)
model.Y = Var(domain=NonNegativeIntegers)

# Restricciones
model.res1 = Constraint(expr=model.X + model.Y >= model.O_min)
model.res2 = Constraint(expr=model.Y >= 2*model.X)
model.res3 = Constraint(expr=model.X + model.Y <= model.CAP)
model.res4 = Constraint(expr=model.Prod_ch * model.X >= model.D_l)
model.res5 = Constraint(expr=model.Prod_m * model.Y >= model.D_c)

# Función objetivo
model.obj = Objective(expr=model.C_ch * model.X + model.C_m * model.Y, sense=minimize)

Podemos apreciar un componente nuevo que no aparecía en el anterior tutorial: el componente Param. Este componente de Pyomo sirve para definir constantes numéricas, es decir, valores conocidos antes de resolver el problema. Al no ser variables, pueden multiplicar componentes de tipo Var() sin romper la condición de linealidad.

En este caso, al ser parámetros escalares, para asignarles un valor podemos especificárselo a través del parámetro default, que es el valor que toman si no se especifica nada más. Más adelante veremos que para estructuras de datos más complicadas se procede de otra manera.

Realmente, no tienen mucha aplicación en este tipo de problema ya que podríamos conseguir el mismo resultado simplemente con una variable en como C_ch = 10 pero así vamos introduciendo el concepto. Además, siempre es buena práctica que todo lo referente al modelo quede reflejado dentro de su estructura, ya que nos facilitará la vida a la hora de depurar.

Complicando el modelo

Vale bien, este era fácil verdad? Como nos ha ido tan bien planificando nuestro rebaño, queremos programar la compra de churras y merinas de los años venideros. El resultado obtenido hasta ahora nos vale para un año, cómo resolvemos para por ejemplo, los siguientes 10 años?

— Fácil no? Resolvemos el mismo modelo 10 veces cambiando las condiciones de contorno y ya.’

— Gracias, voz interior, pero buscaremos una solución un poco más práctica

Obviamente no vamos a solucionar el modelo 10 veces revisando las condiciones cada vez. Vamos a introducir un concepto nuevo, el componente Set.

En el ejemplo anterior, cada uno de los componentes del modelo era escalar. Es decir, tenían un único valor, no vectores o matrices. Esto también incluye a las restricciones. En problemas complicados, es común tener vectores de variables y restricciones. Esto es gestionado por Pyomo a través de los componentes indexados.

Aquí es donde entran los sets. Un Set es una colección de datos, numéricos o no (por ejemplo, strings) que se usa para especificar los índices válidos de los componentes indexados. Hay varias formas de crear un Set en Pyomo:

- Set: Componente genérico para declarar sets a partir de datos
- RangeSet: Genera un set a partir de un rango de números (1, 2, 3…)

Para nuestro caso, vamos a generar un índice de tiempos para indicar el año correspondiente a cada variable:

model.T = RangeSet(10)

Con esto creamos un índice de tiempos que varía desde 0 a 9, representando el año en el índice. También podemos asignar el set declarando previamente una estructura de datos. En el parámetro initialize podemos especificar un `iterable`, como una lista, un set o una tupla:

years = [x for x in range(0, 9)]
model.T = Set(initialize=years)

Un pequeño detalle: cuando Pyomo instancia un RangeSet, utiliza índices empezando desde 1. Esto se debe a que internamente utiliza parejas de key-values para acceder a los índices, lo que ocasiona problemas al intentar pasar una lista (ya que el índice 0 no existe en el diccionario). Esto se puede corregir especificando el índice 0 como start en el RangeSet o empleando diccionarios para inicializar los sets.

Con el set ya creado podemos refactorizar nuestras variables para que sean indexadas. Para crear variables indexadas, hay que pasar el set (o sets) en la declaración de la variable:

model.x = Var(model.T, domain=NonNegative)
model.y = Var(model.T, domain=NonNegative)

Con esto conseguimos que las variables X e Y tengan 10 elementos en su interior, un valor para cada año (x1, x2, x3…).

Ahora hay que reformular las restricciones para tener en cuenta el índice temporal. En este caso no es tan inmediato. Por ejemplo, para reescribir la primera restricción tenemos que pasar el set correspondiente como primer parámetro y adaptar la expresión:

# Method 1
def restriction1_rule(m, t):
  return m.x[t] + m.y[t] >= m.O_min

model.res1 = Constraint(model.T, rule=restriction1_rule)

# Method 2
@model.Constraint(model.T)
def restriction1(m, t):
  return m.x[t] + m.y[t] >= m.O_min

He de decir que soy más fan del método 2, los decoradores ayudan a que todo quede más limpio. El primer argumento de la restricción siempre es el modelo, independientemente de que esta sea indexada o no. Los siguientes argumentos son los sets que corresponden a los índices de las variables con las que queremos trabajar. En nuestro caso t es el índice temporal.

Con este método generaremos una restricción distinta para cada índice temporal. Por simplificar las cosas, tal y como está ahora la restricción representa la cantidad mínima de ovejas para cada año, sin tener en cuenta el estado anterior, pero este no es un caso realista.

¿Qué pasa por ejemplo si compro ovejas de más durante un año? Obviamente las seguimos conservando para el siguiente periodo, en el que no partimos de cero. Podríamos reformular la restricción teniendo en cuenta estas relaciones entre distintos índices temporales:

@model.Constraint(model.T)
def restriction1(m, t):
  if t == 0:
    o_init = 0
  else:
    o_init = m.x[t-1] + m.y[t-1]
  return o_init + m.x[t] + m.y[t] >= m.O_min

Vamos, que se puede complicar todo lo que queramos. Podemos añadir lógica según el índice temporal, relacionar periodos, etc.

También podríamos añadir la componente temporal a algún parámetro para darle un poco más de realismo al problema. Por ejemplo, que la demanda de productos sea un array de valores (dependiendo del año)

Con esto ya podemos reformular el modelo completo. Los detalles del nuevo modelo pueden consultarse con model.pprint(). Hay que cambiar ligeramente la forma de acceder a la solución ya que nuestras variables ahora son arrays:

Solución óptima: 
X: [10.0, 15.0, 10.0, 20.0, 10.0, 13.0, 5.0, 8.0, 6.0, 10.0]
Y: [20.0, 30.0, 20.0, 40.0, 20.0, 26.0, 10.0, 16.0, 12.0, 20.0]
Coste total: 4280.0

Resumen

En este artículo, hemos profundizado en el modelado de problemas de optimización lineal más complejos con Pyomo. Para ello, hemos introducido dos nuevos componentes:

1. Parámetros: Permiten definir constantes numéricas en el modelo. Son útiles para separar la lógica del problema de los valores específicos, facilitando la lectura y el mantenimiento del código.
2. Sets: Permiten definir conjuntos de datos que se utilizan para indexar variables y restricciones. Esto facilita el modelado de problemas con múltiples escenarios, elementos, periodos de tiempo, etc.

A través de un ejemplo práctico, hemos visto cómo usar estos componentes para:

• Modelar variables indexadas: Las variables ahora pueden tener un valor para cada elemento del set, permitiendo representar diferentes entidades o momentos temporales.
• Añadir lógica temporal a las restricciones: Las restricciones pueden tener en cuenta el estado del sistema en diferentes periodos, lo que permite modelar problemas más realistas.

Siguientes pasos

En el siguiente artículo, exploraremos el modelado de problemas de programación lineal entera (IP) y mixta (MILP) con Pyomo. Estos tipos de problemas son aún más complejos que la programación lineal tradicional, pero pueden usarse para modelar una amplia gama de problemas del mundo real.

Recursos adicionales

• Documentación oficial de Pyomo:
• Ejemplos de modelos Pyomo: https://jckantor.github.io/CBE30338/06.99-Pyomo-Examples.html
• Foro de Pyomo
• Williams, H. P. (2013). Model Building in Mathematical Programming, 5th Edition | Wiley
• Pyomo — Optimization Modeling in Python | SpringerLink