Solucionar problemas de optimización es una de las claves de toda persona que se dedique al mundo de los datos: encontrar rutas óptimas que minimicen el recorrido, maximizar el margen de una producción industrial… Son muchos los casos en los que aplicar modelos de optimización es fundamental. Así pues, en este post te voy a explicar cómo puedes crear modelos de optimización en Python con la librería Pyomo. En mi opinión, Pyomo es una de las librerías de optimización por distintos motivos:
- Sirve para para solucionar problemas con programación lineal y no lineal.
- Permite solucionar problemas de programación no lineal con variables binarias.
- La forma de definir el código es muy Pythonic, lo cual es muy interesante.
- Es una librería muy utilizada y con mucho soporte. Por lo que si tienes cualquier duda, es muy probable que encuentres una solución o alguien pueda ayudarte.
Dicho esto, vamos a empezar con el tutorial de cómo realizar optimización en Python con Pyomo. ¡Empecemos con la instalación!
Cómo instalar Pyomo
La instalación de la librería de optimización Pyomo la podemos dividir en dos partes: por un lado la librería y, por otro lado, el solver.
Instalar la librería pyomo es relativamente sencillo. Simplemente tendremos que abrir una terminal e instalar la librería usando pip, como cualquier otra librería. Ejemplo:
pip insall pyomo
Hecho esto, ahora vamos a la segunda cuestión, instalar el solver que vayamos a utilizar.
Instalación del solver usado por Pyomo
Pyomo acepta el uso de distintos tipos de solvers. Puedes obtener el listado de todos los solvers aceptados ejecutando el siguiente comando (deberás tener pyomo instalado previamente):
pyomo help --solvers
En este sentido, voy a explicar cómo puedes instalar dos de los principales solvers:
CBC
- Debian
sudo apt-get install coinor-cbc coinor-libcbc-dev
- conda
conda install coin-or-cbc
- Mac (suponiendo que tienes breweri instalado)
brew tap coin-or-tools/coinor
brew install coin-or-tools/coinor/cbc
GLPK
- Ubuntu
sudo apt-get install python-glpk
sudo apt-get install glpk-utils
- Mac (suponiendo que tienes Homebrew instalado)
brew install glpk
- Windows: puedes seguir este tutorial.
Perfecto, ahora que ya sabes cómo puedes instalar Pyomo y su solver, sigamos con el tutorial viendo cuáles son los componentes básicos de Pyomo.
Componentes Básicos de Pyomo
Modelo
El modelo es donde definimos las variable objetivo, restricciones y todo lo que tiene que ver con el problema a resolver. Los modelos pueden ser de dos tipos:
- Concretos: son aquellos casos en los que queremos crear un optimizador para solucionar un problema con unos datos concretos. Si queremos obtener los resultados para otros datos, deberemos cambiar el modelo.
- Abstractos: sirve para aquellos casos en los que queremos solucionar un problema sean cuales sean los datos.
Si quieres solucionar un problema de forma puntual, usarás el modelo concreto. Si crees que el problema puede volver a darse en el futuro y te gustaría solucionarlo con otros datos, usarás un modelo abstracto.
Puedes definir un modelo concreto de la siguiente manera:
import pyomo.environ as pyo
model = pyo.ConcreteModel()
Asimismo, puedes definir un modelo abstracto de la siguiente manera:
import pyomo.environ as pyo
model = pyo.AbstractModel()
Asimismo, otra cuestión importante es, una vez hemos definido todos los componentes, resolver el modelo. En este sentido, los modelos concretos y los abstractos se resuelven de formas diferentes.
Por un lado, un modelo concreto lo podemos solucionar directamente aplicando el método solve de la clase SolverFactory, tal como se muestra en el siguiente ejemplo:
SolverFactory('glpk').solve(model).write()
Sin embargo, los modelos abstractos son algo más complejos, ya que el modelo como tal no incluye ningún tipo de dato. Así pues, deberemos:
- Definir los datos que utilizará el modelo.
- Crear una instancia con dichos datos.
- Encontrar el óptimo en dicha instancia.
Veamos un ejemplo:
data = {
None : {
"products": ["Product A", "Product B"],
"labor": {
"Product A": 2,
"Product B": 4
},
"material": {
"Product A": 3,
"Product B": 2
},
"product_profit": {
"Product A": 10,
"Product B": 15
},
"total_labor": {None: 500},
"total_material": {None: 800}
}
}
instance = model.create_instance(data)
optimizer = pyo.SolverFactory("glpk")
optimizer.solve(instance)
instance.display()
Sí no entiendes los datos, no te preocupes, lo entenderás en cuanto lleguemos a los ejercicios. Por ahora lo que más importa es que recuerdes que los modelos abstractos y los concretos se resuelven de forma diferente.
Variable
Este componente sirve para definir las variables cuyos valores el modelo de optimización puede cambiar para solucionar el problema.
Por ejemplo, si quisíeramos obtener el número de unidades óptimas a producir de un producto A y un producto B para maximizar el margen, tendríamos dos opciones:
- Si usamos un modelo concreto, deberíamos crear dos variables, una indicando el número de unidades producidas del producto A y otra indicando el número de unidades producidas del producto B.
- Si usamos un modelo abstracto, deberíamos crear una única variable, la cual indica cuántas unidades deberíamos producir de cada producto.
Asismimo, cuando definimos una variable debemos indicar también qué tipo de valores puede tomar esa variable. Esto lo hacemos pasando distintas funciones al parámetro domain. Las opciones son las siguientes:
NonNegativeReals: incluye valores no negativos, incluyendo el 0.NegativeIntegers:Any: todos los posibles valores.Reals: valores con coma.PositiveReals: valores positivos (excluido el 0) con decimalesNonPositiveReals: valores negativos (incluido el 0) con decimales.NegativeReals: valores negativos (excluido el 0) con decimales.NonNegativeReals: valores positivos (incluido el 0) con decimales.PercentFraction: valores con decimales entre 0 y 1.Integers: cualquier valor entero.PositiveIntegers: valor entero positivo.NonPositiveIntegers: valor entero negativo (incluido el 0).NegativeIntegers: valor entero negativo.NonNegativeIntegers: valor entero no negativo (incluido el 0).Boolean: valores booleanas que pueden definirse comoTrueyFalseo 1 y 0.
Sabiendo esto, veamos cómo definiríamos las variables en el caso del modelo concreto:
import pyomo.environ as pyo
# Define the model
model = pyo.ConcreteModel()
# Define the variables
model.productA = Var(domain = pyo.NonNegativeIntegers)
model.productB = Var(domain = pyo.NonNegativeIntegers)
Perfecto, ahora que conoces más sobre las variables, vayamos al siguiente punto: los sets.
Set
Hemos visto cómo se definen las variables en los casos de los modelos concretos, pero, ¿y en el caso de los modelos abstractos? ¿Cómo sabrá el modelo que tengo 2 productos diferentes a fabricar? Y aunque lo pudiese inferir del resto de datos, ¿cómo sé yo cuál es cada producto? Para solucionar estas dudas es donde entran los sets.
Un set no es más que un objeto que nos permite definir cuántos «miembros» hay y qué nombre recibe cada miembro. Por ejemplo, si queremos definir un modelo abstracto para la fabricación de productos, nosotros de antemano no sabemos cuántos productos hay ni cómo se llaman cada uno de estos productos.
Así pues, crearemos un set de productos. Si, más adelante debemos incluir un parámetro (más adelante verás qué es un parámetro) por producto, como por ejemplo, el consumo de materias primas, podemos crear un parámetro en base al set. Así nos aseguraremos de que debe haber un parámetro definido para cada miembro del set.
En este sentido, un set puede tener una cantidad de dimensiones ya definidas (por ejemplo, indico que solo puede aceptar 2 productos), o indefinidas. En el primer caso, el modelo abstracto solo funcionaría con dos productos mientras que, con el segundo, el optimizador funcionaría con una cantidad ilimitada de productos (2, 3, 10, etc.).
Si no te ha quedado 100% claro qué es un set, no te preocupes, creo que se entiende mejor cuando veamos la práctica. De momento, lo que me interesa es que te quedes con que el set nos permite definir cuántos miembros hay y cómo se llama cada uno de ellos. Asimismo, quédate también con que un set se puede definir de la siguiente manera:
import pyomo.environ as pyo
model = pyo.AbstractModel()
# Define the variables
model.products = pyo.Set()
Ahora que hemos visto los sets, sigamos con nuestro tutorial de optimización en Python con Pyomo, viendo qué son los parámetros:
Parámetros
Los parámetros permiten pasar a Pyomo los datos que va a necesitar para solucionar el problema. Por ejemplo, si queremos optimizar la producción que tiene uso de cierta cantidad de materia prima, el número de unidades de materia prima que consume cada producto sería un parámetro.
La realidad es que, si trabajas con un modelo concreto, no es estrictamente necesario que los incluyas (aunque sí recomendable), ya que puedes meter esta información directamente a mano. En el caso de los modelos abstractos, en cambio, sí es obligatorio definir los parámetros del modelo.
Asimismo, es interesante remarcar que puede darse el caso de que haya un parámetro (valor) por cada campo del set (o por varios sets), por ejemplo, que cada producto consuma una cantidad de materias primas determinada. En este caso, deberemos vincular el parámetro con el Set, de tal forma que sea obligatorio que haya definido un parámetro por set. Ejemplo:
import pyomo.environ as pyo
model = pyo.AbstractModel()
# Define the variables
model.products = pyo.Set()
# Defin the parameters (data)
model.material = pyo.Param(model.products)
También puede darse el caso de que haya un parámetro para cada combinación de sets. Supongamos por ejemplo que tenemos 3 fábricas que pueden producir 3 productos, cada una con un coste de transporte diferente. Tanto la fábrica como el producto cada uno de ellos sería un set. Así pues, nuestro parámetro coste de transporte tendrá que tener un valor para cada combinación de fábrica y producto. Esto lo podemos definir de la siguiente manera:
import pyomo.environ as pyo
model = pyo.AbstractModel()
# Define the sets of values
model.factories = pyo.Set()
model.clients = pyo.Set()
# Define the input data
model.transportation_cost = pyo.Param(model.factories, model.clients)
Sí, sé que son muchas cosas las que hay que tener en cuenta en el optimizador, pero no te preocupes, solo nos quedan dos cosas más: los objetivos y los constraints. Veamos qué son los objetivos.
Objetivos
Los obejtivos se creanc con el objeto pyo.Objective y nos permiten definir qué es lo que queremos conseguir maximizar o minimizar utilizando las variables y los parámetros disponibles. Si queremos maximizar el beneficio, por ejemplo, el objetivo debe ser la fórmula para calcular dicho beneficio mientras que, si queremos minimizar el tiempo de transporte, el objetivo debe ser la fórmula para calcular el tiempo total de transporte.
Así pues, poodemos ver como a simple vista se distinguen dos grandes componentes del objetivo: la función objetivo y el tipo de operación que debemos aplicar (maximizar o minimizar).
Para definir si queremos maximizar una función objetivo, simplemente deberemos pasar la función pyo.maximize al parámetro sense y pyo.minimize si queremos minimizar la función.
Por otro lado, la función obejtivo hay dos formas de definirla. Por un lado, podemos definir la función como una función en una única línea, definiendo la función en el parámetro expr. Esta forma es más común en los modelos de optimización concretos y que son relativamente simples.
Sin embargo, cuando el modelo es abstracto o siendo concreto ya es algo más complejo, se tiende a crear una función que recibe el modelo que hemos creado y calcular con ello el valor objetivo y, posteriormente, pasar dicha función al parámetro rule .
Una vez más, creo que el funcionamiento se entenderá mejor cuando, más adelante, veamos unos ejemplos.
Por último, vamos a terminar la sección de componentes básicos de este tutorial de optimización en Python con Pyomo viendo cómo definir los constraints o las limitaciones de nuestro modelo.
Constraints
Los constraints son aquellas limitaciones que debemos aplicar en los proyectos de optimización. Por ejemplo, si queremos fabricar los productos A y B y para ello contamos con la matería prima X, hemos de tener en cuenta que no contamos con una cantidad infinita de materia prima, sino que esta será limitada (digamos, 100). Aquí es donde entran en juego los constraints.
Parecido al caso de los objetivos, hay dos formas de definir los constraints:
-
Definiendo el constraint como una función de una sola línea. Esto es aplicable en los modelos concretos con pocos parámetros. Para ello, debemos definir la función en el parámetro
expr.Ejemplo:
model.labor = Constraint(expr = ((model.productA * 2) + (model.productB*4)) <= 500) -
Crear una función que tomando el modelo, calcule el constraint. Esta función deberemos pasarla al parámetro
rule. Ejemplo:
def total_labor_rule(model):
return (0, sum(model.labor[prod] * model.production[prod] for prod in model.products), model.total_labor)
model.labor_consumed = pyo.Constraint(rule = total_labor_rule)
Tal como vemos, las funciones de constraint pueden devolver o bien un booleano que compruebe si la condición se cumple o no (como en el primer ejemplo), o bien una tupla con tres valores: el límite inferior, el valor obtenido y el límite superior, tal como he hecho en el segundo ejemplo.
Asimismo, puede darse el caso en el que contemos queramos aplicar una restricción por set. En este caso, podemos pasar el set a la hora de definir el constraint y Pyomo se encargará de iterar para cada uno de los valores del set. Si no lo entiendes, no te preocupes, lo veremos en el segundo ejemplo más adelante.
Perfecto, con esto ya conoces cuáles son los componentes básicos de un modelo de optimización en Python con Pyomo. Sé que seguramente tengas dudas, pero no te preocupes, ahora entraremos en distintos ejemplos de optimización usando Pyomo para que puedas entender mejor cómo funcionan. ¿Te suena interesante? ¡Vamos con ello!
Ejemplo 1: Optimización de Fabricación con Pyomo
Contexto del problema a resolver
Supongamos que tenemos una empresa que fabrica dos productos: Producto A y Producto B. La producción de cada producto requiere una cierta cantidad de mano de obra y materias primas. En este sentido, la empresa dispone de 500 unidades de mano de obra y 800 unidades de materia prima. Además, cada producto ofrece diferente ganancia: la ganancia por unidad del Producto A es de $10 y la ganancia por unidad del Producto B es de $15.
La empresa quiere maximizar su beneficio total. Los requisitos de producción y la disponibilidad de recursos son los siguientes:
- El producto A requiere 2 unidades de mano de obra y 3 unidades de materias primas por unidad.
- El producto B requiere 4 unidades de mano de obra y 2 unidades de materias primas por unidad.
Con esta información, deberemos montar un sistema de optimización que nos indique cuánto debemos de producir de cada producto para maximizar el beneficio.
Dicho esto, veámos cómo enfocarlo tanto con un modelos concreto como con un modelo abstracto.
Solución con un modelo concreto
En mi opinión, la mejor forma de solucionar un problema de optimización (ya sea con Pyomo o cualquier otro framework) es analizando poco a poco cada una de los elementos que incluye en framework.
En este caso, hemos dicho que queremos solucionar el problema con un modelo concreto, así pues, lo primero de todo es crear el modelo:
import pyomo.environ as pyo
model = pyo.ConcreteModel()
Ahora que tenemos el modelo creado, vamos al siguiente paso: definir las variables. Si recuerdas, las variables son los campos que el optimizador puede modificar para llegar al resultado.
En nuestro caso, contamos con dos variables: el número de unidades fabricadas de producto A y el número de unidades fabricadas de producto B.
En ambos casos, dicha variable deberá ser no negativa (no puedes producir en negativo, pero sí 0) y con valores enteros, puesto que tampoco puedes crear medios productos. Así pues, indicaremos que el domain sea NonNegativeIntegers.
# Define the variables
model.productA = pyo.Var(domain = pyo.NonNegativeIntegers)
model.productB = pyo.Var(domain = pyo.NonNegativeIntegers)
El siguiente paso es definir los parámetros (datos) del modelo. Tal como he comentado previamente, en el caso de los modelos concretos no es algo obligatorio, pero sí recomendable. Nosotros lo haremos para seguir las buenas prácticas.
Así pues, disponemos de los siguientes datos:
- Coste de mano de obra del producto A: 2.
- Coste de mano de obra del producto B: 4.
- Coste de materia prima del producto A: 3.
- Coste de materia prima del producto B: 2.
- Beneficio del producto A: 10$.
- Beneficio del producto B: 15$.
Así pues, vamos a crear dichos parámetros:
model.laborA = pyo.Param(initialize = 2)
model.laborB = pyo.Param(initialize = 4)
model.materialA = pyo.Param(initialize = 3)
model.materialB = pyo.Param(initialize = 2)
model.profitA = pyo.Param(initialize = 10)
model.profitB = pyo.Param(initialize = 15)
Perfecto, lo siguiente sería definir cuál es el objetivo que queremos obtener. En nuestro caso, queremos maximizar el beneficio. El beneficio se calculará como la suma del número de unidades vendidas por producto multiplicado por su margen correspondiente.
Asimismo, vamos a definir dicho objetivo en una sola línea de código (expresión):
# Define the objective
model.profit = pyo.Objective(
expr = ((model.productA * model.profitA) + (model.productB * model.profitB)),
sense = pyo.maximize
)
Aunque pueda parecer que con esto ya estaría, no es así. Nos faltan incluir las limitaciones de nuestro modelo de optimización. Así pues, si recordamos el problema contamos con dos limitaciones:
- Contamos con 500 unidades de mano de obra.
- Contamos con 800 unidades de materia prima.
Así pues, el consumo de estos recursos debe ser inferior a estos valores. Lo incluimos:
# Define the constraints
model.labor = pyo.Constraint(expr = ((model.productA * model.laborA) + (model.productB* model.laborB)) <= 500)
model.material = pyo.Constraint(expr = ((model.productA * model.materialA) + (model.productB * model.materialB)) <= 800)
¡Ya lo tenemos! Por último, vamos a resolver nuestra problema a ver qué resultados nos ofrece:
pyo.SolverFactory('glpk').solve(model).write()
# Print the results
print(f"Profit = {model.profit()}")
print(f"Units of Product A = {model.productA()}")
print(f"Units of Product B = {model.productB()}")
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 2500.0
Upper bound: 2500.0
Number of objectives: 1
Number of constraints: 2
Number of variables: 2
Number of nonzeros: 4
Sense: maximize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 1
Number of created subproblems: 1
Error rc: 0
Time: 0.038468360900878906
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
Profit = 2500.0
Units of Product A = 250.0
Units of Product B = 0.0
Como podemos ver, la solución pasa por fabrican 250 unidades de producto A y ninguna unidad de producto B.
Ahora, vamos a seguir con el tutorial de optimización con Pyomo, solucionando este mismo problema, pero mediante un modelo abstracto. ¡Vamos a ello!
Solución con un modelo abstracto
Para solucionar el problema previamente planteado con un modelo abstracto, lo primero de todo será definir dicho modelo abstracto:
import pyomo.environ as pyo
model = pyo.AbstractModel()
Asimismo, al tratarse de un modelo abstracto, sabemos que hablaremos de productos, pero necesitaremos saber cuál es cada producto. Así pues, crearemos un Set para poder definir cuántos productos hay y cómo se llama cada uno de los productos:
# Define the sets
model.products = pyo.Set()
Por otro lado, tendremos que definir las variables del modelo, es decir, cuántos productos se fabricarán por producto. En este caso, contamos únicamente con una variable (las unidades a fabricar), las cuales dependerán del set. Así pues, procedemos a definir la variable:
# Product production
model.production = pyo.Var(model.products, domain=pyo.NonNegativeReals)
Asimismo, deberemos indicar también los parámetros del modelo. Al tratarse de un modelo abstracto, no sabemos cuántos productos habrá. En cualquier caso, si sabemos que vamos a necesitar los siguientes datos:
- Consumo de materias primas por producto.
- Consumo de unidades laborales por producto.
- Total de unidades de materias primas disponibles.
- Total de unidades laborales disponibles.
- Beneficio que genera cada producto.
Así pues, crearemos dichos parámetros. En caso de que deba haber un valor por producto, lo indicaremos pasando el set:
# Defin the parameters (data)
model.labor = pyo.Param(model.products)
model.material = pyo.Param(model.products)
model.product_profit = pyo.Param(model.products)
model.total_labor = pyo.Param()
model.total_material = pyo.Param()
Por otro lado, tenemos que definir también el objetivo de nuestro modelo. En este caso, el objetivo es maximizar el beneficio. A diferencia del modelo concreto, esta vez crearemoos una función que obtenga el modelo como input y calcule el valor objetivo. Pasaremos dicha función al parámetro rule para que se utilce como objetivo.
Al usar el modelo en la función, lo más común suele ser iterar por el set (o los sets) que hagan falta para calcular el objetivo, tal como realizo en el ejemplo.
# Define the objective function
def total_profit_rule(model):
return sum(model.product_profit[prod] * model.production[prod] for prod in model.products)
model.total_profit = pyo.Objective(
rule = total_profit_rule,
sense = pyo.maximize
)
Genial, ya tenemos el objetivo, por lo que únicamente nos queda definir las restricciones de nuestro modelo. Al igual que con el objetivo, definiremos las restricciones como funciones a parte, las cuales pasaremos al parámetro rule, tal como se muestra a continuación:
# Constraints
def total_labor_rule(model):
return (0, sum(model.labor[prod] * model.production[prod] for prod in model.products), model.total_labor)
def total_material_rule(model):
return (0, sum(model.material[prod] * model.production[prod] for prod in model.products), model.total_material)
model.labor_consumed = pyo.Constraint(rule = total_labor_rule)
model.material_consumed = pyo.Constraint(rule = total_material_rule)
Ya tenemos el problema definido. Solo queda crear una instancia con unos datos dados, de tal forma que el algoritmo pueda encontrar la solución óptima. Existen diferentes formas de definir los datos. En esta página puedes encontrar todas las opciones posibles. Yo personalmente te explicaré la más común: el uso de diccionarios.
Como digo, la forma más común de pasar datos a Pyomo es mediante diccionarios de Python. En este sentido, debemos crear un diccionario cuya clave es None y cuyo valor es, a su vez, un diccionario donde definiremos los datos. En este segundo diccionario, las claves deben ser los nombres de los atributos que hemos definido en el modelo, en mi caso labor, material, etc.
Si un objeto puede recibir varios valores (por ejemplo, el parámetro labor, el cual debe recibir un valor para cada objeto), deberemos crear un diccionario con el nombre del producto como clave y el valor de dicho parámetro u objeto como valor. Ejemplo:
...
"labor": {
"Product A": 2,
"Product B": 4
},
...
Por el contrario, si un parámetro es global y no depende del set (como el número total de horas de trabajo), la clave debe ser None y el valor el valor que queramos pasar al parámetro u objeto. Ejemplo:
...
"total_labor": {None: 500}
...
Así pues, sabiendo eso, los datos que pasaremos al modelo son los siguientes:
data = {
None : {
"products": ["Product A", "Product B"],
"labor": {
"Product A": 2,
"Product B": 4
},
"material": {
"Product A": 3,
"Product B": 2
},
"product_profit": {
"Product A": 10,
"Product B": 15
},
"total_labor": {None: 500},
"total_material": {None: 800}
}
}
Ahora que tenemos los datos, podemos crear una instancia del optimizador con estos datos y solucionar el problema, tal como muestro a continuación:
instance = model.create_instance(data)
optimizer = pyo.SolverFactory("glpk")
optimizer.solve(instance)
instance.display()
Model unknown
Variables:
production : Size=2, Index=products
Key : Lower : Value : Upper : Fixed : Stale : Domain
Product A : 0 : 250.0 : None : False : False : NonNegativeReals
Product B : 0 : 0.0 : None : False : False : NonNegativeReals
Objectives:
total_profit : Size=1, Index=None, Active=True
Key : Active : Value
None : True : 2500.0
Constraints:
labor_consumed : Size=1
Key : Lower : Body : Upper
None : 0.0 : 500.0 : 500
material_consumed : Size=1
Key : Lower : Body : Upper
None : 0.0 : 750.0 : 800
Como puedes ver, en este caso el resultado del optimizador es el mismo (producir 250 unidades de producto A). La principal diferencia es que, si los valores del modelo cambiasen, ya sea que hubiera más productos o, simplemente, cambiasen los datos, únicamente deberíamos cargar dichos valores en un diccionario y deberíamos volver a ejecutar las 5 últimaas líneas que he mostrado previamente.
Hecho esto, este ha sido un ejemplo bastante sencillo. Ahora veamos un ejemplo algo más complejo para que veas Pyomo funcionando cuando los parámetros dependen de dos sets de datos.
Ejemplo 2: Optimización de Fabricación y transporte a clientes con Pyomo
Eres el gerente de logística responsable de optimizar la distribución de productos de múltiples fábricas a varias ubicaciones de clientes. Tiene un conjunto de fábricas F = {Factory1, Factory2, Factory3} y un conjunto de ubicaciones de clientes C = {Client1, Client2, Client3, Client4}.
Los costes de transporte (en dólares) por unidad de producto desde cada fábrica hasta cada cliente se dan en la siguiente tabla:
| Client1 | Client2 | Client3 | Client4 |
|---|---|---|---|
| Factory1 | 10 | 15 | 20 |
| Factory2 | 18 | 12 | 10 |
| Factory3 | 16 | 19 | 13 |
Las fábricas tienen capacidades de producción de 100, 120 y 80 unidades por día, respectivamente. Además, la demanda de los clientes para cada ubicación es de 80, 70, 60 y 90 unidades por día, respectivamente.
El objetivo es determinar la asignación óptima de cantidades de productos de cada fábrica a cada ubicación de cliente, minimizando el costo total de transporte y satisfaciendo las capacidades de producción y las demandas de los clientes.
Solución con modelo concreto
Si queremos crear un modelo concreto, veremos que la variable es el número de unidades de producto que se transportarán de cada fábrica a cada cliente, es decir, que tendremos 12 variables a definir. Asimismo, contaremos con dos restricciones: la producción por fábrica (el cual define el valor máximo) y la demanda por cliente (la cual se debe cumplir, es decir, que sea igual a ese valor).
Por último, el objetivo es minimizar el coste, por lo que tendremos que calcular el coste total, como multiplicación de las unidades enviadas de cada fábrica a cada cliente.
Así pues, la solución será la siguiente:
model = pyo.ConcreteModel()
# Define the production in each factory
model.factory1_client1 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory1_client2 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory1_client3 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory1_client4 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory2_client1 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory2_client2 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory2_client3 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory2_client4 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory3_client1 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory3_client2 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory3_client3 = pyo.Var(within = pyo.NonNegativeIntegers)
model.factory3_client4 = pyo.Var(within = pyo.NonNegativeIntegers)
# Define the
# Set the maximum production of each factory
model.factory1_production = pyo.Constraint(
expr = (model.factory1_client1 + model.factory1_client2 +
model.factory1_client3 + model.factory1_client4) <= 100
)
model.factory2_production = pyo.Constraint(
expr = (model.factory2_client1 + model.factory2_client2 +
model.factory2_client3 + model.factory2_client4) <= 120
)
model.factory3_production = pyo.Constraint(
expr = (model.factory3_client1 + model.factory3_client2 +
model.factory3_client3 + model.factory3_client4) <= 80
)
# Meet customer demand
model.client1_demand = pyo.Constraint(
expr = (model.factory1_client1 + model.factory2_client1 +
model.factory3_client1) == 80
)
model.client2_demand = pyo.Constraint(
expr = (model.factory1_client2 + model.factory2_client2 +
model.factory3_client2) == 70
)
model.client3_demand = pyo.Constraint(
expr = (model.factory1_client3 + model.factory2_client3 +
model.factory3_client3) == 60
)
model.client4_demand = pyo.Constraint(
expr = (model.factory1_client4 + model.factory2_client4 +
model.factory3_client4) == 90
)
def total_cost(model):
cost = (model.factory1_client1 * 10) + (model.factory1_client2 * 15) \
+ (model.factory1_client3 * 20) + (model.factory1_client4 * 12) \
+ (model.factory2_client1 * 18) + (model.factory2_client2 * 12) \
+ (model.factory2_client3 * 10) + (model.factory2_client4 * 14) \
+ (model.factory3_client1 * 16) + (model.factory3_client2 * 19) \
+ (model.factory3_client3 * 13) + (model.factory3_client4 * 10)
return cost
model.total_transportation_cost = pyo.Objective(
rule = total_cost,
sense = pyo.minimize
)
pyo.SolverFactory('glpk').solve(model).write()
# Print the results
print(f"Total Return = {model.total_transportation_cost()}")
print(f"Factory 1, Client 1 = {model.factory1_client1()}")
print(f"Factory 1, Client 2 = {model.factory1_client2()}")
print(f"Factory 1, Client 3 = {model.factory1_client3()}")
print(f"Factory 1, Client 4 = {model.factory1_client4()}")
print(f"Factory 2, Client 1 = {model.factory2_client1()}")
print(f"Factory 2, Client 2 = {model.factory2_client2()}")
print(f"Factory 2, Client 3 = {model.factory2_client3()}")
print(f"Factory 2, Client 4 = {model.factory2_client4()}")
print(f"Factory 3, Client 1 = {model.factory3_client1()}")
print(f"Factory 3, Client 2 = {model.factory3_client2()}")
print(f"Factory 3, Client 3 = {model.factory3_client3()}")
print(f"Factory 3, Client 4 = {model.factory3_client4()}")
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 3190.0
Upper bound: 3190.0
Number of objectives: 1
Number of constraints: 7
Number of variables: 12
Number of nonzeros: 24
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 1
Number of created subproblems: 1
Error rc: 0
Time: 0.053026676177978516
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
Total Return = 3190.0
Factory 1, Client 1 = 80.0
Factory 1, Client 2 = 10.0
Factory 1, Client 3 = 0.0
Factory 1, Client 4 = 10.0
Factory 2, Client 1 = 0.0
Factory 2, Client 2 = 60.0
Factory 2, Client 3 = 60.0
Factory 2, Client 4 = 0.0
Factory 3, Client 1 = 0.0
Factory 3, Client 2 = 0.0
Factory 3, Client 3 = 0.0
Factory 3, Client 4 = 80.0
Como ves, solucionar problemas algo más complejos mediante modelos concretos, es sencillo, pero bastante laborioso. En mi opinión, aunque sea algo más complejo, suele ser recomendable crear modelos abstractos. Veamos ahora cómo definir el modelo abstracto para este caso de optimización con Pyomo.
Solución mediante modelo abstracto
En este caso, crear el modelo abstracto es algo más complejo, así que lo iremos creando poco a poco.
Para empezar, hay que tener en cuenta que trabajar con dos sets que desconocemos a priori: un set de fábricas y otro set de clientes. Así pues, comenzamos creando dichos sets:
import pyomo.environ as pyo
model = pyo.AbstractModel()
# Define the sets of values
model.factories = pyo.Set()
model.clients = pyo.Set()
Respecto a los parámetros (datos) que recibe el modelo, debemos indicar:
- Una cantidad máxima de fabricación por fábrica.
- Una demanda máxima por cliente.
- Un coste de transporte de cada fábrica a cada cliente.
Si bien los dos primeros puntos es algo que ya hemos visto previamente, en el segundo caso deberemos indicar que debe haber un valor para cada combinación de fábrica y cliente. Esto lo hacemos pasando ambos parámetros, tal como se muestra a continuación:
# Define the input data
model.transportation_cost = pyo.Param(model.factories, model.clients)
model.factory_production = pyo.Param(model.factories)
model.customer_demand = pyo.Param(model.clients)
Asimismo, con la variable ocurre lo mismo, debemos definir cuánto queremos producir y enviar de cada fábrica a cada cliente. Así pues, debe haber una cantidad de fabricación por cada combinación posible de fábrica y cliente. En este caso, además, hay que tener en cuenta que la fabricación no puede ser negativa y que el número de unidades enviadas debe ser un número entero:
# Define the model variable
model.supply_factory_client = pyo.Var(model.factories, model.clients, within = pyo.NonNegativeIntegers)
Por otro lado, debemos definir las restricciones del modelo. Aquí introducimos un concepto algo más complejo. Y es que, en el caso de las unidades fabricadas, deberemos introducir una restricción por fábrica (con la demanda de clientes ocurre lo mismo).
Para ello, podemos crear una función, que reciba un parámetro factory, el cual indicará sobre qué fábrica vamos a fijar el constraint. De esta forma, si pasamos el set factories al constraint, automáticamente iterará para cada fábrica creando, para cada una de ellas, su constraint correspondiente:
def maximum_production_rule(model, factory):
return (0, sum(model.supply_factory_client[factory, client] for client in model.clients), model.factory_production[factory])
model.maximum_production_rule = pyo.Constraint(model.factories, rule=maximum_production_rule)
Ahora podemos aplicar esta misma idea pero para calcular la demanda del cliente:
def meet_demand_rule(model, client):
return sum(model.supply_factory_client[factory, client] for factory in model.factories) == model.customer_demand[client]
model.meet_demand_rule = pyo.Constraint(model.clients, rule=meet_demand_rule)
Por último, deberemos definir el objetivo como la suma de costes de transporte de cada fábrica a cada cliente. Para ello, tenemos que tener en cuenta que el coste de una fábrica a un cliente se accede de la siguiente manera: model.transportation_cost[factory, client].
Sabiendo esto, definimos el objetivo:
def total_cost(model):
total_cost = 0
for client in model.clients:
for factory in model.factories:
total_cost =+ model.transportation_cost[factory, client]
return total_cost
model.total_transportation_cost = pyo.Objective(
rule = total_cost,
sense = pyo.minimize
)
Con esto ya tendríamos el modelo abstracto creado. Solo faltaría crear una instancia a partir de los datos y solucionarla. En este sentido, comentar que si tenemos un parámetro que dependa de varios sets (como el coste de transporte), a la hora de definirlo deberemos pasar una tupla con los distintos valores de cada set, cada uno en la posición que se haya definido.
Veamos pues, cómo definir nuestros datos:
data = {
None: {
"factories": ["Factory1", "Factory2", "Factory3"],
"clients": ["Customer1", "Customer2", "Customer3", "Customer4"],
"transportation_cost": {
("Factory1", "Customer1"): 10,
("Factory1", "Customer2"): 15,
("Factory1", "Customer3"): 20,
("Factory1", "Customer4"): 12,
("Factory2", "Customer1"): 18,
("Factory2", "Customer2"): 12,
("Factory2", "Customer3"): 10,
("Factory2", "Customer4"): 14,
("Factory3", "Customer1"): 16,
("Factory3", "Customer2"): 19,
("Factory3", "Customer3"): 13,
("Factory3", "Customer4"): 10
},
"customer_demand": {
"Customer1": 80,
"Customer2": 70,
"Customer3": 60,
"Customer4": 90
},
"factory_production": {
"Factory1" : 100,
"Factory2": 120,
"Factory3": 80
}
}
}
Finalmente, ahora que tneemos los datos definidos podemos solucionar el problema:
instance = model.create_instance(data)
optimizer = pyo.SolverFactory("glpk")
optimizer.solve(instance)
instance.display()
Model unknown
Variables:
supply_factory_client : Size=12, Index=supply_factory_client_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('Factory1', 'Customer1') : 0 : 80.0 : None : False : False : NonNegativeIntegers
('Factory1', 'Customer2') : 0 : 10.0 : None : False : False : NonNegativeIntegers
('Factory1', 'Customer3') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory1', 'Customer4') : 0 : 10.0 : None : False : False : NonNegativeIntegers
('Factory2', 'Customer1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory2', 'Customer2') : 0 : 60.0 : None : False : False : NonNegativeIntegers
('Factory2', 'Customer3') : 0 : 60.0 : None : False : False : NonNegativeIntegers
('Factory2', 'Customer4') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory3', 'Customer1') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory3', 'Customer2') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory3', 'Customer3') : 0 : 0.0 : None : False : False : NonNegativeIntegers
('Factory3', 'Customer4') : 0 : 80.0 : None : False : False : NonNegativeIntegers
Objectives:
total_transportation_cost : Size=1, Index=None, Active=True
Key : Active : Value
None : True : 10
Constraints:
maximum_production_rule : Size=3
Key : Lower : Body : Upper
Factory1 : 0.0 : 100.0 : 100.0
Factory2 : 0.0 : 120.0 : 120.0
Factory3 : 0.0 : 80.0 : 80.0
meet_demand_rule : Size=4
Key : Lower : Body : Upper
Customer1 : 80.0 : 80.0 : 80.0
Customer2 : 70.0 : 70.0 : 70.0
Customer3 : 60.0 : 60.0 : 60.0
Customer4 : 90.0 : 90.0 : 90.0
Conclusión y próximos pasos
En mi opinión, el mundo de la optimización es un mundo muy interesante que aporta mucho valor a las empresas. En este post sobre cómo crear modelos de optimización en Python con Pyomo, me he centrado en Pyomo, pero no es el único framework. Si consideras que la optimización de procesos es algo que deberías añadir a tu lista de herramientas, te recomiendaría no probar diferentes frameworks y quedarte con aquel que más te guste. De hecho, te recomendaría también probar otro tipo de enfoques, como el de los algoritmos genéticos (de los cuáles hablé en este post).
Por otro lado, si Pyomo te ha gustado y te ha resultado interesante, te animo a seguir practicando con diferentes ejercicios. En este sentido, ChatGPT puede ser de gran ayuda para generar nuevos problemas que solucionar. De hecho, las ideas de los ejemplos que he utilizado en este post me las dió ChatGPT.
En cualquier caso, espero que este post te haya gustado. ¡Nos vemos en el siguiente!