notebook.community

Edit and run

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 [1]:



    


# 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 [8]:



    


sym.var('x', real = True)
f = x**2
f



    


















    
Out[8]:







$$x^{2}$$
















In [12]:



    


df = sym.diff(f, x)
df



    


















    
Out[12]:







$$2 x$$
















In [14]:



    


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



    


















    
Out[14]:







$$0$$

Veamos la gráfica...





In [16]:



    


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



    











In [18]:



    


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 [21]:



    


def f(x):
    return 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 [24]:



    


def g(x):
    return x**3



    











In [26]:



    


dg = sym.diff(g(x), x)
x_c = sym.solve(dg, x)
x_c



    


















    
Out[26]:







$$\left [ 0\right ]$$
















In [27]:



    


x_vec = np.linspace(-2,2,100)
plt.figure(figsize=(8,6))
plt.plot(x_vec, g(x_vec))
plt.xlabel('$x$')
plt.ylabel('$x^3$')
plt.show()

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 [29]:



    


f = x**2
#d2f = sym.diff(f, x, x)
d2f = sym.diff(f, x, 2)
d2f



    


















    
Out[29]:







$$2$$
















In [30]:



    


d2f>0



    


















    
Out[30]:







$$\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 [31]:



    


g(x)



    


















    
Out[31]:







$$x^{3}$$
















In [32]:



    


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



    


















    
Out[32]:







$$6 x$$
















In [33]:



    


d2g.subs(x, 0)



    


















    
Out[33]:







$$0$$

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 [34]:



    


f = x**2-6*x
f



    


















    
Out[34]:







$$x^{2} - 6 x$$
















In [35]:



    


df = sym.diff(f, x)
df



    


















    
Out[35]:







$$2 x - 6$$
















In [36]:



    


x_c = sym.solve(df, x)
x_c



    


















    
Out[36]:







$$\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 [37]:



    


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 [38]:



    


def h(x):
    return x**3-3*x



    











In [39]:



    


dh = sym.diff(h(x), x)
x_c = sym.solve(dh, x)
x_c



    


















    
Out[39]:







$$\left [ -1, \quad 1\right ]$$
















In [40]:



    


h(-2.2), h(x_c[0]), h(x_c[1]), h(1.8)



    


















    
Out[40]:







$$\left ( -4.048000000000003, \quad 2, \quad -2, \quad 0.4320000000000004\right )$$
















In [41]:



    


x_vec = np.linspace(-2.2, 1.8, 100)

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

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 [42]:



    


sym.var('x y')
x, y



    


















    
Out[42]:







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
















In [43]:



    


def f(x, y):
    return x**2 + y**2



    











In [44]:



    


dfx = sym.diff(f(x,y), x)
dfy = sym.diff(f(x,y), y)
dfx, dfy



    


















    
Out[44]:







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
















In [45]:



    


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



    


















    
Out[45]:







$$\left \{ x : 0, \quad y : 0\right \}$$
















In [46]:



    


x_c, y_c = xy_c[x], xy_c[y]
x_c



    


















    
Out[46]:







$$0$$
















In [49]:



    


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[49]:







$$\left \{ 2 : 2\right \}$$
















In [50]:



    


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D



    











In [52]:



    


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([x_c], [y_c], [f(x_c,y_c)], '*r')



    


















    
Out[52]:








[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7fcb03c08470>]

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.

Content source: cesarcontre/Simulacion2017

Similar notebooks:

notebook.community | gallery | about