Funciones de Lyapunov: demostrar estabilidad

Estudiar la estabilidad sin resolver todas las trayectorias

Los capítulos 33 a 36 construyeron varias rutas para comprender la dinámica alrededor de un equilibrio. En una dimensión, el signo del flujo clasifica los puntos estacionarios. En los sistemas lineales, los valores propios organizan los modos de evolución. En los no lineales, el jacobiano da una aproximación local cuando el equilibrio es hiperbólico, mientras que las nulclinas y las regiones invariantes añaden restricciones geométricas que pueden extenderse lejos de ese punto (Hirsch et al., 2013).

Queda otra pregunta. Supongamos que el sistema no puede resolverse, que la linealización es insuficiente o que se busca una afirmación de estabilidad válida en una región explícita. Interesa entonces saber si puede demostrarse que las trayectorias permanecen cerca de un equilibrio, o incluso que terminan acercándose a él, sin conocer sus fórmulas temporales.

El método directo de Lyapunov responde mediante una función escalar del estado. Para un sistema autónomo

con el equilibrio trasladado al origen, se busca una función $V:\mathbb R^n\rightarrow\mathbb R$ que asigne a cada estado un número capaz de ordenar la evolución. La idea recuerda a una energía que disminuye por disipación, aunque ambas cosas no deben identificarse. Una función de Lyapunov puede coincidir con una energía física, ser una combinación artificial de variables o surgir de una ecuación matricial sin interpretación energética independiente. Su papel es matemático, porque se trata de construir una cantidad cuya geometría y cuya variación restringen las trayectorias (Khalil, 2002; LaSalle y Lefschetz, 1961).

La ventaja es importante, porque una sola desigualdad sobre $\dot V$ puede controlar una familia continua de condiciones iniciales. La dificultad se desplaza. En lugar de resolver todas las órbitas, hay que encontrar una función adecuada y demostrar qué propiedades tiene realmente.

Estabilidad, atracción y estabilidad asintótica

Consideremos el equilibrio $z=0$. Es estable en el sentido de Lyapunov si para todo $\varepsilon>0$ existe $\delta>0$ tal que

para todo $t\ge0$ mientras la solución esté definida. La definición no exige que el estado se acerque al equilibrio. Exige que una perturbación suficientemente pequeña permanezca controlada, de modo que un centro lineal ideal puede ser estable sin ser atractivo.

El origen es localmente atractivo si existe una vecindad desde la cual $z(t)\rightarrow0$ cuando $t\rightarrow\infty$. Es asintóticamente estable si es estable y atractivo a la vez. La conjunción importa, porque converger desde algunos estados no basta para establecer la estabilidad en el sentido anterior.

Hablaremos de estabilidad asintótica global cuando la atracción y las condiciones de existencia pertinentes se extiendan a todo el espacio de estados considerado. El adjetivo global no debe añadirse porque una simulación desde muchas condiciones parezca converger. Requiere un argumento que abarque el dominio completo, o la región que se haya especificado como dominio del modelo. Las funciones de Lyapunov demuestran estas propiedades comparando el valor de $V$, que mide de algún modo la distancia al equilibrio, con la derivada de $V$ a lo largo del campo.

Funciones positivas definidas

Una función continua $V(z)$ es positiva definida alrededor del origen si $V(0)=0$ y $V(z)>0$ para todo $z\neq0$ en la región considerada. La condición no exige que $V$ sea una distancia métrica, ya que puede contener términos cruzados, potencias distintas o escalas diferentes. Lo que importa es que el origen sea su mínimo estricto.

La función $V(x,y)=x^2+3y^2$ es positiva definida. También lo es

aunque el término $xy$ impide reconocerlo a simple vista. Para una forma cuadrática $V(z)=z^{\mathsf T}Pz$ con $P=P^{\mathsf T}$, la positividad definida equivale a que todos los valores propios de $P$ sean positivos. Si $\lambda_{\min}(P)$ y $\lambda_{\max}(P)$ son sus valores propios extremos, entonces

Estas desigualdades conectan los niveles de $V$ con las bolas ordinarias del espacio de estados. Esa conexión es la pieza geométrica que permite traducir una disminución de $V$ en una afirmación sobre estabilidad.

Otra propiedad útil es la no acotación radial, también llamada coercividad en este contexto, según la cual $V(z)\rightarrow\infty$ cuando $\|z\|\rightarrow\infty$. Una función positiva definida puede carecer de ella. Para las conclusiones globales la diferencia puede ser decisiva, porque hay que impedir que un nivel finito de $V$ se prolongue hasta estados arbitrariamente lejanos.

Derivada a lo largo de las trayectorias

Si $V$ es diferenciable y $z(t)$ satisface $\dot z=F(z)$, la regla de la cadena da

cantidad que suele abreviarse como $\dot V(z)=\nabla V(z)^{\mathsf T}F(z)$. Aquí no se deriva $V$ respecto de una variable temporal explícita. Se calcula cómo cambia su valor cuando el estado se mueve según el campo. En los sistemas no autónomos habría que añadir el término parcial $\partial V/\partial t$.

Si $\dot V(z)<0$ fuera del equilibrio, cada trayectoria reduce estrictamente el valor de $V$ mientras permanece en la región donde vale la desigualdad. Las curvas de nivel de $V$ se comportan entonces como superficies anidadas que la evolución atraviesa hacia niveles menores. El razonamiento es más fuerte que observar una serie temporal decreciente para una condición inicial, porque la desigualdad es una propiedad del campo en una región completa.

El teorema directo básico y por qué demuestra estabilidad

Una versión local del teorema directo puede formularse así. Supongamos que el origen es un equilibrio de un sistema autónomo con soluciones únicas y que existe una función continuamente diferenciable $V$ en una vecindad del origen tal que $V(0)=0$, $V(z)>0$ para $z\neq0$ y $\dot V(z)\le0$. Entonces el origen es estable en el sentido de Lyapunov.

El argumento puede verse geométricamente. Tomemos una esfera pequeña de radio $\varepsilon$ alrededor del origen. Como $V$ es positiva en esa frontera compacta, alcanza en ella un mínimo positivo,

Por continuidad y porque $V(0)=0$, puede elegirse $\delta<\varepsilon$ de manera que $\|z\|<\delta$ implique $V(z)<m_\varepsilon$. Si una trayectoria comienza dentro de esa bola y $\dot V\le0$, entonces $V(z(t))\le V(z(0))<m_\varepsilon$. La trayectoria no puede alcanzar la esfera $\|z\|=\varepsilon$, donde $V\ge m_\varepsilon$, de modo que queda confinada dentro de ella. Esa es precisamente la estructura $\varepsilon$-$\delta$ de la estabilidad.

Si la condición se refuerza a $\dot V(z)<0$ para $z\neq0$, de modo que $-\dot V$ sea positiva definida, se obtiene bajo las hipótesis regulares habituales la estabilidad asintótica local. La disminución estricta impide que una trayectoria no estacionaria permanezca indefinidamente sobre un mismo nivel de $V$ (Slotine y Li, 1991; Teschl, 2012).

Hay que evitar extender la conclusión más allá de lo demostrado. Si las propiedades de $V$ se probaron solo dentro de una región $D$, el teorema da una conclusión local asociada a esa región. Afirmar la estabilidad asintótica global exige, además, que las soluciones puedan prolongarse y que la geometría de $V$ controle estados arbitrariamente lejanos, por ejemplo mediante una función radialmente no acotada.

Una función cuadrática que no es una energía física

Retomemos un sistema lineal bidimensional del tipo estudiado en el capítulo 34,

Sus valores propios son $-1$ y $-2$, de modo que el análisis modal ya anticipa la estabilidad asintótica del origen. El sistema servirá para aprender otra ruta de demostración.

La función euclídea más obvia sería $V_0(x,y)=x^2+y^2$. Su derivada es

La expresión no es negativa definida, porque el término cruzado puede volverla positiva en ciertos estados. En $(1;\ -0{,}1)$, por ejemplo, vale $0{,}14$. El sistema es estable, pero esta candidata no lo demuestra. Un candidato fallido deja la estabilidad sin refutar y muestra solo que la función elegida no da el certificado buscado.

Para los sistemas lineales de Hurwitz hay una construcción sistemática. Se elige una matriz simétrica positiva definida $Q$ y se busca una matriz simétrica $P$ que satisfaga la ecuación de Lyapunov continua,

Con $Q=I$ y $P$ simétrica, las tres ecuaciones independientes dan

cuyos valores propios, $\lambda_{\pm}=(3\pm\sqrt5)/4$, son positivos. Por tanto,

es positiva definida y, a lo largo de las trayectorias,

La derivada es negativa definida, de modo que la estabilidad asintótica queda demostrada sin integrar una sola condición inicial (figura 37.1).

La función no se dedujo de una energía física del sistema, porque se construyó para certificar una propiedad dinámica. La distinción es especialmente importante fuera de la mecánica. Llamar «energía» a cualquier función decreciente añade una interpretación que la matemática no contiene.

Por qué la ecuación matricial funciona para cualquier sistema lineal de Hurwitz

El ejemplo anterior no depende de una coincidencia particular. Sea $\dot z=Az$ un sistema lineal cuyo espectro está estrictamente en el semiplano izquierdo. Para cualquier matriz simétrica positiva definida $Q$, definamos

La integral converge porque $e^{At}$ decae exponencialmente cuando $A$ es de Hurwitz. Para todo vector no nulo $v$,

de modo que $P$ es positiva definida. La identidad que satisface se obtiene al observar que

e integrar desde cero hasta infinito,

La función $V(z)=z^{\mathsf T}Pz$ satisface, por tanto, $\dot V=-z^{\mathsf T}Qz<0$ fuera del origen. Para los sistemas lineales continuos, la estabilidad espectral del capítulo 34 y la existencia de una función cuadrática estricta son dos maneras estrechamente relacionadas de expresar la misma estructura.

La dirección recíproca también es informativa. Si existen matrices $P>0$ y $Q>0$ con $A^{\mathsf T}P+PA=-Q$, ningún modo puede tener parte real no negativa, porque la función cuadrática decrece estrictamente para todo estado no nulo. La ecuación de Lyapunov es así un criterio algebraico de estabilidad para los sistemas lineales, además de un procedimiento para fabricar funciones.

La equivalencia desaparece en su forma simple al pasar a los sistemas no lineales. En ellos una función de Lyapunov puede existir sin ser cuadrática. Una función cuadrática útil cerca del equilibrio puede, además, dejar de controlar la dinámica lejos de él. Por eso hay que seguir separando las conclusiones locales de las globales.

De la disminución cualitativa a una cota de velocidad

El mismo ejemplo permite obtener algo más que la convergencia. De $\lambda_{\min}\|z\|^2\le V(z)\le\lambda_{\max}\|z\|^2$ y $\dot V=-\|z\|^2$ se sigue $\|z\|^2\ge V/\lambda_{\max}$, de modo que

Al integrar esta desigualdad escalar se obtiene $V(t)\le V(0)\exp(-t/\lambda_{\max})$, con $\lambda_{\max}=(3+\sqrt5)/4\approx1{,}3090$. Al combinar de nuevo las desigualdades espectrales resulta

donde el factor $\sqrt{\lambda_{\max}/\lambda_{\min}}$ vale unos $2{,}618$ y la tasa garantizada es de unos $0{,}382$.

La cota no coincide necesariamente con la trayectoria exacta ni con la tasa dominante del sistema. El modo más lento decae como $e^{-t}$, mucho más rápido que la tasa garantizada. La cota es un límite construido a partir de la función escogida. Distintas funciones de Lyapunov pueden dar cotas más o menos conservadoras para la misma dinámica (figura 37.2).

La distinción entre demostrar la convergencia y estimar su velocidad será recurrente. Una función que prueba la estabilidad puede ser excelente para una conclusión cualitativa y mediocre para cuantificar tiempos de relajación.

Un sistema no lineal con certificado global sencillo

El método no se limita a los sistemas lineales. Consideremos

cuyo único equilibrio es el origen. Tomemos $V(x,y)=\tfrac12(x^2+y^2)$, positiva definida y radialmente no acotada. Su derivada es

negativa para todo $(x,y)\neq(0,0)$. En consecuencia, el origen es asintóticamente estable. Los subniveles $\Omega_c=\{(x,y):V(x,y)\le c\}$ son, además, compactos y positivamente invariantes, porque una trayectoria que comienza en $\Omega_c$ no puede pasar a un nivel mayor si $V$ no aumenta. Como hay subniveles de cualquier tamaño y el campo es suave en todo el plano, la función da un argumento global para la convergencia al origen (figura 37.3).

La conclusión pertenece al modelo matemático completo y no a su linealización. Cerca del origen el término $-x^3$ es de orden superior y el jacobiano también predice estabilidad, mientras que la función de Lyapunov extiende el control más allá de una vecindad infinitesimal.

Derivada semidefinida y convergencia

El signo no estricto requiere cuidado. En el sistema $\dot x=-x$, $\dot y=0$, con $V(x,y)=\tfrac12(x^2+y^2)$, se obtiene $\dot V=-x^2\le0$. La función nunca aumenta, así que el origen es estable. Sin embargo, $y(t)=y_0$ permanece constante, de modo que la trayectoria no converge al origen si $y_0\neq0$. La condición $\dot V\le0$ no permite concluir la estabilidad asintótica.

El conjunto donde se anula la derivada es $\{(x,y):x=0\}$, que en este ejemplo contiene toda una familia de equilibrios. En otros sistemas, el conjunto donde $\dot V=0$ puede contener puntos que no son invariantes, de los que las trayectorias salen. Determinar qué parte de ese conjunto puede contener la dinámica límite exige una herramienta adicional, el principio de invariancia de LaSalle, tema del capítulo 38. El método directo con derivada negativa definida ya permite demostrar mucho. LaSalle amplía lo que puede concluirse cuando la disminución no es estricta.

Regiones de atracción mediante conjuntos de nivel

En un sistema no lineal, una función de Lyapunov puede ser útil aunque sus propiedades se conozcan solo localmente. Supongamos que $V$ es positiva definida en una región $D$ y que $\dot V<0$ en $D\setminus\{0\}$. Supongamos además que el campo es localmente Lipschitz en el dominio abierto $D$ que contiene al origen y que $V$ es continuamente diferenciable. Si existe $c>0$ tal que el subnivel $\Omega_c=\{z\in D:V(z)\le c\}$ es compacto y está contenido en $D$, entonces es positivamente invariante. En su frontera $V=c$, la derivada negativa impide que una trayectoria pase a $V>c$. La compacidad permite prolongar las soluciones hacia adelante y, junto con la disminución estricta de $V$ fuera del origen, garantiza la convergencia al equilibrio para los estados iniciales de $\Omega_c$.

Se obtiene así una estimación interna de la cuenca de atracción, que no tiene por qué ser la cuenca completa. Una trayectoria que comienza fuera de $\Omega_c$ puede converger también al equilibrio, aunque la prueba no lo establezca. El teorema informa dónde se demostró una propiedad, lo que es distinto de indicar dónde deja de ser cierta. Por eso buscar mejores funciones de Lyapunov puede ampliar la región certificada sin cambiar el sistema, ya que la geometría de la prueba no coincide con la geometría completa de la dinámica.

Cómo se encuentra una función de Lyapunov

No existe un algoritmo universal que produzca una función útil para todo sistema no lineal, aunque en ejemplos sencillos funcionan varias estrategias. Puede usarse una cantidad físicamente conservada o disipada, cuando realmente existe y su interpretación está justificada. En un oscilador amortiguado, por ejemplo, la energía mecánica es una candidata natural.

Cerca de un equilibrio pueden proponerse formas cuadráticas $V(z)=z^{\mathsf T}Pz$. En los sistemas lineales de Hurwitz, la ecuación de Lyapunov da una construcción sistemática para cualquier $Q>0$. En los no lineales, la misma función puede servir localmente si la parte lineal domina los términos de orden superior en una vecindad suficientemente pequeña. También puede aprovecharse la estructura algebraica del modelo, mediante sumas de cuadrados, términos que cancelan productos cruzados, funciones logarítmicas en sistemas positivos o integrales construidas a partir del campo.

Lo que no debe hacerse es elegir una función después de mirar una sola trayectoria y tratar su descenso aparente como una demostración general. La prueba requiere establecer la positividad y el signo de $\dot V$ en toda una región, más allá de los puntos que muestrea una simulación. En problemas de mayor dimensión hay métodos computacionales para buscar certificados, entre ellos la programación semidefinida y las técnicas de sumas de cuadrados para sistemas polinomiales. Esas herramientas aparecerán más adelante, cuando se disponga del marco de optimización necesario. Su automatización no cambia el principio lógico. El resultado útil sigue siendo una desigualdad certificada bajo hipótesis explícitas.

Verificación computacional sin sustituir la demostración

El bloque siguiente comprueba numéricamente la construcción del ejemplo lineal, la cota exponencial y el signo de la derivada en el ejemplo no lineal.

import numpy as np
from scipy.linalg import expm, solve_continuous_lyapunov

A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
P = solve_continuous_lyapunov(A.T, -np.eye(2))
print(P)
print(np.linalg.eigvalsh(P))

lmax = np.linalg.eigvalsh(P).max()
peor = 0.0
for ang in np.linspace(0.0, 2*np.pi, 13)[:-1]:
    z0 = 2.0 * np.array([np.cos(ang), np.sin(ang)])
    for t in np.linspace(0.0, 6.0, 121):
        z = expm(A * t) @ z0
        peor = max(peor, (z @ P @ z) / ((z0 @ P @ z0) * np.exp(-t / lmax)))
print(peor)

x, y = np.meshgrid(np.linspace(-3, 3, 601), np.linspace(-3, 3, 601))
derivada = x*(-x**3 + y) + y*(-x - y)
print(np.max(derivada[(x != 0) | (y != 0)]))

La salida reproduce la matriz $P$ y sus valores propios, $0{,}191$ y $1{,}309$. El cociente máximo entre $V(t)$ y su cota vale 1, alcanzado en $t=0$, de modo que ninguna de las doce trayectorias supera la cota. En la malla del ejemplo no lineal, la derivada es negativa en todos los puntos distintos del origen. Estas comprobaciones verifican, sin fundamentar lógicamente el resultado. Una malla finita no demuestra una desigualdad para todos los puntos del plano. Aquí la demostración proviene del álgebra, mientras que la simulación controla que el código, las figuras y las cifras sean coherentes con ella.

Límites conceptuales y errores frecuentes

Una función de Lyapunov demuestra propiedades del sistema especificado, sin demostrar que el sistema real tenga el mecanismo descrito por el modelo. Si una función decrece en un modelo de síntomas, por ejemplo, eso no revela ninguna «energía psicopatológica» ni un principio físico subyacente.

Tampoco debe confundirse una función de Lyapunov con un exponente de Lyapunov. Una función de Lyapunov es una función escalar construida para demostrar la estabilidad o acotar la evolución. Los exponentes, que estudiaremos en el bloque dedicado al caos, cuantifican tasas asintóticas de crecimiento o de contracción de perturbaciones en las direcciones del espacio tangente. Ambos conceptos comparten un apellido histórico y definiciones distintas.

Que $\dot V<0$ no significa que cada variable decrezca. En el sistema lineal central, una componente puede crecer temporalmente en magnitud mientras la combinación cuadrática $V$ disminuye, porque el certificado controla una propiedad agregada del estado sin garantizar la monotonía de cada componente. El método es, por último, suficiente, aunque no necesario en la forma concreta que se haya intentado. Fracasar al buscar una función simple no demuestra inestabilidad. Los teoremas conversos establecen, bajo hipótesis apropiadas, que muchos sistemas asintóticamente estables admiten funciones de Lyapunov, aunque encontrarlas de forma útil puede ser difícil.

Consolidación

El capítulo 33 definió con precisión la estabilidad y la atracción. Los capítulos 34 y 35 mostraron cómo los valores propios y la linealización responden localmente a esas preguntas. El método directo de Lyapunov introduce otra estrategia, que consiste en construir una función escalar cuyos niveles se comportan como barreras geométricas y cuya derivada sobre el campo ordena la evolución.

Una función positiva definida con $\dot V\le0$ demuestra la estabilidad. Si $\dot V$ es negativa definida, se obtiene, bajo las hipótesis regulares correspondientes, la estabilidad asintótica. Extender el resultado globalmente exige controlar también la existencia de las soluciones y el comportamiento de $V$ lejos del equilibrio. La no acotación radial es una condición suficiente frecuente para ello.

El ejemplo lineal mostró que la función útil no tiene por qué ser la norma euclídea ni una energía física. La ecuación matricial de Lyapunov produjo una forma cuadrática cuya derivada es exactamente $-\|z\|^2$. El ejemplo no lineal mostró que una candidata sencilla puede certificar la estabilidad global aunque el campo contenga términos de orden superior. La derivada semidefinida deja abierta la pregunta sobre lo que ocurre dentro del conjunto donde $\dot V=0$. El capítulo 38 la responderá mediante el principio de invariancia de LaSalle y usará los conjuntos de nivel para estimar regiones de atracción de mayor alcance.

Ejercicios

Ejercicio 37.1. Para $A=\begin{pmatrix}-1&0\\1&-2\end{pmatrix}$, (a) resuelva la ecuación de Lyapunov $A^{\mathsf T}P+PA=-I$. (b) Compruebe que $P$ es positiva definida. (c) Escriba la función de Lyapunov resultante y su derivada.

Ejercicio 37.2. Para el oscilador amortiguado $\dot x=y$, $\dot y=-x-0{,}5y$, (a) calcule la derivada de la energía $E=\tfrac12(x^2+y^2)$. (b) ¿Qué permite concluir el método directo con esa función? (c) Use la ecuación de Lyapunov con $Q=I$ para obtener una función con derivada negativa definida.

Ejercicio 37.3. En el sistema lineal central, (a) encuentre un estado donde la derivada de $V_0=x^2+y^2$ sea positiva. (b) Explique por qué ese resultado no contradice la estabilidad asintótica del origen.

Ejercicio 37.4. Con los valores propios de la matriz $P$ del capítulo, (a) calcule la constante y la tasa de la cota exponencial para $\|z(t)\|$. (b) Compárelas con el decaimiento real del modo más lento. (c) ¿Qué propiedad de $P$ hace que la cota sea conservadora?

Ejercicio 37.5. Para la ecuación escalar $\dot x=-x+x^3$, (a) use $V=x^2/2$ para encontrar la región donde $\dot V<0$. (b) Obtenga la mayor estimación de la cuenca de atracción del origen que da ese argumento. (c) Compárela con la cuenca exacta.

Ejercicio 37.6. Para $\dot x=-x+y^2$, $\dot y=-y$, (a) calcule la derivada de $V=\tfrac12(x^2+y^2)$. (b) Muestre que $\dot V<0$ fuera del origen en la franja $|x|<1$ y obtenga la región certificada por los conjuntos de nivel. (c) Compruebe que la cuenca de atracción real es todo el plano.

Referencias

Hirsch, M. W., Smale, S. y Devaney, R. L. (2013). Differential equations, dynamical systems, and an introduction to chaos (3.ª ed.). Academic Press.

Khalil, H. K. (2002). Nonlinear systems (3.ª ed.). Prentice Hall.

LaSalle, J. P. y Lefschetz, S. (1961). Stability by Liapunov's direct method with applications. Academic Press.

Slotine, J.-J. E. y Li, W. (1991). Applied nonlinear control. Prentice Hall.

Teschl, G. (2012). Ordinary differential equations and dynamical systems. American Mathematical Society.