In the previous lessons of this chapter, we described how to estimate dynamic parameters of a manipulator using least squares on the regressor form of the equations of motion. Recall that for each sampled time instant \( k = 1,\dots,N \), the joint torque vector satisfies the linear-in-parameters structure

\[ \boldsymbol{\tau}_k = \mathbf{Y}_k \boldsymbol{\pi} + \boldsymbol{\varepsilon}_k, \quad \boldsymbol{\tau}_k \in \mathbb{R}^m, \ \mathbf{Y}_k \in \mathbb{R}^{m \times p}, \ \boldsymbol{\pi} \in \mathbb{R}^p, \]

where \( m \) is the number of joints and \( p \) the number of dynamic parameters (masses, inertia components, friction coefficients, etc.). Collecting all samples into stacked vectors and matrices gives

\[ \boldsymbol{\tau} = \mathbf{Y} \boldsymbol{\pi} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\tau} \in \mathbb{R}^{Nm},\; \mathbf{Y} \in \mathbb{R}^{Nm \times p}. \]

Least squares estimation yields \( \hat{\boldsymbol{\pi}} \). Model validation then asks: Is the pair \( (\mathbf{Y},\hat{\boldsymbol{\pi}}) \) an adequate description of the robot's dynamics for new trajectories? To answer this, we define validation metrics based on prediction errors (residuals) computed on a validation data set that was not used for estimation.

Let \( \hat{\boldsymbol{\tau}}_k \) denote the torque predicted by the identified model for sample \(k\). The basic residual is

\[ \mathbf{e}_k = \boldsymbol{\tau}_k^{\text{meas}} - \hat{\boldsymbol{\tau}}_k, \quad k = 1,\dots,N. \]

Model validation metrics are functionals \( J(\{\mathbf{e}_k\}, \{\boldsymbol{\tau}_k^{\text{meas}}\}) \) that summarize residual magnitude, normalization, and statistical structure. They should be:

• Scale-aware: Do not unfairly favor models just because torques are small.
• Statistically interpretable: Connect to noise assumptions (e.g. Gaussian).
• Comparable: Allow comparing different model structures (e.g. with or without friction).

flowchart TD
  D["Experimental data: q, qdot, qddot, tau_meas"] --> ID["Parameter estimation (e.g. least squares)"]
  ID --> M["Identified model: tau_hat = Y * pi_hat"]
  M --> VAL["Apply model on validation trajectories"]
  VAL --> RES["Compute residuals: e = tau_meas - tau_hat"]
  RES --> MET["Compute metrics: RMS, VAF, R2, AIC/BIC, tests"]
  MET --> DEC["Decision: accept model, refine, or redesign experiment"]
        

For simplicity, consider a single joint \( m=1 \) and denote scalar torques \( \tau_k^{\text{meas}} \) and \( \hat{\tau}_k \). The residual sequence is

\[ e_k = \tau_k^{\text{meas}} - \hat{\tau}_k, \quad k = 1,\dots,N. \]

Fundamental norm-based error measures include:

• Sum of squared errors (SSE):

  \[ \mathrm{SSE} = \sum_{k=1}^N e_k^2. \]

  This is the quantity minimized by ordinary least squares.
• Root-mean-square error (RMSE):

  \[ \mathrm{RMSE} = \sqrt{\frac{1}{N} \sum_{k=1}^N e_k^2 }. \]

  RMSE has the same units as torque and is widely used in robotics.
• Maximum absolute error:

  \[ \mathrm{MAE}_\infty = \max_{1 \leq k \leq N} |e_k|. \]

  This highlights worst-case mismatch between model and data.

For multi-joint manipulators, residuals become vectors \( \mathbf{e}_k \in \mathbb{R}^m \) and we typically use Euclidean norms:

\[ \mathrm{RMSE}_{\text{multi}} = \sqrt{\frac{1}{N m} \sum_{k=1}^N \left\| \mathbf{e}_k \right\|_2^2 }. \]

A crucial limitation of raw RMSE is its dependence on the scale of the torques; a robot whose joints routinely produce hundreds of N·m will naturally exhibit larger RMSE than a light planar arm, even for equally good relative fits. This motivates normalized metrics.

Normalized metrics compare the residual variance to the variance of the measured signal. They are invariant to uniform scaling of torques and enable fair comparisons across joints and robots.

3.1 Variance Accounted For (VAF)

For a scalar sequence, define sample variance of measured torques and residuals as

\[ \operatorname{var}(\tau^{\text{meas}}) = \frac{1}{N-1} \sum_{k=1}^N \left(\tau_k^{\text{meas}} - \bar{\tau}\right)^2, \quad \bar{\tau} = \frac{1}{N}\sum_{k=1}^N \tau_k^{\text{meas}}, \]

\[ \operatorname{var}(e) = \frac{1}{N-1} \sum_{k=1}^N \left(e_k - \bar{e}\right)^2, \quad \bar{e} = \frac{1}{N}\sum_{k=1}^N e_k. \]

The VAF is then

\[ \mathrm{VAF} = 100 \left( 1 - \frac{\operatorname{var}(e)} {\operatorname{var}(\tau^{\text{meas}})} \right) \% . \]

If the model is perfect up to measurement noise, \( \operatorname{var}(e) \) is small and VAF approaches \(100\%\). Poor models may yield negative VAF, which indicates that predicting a constant mean torque \( \bar{\tau} \) is better than using the model.

3.2 MATLAB-style Fit Percentage

A widely used metric (mirroring the one in MATLAB's System Identification Toolbox) is the fit percentage:

\[ \mathrm{Fit} = 100 \left( 1 - \frac{ \left\| \tau^{\text{meas}} - \hat{\tau} \right\|_2 }{ \left\| \tau^{\text{meas}} - \bar{\tau}\mathbf{1} \right\|_2 } \right) \% , \]

where \( \tau^{\text{meas}} \) and \( \hat{\tau} \) are length-\(N\) vectors and \( \mathbf{1} \) is the all-ones vector. This measure compares the model to a trivial constant model. When least squares is used with a constant term (intercept), one can show that \( \mathrm{Fit} \) is closely related to the classical coefficient of determination \( R^2 \).

3.3 Coefficient of Determination \( R^2 \)

Define the total sum of squares (TSS) and residual sum of squares (RSS):

\[ \mathrm{TSS} = \sum_{k=1}^N \left(\tau_k^{\text{meas}} - \bar{\tau}\right)^2, \quad \mathrm{RSS} = \sum_{k=1}^N e_k^2. \]

Then

\[ R^2 = 1 - \frac{\mathrm{RSS}}{\mathrm{TSS}}. \]

When least squares is used and the regression includes a constant term, the identity \( \mathrm{TSS} = \mathrm{SSR} + \mathrm{RSS} \) holds, where SSR is the regression sum of squares. Since both SSR and RSS are non-negative, it follows that \( 0 \leq R^2 \leq 1 \).

Sketch of proof: Let \( \mathbf{y} = \tau^{\text{meas}} \), \( \hat{\mathbf{y}} = \hat{\tau} \), and \( \mathbf{1} \) the all-ones vector. The least-squares projection of \( \mathbf{y} \) onto the column space of the regressor matrix (which contains \( \mathbf{1} \)) decomposes as

\[ \mathbf{y} - \bar{\tau}\mathbf{1} = (\hat{\mathbf{y}} - \bar{\tau}\mathbf{1}) + (\mathbf{y} - \hat{\mathbf{y}}), \]

where the two right-hand terms are orthogonal. Taking squared norms yields \( \mathrm{TSS} = \mathrm{SSR} + \mathrm{RSS} \), and therefore \( R^2 = \mathrm{SSR}/\mathrm{TSS} \in [0,1] \).

Good normalized fit does not guarantee that the model captures all relevant dynamics. Residuals might still show temporal structure (e.g. periodic patterns) indicating unmodeled inertia, friction, or flexibility. Statistical residual analysis therefore complements scalar metrics.

4.1 Residual Mean and Variance

A basic check is that residuals have zero mean with variance consistent with measurement noise:

\[ \hat{\mu}_e = \frac{1}{N} \sum_{k=1}^N e_k, \quad \hat{\sigma}_e^2 = \frac{1}{N-1} \sum_{k=1}^N (e_k - \hat{\mu}_e)^2. \]

A significantly non-zero \( \hat{\mu}_e \) indicates systematic bias in the model (e.g. missing constant friction).

4.2 Autocorrelation and Whiteness

For an adequate model with independent measurement noise, residuals should be approximately white: uncorrelated over time. The sample autocovariance at lag \( \ell \) is

\[ \hat{R}_e(\ell) = \frac{1}{N-\ell} \sum_{k=\ell+1}^N (e_k - \hat{\mu}_e) (e_{k-\ell} - \hat{\mu}_e), \quad \ell = 0,1,\dots,L. \]

The normalized autocorrelation is

\[ \hat{\rho}_e(\ell) = \frac{\hat{R}_e(\ell)}{\hat{R}_e(0)}. \]

Under the hypothesis of white noise residuals with large \(N\), it can be shown that, approximately, \( \hat{\rho}_e(\ell) \sim \mathcal{N}(0,1/N) \) for \( \ell \geq 1 \). Thus a heuristic \(95\%\) confidence band is

\[ |\hat{\rho}_e(\ell)| \lesssim \frac{1.96}{\sqrt{N}}. \]

Large values outside this band at specific lags often indicate unmodeled dynamics (for example, missing elastic components leading to oscillatory structure in \( e_k \)).

4.3 Correlation with Regressors

A more powerful test considers correlations between residuals and regressor columns. Let \( \mathbf{Y}_k = [Y_{k,1},\dots,Y_{k,p}] \) be the regressor row at time \(k\). For each parameter \(j\) and lag \( \ell \geq 0 \), one can compute

\[ \hat{R}_{eY_j}(\ell) = \frac{1}{N-\ell} \sum_{k=\ell+1}^N (e_k - \hat{\mu}_e) (Y_{k-\ell,j} - \bar{Y}_j), \]

\[ \hat{\rho}_{eY_j}(\ell) = \frac{\hat{R}_{eY_j}(\ell)} {\sqrt{\hat{R}_e(0)\hat{R}_{Y_j}(0)}}, \]

where \( \hat{R}_{Y_j}(0) \) is the sample variance of regressor component \(Y_{k,j}\). For an adequate model, these correlations should be statistically indistinguishable from zero. Significant correlation at lag \( \ell \) suggests that information in \( Y_{k-\ell,j} \) has not been fully exploited by the model, hinting at missing dynamics or nonlinearities.

flowchart TD
  R["Residuals e(k)"] --> M0["Check mean(e) approx 0"]
  R --> AC["Compute autocorrelation rho_e(l)"]
  R --> RC["Compute cross-corr with \nregressors rho_eY(l)"]
  M0 --> D0["Bias → \nadd offsets/friction terms"]
  AC --> D1["Oscillatory structure → \nmissing flexible dynamics"]
  RC --> D2["Significant corr with Y → \nmissing nonlinear/structured terms"]
        

In robot dynamics identification, several model structures are often considered: with or without Coulomb friction, with reduced inertia parameter sets, etc. Metrics like RMSE or VAF alone tend to favor more complex models (larger \(p\)). Information criteria such as AIC and BIC address this by penalizing complexity.

Suppose we have \( N \) samples of an \(m\)-dimensional torque vector and a model with \( p \) estimated parameters. Let

\[ \hat{\sigma}_e^2 = \frac{1}{Nm} \sum_{k=1}^N \left\| \mathbf{e}_k \right\|_2^2 \]

be the empirical noise variance. Under the standard Gaussian noise assumption, the negative log-likelihood is (up to constants)

\[ -\log L(\hat{\boldsymbol{\pi}}) = \frac{Nm}{2} \log \hat{\sigma}_e^2 + \text{const}. \]

The Akaike information criterion (AIC) and Bayesian information criterion (BIC) are then defined as

\[ \mathrm{AIC} = 2p + Nm \log \hat{\sigma}_e^2, \quad \mathrm{BIC} = p \log(Nm) + Nm \log \hat{\sigma}_e^2. \]

Between two candidate models, the one with smaller AIC or BIC is preferred. BIC penalizes complexity more strongly than AIC, and is therefore more conservative in adding parameters (e.g. extra friction terms).

In practice, AIC/BIC are used together with residual analysis: a model with slightly higher AIC but clearly uncorrelated, white residuals may be preferable to a model with marginally lower AIC but structured residuals.

We now implement the above metrics in Python using numpy and scipy. Assume that we already have validation data and a function that evaluates the regressor \( \mathbf{Y}_k \) for given joint states and accelerations, as introduced in earlier lessons of this chapter.

import numpy as np
from dataclasses import dataclass

@dataclass
class ValidationMetrics:
    rmse: float
    vaf: float
    fit_percent: float
    r2: float
    aic: float
    bic: float

def compute_residuals(tau_meas: np.ndarray,
                      tau_pred: np.ndarray) -> np.ndarray:
    """
    tau_meas: shape (N, m)
    tau_pred: shape (N, m)
    returns residuals e = tau_meas - tau_pred, shape (N, m)
    """
    return tau_meas - tau_pred

def compute_basic_metrics(tau_meas: np.ndarray,
                          tau_pred: np.ndarray) -> ValidationMetrics:
    tau_meas = np.asarray(tau_meas)
    tau_pred = np.asarray(tau_pred)
    assert tau_meas.shape == tau_pred.shape
    N, m = tau_meas.shape

    e = tau_meas - tau_pred

    # RMSE over all joints and samples
    rmse = np.sqrt(np.mean(np.sum(e**2, axis=1)))

    # Flatten for convenience
    y = tau_meas.reshape(-1)
    y_hat = tau_pred.reshape(-1)
    e_flat = e.reshape(-1)

    # Means and variance
    y_mean = np.mean(y)
    e_mean = np.mean(e_flat)
    var_y = np.var(y, ddof=1)
    var_e = np.var(e_flat, ddof=1)

    # VAF (in percent)
    vaf = 100.0 * (1.0 - var_e / var_y)

    # Fit percentage
    num = np.linalg.norm(y - y_hat)
    den = np.linalg.norm(y - y_mean)
    fit_percent = 100.0 * (1.0 - num / den)

    # R^2
    tss = np.sum((y - y_mean)**2)
    rss = np.sum((y - y_hat)**2)
    r2 = 1.0 - rss / tss

    # AIC and BIC
    # Assume effective number of parameters p is known:
    p = 30  # example; use actual number of estimated dynamic parameters
    sigma2_hat = np.mean(e_flat**2)
    Nm = N * m
    aic = 2.0 * p + Nm * np.log(sigma2_hat)
    bic = p * np.log(Nm) + Nm * np.log(sigma2_hat)

    return ValidationMetrics(rmse, vaf, fit_percent, r2, aic, bic)

def residual_autocorrelation(e: np.ndarray, L: int) -> np.ndarray:
    """
    e: residuals, shape (N,) for a single joint
    L: maximum lag
    returns normalized autocorrelation rho_e(l) for l = 0..L
    """
    e = np.asarray(e).reshape(-1)
    N = e.size
    e_mean = np.mean(e)
    e0 = e - e_mean
    R0 = np.dot(e0, e0) / N
    rho = np.zeros(L + 1)
    rho[0] = 1.0
    for ell in range(1, L + 1):
        num = np.dot(e0[ell:], e0[:-ell]) / (N - ell)
        rho[ell] = num / R0
    return rho

# Example usage:
# tau_meas_val, tau_pred_val: arrays of shape (N, m) on validation trajectories
# metrics = compute_basic_metrics(tau_meas_val, tau_pred_val)
# rho_joint1 = residual_autocorrelation((tau_meas_val - tau_pred_val)[:, 0], L=50)
      

In a complete identification workflow, the same code can be applied to several candidate models and validation trajectories, selecting the one that best balances low RMSE, high VAF/fit, low AIC/BIC, and white residuals.

In C++, we can use the Eigen library for linear algebra. The code below assumes that torque data are stored as matrices Eigen::MatrixXd of size N x m.

#include <Eigen/Dense>
#include <cmath>
#include <cstddef>

struct ValidationMetricsCpp {
    double rmse;
    double vaf;
    double fit_percent;
    double r2;
    double aic;
    double bic;
};

ValidationMetricsCpp computeMetricsCpp(const Eigen::MatrixXd& tau_meas,
                                       const Eigen::MatrixXd& tau_pred,
                                       std::size_t p_params)
{
    using Eigen::MatrixXd;
    using Eigen::VectorXd;

    const std::ptrdiff_t N = tau_meas.rows();
    const std::ptrdiff_t m = tau_meas.cols();
    MatrixXd e = tau_meas - tau_pred;

    // RMSE (flattened)
    double rmse = std::sqrt((e.array().square().rowwise().sum()).mean());

    // Flatten to a single vector
    VectorXd y = Eigen::Map<const VectorXd>(tau_meas.data(), N * m);
    VectorXd y_hat = Eigen::Map<const VectorXd>(tau_pred.data(), N * m);
    VectorXd e_flat = y - y_hat;

    double y_mean = y.mean();
    double e_mean = e_flat.mean();

    VectorXd y0 = y.array() - y_mean;
    VectorXd e0 = e_flat.array() - e_mean;

    double var_y = (y0.array().square().sum()) / (y.size() - 1);
    double var_e = (e0.array().square().sum()) / (e_flat.size() - 1);

    double vaf = 100.0 * (1.0 - var_e / var_y);

    double num = (y - y_hat).norm();
    double den = (y.array() - y_mean).matrix().norm();
    double fit_percent = 100.0 * (1.0 - num / den);

    double tss = (y.array() - y_mean).square().sum();
    double rss = (y - y_hat).squaredNorm();
    double r2 = 1.0 - rss / tss;

    double sigma2_hat = e_flat.squaredNorm() / static_cast<double>(N * m);
    double Nm = static_cast<double>(N * m);
    double aic = 2.0 * static_cast<double>(p_params) + Nm * std::log(sigma2_hat);
    double bic = static_cast<double>(p_params) * std::log(Nm) + Nm * std::log(sigma2_hat);

    ValidationMetricsCpp out { rmse, vaf, fit_percent, r2, aic, bic };
    return out;
}
      

This function can be integrated into a robot dynamics identification pipeline written in C++, where tau_pred is obtained by evaluating the regressor Y_k and parameter vector pi_hat for each sample.

In Java, we can implement similar computations using plain arrays or matrix libraries such as EJML. Here we illustrate a simple static utility for scalar-joint data.

public final class ValidationMetricsJava {

    public static class Metrics {
        public final double rmse;
        public final double vaf;
        public final double fitPercent;
        public final double r2;
        public final double aic;
        public final double bic;

        public Metrics(double rmse, double vaf, double fitPercent,
                       double r2, double aic, double bic) {
            this.rmse = rmse;
            this.vaf = vaf;
            this.fitPercent = fitPercent;
            this.r2 = r2;
            this.aic = aic;
            this.bic = bic;
        }
    }

    public static Metrics compute(double[] tauMeas, double[] tauPred, int pParams) {
        int N = tauMeas.length;
        if (N != tauPred.length) {
            throw new IllegalArgumentException("Array sizes must match");
        }

        double[] e = new double[N];
        double sumESq = 0.0;
        double sumTau = 0.0;

        for (int k = 0; k < N; ++k) {
            e[k] = tauMeas[k] - tauPred[k];
            sumESq += e[k] * e[k];
            sumTau += tauMeas[k];
        }

        double rmse = Math.sqrt(sumESq / N);
        double tauMean = sumTau / N;

        double tss = 0.0;
        for (int k = 0; k < N; ++k) {
            double d = tauMeas[k] - tauMean;
            tss += d * d;
        }

        // Variances
        double eMean = 0.0;
        for (int k = 0; k < N; ++k) {
            eMean += e[k];
        }
        eMean /= N;

        double varY = 0.0;
        double varE = 0.0;
        for (int k = 0; k < N; ++k) {
            double dy = tauMeas[k] - tauMean;
            double de = e[k] - eMean;
            varY += dy * dy;
            varE += de * de;
        }
        varY /= (N - 1);
        varE /= (N - 1);

        double vaf = 100.0 * (1.0 - varE / varY);

        double num = 0.0;
        double den = 0.0;
        for (int k = 0; k < N; ++k) {
            double d1 = tauMeas[k] - tauPred[k];
            double d2 = tauMeas[k] - tauMean;
            num += d1 * d1;
            den += d2 * d2;
        }
        num = Math.sqrt(num);
        den = Math.sqrt(den);
        double fitPercent = 100.0 * (1.0 - num / den);

        double rss = sumESq;
        double r2 = 1.0 - rss / tss;

        double sigma2Hat = rss / N;
        double aic = 2.0 * pParams + N * Math.log(sigma2Hat);
        double bic = pParams * Math.log(N) + N * Math.log(sigma2Hat);

        return new Metrics(rmse, vaf, fitPercent, r2, aic, bic);
    }

    private ValidationMetricsJava() {
        // utility
    }
}
      

For multi-joint validation, call this method per joint or extend it to two-dimensional arrays.

MATLAB is widely used for robot dynamics modeling. After estimating parameters, we can simulate the identified model in Simulink and compare the simulated torques with measured torques logged from the robot.

function metrics = validate_dynamics_model(tau_meas, tau_pred, p_params)
% tau_meas, tau_pred: N-by-m matrices
% p_params: number of estimated dynamic parameters

[N, m] = size(tau_meas);
e = tau_meas - tau_pred;

% RMSE over all joints
rmse = sqrt(mean(sum(e.^2, 2)));

y = tau_meas(:);
y_hat = tau_pred(:);
e_flat = y - y_hat;

y_mean = mean(y);
e_mean = mean(e_flat);

var_y = var(y, 1);   % population form (1/N)
var_e = var(e_flat, 1);

vaf = 100 * (1 - var_e / var_y);

num = norm(y - y_hat);
den = norm(y - y_mean);
fit_percent = 100 * (1 - num / den);

tss = sum((y - y_mean).^2);
rss = sum((y - y_hat).^2);
r2 = 1 - rss / tss;

sigma2_hat = mean(e_flat.^2);
Nm = N * m;
aic = 2 * p_params + Nm * log(sigma2_hat);
bic = p_params * log(Nm) + Nm * log(sigma2_hat);

metrics = struct('rmse', rmse, ...
                 'vaf', vaf, ...
                 'fit_percent', fit_percent, ...
                 'r2', r2, ...
                 'aic', aic, ...
                 'bic', bic);
end
      

In a Simulink workflow, one would typically:

1. Create a Simulink model implementing the rigid-body dynamics using \( \hat{\boldsymbol{\pi}} \).
2. Feed measured joint trajectories \( q, \dot{q}, \ddot{q} \) as inputs.
3. Use To Workspace blocks to log simulated torques tau_pred.
4. Call validate_dynamics_model with the logged data and measured torques.

Wolfram Mathematica offers convenient symbolic and numeric tools for validation experiments. The snippet below assumes scalar joint data.

(* tauMeas and tauPred are lists of equal length *)
validateDynamicsModel[tauMeas_List, tauPred_List, pParams_Integer] := Module[
  {
    e, n, rmse, y, yMean, eMean, varY, varE, vaf,
    num, den, fit, tss, rss, r2, sigma2Hat, aic, bic
  },
  n = Length[tauMeas];
  If[n != Length[tauPred],
    Return[$Failed, Module]
  ];
  e = tauMeas - tauPred;

  rmse = Sqrt[Mean[e^2]];

  y = tauMeas;
  yMean = Mean[y];
  eMean = Mean[e];

  varY = Variance[y, WorkingPrecision -> MachinePrecision];
  varE = Variance[e, WorkingPrecision -> MachinePrecision];

  vaf = 100.0 (1.0 - varE/varY);

  num = Norm[tauMeas - tauPred];
  den = Norm[tauMeas - ConstantArray[yMean, n]];
  fit = 100.0 (1.0 - num/den);

  tss = Total[(y - yMean)^2];
  rss = Total[e^2];
  r2 = 1.0 - rss/tss;

  sigma2Hat = Mean[e^2];
  aic = 2.0 pParams + n Log[sigma2Hat];
  bic = pParams Log[n] + n Log[sigma2Hat];

  <<Association[
    "RMSE" -> rmse,
    "VAF" -> vaf,
    "FitPercent" -> fit,
    "R2" -> r2,
    "AIC" -> aic,
    "BIC" -> bic
  ]
];

(* Example of residual autocorrelation up to lag L *)
residualAutocorrelation[e_List, L_Integer] := Module[{n, e0, r0, rho},
  n = Length[e];
  e0 = e - Mean[e];
  r0 = Total[e0^2]/n;
  rho = Table[
    If[l == 0, 1.0,
      Total[Take[e0, {l + 1, n}] * Take[e0, {1, n - l}]] / ((n - l) r0)
    ],
    {l, 0, L}
  ];
  rho
];
      

Mathematica's plotting capabilities (e.g. ListPlot) can be used to visualize residual autocorrelation and confidence bands.

This lesson formalized the role of validation metrics in robot dynamics identification. Starting from residuals between measured and predicted torques, we defined norm-based errors (RMSE, maximum error), normalized fit measures (VAF, fit percentage, \( R^2 \)), and statistical residual diagnostics (means, variances, autocorrelation, and cross-correlation with regressors). We then introduced AIC and BIC to balance model fit against complexity, and implemented these quantities in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. Together, these tools allow a control engineer to systematically decide whether an identified dynamic model is adequate for subsequent tasks such as simulation, planning, and (in later chapters) control design.