Programación lineal

La clase pasada vimos que los problemas de programación lineal eran problemas de optimización de funciones lineales sujetos a restricciones también lineales. Motivamos el estudio de esta disciplina con el concepto de complejidad computacional.

Adicionalmente vimos cómo llevar los problemas de programación lineal a una forma determinada que nos será bastante útil en esta clase.

Seguir instrucciones de acá: (ejecutar anaconda prompt como administrador)

https://stackoverflow.com/questions/19928878/installing-pyomo-on-windows-with-anaconda-python

1. Problemas de programación lineal

De acuerdo a lo visto la clase pasada, un problema de programación lineal puede escribirse en la siguiente forma:

\begin{equation} \begin{array}{ll} \min_{x_1,\dots,x_n} & f_1x_1+\dots+f_nx_n \\ \text{s. a. } & a^{eq}_{j,1}x_1+\dots+a^{eq}_{j,n}x_n=b^{eq}_j \text{ para } 1\leq j\leq m_1 \\ & a_{k,1}x_1+\dots+a_{k,n}x_n\leq b_k \text{ para } 1\leq k\leq m_2, \end{array} \end{equation}

donde:

$x_i$ para $i=1,\dots,n$ son las incógnitas o variables de decisión,
$f_i$ para $i=1,\dots,n$ son los coeficientes de la función a optimizar,
$a^{eq}_{j,i}$ para $j=1,\dots,m_1$ e $i=1,\dots,n$, son los coeficientes de la restricción de igualdad,
$a_{k,i}$ para $k=1,\dots,m_2$ e $i=1,\dots,n$, son los coeficientes de la restricción de desigualdad,
$b^{eq}_j$ para $j=1,\dots,m_1$ son valores conocidos que deben ser respetados estrictamente, y
$b_k$ para $k=1,\dots,m_2$ son valores conocidos que no deben ser superados.

Equivalentemente, el problema puede escribirse como

\begin{equation} \begin{array}{ll} \min_{\boldsymbol{x}} & \boldsymbol{f}^T\boldsymbol{x} \\ \text{s. a. } & \boldsymbol{A}_{eq}\boldsymbol{x}=\boldsymbol{b}_{eq} \\ & \boldsymbol{A}\boldsymbol{x}\leq\boldsymbol{b}, \end{array} \end{equation}

donde:

$\boldsymbol{x}=\left[x_1\quad\dots\quad x_n\right]^T$,
$\boldsymbol{f}=\left[f_1\quad\dots\quad f_n\right]^T$,
$\boldsymbol{A}_{eq}=\left[\begin{array}{ccc}a^{eq}_{1,1} & \dots & a^{eq}_{1,n}\\ \vdots & \ddots & \vdots\\ a^{eq}_{m_1,1} & \dots & a^{eq}_{m_1,n}\end{array}\right]$,
$\boldsymbol{A}=\left[\begin{array}{ccc}a_{1,1} & \dots & a_{1,n}\\ \vdots & \ddots & \vdots\\ a_{m_2,1} & \dots & a_{m_2,n}\end{array}\right]$,
$\boldsymbol{b}_{eq}=\left[b^{eq}_1\quad\dots\quad b^{eq}_{m_1}\right]^T$, y
$\boldsymbol{b}=\left[b_1\quad\dots\quad b_{m_2}\right]^T$.

Nota: el problema $\max_{\boldsymbol{x}}\boldsymbol{g}(\boldsymbol{x})$ es equivalente a $\min_{\boldsymbol{x}}-\boldsymbol{g}(\boldsymbol{x})$.

2. Explicar la función pyomo_utilities.py

Para trabajar con pyomo:

pyomo a parte de contener funciones para optimización, es en sí un lenguaje de programación. No tiene solucionadores instalados, utiliza solucionadores externos (glpk, ipopt).
Modelos concretos y modelos abstractos.
Forma de ingresarle los parámetros para un modelo abstracto.
Resultados.

Sin embargo, ya tenemos una función que hace todo el trabajo por nosotros. Únicamente debemos proporcionar los parámetros $\boldsymbol{f}$, $\boldsymbol{A}$ y $\boldsymbol{b}$ ($\boldsymbol{A}_{eq}$ y $\boldsymbol{b}_{eq}$, de ser necesario).

3. Ejemplos de la clase pasada

3.1

Una compañía produce dos productos ($X_1$ y $X_2$) usando dos máquinas ($A$ y $B$). Cada unidad de $X_1$ que se produce requiere 50 minutos en la máquina $A$ y 30 minutos en la máquina $B$. Cada unidad de $X_2$ que se produce requiere 24 minutos en la máquina $A$ y 33 minutos en la máquina $B$.

Al comienzo de la semana hay 30 unidades de $X_1$ y 90 unidades de $X_2$ en inventario. El tiempo de uso disponible de la máquina $A$ es de 40 horas y el de la máquina $B$ es de 35 horas.

La demanda para $X_1$ en la semana actual es de 75 unidades y de $X_2$ es de 95 unidades. La política de la compañía es maximizar la suma combinada de unidades de $X_1$ y $X_2$ en inventario al finalizar la semana.

Formular el problema de decidir cuánto hacer de cada producto en la semana como un problema de programación lineal.

Solución

Sean:

$x_1$ la cantidad de unidades de $X_1$ a ser producidas en la semana, y
$x_2$ la cantidad de unidades de $X_2$ a ser producidas en la semana.

Notar que lo que se quiere es maximizar $(x_1-75+30)+(x_2-95+90)=x_1+x_2-50$. Equivalentemente $x_1+x_2$.

Restricciones:

El tiempo de uso disponible de la máquina $A$ es de 40 horas: $50x_1+24x_2\leq 40(60)\Rightarrow 50x_1+24x_2\leq 2400$.
El tiempo de uso disponible de la máquina $B$ es de 35 horas: $30x_1+33x_2\leq 35(60)\Rightarrow 30x_1+33x_2\leq 2100$.
La demanda para $X_1$ en la semana actual es de 75 unidades: $x_1+30\geq 75\Rightarrow x_1\geq 45\Rightarrow -x_1\leq -45$.
La demanda para $X_2$ en la semana actual es de 95 unidades: $x_2+90\geq 95\Rightarrow x_2\geq 5\Rightarrow -x_2\leq -5$.

Finalmente, el problema puede ser expresado en la forma explicada como: \begin{equation} \begin{array}{ll} \min_{x_1,x_2} & -x_1-x_2 \\ \text{s. a. } & 0x_1+0x_2=0 \\ & 50x_1+24x_2\leq 2400 \\ & 30x_1+33x_2\leq 2100 \\ & -x_1\leq -45 \\ & -x_2\leq -5, \end{array} \end{equation}

o, eqivalentemente \begin{equation} \begin{array}{ll} \min_{\boldsymbol{x}} & \boldsymbol{f}^T\boldsymbol{x} \\ \text{s. a. } & \boldsymbol{A}_{eq}\boldsymbol{x}=\boldsymbol{b}_{eq} \\ & \boldsymbol{A}\boldsymbol{x}\leq\boldsymbol{b}, \end{array} \end{equation} con

$\boldsymbol{f}=\left[-1 \quad -1\right]^T$,
$\boldsymbol{A}_{eq}=\left[0\quad 0\right]$,
$\boldsymbol{A}=\left[\begin{array}{cc}50 & 24 \\ 30 & 33\\ -1 & 0\\ 0 & -1\end{array}\right]$,
$\boldsymbol{b}_{eq}=0$, y
$\boldsymbol{b}=\left[2400\quad 2100\quad -45\quad -5\right]^T$.



In [1]:

    
import numpy as np



In [2]:

    
f = np.array([-1, -1])
A = np.array([[50, 24], [30, 33], [-1, 0], [0, -1]])
b = np.array([2400, 2100, -45, -5])



In [3]:

    
import pyomo_utilities



In [5]:

    
x, obj = pyomo_utilities.linprog(f, A, b)









    



# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: -51.25
  Upper bound: -51.25
  Number of objectives: 1
  Number of constraints: 5
  Number of variables: 3
  Number of nonzeros: 7
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.14040040969848633
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0



In [6]:

    
x









    Out[6]:





array([ 45.  ,   6.25])



In [16]:

    
obj









    Out[16]:





-51.25000000000002



In [7]:

    
obj_real = x.sum()-50
obj_real.round(2)









    Out[7]:





1.25

3.2

Mónica hace aretes y cadenitas de joyería. Es tan buena, que todo lo que hace lo vende.

Le toma 30 minutos hacer un par de aretes y una hora hacer una cadenita, y como Mónica también es estudihambre, solo dispone de 10 horas a la semana para hacer las joyas. Por otra parte, el material que compra solo le alcanza para hacer 15 unidades (el par de aretes cuenta como unidad) de joyas por semana.

La utilidad que le deja la venta de las joyas es \$15 en cada par de aretes y \$20 en cada cadenita.

¿Cuántos pares de aretes y cuántas cadenitas debería hacer Mónica para maximizar su utilidad?

Formular el problema en la forma explicada y obtener la solución gráfica.

Solución

Sean:

$x_1$ la cantidad de pares de aretes que hace Mónica.
$x_2$ la cantidad de cadenitas que hace Mónica.

Notar que lo que se quiere es maximizar $15x_1+20x_2$.

Restricciones:

Dispone de 10 horas semanales: $0.5x_1+x_2\leq 10$.
Material disponible para 15 unidades: $x_1+x_2\leq 15$.
No negatividad: $x_1,x_2\geq 0\Rightarrow -x_1,-x_2\leq 0$.

Finalmente, el problema puede ser expresado en la forma explicada como: \begin{equation} \begin{array}{ll} \min_{x_1,x_2} & -15x_1-20x_2 \\ \text{s. a. } & 0x_1+0x_2=0 \\ & 0.5x_1+x_2\leq 10 \\ & x_1+x_2\leq 15 \\ & -x_1\leq 0 \\ & -x_2\leq 0, \end{array} \end{equation}

o, eqivalentemente

con

$\boldsymbol{f}=\left[-15 \quad -20\right]^T$,
$\boldsymbol{A}_{eq}=\left[0\quad 0\right]$,
$\boldsymbol{A}=\left[\begin{array}{cc}0.5 & 1 \\ 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{array}\right]$,
$\boldsymbol{b}_{eq}=0$, y
$\boldsymbol{b}=\left[10\quad 15 \quad 0 \quad 0\right]^T$.



In [8]:

    
f = np.array([-15, -20])
A = np.array([[0.5, 1], [1, 1], [-1, 0], [0, -1]])
b = np.array([10, 15, 0, 0])



In [9]:

    
x, obj = pyomo_utilities.linprog(f, A, b)









    



# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: -250.0
  Upper bound: -250.0
  Number of objectives: 1
  Number of constraints: 5
  Number of variables: 3
  Number of nonzeros: 7
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.2652006149291992
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0



In [27]:

    
x









    Out[27]:





array([ 10.,   5.])



In [28]:

    
obj









    Out[28]:





-250.0



In [10]:

    
obj_real = -obj
obj_real









    Out[10]:





250.0

4. Problema de transporte

Referencia: https://es.wikipedia.org/wiki/Programaci%C3%B3n_lineal

Este es un caso curioso, con solo 6 variables (un caso real de problema de transporte puede tener fácilmente más de 1.000 variables) en el cual se aprecia la utilidad de este procedimiento de cálculo.

Existen tres minas de carbón cuya producción diaria es:

la mina "a" produce 40 toneladas de carbón por día;
la mina "b" produce 40 t/día; y,
la mina "c" produce 20 t/día.

En la zona hay dos centrales termoeléctricas que consumen:

la central "d" consume 40 t/día de carbón; y,
la central "e" consume 60 t/día.

Los costos de mercado, de transporte por tonelada son:

de "a" a "d" = 2 monedas;
de "a" a "e" = 11 monedas;
de "b" a "d" = 12 monedas;
de "b" a "e" = 24 monedas;
de "c" a "d" = 13 monedas; y,
de "c" a "e" = 18 monedas.

Si se preguntase a los pobladores de la zona cómo organizar el transporte, tal vez la mayoría opinaría que debe aprovecharse el precio ofrecido por el transportista que va de "a" a "d", porque es más conveniente que los otros, debido a que es el de más bajo precio.

En este caso, el costo total del transporte es:

transporte de 40 t de "a" a "d" = 80 monedas;
transporte de 20 t de "c" a "e" = 360 monedas; y,
transporte de 40 t de "b" a "e" = 960 monedas,

Para un total 1.400 monedas.

Sin embargo, formulando el problema para ser resuelto por la programación lineal con

$x_1$ toneladas transportadas de la mina "a" a la central "d"
$x_2$ toneladas transportadas de la mina "a" a la central "e"
$x_3$ toneladas transportadas de la mina "b" a la central "d"
$x_4$ toneladas transportadas de la mina "b" a la central "e"
$x_5$ toneladas transportadas de la mina "c" a la central "d"
$x_6$ toneladas transportadas de la mina "c" a la central "e"

se tienen las siguientes ecuaciones:

Restricciones de la producción:

$x_1 + x_2 \leq 40$
$x_3 + x_4 \leq 40$
$x_5 + x_6 \leq 20$

Restricciones del consumo:

$x_1 + x_3 + x_5 \geq 40$
$x_2 + x_4 + x_6 \geq 60$

La función objetivo será:

$$\min_{x_1,\dots,x_6}2x_1 + 11x_2 + 12x_3 + 24x_4 + 13x_5 + 18x_6$$



In [13]:

    
f = np.array([2, 11, 12, 24, 13, 18])
A = np.array([[1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1], [-1, 0, -1, 0, -1, 0], [0, -1, 0, -1, 0, -1]])
A = np.concatenate((A, -np.eye(6)), axis = 0)
b = np.array([40, 40, 20, -40, -60])
b = np.concatenate((b, np.zeros((6,))))



In [14]:

    
x, obj = pyomo_utilities.linprog(f, A, b)









    



# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 1280.0
  Upper bound: 1280.0
  Number of objectives: 1
  Number of constraints: 12
  Number of variables: 7
  Number of nonzeros: 19
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.15600037574768066
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0



In [15]:

    
x









    Out[15]:





array([ -0.,  40.,  40.,  -0.,  -0.,  20.])



In [17]:

    
obj









    Out[17]:





1280.0

La solución de costo mínimo de transporte diario resulta ser:

$x_2 = 40$ resultando un costo de $11(40) = 480$ monedas
$x_3 = 40$ resultando un costo de $12(40) = 440$ monedas
$x_6 = 20$ resultando un costo de $18(20) = 360$ monedas

para un total de $1280$ monedas, $120$ monedas menos que antes.

5. Problema de transporte más interesante (Actividad)

Referencia: https://relopezbriega.github.io/blog/2017/01/18/problemas-de-optimizacion-con-python/

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:



In [42]:

    
import pandas as pd



In [43]:

    
info = pd.DataFrame({'Bar1': [2, 3], 'Bar2': [4, 1], 'Bar3': [5, 3], 'Bar4': [2, 2], 'Bar5': [1, 3]}, index = ['CerveceriaA', 'CerveceriaB'])

info









    Out[43]:







  
    
      
      Bar1
      Bar2
      Bar3
      Bar4
      Bar5
    
  
  
    
      CerveceriaA
      2
      4
      5
      2
      1
    
    
      CerveceriaB
      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 (en cajas de cerveza) y cada bar, las cuales se detallan en el gráfico de más arriba.

Solución

Sean:

$x_i$ cajas transportadas de la cervecería A al Bar $i$,
$x_{i+5}$ cajas transportadas de la cervecería B al Bar $i$.

La actividad consiste en plantear el problema y resolverlo con linprog...



In [45]:

    
f = np.array([2, 4, 5, 2, 1, 3, 1, 3, 2, 3])
A = np.array([[1, 1, 1, 1, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 1, 1, 1, 1], [-1, 0, 0, 0, 0, -1, 0, 0, 0, 0], [0, -1, 0, 0, 0, 0, -1, 0, 0, 0], [0, 0, -1, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, -1, 0, 0, 0, 0, -1, 0], [0, 0, 0, 0, -1, 0, 0, 0, 0, -1]])
A = np.concatenate((A, -np.eye(10)), axis = 0)
b = np.array([1000, 4000, -500, -900, -1800, -200, -700])
b = np.concatenate((b, np.zeros((10,))))



In [46]:

    
x, obj = pyomo_utilities.linprog(f, A, b)









    



# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 8600.0
  Upper bound: 8600.0
  Number of objectives: 1
  Number of constraints: 18
  Number of variables: 11
  Number of nonzeros: 31
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.20701193809509277
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0



In [47]:

    
x









    Out[47]:





array([  300.,    -0.,    -0.,    -0.,   700.,   200.,   900.,  1800.,
         200.,    -0.])



In [48]:

    
obj









    Out[48]:





8600.0