M5-T1 — Lyapunov Stability

Learning objectives

After this lesson, you should be able to:

• State Lyapunov's theorem and the invariance principle and say what each proves.
• Construct a quadratic Lyapunov function for a stable linear system by solving the Lyapunov equation.
• Estimate a region of attraction as the largest level set inside the region where a bound on V̇ holds.
• Compare a guaranteed region with the true region found by simulation and explain the gap.

The intuition

UAV M5-T5: Attitude control used a nonlinear attitude law and the Navigation track's control barrier functions are a cousin of the method here. Neither says where the guarantee ends. This lesson supplies the tool: an energy argument that holds for the nonlinear system as it is, and says how big the starting error may be.

The ball in the bowl

A ball in a bowl with friction comes to rest at the bottom. You can prove it without solving the motion: the energy of the ball falls at every instant that it moves, and it can only stop at the bottom. A Lyapunov function is any quantity of the state that plays the part of energy: positive away from the equilibrium, zero at it, and falling along every trajectory. If one exists the equilibrium is stable, and if the fall is strict the system returns to it. The function does not have to be physical energy, and finding one is the art of the method.

For a linear system there is always one

If $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}$ is stable, a quadratic $V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}$ works, with $\mathbf{P}$ found by solving a linear equation. For Hopper's pitch loop it is the matrix whose level sets are tilted ellipses in the plane of pitch and pitch rate. The same function can then be tried on the nonlinear loop, where it keeps working inside the region where the nonlinearity has not yet acted.

A guaranteed region, and its price

The gimbal saturates at 0.15 rad. As long as the commanded gimbal angle stays within it the loop is exactly the linear loop, and V falls. The largest level set of V that fits inside that region is a region from which the state is guaranteed to return: a region of attraction estimate. For Hopper it is small: pitch up to 0.106 rad. The loop actually recovers from a start of 6.0 rad, almost sixty times more. The theorem is right and the estimate honest, but a Lyapunov function that was not built for the problem gives a conservative answer, and that is worth knowing before trusting one.

A Lyapunov function proves stability of the nonlinear loop as it is, with no linearization, and a level set of it is a guaranteed region of attraction. The region is usually much smaller than the true one.

Definitions

Equilibrium

A state $\mathbf{x}^*$ with $f(\mathbf{x}^*) = 0$. The origin is assumed, after a shift of coordinates.

Positive definite $V$

$V(0) = 0$ and $V(\mathbf{x}) > 0$ for every other state in the region considered.

Derivative along trajectories $\dot{V}$

$\dot{V} = \nabla V^{\mathsf{T}} f(\mathbf{x})$: the rate at which V changes along a solution, computed without solving.

Region of attraction

The set of initial states from which the state tends to the equilibrium.

Level set

$\Omega_c = \{\mathbf{x} : V(\mathbf{x}) \le c\}$. If $\dot{V} < 0$ on $\Omega_c$ minus the origin, no trajectory can leave it.

The direct method is a pair of statements.

When $\dot{V}$ is only nonpositive, LaSalle's invariance principle still gives convergence if no trajectory other than the origin can stay forever in the set where $\dot{V} = 0$.

Derivation: why the Lyapunov equation works, and the region estimate

For the linear loop and then for a saturated one.

1. $$V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}, \qquad \dot{V} = \dot{\mathbf{x}}^{\mathsf{T}}\mathbf{P}\mathbf{x} + \mathbf{x}^{\mathsf{T}}\mathbf{P}\dot{\mathbf{x}}$$

   The product rule for a quadratic form.

2. $$\dot{V} = \mathbf{x}^{\mathsf{T}}(\mathbf{A}^{\mathsf{T}}\mathbf{P} + \mathbf{P}\mathbf{A})\mathbf{x}$$

   Substitute $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}$.

3. $$\dot{V} = -\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x} < 0 \text{ for } \mathbf{x} \ne 0$$

   By the Lyapunov equation. The existence of P for every stable A comes from the integral $\mathbf{P} = \int_0^\infty e^{\mathbf{A}^{\mathsf{T}}t}\mathbf{Q}e^{\mathbf{A}t}dt$, the cost of a trajectory.

The region estimate for the saturated loop

1. $$u = -\mathbf{K}\mathbf{x} \text{ is the loop whenever } |\mathbf{K}\mathbf{x}| \le u_{\max}$$

   Inside the strip where the commanded gimbal is not saturated the closed loop is exactly linear, and $\dot{V} = -\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x}$ holds.

2. $$\max_{V(\mathbf{x}) \le c} |\mathbf{K}\mathbf{x}|^2 = c\;\mathbf{K}\mathbf{P}^{-1}\mathbf{K}^{\mathsf{T}}$$

   The largest value of a linear function over an ellipsoid $\mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x} \le c$ (Cauchy-Schwarz in the metric of P).

3. $$c^* = \dfrac{u_{\max}^2}{\mathbf{K}\mathbf{P}^{-1}\mathbf{K}^{\mathsf{T}}}$$

   The largest level set inside the strip. On it the loop is linear, V falls, and no trajectory can leave, so every start in it converges: a guaranteed region of attraction.

A level set of V on which V̇ is negative cannot be left. The estimate of the region of attraction is the biggest such level set; it is a sufficient condition, and how far it falls short of the true region depends on the choice of V.

A worked example

Hopper's pitch loop with the PD gains of module 3 ($K_p = 0.8405$, $K_d = 0.2854$) and the gimbal limit 0.15 rad.

1. The loop. $\ddot{\theta} = 0.49\,\theta + 19.62\,\mathrm{sat}(-K_p\theta - K_d q)$. In the linear region its poles are −2.80 ± 2.86j.
2. The function. Solve $\mathbf{A}^{\mathsf{T}}\mathbf{P} + \mathbf{P}\mathbf{A} = -\mathbf{I}$: $\mathbf{P} = \begin{bmatrix}1.693 & 0.031\\ 0.031 & 0.095\end{bmatrix}$.
3. The region. The largest level set inside the strip $|\mathbf{K}\mathbf{x}| \le 0.15$ is $V \le 0.01891$. On the axis of pitch alone it contains up to 0.106 rad; on the axis of pitch rate alone, up to 0.446 rad/s.
4. Check. On a grid of states inside the level set the largest value of the ratio $\dot{V}/|\mathbf{x}|^2$ is −1, as the theory says.
5. The truth. Simulating the saturated loop, it recovers from a pitch of up to 6.00 rad (the gimbal at its limit can hold off the instability up to $19.62 \cdot 0.15 / 0.49 = 6.0$ rad) and from a rate of up to 4.20 rad/s. The guaranteed region is 57 times smaller in pitch than the real one.

import math
import gnckit as g

lam2, b, dmax = 0.49, 19.62, 0.15
kp, kd = 0.8405, 0.2854
A = [[0.0, 1.0], [lam2 - b * kp, -b * kd]]
print('closed-loop poles of the linear pitch loop:', [f'{complex(p):.3f}' for p in g.mat_eigvals(A)])

# Lyapunov function V = x'Px with A'P + PA = -I
P = g.lyap(A, g.mat_eye(2))
print(f'P = [[{P[0][0]:.4f}, {P[0][1]:.4f}], [{P[1][0]:.4f}, {P[1][1]:.4f}]]')

# the linear law is valid while |K x| <= dmax; the largest level set V <= c inside that strip
Pi = g.mat_inv(P)
K = [kp, kd]
kPk = sum(K[i] * Pi[i][j] * K[j] for i in range(2) for j in range(2))
c = dmax ** 2 / kPk
print(f'level set V <= {c:.5f} stays inside the linear region |K x| <= {dmax}')
print(f'on q = 0 it contains |theta| <= {math.sqrt(c / P[0][0]):.3f} rad')

def f(t, s):
    u = max(-dmax, min(dmax, -(kp * s[0] + kd * s[1])))
    return [s[1], lam2 * s[0] + b * u]

def recovers(th0, q0):
    xs = g.integrate(f, [th0, q0], 15.0, 0.002)[1]
    return abs(xs[-1][0]) < 1e-3

# Vdot = -|x|^2 inside the level set for the linear law; check on a grid of states with V <= c
worst = -1e9
for i in range(-60, 61):
    for j in range(-60, 61):
        x = [i * 0.01, j * 0.02]
        v = sum(x[a] * P[a][bb] * x[bb] for a in range(2) for bb in range(2))
        if 0 < v <= c:
            d = g.mat_vec(A, x)
            vd = 2 * sum(x[a] * P[a][bb] * d[bb] for a in range(2) for bb in range(2))
            worst = max(worst, vd / sum(t * t for t in x))
print(f'largest Vdot / |x|^2 over the level set: {worst:.4f} (the theory says -1)')

# true region: largest initial pitch on q = 0 and largest initial pitch rate on theta = 0 that recover
lo, hi = 0.0, 8.0
for _ in range(40):
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if recovers(mid, 0.0) else (lo, mid)
print(f'simulation with the 0.15 rad limit: recovers from theta0 up to {lo:.3f} rad at q0 = 0')
lo, hi = 0.0, 8.0
for _ in range(40):
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if recovers(0.0, mid) else (lo, mid)
print(f'                                    recovers from q0 up to {lo:.3f} rad/s at theta0 = 0')
th = math.sqrt(c / P[0][0]); qq = math.sqrt(c / P[1][1])
print(f'the Lyapunov estimate guarantees theta0 up to {th:.3f} rad and q0 up to {qq:.3f} rad/s')

closed-loop poles of the linear pitch loop: ['-2.800-2.857j', '-2.800+2.857j']
P = [[1.6930, 0.0312], [0.0312, 0.0949]]
level set V <= 0.01891 stays inside the linear region |K x| <= 0.15
on q = 0 it contains |theta| <= 0.106 rad
largest Vdot / |x|^2 over the level set: -1.0000 (the theory says -1)
simulation with the 0.15 rad limit: recovers from theta0 up to 6.002 rad at q0 = 0
                                    recovers from q0 up to 4.202 rad/s at theta0 = 0
the Lyapunov estimate guarantees theta0 up to 0.106 rad and q0 up to 0.446 rad/s

In robotics

Lyapunov arguments are how nonlinear flight laws are justified, and how their limits are stated.

• Proving a nonlinear law. The attitude and trajectory laws of the next lessons (dynamic inversion, sliding mode, adaptation) are all stability-proven by a V built from the error. Writing the proof is how one finds out which assumption the law needs.
• Saturation. The gimbal limit makes every real loop nonlinear. The quadratic estimate above is the simplest way to bound the start from which the unsaturated analysis applies, and a flight rule for the largest permitted attitude error can be taken from it.
• Better functions. The conservatism can be reduced by searching for V by computer (sums of squares, or the linear-matrix-inequality methods in the references), and by allowing the level set to extend into the saturated region with a sector bound on the nonlinearity.
• Adaptive laws. The adaptation law of M5-T5: Adaptive Control and Its Risks is derived by adding a term in the parameter error to V and choosing the update so that $\dot{V}$ stays nonpositive.
• Limits of the method. It is sufficient, never necessary: failing to find a V proves nothing. And it has no notion of margin or noise: a system with V can be as fragile in gain as the linear one.

Matrix definiteness, which the equation relies on, is covered in Linear Algebra M6-T4: Positive Definite Matrices; the lab M5-L1: Schedule the Gains Through the Burn uses it to hold margins through the burn.

Common misconceptions

Key takeaways

1. A Lyapunov function is positive away from the equilibrium and falls along trajectories; it proves stability without solving the motion.
2. For a stable linear system $V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}$ with $\mathbf{A}^{\mathsf{T}}\mathbf{P} + \mathbf{P}\mathbf{A} = -\mathbf{Q}$ always works.
3. The invariance principle extends the result to V̇ ≤ 0 when no trajectory can stay where V̇ = 0.
4. The largest level set of V inside the unsaturated strip is a guaranteed region of attraction: $c^* = u_{\max}^2/\mathbf{K}\mathbf{P}^{-1}\mathbf{K}^{\mathsf{T}}$.
5. For Hopper the guaranteed pitch is 0.106 rad against a true 6.0 rad: the estimate is safe and conservative.
6. The method gives sufficient conditions only, and nothing about margins.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• K. J. Åström and R. M. Murray, Feedback Systems, 2nd edition, Princeton University Press, 2021 , the chapter on Lyapunov stability. Free to read as the authors' PDF; the book is copyrighted.
• R. Tedrake, Underactuated Robotics, the chapter on Lyapunov analysis , with region-of-attraction estimates by sums of squares. Free to read; licence CC BY-NC-ND 4.0.
• MIT 6.832 Underactuated Robotics (OpenCourseWare) , the lectures on Lyapunov functions. Free; licence CC BY-NC-SA 4.0.
• H. K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall, 2002, the chapters on Lyapunov stability and the region of attraction. Paid.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

Other courses on this site treat neighbouring topics from their own angle:

• Lab: Build a CBF Safety Filter
• Navigation M7-T2: Safety Filters: Control Barrier Functions
• UAV M5-T5: Attitude control
• UAV M5-T16: Disturbance rejection and wind compensation
• UAV M1-T6: Roll, pitch, yaw and thrust generation
• Linear Algebra M6-T7: The Moore-Penrose Pseudoinverse
• Linear Algebra M8-T9: Null-Space Control
• Robot Learning M9-T4: Hybrid Learning and Control Systems

This course leaves the following to the pages that already teach them:

• Control barrier functions and safety filters
• Geometric attitude control of a multirotor on SO(3)
• Disturbance observers for wind
• Learned controllers and residual policies