By Chapter 13 you have seen how ideal rigid-body dynamics of an n-DOF robot can be written as

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

where \(q,\dot{q},\ddot{q} \in \mathbb{R}^n\) are joint position, velocity, acceleration, \(\mathbf{M}(q)\) is the inertia matrix, \(\mathbf{C}\) is the Coriolis/centrifugal matrix, \(\mathbf{g}\) gravity, \(\mathbf{f}\) non-ideal effects (e.g. friction), and \(\mathbf{d}(t)\) represents external disturbances.

Typical sources of uncertainty are:

• Parametric uncertainty: in link masses, centers of mass, inertia tensors, friction coefficients, gear ratios, etc.
• Structural/unmodeled dynamics: flexible links, joint backlash, high-frequency actuator dynamics that are not included in the model.
• Exogenous disturbances: contact forces, payload variations, cable drag, aerodynamic drag.
• State estimation errors: encoder quantization, sensor bias and noise that corrupt \(q,\dot{q}\).

Our goal is not to eliminate uncertainty but to model it explicitly so that we can attach rigorous error bounds to quantities such as predicted joint torques, accelerations, or end-effector motion.

flowchart TD
  SRC["Physical robot"] --> U1["Identify sources: parameters, unmodeled, disturbances"]
  U1 --> C1["Classify: parametric vs additive vs measurement"]
  C1 --> R1["Choose representation: interval / norm-bounded / stochastic"]
  R1 --> BND["Derive error bounds on tau, qddot"]
  BND --> V1["Validate bounds with experiments or simulation"]
  V1 --> REF["Refine model and uncertainty sets"]
        

In practice we only have access to a nominal model, obtained, for example, from CAD data and simple experiments. Denote the nominal model by

\[ \hat{\mathbf{M}}(q)\ddot{q} + \hat{\mathbf{C}}(q,\dot{q})\dot{q} + \hat{\mathbf{g}}(q) + \hat{\mathbf{f}}(\dot{q}) = \boldsymbol{\tau}. \]

The true dynamics can be written as

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

Define parameterized errors \( \Delta\mathbf{M} = \mathbf{M} - \hat{\mathbf{M}} \), \( \Delta\mathbf{C} = \mathbf{C} - \hat{\mathbf{C}} \), \( \Delta\mathbf{g} = \mathbf{g} - \hat{\mathbf{g}} \), \( \Delta\mathbf{f} = \mathbf{f} - \hat{\mathbf{f}} \). Subtracting the nominal from the true model yields the dynamic mismatch

\[ \Delta\mathbf{M}(q)\ddot{q} + \Delta\mathbf{C}(q,\dot{q})\dot{q} + \Delta\mathbf{g}(q) + \Delta\mathbf{f}(\dot{q}) = \mathbf{d}(t). \]

If we group all modeling errors and disturbances into a single term \( \mathbf{w}(t) \in \mathbb{R}^n \), we can write

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

where \( \mathbf{w}(t) \) is an additive uncertainty term. A common modeling assumption in dynamics and control analysis is that \( \mathbf{w}(t) \) is bounded:

\[ \|\mathbf{w}(t)\| \le \bar{w}, \quad \forall t, \]

for some known constant \( \bar{w} > 0 \). This leads to robust guarantees that hold for all disturbances of magnitude at most \( \bar{w} \).

From Chapter 13, the rigid-body manipulator dynamics can be written in a linearly parameterized form

\[ \boldsymbol{\tau}(t) = \mathbf{Y}\big(q(t),\dot{q}(t),\ddot{q}(t)\big)\,\boldsymbol{\theta}, \]

where \( \boldsymbol{\theta} \in \mathbb{R}^p \) collects inertial and friction parameters and \( \mathbf{Y} \in \mathbb{R}^{n \times p} \) is the regressor matrix, known from kinematics and desired motion.

In identification we estimate \( \hat{\boldsymbol{\theta}} \) so that the nominal model becomes

\[ \hat{\boldsymbol{\tau}}(t) = \mathbf{Y}\big(q(t),\dot{q}(t),\ddot{q}(t)\big)\,\hat{\boldsymbol{\theta}}. \]

Define the parameter error \( \tilde{\boldsymbol{\theta}} = \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \). Then the torque prediction error is

\[ \boldsymbol{\tau}(t) - \hat{\boldsymbol{\tau}}(t) = \mathbf{Y}\big(q,\dot{q},\ddot{q}\big)\,\tilde{\boldsymbol{\theta}}. \]

If we have a norm bound on the parameter error \( \|\tilde{\boldsymbol{\theta}}\| \le \bar{\theta} \), then the induced torque error is bounded by the operator norm inequality

\[ \|\boldsymbol{\tau} - \hat{\boldsymbol{\tau}}\| \le \|\mathbf{Y}\|\,\|\tilde{\boldsymbol{\theta}}\| \le \|\mathbf{Y}\|\,\bar{\theta}, \]

where \( \|\mathbf{Y}\| \) is any sub-multiplicative matrix norm (e.g., spectral or Frobenius norm). This provides a direct and computationally cheap way to compute a torque error bound given parameter uncertainties.

flowchart TD
  QDD["q, dq, qdd (trajectory)"] --> YMAT["Compute regressor Y(q,dq,qdd)"]
  YMAT --> THETA["True parameters theta"]
  YMAT --> THETAH["Estimated parameters theta_hat"]
  THETA --> TAU["True tau = Y * theta"]
  THETAH --> TAUH["Nominal tau_hat = Y * theta_hat"]
  TAU --> ERR["Error = tau - tau_hat"]
  TAUH --> ERR
        

Consider the simplified case where the only uncertainty is in the inertia matrix:

\[ \mathbf{M}(q)\ddot{q} = \boldsymbol{\tau}, \quad \hat{\mathbf{M}}(q)\ddot{q}_{\text{nom}} = \boldsymbol{\tau}. \]

The true and nominal accelerations are \( \ddot{q} = \mathbf{M}^{-1}\boldsymbol{\tau} \), \( \ddot{q}_{\text{nom}} = \hat{\mathbf{M}}^{-1}\boldsymbol{\tau} \). Let \( \Delta\mathbf{M} = \mathbf{M}-\hat{\mathbf{M}} \). Assume \( \hat{\mathbf{M}} \) is positive definite and \( \|\hat{\mathbf{M}}^{-1}\Delta\mathbf{M}\| < 1 \), so that \( \mathbf{M} \) is also invertible. A Neumann series argument yields

\[ \mathbf{M}^{-1} = \big(\hat{\mathbf{M}} + \Delta\mathbf{M}\big)^{-1} = \hat{\mathbf{M}}^{-1} \big(\mathbf{I} + \hat{\mathbf{M}}^{-1}\Delta\mathbf{M}\big)^{-1}. \]

Hence the acceleration error is

\[ \ddot{q} - \ddot{q}_{\text{nom}} = \big(\mathbf{M}^{-1} - \hat{\mathbf{M}}^{-1}\big)\boldsymbol{\tau} = \hat{\mathbf{M}}^{-1} \Big[\big(\mathbf{I} + \hat{\mathbf{M}}^{-1}\Delta\mathbf{M}\big)^{-1} - \mathbf{I}\Big] \boldsymbol{\tau}. \]

Using the identity \( (\mathbf{I}+\mathbf{X})^{-1} - \mathbf{I} = -\mathbf{X}(\mathbf{I}+\mathbf{X})^{-1} \) with \( \mathbf{X} = \hat{\mathbf{M}}^{-1}\Delta\mathbf{M} \), we obtain

\[ \ddot{q} - \ddot{q}_{\text{nom}} = -\hat{\mathbf{M}}^{-1} \hat{\mathbf{M}}^{-1}\Delta\mathbf{M} \big(\mathbf{I} + \hat{\mathbf{M}}^{-1}\Delta\mathbf{M}\big)^{-1} \boldsymbol{\tau}. \]

Therefore the acceleration error norm can be bounded as

\[ \|\ddot{q} - \ddot{q}_{\text{nom}}\| \le \|\hat{\mathbf{M}}^{-1}\|^2 \frac{\|\Delta\mathbf{M}\|}{1-\|\hat{\mathbf{M}}^{-1}\Delta\mathbf{M}\|} \,\|\boldsymbol{\tau}\|. \]

In practice we may know a bound \( \|\Delta\mathbf{M}(q)\| \le \bar{M}_{\Delta} \) over a region of interest in configuration space, which leads to a global bound on the acceleration error for all admissible torques.

An alternative viewpoint is interval uncertainty on individual parameters. Suppose each parameter \( \theta_i \) satisfies \( \theta_i \in [\underline{\theta}_i,\overline{\theta}_i] \). Then the admissible parameter set is an axis-aligned box \(\Theta \subset \mathbb{R}^p\).

For any fixed regressor row \( \mathbf{y}_j^{\top} \) (corresponding to joint \( j \)), the torque component is \( \tau_j = \mathbf{y}_j^{\top}\boldsymbol{\theta} \). The worst-case deviation of \( \tau_j \) from its nominal value \( \hat{\tau}_j = \mathbf{y}_j^{\top}\hat{\boldsymbol{\theta}} \) is

\[ \big|\tau_j - \hat{\tau}_j\big| = \Big|\sum_{i=1}^p y_{ji}\,\tilde{\theta}_i\Big| \le \sum_{i=1}^p |y_{ji}|\,\delta_i, \]

where \( \delta_i = \max\big( |\overline{\theta}_i - \hat{\theta}_i|, |\underline{\theta}_i - \hat{\theta}_i| \big) \).

This bound is tight if we can choose parameter combinations independently inside their intervals. In robotics, dependence between parameters (e.g., inertia tensor constraints) can be captured by replacing the box set \( \Theta \) with a polytope or ellipsoid and solving small optimization problems to find worst-case torques.

For a 1-DOF pendulum with dynamics \( I\ddot{q} + b\dot{q} + mgl\sin q = \tau \) and parameter intervals \( I \in [\underline{I},\overline{I}] \), \( b \in [\underline{b},\overline{b}] \), \( (mgl) \in [\underline{k},\overline{k}] \), the regressor is \( \mathbf{y} = [\ddot{q},\dot{q},\sin q] \). The worst-case torque deviation at a given motion \( (q,\dot{q},\ddot{q}) \) is

\[ \big|\tau - \hat{\tau}\big| \le |\ddot{q}|\,\delta_I + |\dot{q}|\,\delta_b + |\sin q|\,\delta_k, \]

where \( \delta_I,\delta_b,\delta_k \) are the scalar parameter interval radii.

When parameter estimates are obtained from noisy experimental data, it is natural to model \( \boldsymbol{\theta} \) as a random vector with mean \( \hat{\boldsymbol{\theta}} \) and covariance \( \boldsymbol{\Sigma}_{\theta} \succ 0 \). Assuming an unbiased estimator, the torque becomes a random quantity

\[ \boldsymbol{\tau} = \mathbf{Y}\boldsymbol{\theta}, \quad \mathbb{E}[\boldsymbol{\tau}] = \mathbf{Y}\hat{\boldsymbol{\theta}}. \]

Since the mapping is linear, the covariance of the torque is

\[ \boldsymbol{\Sigma}_{\tau} = \mathbf{Y}\,\boldsymbol{\Sigma}_{\theta}\,\mathbf{Y}^{\top}. \]

If, in addition, \( \boldsymbol{\theta} \) is Gaussian, then \( \boldsymbol{\tau} \) is also Gaussian. Even without Gaussian assumptions, the covariance expression is valid and can be used to construct probabilistic bounds of the form

\[ \mathbb{P}\Big( \big\|\boldsymbol{\tau} - \mathbf{Y}\hat{\boldsymbol{\theta}}\big\| \le \gamma \Big) \approx 1-\alpha, \]

where \( \gamma \) can be computed from the eigenvalues of \( \boldsymbol{\Sigma}_{\tau} \) and an assumed confidence level \( 1-\alpha \). This connects parameter estimation uncertainty to confidence regions on predicted torques.

We illustrate interval and norm-bounded error computation for the 1-DOF pendulum dynamics \( I\ddot{q} + b\dot{q} + k\sin q = \tau \) with \( k = mgl \). This can be seen as a special case of a single-joint link from a more complex manipulator.

import numpy as np

def pendulum_regressor(q, dq, ddq):
    """
    Regressor for I*ddq + b*dq + k*sin(q) = tau.
    Returns Y (1x3) so that tau = Y @ theta, theta = [I, b, k].
    """
    return np.array([[ddq, dq, np.sin(q)]], dtype=float)

# Nominal parameters
theta_hat = np.array([0.05, 0.01, 0.5])  # [I, b, k]

# Interval radii (absolute)
delta = np.array([0.005, 0.002, 0.05])   # +-10% for I, 20% for b, 10% for k

# Example trajectory point
q = np.deg2rad(30.0)
dq = 0.5
ddq = 1.0
Y = pendulum_regressor(q, dq, ddq)

# Nominal torque
tau_hat = float(Y @ theta_hat)

# Worst-case interval bound (element-wise)
tau_interval_radius = float(np.sum(np.abs(Y) * delta))
print("Nominal tau:", tau_hat)
print("Interval bound on |tau - tau_hat|:", tau_interval_radius)

# Norm-bounded uncertainty: treat theta error vector with 2-norm bound
theta_bar = np.linalg.norm(delta, ord=2)  # crude upper bound
Y_norm = np.linalg.norm(Y, ord=2)
tau_norm_bound = Y_norm * theta_bar
print("2-norm based bound on ||tau - tau_hat||:", tau_norm_bound)

# Monte Carlo sampling inside the parameter box to empirically check bounds
rng = np.random.default_rng(0)
N = 10000
errors = []
for _ in range(N):
    # Sample each parameter uniformly inside [theta_hat - delta, theta_hat + delta]
    theta_sample = theta_hat + rng.uniform(low=-delta, high=delta)
    tau_sample = float(Y @ theta_sample)
    errors.append(abs(tau_sample - tau_hat))

errors = np.array(errors)
print("Max error over samples:", errors.max())
print("Interval bound is valid?:", errors.max() <= tau_interval_radius + 1e-9)
      

In larger manipulators, the same pattern applies with \( \mathbf{Y} \in \mathbb{R}^{n \times p} \) obtained, for example, using roboticstoolbox-python to compute dynamics and then extracting the regressor for a given trajectory.

In C++, libraries such as RBDL or Pinocchio can compute \( \mathbf{M}(q) \), \( \mathbf{C}(q,\dot{q}) \), and \( \mathbf{g}(q) \). Here we show a minimal example for the 1-DOF pendulum using only the regressor idea (no external library).

#include <iostream>
#include <cmath>
#include <array>

using Vec3 = std::array<double, 3>;

Vec3 pendulumRegressor(double q, double dq, double ddq) {
    // Y = [ddq, dq, sin(q)]
    return {ddq, dq, std::sin(q)};
}

double dot3(const Vec3& a, const Vec3& b) {
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}

int main() {
    // Nominal parameters theta_hat = [I, b, k]
    Vec3 theta_hat{0.05, 0.01, 0.5};
    // Interval radii delta
    Vec3 delta{0.005, 0.002, 0.05};

    double q = 30.0 * M_PI / 180.0;
    double dq = 0.5;
    double ddq = 1.0;

    Vec3 Y = pendulumRegressor(q, dq, ddq);
    double tau_hat = dot3(Y, theta_hat);

    // Interval bound: sum |Y_i| * delta_i
    double tau_interval_radius =
        std::fabs(Y[0]) * delta[0] +
        std::fabs(Y[1]) * delta[1] +
        std::fabs(Y[2]) * delta[2];

    std::cout << "Nominal tau: " << tau_hat << "\n";
    std::cout << "Interval bound on |tau - tau_hat|: "
              << tau_interval_radius << "\n";

    return 0;
}
      

When using RBDL or Pinocchio, one can obtain the nominal dynamics and then perturb the inertial parameters inside the model structure to empirically check the theoretical bounds derived earlier.

Java is less common for low-level robotics dynamics, but libraries like EJML provide convenient matrix operations. The following snippet performs Monte Carlo sampling of parameter uncertainty for the pendulum example.

import org.ejml.simple.SimpleMatrix;
import java.util.Random;

public class PendulumUncertaintyJava {

    public static SimpleMatrix regressor(double q, double dq, double ddq) {
        double[] data = {ddq, dq, Math.sin(q)};
        return new SimpleMatrix(1, 3, true, data);
    }

    public static void main(String[] args) {
        double[] thetaHat = {0.05, 0.01, 0.5};
        double[] delta = {0.005, 0.002, 0.05};

        SimpleMatrix thetaHatMat =
            new SimpleMatrix(3, 1, true, thetaHat);

        double q = Math.toRadians(30.0);
        double dq = 0.5;
        double ddq = 1.0;

        SimpleMatrix Y = regressor(q, dq, ddq);
        double tauHat = Y.mult(thetaHatMat).get(0, 0);

        // Analytic interval bound
        double tauIntervalRadius =
            Math.abs(Y.get(0, 0)) * delta[0] +
            Math.abs(Y.get(0, 1)) * delta[1] +
            Math.abs(Y.get(0, 2)) * delta[2];

        System.out.println("Nominal tau: " + tauHat);
        System.out.println("Interval bound |tau - tauHat|: " + tauIntervalRadius);

        // Monte Carlo check
        Random rng = new Random(0);
        int N = 10000;
        double maxErr = 0.0;
        for (int k = 0; k < N; ++k) {
            double[] thetaSample = new double[3];
            for (int i = 0; i < 3; ++i) {
                double u = rng.nextDouble(); // in [0,1]
                thetaSample[i] = thetaHat[i] + (2.0 * u - 1.0) * delta[i];
            }
            SimpleMatrix thetaSampleMat =
                new SimpleMatrix(3, 1, true, thetaSample);
            double tauSample = Y.mult(thetaSampleMat).get(0, 0);
            double err = Math.abs(tauSample - tauHat);
            if (err > maxErr) {
                maxErr = err;
            }
        }
        System.out.println("Max error over samples: " + maxErr);
    }
}
      

Such Monte Carlo experiments complement analytic bounds and help build intuition about the conservatism of interval and norm-bounded models.

10.1 MATLAB/Simulink

MATLAB with the Robotics System Toolbox can represent the manipulator dynamics, while Robust Control Toolbox provides explicit uncertain objects. For the pendulum we can use ureal to encode parameter intervals.

% Pendulum parameters as uncertain reals (10% relative uncertainty)
I = ureal('I', 0.05, 'Percentage', 10);
b = ureal('b', 0.01, 'Percentage', 20);
k = ureal('k', 0.5,  'Percentage', 10);

% Regressor at a given state (q, dq, ddq)
q   = deg2rad(30);
dq  = 0.5;
ddq = 1.0;
Y   = [ddq, dq, sin(q)];

theta_hat = [0.05; 0.01; 0.5];
tau_hat   = Y * theta_hat;

% Uncertain parameter vector
theta_unc = [I; b; k];
tau_unc   = Y * theta_unc;   % This is an uncertain scalar (uss)

% Sample the uncertain parameters and compute tau
Ns = 1000;
samples = usample(theta_unc, Ns);  % Ns-by-1 struct of samples
tau_samples = zeros(Ns, 1);
for i = 1:Ns
    th = [samples(i).I; samples(i).b; samples(i).k];
    tau_samples(i) = Y * th;
end

max_err = max(abs(tau_samples - tau_hat));
fprintf('Max |tau - tau_hat| over samples: %g\n', max_err);

% Simulink modeling idea (conceptual):
%  - Build a Simulink block diagram for q, dq, ddq dynamics.
%  - Implement the dynamics I*ddq + b*dq + k*sin(q) = tau using Gain and Sum blocks.
%  - Use the Uncertain State Space blocks or mask the gains with uncertain
%    parameters I, b, k to perform worst-case and Monte Carlo simulations.
      

10.2 Wolfram Mathematica

Mathematica is well-suited for symbolic expressions and interval arithmetic. The following code computes symbolic torque expressions and interval bounds for the pendulum.

(* Symbolic parameters and states *)
ClearAll["Global`*"];
I   = Symbol["I"];
b   = Symbol["b"];
k   = Symbol["k"];
q   = Symbol["q"];
dq  = Symbol["dq"];
ddq = Symbol["ddq"];

tauExpr = I*ddq + b*dq + k*Sin[q];

(* Nominal and interval values *)
thetaHat = {I -> 0.05, b -> 0.01, k -> 0.5};
delta    = {0.005, 0.002, 0.05};

qVal   = 30 Degree;
dqVal  = 0.5;
ddqVal = 1.0;

tauHat = tauExpr /. thetaHat /. {q -> qVal, dq -> dqVal, ddq -> ddqVal} // N

(* Interval parameters *)
Iint = Interval[{0.05 - delta[[1]], 0.05 + delta[[1]]}];
bint = Interval[{0.01 - delta[[2]], 0.01 + delta[[2]]}];
kint = Interval[{0.5  - delta[[3]], 0.5  + delta[[3]]}];

tauInterval =
  (Iint*ddqVal + bint*dqVal + kint*Sin[qVal]) // Simplify

(* Extract numeric bounds and radius *)
{tauMin, tauMax} = IntervalBounds[tauInterval];
tauRadius = (tauMax - tauMin)/2.0
      

The interval radius computed symbolically corresponds directly to the worst-case deviation bounds derived analytically in Section 5.

In this lesson we formulated uncertainty in robot dynamics at the equation-of-motion level, distinguishing additive disturbances from parametric uncertainty and exploiting the linear parameterization of manipulator dynamics. We derived norm-based and interval-based bounds on torque and acceleration errors, and related stochastic parameter covariance to torque covariance. Implementation examples in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica demonstrated how these concepts can be realized computationally for simple systems such as a 1-DOF pendulum, providing a template for more complex manipulators.