For an \( n \)-DOF serial manipulator with generalized coordinates \( \mathbf{q} \in \mathbb{R}^n \), the joint-space dynamics (neglecting friction and external contacts) have the generic form

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

where \( \mathbf{M}(\mathbf{q}) \) is the symmetric positive-definite inertia matrix, \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) captures centrifugal and Coriolis effects, and \( \mathbf{g}(\mathbf{q}) \) is the gravity vector.

From Chapters 10–12 you know how to derive these matrices using Lagrange–Euler or Newton–Euler methods. In this lesson we view the dynamics as a linear regression in a vector of inertial parameters:

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

where \( \mathbf{Y} \) is the regressor matrix and \( \boldsymbol{\pi} \) is a stacked vector of scalar inertial parameters (mass, products of inertia, first moments, etc.).

We distinguish:

• Full parameter set: a parameter vector \( \boldsymbol{\pi} \in \mathbb{R}^p \) that contains all standard link inertial parameters (e.g. \(10\) per link), constructed so that the dynamics are linear in \( \boldsymbol{\pi} \).
• Minimal (base) parameter set: a shortest-length parameter vector \( \boldsymbol{\beta} \in \mathbb{R}^r \), with \( r \leq p \), such that the same dynamics can be written as \( \boldsymbol{\tau} = \mathbf{Y}_b(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) \boldsymbol{\beta} \) with \( \mathbf{Y}_b \) having full column rank \( r \).

Minimal parameterization is crucial in identification (Chapter 19), model-based control, and numerical robustness, because redundant parameters lead to poor conditioning and non-uniqueness.

flowchart TD
  A["Physical robot links (m, com, inertia)"]
    --> B["Full parameter vector pi (10*n scalars)"]
  B --> C["Dynamics: tau = Y(q,qd,qdd) * pi"]
  C --> D["Detect dependencies in Y"]
  D --> E["Transform to base parameters beta"]
  E --> F["Dynamics: tau = Yb(q,qd,qdd) * beta (full rank)"]
        

Consider link \( i \) with a body-fixed frame \( \{O_i\} \). Let:

• \( m_i \): mass of link \( i \).
• \( \mathbf{c}_i = [c_{ix}, c_{iy}, c_{iz}]^\top \): position of the center of mass expressed in frame \( \{O_i\} \).
• \( \bar{\mathbf{I}}_i \): \(3\times 3\) inertia tensor about \( O_i \) expressed in \( \{O_i\} \), with components \( I_{ixx}, I_{iyy}, I_{izz}, I_{ixy}, I_{ixz}, I_{iyz} \) .

The spatial inertia of link \( i \) in frame \( \{O_i\} \) is

\[ \mathbf{I}_i = \begin{bmatrix} \bar{\mathbf{I}}_i & m_i[\mathbf{c}_i]_{\times} \\ m_i[\mathbf{c}_i]_{\times}^\top & m_i\mathbf{I}_3 \end{bmatrix}, \]

where \( [\mathbf{c}_i]_{\times} \) is the skew-symmetric cross-product matrix and \( \mathbf{I}_3 \) is the \(3\times 3\) identity matrix. All elements of \( \mathbf{I}_i \) are linear functions of the following standard inertial parameter vector for link \( i \):

\[ \boldsymbol{\pi}_i = \begin{bmatrix} m_i \\ m_i c_{ix} \\ m_i c_{iy} \\ m_i c_{iz} \\ I_{ixx} \\ I_{iyy} \\ I_{izz} \\ I_{ixy} \\ I_{ixz} \\ I_{iyz} \end{bmatrix} \in \mathbb{R}^{10}. \]

For an \( n \)-link serial chain, the full parameter vector is obtained by stacking all link parameters:

\[ \boldsymbol{\pi} = \begin{bmatrix} \boldsymbol{\pi}_1 \\ \boldsymbol{\pi}_2 \\ \vdots \\ \boldsymbol{\pi}_n \end{bmatrix} \in \mathbb{R}^{10n}. \]

Because each spatial inertia matrix and hence the kinetic and potential energies are linear in \( \boldsymbol{\pi} \), the resulting equations of motion are linear in \( \boldsymbol{\pi} \) as well. This is the foundation for writing the dynamics as \( \boldsymbol{\tau} = \mathbf{Y}\boldsymbol{\pi} \).

Theorem. For any rigid serial manipulator with known kinematics and spatial inertias \( \mathbf{I}_i \) parameterized linearly by \( \boldsymbol{\pi} \), the joint dynamics can be written in the form \( \boldsymbol{\tau}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) = \mathbf{Y}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) \boldsymbol{\pi} \).

Sketch of proof.

1. For link \( i \) with spatial twist \( \mathbf{v}_i \in \mathbb{R}^6 \), the kinetic energy is

   \[ T_i = \tfrac{1}{2}\mathbf{v}_i^\top \mathbf{I}_i \mathbf{v}_i. \]

   Because \( \mathbf{I}_i \) is linear in \( \boldsymbol{\pi}_i \), we can write \( T_i = \boldsymbol{\phi}_i(\mathbf{q},\dot{\mathbf{q}})^\top \boldsymbol{\pi}_i \) for some vector of scalar functions \( \boldsymbol{\phi}_i \).
2. The total kinetic energy is

   \[ T(\mathbf{q},\dot{\mathbf{q}}) = \sum_{i=1}^n T_i = \boldsymbol{\phi}(\mathbf{q},\dot{\mathbf{q}})^\top \boldsymbol{\pi}, \]

   where \( \boldsymbol{\phi} = [\boldsymbol{\phi}_1^\top,\dots, \boldsymbol{\phi}_n^\top]^\top \) .
3. Gravitational potential energy for link \( i \) is

   \[ U_i(\mathbf{q}) = m_i \mathbf{g}_0^\top \mathbf{p}_{C_i}(\mathbf{q}) = \boldsymbol{\psi}_i(\mathbf{q})^\top \boldsymbol{\pi}_i, \]

   which is linear in \( \boldsymbol{\pi}_i \) because \( m_i \) and \( m_i\mathbf{c}_i \) appear linearly.
4. Summing over links gives

   \[ U(\mathbf{q}) = \boldsymbol{\psi}(\mathbf{q})^\top \boldsymbol{\pi}, \]

   with \( \boldsymbol{\psi} = [\boldsymbol{\psi}_1^\top,\dots,\boldsymbol{\psi}_n^\top]^\top \).
5. The Lagrangian is

   \[ \mathcal{L}(\mathbf{q},\dot{\mathbf{q}}) = T(\mathbf{q},\dot{\mathbf{q}}) - U(\mathbf{q}) = \boldsymbol{\ell}(\mathbf{q},\dot{\mathbf{q}})^\top \boldsymbol{\pi}, \]

   with \( \boldsymbol{\ell} = \boldsymbol{\phi} - \boldsymbol{\psi} \) , still linear in \( \boldsymbol{\pi} \).
6. Lagrange's equations

   \[ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}}{\partial \dot{q}_j} \right) - \frac{\partial \mathcal{L}}{\partial q_j} = \tau_j ,\quad j=1,\dots,n, \]

   involve linear operators (partial derivatives and time derivative) acting on a function that is linear in \( \boldsymbol{\pi} \). Linearity is preserved, so each joint equation is a linear combination of the entries of \( \boldsymbol{\pi} \).

Stacking the \( n \) joint equations yields the regressor form with \( \mathbf{Y}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) \) of size \( n\times 10n \).

In practice, many columns of the regressor \( \mathbf{Y} \) are linearly dependent because of kinematic constraints, choice of link frames, and dynamic symmetries. Consequently, not all components of \( \boldsymbol{\pi} \) affect the joint torques independently.

For some finite collection of samples \( \{(\mathbf{q}_k,\dot{\mathbf{q}}_k,\ddot{\mathbf{q}}_k)\}_{k=1}^N \) , define the stacked regressor

\[ \mathbf{W} = \begin{bmatrix} \mathbf{Y}(\mathbf{q}_1,\dot{\mathbf{q}}_1,\ddot{\mathbf{q}}_1) \\ \mathbf{Y}(\mathbf{q}_2,\dot{\mathbf{q}}_2,\ddot{\mathbf{q}}_2) \\ \vdots \\ \mathbf{Y}(\mathbf{q}_N,\dot{\mathbf{q}}_N,\ddot{\mathbf{q}}_N) \end{bmatrix} \in \mathbb{R}^{(Nn)\times p}, \]

where \( p=10n \). Ideally, the rank of \( \mathbf{W} \) for a generic persistently exciting data set equals the structural rank of the dynamics.

Define the null space of the regressor:

\[ \mathcal{N} = \left\{ \Delta\boldsymbol{\pi} \in \mathbb{R}^p \,\middle|\, \mathbf{Y}(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) \Delta\boldsymbol{\pi} = \mathbf{0} \text{ for all admissible } (\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}) \right\}. \]

Any two parameter vectors \( \boldsymbol{\pi}_1,\boldsymbol{\pi}_2 \) with \( \boldsymbol{\pi}_1 - \boldsymbol{\pi}_2 \in \mathcal{N} \) generate identical joint torques for all trajectories and are dynamically indistinguishable.

Let \( r = \dim(\mathrm{im}(\mathbf{Y})) \) denote the rank and \( p-r = \dim(\mathcal{N}) \) the dimension of the null space. We can choose a matrix \( \mathbf{V}_1 \in \mathbb{R}^{p\times r} \) whose columns form a basis for the orthogonal complement of the null space. Then the base parameter vector is

\[ \boldsymbol{\beta} = \mathbf{V}_1^\top \boldsymbol{\pi} \in \mathbb{R}^r, \]

and there exists a reduced regressor \( \mathbf{Y}_b \) such that

\[ \boldsymbol{\tau} = \mathbf{Y} \boldsymbol{\pi} = \mathbf{Y} \mathbf{V}_1 \boldsymbol{\beta} = \mathbf{Y}_b \boldsymbol{\beta}, \quad \mathrm{rank}(\mathbf{Y}_b)=r. \]

By construction, \( \boldsymbol{\beta} \) is a minimal parameter set: no component can be removed without losing the ability to reproduce all possible joint torques.

Numerically, we approximate the structural rank by computing the singular value decomposition (SVD) of the stacked regressor \( \mathbf{W} \). Let

\[ \mathbf{W} = \mathbf{U} \begin{bmatrix} \boldsymbol{\Sigma}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \mathbf{V}^\top = \mathbf{U} \begin{bmatrix} \boldsymbol{\Sigma}_1 & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{V}_1^\top \\ \mathbf{V}_2^\top \end{bmatrix}, \]

where \( \boldsymbol{\Sigma}_1 \) contains the singular values above a numerical tolerance, and \( \mathbf{V} = [\mathbf{V}_1\ \mathbf{V}_2] \) is orthogonal.

• Columns of \( \mathbf{V}_2 \) span the numerical null space, i.e., combinations of parameters that do not affect joint torques.
• Columns of \( \mathbf{V}_1 \) define the base parameters \( \boldsymbol{\beta} = \mathbf{V}_1^\top \boldsymbol{\pi} \).

In practice we:

1. Generate a sufficiently rich set of configurations, velocities, and accelerations.
2. Compute the regressor \( \mathbf{Y} \) for each sample and stack into \( \mathbf{W} \).
3. Compute the SVD and inspect singular values to determine numerical rank.
4. Form \( \mathbf{V}_1 \) from the columns corresponding to nonzero singular values.

flowchart TD
  S["Sample trajectories (q, qd, qdd)"]
    --> Ymat["Compute regressors Y_k"]
  Ymat --> Wmat["Stack into big matrix W"]
  Wmat --> SVD["SVD: W = U * Sigma * V^T"]
  SVD --> RANK["Choose rank r by singular values"]
  RANK --> V1["Extract V1 (first r columns of V)"]
  V1 --> BETA["Define base params beta = V1^T * pi"]
  BETA --> YB["Reduced regressor Yb = Y * V1"]
        

Consider a single-link planar pendulum (revolute joint at the base) with generalized coordinate \( q \), link length \( L \), COM distance \( c \) from the joint, mass \( m \), and inertia about the joint \( I_O \). Take gravity pointing downward and measure \( q \) from the downward vertical.

The kinetic and potential energies are

\[ T = \tfrac{1}{2} \left( I_O + m c^2 \right) \dot{q}^2, \quad U = m g c (1 - \cos q). \]

The Lagrangian \( \mathcal{L} = T-U \) yields the equation of motion

\[ \left( I_O + m c^2 \right)\ddot{q} + m g c \sin q = \tau. \]

The physical parameters are \( I_O \), \( m \), and \( c \). However, the dynamics depend on them only through the combinations

\[ \theta_1 = I_O + m c^2, \qquad \theta_2 = m c. \]

We can rewrite the dynamics as

\[ \tau = \underbrace{ \begin{bmatrix} \ddot{q} & g \sin q \end{bmatrix} }_{\mathbf{Y}(q,\dot{q},\ddot{q})} \underbrace{ \begin{bmatrix} \theta_1 \\ \theta_2 \end{bmatrix} }_{\boldsymbol{\beta}}. \]

Thus the base parameter vector is \( \boldsymbol{\beta} = [\theta_1,\theta_2]^\top \) of dimension 2, while the full parameter vector \( \boldsymbol{\pi}_{\text{phys}} = [I_O, m, c]^\top \) has dimension 3. The mapping

\[ \boldsymbol{\beta} = \begin{bmatrix} 1 & c^2 & 0 \\ 0 & c & m \end{bmatrix} \boldsymbol{\pi}_{\text{phys}} \]

shows explicitly that only two independent combinations of \( I_O, m, c \) affect the joint torque. Any choice of physical parameters that keeps \( \theta_1,\theta_2 \) fixed produces identical dynamics.

In Python we can use numpy and scipy or sympy together with robotics libraries such as roboticstoolbox (Corke) or pinocchio bindings to construct regressors. Below we create a synthetic regressor for the 1-DOF pendulum example and recover the 2D base parameter space from a 4D over-parameterized vector.

Suppose we choose the full parameter vector \( \boldsymbol{\pi} = [I_O,\; m,\; m c^2,\; m c]^\top \) so that \( I_O \), \( m \), \( m c^2 \), and \( m c \) appear linearly. The regressor row for one sample is then \( [\ddot{q},\; 0,\; \ddot{q},\; g\sin q] \).

import numpy as np

g = 9.81

def pendulum_regressor_row(q, qdd):
    """
    Regressor row for the over-parameterized pendulum model:
    tau = [qdd, 0, qdd, g*sin(q)] @ [I0, m, m*c**2, m*c]
    """
    return np.array([qdd, 0.0, qdd, g * np.sin(q)], dtype=float)

# Generate random samples of (q, qdd)
rng = np.random.default_rng(0)
N = 50
qs = rng.uniform(low=-np.pi, high=np.pi, size=N)
qdds = rng.uniform(low=-5.0, high=5.0, size=N)

# Stack into big regressor W (size N x 4)
W = np.vstack([pendulum_regressor_row(q, qdd) for q, qdd in zip(qs, qdds)])

# Compute SVD of W
U, S, Vt = np.linalg.svd(W, full_matrices=False)

print("Singular values:", S)

# Numerical rank: count singular values above a tolerance
tol = 1e-8
r = np.sum(S > tol)
print("Numerical rank r =", r)

V1 = Vt[:r, :].T   # first r right-singular vectors (columns)
V2 = Vt[r:, :].T   # remaining p-r columns (null space basis)

print("V1 (basis of base parameter space):\n", V1)
print("V2 (basis of null space):\n", V2)

# Transform a sample full parameter vector pi into base parameters beta:
pi_full = np.array([0.1, 2.0, 2.0 * 0.25, 2.0 * 0.5])  # I0, m, m*c**2, m*c
beta = V1.T @ pi_full
print("Base parameter coordinates beta:", beta)
      

The numerical rank is \( r=2 \), confirming that only two independent combinations of the four parameters matter. In a more realistic multi-DOF manipulator implementation, we would replace pendulum_regressor_row by calls to a robotics toolbox that returns the regressor for given joint states.

In C++, rigid-body dynamics libraries such as RBDL or Pinocchio provide efficient algorithms for forward dynamics and inverse dynamics. Although not all expose a regressor directly, you can construct it by evaluating joint torques for carefully chosen parameter vectors, or by using symbolic tools offline.

The following snippet shows how to compute the SVD of a stacked regressor using Eigen and transform a full parameter vector into base parameters:

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

int main() {
    using Eigen::MatrixXd;
    using Eigen::VectorXd;

    const int N = 50;   // number of samples
    const int p = 4;    // number of full parameters in this example

    MatrixXd W(N, p);

    // TODO: fill W with regressor rows from your dynamics library.
    // Here we simply fake something with random values to illustrate:
    W.setRandom();

    // Compute thin SVD: W = U * S * V^T
    Eigen::JacobiSVD<MatrixXd> svd(W, Eigen::ComputeThinU | Eigen::ComputeThinV);
    VectorXd S = svd.singularValues();
    MatrixXd V = svd.matrixV();

    double tol = 1e-8;
    int r = 0;
    for (int i = 0; i < S.size(); ++i) {
        if (S(i) > tol) { ++r; }
    }
    std::cout << "Rank r = " << r << std::endl;

    MatrixXd V1 = V.leftCols(r);         // base parameter subspace
    MatrixXd V2 = V.rightCols(p - r);    // null-space basis

    // Example full parameter vector
    VectorXd pi_full(p);
    pi_full << 0.1, 2.0, 0.5, 1.0;

    // Base parameters beta = V1^T * pi
    VectorXd beta = V1.transpose() * pi_full;

    std::cout << "beta =\n" << beta << std::endl;
    std::cout << "Null-space basis V2 =\n" << V2 << std::endl;

    return 0;
}
      

For a real manipulator, W would be filled by calling, for example, RBDL::InverseDynamics on unit parameter variations or by an explicit regressor routine if available. The SVD step is identical.

In Java, popular linear algebra libraries like EJML or Apache Commons Math support SVD and matrix operations. Robot-specific functionality can be layered on top (e.g., your own Newton–Euler implementation or bindings to C++ libraries).

import org.ejml.simple.SimpleMatrix;

public class BaseParameterDemo {
    public static void main(String[] args) {
        int N = 50;  // samples
        int p = 4;   // full parameters

        // Build a random regressor W (replace with real robot regressor)
        SimpleMatrix W = SimpleMatrix.random_DDRM(N, p, -1.0, 1.0, null);

        // Compute SVD: W = U * S * V^T
        var svd = W.svd();
        SimpleMatrix V = svd.getV();
        double[] s = svd.getW().extractDiag().getDDRM().getData();

        // Determine numerical rank
        double tol = 1e-8;
        int r = 0;
        for (double value : s) {
            if (value > tol) {
                r++;
            }
        }
        System.out.println("Rank r = " + r);

        // Extract V1 (first r columns) and V2 (remaining columns)
        SimpleMatrix V1 = V.extractMatrix(0, SimpleMatrix.END, 0, r);
        SimpleMatrix V2 = V.extractMatrix(0, SimpleMatrix.END, r, p);

        // Example full parameter vector
        SimpleMatrix piFull = new SimpleMatrix(p, 1, true,
                                               0.1, 2.0, 0.5, 1.0);
        // Base parameters beta = V1^T * pi
        SimpleMatrix beta = V1.transpose().mult(piFull);

        System.out.println("beta = ");
        beta.print();
        System.out.println("Null-space basis V2 = ");
        V2.print();
    }
}
      

The core pattern is identical to the Python and C++ examples: build a stacked regressor W, compute its SVD, and obtain base parameters as coordinates in the column space of V1.

MATLAB, together with the Robotics System Toolbox and Symbolic Math Toolbox, is a natural environment for deriving and reducing inertial parameter sets. Simulink can then be used to simulate dynamics using the minimal parameter model.

The following script constructs the pendulum regressor, computes its SVD, and forms base parameters:

g = 9.81;

% Regressor row for over-parameterized pendulum
% tau = [qdd, 0, qdd, g*sin(q)] * [I0; m; m*c^2; m*c]
pendulumRegressorRow = @(q, qdd) [qdd, 0, qdd, g * sin(q)];

% Generate samples
N  = 50;
qs   = (rand(N,1) - 0.5) * 2*pi;
qdds = (rand(N,1) - 0.5) * 10;

W = zeros(N, 4);
for k = 1:N
    W(k,:) = pendulumRegressorRow(qs(k), qdds(k));
end

% SVD
[U, S, V] = svd(W, "econ");
singularVals = diag(S);
tol = 1e-8;
r = sum(singularVals > tol);

V1 = V(:,1:r);
V2 = V(:,r+1:end);

disp("Rank r = " + r);
disp("V1 (base parameter basis):");
disp(V1);
disp("V2 (null-space basis):");
disp(V2);

% Example full parameter vector
pi_full = [0.1; 2.0; 0.5; 1.0];
beta = V1.' * pi_full;
disp("Base parameters beta:");
disp(beta);
      

In Simulink, you can implement two variants of the pendulum model: one using the full parameter vector and one using the base parameters \( \boldsymbol{\beta} \), and verify that the simulated joint torques match for all trajectories.

Mathematica is very convenient for deriving symbolic regressors and separating parameter combinations. The next code derives the pendulum equation of motion and extracts its regressor with respect to an over-parameterized vector \( [I_O,\; m,\; m c^2,\; m c]^\top \).

Clear[q, t, I0, m, c, g];

q[t_] := q[t];

(* Kinetic and potential energy *)
T = 1/2 (I0 + m c^2) (D[q[t], t])^2;
U = m g c (1 - Cos[q[t]]);
L = T - U;

(* Lagrange equation *)
eom = D[D[L, D[q[t], t]], t] - D[L, q[t]] // Simplify;

(* Replace derivatives by symbols qdd, qd, and q *)
eomSub = eom /. {q[t] -> q, D[q[t], t] -> qd, D[q[t], {t, 2}] -> qdd};

(* Over-parameterized vector *)
params = {I0, m, m c^2, m c};

(* Collect eom w.r.t parameters and extract coefficient matrix *)
eomCollected = Collect[eomSub, params];
{Yrow, piVec} = CoefficientArrays[eomCollected, params];

Yrow // MatrixForm
piVec // MatrixForm
      

The resulting row Yrow coincides with \( [\ddot{q},\; 0,\; \ddot{q},\; g\sin q] \), and the symbolic machinery can be extended to multi-DOF manipulators. You can then export the regressor to other languages or use Mathematica's own linear algebra tools to compute base parameters.

In this lesson we introduced the notion of full and minimal inertial parameter sets for rigid manipulators. Starting from spatial inertias, we showed that robot dynamics are linear in a standard 10-parameter description per link and can be written as \( \boldsymbol{\tau} = \mathbf{Y}\boldsymbol{\pi} \). We then used linear algebra to characterize the null space of the regressor, which captures dynamically indistinguishable parameter variations, and defined the base parameter vector \( \boldsymbol{\beta} \) as a minimal set of independent combinations. An SVD-based algorithm provides a practical way to compute base parameters numerically. Finally, we illustrated the concepts with a 1-DOF pendulum and implemented the reduction procedure in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.