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

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).
  2. Distribución óptima de recursos.
    • Distribución de rutas de vuelo.
    • Encontrar una dieta óptima.
  3. Optimización financiera
    • Administración de riesgos.

En esta clase veremos aspectos básicos de optimización. En específico, veremos cómo obtener máximos y mínimos de una función escalar de una variable (como en cálculo diferencial).

Basamos todos los resultados en los siguientes teoremas:

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$ tiene un mínimo global en $x=0$, pues

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

In [37]:
# Librería de cálculo simbólico
import sympy as sym
# Para imprimir en formato TeX
from sympy import init_printing; init_printing(use_latex='mathjax')

In [38]:
sym.var('x', real = True)
f = x**2
f


Out[38]:
$$x^{2}$$

In [39]:
df = sym.diff(f, x)
df


Out[39]:
$$2 x$$

In [40]:
x_c = sym.solve(df, x)
x_c[0]


Out[40]:
$$0$$

Veamos la gráfica...


In [41]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [42]:
f_num = sym.lambdify([x], f, 'numpy')
x_vec = np.linspace(-5, 5, 100)

plt.plot(x_vec, f_num(x_vec))
plt.xlabel('$x$')
plt.ylabel('$x^2$')
plt.show()


Otra manera de hacer lo anterior


In [44]:
def f(x):
    return x**2

In [45]:
f_sym = f(x)
f_sym


Out[45]:
$$x^{2}$$

In [47]:
df = sym.diff(f(x), x)
df


Out[47]:
$$2 x$$

In [48]:
x_c = sym.solve(df, x)
x_c[0]


Out[48]:
$$0$$

El converso del teorema anterior no es cierto.

Actividad

Considere $g(x)=x^3$.

  • Usando sympy, muestre que $g'(0)=0$.
  • Sin embargo, descartar que $x=0$ es un extremo de $g(x)$ viendo su gráfica.

In [ ]:


In [ ]:


In [ ]:

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$ tiene un mínimo relativo en $x=0$.

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


In [7]:
f = x**2
#d2f = sym.diff(f, x, x)
d2f = sym.diff(f, x, 2)
d2f


Out[7]:
$$2$$

In [8]:
d2f>0


Out[8]:
$$\mathrm{True}$$

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

Actividad

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


In [ ]:

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

  • Evaluar $f$ en los extremos $x=a$ y $x=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.
  • 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-6x$ en $\left[0,5\right]$.

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


In [11]:
f = x**2-6*x
f


Out[11]:
$$x^{2} - 6 x$$

In [17]:
df = sym.diff(f, x)
df


Out[17]:
$$2 x - 6$$

In [18]:
x_c = sym.solve(df, x)
x_c


Out[18]:
$$\left [ 3\right ]$$

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


In [21]:
f.subs(x, 0), f.subs(x, 5), f.subs(x, x_c[0])


Out[21]:
$$\left ( 0, \quad -5, \quad -9\right )$$

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


In [29]:
f_num = sym.lambdify([x], f, 'numpy')
x_vec = np.linspace(0, 5, 100)

plt.figure(figsize=(8,6))
plt.plot(x_vec, f_num(x_vec), 'k', label = '$y=f(x)$')
plt.plot([0], [0], '*r', label = '$(0,0=\max_{0\leq x\leq 5} f(x))$')
plt.plot([3], [-9], '*b', label = '$(3,-9=\min_{0\leq x\leq 5} f(x))$')
plt.legend(loc='best')
plt.xlabel('x')
plt.show()


Actividad

Determinar los valores extremos absolutos de $h(x)=x^3-3x$ en $\left[-2.2,1.8\right]$, usando sympy. Mostrar en una gráfica.


In [ ]:


In [ ]:


In [ ]:

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.


In [49]:
sym.var('x y')
x, y


Out[49]:
$$\left ( x, \quad y\right )$$

In [50]:
def f(x, y):
    return x**2 + y**2

In [52]:
dfx = sym.diff(f(x,y), x)
dfy = sym.diff(f(x,y), y)
dfx, dfy


Out[52]:
$$\left ( 2 x, \quad 2 y\right )$$

In [62]:
xy_c = sym.solve([dfx, dfy], [x, y])
xy_c


Out[62]:
$$\left \{ x : 0, \quad y : 0\right \}$$

In [61]:
x_c, y_c = xy_c[x], xy_c[y]

In [76]:
d2fx = sym.diff(f(x,y), x, 2)
d2fy = sym.diff(f(x,y), y, 2)
dfxy = sym.diff(f(x,y), x, y)

Jf = sym.Matrix([[d2fx, dfxy], [dfxy, d2fy]])
Jf.eigenvals()


Out[76]:
$$\left \{ 2 : 2\right \}$$

In [77]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

In [79]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x = np.linspace(-2, 2, 100)
y = x
X, Y = np.meshgrid(x, y)

ax.plot_surface(X, Y, f(X, Y))
ax.plot([0], [0], [0], '*r')


Out[79]:
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7fd179f4cbe0>]

Tarea (para el martes 31 de Octubre a las 23:00).

Elaborar un algoritmo, usando sympy, que devuelva el máximo y el mínimo absoluto de una función en un intervalo dado (finito), y que además dibuje una gráfica de la función en dicho intervalo y señalando los puntos máximo y mínimo absolutos.

Recordar: examen el viernes 27 de octubre. Lo entregan el miércoles 1 de noviembre antes de las 23:00.

Para este viernes 27 de octubre: subir los avances que tengan del proyecto.

Created with Jupyter by Esteban Jiménez Rodríguez.