text:
###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth. Partly based on David Ketcheson's pendulum lesson, also under CC-BY.

# Phugoid Oscillation

Welcome back! This is the second IPython Notebook of the series _"The phugoid model of glider flight"_, the first learning module of the course [**"Practical Numerical Methods with Python."**](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about)

In the first notebook, _"Phugoid Motion"_, we described the physics of an aircraft's oscillatory trajectory, seen as an exchange of kinetic and potential energy. This analysis goes back to Frederick Lanchester, who published his book _"Aerodonetics"_ on aircraft stability in 1909. We concluded that first exposure to our problem of interest by plotting the flight paths predicted by Lanchester's analysis, known as _phugoids_.

Here, we will look at the situation when an aircraft is initially moving on the straight-line phugoid (obtained with the parameters $C=2/3$, $\cos\theta=1$, and $z=z_t$ in the previous analysis), and experiences a small upset, a wind gust that slightly perturbs its path. It will then enter into a gentle oscillation around the previous straight-line path: a _phugoid oscillation_.

If the aircraft experiences an upward acceleration of $-d^2z/dt^2$, and we assume that the perturbation is small, then $\cos\theta=1$ is a good approximation and Newton's second law in the vertical direction is:

$$\begin{equation}
L - W = - \frac{W}{g}\frac{d^2 z}{dt^2}
\end{equation}$$

In the previous notebook, we saw that the following relation holds for the ratio of lift to weight, in terms of the trim velocity $v_t$:

$$\begin{equation}
\frac{L}{W}=\frac{v^2}{v_t^2}
\end{equation}$$

This will be useful: we can divide Equation (1) by the weight and use Equation (2) to replace $L/W$. Another useful relation from the previous notebook expressed the conservation of energy (per unit mass) as $v^2 = 2 gz$. With this, Equation (1) is rearranged as:

$$\begin{equation}
\frac{d^2z}{dt^2} + \frac{gz}{z_t} = g
\end{equation}$$

Look at Equation (3) for a moment. Does it ring a bell? Do you recognize it?

If you remember from your physics courses the equation for _simple harmonic motion_, you should see the similarity!

Take the case of a simple spring. Hooke's law is $F=-kx$, where $F$ is a restoring force, $x$ the displacement from a position of equilibrium and $k$ the spring constant. This results in the following ordinary differential equation for the displacement:

$$\begin{equation}
m \frac{d^2 x}{dt^2}= -\frac{k}{m}x
\end{equation}$$

which has the solution $x(t) = A \cos(\omega t- \phi)$, representing simple harmonic motion with an angular frequency $\omega=\sqrt{k/m}=2\pi f$ and phase angle $\phi$.

Now look back at Equation (3): it has the same form and it represents simple harmonic motion with angular frequency $\omega=\sqrt{g/z_t}$.

Think about this for a moment ... we can immediately say what the period of the oscillation is: exactly $2 \pi \sqrt{z_t/g}$ — or, in terms of the trim velocity, $\pi \sqrt{2} v_t/g$.

_This is a remarkable result!_ Think about it: we know nothing about the aircraft, or the flight altitude, yet we can obtain the period of the phugoid oscillation simply as a function of the trim velocity. For example, if trim velocity is 200 knots, we get a phugoid period of about 47 seconds—over that time, you really would not notice anything if you were flying in that aircraft.

Next, we want to be able to compute the trajectory of the aircraft for a given initial perturbance. We will do this by numerically integrating the equation of motion.

## Prepare to integrate

We want to integrate the differential equation and plot the trajectory of the aircraft. Are you ready?

The equation for the phugoid oscillation is a second-order, ordinary differential equation (ODE). Let's represent the time derivative with a prime, and write it like this:

$$\begin{equation}
z(t)'' + \frac{g \,z(t)}{z_t}=g
\end{equation}$$

There's a convenient trick when we work with ODEs: we can turn this 2nd-order equation into a system of two 1st-order equations. Like this:

$$\begin{eqnarray}
z'(t) &=& v(t)\\
v'(t) &=& g\left(1-\frac{z(t)}{z_t}\right)
\end{eqnarray}$$

Are you following? Make sure you are following the derivations, even if it means writing the equations down in your own notes! (Yes, the old-fashioned paper way.)

Another way to look at a system of two 1st-order ODEs is by using vectors. You can make a vector with your two independent variables,

$$\begin{equation}
\vec{u}  = \begin{pmatrix} z \\ v \end{pmatrix}
\end{equation}$$

and write the differential system as a single vector equation:

\begin{equation}
\vec{u}'(t)  = \begin{pmatrix} v \\ g-g\frac{z(t)}{z_t} \end{pmatrix}
\end{equation}

If you call the right-hand-side $\vec{f}(\vec{u})$, then the equation is very short: $\vec{u}'(t) = \vec{f}(\vec{u})$—but let's drop those arrows to denote vectors from now on, as they are a bit cumbersome: just remember that $u$ and $f$ are vectors in the phugoid equation of motion.

Next, we'll prepare to solve this problem numerically.

## Initial value problems

Let's step back for a moment. Suppose we have a first-order ODE $u'=f(u)$. You know that if we were to integrate this, there would be an arbitrary constant of integration. To find its value, we do need to know one point on the curve $(t, u)$. When the derivative in the ODE is with respect to time, we call that point the _initial value_ and write something like this:

$$u(t=0)=u_0$$

In the case of a second-order ODE, we already saw how to write it as a system of first-order ODEs, and we would need an initial value for each equation: two conditions are needed to determine our constants of integration. The same applies for higher-order ODEs: if it is of order $n$, we can write it as $n$ first-order equations, and we need $n$ known values. If we have that data, we call the problem an _initial value problem_.

Remember the definition of a derivative? The derivative represents the slope of the tangent at a point of the curve $u=u(t)$, and the definition of the derivative $u'$ for a function is:

$$u'(t) = \lim_{\Delta t\rightarrow 0} \frac{u(t+\Delta t)-u(t)}{\Delta t}$$

If the step $\Delta t$ is already very small, we can _approximate_ the derivative by dropping the limit. We can write:

$$\begin{equation}
u(t+\Delta t) \approx u(t) + u'(t) \Delta t
\end{equation}$$

With this equation, and because we know $u'(t)=f(u)$, if we have an initial value, we can step by $\Delta t$ and find the value of $u(t+\Delta t)$, then we can take this value, and find $u(t+2\Delta t)$, and so on: we say that we _step in time_, numerically finding the solution $u(t)$ for a range of values: $t_1, t_2, t_3 \cdots$, each separated by $\Delta t$. The numerical solution of the ODE is simply the table of values $t_i, u_i$ that results from this process.

## Discretization

In order to execute the process described above and find the numerical solution of the ODE, we start by choosing the values $t_1,t_2,t_3 \cdots t_n$—we call these values our *grid* in time. The first point of the grid is given by our _initial value_, and the small difference between two consecutive times is called the _time step_, denoted by $\Delta t$.  The solution value at time $t_n$ is denoted by $u_n$.

Let's build a time grid for our problem. We first choose a final time $T$ and the time step $\Delta t$. In code, we'll use readily identifiable variable names: `T` and `dt`, respectively. With those values set, we can calculate the number of time steps that will be needed to reach the final time; we call that variable `N`.

Let's write some code. The first thing we do in Python is load our favorite libraries: NumPy for array operations, and the Pyplot module in Matplotlib, to later on be able to plot the numerical solution. The line `%matplotlib inline` tells IPython Notebook to show the plots inline.

```
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
```

Now, initialize `T` and `dt`, calculate `N` and build a NumPy array with all the values of time that make up the grid.

```
T = 100.0
dt = 0.01
N = int(T/dt)
t = numpy.linspace(0.0, T, N)
```

We have our grid! Now it's time to apply the numerical time stepping represented by Equation (10).

## Euler's method

The approximate solution at time $t_n$ is $u_n$, and the numerical solution of the differential equation consists of computing a sequence of approximate solutions by the following formula, based on Equation (10):

$$u_{n+1} = u_n + \Delta t \,f(u_n).$$

This formula is called **Euler's method**.

For the equations of the phugoid oscillation, Euler's method gives the following algorithm that we need to implement in code:

\begin{align}
z_{n+1} & = v_n + \Delta t \, z_n \\
v_{n+1} & = g - \frac{g}{z_t} \, z_n + \Delta t \, v_n.
\end{align}

### And solve!

To apply the numerical solution method, we need to set things up in code: define the parameter values needed in the model, initialize a NumPy array to hold the discrete solution values, and initialize another array for the angular values.

\begin{align}
z_{n+1} & = v_n + \Delta t \, z_n \\
v_{n+1} & = g - \frac{g}{z_t} \, z_n + \Delta t \, v_n.
\end{align}

```
# initial conditions
z0 = 100.  #altitude
v  = 10.   #upward velocity resulting from gust
zt = 100.
g  = 9.81

u = numpy.array([z0, v])

# initialize an array to hold the changing angle values
z = numpy.zeros(N)
z[0] = z0
```

You should pay attention to a couple of things: (1) See how there is a dot after the numbers used to define our parameters? We just want to be explicit (as a good habit) that these variables are real numbers, called "floats." (2) We both _created_ and _initialized_ with zeros everywhere the solution vector `z`. Look up the documentation for the handy NumPy function `zeros()`, if you need to. (3) In the last line above, we assign the _initial value_ to the first element of the solution vector: `z[0]`.

Now we can step in time using Euler's method. Notice how we are time stepping the two independend variables at once in the time iterations.

```
# time-loop using Euler's method
for n in range(N-1):
    u = u + dt*numpy.array([u[1], g*(1-u[0]/zt)])
    z[n] = u[0]
```

Make sure you understand what this code is doing. This is a basic pattern in numerical methods: iterations in a time variable that apply a numerical scheme at each step.

## Plot the solution

If the code is correct, we have stored in the array `z` the position of the glider at each time. Let's use Matplotlib to examine the flight path of the arcraft.

You should explore the [Matplotlib tutorial](http://matplotlib.org/users/pyplot_tutorial.html) (if you need to) and familiarize yourself with the command-style functions that control the size, labels, line style, and so on. Creating good plots is a useful skill: it is about communicating your results effectively.

Here, we set the figure size, the limits of the vertical axis, the format of tick-marks, and axis labels. The final line actually produces the plot, with our chosen line style (continuous black line).

```
plt.figure(figsize=(10,4))   #set plot size
plt.ylim(40,150)             #y-axis plot limits
plt.tick_params(axis='both', labelsize=14) #increase font size for ticks
plt.xlabel('t', fontsize=14) #x label
plt.ylabel('z', fontsize=14) #y label
plt.plot(t,z, 'k-');
```

That sure looks like a bumpy ride!

#####Explore and think

Try changing the value of `v` in the initial conditions.

* What happens when you have a larger gust?
* What about a smaller gust?
* What happens if there isn't a gust (`v = 0`)?

## Exact solution

The equation for phugoid oscillations is a 2nd-order, linear ODE and it has an exact solution of the following form:

$$
\begin{equation}
z(t) = A \sin \left(\sqrt{\frac{g}{z_t}} t \right) + B \cos \left(\sqrt{\frac{g}{z_t}} t \right) + z_t
\end{equation}
$$

where $A$ and $B$ are constants that we solve for using initial conditions.

Our numerical solution used the initial conditions:

$$\begin{eqnarray}
z(0) = z_0 \\
v(0) = v_0
\end{eqnarray}$$

Applying these to the exact solution and solving for $A$ and $B$, we get:

$$\begin{equation}
z(t) = v_0 \sqrt{\frac{z_t}{g}} \sin \left(\sqrt{\frac{g}{z_t}} t \right) + (z_0-z_t) \cos \left(\sqrt{\frac{g}{z_t}} t \right) + z_t
\end{equation}$$

We already defined all of these variables for our numerical solution, so we can immediately compute the exact solution.

#####Pro tip:

The expression is a  bit long —if you don't feel like scrolling left and right, you can insert a line break using a backslash '\' and Python will treat the next line as a continuation of the first.

```
z_exact = v*(zt/g)**.5*numpy.sin((g/zt)**.5*t)+\
            (z0-zt)*numpy.cos((g/zt)**.5*t)+zt
```

## Plot the exact solution

Now we can plot our exact solution!  Even better, we can plot _both_ the numerical solution *and* the exact solution to see how well Euler's method approximated the phugoid oscillations.

To add another curve to a plot, simply type a second `plt.plot()` statement.  We also added a legend using the format

```Python
plt.legend(["line 1","line 2"])
```

```
plt.figure(figsize=(10,4))
plt.ylim(40,150)             #y-axis plot limits
plt.tick_params(axis='both', labelsize=14) #increase font size for ticks
plt.xlabel('t', fontsize=14) #x label
plt.ylabel('z', fontsize=14) #y label
plt.plot(t,z)
plt.plot(t, z_exact)
plt.legend(['Numerical Solution','Analytical Solution']);
```

That looks like pretty good agreement, but what's happening towards the end?

Euler's method, like all numerical methods, is an approximation.  If the method is *convergent*, the approximation will get closer and closer to the exact solution as we reduce the size of the step, $\Delta t$. If the method is convergent, the error in our numerical method will tend to zero.

## Convergence

To compare the two solutions, we need to use a **norm** of the difference, like the $L_1$ norm, for example.

$$ E = \Delta t \sum_{n=0}^N \left|z(t_n) - z_n\right|$$

The $L_1$ norm is the sum of the individual differences between the exact and the numerical solutions, at each mesh point.

We check for convergence by calculating the numerical solution using progressively smaller values of `dt`. We already have most of the code that we need.  We just need to add an extra loop and an array of different $\Delta t$ values to iterate through.

You should read the documentation for the function [enumerate](https://docs.python.org/2/library/functions.html#enumerate) and make sure you understand how we're using it below.

#####Warning:

The cell below can take a little while to finish (the last $\Delta t$ value alone requires 1 million iterations!).  If the cell is still running, the input label will say `In [*]`.  When it finishes, the `*` will be replaced by a number.

```
# time-increment array
dt_values = numpy.array([0.1, 0.05, 0.01, 0.005, 0.001, 0.0001])

# array that will contain solution of each grid
z_values = numpy.empty_like(dt_values, dtype=numpy.ndarray)

for i, dt in enumerate(dt_values):
    N = int(T/dt)    # number of time-steps
    ### discretize the time using numpy.linspace() ###
    t = numpy.linspace(0.0, T, N)

    # initial conditions
    u = numpy.array([z0, v])
    z = numpy.empty_like(t)
    z[0] = z0

    # time loop - Euler method
    for n in range(N-1):
        ### compute next solution using Euler method ###
        u[:] = u[:] + dt*numpy.array([u[1], g*(1-u[0]/zt)])
        z[n] = u[0]   # store the elevation at time-step n+1

    z_values[i] = z.copy()    # store the total elevation calculation grid i
```

### Calculate the error

We now have numerical solutions for each $\Delta t$ in the array `z_values`.  To calculate the error corresponding to each $\Delta t$, we can write a function!

```
def get_error(z, dt):
    """Returns the error relative to analytical solution using L-1 norm.

    Parameters
    ----------
    z : list of float
        numerical solution.
    z_exact : list of float
        analytical solution.
    dt : float
        time increment.

    Returns
    -------
    err : float
        L_{1} norm of the error with respect to the exact solution.
    """
    N = len(z)
    t = numpy.linspace(0.0, T, N)

    z_exact = v*(zt/g)**.5*numpy.sin((g/zt)**.5*t)+\
                (z0-zt)*numpy.cos((g/zt)**.5*t)+zt

    return dt * numpy.sum(numpy.abs(z-z_exact))
```

**Note**: in the last line of the function below, we perform an 'array operation':

`z - z_exact`

We are *not* subtracting one value from another.  Instead, we are taking the difference between elements at each corresponding index in both arrays.  Here is a quick example:

```
a = numpy.array([1, 2, 3])
b = numpy.array([4, 4, 4])

b - a
```

    array([3, 2, 1])

Now, we iterate through each $\Delta t$ value and calculate the corresponding error.

```
error_values = numpy.empty_like(dt_values)

for i, dt in enumerate(dt_values):
    ### call the function get_error() ###
    error_values[i] = get_error(z, dt)
```

Remember, *if* the method is convergent then the error should get smaller as  $\Delta t$ gets smaller.  To visualize this, let's plot $\Delta t$ vs. error.  If you use `plt.plot` you won't get a very useful result.  Instead, use `plt.loglog` to create the same plot with a log-log scale.  This is what we do almost always to assess the errors of a numerical scheme graphically.

```
plt.figure(figsize=(10, 6))
plt.tick_params(axis='both', labelsize=14) #increase tick font size
plt.grid(True)                         #turn on grid lines
plt.xlabel('$\Delta t$', fontsize=16)  #x label
plt.ylabel('Error', fontsize=16)       #y label
plt.loglog(dt_values, error_values, 'ko-')  #log-log plot
plt.axis('equal')                      #make axes scale equally;
```

This is the kind of result we like to see!  As $\Delta t$ shrinks, the error gets smaller and smaller, like it should.

##### Challenge!

We calculated the error for several different timestep sizes using two nested `for` loops.  That worked, but whenever possible, we like to re-use code (and not just copy and paste it!).

Create a function that implements Euler's method and re-write the error analysis cell to use your function.

---

######The cell below loads the style of this notebook.

```
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
```

<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>
<style>

@font-face {
    font-family: "Computer Modern";
    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

#notebook_panel { /* main background */
    background: rgb(245,245,245);
}

div.cell { /* set cell width */
    width: 750px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    width: 1000px;
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running {
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    background-color: rgb(256,256,256);
    border-radius: 0px;
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Alegreya Sans' sans-serif;
    line-height: 140%;
    font-size: 125%;
    font-weight: 400;
    width:600px;
    margin-left:auto;
    margin-right:auto;
}

/* Formatting for header cells */
.text_cell_render h1 {
    font-family: 'Nixie One', serif;
    font-style:regular;
    font-weight: 400;
    font-size: 45pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}
.text_cell_render h2 {
    font-family: 'Nixie One', serif;
    font-weight: 400;
    font-size: 30pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}

.text_cell_render h3 {
    font-family: 'Nixie One', serif;
    margin-top:16px;
	font-size: 22pt;
    font-weight: 600;
    margin-bottom: 3px;
    font-style: regular;
    color: rgb(102,102,0);
}

.text_cell_render h4 {    /*Use this for captions*/
    font-family: 'Nixie One', serif;
    font-size: 14pt;
    text-align: center;
    margin-top: 0em;
    margin-bottom: 2em;
    font-style: regular;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-family: 'Nixie One', sans-serif;
    font-weight: 400;
    font-size: 16pt;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'PT Mono', sans-serif;
    font-weight: 300;
    font-size: 9pt;
    line-height: 100%;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
        font-family: "PT Mono";
        font-size: 90%;
}

</style>

max_stars_repo_name: mojutijatir/numerical-mooc

max_issues_repo_name: mojutijatir/numerical-mooc

max_forks_repo_name: mojutijatir/numerical-mooc


text:
#
Fracciones y divisibilidad.

Para realizar la mayor parte de los cálculos matemáticos de modo exacto utilizaremos la biblioteca **SYMPY**, que no viene de serie con Python, aunque si viene incluida en las distribuciones de Python para cálculo científico. Para ello se debe importar la biblioteca, y en general lo haremos con la sentencia **from sympy import ***.

Además presentaremos los resultados en formato $\LaTeX$. Para ello utilizamos la función **init_printing** (escribimos "init" y pulsamos el tabulador). Estas dos expresiones deben aparecer siempre al comienzo del notebook.

Para escribir fracciones utilizaremos la función **Rational** (es imprescindible escribir la primera letra con mayúsculas). Con las fracciones se pueden realizar todas las operaciones habituales y los resultados se presentan también en forma de fracción y no en formato decimal. Pero si en la operación introducimos un número decimal, la operación la realiza ya en decimales.

###
Realiza las siguientes operaciones con fracciones.

* Simplifica la fracción $\displaystyle \frac{46}{50}$

* $\displaystyle \frac{2}{3} +\frac{6}{7} + 1$

* $\displaystyle \frac{2}{5} \times \frac{7}{6}$

* $\displaystyle \frac{3}{4} + 1.0$

* Calcula $\displaystyle \left(\frac{2}{3}\right)^{100}$ y da el resultado en decimales.

Sympy puede factorizar números enteros de tamaño razonable con la función **factorint**. El resultado es lo que en Python se denomina un diccionario, donde aparecen el factor primo y el exponente separador por dos puntos. Si queremos que aparezca en el formato estandar, debemos añadir **visual = True**.

###
Factoriza los números 612 y 5292 y comprueba el resultado.

Las funciones para calcular el máximo común denominador y el mínimo común múltiplo son **gcd** y **lcm**, que son las iniciales inglesas que corresponden a nuestros **mcd** y **mcm**.

###
Calcula el mcd y el mcm de los números anteriores.

Para comprobar si un número es primo se utiliza la función **isprime**. Para números "pequeños" utiliza algoritmos exactos para su comprobación, pero para números grandes utiliza métodos probabilísticos. Parece que  esta función no es necesaria, pues con la factorización se puede saber si un número es primo o no lo es. Sin embargo los métodos de factorización son muy costosos computacionalmente y los test de primalidad no lo son tanto.

###
Comprueba si $5673$ es primo. Lo mismo con $2^{32}+ 1$. Factorízalos.

Para trabajar con divisores Sympy tiene dos funciones: **divisors**, que devuelve una lista con todos los divisores, y **divisor_count**, que únicamente cuenta el número de divisores.

###
Calcula todos los divisores de 180.

La función **sum** suma todos los elementos de una lista.

###
Demuestra que 496 es un número perfecto.

Para calcular el primo anterior y el primo siguiente se utilizan las funciones **nextprime** y **prevprime**. Para calcular el número de primos menores que una cantidad dada se utiliza **primepi**.

###
Calcula el primo anterior y posterior a 1000.

###
Calcula $\pi(1000)$ (la cantidad de números primos menores que 1000).

Para calcular con decimales cualquier expresión se utiliza el método **evalf**. Para utilizar un método, primero debemos escribir la expresión, después un punto y finalmente el nombre del método. Muchas veces lo más cómodo es guardar la expresión en una variable y emplear el método sobre dicha variable. En este caso el método **evalf** tiene un argumento opcional, que es el número de cifras significativas.

Una de las "funciones mágicas" de Ipython es **%%time**, que nos permite calcular el tiempo de cómputo de la celda en cuestión. Esto es una función propia de IPython, no disponible en Python.

###
Calcula el número $\pi$ y $\sqrt{2}$ con 50 cifras significativas. Calcula algunos tiempos de cómputo.

ext: ipynb

max_stars_repo_path: notebooks/03.- Fracciones y divisibilidad.ipynb

max_stars_repo_name: disoftw/python-for-maths

max_issues_repo_name: disoftw/python-for-maths

max_forks_repo_name: disoftw/python-for-maths


text:
# 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>


text:
# Section 1.2 $\quad$ Matrices

## Definitions

- A rectangular array of $m\times n$ real or complex numbers arranged in $m$ horizontal **rows** and $n$ vertical **columns**:
$$
    A =
    \left[
      \begin{array}{cccc}
          a_{11} & a_{12} & \cdots & a_{1n} \\
          a_{21} & a_{22} & \cdots & a_{2n} \\
          \vdots & \vdots & \ddots & \vdots \\
          a_{m1} & a_{m2} & \cdots & a_{mn} \\
      \end{array}
    \right]\qquad\qquad (1)
$$
    is called an $\underline{\hspace{1.5in}}$. We say that the matrix $A$ is $\underline{\hspace{1in}}$.

- The $i$**th row** and the $j$**th column** of $A$ are

- If $m=n$, we say that $A$ is a $\underline{\hspace{2in}}$, and that the numbers $a_{11}$, $a_{22}$, $\cdots$, $a_{nn}$ form the $\underline{\hspace{1in}}$ of $A$.

- The number $a_{ij}$ is called the $\underline{\hspace{1.5in}}$ of $A$, or the $\underline{\hspace{1.5in}}$ of $A$. We can write the matrix $A$ as

- An $n\times 1$ matrix is called an $\underline{\hspace{1.5in}}$ or a $\underline{\hspace{1.5in}}$ when $n$ is understood.

- Two $m\times n$ matrices $A = [a_{ij}]$ and $B = [b_{ij}]$ are **equal** if

### Example 1

Examples of matrices
$$
  A = \left[
     \begin{array}{ccc}
       1 & 2 & 3 \\
     -1 & 0 & 1 \\
     \end{array}
   \right],\quad
  B = \left[
     \begin{array}{cc}
       1+i & 4i \\
       2-3i & -3 \\
     \end{array}
   \right],\quad
  C = \left[
     \begin{array}{c}
       1\\
       -1\\
       2 \\
     \end{array}
   \right].
$$

### Example 2

If
$$
  \left[
    \begin{array}{cc}
      a+2b & 2a-b \\
      2c+d & c-2d \\
    \end{array}
  \right] =
  \left[
    \begin{array}{cc}
      4 & -2 \\
      4 & -3 \\
    \end{array}
  \right]
$$
find $a,b,c$, and $d$.

```python
from sympy import *

a, b, c, d = symbols('a b c d');

eq1 = a + 2*b - 4;
eq2 = 2*a - b + 2;
eq3 = 2*c + d - 4;
eq4 = c - 2*d + 3;

solve([eq1, eq2, eq3, eq4], (a, b, c, d))
```

    {a: 0, d: 2, c: 1, b: 2}

## Matrix Operations

**Addition**: If $A = [a_{ij}]$ and $B = [b_{ij}]$ are both $m\times n$ matrices, then the **sum** $A+B$ is an $m\times n$ matrix $C = [c_{ij}]$ defined by

### Example 3

Let
$$
  A =
  \left[
    \begin{array}{ccc}
      1 & -2 & 3 \\
      2 & -1 & 4 \\
    \end{array}
  \right], ~~~~
  B =
  \left[
    \begin{array}{ccc}
      0 & 2 & 1 \\
      1 & 3 & -4 \\
    \end{array}
  \right].
$$
Then
$$
  A + B = \hspace{4in}
$$

```python
from sympy import *

A = Matrix([[1, -2, 3], [2, -1, 4]]);
B = Matrix([[0, 2, 1], [1, 3, -4]]);

A + B
```

    Matrix([
    [1, 0, 4],
    [3, 2, 0]])

**Remark:** the sum of $A$ and $B$ is defined only when $\underline{\hspace{3in}}$.

## Scalar Multiplication

If $A = [a_{ij}]$ is an $m\times n$ matrix and $r$ is a real number, then the **scalar multiple** of $A$ by $r$ is an $m\times n$ matrix $C = [c_{ij}]$ where

### Example 4

Let
$$
  A =
  \left[
    \begin{array}{ccc}
      1 & -2 & 3 \\
      2 & -1 & 4 \\
    \end{array}
  \right], ~~~~
  B =
  \left[
    \begin{array}{ccc}
      0 & 2 & 1 \\
      1 & 3 & -4 \\
    \end{array}
  \right].
$$
Then
$$
  5A = \hspace{4in}
$$
$$
  A-B = \hspace{4in}
$$

```python
from sympy import *

A = Matrix([[1, -2, 3], [2, -1, 4]]);
B = Matrix([[0, 2, 1], [1, 3, -4]]);

5*A
```

    Matrix([
    [ 5, -10, 15],
    [10,  -5, 20]])

```python
from sympy import *

A = Matrix([[1, -2, 3], [2, -1, 4]]);
B = Matrix([[0, 2, 1], [1, 3, -4]]);

A - B
```

    Matrix([
    [1, -4, 2],
    [1, -4, 8]])

## More Definitions

- If $A_1,A_2,\cdots,A_k$ are matrices of the same dimension, and $c_1,c_2,\cdots,c_k$ are real numbers, then
$$
  c_1A_1 + c_2A_2 + \cdots + c_kA_k \qquad\qquad (2)
$$
    is called a $\underline{\hspace{2.5in}}$ of $A_1,A_2,\cdots,A_k$.

    The numbers $c_1,c_2,\cdots,c_k$ are called $\underline{\hspace{1.5in}}$.

- The linear combination (2) can be expressed using summation notation:

- If $A = [a_{ij}]$ is an $m\times n$ matrix, then the $\underline{\hspace{1.5in}}$ of $A$ is an $n\times m$ matrix $A^T = [a_{ij}^T]$ defined by
$$
        a_{ij}^T = \qquad\qquad\qquad\qquad
$$

### Example 5

Let
$$
  A =
  \left[
    \begin{array}{ccc}
      1 & -2 & 3 \\
      2 & -1 & 4 \\
    \end{array}
  \right]
$$
Then
$$
  A^T = \hspace{4in}
$$

```python
from sympy import *

A = Matrix([[1, -2, 3], [2, -1, 4]]);

A.T
```

    Matrix([
    [ 1,  2],
    [-2, -1],
    [ 3,  4]])

### Example 6

Suppose $A$ and $B$ are both $m\times n$ matrices. Then
$$
  (2A + 3B)^T = \hspace{4in}
$$

$$
  (A^T)^T = \hspace{4in}
$$

$$
  (A - 2B^T)^T = \hspace{4in}
$$