Recall from previous chapters that the dynamics of an \( n \)-DOF rigid robot manipulator can be written as

\[ \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}, \quad \mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}\in\mathbb{R}^n. \]

Because inertial, gravitational and many friction parameters enter linearly, we can rewrite the inverse dynamics in the standard linear-in-parameters form

\[ \boldsymbol{\tau}(t) = \mathbf{Y}\!\big(\mathbf{q}(t),\dot{\mathbf{q}}(t),\ddot{\mathbf{q}}(t)\big) \,\boldsymbol{\theta}, \quad \boldsymbol{\theta}\in\mathbb{R}^p, \]

where \( \mathbf{Y} \) is the dynamic regressor and \( \boldsymbol{\theta} \) collects dynamic parameters (link masses, centers of mass, inertia components, friction coefficients, etc.). If we sample along a trajectory at times \( t_1,\dots,t_N \), we obtain stacked equations

\[ \mathbf{b} = \begin{bmatrix} \boldsymbol{\tau}(t_1)\\ \vdots\\ \boldsymbol{\tau}(t_N) \end{bmatrix} \in\mathbb{R}^{Nn}, \quad \boldsymbol{\Phi} = \begin{bmatrix} \mathbf{Y}(t_1)\\ \vdots\\ \mathbf{Y}(t_N) \end{bmatrix} \in\mathbb{R}^{Nn\times p}. \]

The excitation trajectory \( \mathbf{q}(t) \) determines \( \boldsymbol{\Phi} \). For parameter identification, what matters is \( \boldsymbol{\Phi} \): if it is poorly conditioned or rank-deficient, some parameters (or parameter combinations) cannot be reliably estimated.

flowchart TD
  Q["Trajectory q(t), qdot(t), qddot(t)"] --> Y["Regressor Y(q, qdot, qddot)"]
  Y --> PHI["Stacked matrix 'Phi' from samples"]
  PHI --> INFO["Information matrix 'Phi^T Phi'"]
  INFO --> QUAL["Parameter quality (rank, condition, variance)"]
  QUAL --> DESIGN["Adjust trajectory parameters and repeat if needed"]
        

Goal of excitation trajectory design: choose a feasible joint trajectory \( \mathbf{q}(t) \) that makes \( \boldsymbol{\Phi} \) as informative as possible, subject to physical and safety constraints (joint limits, velocity/acceleration bounds, actuator limits, collision avoidance, etc.).

We first relate the structure of \( \boldsymbol{\Phi} \) to identifiability of the parameters.

Proposition 1 (Uniqueness without noise).

Assume perfect measurements, so that for some unknown \( \boldsymbol{\theta}^{\star} \) we have

\[ \mathbf{b} = \boldsymbol{\Phi}\,\boldsymbol{\theta}^{\star}. \]

Then \( \boldsymbol{\theta}^{\star} \) is unique if and only if \( \boldsymbol{\Phi} \) has full column rank \( p \).

Proof.

(Only if): If \( \operatorname{rank}(\boldsymbol{\Phi}) < p \) then the null space \( \mathcal{N}(\boldsymbol{\Phi}) = \{\mathbf{z}\neq\mathbf{0} : \boldsymbol{\Phi}\mathbf{z}=\mathbf{0}\} \) is non-trivial. For any \( \mathbf{z}\in\mathcal{N}(\boldsymbol{\Phi}) \) and scalar \( \lambda \),

\[ \boldsymbol{\Phi}(\boldsymbol{\theta}^{\star} + \lambda\mathbf{z}) = \boldsymbol{\Phi}\boldsymbol{\theta}^{\star} + \lambda\boldsymbol{\Phi}\mathbf{z} = \mathbf{b}, \]

so there are infinitely many parameter vectors consistent with the data.

(If): If \( \operatorname{rank}(\boldsymbol{\Phi}) = p \) and \( \boldsymbol{\Phi}\boldsymbol{\theta}_1 = \boldsymbol{\Phi}\boldsymbol{\theta}_2 \) then

\[ \boldsymbol{\Phi}(\boldsymbol{\theta}_1-\boldsymbol{\theta}_2)=\mathbf{0} \;\Rightarrow\; \boldsymbol{\theta}_1-\boldsymbol{\theta}_2=\mathbf{0}, \]

because the only vector in the null space of a full-column-rank matrix is the zero vector. Hence \( \boldsymbol{\theta}_1=\boldsymbol{\theta}_2 \). \(\square\)

In practice, joint measurements and torques are noisy, and we will use least-squares estimation (next lesson). For now, recall the result from basic linear regression: if the measurement noise has covariance \( \sigma^2\mathbf{I} \), the least-squares estimate \( \hat{\boldsymbol{\theta}} \) has covariance

\[ \operatorname{Cov}(\hat{\boldsymbol{\theta}}) \approx \sigma^2 \big(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}\big)^{-1}. \]

Thus, \( \boldsymbol{\Phi}^{\top}\boldsymbol{\Phi} \) plays the role of an information matrix. Typical design criteria include:

• D-optimality: maximize \( \det(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}) \) (equivalently, minimize \( \det\big((\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi})^{-1}\big) \)).
• A-optimality: minimize \( \operatorname{tr}\big((\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi})^{-1}\big) \) (average variance of parameters).
• Condition number: minimize \( \kappa(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}) = \lambda_{\max}/\lambda_{\min} \) to avoid nearly unidentifiable parameter combinations.

The excitation trajectory design problem is thus an optimization over trajectories that seeks a good compromise between these criteria and robot constraints.

In continuous time, persistent excitation is a property of the regressor that guarantees that all parameter directions are sufficiently excited over any sliding time window. It is fundamental in adaptive control and on-line identification, but conceptually it is the same requirement we saw above: no parameter direction may remain “invisible”.

Definition (continuous-time PE).

The regressor \( \mathbf{Y}(t)\in\mathbb{R}^{n\times p} \) is persistently exciting of order \( p \) if there exist constants \( T > 0 \) and \( \alpha > 0 \) such that for all \( t_0 \)

\[ \int_{t_0}^{t_0+T} \mathbf{Y}(t)^{\top}\mathbf{Y}(t)\,\mathrm{d}t \succeq \alpha\,\mathbf{I}_p. \]

In discrete time, with samples \( \mathbf{Y}_k = \mathbf{Y}(t_k) \), we say the signal is PE if there exist \( N\in\mathbb{N} \) and \( \alpha > 0 \) such that for all \( k_0 \):

\[ \sum_{k=k_0}^{k_0+N-1} \mathbf{Y}_k^{\top}\mathbf{Y}_k \succeq \alpha\,\mathbf{I}_p. \]

Intuitively, this means that over any window of length \(N\) (or \(T\)), the Gram matrix of regressors has all eigenvalues bounded away from zero. This is a uniform version of the full-rank requirement.

Remark.

For linear parameter adaptation laws of the form \( \dot{\boldsymbol{\theta}} = -\Gamma\,\mathbf{Y}^{\top}e \) (with suitable error signal \( e \) and \( \Gamma \)), PE guarantees convergence of the parameter error to zero. While we do not study adaptive control here, the same condition explains why some trajectories (e.g., constant joint positions, or motions along a single joint) are fundamentally insufficient for identification.

flowchart TD
  A["Pick parametric form for q(t) (e.g. multi-sine)"] --> B["Simulate Y(t) for all joints"]
  B --> C["Compute sampled matrix 'Phi'"]
  C --> D["Evaluate 'Phi^T Phi' eigenvalues over windows"]
  D --> E{"Small eigenvalues?"}
  E -->|yes| TUNE["Adjust amplitudes/frequencies/phases"]
  E -->|no| ACCEPT["Trajectory is sufficiently exciting"]
        

In practice, checking PE exactly is impossible with finite data. Instead, we numerically approximate \( \boldsymbol{\Phi}^{\top}\boldsymbol{\Phi} \) from a simulated trajectory, and require:

• full column rank \( p \),
• minimum eigenvalue \( \lambda_{\min}(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}) \) above a threshold,
• condition number \( \kappa(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}) \) below a threshold.

We now describe trajectory families commonly used in robot dynamics identification. A typical choice is a finite Fourier/multi-sine parameterization for each joint \( i = 1,\dots,n \):

\[ q_i(t) = q_{i0} + \sum_{m=1}^M a_{i,m}\, \sin\!\big(\omega_{i,m} t + \varphi_{i,m}\big). \]

The corresponding velocity and acceleration are

\[ \dot{q}_i(t) = \sum_{m=1}^M a_{i,m}\omega_{i,m} \cos\!\big(\omega_{i,m} t + \varphi_{i,m}\big), \quad \ddot{q}_i(t) = -\sum_{m=1}^M a_{i,m}\omega_{i,m}^2 \sin\!\big(\omega_{i,m} t + \varphi_{i,m}\big). \]

We must respect bounds for each joint:

\[ q_i^{\min} \le q_i(t) \le q_i^{\max},\quad |\dot{q}_i(t)| \le \dot{q}_i^{\max},\quad |\ddot{q}_i(t)| \le \ddot{q}_i^{\max}. \]

Simple sufficient (though conservative) inequalities follow from triangle inequalities, e.g.

\[ \max_t |\dot{q}_i(t)| \le \sum_{m=1}^M |a_{i,m}\omega_{i,m}|, \quad \max_t |\ddot{q}_i(t)| \le \sum_{m=1}^M |a_{i,m}\omega_{i,m}^2|. \]

Let \( \mathbf{p} \) collect all amplitudes, frequencies and phases for all joints. A generic excitation trajectory design problem can then be written as:

\[ \begin{aligned} \min_{\mathbf{p}} \quad & J(\mathbf{p}) = -\log\det\big( \boldsymbol{\Phi}(\mathbf{p})^{\top}\boldsymbol{\Phi}(\mathbf{p})\big) \\ \text{s.t.}\quad & q_i^{\min} \le q_i(t_k,\mathbf{p}) \le q_i^{\max}, \\ & |\dot{q}_i(t_k,\mathbf{p})| \le \dot{q}_i^{\max}, \\ & |\ddot{q}_i(t_k,\mathbf{p})| \le \ddot{q}_i^{\max}, \quad \forall i,k, \end{aligned} \]

where \( t_k \) are discretization times. Other criteria such as minimizing the condition number or maximizing the minimum singular value of \( \boldsymbol{\Phi}(\mathbf{p}) \) are also common. Solving this nonconvex problem typically requires numerical optimization (gradient-based or evolutionary algorithms).

Alternative parameterizations used in the literature include B-splines, fifth-order polynomials patched to ensure continuity of \(q,\dot{q},\ddot{q}\), and PRBS-like inputs filtered through joint actuators. The mathematical treatment is similar: the regressor remains linear in parameters, and the information matrix depends on a finite vector of trajectory parameters.

We illustrate a minimal Python workflow using pinocchio for dynamics and numpy/scipy for linear algebra. The idea: (i) define a multi-sine trajectory for each joint, (ii) sample it, (iii) compute the joint torque regressor \( \mathbf{Y} \) at each sample via pinocchio.computeJointTorqueRegressor, and (iv) evaluate the condition number of \( \boldsymbol{\Phi}^{\top}\boldsymbol{\Phi} \).

import numpy as np
import pinocchio as pin

# 1) Load a robot model (URDF path and reference frame adapted to your setup)
model = pin.buildModelFromUrdf("path/to/your_robot.urdf")
data = model.createData()
nq = model.nq
nv = model.nv

# 2) Multi-sine parameterization for each joint
M = 3                      # number of harmonics per joint
T = 5.0                    # experiment duration [s]
dt = 0.002                 # sampling period [s]
t = np.arange(0.0, T, dt)
N = t.size

# Example parameters (to be later optimized)
q0   = np.zeros(nq)
A    = 0.4 * np.ones((nq, M))      # amplitudes
W    = np.linspace(1.0, 5.0, M)    # base frequencies [rad/s]
PHI0 = np.linspace(0.0, np.pi/2, M)

def multi_sine_joint_traj(t, q0, A_row, W, PHI0):
    """
    Single-joint multi-sine q(t), qdot(t), qddot(t).
    A_row: shape (M,)
    """
    q = np.full_like(t, q0, dtype=float)
    qd = np.zeros_like(t, dtype=float)
    qdd = np.zeros_like(t, dtype=float)
    for m in range(len(W)):
        w = W[m]
        a = A_row[m]
        ph = PHI0[m]
        s = np.sin(w * t + ph)
        c = np.cos(w * t + ph)
        q += a * s
        qd += a * w * c
        qdd += -a * (w ** 2) * s
    return q, qd, qdd

# 3) Build joint-space trajectories (q, qdot, qddot) for all joints
q_traj   = np.zeros((N, nq))
qd_traj  = np.zeros((N, nq))
qdd_traj = np.zeros((N, nq))

for i in range(nq):
    q_i, qd_i, qdd_i = multi_sine_joint_traj(t, q0[i], A[i, :], W, PHI0)
    q_traj[:, i]   = q_i
    qd_traj[:, i]  = qd_i
    qdd_traj[:, i] = qdd_i

# 4) For each sample, compute the dynamic regressor using Pinocchio
Y_list = []
for k in range(N):
    qk = q_traj[k, :]
    qdk = qd_traj[k, :]
    qddk = qdd_traj[k, :]
    # Ensure correct vector types (column vectors)
    qk = np.asarray(qk).reshape(-1, 1)
    qdk = np.asarray(qdk).reshape(-1, 1)
    qddk = np.asarray(qddk).reshape(-1, 1)

    # Compute regressor: Y_k so that tau = Y_k * theta
    Yk = pin.computeJointTorqueRegressor(model, data, qk, qdk, qddk)
    # Yk has shape (nv, p)
    Y_list.append(Yk)

# Stack into Phi
Phi = np.vstack(Y_list)   # shape (N*nv, p)

# 5) Information matrix and condition number
F = Phi.T @ Phi
eigs = np.linalg.eigvalsh(F)
lambda_min = np.min(eigs)
lambda_max = np.max(eigs)
cond = lambda_max / lambda_min
print("lambda_min =", lambda_min)
print("cond(F)    =", cond)
      

By wrapping the above in an outer optimization loop over the multi-sine parameters \( A,W,\Phi_0 \), one can minimize cond or maximize lambda_min to obtain a more exciting trajectory. This is precisely the numerical counterpart of the abstract criteria defined earlier.

The same workflow can be implemented in C++ using Pinocchio’s computeJointTorqueRegressor. Below is a sketch focusing on the trajectory and information-matrix aspects (error handling and URDF paths omitted for brevity).

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

#include <pinocchio/parsers/urdf.hpp>
#include <pinocchio/algorithm/regressor.hpp>

using Eigen::MatrixXd;
using Eigen::VectorXd;

struct MultiSineParams {
  MatrixXd A;      // (nq x M)
  Eigen::VectorXd W;
  Eigen::VectorXd PHI0;
  VectorXd q0;
};

void multiSineJoint(const Eigen::VectorXd &t,
                    double q0,
                    const Eigen::VectorXd &A_row,
                    const Eigen::VectorXd &W,
                    const Eigen::VectorXd &PHI0,
                    Eigen::VectorXd &q,
                    Eigen::VectorXd &qd,
                    Eigen::VectorXd &qdd)
{
  const int N = static_cast<int>(t.size());
  const int M = static_cast<int>(W.size());
  q  = VectorXd::Constant(N, q0);
  qd = VectorXd::Zero(N);
  qdd = VectorXd::Zero(N);
  for (int m = 0; m < M; ++m) {
    double w = W[m];
    double a = A_row[m];
    double ph = PHI0[m];
    for (int k = 0; k < N; ++k) {
      double arg = w * t[k] + ph;
      double s = std::sin(arg);
      double c = std::cos(arg);
      q[k]  += a * s;
      qd[k] += a * w * c;
      qdd[k] += -a * w * w * s;
    }
  }
}

int main()
{
  // 1) Load model
  pinocchio::Model model;
  pinocchio::urdf::buildModel("path/to/your_robot.urdf", model);
  pinocchio::Data data(model);
  const int nq = static_cast<int>(model.nq);
  const int nv = static_cast<int>(model.nv);

  // 2) Time grid
  double T = 5.0;
  double dt = 0.002;
  int N = static_cast<int>(T / dt);
  VectorXd t(N);
  for (int k = 0; k < N; ++k) t[k] = k * dt;

  // 3) Example parameters
  int M = 3;
  MultiSineParams P;
  P.A = MatrixXd::Constant(nq, M, 0.4);
  P.W = VectorXd::LinSpaced(M, 1.0, 5.0);
  P.PHI0 = VectorXd::LinSpaced(M, 0.0, M_PI / 2.0);
  P.q0 = VectorXd::Zero(nq);

  // 4) Build trajectories
  MatrixXd q_traj(N, nq), qd_traj(N, nq), qdd_traj(N, nq);
  for (int i = 0; i < nq; ++i) {
    Eigen::VectorXd q, qd, qdd;
    multiSineJoint(t, P.q0[i], P.A.row(i).transpose(), P.W, P.PHI0,
                   q, qd, qdd);
    q_traj.col(i) = q;
    qd_traj.col(i) = qd;
    qdd_traj.col(i) = qdd;
  }

  // 5) Stack regressor
  // p is dynamic parameter dimension (query from a first call)
  // Here we create Phi with unknown p: start with first sample
  VectorXd q0_vec = q_traj.row(0).transpose();
  VectorXd qd0_vec = qd_traj.row(0).transpose();
  VectorXd qdd0_vec = qdd_traj.row(0).transpose();
  MatrixXd Y0 = pinocchio::computeJointTorqueRegressor(
      model, data, q0_vec, qd0_vec, qdd0_vec);
  int p = static_cast<int>(Y0.cols());

  MatrixXd Phi(N * nv, p);
  Phi.block(0, 0, nv, p) = Y0;

  for (int k = 1; k < N; ++k) {
    VectorXd qk = q_traj.row(k).transpose();
    VectorXd qdk = qd_traj.row(k).transpose();
    VectorXd qddk = qdd_traj.row(k).transpose();
    MatrixXd Yk = pinocchio::computeJointTorqueRegressor(
        model, data, qk, qdk, qddk);
    Phi.block(k * nv, 0, nv, p) = Yk;
  }

  // 6) Information matrix and condition number
  MatrixXd F = Phi.transpose() * Phi;
  Eigen::SelfAdjointEigenSolver<MatrixXd> es(F);
  double lambda_min = es.eigenvalues().minCoeff();
  double lambda_max = es.eigenvalues().maxCoeff();
  double cond = lambda_max / lambda_min;
  std::cout << "lambda_min = " << lambda_min << std::endl;
  std::cout << "cond(F)    = " << cond << std::endl;

  return 0;
}
      

Such a C++ utility can be used offline to scan many candidate parameter sets and select those that provide numerically well-conditioned information matrices before running physical identification experiments.

Although robotics libraries for Java are less standardized, Java is convenient for trajectory generation and matrix computations, e.g. using EJML. Suppose we have a function computeRegressor(q, qd, qdd) that returns \( \mathbf{Y} \) for a planar 2-DOF arm (derived symbolically in earlier chapters). We can still study excitation quality in pure Java:

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;

public class ExcitationDesign {

    static double[][] multiSineJoint(double[] t, double q0,
                                     double[] A, double[] W, double[] PHI0) {
        int N = t.length;
        int M = W.length;
        double[] q = new double[N];
        double[] qd = new double[N];
        double[] qdd = new double[N];
        for (int k = 0; k < N; ++k) {
            q[k] = q0;
            qd[k] = 0.0;
            qdd[k] = 0.0;
        }
        for (int m = 0; m < M; ++m) {
            double w = W[m];
            double a = A[m];
            double ph = PHI0[m];
            for (int k = 0; k < N; ++k) {
                double arg = w * t[k] + ph;
                double s = Math.sin(arg);
                double c = Math.cos(arg);
                q[k]  += a * s;
                qd[k] += a * w * c;
                qdd[k] += -a * w * w * s;
            }
        }
        return new double[][]{q, qd, qdd};
    }

    // Placeholder for symbolic regressor of a 2R arm: Y(q, qd, qdd) in R^(2 x p)
    static DMatrixRMaj computeRegressor(double q1, double q2,
                                        double qd1, double qd2,
                                        double qdd1, double qdd2) {
        int p = 6; // example minimal parameter dimension
        DMatrixRMaj Y = new DMatrixRMaj(2, p);
        // Fill Y with appropriate basis functions of q, qd, qdd
        // (e.g., cos(q2), 2*qd1*qd2, qdd1, etc.)
        // ...
        return Y;
    }

    public static void main(String[] args) {
        double T = 5.0;
        double dt = 0.002;
        int N = (int) (T / dt);
        double[] t = new double[N];
        for (int k = 0; k < N; ++k) t[k] = k * dt;

        int M = 3;
        double[] A1 = {0.4, 0.3, 0.2};
        double[] W = {1.0, 2.5, 4.0};
        double[] PHI0 = {0.0, Math.PI / 3.0, Math.PI / 2.0};

        double[][] joint1 = multiSineJoint(t, 0.0, A1, W, PHI0);
        double[][] joint2 = multiSineJoint(t, 0.0, A1, W, PHI0); // copy for simplicity

        int p = 6;
        DMatrixRMaj Phi = new DMatrixRMaj(2 * N, p);
        DMatrixRMaj Yk;
        for (int k = 0; k < N; ++k) {
            double q1 = joint1[0][k], qd1 = joint1[1][k], qdd1 = joint1[2][k];
            double q2 = joint2[0][k], qd2 = joint2[1][k], qdd2 = joint2[2][k];
            Yk = computeRegressor(q1, q2, qd1, qd2, qdd1, qdd2);
            // copy into Phi at rows 2k and 2k+1
            for (int r = 0; r < 2; ++r) {
                for (int c = 0; c < p; ++c) {
                    Phi.set(2 * k + r, c, Yk.get(r, c));
                }
            }
        }

        // F = Phi^T Phi
        DMatrixRMaj F = new DMatrixRMaj(p, p);
        CommonOps_DDRM.multTransA(Phi, Phi, F);

        // eigenvalues (EJML has symmetric eigensolvers)
        org.ejml.dense.row.decomposition.eig.SymmetricQRAlgorithmDecomposition_DDRM eig =
                new org.ejml.dense.row.decomposition.eig.SymmetricQRAlgorithmDecomposition_DDRM(p, true);
        eig.decompose(F);
        double lambdaMin = Double.POSITIVE_INFINITY;
        double lambdaMax = 0.0;
        for (int i = 0; i < p; ++i) {
            double val = eig.getEigenvalue(i).getReal();
            lambdaMin = Math.min(lambdaMin, val);
            lambdaMax = Math.max(lambdaMax, val);
        }
        double cond = lambdaMax / lambdaMin;
        System.out.println("lambda_min = " + lambdaMin);
        System.out.println("cond(F)    = " + cond);
    }
}
      

This Java code emphasizes the structure of trajectory design and information-matrix evaluation; the symbolic expression for computeRegressor would come from previous analytical derivations for a specific arm.

In MATLAB, dynamic parameter identification is often carried out with Simulink models of the robot and a separate script that generates the excitation trajectory and passes it via workspace signals.

% Parameters
T  = 5.0;          % duration [s]
dt = 0.002;
t  = (0:dt:T).';
N  = numel(t);
nq = 6;            % 6-DOF arm
M  = 4;            % harmonics per joint

q0   = zeros(1, nq);
A    = 0.4 * ones(nq, M);
W    = linspace(1.0, 6.0, M);
PHI0 = linspace(0.0, pi/2, M);

q   = zeros(N, nq);
qd  = zeros(N, nq);
qdd = zeros(N, nq);

for i = 1:nq
    qi  = q0(i) * ones(N, 1);
    qdi = zeros(N, 1);
    qddi = zeros(N, 1);
    for m = 1:M
        w  = W(m);
        a  = A(i, m);
        ph = PHI0(m);
        s  = sin(w * t + ph);
        c  = cos(w * t + ph);
        qi  = qi  + a * s;
        qdi = qdi + a * w * c;
        qddi = qddi - a * w^2 * s;
    end
    q(:, i)   = qi;
    qd(:, i)  = qdi;
    qdd(:, i) = qddi;
end

% Export as Simulink timeseries for joint commands
excitation.q   = timeseries(q, t);
excitation.qd  = timeseries(qd, t);
excitation.qdd = timeseries(qdd, t);

% (In Simulink)
%  - Use "From Workspace" blocks to feed excitation.q to a joint trajectory
%    controller that tracks q(t).
%  - Measure joint torques tau_meas and joint positions.
%  - After simulation, build the regressor matrix Phi and evaluate
%    F = Phi' * Phi to check excitation quality.

% Example: compute information matrix from offline regressor function
Phi = [];
for k = 1:N
    qk   = q(k, :).';
    qdk  = qd(k, :).';
    qddk = qdd(k, :).';
    Yk = myRobotRegressor(qk, qdk, qddk);   % user-defined
    Phi = [Phi; Yk]; %#ok<AGROW>
end

F = Phi' * Phi;
lambda = eig(F);
lambda_min = min(lambda);
condF = max(lambda) / lambda_min;
fprintf('lambda_min = %.3e, cond(F) = %.3e\n', lambda_min, condF);
      

The symbolic function myRobotRegressor can be generated using MATLAB’s Symbolic Toolbox from the Lagrange or Newton–Euler derivations, or numerically via existing robotics toolboxes.

Mathematica is well suited for symbolic manipulations of both the robot dynamics and the excitation trajectory. For instance, we can derive the Gram matrix analytically for a simplified 1-DOF pendulum and optimize the Fourier coefficients under constraints.

(* Time grid *)
T  = 5.0;
dt = 0.01;
tList = Range[0.0, T, dt];

(* Single-joint Fourier trajectory *)
M = 3;
a = Array[a, M];
w = {1.0, 2.0, 3.0};
phi0 = {0.0, Pi/3, 2 Pi/3};
q0 = 0.0;

q[t_] := q0 + Sum[a[[m]]*Sin[w[[m]]*t + phi0[[m]]], {m, 1, M}];
qd[t_] := D[q[t], t];
qdd[t_] := D[q[t], {t, 2}];

(* Example scalar regressor for a 1-DOF pendulum:
   Y(t) = [qdd(t), Cos[q(t)]] corresponding to [I, m*g*l] *)
Y[t_] := {qdd[t], Cos[q[t]]};

(* Sampled Phi *)
Phi = Table[Y[t], {t, tList}];
PhiMat = N[Phi]; (* numeric matrix, rows are samples *)

(* Information matrix and determinant as a function of a1,a2,a3 *)
F[a1_?NumericQ, a2_?NumericQ, a3_?NumericQ] := Module[
  {PhiNum, FF},
  PhiNum = PhiMat /. {a[1] -> a1, a[2] -> a2, a[3] -> a3};
  FF = Transpose[PhiNum].PhiNum;
  Det[FF]
];

(* Simple optimization: maximize log det(F) with amplitude bounds *)
sol = NMaximize[
  {Log[F[a1, a2, a3]],
   -0.8 <= a1 <= 0.8 && -0.8 <= a2 <= 0.8 && -0.8 <= a3 <= 0.8},
  {a1, a2, a3}
];

sol
      

This symbolic-numeric mix can be extended to multi-DOF manipulators, though the expressions for \( \mathbf{Y} \) grow quickly in size. Still, Mathematica provides an excellent environment to prototype and analyze theoretical excitation criteria.

In this lesson we connected the linear-in-parameters representation of robot dynamics with the design of informative excitation trajectories. The key object is the stacked regressor \( \boldsymbol{\Phi} \), whose rank and conditioning determine which parameter combinations can be identified and with what accuracy. Persistent excitation provides a mathematical condition that guarantees that all parameter directions are sufficiently excited over time. We formulated trajectory design as a constrained optimization problem over parametric families (multi-sine, Fourier series, B-splines) and showed how to evaluate excitation quality via the information matrix \( \boldsymbol{\Phi}^{\top}\boldsymbol{\Phi} \).

Through Python, C++, Java, MATLAB/Simulink and Mathematica examples, we emphasized that, regardless of language or specific library, the core design loop is the same: parameterize trajectories, simulate or measure the regressor, compute information metrics, and iterate until the excitation is adequate and all physical constraints are satisfied. In the next lesson, we will formally develop least-squares methods that exploit these trajectories to estimate the dynamic parameters.