Real actuators (electric motors, hydraulic cylinders, etc.) cannot generate arbitrarily large torque, force, velocity, or power. At the dynamics level we mostly care about torque limits, expressed as bounds on joint torques \( \boldsymbol{\tau} \in \mathbb{R}^n \). For a manipulator with \( n \) joints and dynamics

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

saturation is typically modeled as box constraints on the actual actuator torque:

\[ \boldsymbol{\tau}_{\min} \;\leq\; \boldsymbol{\tau} \;\leq\; \boldsymbol{\tau}_{\max}, \]

with the inequalities interpreted componentwise. The commanded torque \( \boldsymbol{\tau}^{\mathrm{cmd}} \) is transformed by a saturation nonlinearity of the actuators and drive electronics to produce the realized torque \( \boldsymbol{\tau} \).

At a given configuration and velocity, these torque bounds imply a feasible set of accelerations, and through the Jacobian a feasible set of task-space wrenches. This view is crucial when dynamics is later combined with control, trajectory generation, or safety analysis.

flowchart TD
  DYN["Dynamics: M(q) qddot + h(q,qd) = tau + tau_ext"]
    --> BOX["Torque limits: tau_min <= tau <= tau_max"]
  BOX --> ACC["Feasible joint accelerations qddot"]
  ACC --> WR["Feasible task wrenches via J^T F = tau"]
  WR --> IMP["Impact on reachable motions & performance"]
        

For a single scalar actuator (one joint), the static saturation map is the piecewise-affine function

\[ s(u) = \begin{cases} u_{\max}, & u \geq u_{\max}, \\ u, & u_{\min} \leq u \leq u_{\max}, \\ u_{\min}, & u \leq u_{\min}, \end{cases} \]

where \( u \) is the commanded actuator input (e.g., torque command), and \( u_{\min} < 0 < u_{\max} \) are the minimum and maximum achievable torques.

In vector form, with \( \boldsymbol{\tau}^{\mathrm{cmd}} \), \( \boldsymbol{\tau}_{\min} \), \( \boldsymbol{\tau}_{\max} \) in \( \mathbb{R}^n \), the saturation map is applied componentwise:

\[ \boldsymbol{\tau} = \operatorname{sat} \!\left( \boldsymbol{\tau}^{\mathrm{cmd}}; \boldsymbol{\tau}_{\min}, \boldsymbol{\tau}_{\max} \right), \quad \bigl(\operatorname{sat}(\cdot)\bigr)_i = \min\bigl\{ \tau_{\max,i}, \max\{ \tau_{\min,i}, \tau^{\mathrm{cmd}}_i \} \bigr\}. \]

This is a memoryless (static) nonlinearity, monotone, and non-expansive in the Euclidean norm, i.e. \( \|\operatorname{sat}(\mathbf{u}) - \operatorname{sat}(\mathbf{v})\| \leq \|\mathbf{u}-\mathbf{v}\| \). The proof is componentwise: for each scalar actuator, the derivative of \( s(u) \) with respect to \( u \) lies in \( [0,1] \), hence the Lipschitz constant is at most 1.

Recall the standard manipulator dynamics in joint space:

\[ \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) = \boldsymbol{\tau} + \boldsymbol{\tau}_{\mathrm{ext}}, \]

where \( \mathbf{M}(\mathbf{q}) \) is the positive definite inertia matrix and \( \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) \) collects Coriolis/centrifugal, gravity, and friction terms introduced earlier in this chapter.

Solving for the acceleration gives

\[ \ddot{\mathbf{q}} = \mathbf{M}(\mathbf{q})^{-1} \left( \boldsymbol{\tau} + \boldsymbol{\tau}_{\mathrm{ext}} - \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) \right). \]

Imposing torque bounds \( \boldsymbol{\tau}_{\min} \leq \boldsymbol{\tau} \leq \boldsymbol{\tau}_{\max} \) is equivalent to inequality constraints on accelerations:

\[ \boldsymbol{\tau}_{\min} \;\leq\; \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) - \boldsymbol{\tau}_{\mathrm{ext}} \;\leq\; \boldsymbol{\tau}_{\max}. \]

For fixed \( \mathbf{q}, \dot{\mathbf{q}} \) and \( \boldsymbol{\tau}_{\mathrm{ext}} \), the acceleration \( \ddot{\mathbf{q}} \) lies in a convex polyhedron:

\[ \mathcal{A}(\mathbf{q},\dot{\mathbf{q}}) = \left\{ \ddot{\mathbf{q}} \in \mathbb{R}^n \;\middle|\; \mathbf{A}(\mathbf{q}) \ddot{\mathbf{q}} \leq \mathbf{b}(\mathbf{q},\dot{\mathbf{q}}) \right\}, \]

where the inequality matrix and vector can be written as

\[ \mathbf{A}(\mathbf{q}) = \begin{bmatrix} \mathbf{M}(\mathbf{q}) \\[0.3em] -\mathbf{M}(\mathbf{q}) \end{bmatrix}, \quad \mathbf{b}(\mathbf{q},\dot{\mathbf{q}}) = \begin{bmatrix} \boldsymbol{\tau}_{\max} - \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) + \boldsymbol{\tau}_{\mathrm{ext}} \\ -\boldsymbol{\tau}_{\min} + \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) - \boldsymbol{\tau}_{\mathrm{ext}} \end{bmatrix}. \]

Thus, saturation does not merely clip torques; it carves out a polyhedral cone or polytope in acceleration space, depending on whether the bounds are symmetric and whether gravity is canceled.

Let \( \mathbf{J}(\mathbf{q}) \) be the manipulator Jacobian mapping joint velocities to end-effector twists. For a quasi-static mapping between end-effector wrench \( \mathbf{F} \in \mathbb{R}^m \) and joint torque we have

\[ \boldsymbol{\tau} = \mathbf{J}(\mathbf{q})^\top \mathbf{F} + \boldsymbol{\tau}_{\mathrm{bias}}, \]

where \( \boldsymbol{\tau}_{\mathrm{bias}} \) represents any gravity or friction compensation torques already applied. Imposing torque limits gives

\[ \boldsymbol{\tau}_{\min} \leq \mathbf{J}(\mathbf{q})^\top \mathbf{F} + \boldsymbol{\tau}_{\mathrm{bias}} \leq \boldsymbol{\tau}_{\max}. \]

The set of feasible wrenches at a given configuration is therefore

\[ \mathcal{W}(\mathbf{q}) = \left\{ \mathbf{F} \in \mathbb{R}^m \;\middle|\; \boldsymbol{\tau}_{\min} - \boldsymbol{\tau}_{\mathrm{bias}} \leq \mathbf{J}(\mathbf{q})^\top \mathbf{F} \leq \boldsymbol{\tau}_{\max} - \boldsymbol{\tau}_{\mathrm{bias}} \right\}. \]

This is again a convex polytope (intersection of half-spaces) in wrench space. Its shape and volume depend on both actuator limits and the Jacobian, hence on configuration. As the robot approaches kinematic singularities, the geometry of \( \mathcal{W}(\mathbf{q}) \) reflects the loss of directional force capability.

flowchart TD
  Q["Joint space q"] --> J["Jacobian J(q)"]
  J --> FSET["Torque-limited wrench set W(q)"]
  FSET --> TASK["Task-level capabilities (force directions, magnitudes)"]
        

A single revolute joint with inertia \( I > 0 \) and viscous friction coefficient \( b \geq 0 \) has dynamics

\[ I \ddot{q} + b \dot{q} + g(q) = \tau, \quad \tau_{\min} \leq \tau \leq \tau_{\max}, \]

where \( g(q) \) is the gravity torque. Solving for acceleration:

\[ \ddot{q} = \frac{1}{I} \left( \tau - b \dot{q} - g(q) \right). \]

The acceleration bounds implied by torque limits are

\[ \ddot{q}_{\min} = \frac{1}{I} \left( \tau_{\min} - b \dot{q} - g(q) \right), \quad \ddot{q}_{\max} = \frac{1}{I} \left( \tau_{\max} - b \dot{q} - g(q) \right). \]

Even this simple example shows that actuator saturation translates into state-dependent acceleration limits. When \( \dot{q} \) is large, friction consumes some of the available torque budget, shrinking the feasible acceleration interval.

In Python, robotics modeling is commonly supported by packages such as numpy, scipy, and roboticstoolbox-python. Below we implement a small saturation utility and a function that computes feasible accelerations for a 1-DOF joint, as in the previous section.

import numpy as np

def sat_scalar(u, u_min, u_max):
    """
    Scalar saturation.
    """
    return max(u_min, min(u, u_max))

def sat_vector(u, u_min, u_max):
    """
    Componentwise saturation of vectors (numpy arrays).
    """
    u = np.asarray(u, dtype=float)
    u_min = np.asarray(u_min, dtype=float)
    u_max = np.asarray(u_max, dtype=float)
    return np.minimum(u_max, np.maximum(u_min, u))

def feasible_acc_1d(q, qd, tau_min, tau_max, I, b, g_func):
    """
    Compute the min and max acceleration consistent with torque limits
    for a 1-DOF joint.
    """
    g = g_func(q)
    ddq_min = (tau_min - b * qd - g) / I
    ddq_max = (tau_max - b * qd - g) / I
    return ddq_min, ddq_max

# Example usage:
I = 0.2      # kg m^2
b = 0.05     # N m s/rad
tau_min = -5.0
tau_max = 5.0

def g_func(q):
    # Simple gravity torque model, e.g. for a pendulum-like link
    m = 1.0
    l = 0.5
    g0 = 9.81
    return m * g0 * l * np.sin(q)

q = 0.5       # rad
qd = 0.0      # rad/s

# Commanded torque (from a controller, not modeled here)
tau_cmd = 8.0
tau_act = sat_scalar(tau_cmd, tau_min, tau_max)

ddq_min, ddq_max = feasible_acc_1d(q, qd, tau_min, tau_max, I, b, g_func)

print("Actual torque:", tau_act)
print("Feasible acceleration interval:", ddq_min, ddq_max)
      

In a full multi-DOF robot model (e.g. using roboticstoolbox-python), the same ideas extend by substituting symbolic or numeric inertia matrices and computing the polyhedral set \( \mathcal{A}(\mathbf{q},\dot{\mathbf{q}}) \).

In C++, robotics dynamics are often implemented with Eigen-based libraries such as RBDL or Pinocchio. Here we show a simple saturation in plain Eigen vectors, which can be wrapped around the torque output of any dynamics or control module.

#include <Eigen/Dense>

Eigen::VectorXd saturateTorque(
    const Eigen::VectorXd& tau_cmd,
    const Eigen::VectorXd& tau_min,
    const Eigen::VectorXd& tau_max)
{
    Eigen::VectorXd tau = tau_cmd;
    for (int i = 0; i < tau.size(); ++i) {
        if (tau(i) > tau_max(i)) {
            tau(i) = tau_max(i);
        } else if (tau(i) < tau_min(i)) {
            tau(i) = tau_min(i);
        }
    }
    return tau;
}

// Example: wrap inverse dynamics (not shown) with saturation
Eigen::VectorXd computeSaturatedTorque(
    const Eigen::VectorXd& q,
    const Eigen::VectorXd& qd,
    const Eigen::VectorXd& qdd_des,
    const Eigen::VectorXd& tau_min,
    const Eigen::VectorXd& tau_max)
{
    // Suppose invDynamics(q, qd, qdd_des) is provided by a dynamics library
    Eigen::VectorXd tau_cmd = invDynamics(q, qd, qdd_des);
    Eigen::VectorXd tau_act = saturateTorque(tau_cmd, tau_min, tau_max);
    return tau_act;
}
      

When integrating this into a full dynamic simulator, one must be aware that saturation introduces nonlinearity and can cause loss of passivity if not handled carefully, especially when combined with idealized controllers (covered in later courses).

Java-based robotics stacks (for example, some ROS Java bindings or custom simulation tools) often use numerical libraries such as EJML or Apache Commons Math. The core saturation logic, however, is simple array manipulation:

public class TorqueSaturation {

    public static double saturateScalar(double u, double uMin, double uMax) {
        if (u > uMax) {
            return uMax;
        } else if (u < uMin) {
            return uMin;
        } else {
            return u;
        }
    }

    public static double[] saturateVector(
            double[] tauCmd,
            double[] tauMin,
            double[] tauMax) {

        int n = tauCmd.length;
        double[] tau = new double[n];
        for (int i = 0; i < n; i++) {
            double u = tauCmd[i];
            double uMin = tauMin[i];
            double uMax = tauMax[i];
            if (u > uMax) {
                tau[i] = uMax;
            } else if (u < uMin) {
                tau[i] = uMin;
            } else {
                tau[i] = u;
            }
        }
        return tau;
    }
}
      

A Java dynamics engine would call saturateVector on the torques computed from its equations of motion before applying them to the numerical integrator that advances joint states.

In MATLAB with the Robotics System Toolbox, joint dynamics are often simulated either via scripts or Simulink. Simulink offers a built-in Saturation block that can be placed on the torque signal line. The following MATLAB function implements scalar and vector saturation for use in script-based simulations:

function tau = saturateTorque(tau_cmd, tau_min, tau_max)
%SATURATETORQUE Componentwise saturation of joint torques.
%
% tau_cmd, tau_min, tau_max: column vectors (n x 1)

tau = tau_cmd;
for i = 1:length(tau_cmd)
    if tau(i) > tau_max(i)
        tau(i) = tau_max(i);
    elseif tau(i) < tau_min(i)
        tau(i) = tau_min(i);
    end
end
end

% Example usage in a time-stepping loop:
%
% for k = 1:N
%   tau_cmd = M(q) * qdd_des(:,k) + h(q, qd);
%   tau_act = saturateTorque(tau_cmd, tau_min, tau_max);
%   qdd = M(q) \ (tau_act - h(q, qd));
%   % Integrate qd, q ...
% end
      

In Simulink, one typically connects the output of an inverse dynamics block (built from RigidBodyTree models) to a Saturation block whose upper and lower limits are given by \( \boldsymbol{\tau}_{\max} \) and \( \boldsymbol{\tau}_{\min} \). The saturated torque then feeds the joint integrator subsystem.

Wolfram Mathematica is well-suited for symbolic manipulator dynamics. Below we define a saturation function and derive acceleration bounds for the 1-DOF example symbolically.

(* Scalar saturation *)
sat[u_, umin_, umax_] :=
  Piecewise[{
    {umax, u >= umax},
    {u, umin <= u <= umax},
    {umin, u <= umin}
  }];

(* 1-DOF dynamics parameters *)
Clear[I, b, q, qd, gq, tau, tauMin, tauMax];
gq = Symbol["g"][q];

(* Acceleration in terms of torque *)
ddq[tau_] := (tau - b*qd - gq)/I;

(* Symbolic acceleration bounds *)
ddqMin = Simplify[ddq[tauMin]];
ddqMax = Simplify[ddq[tauMax]];

{ddqMin, ddqMax}
      

Using Mathematica's RegionPlot and related functions, one can visualize feasible acceleration regions or wrench polytopes for low-dimensional subsystems and investigate how they vary with configuration.

In this lesson we treated actuator saturation as a formal constraint in the dynamics of robot manipulators. Starting from torque bounds, we derived convex polyhedral sets of feasible accelerations and task-space wrenches, showing how saturation is a structured but nonlinear effect that limits performance. We also implemented generic saturation utilities and basic dynamics-aware computations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, preparing the way for later courses where these constraints interact with control design and trajectory optimization.