# Optimización de funciones escalares diferenciables con `SymPy`



- Mediante optimización se obtienen soluciones elegantes tanto en teoría como en ciertas aplicaciones. 
- La teoría de optimización usa elementos comenzando con cálculo elemental y álgebra lineal básica, y luego se extiende con análisis funcional y convexo.
- Las aplicaciones en optimización involucran ciencia, ingeniería, economía, finanzas e industria.
- El amplio y creciente uso de la optimización lo hace escencial para estudiantes y profesionales de cualquier rama de la ciencia y la tecnología.

**Referencia:**
- http://www.math.uwaterloo.ca/~hwolkowi//henry/reports/talks.d/t06talks.d/06msribirs.d/optimportance.shtml





### Algunas aplicaciones son:

1. **Ingeniería:**
 - Encontrar la composición de equilibrio de una mezcla de diferentes átomos.
 - Planeación de ruta para un robot (o vehículo aéreo no tripulado).
 - Planeación de la mano de obra óptima en una construcción o planta de producción.
 - Análisis conformacional - Isomería conformacional en modelado molecular
 - Desempeño en fermentación de productos.
 - Estimación de parámetros en bioquímica e ingeniería metabólica.
 
 
2. **Distribución óptima de recursos:**
 - Distribución de rutas de vuelo.
 - Encontrar una dieta óptima.
 - Planeación de ruta óptima.
 
 
3. **Optimización financiera:**
 - Administración de riesgos.
 - Portafolios de inversión.
 - Problema de maximización de la utilidad.
 - 

### <p style="color:green;">En esta clase veremos aspectos básicos de optimización. Específicamente, veremos cómo obtener máximos y mínimos de una función escalar de una variable (como en cálculo diferencial).</p>



___
## 0. Librerías que usaremos


Como ya dijimos en la primer clase `python` es el lenguaje de programación (el cual es de alto nivel). Sin embargo, `python` solo tiene unos pocos comandos primitivos y para hacer más fácil su uso en nuestras actividades de simulación en ingeniería, otras personas ya han escrito ciertas librerías por nosotros.

### 0.1 `NumPy`



`NumPy` (Numerical Python) es la librería fundamental para computación científica (numérica) con `Python`. Contiene, entre otras cosas:
- un objeto tipo arreglo N-dimensional muy poderoso
- funciones sofisticadas
- funciones de álgebra lineal, transformada de Fourier y números aleatorios.

Por lo anterior, `NumPy` es de amplio uso entre la comunidad científica e ingenieril (por su manejo de cantidades vectoriales). De la misma manera, se usa para guardar datos. Para nuestros propósitos, se puede usar libremente.

**Referencia:**
- http://www.numpy.org/

`NumPy` ya viene incluido en la instalación estándar de Anaconda por defecto. Para comenzar a usarlo, solo debemos de importarlo:


```python
# Importar librería numpy
import numpy as np
```

### 0.2 `SymPy`



`SymPy` (Symbolic Python) es una librería de `Python` para matemáticas simbólicas. Su objetivo es convertirse en un sistema de álgebra computacional con las mejores características, manteniendo el código lo más simple posible para que sea comprensible.

**Referencia:**
- http://www.sympy.org/en/index.html

`SymPy` ya viene incluido en la instalación estándar de Anaconda por defecto. Para comenzar a usarlo, solo debemos de importarlo:


```python
#Importarmos la librería sympy
import sympy as sym

# Imprimir datos en formato Latex
sym.init_printing(use_latex = 'mathjax')
```

La funcionalidad de imprimir en formato LaTex que nos da `SymPy` mediante el proyecto `mathjax` hace de `SymPy` una herramienta muy atractiva...

Notar que en `SymPy` y en `NumPy` existen funciones con el mismo nombre, pero reciben tipos de datos diferentes...


```python
sym.var('x')
```




$$x$$




```python
sym.cos(x), np.cos(3.1416)
```




$$\left ( \cos{\left (x \right )}, \quad -0.9999999999730151\right )$$



Explicar el uso de la sintaxis `from numpy import *` y sus peligros (no recomendable).


```python
# importar con * y ver que pasa
#from numpy import *
#from sympy import *
# No recomendado
#b=array([3])
#b
#cos(x), 
#cos(3.1416)
```

### 0.3 `PyPlot` de `matplotlib`



El módulo `PyPlot` de la librería `matplotlib` contiene funciones que nos permite generar una gran cantidad de gráficas rápidamente. Las funciones de este módulo están escritas con el mismo nombre que las funciones para graficar en `Matlab`.

**Referencia:**
- https://matplotlib.org/api/pyplot_summary.html


```python
# importar matplotlib.pyplot
import matplotlib.pyplot as plt
# comando para que las gráficas salgan en la misma ventana
%matplotlib inline
```

Ya que revisamos todas las librerías que usaremos, empecemos con la clase como tal...

___
<p style="color:red;">Todos los resultados que veremos son basados en el siguiente teorema:</p>

## 1. Teorema de Fermat (análisis)

Si una función $f(x)$ alcanza un máximo o mínimo local en $x=c$, y si la derivada $f'(c)$ existe en el punto $c$, entonces $f'(c) = 0$.

### Ejemplo

Sabemos que la función $f(x)=(x-2)^2$ tiene un mínimo global en $x=2$, pues

$$f(x)=(x-2)^2\geq0,\qquad\text{y}\qquad f(x)=(x-2)^2=0 \qquad\text{si y solo si}\qquad x=2.$$


```python
#Declarar variable x como real
sym.var('x', real=True)
```




$$x$$




```python
# Declarar función f y mostrar
f = (x-2)**2
f
```




$$\left(x - 2\right)^{2}$$




```python
# Encontrar la primera derivada de la función f
df = sym.diff(f,x)
df
```




$$2 x - 4$$




```python
#Resolver f'(x)= 0 y mostrar punto
xc = sym.solve(df,x)
xc
```




$$\left [ 2\right ]$$



Veamos la gráfica...


```python
# Convertir función simbólica f a una función numérica
fnum = sym.lambdify(x,f,modules='numpy')
```


```python
# Definir coordenas x
xnum = np.linspace(-4,8,200)
ynum = fnum(xnum)
```


```python
# Graficar
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum,ynum,'k',label = '$y = (x-2)^2$')
plt.plot([2],[0], '*y', ms = 15, label = '$(2,0)$')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$f(x) = (x-2)^2$')
plt.legend(loc = 'best')
plt.grid()
```

Ver diferencias entre f y f_num


```python

```

### **Otra manera de hacer lo anterior**

Revisemos el concepto de función en `Python`


```python
#Definimos una función en Python
def nuestra_funcion_f(x):
    return (x-2)**2
```


```python
x
```




$$x$$




```python
nuestra_funcion_f(x)
```




$$\left(x - 2\right)^{2}$$




```python
df = sym.diff(nuestra_funcion_f(x),x)
df
```




$$2 x - 4$$




```python
xc = sym.solve(df,x)
xc
```




$$\left [ 2\right ]$$




```python
# Graficar
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum,nuestra_funcion_f(xnum),'k',label = '$y = (x-2)^2$')
plt.plot([2],[0], '*y', ms = 15, label = '$(2,0)$')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$f(x) = (x-2)^2$')
plt.legend(loc = 'best')
plt.grid()
```

**El converso del teorema anterior no es cierto.**

### Actividad
Considere $g(x)=(x-2)^3$. 
- Usando `sympy`, muestre que $g'(2)=0$.
- Sin embargo, descartar que $x=2$ es un extremo de $g(x)$ viendo su **gráfica**.


```python
def g(x):
    return (x-2)**3
```


```python
x
```




$$x$$




```python
df = sym.diff(g(x),x)
df
```




$$3 \left(x - 2\right)^{2}$$




```python
xc = sym.solve(df,x)
xc
```




$$\left [ 2\right ]$$




```python
# Graficar
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum,g(xnum),'k',label = '$y = (x-2)^3$')
plt.plot([2],[0], '*y', ms = 15, label = '$(2,0)$')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$g(x) = (x-2)^3$')
plt.legend(loc = 'best')
plt.grid()
```

## 2. Criterio de la segunda derivada

Sea $f(x)$ una función tal que $f’(c)=0$ y cuya segunda derivada existe en un intervalo abierto que contiene a $c$.
- Si $f’’(c)>0$, entonces $f(c)$ es un mínimo relativo.
- Si $f’’(c)<0$, entonces $f(c)$ es un máximo relativo.
- Si $f’’(c)=0$, entonces el criterio no decide. 

### Ejemplo

Mostrar, usando `sympy`, que la función $f(x)=(x-2)^2$ tiene un mínimo relativo en $x=2$.

Ya vimos que $f'(2)=0$. Notemos que:


```python
def f(x):
    return (x-2)**2
```


```python
d2f = sym.diff(f(x),x,2)
d2f
```




$$2$$




```python
d2f>0
```




$$\mathrm{True}$$



Por tanto, por el criterio de la segunda derivada, $f(2)=0$ es un mínimo relativo (en efecto, el mínimo global).

### Ejemplo
¿Qué pasa con $g(x)=(x-2)^3$ al intentar utilizar el criterio de la segunda derivada? (usar `sympy`).


```python
d2g = sym.diff(g(x),x,2)
d2g
```




$$6 \left(x - 2\right)$$




```python
d2g.subs(x,2)
```




$$0$$




```python
# Graficar
#Se crea nuevo eje de las x
xnum1 = np.linspace(0,4,200)
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum1,h(xnum1),'k',label = '$y = (x-2)^3$')
plt.plot([2],[0], '*y', ms = 15, label = '$(2,0)$')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$g(x) = (x-2)^3$')
plt.legend(loc = 'best')
plt.grid()
```

Como $g''(2)=0$ entonces el criterio de la segunda derivada no concluye.

### Actividad

¿Qué pasa con $h(x)=(x-2)^4$ al intentar utilizar el criterio de la segunda derivada?.


```python
def h(x):
    return (x-2)**4
```


```python
x
```




$$x$$




```python
d2h = sym.diff(h(x),x,2)
d2h
```




$$12 \left(x - 2\right)^{2}$$




```python
d2h.subs(x,2)
```




$$0$$




```python
# Graficar
#Se crea nuevo eje de las x
xnum1 = np.linspace(0,4,200)
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum1,h(xnum1),'g',label = '$y = (x-2)^3$')
plt.plot([2],[0], '*y', ms = 15, label = '$(2,0)$')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$g(x) = (x-2)^3$')
plt.legend(loc = 'best')
plt.grid()
```

## 3. Método para determinar extremos absolutos de una función continua y=f(x) en [a,b] 

- Determinar todos los valores críticos $c_1, c_2, c_3, \dots, c_n$ en $(a,b)$.
- Evaluar $f$ en todos los valores críticos y en los extremos $x=a$ y $x=b$.
- El más grande y el más pequeño de los valores de la lista $f(a), f(b), f(c_1), f(c_2), \dots, f(c_n)$ son el máximo absoluto y el mínimo absoluto, respectivamente, de $f$ en el intervalo $[a,b]$.

### Ejemplo

Determinar los extremos absolutos de $f(x)=x^2-2x+5$ en $\left[0,3\right]$.

Obtenemos los puntos críticos de $f$ en $\left[0,3\right]$:


```python
def f(x):
    return x**2-2*x+5
```


```python
df = sym.diff(f(x),x)
df
```




$$2 x - 2$$




```python
xc = sym.solve(df,x)
xc
```




$$\left [ 1\right ]$$



Evaluamos $f$ en los extremos y en los puntos críticos:


```python
f(0), f(1), f(3)
```




$$\left ( 5, \quad 4, \quad 8\right )$$



Concluimos que el máximo absoluto de $f$ en $\left[0,3\right]$ es $8$ y se alcanza en $x=3$, y que el mínimo absoluto es $4$ y se alcanza en $x=1$.


```python
xnum = np.linspace(0,3,480)
#Grafica
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum,f(xnum),'b',label = '$y = x^2-2x+5$')
plt.plot([3],[8], '*y', ms = 15, label = '$(3,8)$, máximo absoluto en el intervalo')
plt.plot([1],[4], '*y', ms = 15, label = '$(1,4)$, mínimo absoluto en el intervalo')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$f(x) = x^2-2x+5$')
plt.legend(loc = 'best')
plt.grid()

```

### Actividad
Determinar los valores extremos absolutos de  $f(x) = (x-3)^2(x+3)^2$ en el intervalo $[1,\,4]$, usando `sympy`. Mostrar en una gráfica.


```python
def f1(x):
    return (x-3)**2*(x+3)**2
```


```python
df1 = sym.diff(f1(x),x)
df1
```




$$\left(x - 3\right)^{2} \left(2 x + 6\right) + \left(x + 3\right)^{2} \left(2 x - 6\right)$$




```python
xc1 = sym.solve(df1,x)
xc1
```




$$\left [ -3, \quad 0, \quad 3\right ]$$




```python
f1(1), f1(3), f1(4)
```




$$\left ( 64, \quad 0, \quad 49\right )$$




```python
xnum = np.linspace(1,4,480)
#Grafica
# Crear nueva ventana de gráficos y se da medida a la ventana
plt.figure(figsize = (8,6))

plt.plot(xnum,f1(xnum),'b',label = '$y = (x-3)^2(x+3)^2$')
#plt.plot([3],[8], '*y', ms = 15, label = '$(3,8)$, máximo absoluto en el intervalo')
#plt.plot([1],[4], '*y', ms = 15, label = '$(1,4)$, mínimo absoluto en el intervalo')

#Nombrar los ejes de la gráfica
plt.xlabel('$x$')
plt.ylabel('$f(x) = (x-3)^2(x+3)^2$')
plt.legend(loc = 'best')
plt.grid()
```


```python

```


```python

```

### En varias variables...

El procedimiento es análogo. 

Si una función $f:\mathbb{R}^n\to\mathbb{R}$ alcanza un máximo o mínimo local en $\boldsymbol{x}=\boldsymbol{c}\in\mathbb{R}^n$, y $f$ es diferenciable en el punto $\boldsymbol{x}=\boldsymbol{c}$, entonces $\left.\frac{\partial f}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}=\boldsymbol{c}}=\boldsymbol{0}$ (todas las derivadas parciales en el punto $\boldsymbol{x}=\boldsymbol{c}$ son cero).

**Criterio de la segunda derivada:** para ver si es máximo o mínimo, se toma la segunda derivada (matriz jacobiana) y se verifica definición negativa o positiva, respectivamente.

Si se restringe a cierta región, hay ciertas técnicas. La más general, pero también la más compleja es la de **multiplicadores de Lagrange**.

**Ejemplo:** Determinar los valores extremo de $z(x,y)= (x-3)^2 + (y+5)^2$ usando el criterio de la segunda derivada y corroborarlo con la gráfica.


```python
sym.var('x y')
```




$$\left ( x, \quad y\right )$$




```python
def z(x,y):
    return (x-3)**2 + (y+5)**2
```


```python
dzx = sym.diff(z(x,y),x) 
dzy = sym.diff(z(x,y),y)
dzx,dzy
```




$$\left ( 2 x - 6, \quad 2 y + 10\right )$$




```python
xy_c = sym.solve([dzx,dzy],[x,y])
xy_c
```




$$\left \{ x : 3, \quad y : -5\right \}$$




```python
d2zx = sym.diff(z(x,y),x,2) 
d2zy = sym.diff(z(x,y),y,2)
dzxy = sym.diff(z(x,y),x,y)

Jz = sym.Matrix([[d2zx,dzxy],[dzxy,d2zy]])
Jz
```




$$\left[\begin{matrix}2 & 0\\0 & 2\end{matrix}\right]$$




```python
Jz.eigenvals()
Jz.det()
```




$$4$$



### Tarea.

1. Obtener máximos y mínimos absolutos, usando `sympy`, de las funciones dadas en los intervalos dados, y graficar la función en dicho intervalo señalando los puntos máximo y mínimo.

    + $f(x)=-3x^3+9x+2$ en $[-1.8,\,1.5]$
    + $g(x) = x^3-6x^2+9x-4$ en el intervalo $[0, \,2.5]$



<footer id="attribution" style="float:right; color:#808080; background:#fff;">
Created with Jupyter by Cristian Camilo Zapata Zuluaga.
</footer>

ee307105c28736c78ce809425437caeea01b997f
ipynb
Jupyter Notebook
Modulo1/Clase4_OptimizacionSympy.ipynb
ArellanoMCarlos/SimMat2019-1
c84b92a581916572352615806d31961a468da3d9
MIT
Modulo1/Clase4_OptimizacionSympy.ipynb
ArellanoMCarlos/SimMat2019-1
c84b92a581916572352615806d31961a468da3d9
MIT
Modulo1/Clase4_OptimizacionSympy.ipynb
ArellanoMCarlos/SimMat2019-1
c84b92a581916572352615806d31961a468da3d9
MIT
Qwen/Qwen-72B
1. YES
2. YES