Consider a nonlinear state-space model in joint-space or task-space coordinates:

\[ \dot{x} = f(x,u), \qquad y = h(x,u), \]

where \( x \in \mathbb{R}^n \) collects robot states (e.g., joint positions and velocities), \( u \in \mathbb{R}^m \) is the control input, and \( y \) is the output (joint positions, end-effector pose, etc.). An equilibrium \( (x^\star, u^\star) \) satisfies

\[ f(x^\star, u^\star) = 0. \]

By shifting coordinates (which you have seen in the linearization lesson), we can usually assume \( x^\star = 0, u^\star = 0 \) and study the origin:

\[ \dot{x} = f(x), \qquad f(0) = 0. \]

We say that the origin is:

• Stable in the sense of Lyapunov if for every \( \varepsilon > 0 \) there exists a \( \delta(\varepsilon) > 0 \) such that \( \|x(0)\| < \delta(\varepsilon) \) implies \( \|x(t)\| < \varepsilon \) for all \( t \geq 0 \).
• Asymptotically stable if it is Lyapunov stable and, in addition, \( \lim_{t \to \infty} \|x(t)\| = 0 \) for all initial states in some neighborhood of the origin.
• Exponentially stable if there exist constants \( M \geq 1 \) and \( \alpha > 0 \) such that

\[ \|x(t)\| \leq M e^{-\alpha t} \|x(0)\| \quad \text{for all } t \geq 0 \text{ and all sufficiently small } x(0). \]

For robot control, we normally want asymptotic or even exponential stability of tracking errors, because these directly connect to convergence speed and robustness to small disturbances.

flowchart TD
  A["Robot dynamics model: dot(x)=f(x,u)"] --> B["Identify equilibrium (x*,u*)"]
  B --> C["Shift to error coordinates: e = x - x*"]
  C --> D["Choose notion: internal, BIBO, ISS"]
  D --> E["Apply appropriate test (Lyapunov, transfer function, ISS Lyapunov)"]
        

Lyapunov's direct method provides a powerful, coordinate-independent tool to verify stability without explicitly solving the differential equations. For the autonomous system \( \dot{x} = f(x) \), a Lyapunov function candidate is a scalar function \( V:\mathbb{R}^n \to \mathbb{R} \) such that:

• Positive definiteness: \( V(0) = 0 \) and \( V(x) > 0 \) for all \( x \neq 0 \).
• Radial unboundedness (optional but useful): \( V(x) \to \infty \) as \( \|x\| \to \infty \).

The derivative of \( V \) along trajectories is

\[ \dot{V}(x) = \nabla V(x)^{\top} f(x). \]

The classical Lyapunov theorems state:

• If \( V(x) \) is positive definite and \( \dot{V}(x) \leq 0 \) in a neighborhood of the origin, then the origin is Lyapunov stable.
• If \( V(x) \) is positive definite and \( \dot{V}(x) \) is negative definite (strictly negative for all \( x \neq 0 \)) in some neighborhood, then the origin is asymptotically stable.
• If, in addition, there exist constants \( c_1, c_2, c_3, c_4 > 0 \) such that

\[ c_1 \|x\|^2 \leq V(x) \leq c_2 \|x\|^2,\quad \dot{V}(x) \leq -c_3 \|x\|^2,\quad \|\nabla V(x)\| \leq c_4 \|x\| \]

then the origin is exponentially stable.

In robotics, a classical Lyapunov candidate is the total energy of a manipulator (kinetic plus potential energy) or its shaped version under feedback. For a single-DOF rigid joint with inertia \( I > 0 \), position \( q \), and velocity \( \dot{q} \), consider a PD control law (already sketched in earlier lessons):

\[ I\ddot{q} + k_d \dot{q} + k_p q = 0, \quad k_p > 0,\; k_d > 0. \]

A Lyapunov candidate is

\[ V(q,\dot{q}) = \tfrac{1}{2} I \dot{q}^2 + \tfrac{1}{2} k_p q^2. \]

Its derivative along trajectories is

\[ \dot{V} = I \dot{q} \ddot{q} + k_p q \dot{q} = \dot{q} (I\ddot{q} + k_p q) = \dot{q} (-k_d \dot{q}) = -k_d \dot{q}^2 \leq 0. \]

This shows that the origin \( (q,\dot{q}) = (0,0) \) is asymptotically stable (and in fact exponentially stable, which can be seen from the linear form of the dynamics).

Consider the linearized robot error dynamics from Lesson 3:

\[ \dot{x} = A x, \]

where \( A \in \mathbb{R}^{n \times n} \). The equilibrium \( x = 0 \) is asymptotically stable if and only if the matrix \( A \) is Hurwitz (all eigenvalues have strictly negative real part).

Lyapunov's second theorem gives the following equivalent condition: For any symmetric positive-definite matrix \( Q = Q^{\top} > 0 \), there exists a unique symmetric positive-definite matrix \( P \) such that

\[ A^{\top} P + P A = -Q. \]

If such \( P > 0 \) exists, then \( V(x) = x^{\top} P x \) is a Lyapunov function and the origin is exponentially stable.

For robotic manipulators linearized around a configuration, the matrix \( A \) depends on the inertia matrix, Coriolis terms, stiffness, and damping. Controller design (e.g., PD gains) is often chosen specifically to make the closed-loop \( A \) Hurwitz. The Lyapunov equation then certifies exponential convergence of joint or task-space error.

Lyapunov stability is an internal property: it concerns state trajectories. For controller design, we also care about how outputs react to bounded inputs. For LTI systems with transfer function \( G(s) \), BIBO stability is defined as follows:

The system is BIBO stable if for every input \( u(t) \) satisfying \( |u(t)| \leq M \) for all \( t \geq 0 \) (bounded input), there exists some constant \( N(M) > 0 \) such that the output satisfies \( |y(t)| \leq N(M) \) for all \( t \geq 0 \).

For a single-input single-output LTI system with impulse response \( h(t) \), the BIBO stability condition is:

\[ \int_{0}^{\infty} |h(t)| \, dt < \infty \quad \Longleftrightarrow \quad \text{BIBO stable}. \]

In the strictly proper, minimal, LTI case, BIBO stability is equivalent to all poles of \( G(s) \) lying strictly in the open left half-plane, which coincides with exponential internal stability.

For robot control, BIBO stability ensures that bounded joint torque commands produce bounded joint positions and velocities (and end-effector coordinates), provided we have a well-posed model and sensors. In multi-DOF robots, BIBO stability is usually studied by examining the poles of loop transfer functions around each joint or around a task-space channel.

Input-to-state stability generalizes Lyapunov stability to systems with exogenous inputs (disturbances, unmodeled dynamics). Consider

\[ \dot{x} = f(x,w), \qquad x(0) = x_0, \]

where \( w(t) \) is an input/disturbance (e.g., external forces on a manipulator, unmodeled friction). The system is ISS if there exist functions \( \beta \) of class \( \mathcal{KL} \) and \( \gamma \) of class \( \mathcal{K} \) such that for all \( t \geq 0 \)

\[ \|x(t)\| \leq \beta(\|x_0\|, t) + \gamma \!\Big( \sup_{0 \leq \sigma \leq t} \|w(\sigma)\| \Big). \]

Intuitively, ISS means:

• With no disturbance (\( w \equiv 0 \)), the state converges to the origin (like asymptotic stability).
• With bounded disturbance, the state remains bounded and its bound grows at most with the size of the disturbance.

An ISS-Lyapunov function is a function \( V(x) \) that is positive definite and radially unbounded and for which there exist class \( \mathcal{K} \) functions \( \alpha_1,\alpha_2,\alpha_3 \) and \( \gamma \) such that

\[ \alpha_1(\|x\|) \leq V(x) \leq \alpha_2(\|x\|), \]

\[ \dot{V}(x,w) \leq -\alpha_3(\|x\|) + \gamma(\|w\|). \]

For many robot controllers, we can show that the closed-loop dynamics are ISS with respect to bounded disturbances in torques, model errors, or contact perturbations. This justifies that small modeling inaccuracies do not destabilize the robot but only cause small steady-state errors.

flowchart TD
  S["Closed-loop robot model: dot(x)=f(x,w)"] --> V["Pick candidate V(x)"]
  V --> D["Compute dV/dt along trajectories"]
  D --> J["Check inequality: dV/dt <= -alpha3(||x||) + gamma(||w||)"]
  J --> R["If satisfied, conclude ISS and robustness to disturbances"]
        

In Python, we can use numpy and scipy (for Lyapunov equations) together with robotics libraries such as pinocchio or python-control to analyze stability of linearized joint-space models.

Below, we construct a simple 2-DOF linearized joint error model \( \dot{x} = A x \) and check exponential stability by solving the Lyapunov equation with a chosen positive-definite matrix \( Q \).

import numpy as np
from numpy.linalg import eigvals
from scipy.linalg import solve_continuous_lyapunov

# Example 2-DOF joint error dynamics: x = [q1, q2, dq1, dq2]
# A could come from linearizing the full rigid-body dynamics
A = np.array([
    [0.0, 0.0, 1.0, 0.0],
    [0.0, 0.0, 0.0, 1.0],
    [-50.0, 0.0, -10.0, 0.0],
    [0.0, -40.0, 0.0, -12.0],
])

# Check eigenvalues (Hurwitz if all real parts < 0)
lam = eigvals(A)
print("Eigenvalues of A:", lam)

# Choose Q > 0 (e.g., identity)
Q = np.eye(4)

# Solve A^T P + P A = -Q for P
P = solve_continuous_lyapunov(A, -Q)

# Check symmetry and definiteness of P
print("P (Lyapunov solution):\n", P)
print("Eigenvalues of P:", eigvals(P))

def V(x):
    return float(x.T @ P @ x)

def dVdt(x):
    # For linear system, dV/dt = x^T (A^T P + P A) x = -x^T Q x
    return float(-x.T @ Q @ x)

# Example trajectory simulation
def simulate(x0, dt=0.001, T=2.0):
    x = x0.copy()
    xs = [x.copy()]
    ts = [0.0]
    t = 0.0
    while t < T:
        xdot = A @ x
        x = x + dt * xdot
        t += dt
        xs.append(x.copy())
        ts.append(t)
    return np.array(ts), np.array(xs)

x0 = np.array([0.1, -0.1, 0.0, 0.0])
ts, xs = simulate(x0)

print("Initial V:", V(x0))
print("Final V:", V(xs[-1]))
      

In a robotics stack, the matrix A can be obtained by:

• Using pinocchio.computeABADerivatives to get linearization of the dynamics.
• Using python-control to create a StateSpace model and then calling Eigen or SciPy routines for Lyapunov analysis.

In C++, the combination of Eigen (for linear algebra) and typical robot dynamics libraries such as RBDL or Orocos KDL allows you to assemble linearized models and perform stability checks in real-time control software (e.g., in a ROS controller).

For a discrete-time system \( x_{k+1} = A_d x_k \), exponential stability is equivalent to all eigenvalues of \( A_d \) lying strictly inside the unit circle. An alternative is to find \( P > 0 \) such that

\[ A_d^{\top} P A_d - P = -Q, \quad Q > 0. \]

Below is a minimal C++ snippet (pseudo-code style) illustrating how to: (i) compute eigenvalues and (ii) check Lyapunov decrease using Eigen.

#include <iostream>
#include <Eigen/Dense>

int main() {
    using namespace Eigen;

    // Example: discrete-time joint error dynamics matrix Ad (4x4)
    Matrix4d Ad;
    Ad <<
        0.98, 0.0,  0.01, 0.0,
        0.0,  0.97, 0.0,  0.01,
       -0.5,  0.0,  0.90, 0.0,
        0.0, -0.4,  0.0,  0.88;

    // 1. Spectral radius check
    EigenSolver<Matrix4d> es(Ad);
    std::cout << "Eigenvalues of Ad:\n" << es.eigenvalues() << std::endl;

    // 2. Lyapunov function candidate V(x) = x^T P x
    Matrix4d Q = Matrix4d::Identity();
    Matrix4d P = Matrix4d::Identity(); // here we just choose P = I for illustration

    // Check that V decreases along one time-step for a sample state x
    Vector4d x;
    x << 0.1, -0.1, 0.0, 0.0;

    double Vx  = x.transpose() * P * x;
    Vector4d x_next = Ad * x;
    double Vxnext = x_next.transpose() * P * x_next;

    std::cout << "V(x) = " << Vx << ", V(x_next) = " << Vxnext << std::endl;

    if (Vxnext < Vx) {
        std::cout << "Lyapunov function decreases: system looks stable for this x." << std::endl;
    } else {
        std::cout << "No decrease for this x: adjust Ad or P." << std::endl;
    }

    return 0;
}
      

In a ROS-based robot controller, Ad could be constructed from a discretized continuous-time closed-loop model using the sampling time of the joint controller. Libraries like RBDL can compute the continuous-time dynamics, while the controller node uses Eigen to test stability properties online or offline.

Java is often used in educational and competition robotics (e.g., FRC's WPILib). We can simulate joint dynamics and verify Lyapunov-like behavior (energy decrease) numerically. For matrix operations, a common choice is EJML.

Below, we simulate the scalar PD-controlled joint: \( I \ddot{q} + k_d \dot{q} + k_p q = 0 \) using explicit Euler integration and inspect the decay of \( V = \tfrac{1}{2} I \dot{q}^2 + \tfrac{1}{2} k_p q^2 \).

public class PdStabilityDemo {
    public static void main(String[] args) {
        double I = 1.0;
        double kp = 50.0;
        double kd = 10.0;

        double q = 0.2;     // initial position error
        double dq = 0.0;    // initial velocity
        double dt = 0.001;  // time step
        double T  = 2.0;    // total simulation time

        double t = 0.0;
        while (t < T) {
            // Dynamics: I * ddq + kd * dq + kp * q = 0
            double ddq = -(kd * dq + kp * q) / I;

            // Integrate
            dq += dt * ddq;
            q  += dt * dq;
            t  += dt;

            // Lyapunov function (energy-like)
            double V = 0.5 * I * dq * dq + 0.5 * kp * q * q;
            if ((int)(t * 1000) % 100 == 0) {
                System.out.println("t=" + t + " q=" + q + " dq=" + dq + " V=" + V);
            }
        }
    }
}
      

In a real Java-based robot project (e.g., using WPILib), the same dynamics can be wrapped into a periodic control loop. The observation that V decreases over time for small step sizes and suitable gains is a discrete-time analog of the continuous-time Lyapunov argument.

9.1 MATLAB/Simulink — Lyapunov Equation and Simulation

MATLAB, with the Control System Toolbox and Robotics System Toolbox, is widely used for robot control design. The continuous-time Lyapunov equation can be solved with lyap.

% Example: 2-DOF linearized joint error dynamics
A = [0 0 1 0;
     0 0 0 1;
    -50 0 -10 0;
     0 -40 0 -12];

% Choose Q > 0
Q = eye(4);

% Solve A' P + P A = -Q
P = lyap(A', Q);

% Check eigenvalues of A (Hurwitz condition)
lambda = eig(A);

disp('Eigenvalues of A:');
disp(lambda);
disp('Lyapunov matrix P:');
disp(P);

% Simulink:
%  - Create a State-Space block with A, B = zeros(4,1), C = eye(4), D = 0.
%  - Feed an initial condition x0 in an IC block.
%  - Scope the states to visualize convergence.
      

9.2 Wolfram Mathematica — Symbolic Lyapunov Analysis

Mathematica supports symbolic and numeric Lyapunov analysis via LyapunovSolve. This is useful when studying parametric robot models.

(* 2-DOF linear joint error dynamics *)
A = { {0, 0, 1, 0},
     {0, 0, 0, 1},
     {-50, 0, -10, 0},
     {0, -40, 0, -12} };

Q = IdentityMatrix[4];

(* Solve continuous-time Lyapunov equation A^T P + P A == -Q *)
P = LyapunovSolve[Transpose[A], -Q];

Eigenvalues[A]
Eigenvalues[P]

(* Define V(x) = x^T P x and check its derivative symbolically *)
x = Array[x, 4];
V = x.P.x;

(* For linear systems, dV/dt = x^T (A^T P + P A) x = -x^T Q x *)
Simplify[Transpose[x].(Transpose[A].P + P.A).x + x.Q.x]
      

In more advanced lessons (e.g., robust and adaptive control), Mathematica can be used to derive symbolic ISS-Lyapunov conditions for simplified robot models.

In this lesson, we introduced several key stability notions for robot control: Lyapunov (internal) stability of equilibria, BIBO stability for input-output behavior, and ISS for robustness to disturbances. We saw how Lyapunov functions translate robot energy or quadratic forms into rigorous stability proofs, and how linear Lyapunov equations provide a computational test for Hurwitz matrices derived from linearized robot dynamics. We also connected BIBO stability to the location of poles of transfer functions and introduced ISS and ISS-Lyapunov functions as a unifying framework for robustness analysis. These tools will be used repeatedly in later chapters (e.g., model-based nonlinear control, robust control, and adaptive control).