Esta notebook fue creada originalmente como un blog post por Raúl E. López Briega en Matemáticas, análisis de datos y python. El contenido esta bajo la licencia BSD.
La optimización es fundamental para cualquier problema relacionado con la toma de decisiones, ya sea en ingeniería o en ciencias económicas. La tarea de tomar decisiones implica elegir entre varias alternativas. Esta opción va a estar gobernada por nuestro deseo de tomar la "mejor" decisión posible. Que tan buena va a ser cada una de las alternativas va a estar descripta por una función objetivo o índice de rendimiento. La teoría y los métodos de optimización nos van a ayudar a seleccionar la mejor alternativa de acuerdo a esta función objetivo dada. El área de optimización ha recibido gran atención en los últimos años, principalmente por el rápido desarrollo de las ciencias de computación, incluido el desarrollo y la disponibilidad de herramientas de software sumamente amigables, procesadores paralelos de alta velocidad, y redes neuronales artificiales. El poder de los métodos de optimización reside en la posibilidad de determinar la solución óptima sin realmente tener que probar todos los casos posibles. Para logar esto, se utiliza un nivel modesto de Matemáticas y se realizan cálculos numéricos iterativos utilizando procedimientos lógicos claramente definidos o algoritmos implementados en computadoras.
Un problema de optimización comienza con un conjunto de variables independientes o parámetros, y a menudo incluye condiciones o restricciones que definen los valores aceptables de estas variables. Tales restricciones se denominan las limitaciones del problema. El otro componente esencial de un problema de optimización es una medida única de "bondad", denominada función objetivo, la cual va a depender también de las variables independientes. La solución al problema de optimización va a estar dada por un conjunto de valores permitidos para las variables independientes, de acuerdo con las limitaciones del problema, en los que la función objetivo asume un valor óptimo. En términos matemáticos, la optimización implica normalmente maximizar o minimizar la función objetivo.
Para aplicar los resultados matemáticos y técnicas numéricas de la teoría de optimización a problemas concretos, es necesario delinear claramente los límites del sistema a optimizar, definir los parámetros cuantitativos que se utilizaran como criterio en base al cual serán clasificadas las alternativas para determinar la "mejor" opción, seleccionar las variables que se utilizarán para caracterizar o identificar alternativas; y definir el modelo que exprese la forma en que las variables estarán relacionadas. Una buena formulación del problema es la clave para el éxito de un problema de optimización y es en alto grado un arte. Se aprende a través de la práctica y el estudio de aplicaciones exitosas y se basa en el conocimiento de las fortalezas, debilidades y particularidades de las técnicas proporcionadas por los distintos métodos de optimización.
Un paso importante en cualquier problema de optimización es clasificar nuestro modelo, ya que los algoritmos para resolver problemas de optimización generalmente están diseñados para atacar un tipo de problema en particular. Algunas de las principales clasificaciones que vamos a poder encontrar son las siguientes:
Matemáticamente, un problema de optimización con restricciones asume la siguiente forma:
$$ \begin{array}{ll} \min \hspace{1cm} f_0(x)\\ \mbox{sujeto a } \ f_i (x) \leq b_i, \hspace{1cm} i=1, \dots, m. \end{array} $$En donde el vector $x = (x_1, \dots, x_n)$ es la variable de optimización del problema, la función $f_0: \mathbb{R}^n \rightarrow \mathbb{R}$ es la función objetivo, las funciones $f_i: \mathbb{R}^n \rightarrow \mathbb{R}, i=1, \dots, m$; son las funciones de restricciones de desigualdad; y las constantes $b_1, \dots, m$ son los límites de las restricciones. Dentro de este marco, un caso importante es el de la programación lineal, en el cual la función objetivo y las restricciones son lineales. El objetivo de la programación lineal es determinar los valores de las variables de decisión que maximizan o minimizan una función objetivo lineal, y en donde las variables de decisión están sujetas a restricciones lineales. En general, el objetivo es encontrar un punto que minimice la función objetivo al mismo tiempo que satisface las restricciones.
Para resolver un problema de programación lineal, debemos seguir los siguientes pasos:
Elegir las incógnitas o variables de decisión.
Escribir la función objetivo en función de los datos del problema.
Escribir las restricciones en forma de sistema de inecuaciones.
Averiguar el conjunto de soluciones factibles representando gráficamente las restricciones.
Calcular las coordenadas de los vértices del recinto de soluciones factibles (si son pocos).
Calcular el valor de la función objetivo en cada uno de los vértices para ver en cuál de ellos presenta el valor máximo o mínimo según nos pida el problema (hay que tener en cuenta aquí la posible no existencia de solución).
Uno de los algoritmos más eficientes para resolver problemas de programación lineal, es el método simplex.
Otro caso importante de optimización que debemos destacar es el de la optimización convexa, en el cual la función objetivo y las restricciones son convexas. En realidad, la programación lineal que vimos anteriormente, no es más que un caso especial de optimización convexa.
Podemos decir que un un conjunto es convexo si se puede ir de cualquier punto a cualquier otro en línea recta, sin salir del mismo. El concepto de convexidad es el opuesto de concavidad. LLevado a las funciones, podemos decir que una función es convexa en un intervalo (a,c), si para todo punto b del intervalo la recta tangente en ese punto queda por debajo de la función.
Los conjuntos y funciones convexas tienen algunas propiedades que los hacen especiales para problemas de optimización, como ser:
Una función convexa no tiene mínimos locales que no sean globales.
Un conjunto convexo tiene un interior relativo no vacío.
Un conjunto convexo está conectado y tiene direcciones factibles en cualquier punto.
Una función no convexa puede ser convexificada manteniendo al mismo tiempo lo óptima de sus mínimos globales.
Una función convexa es continua dentro del interior de su dominio, y tiene buenas propiedades de diferenciación.
entre otras
Como es de esperar, en el ecosistema científico de Python podemos encontrar varias librerías que nos van a ayudar a enfrentar los problemas de optimización, entre las que podemos destacar:
Debemos destacar que tanto PuLP como Pyomo requieren la instalación adicional de diferentes solvers para poder resolver los problemas de optimización. Algunos de los solvers que soportan son: GLPK, COIN CLP/CBC, CPLEX y GURUBI, entre otros.
Bien, ahora llegó el momento de ver como podemos resolver algunos problemas de optimización con la ayuda de las librerías antes mencionadas.
En general, podemos resolver problemas de Mínimos cuadrados utilizando un poco de álgebra lineal, pero los Mínimos cuadrados también pueden ser vistos como un problema de optimización; ya que como su nombre lo indica, debemos minimizar la suma de los cuadrados de las diferencias entre los puntos generados por la función elegida y los correspondientes valores en los datos. scipy.optimize nos ofrece algunos métodos para resolver este tipo de problemas utilizando técnicas de optimización, como por ejemplo el algoritmo Gauss-Newton. Estos métodos pueden sernos de mucha utilizar sobre todo si la función tiene componentes no lineales. Veamos un ejemplo
In [1]:
# <!-- collapse=True -->
# importando modulos necesarios
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
import cvxopt
import pulp
from pyomo.environ import *
from pyomo.opt import SolverFactory
import pyomo.environ
np.random.seed(1984) #replicar random
%matplotlib inline
In [2]:
# Ejemplo mínimos cuadrados no lineales utilizando scipy.optimize
beta = (0.25, 0.75, 0.5)
# funcion modelo
def f(x, b0, b1, b2):
return b0 + b1 * np.exp(-b2 * x**2)
# datos aleatorios para simular las observaciones
xdata = np.linspace(0, 5, 50)
y = f(xdata, *beta)
ydata = y + 0.05 * np.random.randn(len(xdata))
# función residual
def g(beta):
return ydata - f(xdata, *beta)
# comenzamos la optimización
beta_start = (1, 1, 1)
beta_opt, beta_cov = optimize.leastsq(g, beta_start)
beta_opt
Out[2]:
In [3]:
# graficamos
fig, ax = plt.subplots(figsize=(10,8))
ax.scatter(xdata, ydata)
ax.plot(xdata, y, 'r', lw=2)
ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2)
ax.set_xlim(0, 5)
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)
ax.set_title('Mínimos cuadrados no lineales')
plt.show()
Las restricciones añaden otro nivel de complejidad a los problemas de optimización. En su forma más simple, simplemente consiste en poner algunos límites sobre los valores que las variables pueden alcanzar. Por ejemplo, podrías intentar resolver el siguiente problema:
$$\min \ f(x)= (x_1 -1)^2 - (x_2 -1)^2 \hspace{1cm} \mbox{sujeto a } \ 2 \leq x_1 \leq 3 \mbox{ y } 0 \leq x_2 \leq 2 $$Este tipo de problema se puede resolver utilizando el método L-BFGS-B que nos ofrece scipy.optimize.
In [4]:
# Ejemplo de optimización con restricciones scipy
# defino una funcion de ayuda para subregion en el gráfico
def func_X_Y_to_XY(f, X, Y):
"""
Wrapper for f(X, Y) -> f([X, Y])
"""
s = np.shape(X)
return f(np.vstack([X.ravel(), Y.ravel()])).reshape(*s)
# función a minimizar
def f(X):
x, y = X
return (x - 1)**2 + (y - 1)**2
# minimizo la función si restricciones
x_opt = optimize.minimize(f, (1, 1), method='BFGS').x
# el mínimo para las restricciones
bnd_x1, bnd_x2 = (2, 3), (0, 2)
x_cons_opt = optimize.minimize(f, np.array([1, 1]), method='L-BFGS-B',
bounds=[bnd_x1, bnd_x2]).x
# graficando la solución
fig, ax = plt.subplots(figsize=(10, 8))
x_ = y_ = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 50)
ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15)
ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15)
bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]),
bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0],
facecolor="grey")
ax.add_patch(bound_rect)
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
plt.colorbar(c, ax=ax)
ax.set_title('Optimización con restricciones')
plt.show()
Las restricciones que se definen por igualdades o desigualdades que incluyen más de una variable son más complicadas de tratar. Sin embargo, existen técnicas generales que también podemos utilizar para este tipo de problemas. Volviendo al ejemplo anterior, cambiemos la restricción por una más compleja, como ser: $g(x) = x_1 - 1.75 - ( x_0 - 0.75 )^4 \geq 0$. Para resolver este problema scipy.optimize nos ofrece un método llamado programación secuencial por mínimos cuadrados, o SLSQP por su abreviatura en inglés.
In [5]:
# Ejemplo scipy SLSQP
# funcion de restriccion
def g(X):
return X[1] - 1.75 - (X[0] - 0.75)**4
# definimos el diccionario con la restricción
restriccion = dict(type='ineq', fun=g)
# resolvemos
x_opt = optimize.minimize(f, (0, 0), method='BFGS').x
x_cons_opt = optimize.minimize(f, (0, 0), method='SLSQP',
constraints=restriccion).x
# graficamos
ig, ax = plt.subplots(figsize=(10, 8))
x_ = y_ = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 50)
ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15)
ax.plot(x_, 1.75 + (x_-0.75)**4, 'k-', markersize=15)
ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color='grey')
ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15)
ax.set_ylim(-1, 3)
ax.set_xlabel(r"$x_0$", fontsize=18)
ax.set_ylabel(r"$x_1$", fontsize=18)
plt.colorbar(c, ax=ax)
plt.show()
Ahora veamos un ejemplo de programación lineal utilizando los métodos de optimización convexa que nos ofrece la librería CVXopt. Recordemos que la programación lineal es un caso especial de la optimización convexa y el principal algoritmo que se utiliza es el método simplex. Para este pequeño ejemplo vamos a maximizar la siguiente función objetivo: $$f(x_1,x_2) = 50x_1 + 40x_2$$
con las siguientes restricciones:
$$x_{1} + 1.5x_{2} \leq 750 \\ 2x_1 + x_2 \leq 1000 \\ x_1 \geq 0 \\ x_2 \geq 0 $$
In [6]:
# Ejemplo programación lineal con CVXopt
# Resolviendo el problema con cvxopt
A = cvxopt.matrix([[-1., -2., 1., 0.], # columna de x1
[-1.5, -1., 0., 1.]]) # columna de x2
b = cvxopt.matrix([750., 1000., 0., 0.]) # resultados
c = cvxopt.matrix([50., 40.]) # funcion objetivo
# resolviendo el problema
sol=cvxopt.solvers.lp(c,A,b)
In [7]:
# imprimiendo la solucion.
print('{0:.2f}, {1:.2f}'.format(sol['x'][0]*-1, sol['x'][1]*-1))
In [8]:
# Resolviendo la optimizacion graficamente.
x_vals = np.linspace(0, 800, 10) # 10 valores entre 0 y 800
y1 = ((750 - x_vals)/1.5) # x1 + 1.5x2 = 750
y2 = (1000 - 2*x_vals) # 2x1 + x2 = 1000
plt.figure(figsize=(10,8))
plt.plot(x_vals, y1, label=r'$x_1 + 1.5x_2 \leq 750$')
plt.plot(x_vals, y2, label=r'$2x_1 + x_2 \leq 1000$') #
plt.plot(375, 250, 'b*', markersize=15)
# Región factible
y3 = np.minimum(y1, y2)
plt.fill_between(x_vals, 0, y3, alpha=0.15, color='b')
plt.axis(ymin = 0)
plt.title('Optimización lineal')
plt.legend()
plt.show()
Como podemos ver, tanto la solución utilizando CVXopt, como la solución gráfica; nos devuelven el mismo resultado $x_1=375$ y $x_2=250$.
El problema de transporte es un problema clásico de programación lineal en el cual se debe minimizar el costo del abastecimiento a una serie de puntos de demanda a partir de un grupo de puntos de oferta, teniendo en cuenta los distintos precios de envío de cada punto de oferta a cada punto de demanda. Por ejemplo, supongamos que tenemos que enviar cajas de cervezas de 2 cervecerías a 5 bares de acuerdo al siguiente gráfico:
Asimismo, supongamos que nuestro gerente financiero nos informa que el costo de transporte por caja de cada ruta se conforma de acuerdo a la siguiente tabla:
| Bar 1 | Bar 2 | Bar 3 | Bar 4 | Bar 5 | |
|---|---|---|---|---|---|
| Cervecería A | 2 | 4 | 5 | 2 | 1 |
| Cervecería B | 3 | 1 | 3 | 2 | 3 |
Y por último, las restricciones del problema, van a estar dadas por las capacidades de oferta y demanda de cada cervecería y cada bar, las cuales se detallan en el gráfico de más arriba.
Veamos como podemos modelar este ejemplo con la ayuda de PuLP y de Pyomo.
In [9]:
# Ejemplo del problema de transporte de las cervezas utilizando PuLP
# Creamos la variable prob que contiene los datos del problema
prob = pulp.LpProblem("Problema de distribución de cerveza", pulp.LpMinimize)
In [10]:
# Creamos lista de cervecerías o nodos de oferta
cervecerias = ["Cervecería A", "Cervercería B"]
# diccionario con la capacidad de oferta de cada cerveceria
oferta = {"Cervecería A": 1000,
"Cervercería B": 4000}
# Creamos la lista de los bares o nodos de demanda
bares = ["Bar 1", "Bar 2", "Bar 3", "Bar 4", "Bar 5"]
# diccionario con la capacidad de demanda de cada bar
demanda = {"Bar 1":500,
"Bar 2":900,
"Bar 3":1800,
"Bar 4":200,
"Bar 5":700,}
# Lista con los costos de transporte de cada nodo
costos = [ #Bares
#1 2 3 4 5
[2,4,5,2,1],#A Cervecerías
[3,1,3,2,3] #B
]
# Convertimos los costos en un diccionario de PuLP
costos = pulp.makeDict([cervecerias, bares], costos,0)
# Creamos una lista de tuplas que contiene todas las posibles rutas de tranporte.
rutas = [(c,b) for c in cervecerias for b in bares]
In [11]:
# creamos diccionario x que contendrá la candidad enviada en las rutas
x = pulp.LpVariable.dicts("ruta", (cervecerias, bares),
lowBound = 0,
cat = pulp.LpInteger)
# Agregamos la función objetivo al problema
prob += sum([x[c][b]*costos[c][b] for (c,b) in rutas]), \
"Suma_de_costos_de_transporte"
# Agregamos la restricción de máxima oferta de cada cervecería al problema.
for c in cervecerias:
prob += sum([x[c][b] for b in bares]) <= oferta[c], \
"Suma_de_Productos_que_salen_de_cervecerias_%s"%c
# Agregamos la restricción de demanda mínima de cada bar
for b in bares:
prob += sum([x[c][b] for c in cervecerias]) >= demanda[b], \
"Sum_of_Products_into_Bar%s"%b
# Los datos del problema son exportado a archivo .lp
prob.writeLP("problemaDeTransporte.lp")
# Resolviendo el problema.
prob.solve()
# Imprimimos el estado del problema.
print("Status: {}".format(pulp.LpStatus[prob.status]))
# Imprimimos cada variable con su solución óptima.
for v in prob.variables():
print("{0:} = {1:}".format(v.name, v.varValue))
# Imprimimos el valor óptimo de la función objetivo
print("Costo total de transporte = {}".format(prob.objective.value()))
Como vemos, la solución óptima que encontramos con la ayuda de PuLP, nos dice que deberíamos enviar desde la Cervecería A, 300 cajas al Bar 1 y 700 cajas al Bar 5; y que desde la Cervecería B deberíamos enviar 200 cajas al Bar 1, 900 cajas al Bar 2, 1800 cajas al Bar 3 y 200 cajas al Bar 4. De esta forma podemos minimizar el costo de transporte a un total de 8600.
Veamos ahora si podemos formular este mismo problema con Pyomo, así podemos darnos una idea de las diferencias entre las herramientas.
In [12]:
# Ejemplo del problema de transporte de las cervezas utilizando Pyomo
# Creamos el modelo
modelo = ConcreteModel()
# Creamos los nodos de oferta y demanda
modelo.i = Set(initialize=['Cervecería A','Cervecería B'], doc='Cervecerías')
modelo.j = Set(initialize=['Bar 1', 'Bar 2', 'Bar 3', 'Bar 4', 'Bar 5'],
doc='Bares')
# Definimos las capacidades de oferta y demanda
modelo.a = Param(modelo.i, initialize={'Cervecería A':1000,'Cervecería B':4000},
doc='Capacidad de oferta de las cervecerías')
modelo.b = Param(modelo.j, initialize={'Bar 1':500,'Bar 2':900,'Bar 3':1800,
'Bar 4':200, 'Bar 5':700 },
doc='Demanda de cada bar')
# Costo de transporte
costos = {
('Cervecería A', 'Bar 1'): 2,
('Cervecería A', 'Bar 2'): 4,
('Cervecería A', 'Bar 3'): 5,
('Cervecería A', 'Bar 4'): 2,
('Cervecería A', 'Bar 5'): 1,
('Cervecería B', 'Bar 1'): 3,
('Cervecería B', 'Bar 2'): 1,
('Cervecería B', 'Bar 3'): 3,
('Cervecería B', 'Bar 4'): 2,
('Cervecería B', 'Bar 5'): 3
}
modelo.d = Param(modelo.i, modelo.j, initialize=costos,
doc='Costo de transporte')
# definimos el costo de tranporte
def f_costo(modelo, i, j):
return modelo.d[i,j]
modelo.c = Param(modelo.i, modelo.j, initialize=f_costo,
doc='Costo de transporte')
# definimos variable x con las cantidades de cajas enviadas
modelo.x = Var(modelo.i, modelo.j, bounds=(0.0,None),
doc='Cantidad de cajas')
In [13]:
## Definimos las restricciones ##
# Límite de oferta
def f_oferta(modelo, i):
return sum(modelo.x[i,j] for j in modelo.j) <= modelo.a[i]
modelo.oferta = Constraint(modelo.i, rule=f_oferta,
doc='Límites oferta de cada Cervecería')
# Límite de demanda
def f_demanda(modelo, j):
return sum(modelo.x[i,j] for i in modelo.i) >= modelo.b[j]
modelo.demanda = Constraint(modelo.j, rule=f_demanda,
doc='Límites demanda de cada bar')
In [14]:
## Definimos la función objetivo y resolvemos el problema ##
# Función objetivo
def f_objetivo(modelo):
return sum(modelo.c[i,j]*modelo.x[i,j] for i in modelo.i for j in modelo.j)
modelo.objetivo = Objective(rule=f_objetivo, sense=minimize,
doc='Función Objetivo')
# resolvemos el problema e imprimimos resultados
def pyomo_postprocess(options=None, instance=None, results=None):
modelo.x.display()
# utilizamos solver glpk
opt = SolverFactory("glpk")
resultados = opt.solve(modelo)
# imprimimos resultados
print("\nSolución óptima encontrada\n" + '-'*80)
pyomo_postprocess(None, None, resultados)
Aquí concluye este artículo, como vemos la optimización se puede aplicar a un amplio rango de problemas, por lo que es sumamente importe conocer sus principales métodos y como podemos aplicarlos con éxito.
Saludos!
Este post fue escrito utilizando IPython notebook. Pueden descargar este notebook o ver su version estática en nbviewer.