Consider an uncertain closed-loop robot system written as

\[ \dot{\mathbf{x}} = f(\mathbf{x}, t, \Delta), \quad \Delta \in \mathcal{D}, \]

where \( \mathbf{x} \) is the state (typically tracking error and its derivative), and \( \Delta \) collects parametric uncertainty, unmodeled dynamics, and disturbances (as introduced in Lesson 1). Robust stability concerns properties that hold for all admissible uncertainties \( \Delta \in \mathcal{D} \).

A family of systems \( \{\Sigma(\Delta)\}_{\Delta\in\mathcal{D}} \) with equilibrium \( \mathbf{x}=0 \) is:

• Robustly stable if for every \( \Delta \in \mathcal{D} \) and every \( \varepsilon > 0 \) there exists \( \delta(\varepsilon) > 0 \) such that \( \|\mathbf{x}(0)\| < \delta(\varepsilon) \) implies \( \|\mathbf{x}(t)\| < \varepsilon \) for all \( t \ge 0 \).
• Robustly asymptotically stable if, in addition, \( \|\mathbf{x}(t)\| \to 0 \) as \( t \to \infty \) for all \( \Delta \in \mathcal{D} \).

For robot manipulators, \( \Delta \) typically encodes:

• Parametric uncertainty in masses, inertias, and link lengths.
• Unmodeled dynamics such as friction, flexibilities, and actuator dynamics.
• Exogenous disturbances, e.g., contact forces or external torques.

Starting from the standard rigid manipulator model (known from the dynamics course)

\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau} + \mathbf{d}(t), \]

robust stability analysis treats the closed loop \( (\mathbf{q},\dot{\mathbf{q}}) \) under a given feedback law \( \boldsymbol{\tau}=\pi(\cdot) \) and admissible uncertainties \( \Delta \).

flowchart TD
  QD["Robot dynamics M(q) qdd + C(q,qd) qd + g(q) = tau + d"] --> Ctl["Controller pi(.) (PD, computed torque, etc.)"]
  Ctl --> CL["Closed-loop error dynamics"]
  CL --> U["Identify uncertainty set: param, unmodeled, disturbance"]
  U --> Tool["Select analysis tool: Lyapunov / ISS / linearization"]
  Tool --> Cond["Derive robust stability inequality or LMI"]
  Cond --> Check["Check gains vs. uncertainty bounds"]
  Check --> Design["If condition satisfied: robustly stable"]
        

Let \( \mathbf{q}_d(t) \) be a smooth desired joint trajectory with position, velocity, and acceleration \( \mathbf{q}_d, \dot{\mathbf{q}}_d, \ddot{\mathbf{q}}_d \). Define tracking errors

\[ \mathbf{e} = \mathbf{q} - \mathbf{q}_d, \quad \dot{\mathbf{e}} = \dot{\mathbf{q}} - \dot{\mathbf{q}}_d. \]

Consider a nominal model-based controller (e.g. a computed-torque or PD+gravity scheme) designed from nominal dynamics \( \hat{\mathbf{M}}, \hat{\mathbf{C}}, \hat{\mathbf{g}} \). With suitable design, the nominal closed-loop error dynamics can be written as

\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{e}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{e}} + \mathbf{K}_d \dot{\mathbf{e}} + \mathbf{K}_p \mathbf{e} = \mathbf{0}, \]

where \( \mathbf{K}_p \), \( \mathbf{K}_d \) are positive definite gain matrices.

In reality, we have uncertainty and disturbances; grouping all unmodeled terms into \( \mathbf{w} \) gives

\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{e}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{e}} + \mathbf{K}_d \dot{\mathbf{e}} + \mathbf{K}_p \mathbf{e} = \mathbf{w}(\mathbf{q},\dot{\mathbf{q}},\mathbf{e},\dot{\mathbf{e}},t). \]

Typical assumptions (consistent with robot dynamics) are:

\[ \begin{aligned} &\mathbf{M}(\mathbf{q}) = \mathbf{M}^\top(\mathbf{q}),\quad m_{\min}\mathbf{I} \preceq \mathbf{M}(\mathbf{q}) \preceq m_{\max}\mathbf{I}, \\ &\dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \text{ is skew-symmetric}, \\ &\|\mathbf{w}(\cdot)\| \le \rho_0 + \rho_1 \|\mathbf{e}\| + \rho_2 \|\dot{\mathbf{e}}\|, \end{aligned} \]

for known nonnegative constants \( \rho_0,\rho_1,\rho_2 \) that bound parametric errors and disturbances over a given operating region.

We now derive a sufficient condition ensuring that the origin \( (\mathbf{e},\dot{\mathbf{e}}) = (\mathbf{0},\mathbf{0}) \) is robustly ultimately bounded for the uncertain error dynamics. Consider the Lyapunov candidate

\[ V(\mathbf{e},\dot{\mathbf{e}},\mathbf{q}) = \tfrac{1}{2}\dot{\mathbf{e}}^\top \mathbf{M}(\mathbf{q})\dot{\mathbf{e}} + \tfrac{1}{2}\mathbf{e}^\top \mathbf{K}_p \mathbf{e}. \]

Using the bounds on \( \mathbf{M} \) and \( \mathbf{K}_p \), we have

\[ \tfrac{1}{2}m_{\min}\|\dot{\mathbf{e}}\|^2 \le V(\mathbf{e},\dot{\mathbf{e}},\mathbf{q}) \le \tfrac{1}{2}m_{\max}\|\dot{\mathbf{e}}\|^2 + \tfrac{1}{2}\lambda_{\max}(\mathbf{K}_p)\|\mathbf{e}\|^2. \]

Differentiating \( V \) along trajectories and using the skew-symmetry property of \( \dot{\mathbf{M}}-2\mathbf{C} \) gives the standard robot identity

\[ \dot{V} = -\dot{\mathbf{e}}^\top \mathbf{K}_d \dot{\mathbf{e}} + \dot{\mathbf{e}}^\top \mathbf{w}. \]

Applying Cauchy–Schwarz and the bound on \( \mathbf{w} \):

\[ \dot{V} \le -\lambda_{\min}(\mathbf{K}_d)\|\dot{\mathbf{e}}\|^2 + \|\dot{\mathbf{e}}\|\, \bigl(\rho_0 + \rho_1\|\mathbf{e}\| + \rho_2\|\dot{\mathbf{e}}\|\bigr). \]

Rearranging,

\[ \dot{V} \le -\bigl(\lambda_{\min}(\mathbf{K}_d)-\rho_2\bigr)\|\dot{\mathbf{e}}\|^2 + \rho_1\|\mathbf{e}\|\|\dot{\mathbf{e}}\| + \rho_0\|\dot{\mathbf{e}}\|. \]

Using the inequality \( ab \le \tfrac{\varepsilon}{2}a^2 + \tfrac{1}{2\varepsilon}b^2 \) for any \( \varepsilon > 0 \) on the mixed terms, we obtain

\[ \begin{aligned} \rho_1\|\mathbf{e}\|\|\dot{\mathbf{e}}\| &\le \tfrac{\varepsilon_1}{2}\|\dot{\mathbf{e}}\|^2 + \tfrac{\rho_1^2}{2\varepsilon_1}\|\mathbf{e}\|^2, \\ \rho_0\|\dot{\mathbf{e}}\| &\le \tfrac{\varepsilon_2}{2}\|\dot{\mathbf{e}}\|^2 + \tfrac{\rho_0^2}{2\varepsilon_2}. \end{aligned} \]

Choosing small \( \varepsilon_1,\varepsilon_2 > 0 \) and substituting,

\[ \dot{V} \le -\alpha_d\|\dot{\mathbf{e}}\|^2 + \alpha_e\|\mathbf{e}\|^2 + c_0, \]

where

\[ \alpha_d = \lambda_{\min}(\mathbf{K}_d) -\rho_2-\tfrac{\varepsilon_1}{2}-\tfrac{\varepsilon_2}{2}, \quad \alpha_e = \tfrac{\rho_1^2}{2\varepsilon_1}, \quad c_0 = \tfrac{\rho_0^2}{2\varepsilon_2}. \]

Provided we choose gains such that \( \alpha_d > 0 \) (i.e. \( \lambda_{\min}(\mathbf{K}_d) \) sufficiently large compared to \( \rho_2 \)), the derivative is negative definite outside a ball around the origin. Using norm equivalences between \( V \) and \( \|(\mathbf{e},\dot{\mathbf{e}})\|^2 \) we obtain:

Theorem 1 (Lyapunov Robust Stability Condition).

Assume the robot dynamics satisfy the standard manipulator properties, and that the disturbance/uncertainty term \( \mathbf{w} \) obeys the affine bound above. If the derivative gains are chosen such that

\[ \lambda_{\min}(\mathbf{K}_d) > \rho_2, \]

then there exist positive constants \( c_1,c_2,c_3 \) such that for all admissible uncertainties,

\[ \dot{V} \le -c_1 \bigl(\|\mathbf{e}\|^2+\|\dot{\mathbf{e}}\|^2\bigr) + c_2, \]

and the tracking error is uniformly ultimately bounded in a ball of radius \( \sqrt{c_2/c_1} \). Increasing \( \mathbf{K}_d \) shrinks this bound.

This condition is conservative (it uses bounds) but extremely common in robust robot control: it formalizes the intuition that a sufficiently high damping gain \( \mathbf{K}_d \) overcomes bounded modeling errors and disturbances.

From Lesson 1 in Chapter 1, recall that an input-to-state stable (ISS) system satisfies

\[ \|\mathbf{x}(t)\| \le \beta\bigl(\|\mathbf{x}(0)\|, t\bigr) + \gamma\Bigl(\sup_{0 \le s \le t} \|\mathbf{w}(s)\|\Bigr), \]

for some class-\(\mathcal{K}\mathcal{L}\) function \( \beta \) and class-\(\mathcal{K}\) function \( \gamma \). The Lyapunov inequality from the previous section has exactly the ISS form: the derivative is negative apart from a term depending on the disturbance magnitude.

More precisely, using the comparison lemma and norm equivalences between \( V \) and \( \|(\mathbf{e},\dot{\mathbf{e}})\|^2 \), one can show that

\[ \|(\mathbf{e}(t),\dot{\mathbf{e}}(t))\| \le \beta\bigl(\|(\mathbf{e}(0),\dot{\mathbf{e}}(0))\|, t\bigr) + \gamma(\rho_0), \]

where \( \rho_0 \) is the bound on the disturbance offset. Thus the robot under PD/model-based control is ISS with respect to disturbances \( \mathbf{w} \), and its steady-state tracking error scales monotonically with the disturbance bound.

This ISS viewpoint is extremely useful in later chapters (e.g. disturbance observers, learning augmentation), where the disturbance term \( \mathbf{w} \) is itself the output of another dynamical system that must be kept “small” in an input–output sense.

flowchart TD
  E["Error state x = (e, edot)"] --> Sys["Closed-loop robot (nominal + PD)"]
  Sys --> W["Uncertainty / disturbance w"]
  W --> Sys
  Sys --> ISS["Lyapunov inequality: Vdot <= -alpha(||x||) + sigma(||w||)"]
  ISS --> Bound["ISS bound: ||x(t)|| <= beta(||x0||,t) + gamma(sup||w||)"]
        

Around a fixed trajectory or operating point, the joint-space error dynamics can be linearized as

\[ \dot{\mathbf{x}} = \bigl(\mathbf{A} + \Delta\mathbf{A}\bigr)\mathbf{x}, \quad \|\Delta\mathbf{A}\| \le \delta, \]

where \( \mathbf{x} \) stacks \( \mathbf{e} \) and \( \dot{\mathbf{e}} \), and \( \delta \) is a known bound on the linearization error plus parametric mismatch.

For the nominal system \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \), exponential stability is equivalent to the existence of a quadratic Lyapunov function \( V=\mathbf{x}^\top \mathbf{P}\mathbf{x} \) with \( \mathbf{P}\succ 0 \) such that

\[ \mathbf{A}^\top\mathbf{P} + \mathbf{P}\mathbf{A} \preceq -\mathbf{Q}, \quad \mathbf{Q}\succ 0. \]

For the uncertain system, we have

\[ \dot{V} = \mathbf{x}^\top\bigl((\mathbf{A}+\Delta\mathbf{A})^\top\mathbf{P} + \mathbf{P}(\mathbf{A}+\Delta\mathbf{A})\bigr)\mathbf{x}. \]

Using the bound \( \|\Delta\mathbf{A}\| \le \delta \Rightarrow \|\mathbf{P}\Delta\mathbf{A} + \Delta\mathbf{A}^\top\mathbf{P}\| \le 2\delta\|\mathbf{P}\| \), we obtain

\[ \dot{V} \le \mathbf{x}^\top \bigl(\mathbf{A}^\top\mathbf{P} + \mathbf{P}\mathbf{A} + 2\delta\|\mathbf{P}\|\mathbf{I}\bigr) \mathbf{x}. \]

Therefore, the following LMI is a sufficient robust stability condition:

Proposition 2 (Quadratic Robust Stability Condition).

If there exists \( \mathbf{P}\succ 0 \) and \( \mathbf{Q}\succ 0 \) such that

\[ \mathbf{A}^\top\mathbf{P} + \mathbf{P}\mathbf{A} + 2\delta\|\mathbf{P}\|\mathbf{I} \preceq -\mathbf{Q}, \]

then \( \dot{\mathbf{x}}=(\mathbf{A}+\Delta\mathbf{A})\mathbf{x} \) is exponentially stable for all \( \Delta\mathbf{A} \) with \( \|\Delta\mathbf{A}\| \le \delta \).

In practice, one solves a slightly more conservative but convex LMI by bounding \( \|\mathbf{P}\| \) with an auxiliary variable, or uses tools such as MATLAB's Robust Control Toolbox to numerically verify this condition.

Consider a single revolute joint with uncertain inertia \( J \) and viscous friction \( b \):

\[ J\ddot{q} + b\dot{q} = \tau. \]

Let the desired trajectory be constant \( q_d \) (regulation). A simple PD controller is

\[ \tau = -k_p(q-q_d) - k_d(\dot{q}-0), \quad k_p > 0,\; k_d > 0. \]

Define \( e=q-q_d \), \( \dot{e}=\dot{q} \). The closed-loop dynamics are

\[ J\ddot{e} + (b+k_d)\dot{e} + k_p e = 0. \]

Dividing by \( J \) gives

\[ \ddot{e} + a_1\dot{e} + a_0 e = 0, \quad a_1 = \frac{b+k_d}{J},\; a_0 = \frac{k_p}{J}. \]

The characteristic polynomial is \( s^2 + a_1 s + a_0 \). From basic control theory (Routh–Hurwitz for a second-order system), all roots are in the open left half-plane if and only if \( a_1 > 0, a_0 > 0 \). Since \( J > 0 \), \( k_p > 0 \), and \( b+k_d > 0 \), we conclude:

Result. If \( k_p > 0 \) and \( b + k_d > 0 \), the single-joint PD controller is robustly asymptotically stable with respect to any positive inertia \( J \), no matter how large or small (but positive) it is. Thus an entire interval \( J \in [J_{\min},J_{\max}] \subset (0,\infty) \) is covered without changing the controller.

This simple example illustrates robust stability via parameter-insensitive sign conditions on the characteristic polynomial, a recurring theme in robust robot control design.

We now implement a grid-based robust stability check and a simple simulation for the single-joint PD system using Python. For more complex manipulators one would typically combine this with libraries such as roboticstoolbox-python, but here we stay with a scalar joint to keep formulas explicit.

import numpy as np
from numpy.linalg import eigvals
from math import isfinite

# Nominal parameters
J_nom = 1.0     # nominal inertia
b = 0.2         # viscous friction
k_p = 25.0
k_d = 10.0

def A_matrix(J, b, k_p, k_d):
    """
    State matrix for single-joint PD control:
        e_dot   = v
        v_dot   = -(k_p/J) * e - ((b + k_d)/J) * v
    """
    return np.array([[0.0,        1.0],
                     [-k_p / J,  -(b + k_d) / J]])

# Uncertainty range on inertia: J in [0.5, 1.5]
J_grid = np.linspace(0.5, 1.5, 21)
robust_stable = True

for J in J_grid:
    A = A_matrix(J, b, k_p, k_d)
    lam = eigvals(A)
    max_real = np.max(np.real(lam))
    if max_real >= 0.0 or not isfinite(max_real):
        print("Potential instability at J =", J, "eigs =", lam)
        robust_stable = False
        break

print("Grid-based robustness check (eigenvalues real parts < 0):", robust_stable)

# Simple time-domain simulation for several J values
def simulate_joint(J, b, k_p, k_d, e0=0.5, edot0=0.0, dt=1e-3, T=2.0):
    n_steps = int(T / dt)
    e = e0
    edot = edot0
    traj = []
    for k in range(n_steps):
        # PD torque for regulation to e = 0
        tau = -k_p * e - k_d * edot
        # Dynamics: J * eddot + b * edot = tau
        eddot = (tau - b * edot) / J
        # Euler integration (small dt)
        edot = edot + dt * eddot
        e = e + dt * edot
        traj.append((k * dt, e, edot))
    return np.array(traj)

if robust_stable:
    for J in [0.5, 1.0, 1.5]:
        traj = simulate_joint(J, b, k_p, k_d)
        print("Final error for J =", J, "is", traj[-1, 1])
      

The eigenvalue check corresponds directly to verifying that the linearization is asymptotically stable for each admissible inertia value. For nonlinear multi-DOF robots, the same idea can be applied to local linearizations along a trajectory.

In C++, linear algebra in robotics is frequently handled by the Eigen library. The following snippet performs the same grid-based eigenvalue robustness check as the Python code.

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

int main() {
    double b = 0.2;
    double k_p = 25.0;
    double k_d = 10.0;

    auto A_matrix = [&] (double J) {
        Eigen::Matrix2d A;
        A << 0.0,        1.0,
              -k_p / J,  -(b + k_d) / J;
        return A;
    };

    std::vector<double> J_values;
    for (int i = 0; i <= 20; ++i) {
        double J = 0.5 + (1.0 / 20.0) * i; // from 0.5 to 1.5
        J_values.push_back(J);
    }

    bool robust_stable = true;
    for (double J : J_values) {
        Eigen::Matrix2d A = A_matrix(J);
        Eigen::EigenSolver<Eigen::Matrix2d> es(A);
        Eigen::VectorXcd lam = es.eigenvalues();
        double max_real = lam.real().maxCoeff();
        if (max_real >= 0.0) {
            std::cout << "Unstable or marginally stable at J = "
                      << J << std::endl;
            std::cout << "Eigenvalues: " << lam.transpose() << std::endl;
            robust_stable = false;
            break;
        }
    }

    std::cout << "Robust stability on grid: "
              << (robust_stable ? "true" : "false") << std::endl;
    return 0;
}
      

For higher-DOF manipulators one would assemble block matrices representing joint-space linearizations and re-use the same eigenvalue logic; Eigen naturally scales to those dimensions.

In Java-based robotics stacks, one often uses matrix libraries such as EJML for linear algebra. Here we focus on a simple explicit simulation loop that checks decay of the error for several inertia values.

public class RobustSingleJointPD {

    static class TrajectorySample {
        public final double t;
        public final double e;
        public final double edot;

        public TrajectorySample(double t, double e, double edot) {
            this.t = t;
            this.e = e;
            this.edot = edot;
        }
    }

    public static TrajectorySample[] simulate(double J,
                                             double b,
                                             double k_p,
                                             double k_d,
                                             double e0,
                                             double edot0,
                                             double dt,
                                             double T) {
        int nSteps = (int) Math.round(T / dt);
        TrajectorySample[] traj = new TrajectorySample[nSteps];
        double e = e0;
        double edot = edot0;
        for (int k = 0; k < nSteps; ++k) {
            double tau = -k_p * e - k_d * edot;
            double eddot = (tau - b * edot) / J;
            edot = edot + dt * eddot;
            e = e + dt * edot;
            double t = k * dt;
            traj[k] = new TrajectorySample(t, e, edot);
        }
        return traj;
    }

    public static void main(String[] args) {
        double b = 0.2;
        double k_p = 25.0;
        double k_d = 10.0;

        double[] J_values = {0.5, 1.0, 1.5};
        for (double J : J_values) {
            TrajectorySample[] traj =
                simulate(J, b, k_p, k_d,
                         0.5, 0.0, 1e-3, 2.0);
            TrajectorySample last = traj[traj.length - 1];
            System.out.println("J = " + J
                + ", final e = " + last.e
                + ", final edot = " + last.edot);
        }
    }
}
      

Observing that e and edot converge to zero for the tested inertia values provides numerical evidence of robust stability consistent with the analytical conditions.

MATLAB is widely used in robotics for both analysis and implementation (with Simulink for real-time or HIL setups). Below is a script that sweeps the inertia parameter, computes the closed-loop eigenvalues, and can be combined with a Simulink model of the joint dynamics.

b   = 0.2;
k_p = 25.0;
k_d = 10.0;

J_values = linspace(0.5, 1.5, 21);
robustStable = true;

for J = J_values
    A = [0,              1;
        -k_p / J,  -(b + k_d) / J];
    lam = eig(A);
    if max(real(lam)) >= 0
        fprintf('Potential instability at J = %.3f\n', J);
        disp(lam);
        robustStable = false;
        break;
    end
end

fprintf('Robust stability on grid: %d\n', robustStable);

% Optional: connect to a Simulink model "single_joint_pd"
% where J is a workspace parameter used in the dynamics block.
% set_param('single_joint_pd', 'SimulationCommand', 'start');
      

In Simulink, the joint can be represented by a second-order transfer function block parameterized by J and b, with a PD controller implemented using basic gain and summation blocks. One can then batch-run simulations for different values of J and visually confirm robust convergence.

Wolfram Mathematica provides symbolic and numeric tools to verify stability conditions, including symbolic eigenvalues and Lyapunov inequalities.

(* Parameters *)
b   = 0.2;
kP  = 25.0;
kD  = 10.0;

(* State matrix for given J *)
A[J_] := { {0, 1},
          {-kP/J, -(b + kD)/J} };

(* Robust eigenvalue sweep *)
jVals = Range[0.5, 1.5, 0.05];
robustStable = True;

Do[
  lam = Eigenvalues[A[J]];
  maxReal = Max[Re[lam]];
  If[maxReal >= 0,
    Print["Unstable or marginal at J = ", J, " eigenvalues = ", lam];
    robustStable = False;
    Break[];
  ],
  {J, jVals}
];

Print["Robust stability on grid: ", robustStable];

(* Example Lyapunov matrix for a chosen J0 *)
J0 = 1.0;
A0 = A[J0];

(* Solve continuous Lyapunov equation A0' P + P A0 = -Q *)
Q = IdentityMatrix[2];
P = LyapunovSolve[Transpose[A0], -Q];
Print["Lyapunov matrix P at J0 = 1: ", MatrixForm[P]];
      

The Lyapunov equation solution P demonstrates the existence of a quadratic Lyapunov function for the nominal system; combined with the perturbation bounds, this leads back to the matrix inequality conditions in Section 5.

In this lesson we formalized robust stability for uncertain robot manipulators, derived Lyapunov and ISS-based sufficient conditions, and connected them to linearized matrix inequalities. For single-joint PD control, we showed explicitly how Routh–Hurwitz conditions yield parametric robust stability with respect to inertia variation. We then implemented practical robustness checks and simulations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. These robust stability conditions underpin the design of sliding-mode controllers, disturbance observers, and adaptive laws in the subsequent lessons of this chapter.