# Programación lineal



> La programación lineal es el campo de la optimización matemática dedicado a maximizar o minimizar (optimizar) funciones lineales, denominada función objetivo, de tal forma que las variables de dicha función estén sujetas a una serie de restricciones expresadas mediante un sistema de ecuaciones o inecuaciones también lineales.

**Referencias:**
- https://es.wikipedia.org/wiki/Programaci%C3%B3n_lineal
- https://docs.scipy.org/doc/scipy/reference/

## 1. Apuntes históricos



- 1826:	Joseph Fourier anticipa la programación lineal. Carl Friedrich Gauss resuelve ecuaciones lineales por eliminación "gaussiana".
- 1902:	Gyula Farkas concibe un método para resolver sistemas de inecuaciones.
- Es hasta la Segunda Guerra Mundial que se plantea la programación lineal como un modelo matemático para planificar gastos y retornos, de modo que se reduzcan costos de guerra y aumentar pérdidas del enemigo. Secreto hasta 1947 (posguerra).
- 1947:	George Dantzig publica el algoritmo simplex y John von Neumann desarrolló la teoría de la dualidad. Se sabe que Leonid Kantoróvich también formuló la teoría en forma independiente.
- Fue usado por muchas industrias en la planificación diaria.

**Hasta acá, tiempos exponenciales de solución. Lo siguiente, tiempo polinomial.**

- 1979: Leonid Khachiyan, diseñó el llamado Algoritmo del elipsoide, a través del cual demostró que el problema de la programación lineal es resoluble de manera eficiente, es decir, en tiempo polinomial.
- 1984: Narendra Karmarkar introduce el método del punto interior para resolver problemas de programación lineal.

**Mencionar complejidad computacional.**

## 2. Motivación

Ya la clase pasada habíamos mencionado que cuando se quería optimizar una función de varias variables con restricciones, se podía aplicar siempre el método de Multiplicadores de Lagrange. Sin embargo, este método es computacionalmente muy complejo conforme crece el número de variables.

Por tanto, cuando la función a optimizar y las restricciones son de caracter lineal, los métodos de solución que se pueden desarrollar son computacionalmente eficientes, por lo que es útil realizar la distinción.

## 3. Problemas de programación lineal

### 3.1. Ejemplo básico
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$ e $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+x_2$.

Restricciones:
1. 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$.
2. 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$.
3. 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$.
4. 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. }  & 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{c}^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{c}=\left[-1 \quad -1\right]^T$,
- $\boldsymbol{A}=\left[\begin{array}{cc}50 & 24 \\ 30 & 33\\ -1 & 0\\ 0 & -1\end{array}\right]$, y
- $\boldsymbol{b}=\left[2400\quad 2100\quad -45\quad -5\right]^T$.

Preferiremos, en adelante, la notación vectorial/matricial.

### 3.2. En general
De acuerdo a lo descrito anteriormente, un problema de programación lineal puede escribirse en la siguiente forma:

\begin{equation}
\begin{array}{ll}
\min_{x_1,\dots,x_n} & c_1x_1+\dots+c_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,
- $c_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{c}^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{c}=\left[c_1\quad\dots\quad c_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})$.

#### Bueno, y una vez planteado, ¿cómo se resuelve el problema?

Este problema está sencillo pues solo es en dos variables. La solución gráfica es válida.


```python
# Importar el módulo pyplot de la librería matplotlib, y la librería numpy
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
```


```python
# Definir las funciones frontera de la región factible
def f1(x1):
    return (2400 - 50 * x1) / 24
def f2(x1):
    return (2100 - 30 * x1) / 33
```


```python
# Vector independiente x1 y funciones evaluadas
x1 = np.linspace(0, 100)
y1 = f1(x1)
y2 = f2(x1)
```


```python
# Gráfica de la región factible y esquinas
plt.figure(figsize=(6, 4))
plt.plot(x1, y1, 'b', lw=2, label='$50x_1 + 24x_2 = 2400$')
plt.plot(x1, y2, 'r', lw=2, label='$30x_1 + 33x_2 = 2100$')
plt.axvline(x=45, c='y', label='$x_1=45$')
plt.axhline(y=5, c='c', label='$x_2=5$')
plt.text(45, 5, '$(45, 5)$')
plt.text(45, (2400 - 50 * 45) / 24, f'$(45, {(2400 - 50 * 45) / 24})$')
plt.text((2400 - 24 * 5) / 50, 5, f'$({(2400 - 24 * 5) / 50}, 5)$')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.grid()
plt.legend(loc=2)
plt.axis([44, 46, 4, 7])
```


```python
45 + 5, 45 + 6.25, 45.6 + 5
```




    (50, 51.25, 50.6)



De acuerdo a la solución gráfica, debemos producir:

- 45 unidades de $X_1$, y
- 6.25 unidades de $X_2$.

Además, el inventario total de productos $X_1$ y $X_2$ al final de la semana es de 1.25 unidades.

**Actividad.** 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 (puede ser a mano).

**Diez minutos: quien primero lo haga, y me lo mande bien explicado al correo, tiene 20 puntos para el examen final.**

## 5. ¿Cómo se resuelve en python?

### 5.1 Librería `SciPy`



`SciPy` es un software de código abierto basado en `Python` para matemáticas, ciencia e ingeniería. 

En particular, los siguientes son algunos de los paquetes básicos:
- `NumPy`
- **Librería `SciPy`**
- `SymPy`
- `matplotlib`
- `pandas`

La **Librería `SciPy`** es uno de los paquetes principales y provee varias rutinas numéricas eficientes. Entre ellas, para integración numérica y optimización.

En esta clase, y en lo que resta del módulo, estaremos utilizando el módulo `optimize` de la librería `SciPy`.

**Importémoslo**


```python
# Importar el módulo optimize de la librería scipy
# import scipy.optimize as opt
from scipy import optimize as opt
```

El módulo `optimize` que acabamos de importar contiene varias funciones para optimización y búsqueda de raices ($f(x)=0$). Entre ellas se encuentra la función `linprog`


```python
# Función linprog del módulo optimize
help(opt.linprog)
```

    Help on function linprog in module scipy.optimize._linprog:
    
    linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='interior-point', callback=None, options=None, x0=None)
        Linear programming: minimize a linear objective function subject to linear
        equality and inequality constraints.
        
        Linear programming solves problems of the following form:
        
        .. math::
        
            \min_x \ & c^T x \\
            \mbox{such that} \ & A_{ub} x \leq b_{ub},\\
            & A_{eq} x = b_{eq},\\
            & l \leq x \leq u ,
        
        where :math:`x` is a vector of decision variables; :math:`c`,
        :math:`b_{ub}`, :math:`b_{eq}`, :math:`l`, and :math:`u` are vectors; and
        :math:`A_{ub}` and :math:`A_{eq}` are matrices.
        
        Informally, that's:
        
        minimize::
        
            c @ x
        
        such that::
        
            A_ub @ x <= b_ub
            A_eq @ x == b_eq
            lb <= x <= ub
        
        Note that by default ``lb = 0`` and ``ub = None`` unless specified with
        ``bounds``.
        
        Parameters
        ----------
        c : 1D array
            The coefficients of the linear objective function to be minimized.
        A_ub : 2D array, optional
            The inequality constraint matrix. Each row of ``A_ub`` specifies the
            coefficients of a linear inequality constraint on ``x``.
        b_ub : 1D array, optional
            The inequality constraint vector. Each element represents an
            upper bound on the corresponding value of ``A_ub @ x``.
        A_eq : 2D array, optional
            The equality constraint matrix. Each row of ``A_eq`` specifies the
            coefficients of a linear equality constraint on ``x``.
        b_eq : 1D array, optional
            The equality constraint vector. Each element of ``A_eq @ x`` must equal
            the corresponding element of ``b_eq``.
        bounds : sequence, optional
            A sequence of ``(min, max)`` pairs for each element in ``x``, defining
            the minimum and maximum values of that decision variable. Use ``None`` to
            indicate that there is no bound. By default, bounds are ``(0, None)``
            (all decision variables are non-negative).
            If a single tuple ``(min, max)`` is provided, then ``min`` and
            ``max`` will serve as bounds for all decision variables.
        method : {'interior-point', 'revised simplex', 'simplex'}, optional
            The algorithm used to solve the standard form problem.
            :ref:`'interior-point' <optimize.linprog-interior-point>` (default),
            :ref:`'revised simplex' <optimize.linprog-revised_simplex>`, and
            :ref:`'simplex' <optimize.linprog-simplex>` (legacy)
            are supported.
        callback : callable, optional
            If a callback function is provided, it will be called at least once per
            iteration of the algorithm. The callback function must accept a single
            `scipy.optimize.OptimizeResult` consisting of the following fields:
        
                x : 1D array
                    The current solution vector.
                fun : float
                    The current value of the objective function ``c @ x``.
                success : bool
                    ``True`` when the algorithm has completed successfully.
                slack : 1D array
                    The (nominally positive) values of the slack,
                    ``b_ub - A_ub @ x``.
                con : 1D array
                    The (nominally zero) residuals of the equality constraints,
                    ``b_eq - A_eq @ x``.
                phase : int
                    The phase of the algorithm being executed.
                status : int
                    An integer representing the status of the algorithm.
        
                    ``0`` : Optimization proceeding nominally.
        
                    ``1`` : Iteration limit reached.
        
                    ``2`` : Problem appears to be infeasible.
        
                    ``3`` : Problem appears to be unbounded.
        
                    ``4`` : Numerical difficulties encountered.
        
                nit : int
                    The current iteration number.
                message : str
                    A string descriptor of the algorithm status.
        
        options : dict, optional
            A dictionary of solver options. All methods accept the following
            options:
        
                maxiter : int
                    Maximum number of iterations to perform.
                    Default: see method-specific documentation.
                disp : bool
                    Set to ``True`` to print convergence messages.
                    Default: ``False``.
                autoscale : bool
                    Set to ``True`` to automatically perform equilibration.
                    Consider using this option if the numerical values in the
                    constraints are separated by several orders of magnitude.
                    Default: ``False``.
                presolve : bool
                    Set to ``False`` to disable automatic presolve.
                    Default: ``True``.
                rr : bool
                    Set to ``False`` to disable automatic redundancy removal.
                    Default: ``True``.
        
            For method-specific options, see
            :func:`show_options('linprog') <show_options>`.
        
        x0 : 1D array, optional
            Guess values of the decision variables, which will be refined by
            the optimization algorithm. This argument is currently used only by the
            'revised simplex' method, and can only be used if `x0` represents a
            basic feasible solution.
        
        
        Returns
        -------
        res : OptimizeResult
            A :class:`scipy.optimize.OptimizeResult` consisting of the fields:
        
                x : 1D array
                    The values of the decision variables that minimizes the
                    objective function while satisfying the constraints.
                fun : float
                    The optimal value of the objective function ``c @ x``.
                slack : 1D array
                    The (nominally positive) values of the slack variables,
                    ``b_ub - A_ub @ x``.
                con : 1D array
                    The (nominally zero) residuals of the equality constraints,
                    ``b_eq - A_eq @ x``.
                success : bool
                    ``True`` when the algorithm succeeds in finding an optimal
                    solution.
                status : int
                    An integer representing the exit status of the algorithm.
        
                    ``0`` : Optimization terminated successfully.
        
                    ``1`` : Iteration limit reached.
        
                    ``2`` : Problem appears to be infeasible.
        
                    ``3`` : Problem appears to be unbounded.
        
                    ``4`` : Numerical difficulties encountered.
        
                nit : int
                    The total number of iterations performed in all phases.
                message : str
                    A string descriptor of the exit status of the algorithm.
        
        See Also
        --------
        show_options : Additional options accepted by the solvers.
        
        Notes
        -----
        This section describes the available solvers that can be selected by the
        'method' parameter.
        
        :ref:`'interior-point' <optimize.linprog-interior-point>` is the default
        as it is typically the fastest and most robust method.
        :ref:`'revised simplex' <optimize.linprog-revised_simplex>` is more
        accurate for the problems it solves.
        :ref:`'simplex' <optimize.linprog-simplex>` is the legacy method and is
        included for backwards compatibility and educational purposes.
        
        Method *interior-point* uses the primal-dual path following algorithm
        as outlined in [4]_. This algorithm supports sparse constraint matrices and
        is typically faster than the simplex methods, especially for large, sparse
        problems. Note, however, that the solution returned may be slightly less
        accurate than those of the simplex methods and will not, in general,
        correspond with a vertex of the polytope defined by the constraints.
        
        .. versionadded:: 1.0.0
        
        Method *revised simplex* uses the revised simplex method as described in
        [9]_, except that a factorization [11]_ of the basis matrix, rather than
        its inverse, is efficiently maintained and used to solve the linear systems
        at each iteration of the algorithm.
        
        .. versionadded:: 1.3.0
        
        Method *simplex* uses a traditional, full-tableau implementation of
        Dantzig's simplex algorithm [1]_, [2]_ (*not* the
        Nelder-Mead simplex). This algorithm is included for backwards
        compatibility and educational purposes.
        
        .. versionadded:: 0.15.0
        
        Before applying any method, a presolve procedure based on [8]_ attempts
        to identify trivial infeasibilities, trivial unboundedness, and potential
        problem simplifications. Specifically, it checks for:
        
        - rows of zeros in ``A_eq`` or ``A_ub``, representing trivial constraints;
        - columns of zeros in ``A_eq`` `and` ``A_ub``, representing unconstrained
          variables;
        - column singletons in ``A_eq``, representing fixed variables; and
        - column singletons in ``A_ub``, representing simple bounds.
        
        If presolve reveals that the problem is unbounded (e.g. an unconstrained
        and unbounded variable has negative cost) or infeasible (e.g. a row of
        zeros in ``A_eq`` corresponds with a nonzero in ``b_eq``), the solver
        terminates with the appropriate status code. Note that presolve terminates
        as soon as any sign of unboundedness is detected; consequently, a problem
        may be reported as unbounded when in reality the problem is infeasible
        (but infeasibility has not been detected yet). Therefore, if it is
        important to know whether the problem is actually infeasible, solve the
        problem again with option ``presolve=False``.
        
        If neither infeasibility nor unboundedness are detected in a single pass
        of the presolve, bounds are tightened where possible and fixed
        variables are removed from the problem. Then, linearly dependent rows
        of the ``A_eq`` matrix are removed, (unless they represent an
        infeasibility) to avoid numerical difficulties in the primary solve
        routine. Note that rows that are nearly linearly dependent (within a
        prescribed tolerance) may also be removed, which can change the optimal
        solution in rare cases. If this is a concern, eliminate redundancy from
        your problem formulation and run with option ``rr=False`` or
        ``presolve=False``.
        
        Several potential improvements can be made here: additional presolve
        checks outlined in [8]_ should be implemented, the presolve routine should
        be run multiple times (until no further simplifications can be made), and
        more of the efficiency improvements from [5]_ should be implemented in the
        redundancy removal routines.
        
        After presolve, the problem is transformed to standard form by converting
        the (tightened) simple bounds to upper bound constraints, introducing
        non-negative slack variables for inequality constraints, and expressing
        unbounded variables as the difference between two non-negative variables.
        Optionally, the problem is automatically scaled via equilibration [12]_.
        The selected algorithm solves the standard form problem, and a
        postprocessing routine converts the result to a solution to the original
        problem.
        
        References
        ----------
        .. [1] Dantzig, George B., Linear programming and extensions. Rand
               Corporation Research Study Princeton Univ. Press, Princeton, NJ,
               1963
        .. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to
               Mathematical Programming", McGraw-Hill, Chapter 4.
        .. [3] Bland, Robert G. New finite pivoting rules for the simplex method.
               Mathematics of Operations Research (2), 1977: pp. 103-107.
        .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point
               optimizer for linear programming: an implementation of the
               homogeneous algorithm." High performance optimization. Springer US,
               2000. 197-232.
        .. [5] Andersen, Erling D. "Finding all linearly dependent rows in
               large-scale linear programming." Optimization Methods and Software
               6.3 (1995): 219-227.
        .. [6] Freund, Robert M. "Primal-Dual Interior-Point Methods for Linear
               Programming based on Newton's Method." Unpublished Course Notes,
               March 2004. Available 2/25/2017 at
               https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf
        .. [7] Fourer, Robert. "Solving Linear Programs by Interior-Point Methods."
               Unpublished Course Notes, August 26, 2005. Available 2/25/2017 at
               http://www.4er.org/CourseNotes/Book%20B/B-III.pdf
        .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear
               programming." Mathematical Programming 71.2 (1995): 221-245.
        .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear
               programming." Athena Scientific 1 (1997): 997.
        .. [10] Andersen, Erling D., et al. Implementation of interior point
                methods for large scale linear programming. HEC/Universite de
                Geneve, 1996.
        .. [11] Bartels, Richard H. "A stabilization of the simplex method."
                Journal in  Numerische Mathematik 16.5 (1971): 414-434.
        .. [12] Tomlin, J. A. "On scaling linear programming problems."
                Mathematical Programming Study 4 (1975): 146-166.
        
        Examples
        --------
        Consider the following problem:
        
        .. math::
        
            \min_{x_0, x_1} \ -x_0 + 4x_1 & \\
            \mbox{such that} \ -3x_0 + x_1 & \leq 6,\\
            -x_0 - 2x_1 & \geq -4,\\
            x_1 & \geq -3.
        
        The problem is not presented in the form accepted by `linprog`. This is
        easily remedied by converting the "greater than" inequality
        constraint to a "less than" inequality constraint by
        multiplying both sides by a factor of :math:`-1`. Note also that the last
        constraint is really the simple bound :math:`-3 \leq x_1 \leq \infty`.
        Finally, since there are no bounds on :math:`x_0`, we must explicitly
        specify the bounds :math:`-\infty \leq x_0 \leq \infty`, as the
        default is for variables to be non-negative. After collecting coeffecients
        into arrays and tuples, the input for this problem is:
        
        >>> c = [-1, 4]
        >>> A = [[-3, 1], [1, 2]]
        >>> b = [6, 4]
        >>> x0_bounds = (None, None)
        >>> x1_bounds = (-3, None)
        >>> from scipy.optimize import linprog
        >>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])
        
        Note that the default method for `linprog` is 'interior-point', which is
        approximate by nature.
        
        >>> print(res)
             con: array([], dtype=float64)
             fun: -21.99999984082494 # may vary
         message: 'Optimization terminated successfully.'
             nit: 6 # may vary
           slack: array([3.89999997e+01, 8.46872439e-08] # may vary
          status: 0
         success: True
               x: array([ 9.99999989, -2.99999999]) # may vary
        
        If you need greater accuracy, try 'revised simplex'.
        
        >>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds], method='revised simplex')
        >>> print(res)
             con: array([], dtype=float64)
             fun: -22.0 # may vary
         message: 'Optimization terminated successfully.'
             nit: 1 # may vary
           slack: array([39.,  0.]) # may vary
          status: 0
         success: True
               x: array([10., -3.]) # may vary
    


\begin{align}
\min_x \ & c^T x \\
\mbox{such that} \ & A_{ub} x \leq b_{ub},\\
        & A_{eq} x = b_{eq},\\
        & l \leq x \leq u ,
\end{align}

la cual resuelve problemas como los que aprendimos a plantear.

### 5.2 Solución del ejemplo básico con linprog
Ya hicimos la solución gráfica. Contrastemos con la solución que nos da `linprog`...

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


```python
# Importar numpy para crear las matrices
import numpy as np
```


```python
# Crear las matrices para resolver el problema
c = np.array([-1, -1])
A_ineq = np.array([[50, 24],
              [30, 33],
              [-1,  0],
              [0,  -1]])
b_ineq = np.array([2400, 2100, -45, -5])
```


```python
c, A_ineq, b_ineq
```




    (array([-1, -1]),
     array([[50, 24],
            [30, 33],
            [-1,  0],
            [ 0, -1]]),
     array([2400, 2100,  -45,   -5]))




```python
# Resolver utilizando linprog
solucion = opt.linprog(c=c, A_ub=A_ineq, b_ub=b_ineq)
```


```python
# Mostrar el resultado
solucion
```




         con: array([], dtype=float64)
         fun: -51.249999998664535
     message: 'Optimization terminated successfully.'
         nit: 8
       slack: array([1.43763828e-08, 5.43750000e+02, 6.79790446e-10, 1.25000000e+00])
      status: 0
     success: True
           x: array([45.  ,  6.25])




```python
# Extraer el vector solución
solucion.x
```




    array([45.  ,  6.25])



**Conclusión**
- 45 unidades de $X_1$, y
- 6.25 unidades de $X_2$.

Además, el inventario total de productos $X_1$ y $X_2$ al final de la semana es de 1.25 unidades.

**Otra forma:** poner las cotas de las variables a parte


```python
# Escribir matrices y cotas
c = np.array([-1, -1])
A_ineq = np.array([[50, 24],
                   [30, 33]])
b_ineq = np.array([2400, 2100])
# Cotas
bnds = ((45, None),  # cotas para x1: 45 <= x1 <= infinto
        (5, None))   # cotas para x2: 5 <= x2 <= infinto
```


```python
# Resolver
solucion2 = opt.linprog(c=c, A_ub=A_ineq, b_ub=b_ineq, bounds=bnds)
```


```python
# Mostrar el resultado
solucion2
```




         con: array([], dtype=float64)
         fun: -51.249999998664535
     message: 'Optimization terminated successfully.'
         nit: 8
       slack: array([1.43763828e-08, 5.43750000e+02])
      status: 0
     success: True
           x: array([45.  ,  6.25])



**Actividad.** Resolver el ejemplo de Mónica y sus tiliches con `linprog`


```python
# Resolver acá

```


```python

```


```python

```

## 6. Problema de transporte 1
- **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$$

Resolver con `linprog`


```python
# Matrices y cotas
c = np.array([2, 11, 12, 24, 13, 18])
A_ineq = 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]])
b_ineq = np.array([40, 40 ,20 , -40, -60])
# bnds = ((0, None), (0, None), (0, None), (0, None), (0, None), (0, None))
bnds = ((0, None),) * 6
```


```python
# Resolver
solucion = opt.linprog(c=c, A_ub=A_ineq, b_ub=b_ineq, bounds=bnds)
```


```python
opt.linprog(c=c, A_ub=A_ineq, b_ub=b_ineq)
```




         con: array([], dtype=float64)
         fun: 1279.999999830141
     message: 'Optimization terminated successfully.'
         nit: 5
       slack: array([ 5.41090373e-09,  5.60606850e-09,  2.64986610e-09, -5.41140821e-09,
           -8.25542656e-09])
      status: 0
     success: True
           x: array([1.85251940e-09, 4.00000000e+01, 4.00000000e+01, 1.85923309e-09,
           2.01371561e-10, 2.00000000e+01])




```python
solucion
```




         con: array([], dtype=float64)
         fun: 1279.999999830141
     message: 'Optimization terminated successfully.'
         nit: 5
       slack: array([ 5.41090373e-09,  5.60606850e-09,  2.64986610e-09, -5.41140821e-09,
           -8.25542656e-09])
      status: 0
     success: True
           x: array([1.85251940e-09, 4.00000000e+01, 4.00000000e+01, 1.85923309e-09,
           2.01371561e-10, 2.00000000e+01])




```python
# Mostrar resultado
x_min = solucion.x
np.round(x_min)
```




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




```python
c.T.dot(x_min)
```




    1279.999999830141




```python
solucion.fun
```




    1279.999999830141



**Conclusión**

- Se deben transportar 40 ton. de la mina "a" a la central "e" ($x_2$)
- Se deben transportar 40 ton. de la mina "b" a la central "d" ($x_3$)
- Se deben transportar 20 ton. de la mina "c" a la central "e" ($x_6$)

## 7. Optimización de inversión en bonos

**Referencia:**


```python
from IPython.display import YouTubeVideo
YouTubeVideo('gukxBus8lOs')
```










El objetivo de este problema es determinar la mejor estrategia de inversión, dados diferentes tipos de bono, la máxima cantidad que puede ser invertida en cada bono, el porcentaje de retorno y los años de madurez. También hay una cantidad fija de dinero disponible ($\$750,000$). Por lo menos la mitad de este dinero debe ser invertido en bonos con 10 años o más para la madurez. Se puede invertir un máximo del $25\%$ de esta cantidad en cada bono. Finalmente, hay otra restricción que no permite usar más de $35\%$ en bonos de alto riesgo.

Existen seis (6) opciones de inversión con las letras correspondientes $A_i$

1. $A_1$:(Tasa de retorno=$8.65\%$; Años para la madurez=11, Riesgo=Bajo)
1. $A_2$:(Tasa de retorno=$9.50\%$; Años para la madurez=10, Riesgo=Alto)
1. $A_3$:(Tasa de retorno=$10.00\%$; Años para la madurez=6, Riesgo=Alto)
1. $A_4$:(Tasa de retorno=$8.75\%$; Años para la madurez=10, Riesgo=Bajo)
1. $A_5$:(Tasa de retorno=$9.25\%$; Años para la madurez=7, Riesgo=Alto)
1. $A_6$:(Tasa de retorno=$9.00\%$; Años para la madurez=13, Riesgo=Bajo)

Lo que se quiere entonces es maximizar el retorno que deja la inversión.

Este problema puede ser resuelto con programación lineal. Formalmente, puede ser descrito como:

$$\max_{A_1,A_2,...,A_6}\sum^{6}_{i=1} A_iR_i,$$

donde $A_i$ representa la cantidad invertida en la opción, y $R_i$ representa la tasa de retorno respectiva.

Plantear restricciones...


```python
# Matrices y cotas

```


```python
# Resolver

```


```python
# Mostrar resultado

```

Recordar que en el problema minimizamos $-\sum^{6}_{i=1} A_iR_i$. El rendimiento obtenido es entonces:


```python

```

**Conclusión**

- 

## 8. Tarea 

### 1. Diseño de la Dieta Óptima

Se quiere producir comida para gatos de la manera más barata, no obstante se debe también asegurar que se cumplan los datos requeridos de analisis nutricional. Por lo que se quiere variar la cantidad de cada ingrediente para cumplir con los estandares nutricionales. Los requisitos que se tienen es que en 100 gramos, se deben tener por lo menos 8 gramos de proteína y 6 gramos de grasa. Así mismo, no se debe tener más de 2 gramos de fibra y 0.4 gramos de sal.  

Los datos nutricionales se pueden obtener de la siguiente tabla:

|Ingrediente|Proteína|Grasa|Fibra|Sal|
|----|---- |---- |---- |---- |
|Pollo|  10.0%|08.0%|00.1%|00.2%
|Carne|  20.0%|10.0%|00.5%|00.5%
|Cordero|15.0%|11.0%|00.5%|00.7%
|Arroz|  00.0%|01.0%|10.0%|00.2%
|Trigo|  04.0%|01.0%|15.0%|00.8%
|Gel|    00.0%|00.0%|00.0%|00.0%

Los costos de cada producto son:

Ingrediente|Costo por gramo
:----|----
Pollo|$\$$0.013
Carne|$\$$0.008
Cordero|$\$$0.010
Arroz|$\$$0.002
Trigo|$\$$0.005
Gel|$\$$0.001    

Lo que se busca optimizar en este caso es la cantidad de productos que se debe utilizar en la comida de gato, para simplificar la notación se van a nombrar las siguientes variables: 

$x_1:$ Gramos de pollo  
$x_2:$ Gramos de carne  
$x_3:$ Gramos de cordero  
$x_4:$ Gramos de arroz  
$x_5:$ Gramos de trigo  
$x_6:$ Gramos de gel  

Con los datos, se puede plantear la función objetivo, está dada por la siguiente expresión:

$$\min 0.013 x_1 + 0.008 x_2 + 0.010 x_3 + 0.002 x_4 + 0.005 x_5 + 0.001 x_6$$

Las restricciones estarían dadas por el siguiente conjunto de ecuaciones:

$x_1+x_2+x_3+x_4+x_5+x_6=100$  

$(10.0 x_1+ 20.0 x_2+ 15.0 x_3+ 00.0 x_4+ 04.0 x_5+ 00.0 x_6)/100 \geq 8.0$  

$(08.0 x_1+ 10.0 x_2+ 11.0 x_3+ 01.0 x_4+ 01.0 x_5+ 00.0 x_6)/100 \geq 6.0$  

$(00.1 x_1+ 00.5 x_2+ 00.5 x_3+ 10.0 x_4+ 15.0 x_5+ 00.0 x_6)/100 \leq 2.0$  

$(00.2 x_1+ 00.5 x_2+ 00.7 x_3+ 00.2 x_4+ 00.8 x_5+ 00.0 x_6)/100 \leq 0.4$  

La primer condición asegura que la cantidad de productos que se usará cumple con los 100 gramos. Las siguientes sólo siguen los lineamientos planteados para cumplir con los requisitos nutrimentales. 

### 2. Otro problema de transporte

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 (Modelo y Cuauhtémoc Moctezuma) 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:


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




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Bar1</th>
      <th>Bar2</th>
      <th>Bar3</th>
      <th>Bar4</th>
      <th>Bar5</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>CerveceriaA</th>
      <td>2</td>
      <td>4</td>
      <td>5</td>
      <td>2</td>
      <td>1</td>
    </tr>
    <tr>
      <th>CerveceriaB</th>
      <td>3</td>
      <td>1</td>
      <td>3</td>
      <td>2</td>
      <td>3</td>
    </tr>
  </tbody>
</table>
</div>



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.

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 tarea consiste en plantear el problema de minimizar el costo de transporte de la forma vista y resolverlo con `linprog`.

Deben crear un notebook de jupyter (archivo .ipynb) y llamarlo Tarea4_ApellidoNombre, y subirlo a moodle.

**Definir fecha**



<footer id="attribution" style="float:right; color:#808080; background:#fff;">
Created with Jupyter by Esteban Jiménez Rodríguez.
</footer>

ipynb
e0dff2c36bb2b72b230384986d8d225bceb6a172
Jupyter Notebook
1184b25ed310d1f2c328673b95c542c1e45166d5
MIT
eduardoxsalazar/simmatp2021
Modulo1/Clase3_ProgramacionLineal.ipynb
1184b25ed310d1f2c328673b95c542c1e45166d5
MIT
eduardoxsalazar/simmatp2021
Modulo1/Clase3_ProgramacionLineal.ipynb
1184b25ed310d1f2c328673b95c542c1e45166d5
MIT
eduardoxsalazar/simmatp2021
Modulo1/Clase3_ProgramacionLineal.ipynb