ER3.4. Método indirecto de Lyapunov aplicado al modelo no lineal de un péndulo simple

Introducción (planteamiento)

Dado el siguiente modelo no lineal de un péndulo simple con un hilo rígido extremadamente delgado, analizar la estabilidad en cada uno de los puntos de equilibrio:

$\left\{ \begin{aligned} &\dot{x}_1=x_2\\ &\dot{x}_2=-\left( \frac{g}{l} \right) \sin x_1-\left( \frac{f}{m} \right) x_2 + \frac{u}{ml}\\\end{aligned} \right.\\ u \ne 0$

Método (plan de solución)

Pasos: (i) cálculo de los puntos de equilibrio, (ii) linealización, (iii) cálculo de los valores propios en cada punto de equilibrio, (iv) análisis de los resultados.

Resultados (solución)

Puntos de equilibrio:

$x_{20}=0\\ -\left( \frac{g}{l} \right) \sin x_{10}-\left( \frac{f}{m} \right) x_{20}+\frac{1}{ml}u_0=0$

Es decir,

$x_{20}=0\\ \sin x_{10}=\frac{u_0}{mg}$    

Existen puntos de equilibrio siempre y cuando $\left| u_0/(mg) \right|\leqslant 1$. 

Por ejemplo, si la entrada es la mitad del peso $(u_0=0.5mg)$, entonces $\sin x_0=0.5$, $x_0=\pi/6=0.5236~\mathrm{rad}$ o  $x_0=5\pi /6=2.6180~\mathrm{rad}$. Tomando valores arbitrarios de los parámetros, la entrada igual a la mitad del peso y las condiciones iniciales iguales a cero y cerca de $2.6180~\mathrm{rad}$, se obtiene siempre el punto de equilibrio en 0.5236 rad, como se observa de las siguientes respuestas temporales obtenidas a partir de la simulación (se omite el código).

Respuesta temporal partiendo de un ángulo igual a cero:

Partiendo de un ángulo igual a 2.6180 rad:

El retrato de fase con la función pplane de MATLAB ilustra mejor los casos anteriores, donde se debe ubicar en el gráfico la condición inicial y seguir la respectiva trayectoria de estado con su dirección (se invita al lector a observar lo que sucede cuando $x_1(0)$ está por debajo o por encima de 2.618 rad:

Para la estabilidad se calcula el jacobiano:

$\left. \mathbf{A}=\frac{\partial f}{\partial \mathbf{x}} \right|_{\mathbf{x}_0}=\left[ \begin{matrix} \frac{\partial f_1}{\partial x_1}&  \frac{\partial f_1}{\partial x_2}\\ \frac{\partial f_2}{\partial x_1}&  \frac{\partial f_2}{\partial x_2}\\\end{matrix} \right] _{\mathbf{x}_0}=\left[ \begin{matrix} 0&  1\\ -\left( \frac{g}{l} \right) \cos x_1&  -\frac{f}{m}\\\end{matrix} \right] _{\mathbf{x}_0}$

Para $\mathbf{x}_0=(\pi /6,0)$ se tiene:

$\mathbf{A}=\left[ \begin{matrix} 0& 1\\ -\frac{\sqrt{3}g}{2l}&  -\frac{f}{m}\\\end{matrix} \right]$

Cálculo del jacobiano con la función jacobian de MATLAB:

syms x1 x2 g l f m s

syms x1 x2 g l f m s 

J = jacobian([x2; -(g/l)*sin(x1) - (f/m)*x2 ],[x1 x2]); 

x1 = pi/6; 

x2 = 0; 

A = eval(J); 

det(s*eye(2)-A); 

J, A, vp = eig(A)

 

Valores propios:

$\lambda _{1,2}=-\frac{f}{2m}\pm \frac{1}{2}\sqrt{\left( \frac{f}{m} \right) ^2-\frac{8\sqrt{3}g}{l}}$

Debido al segundo término dentro del radical, el radical será siempre menor que $(f/m)$ y las raíces siempre serán negativas, por lo que el sistema es asintóticamente estable en ese punto de equilibrio. 

Para $\mathbf{x}_0=(2\pi /3,0)$ se tiene:

$\mathbf{A}=\left[ \begin{matrix} 0&  1\\ \frac{g}{2l}&  -\frac{f}{m}\\\end{matrix} \right]$

Los valores propios son:

$\lambda _{1,2}=-\frac{f}{2m}\pm \frac{1}{2}\sqrt{\left( \frac{f}{m} \right) ^2+\frac{2g}{l}}$

En este caso hay una raíz en el semiplano derecho y una en el semiplano izquierdo, por lo que el punto de equilibrio es inestable (punto de silla), tal y como se observa en el retrato de fase de arriba.

Si se aplica el método directo de Lyapunov lo más natural (mas no necesario) es utilizar la energía total del péndulo como función definida positiva (localmente, dado que es igual a cero en varios puntos donde $\cos x_1=1$):

$E=E_c+E_p=\frac{ml^2x_{2}^{2}}{2}+mg(l-l\cos x_1)$

Función candidata de Lyapunov (es importante aclarar que esta función no necesariamente debe corresponder a una magnitud física (propiedad física que puede medirse) y ese es el poder del método directo de Lyapunov):

$V(\mathbf{x})=\frac{E}{ml^2}=\left( \frac{g}{l} \right) \left( 1-\cos x_1 \right) +\frac{1}{2}x_{2}^{2}$

Derivando:

$\dot{V}(\mathbf{x})=\left( \frac{g}{l} \right) \dot{x}_1\sin x_1+x_2\dot{x}_2=-\left( \frac{k}{m} \right) x_{2}^{2}\le 0, \forall x_1$

Se puede asegurar sólo que el punto de equilibrio es estable, ya que la función $\dot{V}(\mathbf{x})$ es semidefinida negativa (la función es igual a cero no solo en el origen, sino para todos los valores de $x_1$). La función $V(\mathbf{x})$ falla debido a que el punto de equilibrio realmente es asintóticamente estable (como se vio anteriormente), por lo que se debe buscar otra función.

Si no hay fricción $(f=0)$ entonces las raíces están sobre el eje imaginario y el primer teorema de Lyapunov no puede afirmar nada. El retrato de fase, sin embargo, muestra que se tiene un centro:

Discusión y verificación

Los cálculos matemáticos y la información de la simulación y el retrato de fase muestran una coherencia de los resultados, lo cual da más confianza al análisis. 

Se puede observar que el modelo tiene dos puntos de equilibrio, un foco estable en 30° y un punto de silla (inestable) en 150° (2.6180 rad). Tanto el primer teorema de Lyapunov como el retrato de fase muestran esos comportamientos. 

Lo más interesante es observar lo que sucede cuando el péndulo parte de un ángulo inferior o superior a 150°: cuando es inferior converge al punto de equilibrio estable. 

Si el modelo no tiene fricción se obtienen valores propios sobre el eje imaginario y el primer teorema de Lyapunov no permite obtener el tipo de punto de equilibrio, aunque el retrato de fase si permite ver que se tiene un centro. 

De esta manera, el primer teorema de Lyapunov da las pautas para el análisis de estabilidad de sistemas lineales.