Optimización Entera Mixta No Lineal (MINLP) con R y Pyomo

V Jornadas de Usuarios de R (Diciembre 2013) | R Madrid (Marzo 2014)

Jorge Ayuso Rejas

Optimización de un presupuesto

Problema de partida

Una empresa quiere repartir el presupuesto de una semana destinado a publicidad en los siguientes medios: Exterior, Online, Prensa, Radio, Revistas y Televisión.

Solución propuesta

Plantear un problema de optimización gracias a las siguientes curvas de aportes*:

Para cada medio \(i\) definimos:

\[f_i(x_i) = e^{A_i-B_i/x_i}\]

donde \(A_i,B_i>0\) parámetros a definir y \(x_i\) es la inversión semanal en euros en el medio \(i\).




*Llamamos aporte a cualquier variable de negocio interesante para el cliente: Ventas en volumnes, leads, llamadas...

Curvas de aporte

plot of chunk grafico_aporte

Planteamiento del problema


Problema 1 (NLP)

Maximizar: \[ \sum_{i \in \text{medios}} e^{A_i-B_i/x_i} \] sujeto a:

  • \(\sum_{i \in \text{medios}} x_i \leq\) Presupuesto Total
  • \(min_i\leq x_i \leq max_i\quad \forall i \in\) Medios


Problema 2 (MINLP)

Maximizar: \[ \sum_{i \in \text{medios}} w_i e^{A_i-B_i/x_i} \] sujeto a:

  • \(\sum_{i \in \text{medios}} w_i x_i \leq\) Presupuesto Total
  • \(w_i \in \{0,1\} \quad \forall i \in\) Medios
  • \(min_i\leq x_i \leq max_i\quad \forall i \in\) Medios

Implementación del problema: Búsqueda (I)

Paso 1 Buscamos algún paquete que nos sirva en: Task Views Optimization

Encontramos paquetes para resolver el problema 1 (NLP) pero ninguno para el problema 2 (MINLP).

Implementación del problema: Búsqueda (II)

Paso 2 Buscamos software/algoritmos libres que resuelvan el problema, encontramos:


  • Las tres soluciones son libres y compatibles con la libreria ASL:

The AMPL Solver Library (ASL) which allows to read the nl files and provides the automatic differentiation functionality is open-source. It is used in many solvers to implement AMPL connection.
Fuente: Wikipedia

  • Además permiten ser llamdas desde C así que se podría implementar una conexión en R.

Implementación del problema: Pyomo

(http://software.sandia.gov/trac/coopr/wiki/Pyomo)

The Python Optimization Modeling Objects (Pyomo) package is an open source tool for modeling optimization applications in Python. Pyomo can be used to define symbolic problems, create concrete problem instances, and solve these instances with standard solvers. Pyomo provides a capability that is commonly associated with algebraic modeling languages such as AMPL, AIMMS, and GAMS, but Pyomo's modeling objects are embedded within a full-featured high-level programming language with a rich set of supporting libraries.


  • También dispone de la librería ASL así que podemos conectar con los algoritmos anteriores.
  • Más fácil que implementar una conexión en C.
  • Diferenciación automática al usar ASL.

Pyomo: Ejemplo planteado

Veamos cómo definimos el problema en Pyomo:

# Imports
from coopr.pyomo import *
# ***********************************
model = AbstractModel()
# ***********************************
model.datos = Set()
# *********************************** Definimos los parámetros
model.MIN = Param(model.datos, within=NonNegativeReals)
model.MAX = Param(model.datos, within=NonNegativeReals)
model.A = Param(model.datos, within=NonNegativeReals)
model.B = Param(model.datos, within=NonNegativeReals)
model.Presupuesto = Param(within=NonNegativeReals)
# *********************************** Definimos las variables
def Level_bounds(model, i):
  return (model.MIN[i], model.MAX[i])
model.x = Var(model.datos, bounds=Level_bounds)
model.w = Var(model.datos, within=Binary)
# *********************************** Definimos la función objetivo, suma de las curvas de aporte
def Total_Cost_rule(model):
  return sum([model.w[j] *  exp(model.A[j]-model.B[j]/model.x[j]) for j in model.datos])
model.Total_Cost = Objective(rule=Total_Cost_rule,sense=maximize)
# *********************************** Introducimos la restricción de presupuesto total
def Demand_rule(model):
  return sum([model.w[i] * model.x[i] for i in model.datos]) <= model.Presupuesto
model.Demand = Constraint(rule=Demand_rule)

Pyomo: Convertir los datos al formato adecuado

Escribimos una función data_ampl que nos convierte el data.frame "datos" al formato adecuado para Pyomo:

datos
##                A      B    MIN    MAX
## Television 7.091 270164 183360 421967
## Online     6.572  62165  38488 110009
## Prensa     6.780 105753  71820 165048
## Radio      6.183  64580  43420 102131
## Exterior   4.169  26693  20854  34972
cat(paste(c(data_ampl(datos,name="datos"),data_ampl(570000,name="Presupuesto",tipo="param")),collapse="\n"))
## param: datos : "A" "B" "MIN" "MAX" :=
## "Television" 7.091114787 270163.8599 183359.5768 421967.2358
## "Online" 6.571585223 62164.76381 38488.03122 110009.0971
## "Prensa" 6.779516794 105752.5395 71820.20725 165047.6584
## "Radio" 6.18270522 64579.77595 43419.72799 102131.3559
## "Exterior" 4.168857373 26693.37017 20854.44233 34971.98656;
## param Presupuesto := 570000;

Los detalle de la función y objetos que podemos definir se pueden ver con más detalle en el documento metodológico.

Pyomo: Funcionamiento

Guardamos el modelo que hemos definido en modelo.py y los parámetros en datos.dat.
Podemos ejecutar el modelo con:

system("pyomo modelo.py datos.dat --solver bonmin", intern = TRUE)[1:13]
##  [1] "[    0.00] Setting up Pyomo environment"         
##  [2] "[    0.00] Applying Pyomo preprocessing actions" 
##  [3] "[    0.00] Creating model"                       
##  [4] "[    0.02] Applying solver"                      
##  [5] "[    0.08] Processing results"                   
##  [6] "    Number of solutions: 1"                      
##  [7] "    Solution Information"                        
##  [8] "      Gap: <undefined>"                          
##  [9] "      Status: optimal"                           
## [10] "      Function Value: 1429.08707683"             
## [11] "    Solver results file: results.json"           
## [12] "[    0.08] Applying Pyomo postprocessing actions"
## [13] "[    0.08] Pyomo Finished"

Los resultados se guardan en un archivo json que podemos leer fácilmente con el paquete RJSONIO.

Ejemplo Final

  1. Objetivo
    Queremos poder realizar optimizaciones de maner fácil y sencilla para cualquier usuario.

  2. Solución
    Usando Shiny logramos sin grandes esfuerzos una aplicacion web con el optimizador presentado.

  3. Explicación
    El usuario introduce los parámetros del problema, estos parámetros son procesados y enviados al optimizador. Después, recuperamos el resultado y lo representamos de manera gráfica gracias a los paquetes rCharts y googleVis.

    http://jornadas-r.conento.com

¡Gracias!