Let the robot joint configuration be \( q \in \mathbb{R}^n \), with dynamics already known from kinematics/dynamics:

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

where \( M(q) \) is positive definite, \( C(q,\dot{q}) \) is the Coriolis/centrifugal term, \( g(q) \) is gravity, \( \tau \) is the control torque, and \( d(t) \) aggregates disturbances/unmodeled dynamics.

Define the tracking error for a desired joint trajectory \( q_d(t) \) as \( e = q - q_d \), \( \dot{e} = \dot{q} - \dot{q}_d \). Your controller (possibly including PD, computed-torque, impedance, MPC, etc.) produces \( \tau = u(q,\dot{q},q_d,\dot{q}_d,\ldots) \).

Introduce the state vector \( x = \begin{bmatrix} e \\ \dot{e} \end{bmatrix} \in \mathbb{R}^{2n} \). The closed-loop error dynamics can be written abstractly as

\[ \dot{x} = f(x,t) + G(x,t)\,d(t). \]

For local analysis around an equilibrium or a nominal trajectory \( x^{\star}(t) \), we use a linearization of the form

\[ \dot{\tilde{x}} = A_{\text{cl}}(t)\,\tilde{x} + B_d(t)\,\tilde{d}(t), \quad \tilde{x} = x - x^{\star},\; \tilde{d} = d - d^{\star}. \]

For time-invariant controllers around a fixed equilibrium, \( A_{\text{cl}} \) is constant. Internal (asymptotic) stability of the nominal closed loop is characterized by \( A_{\text{cl}} \) being Hurwitz, i.e. all eigenvalues having strictly negative real parts.

For a linear time-invariant (LTI) closed-loop model \( \dot{x} = A_{\text{cl}} x \), the classical Lyapunov theorem says:

\( A_{\text{cl}} \) is Hurwitz if and only if for every symmetric positive definite matrix \( Q \), there exists a unique symmetric positive definite \( P \) solving

\[ A_{\text{cl}}^{\top} P + P A_{\text{cl}} = -Q. \]

A quadratic Lyapunov function \( V(x) = x^{\top} P x \) then satisfies

\[ \dot{V}(x) = x^{\top} \left( A_{\text{cl}}^{\top} P + P A_{\text{cl}} \right) x = - x^{\top} Q x \leq 0, \]

and in fact \( \dot{V}(x) \leq -\lambda_{\min}(Q)\,\|x\|^2 \), which implies exponential convergence of \( x(t) \) to zero.

An important numerical stability margin is the spectral abscissa \( \alpha(A_{\text{cl}}) = \max_i \operatorname{Re}(\lambda_i(A_{\text{cl}})) \). A negative value indicates stability, and the magnitude \( -\alpha(A_{\text{cl}}) \) measures how far the eigenvalues lie from the imaginary axis. For a complex eigenvalue \( \lambda = \sigma + j \omega \), the approximate damping ratio is

\[ \zeta = -\frac{\sigma}{\sqrt{\sigma^2 + \omega^2}}, \]

which is directly related to overshoot for lightly damped modes. In your capstone project, you should compute \( \alpha(A_{\text{cl}}) \) and the damping ratios of the dominant eigenvalues for each operating point of interest.

For nonlinear robot controllers (e.g., computed-torque, sliding mode, adaptive), Lyapunov analysis proceeds at the level of the original dynamics. For instance, consider a fully actuated manipulator with a simple PD + gravity compensation law (ignoring Coriolis terms for illustration):

\[ M(q)\ddot{e} + K_d \dot{e} + K_p e = 0,\quad K_p = K_p^{\top} > 0,\; K_d = K_d^{\top} > 0. \]

Choose \( V(e,\dot{e}) = \tfrac{1}{2}\dot{e}^{\top} M(q) \dot{e} + \tfrac{1}{2} e^{\top} K_p e \). Using symmetry and positive definiteness of \( M(q) \) and \( K_p \), we have

\[ V(e,\dot{e}) \geq \tfrac{1}{2}\lambda_{\min}(M)\|\dot{e}\|^2 + \tfrac{1}{2}\lambda_{\min}(K_p)\|e\|^2, \]

so \( V \) is radially unbounded in \((e,\dot{e})\). Differentiating along trajectories and using the skew-symmetry property of \( \dot{M}(q) - 2C(q,\dot{q}) \), one obtains

\[ \dot{V}(e,\dot{e}) = -\dot{e}^{\top} K_d \dot{e} \leq -\lambda_{\min}(K_d)\,\|\dot{e}\|^2 \leq 0, \]

implying global asymptotic stability of the tracking error dynamics for the nominal model. This Lyapunov certificate remains a central tool for your capstone analysis.

Suppose your nominal linearized closed loop is \( \dot{x} = A_{\text{cl}} x \), and the true dynamics are

\[ \dot{x} = (A_{\text{cl}} + \Delta A) x, \]

where \( \Delta A \) captures parametric uncertainty (e.g., payload variation, friction errors). Let \( P > 0 \) solve the Lyapunov equation

\[ A_{\text{cl}}^{\top} P + P A_{\text{cl}} = -Q,\quad Q = Q^{\top} > 0. \]

Consider again \( V(x) = x^{\top} P x \). For the perturbed system,

\[ \begin{aligned} \dot{V}(x) &= x^{\top}\left[(A_{\text{cl}}+\Delta A)^{\top}P + P(A_{\text{cl}}+\Delta A)\right]x \\ &= -x^{\top} Q x + x^{\top}\left(\Delta A^{\top}P + P\Delta A\right)x. \end{aligned} \]

Using the induced 2-norm, \( \|x^{\top}\Delta A^{\top}P x\| \leq \|P\|_2\,\|\Delta A\|_2\,\|x\|^2 \) and similarly for \( P\Delta A \), so

\[ \dot{V}(x) \leq -\lambda_{\min}(Q)\,\|x\|^2 + 2\|P\|_2\,\|\Delta A\|_2\,\|x\|^2. \]

Therefore, a sufficient condition for robust exponential stability of the perturbed system for all uncertainties satisfying \( \|\Delta A\|_2 \leq \bar{\delta} \) is

\[ 2\|P\|_2 \,\bar{\delta} < \lambda_{\min}(Q). \]

Equivalently, \( \bar{\delta} < \dfrac{\lambda_{\min}(Q)}{2\|P\|_2} \). In practice, you:

• solve the Lyapunov equation for \( P \) given a convenient \( Q \) (typically identity or diagonal);
• compute \( \lambda_{\min}(Q) \) and \( \|P\|_2 \);
• obtain a conservative bound on admissible model error \( \bar{\delta} \).

This bound can be compared against estimated modeling errors in inertia, payload, and friction obtained earlier in your project.

Many disturbances enter the robot as inputs: joint torque disturbances, force interaction, or sensor noise injected at the controller input. For an LTI disturbance channel \( \dot{x} = A_{\text{cl}}x + B_d d,\; y = Cx \), the induced \( \mathcal{L}_2 \) gain defines a robustness measure:

\[ \|G\|_{\infty} = \sup_{\omega \in \mathbb{R}} \bar{\sigma}\big(G(j\omega)\big),\quad G(s) = C(sI - A_{\text{cl}})^{-1}B_d, \]

where \( \bar{\sigma}(\cdot) \) is the largest singular value. A small \( \|G\|_{\infty} \) indicates good attenuation of worst-case disturbances across all frequencies.

For nonlinear closed loops, input-to-state stability (ISS) formalizes robustness: roughly, the state remains small whenever disturbances are small. An ISS-Lyapunov function \( V(x) \) satisfies, for some class-\( \mathcal{K}_{\infty} \) functions \( \alpha_1,\alpha_2,\alpha_3 \) and class-\( \mathcal{K} \) function \( \gamma \),

\[ \alpha_1(\|x\|) \leq V(x) \leq \alpha_2(\|x\|),\quad \dot{V}(x,d) \leq -\alpha_3(\|x\|) + \gamma(\|d\|). \]

Your previous lessons on robust and adaptive robot control have effectively constructed such \( V \) for specific architectures (e.g., sliding-mode or adaptive computed-torque). In this capstone lesson, you interpret simulation results through this ISS lens: you verify that disturbance energy is “channeled” into bounded state deviations.

For SISO loops, classical gain and phase margins (computed from the loop transfer function) provide intuitive robustness measures, and for MIMO loops, singular-value Bode plots of \( G(j\omega) \) provide analogous information.

A systematic stability/robustness evaluation for your capstone controller can follow the workflow below.

flowchart TD
  M["Closed loop model from project"] --> L["Local linearization A_cl,B_d"]
  L --> EIG["Eigenvalues,damping ratios"]
  L --> LYA["Solve Lyapunov eq A_cl^T P + P A_cl = -Q"]
  EIG --> DEC["Check stability margins and specs"]
  LYA --> ROB["Compute delta_bar from P,Q"]
  DEC --> SIM["Nonlinear simulations with disturbances"]
  ROB --> SIM
  SIM --> PASS["Specs satisfied?"]
  PASS -->|yes| DOC["Document results in report"]
  PASS -->|no| RETUNE["Retune controller / redesign"]
  RETUNE --> M
        

This loop is executed for all relevant operating conditions: free motion, contact, different payloads, and different trajectories.

Analytical bounds based on Lyapunov functions are conservative. To complement them, you should perform Monte Carlo testing over uncertainty and disturbance samples. A typical Monte Carlo campaign consists of:

• Sampling uncertain physical parameters (inertias, mass, friction) within admissible intervals or probability distributions.
• Sampling reference trajectories and contact conditions.
• Injecting disturbances (step torques, sinusoids, noise) and measuring performance metrics (tracking error, control effort, constraint violations).

flowchart TD
  U["Sample uncertainty (mass,friction,delay)"] --> S["Simulate closed loop over horizon"]
  D["Sample disturbances,noise"] --> S
  S --> MET["Compute metrics: max error,overshoot,safety"]
  MET --> THR["Compare with thresholds"]
  THR -->|violated| LOG["Log failure scenario"]
  THR -->|ok| CNT["Count as success"]
        

From these simulations you can estimate empirical probabilities of failure, worst-case performance, and robustness margins that supplement your analytic guarantees.

In Python, you can combine numpy, scipy, the python-control library, and a robotics library such as roboticstoolbox to evaluate stability and robustness. Below is a simplified example for one joint with PD control around an equilibrium, plus a Monte Carlo test over inertia uncertainty.

import numpy as np
from numpy.linalg import eig, norm
from scipy.linalg import solve_continuous_lyapunov as lyap
import control

# Nominal scalar joint model: J * q_dd + Kd * q_d + Kp * q = 0
J_nom = 0.5   # nominal inertia
Kp = 50.0
Kd = 5.0

def A_cl(J):
    return np.array([[0.0, 1.0],
                     [-Kp / J, -Kd / J]])

# Nominal closed-loop matrix and eigenvalues
A = A_cl(J_nom)
evals, _ = eig(A)
print("Nominal eigenvalues:", evals)

# Lyapunov-based robustness bound
Q = np.eye(2)
P = lyap(A.T, Q)  # solves A^T P + P A = -Q

lam_min_Q = np.min(np.linalg.eigvals(Q).real)
P_norm = norm(P, 2)
delta_bar = lam_min_Q / (2.0 * P_norm)
print("Conservative bound on ||Delta A||_2:", delta_bar)

# Monte Carlo: sample inertia uncertainty and test stability numerically
rng = np.random.default_rng(0)
num_samples = 200
unstable_count = 0
worst_alpha = -1e9

for k in range(num_samples):
    # Sample inertia within +-30%
    J = J_nom * (1.0 + 0.3 * (2.0 * rng.random() - 1.0))
    A_k = A_cl(J)
    evals_k, _ = eig(A_k)
    alpha_k = np.max(evals_k.real)
    worst_alpha = max(worst_alpha, alpha_k)
    if alpha_k >= 0.0:
        unstable_count += 1

print("Worst spectral abscissa over samples:", worst_alpha)
print("Number of unstable samples:", unstable_count)

# Input-output robustness: L2 gain of disturbance channel
# Disturbance torque enters acceleration: x_dot = A x + B_d * d
B_d = np.array([[0.0],
                [1.0 / J_nom]])
C = np.eye(2)   # measure states for simplicity
G = control.ss(A, B_d, C, 0.0)
hinf_norm, peak_w = control.hinfnorm(G)
print("H-infinity norm of disturbance-to-state map:", hinf_norm)
print("Frequency at which peak occurs:", peak_w)
      

In your actual project, A and B_d will come from linearizing the full manipulator model (for instance, with roboticstoolbox) around your reference trajectories and plugging your actual control law into the dynamics.

For real-time implementations (e.g., using ROS control), C++ with the Eigen library is standard. Below is a minimal example that computes eigenvalues of the closed-loop matrix for a scalar joint PD controller and prints the spectral abscissa.

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

int main() {
  double J_nom = 0.5;
  double Kp = 50.0;
  double Kd = 5.0;

  Eigen::Matrix2d A;
  A << 0.0, 1.0,
        -Kp / J_nom, -Kd / J_nom;

  Eigen::EigenSolver<Eigen::Matrix2d> es(A);
  Eigen::VectorXcd evals = es.eigenvalues();

  double alpha = -1e9;
  for (int i = 0; i < evals.size(); ++i) {
    double real_part = evals(i).real();
    alpha = std::max(alpha, real_part);
  }

  std::cout << "Eigenvalues: " << evals.transpose() << std::endl;
  std::cout << "Spectral abscissa alpha(A_cl): " << alpha << std::endl;

  if (alpha < 0.0) {
    std::cout << "Closed loop is (locally) asymptotically stable." << std::endl;
  } else {
    std::cout << "Closed loop is unstable or marginally stable." << std::endl;
  }

  return 0;
}
      

In a ROS-based project, you would embed similar checks in offline tools or unit tests that operate on the linearizations of your controller plus the robot model provided by, for example, urdf parsers and dynamics libraries (e.g., KDL, pinocchio).

For Java-based robotics stacks (e.g., using ROSJava), numerical linear algebra can be handled through libraries like EJML or Apache Commons Math. The following example uses Apache Commons Math to compute eigenvalues and a simple Monte Carlo stability check.

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import java.util.Random;

public class JointStabilityCheck {
    private static RealMatrix Acl(double J, double Kp, double Kd) {
        double[][] data = {
            {0.0, 1.0},
            {-Kp / J, -Kd / J}
        };
        return new Array2DRowRealMatrix(data);
    }

    public static void main(String[] args) {
        double J_nom = 0.5;
        double Kp = 50.0;
        double Kd = 5.0;

        RealMatrix A = Acl(J_nom, Kp, Kd);
        EigenDecomposition ed = new EigenDecomposition(A);
        double[] realParts = ed.getRealEigenvalues();
        double alpha = -1e9;
        for (double r : realParts) {
            if (r > alpha) alpha = r;
        }
        System.out.println("Nominal spectral abscissa: " + alpha);

        // Monte Carlo over inertia uncertainty
        Random rng = new Random(0);
        int numSamples = 200;
        int unstable = 0;
        double worstAlpha = -1e9;

        for (int k = 0; k < numSamples; ++k) {
            double delta = 0.3 * (2.0 * rng.nextDouble() - 1.0); // +-30%
            double J = J_nom * (1.0 + delta);
            RealMatrix Ak = Acl(J, Kp, Kd);
            EigenDecomposition edk = new EigenDecomposition(Ak);
            double[] rparts = edk.getRealEigenvalues();
            double alpha_k = -1e9;
            for (double r : rparts) {
                if (r > alpha_k) alpha_k = r;
            }
            worstAlpha = Math.max(worstAlpha, alpha_k);
            if (alpha_k >= 0.0) unstable++;
        }

        System.out.println("Worst spectral abscissa over samples: " + worstAlpha);
        System.out.println("Unstable samples: " + unstable);
    }
}
      

Such utilities integrate naturally with Java-based simulation and visualization environments used for your capstone demonstrations.

MATLAB with the Control System Toolbox and Robotics System Toolbox is particularly convenient for stability and robustness analysis, and for integrating with Simulink models of your capstone controller.

% Joint PD closed-loop example
J_nom = 0.5;
Kp = 50;
Kd = 5;

A = [0 1;
    -Kp/J_nom -Kd/J_nom];
B_d = [0; 1/J_nom];
C = eye(2);
D = [0; 0];

% Eigenvalues and damping
ev = eig(A);
alpha = max(real(ev));

fprintf('Eigenvalues:\n');
disp(ev);
fprintf('Spectral abscissa alpha(A_cl) = %g\n', alpha);

% Lyapunov equation
Q = eye(2);
P = lyap(A', Q);
lam_min_Q = min(eig(Q));
P_norm = norm(P, 2);
delta_bar = lam_min_Q / (2 * P_norm);
fprintf('Conservative bound on ||Delta A||_2: %g\n', delta_bar);

% Input-output robustness (H-infinity norm)
G = ss(A, B_d, C, D);
[hinfNorm, freq] = hinfnorm(G);
fprintf('H-infinity norm: %g at frequency %g rad/s\n', hinfNorm, freq);

% Simulink integration:
% 1. Build or reuse a Simulink model of your robot + controller.
% 2. Use linearize() or slLinearizer to obtain (A,B,C,D) at operating points.
% 3. Call eig, lyap, hinfnorm, and margin on these linearizations.
      

For full robot models, rigidBodyTree from Robotics System Toolbox can provide dynamics, which can be linearized with your controller block using linearize at specified configurations and velocities.

In Wolfram Mathematica, you can use StateSpaceModel, LyapunovSolve, and frequency-domain tools to analyze the same closed-loop systems.

(* Joint PD closed loop *)
Jnom = 0.5;
Kp = 50.;
Kd = 5.;

A = { {0., 1.}, {-Kp/Jnom, -Kd/Jnom} };
Bd = { {0.}, {1./Jnom} };
C = IdentityMatrix[2];
D = { {0.}, {0.} };

sys = StateSpaceModel[{A, Bd, C, D}];

(* Eigenvalues and damping ratios *)
evals = Eigenvalues[A];
alpha = Max[Re[evals]];
Print["Eigenvalues: ", evals];
Print["Spectral abscissa: ", alpha];

(* Lyapunov equation A^T P + P A = -Q *)
Q = IdentityMatrix[2];
P = LyapunovSolve[Transpose[A], -Q];
lambdaMinQ = Min[Eigenvalues[Q]];
pNorm = Norm[P, 2];
deltaBar = lambdaMinQ/(2. pNorm);
Print["Conservative bound on ||Delta A||_2: ", deltaBar];

(* H-infinity norm of disturbance-to-state map *)
hinf = HInfinityNorm[sys];
Print["H-infinity norm: ", hinf];

(* Nyquist or singular value plots for full-loop analysis can be done with:
   NyquistPlot[sys];  or  SingularValuePlot[sys]; *)
      

Mathematica is especially handy for symbolic manipulations of Lyapunov conditions and for verifying stability of simplified low-DOF models before moving to high-DOF numerical analysis.

Problem 1 (Lyapunov Stability of PD-Controlled Joint):
Consider a single-DOF joint with inertia \( J > 0 \) and dynamics \( J \ddot{q} = \tau \). The control law is \( \tau = -K_p (q - q_d) - K_d (\dot{q} - \dot{q}_d) \) with \( K_p > 0, K_d > 0 \) and constant desired position \( q_d \), \( \dot{q}_d = 0 \). Prove global asymptotic stability of the error dynamics.

Solution:
The error dynamics are

\[ e = q - q_d,\quad \dot{e} = \dot{q},\quad J \ddot{e} + K_d \dot{e} + K_p e = 0. \]

Choose \( V(e,\dot{e}) = \tfrac{1}{2} J \dot{e}^2 + \tfrac{1}{2} K_p e^2 \). Clearly, \( V(e,\dot{e}) > 0 \) for all \( (e,\dot{e}) \neq (0,0) \) and \( V(0,0) = 0 \), and \( V \to \infty \) as \( \|(e,\dot{e})\| \to \infty \). Differentiating along trajectories:

\[ \dot{V} = J \dot{e} \ddot{e} + K_p e \dot{e} = \dot{e}(J \ddot{e} + K_p e) = \dot{e}(-K_d \dot{e}) = -K_d \dot{e}^2 \leq 0. \]

The largest invariant set inside \( \{(e,\dot{e}) : \dot{V} = 0\} \) is \( \{(0,0)\} \) because \( \dot{V} = 0 \) implies \( \dot{e} = 0 \), and then from the dynamics \( K_p e = 0 \), so \( e = 0 \). By LaSalle's invariance principle, the equilibrium \( (e,\dot{e}) = (0,0) \) is globally asymptotically stable.

Problem 2 (Robustness Radius from Lyapunov Inequality):
Let \( A_{\text{cl}} \) be Hurwitz and suppose that for some symmetric \( P > 0 \) and \( Q > 0 \), \( A_{\text{cl}}^{\top} P + P A_{\text{cl}} = -Q \). Show that any perturbed system \( \dot{x} = (A_{\text{cl}} + \Delta A)x \) is exponentially stable for all \( \Delta A \) satisfying \( \|\Delta A\|_2 < \lambda_{\min}(Q)/(2\|P\|_2) \).

Solution:
For \( V(x) = x^{\top} P x \), we have

\[ \begin{aligned} \dot{V} &= x^{\top}\big[(A_{\text{cl}}+\Delta A)^{\top}P + P(A_{\text{cl}}+\Delta A)\big]x \\ &= -x^{\top}Qx + x^{\top}(\Delta A^{\top}P + P\Delta A)x. \end{aligned} \]

Using the induced 2-norm and symmetry of \( P \), \( \|x^{\top}\Delta A^{\top}P x\| \leq \|P\|_2\|\Delta A\|_2\|x\|^2 \) and the same for \( x^{\top}P\Delta A x \), hence

\[ \dot{V} \leq -\lambda_{\min}(Q)\|x\|^2 + 2\|P\|_2\|\Delta A\|_2\|x\|^2. \]

If \( \|\Delta A\|_2 < \lambda_{\min}(Q)/(2\|P\|_2) \), then \( \dot{V} \leq -\varepsilon \|x\|^2 \) for some \( \varepsilon > 0 \), and exponential convergence follows from standard Lyapunov arguments.

Problem 3 (ISS-Like Bound for Disturbance Rejection):
Let the error dynamics be \( \dot{x} = f(x) + G d \), and suppose there exists a continuously differentiable \( V(x) \) and positive constants \( c_1,c_2,c_3,c_4 \) such that

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

Show that for any bounded disturbance with \( \|d(t)\| \leq \bar{d} \), the state satisfies \( \|x(t)\| \leq \beta(\|x(0)\|, t) + \gamma(\bar{d}) \) for suitable comparison functions \( \beta, \gamma \).

Solution (sketch):
The inequalities are an ISS-Lyapunov condition. Combining them yields

\[ \dot{V} \leq -\frac{c_3}{c_2} V + \frac{c_4}{c_1}\bar{d}^2. \]

This scalar differential inequality can be solved explicitly to give

\[ V(t) \leq V(0) e^{-\frac{c_3}{c_2} t} + \frac{c_4}{c_1}\bar{d}^2 \left( 1 - e^{-\frac{c_3}{c_2} t} \right). \]

Using \( c_1\|x(t)\|^2 \leq V(t) \) and rearranging yields \( \|x(t)\| \leq \beta(\|x(0)\|,t) + \gamma(\bar{d}) \) with \( \beta(r,t) = \sqrt{\tfrac{c_2}{c_1}} r e^{-\frac{c_3}{2c_2} t} \) and \( \gamma(s) = \sqrt{\tfrac{c_4}{c_1 c_3}}\,s \).

Problem 4 (Empirical Robustness Index):
You perform \( N \) Monte Carlo simulations with independently sampled uncertainty and disturbances. Let \( J_i \) be a scalar performance measure from the \( i \)-th run (e.g., maximum tracking error). Define the empirical worst-case index \( J_{\max} = \max_i J_i \) and the empirical \( 95\% \) quantile \( J_{0.95} \). Explain how you would use \( J_{\max} \) and \( J_{0.95} \) to argue about robustness in your capstone report, assuming you have formal stability guarantees from Lyapunov analysis.

Solution:
The formal Lyapunov analysis guarantees stability for all uncertainties within a specified set (e.g., bounded by \( \|\Delta A\|_2 < \bar{\delta} \)) but does not directly quantify performance. Monte Carlo sampling within and slightly beyond that set gives numerical evidence about typical and extreme performance. The value \( J_{\max} \) provides a conservative upper bound on observed performance degradation under the tested scenarios, while \( J_{0.95} \) characterizes performance in the majority (95%) of cases. If both indices lie below design thresholds (e.g., maximum allowable tracking error, saturation limits), you can argue that not only is the system provably stable, but it is also robust in the sense of meeting quantitative performance specifications under realistic uncertainty and disturbance levels.

In this lesson you connected the theoretical notions of Lyapunov stability, robust stability margins, and input–output robustness to a concrete evaluation pipeline for your capstone controller. You saw how to derive Lyapunov-based robustness radii, compute eigenvalue-based stability and damping margins, and estimate disturbance attenuation via \( \mathcal{L}_2 \)-gains. Complementing these analytic tools, you implemented Monte Carlo robustness testing and stability checks in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, preparing you to make rigorous stability/robustness claims in your final project report.