Chapter 15: Learning-Augmented Control (Non-RL Focus)

Lesson 5: Lab: Learning Residual Dynamics for Better Tracking

This lab implements a learning-augmented controller for a robot joint (or a low-DOF manipulator) by identifying a residual dynamics model on top of an existing model-based controller (e.g., computed-torque). We go from mathematical formulation to a complete simulation and controller code in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

1. Lab Objectives and Setup

We consider a robot manipulator with known nominal rigid-body dynamics and unknown modeling errors (unmodeled friction, flexibilities, actuation imperfections). The true joint-space dynamics are

\[ M(q)\ddot q + C(q,\dot q)\dot q + g(q) + d(q,\dot q,\ddot q) = \tau, \]

where \( q \in \mathbb{R}^n \), \( M(q) \) is the inertia matrix, \( C(q,\dot q)\dot q \) the Coriolis/centrifugal term, \( g(q) \) gravity, and \( d(\cdot) \) is an unknown disturbance / modeling error. We only have a nominal model

\[ M_n(q)\ddot q + C_n(q,\dot q)\dot q + g_n(q) \approx \tau. \]

The goal of this lab is to:

  • Implement a nominal model-based tracking controller.
  • Collect data under this controller and estimate a residual map \( r^\star(z) \) that approximates the missing dynamics, where \( z \) is a feature vector built from \( (q,\dot q,\ddot q_d) \).
  • Close the loop with a learning-augmented controller \( \tau = \tau_{\mathrm{nom}} + \hat r_\theta(z) \) and compare tracking performance.
  • Check (at least qualitatively) that the learned residual preserves the stability guarantees from the nominal design (by bounding the residual).

The overall experimental flow is illustrated below.

flowchart TD
  T["Desired trajectory q_d(t)"] --> N["Nominal controller tau_nom(q, qdot, qd, qd_dot, qd_ddot)"]
  N --> R["Robot (true dynamics)"]
  R --> S["Sensors (q, qdot, tau)"]
  S --> L["Log data and build dataset (z_k, y_k)"]
  L --> F["Fit residual model hat r_theta(z)"]
  F --> C["Updated controller tau = tau_nom + hat r_theta(z)"]
  C --> R
        

We will use a low-DOF case (1-DOF or 2-DOF) in code, but the derivations are valid for general manipulators.

2. Nominal Computed-Torque Tracking and Error Dynamics

Let the desired joint trajectory be \( q_d(t) \), with derivatives \( \dot q_d(t), \ddot q_d(t) \). Define tracking error and its derivative as

\[ e = q - q_d, \quad \dot e = \dot q - \dot q_d. \]

The standard computed-torque controller built on the nominal model is

\[ \tau_{\mathrm{nom}}(t) = M_n(q)\bigl(\ddot q_d - K_d \dot e - K_p e\bigr) + C_n(q,\dot q)\dot q + g_n(q), \]

where \( K_p, K_d \) are positive-definite gain matrices. Substituting \( \tau = \tau_{\mathrm{nom}} \) into the true dynamics, we obtain

\[ \begin{aligned} M(q)\ddot q + C(q,\dot q)\dot q + g(q) + d(q,\dot q,\ddot q) &= M_n(q)\bigl(\ddot q_d - K_d \dot e - K_p e\bigr) + C_n(q,\dot q)\dot q + g_n(q). \end{aligned} \]

Rearranging and expressing in terms of the tracking error dynamics, we get

\[ \begin{aligned} M(q)\ddot e + C(q,\dot q)\dot e + K_d M_n(q)\dot e + K_p M_n(q) e = \Delta(q,\dot q,\ddot q_d), \end{aligned} \]

where the lumped disturbance \( \Delta \) collects both the physical disturbance \( d \) and the nominal model mismatch \( M_n - M, C_n - C, g_n - g \). Under mild assumptions and sufficiently large gains, the nominal controller renders the origin \( e=\dot e=0 \) input-to-state stable (ISS) with respect to \( \Delta \). Thus, if \( \|\Delta\| \) is small, the tracking error remains small; but if modeling errors are large, performance is limited.

3. Residual Dynamics Model and Closed-Loop Structure

We introduce a residual torque map \( r^\star(z) \in \mathbb{R}^n \) that approximates the unknown lumped disturbance \( \Delta \). The feature vector \( z \) is built from the measured state and desired motion:

\[ z = \phi(q, \dot q, \ddot q_d) = \begin{bmatrix} q \\ \dot q \\ \ddot q_d \end{bmatrix} \in \mathbb{R}^{3n}. \]

During learning, we estimate a parameterized model \( \hat r_\theta(z) \) and define the learning-augmented control law

\[ \tau(t) = \tau_{\mathrm{nom}}(t) + \hat r_\theta\bigl(z(t)\bigr). \]

A reasonable target for the residual is the torque error:

\[ y_k = \tau_k - \tau_{\mathrm{nom}}(z_k) \approx r^\star(z_k) \]

obtained from experiments with the baseline controller. After fitting \( \hat r_\theta \) on the dataset \( \{(z_k,y_k)\}_{k=1}^N \), we use it online to augment the control.

To reason about stability, let \( \tilde r(z) = r^\star(z) - \hat r_\theta(z) \). The error dynamics become

\[ \begin{aligned} M(q)\ddot e + C(q,\dot q)\dot e + K_d M_n(q)\dot e + K_p M_n(q) e = \tilde r(z). \end{aligned} \]

Using the Lyapunov candidate

\[ V(e,\dot e) = \tfrac{1}{2}\dot e^\top M(q)\dot e + \tfrac{1}{2} e^\top K_p e, \]

one obtains (using standard manipulator properties and the skew-symmetry of \( \dot M - 2C \)) an inequality of the form

\[ \dot V \leq -\alpha_1 \|\dot e\|^2 - \alpha_2 \|e\|^2 + \beta \|\dot e\| \,\|\tilde r(z)\|, \]

for some positive constants \( \alpha_1, \alpha_2, \beta \). If we enforce a bound on the residual error \( \|\tilde r(z)\| \leq \bar r \) for all reachable \( z \) (e.g., by gain scheduling, filtering, or saturation of \( \hat r_\theta \)), then the origin is ISS with respect to \( \bar r \). This motivates conservative residual usage in the lab.

4. Data Collection Protocol

We now describe the procedure to generate training data for the residual model in simulation (and analogously on real hardware). We assume sensor access to \( q(t), \dot q(t) \) and the applied torque \( \tau(t) \).

  1. Choose a rich set of reference trajectories \( q_d^{(j)}(t) \), for \( j = 1,\dots,J \), that sufficiently excite the dynamics (e.g., sinusoids with different frequencies, chirps).
  2. For each trajectory, run the nominal controller \( \tau_{\mathrm{nom}}(t) \) without residuals and log \( \{q(t_k), \dot q(t_k), \ddot q_d(t_k), \tau(t_k)\} \) at sampling times \( t_k \).
  3. For each sample, construct \( z_k = \phi(q(t_k),\dot q(t_k),\ddot q_d(t_k)) \) and target \( y_k = \tau(t_k) - \tau_{\mathrm{nom}}(z_k) \).
  4. Aggregate the dataset \( \mathcal{D} = \{(z_k,y_k)\}_{k=1}^N \) and use it to fit a regression model \( \hat r_\theta \) (Gaussian Process, kernel ridge regression, small neural network, etc.).
  5. Validate on held-out trajectories and quantify prediction error versus torque magnitude to choose safe bounds for online compensation.
flowchart TD
  A["Run nominal controller"] --> B["Log (q, qdot, qd_ddot, tau)"]
  B --> C["Build features z_k and targets y_k"]
  C --> D["Fit regressor hat r_theta(z)"]
  D --> E["Validate and bound residual error"]
  E --> F["Deploy in closed loop with saturation"]
        

In the remainder of the lab, we will instantiate this procedure first in Python, then sketch implementations in C++, Java, MATLAB/Simulink, and Mathematica.

5. Python Implementation: 1-DOF Residual Learning Lab

We consider a 1-DOF rotary joint with true dynamics

\[ m_{\mathrm{true}} \ddot q + b_{\mathrm{true}} \dot q + g_{\mathrm{true}}\sin(q) + d(q,\dot q) = \tau, \]

and a nominal model that ignores nonlinear friction and uses slightly wrong parameters. We implement:

  • Nominal computed-torque tracking.
  • Data collection and Gaussian Process regression of the residual.
  • Closed-loop evaluation with the learned residual.

import numpy as np

# Optional: scikit-learn for Gaussian Processes
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# True parameters
m_true = 1.2
b_true = 0.3
g_true = 9.81

# Nominal parameters (slightly wrong)
m_nom = 1.0
b_nom = 0.1
g_nom = 9.81

# Unknown disturbance term (unmodeled dynamics)
def disturbance(q, qd):
    # Smooth but nonlinear "extra" friction and flex effect
    return 0.4 * np.sin(2.0 * q) + 0.05 * qd**3

def true_dynamics(q, qd, tau):
    d = disturbance(q, qd)
    qdd = (tau - b_true * qd - g_true * np.sin(q) - d) / m_true
    return qdd, d

def nominal_torque(q, qd, qd_des, qdd_des, Kp, Kd):
    e = q - qd_des
    edot = qd - (qd_des)  # qd_des is dot(q_d); passed in explicitly
    v = qdd_des - Kd * edot - Kp * e
    tau_nom = m_nom * v + b_nom * qd + g_nom * np.sin(q)
    return tau_nom

# Desired trajectory (moderately rich)
def desired_trajectory(t):
    # q_d(t), qd_d(t), qdd_d(t)
    A = 0.7
    w = 0.7
    qd = A * np.sin(w * t)
    qd_dot = A * w * np.cos(w * t)
    qd_ddot = -A * w**2 * np.sin(w * t)
    return qd, qd_dot, qd_ddot

# Simulation parameters
dt = 0.002
T = 12.0
N = int(T / dt)

Kp = 25.0
Kd = 10.0

def rollout(with_residual=False, gpr=None, r_max=2.0, collect_data=False):
    q = 0.0
    qd = 0.0
    qs = []
    qds = []
    qd_des_list = []
    qdd_des_list = []
    taus = []
    Z = []
    Y = []

    for k in range(N):
        t = k * dt
        qd_des, qd_dot_des, qdd_des = desired_trajectory(t)
        tau_nom = nominal_torque(q, qd, qd_des, qdd_des, Kp, Kd)

        z = np.array([q, qd, qd_des, qdd_des], dtype=float)

        if with_residual and gpr is not None:
            # Predict residual and saturate
            r_hat = gpr.predict(z.reshape(1, -1))[0]
            r_hat = float(np.clip(r_hat, -r_max, r_max))
        else:
            r_hat = 0.0

        tau = tau_nom + r_hat

        # True dynamics integration (semi-implicit Euler)
        qdd, d_true = true_dynamics(q, qd, tau)
        qd = qd + dt * qdd
        q = q + dt * qd

        # Store logs
        qs.append(q)
        qds.append(qd)
        qd_des_list.append(qd_des)
        qdd_des_list.append(qdd_des)
        taus.append(tau)

        if collect_data:
            y = tau - tau_nom   # target residual
            Z.append(z)
            Y.append(y)

    result = {
        "q": np.array(qs),
        "qd": np.array(qds),
        "q_d": np.array(qd_des_list),
        "qd_d": np.gradient(np.array(qd_des_list), dt),
        "qdd_d": np.array(qdd_des_list),
        "tau": np.array(taus),
    }
    if collect_data:
        result["Z"] = np.array(Z)
        result["Y"] = np.array(Y)
    return result

# 1) Collect data under nominal controller
nominal_run = rollout(with_residual=False, collect_data=True)
Z_train = nominal_run["Z"]
Y_train = nominal_run["Y"]

# 2) Fit a scalar-output GP residual model
kernel = 1.0 * RBF(length_scale=0.5) + WhiteKernel(noise_level=1e-4)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=True)
gpr.fit(Z_train, Y_train)

print("Fitted GP kernel:", gpr.kernel_)

# 3) Evaluate tracking with and without residual
baseline = rollout(with_residual=False, collect_data=False)
augmented = rollout(with_residual=True, gpr=gpr, r_max=2.0, collect_data=False)

# 4) Compute RMS tracking error
def rms(x):
    return np.sqrt(np.mean(x**2))

err_base = baseline["q"] - baseline["q_d"]
err_aug = augmented["q"] - augmented["q_d"]

print("RMS error baseline :", rms(err_base))
print("RMS error augmented:", rms(err_aug))
      

In a full lab, the students should:

  • Plot \( q(t) \) versus \( q_d(t) \) for both controllers.
  • Plot the learned residual \( \hat r_\theta(z_k) \) against the true lumped disturbance (which is known in simulation).
  • Experiment with different kernels, data lengths, and saturation levels \( r_{\max} \) to observe tracking vs safety trade-offs.

6. C++ Implementation Sketch with Eigen

In C++, we typically rely on a robotics dynamics library (e.g., RBDL, Pinocchio) to compute the nominal torque and use Eigen for linear algebra. Here is a minimal sketch for a 1-DOF joint using an RBF (Radial Basis Function) residual model implemented from scratch.


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

struct RBFResidual {
    // phi(z) = exp(-gamma * ||z - c_j||^2) for j = 1..M
    std::vector<Eigen::VectorXd> centers;
    Eigen::VectorXd weights;
    double gamma;

    double predict(const Eigen::VectorXd& z) const {
        const int M = static_cast<int>(centers.size());
        Eigen::VectorXd phi(M);
        for (int j = 0; j < M; ++j) {
            double r2 = (z - centers[j]).squaredNorm();
            phi(j) = std::exp(-gamma * r2);
        }
        return weights.dot(phi);
    }
};

double nominal_torque(double q, double qd, double qd_des,
                      double qdd_des, double Kp, double Kd,
                      double m_nom, double b_nom, double g_nom)
{
    double e = q - qd_des;
    double edot = qd - qd_des; // qd_des is dot(q_d)
    double v = qdd_des - Kd * edot - Kp * e;
    double tau_nom = m_nom * v + b_nom * qd + g_nom * std::sin(q);
    return tau_nom;
}

double saturate(double x, double r_max) {
    if (x > r_max) return r_max;
    if (x < -r_max) return -r_max;
    return x;
}

// In the control loop (pseudo-code):
// 1. Build feature vector z = [q, qd, qd_des, qdd_des]^T
// 2. Compute tau_nom
// 3. Compute residual r_hat = saturate(rbf.predict(z), r_max)
// 4. Apply tau = tau_nom + r_hat
void controller_step(double q, double qd,
                     double qd_des, double qdd_des,
                     const RBFResidual& rbf,
                     double Kp, double Kd,
                     double m_nom, double b_nom, double g_nom,
                     double r_max,
                     double& tau_out)
{
    Eigen::VectorXd z(4);
    z << q, qd, qd_des, qdd_des;

    double tau_nom = nominal_torque(q, qd, qd_des, qdd_des,
                                    Kp, Kd, m_nom, b_nom, g_nom);
    double r_hat = saturate(rbf.predict(z), r_max);
    tau_out = tau_nom + r_hat;
}
      

In a full implementation, the centers and weights of RBFResidual are learned offline by solving a regularized least-squares problem using the dataset \( \{(z_k,y_k)\} \) collected under the nominal controller.

7. Java Implementation Sketch with EJML

Java is less common for low-level robot control, but libraries like EJML provide matrix support. Below is a conceptual sketch using a simple linear residual model \( \hat r_\theta(z) = w^\top z \) that has been learned offline.


import org.ejml.simple.SimpleMatrix;

public class ResidualController1DOF {
    private final double mNom;
    private final double bNom;
    private final double gNom;
    private final double Kp;
    private final double Kd;
    private final SimpleMatrix w; // 4x1 vector of learned weights
    private final double rMax;

    public ResidualController1DOF(double mNom, double bNom, double gNom,
                                  double Kp, double Kd,
                                  SimpleMatrix w, double rMax) {
        this.mNom = mNom;
        this.bNom = bNom;
        this.gNom = gNom;
        this.Kp = Kp;
        this.Kd = Kd;
        this.w = w;
        this.rMax = rMax;
    }

    private double nominalTorque(double q, double qd,
                                 double qdDes, double qddDes) {
        double e = q - qdDes;
        double edot = qd - qdDes;
        double v = qddDes - Kd * edot - Kp * e;
        return mNom * v + bNom * qd + gNom * Math.sin(q);
    }

    private double saturate(double x) {
        if (x > rMax) return rMax;
        if (x < -rMax) return -rMax;
        return x;
    }

    public double computeTorque(double q, double qd,
                                double qdDes, double qddDes) {
        // z = [q, qd, qdDes, qddDes]^T
        double[] zArr = new double[] { q, qd, qdDes, qddDes };
        SimpleMatrix z = new SimpleMatrix(4, 1, true, zArr);

        double tauNom = nominalTorque(q, qd, qdDes, qddDes);
        double rHat = w.dot(z);
        rHat = saturate(rHat);

        return tauNom + rHat;
    }
}
      

The weight vector w is obtained by solving a ridge regression problem on \( (z_k,y_k) \). In a teaching lab, students can implement the least-squares solution \( w = (Z^\top Z + \lambda I)^{-1} Z^\top y \) using EJML.

8. MATLAB/Simulink Implementation

MATLAB is widely used for control and Simulink is convenient for simulation. We outline:

  • A MATLAB script to simulate the nominal controller and collect data.
  • Residual learning via fitrgp (Gaussian Process regression).
  • Integration of the learned residual into a Simulink model via a MATLAB Function block.

% Parameters
m_true = 1.2; b_true = 0.3; g_true = 9.81;
m_nom  = 1.0; b_nom  = 0.1; g_nom  = 9.81;
Kp = 25; Kd = 10;

dt = 0.002;
T  = 12.0;
t  = 0:dt:T;
N  = numel(t);

q  = 0;
qd = 0;

Z = zeros(N,4);
Y = zeros(N,1);

disturbance = @(q, qd) 0.4 * sin(2*q) + 0.05 * qd.^3;

for k = 1:N
    tk = t(k);
    A = 0.7; w = 0.7;
    qd_des    = A * sin(w * tk);
    qd_dotdes = A * w * cos(w * tk);
    qdd_des   = -A * w^2 * sin(w * tk);

    e    = q - qd_des;
    edot = qd - qd_dotdes;
    v    = qdd_des - Kd * edot - Kp * e;
    tau_nom = m_nom * v + b_nom * qd + g_nom * sin(q);

    d_true = disturbance(q, qd);
    qdd    = (tau_nom - b_true * qd - g_true * sin(q) - d_true) / m_true;

    % Euler integration
    qd = qd + dt * qdd;
    q  = q  + dt * qd;

    Z(k,:) = [q, qd, qd_des, qdd_des];
    Y(k)   = tau_nom - (m_nom * v + b_nom * qd + g_nom * sin(q)); % torque error
end

% Fit a GP residual model
gprMdl = fitrgp(Z, Y, "BasisFunction", "none", ...
                       "KernelFunction", "squaredexponential", ...
                       "Sigma", 1e-3);

% Save the model for Simulink
save("residualGP.mat", "gprMdl");
      

In Simulink:

  1. Create a block diagram with a plant block representing the joint dynamics and a nominal computed-torque controller.
  2. Add a MATLAB Function block that:
    • Loads gprMdl from the workspace.
    • Builds the feature vector \( z \).
    • Calls predict(gprMdl, z') to compute the residual.
    • Saturates the residual and adds it to tau_nom.
  3. Run the simulation with and without enabling the residual for comparison.

9. Wolfram Mathematica Implementation

Mathematica has built-in Gaussian Process and regression capabilities. The following code shows how to fit a residual model from simulated data (generation analogous to the Python or MATLAB script) and use it for closed-loop simulation.


(* Assume data has been generated in some way as lists qList, qdList,
   qdDesList, qddDesList, tauList, tauNomList of equal length. *)

zData = Transpose[{qList, qdList, qdDesList, qddDesList}];
yData = tauList - tauNomList;

trainData = Thread[zData -> yData];

(* Gaussian process regression model *)
gpModel = GaussianProcessRegression[trainData];

(* Residual-augmented torque function *)
tauAugmented[q_, qd_, qdDes_, qddDes_] := Module[
  {z, tauNom, rHat, rMax = 2.0},
  z = {q, qd, qdDes, qddDes};
  tauNom = tauNominal[q, qd, qdDes, qddDes]; (* user-defined nominal torque *)
  rHat = Mean[gpModel[z]];
  rHat = Clip[rHat, {-rMax, rMax}];
  tauNom + rHat
];

(* Example of numerical simulation using NDSolve *)
mTrue = 1.2; bTrue = 0.3; gTrue = 9.81;
disturbance[q_, qd_] := 0.4*Sin[2 q] + 0.05*qd^3;

eqns = {
  q''[t] == (tauAugmented[q[t], q'[t], qdDes[t], qddDes[t]] -
              bTrue*q'[t] - gTrue*Sin[q[t]] -
              disturbance[q[t], q'[t]])/mTrue,
  q[0] == 0, q'[0] == 0
};

sol = NDSolve[eqns, {q}, {t, 0, 12}];
      

Here qdDes[t] and qddDes[t] are symbolic definitions of the desired trajectory and its derivatives, and tauNominal is the nominal computed-torque controller.

10. Problems and Solutions

Problem 1 (Residual Target Derivation): Consider the true dynamics \( M(q)\ddot q + C(q,\dot q)\dot q + g(q) + d(q,\dot q,\ddot q) = \tau \) and the nominal model-based controller \( \tau_{\mathrm{nom}}(z) \) built from \( (q,\dot q,\ddot q_d) \). Show that the ideal residual torque \( r^\star(z_k) \) at time \( t_k \) can be expressed as \( r^\star(z_k) = \tau_k - \tau_{\mathrm{nom}}(z_k) \).

Solution: From the true dynamics at time \( t_k \) we have

\[ M(q_k)\ddot q_k + C(q_k,\dot q_k)\dot q_k + g(q_k) + d(q_k,\dot q_k,\ddot q_k) = \tau_k. \]

Suppose we use the control law \( \tau_k = \tau_{\mathrm{nom}}(z_k) + r^\star(z_k) \), where \( \tau_{\mathrm{nom}} \) is computed from the nominal model. Solving for \( r^\star(z_k) \) gives

\[ r^\star(z_k) = \tau_k - \tau_{\mathrm{nom}}(z_k), \]

which is precisely the quantity we can form from logged data. Thus the supervised learning regression target is the torque discrepancy between true and nominal models.

Problem 2 (ISS Bound with Residual): Using the Lyapunov function \( V(e,\dot e) = \tfrac{1}{2}\dot e^\top M(q)\dot e + \tfrac{1}{2}e^\top K_p e \), show that if the residual modeling error satisfies \( \|\tilde r(z)\| \leq \bar r \) and gains \( K_p, K_d \) are large enough, then the tracking error \( e(t) \) is ultimately bounded by a ball whose radius is proportional to \( \bar r \).

Solution: From Section 3, the closed-loop error dynamics with residual are

\[ M(q)\ddot e + C(q,\dot q)\dot e + K_d M_n(q)\dot e + K_p M_n(q)e = \tilde r(z). \]

Differentiating \( V \) and using standard manipulator identities, we obtain

\[ \dot V \leq -\alpha_1 \|\dot e\|^2 - \alpha_2 \|e\|^2 + \beta \|\dot e\| \,\|\tilde r(z)\|, \]

where \( \alpha_1, \alpha_2, \beta > 0 \) depend on bounds on \( M, M_n, K_p, K_d \). Using Young's inequality \( \beta \|\dot e\| \,\|\tilde r\| \leq \tfrac{\alpha_1}{2}\|\dot e\|^2 + \tfrac{\beta^2}{2\alpha_1}\|\tilde r\|^2 \), we get

\[ \dot V \leq -\tfrac{\alpha_1}{2}\|\dot e\|^2 - \alpha_2\|e\|^2 + \tfrac{\beta^2}{2\alpha_1} \|\tilde r\|^2. \]

If \( \|\tilde r(z)\| \leq \bar r \), then \( \dot V \leq -\lambda V + \gamma \bar r^2 \) for some \( \lambda, \gamma > 0 \), which implies that \( V(t) \) (and hence \( \|e(t)\| \)) converges to a ball whose radius is proportional to \( \bar r \). This is the standard ISS ultimate boundedness result.

Problem 3 (Choice of Features): Explain why using only \( q \) as the feature vector \( z \) is typically insufficient for residual dynamics learning, and argue why including \( \dot q \) and \( \ddot q_d \) is beneficial.

Solution: The disturbance term \( d(q,\dot q,\ddot q) \) often depends on velocities and accelerations (e.g. viscous or Coulomb friction, actuator dynamics). If the residual model only sees \( q \), it cannot distinguish cases with identical positions but different velocities or accelerations, so the learned map will blur together distinct regimes and perform poorly. Including \( \dot q \) captures velocity-dependent effects, and including information related to the commanded acceleration (e.g. \( \ddot q_d \)) captures acceleration-dependent unmodeled dynamics. This leads to a better approximation of \( d \) and thus better compensation.

Problem 4 (Regularization in Residual Learning): Consider linear residual regression \( \hat r_\theta(z) = w^\top z \) with ridge penalty \( \lambda \|w\|^2 \). Write the closed-form solution for \( w \) in terms of the data matrix \( Z \in \mathbb{R}^{N \times p} \) and target vector \( y \in \mathbb{R}^N \), and explain how increasing \( \lambda \) interacts with robustness of the learning-augmented controller.

Solution: The ridge regression problem is

\[ \min_w \|Z w - y\|^2 + \lambda \|w\|^2. \]

The solution is

\[ w^\star = (Z^\top Z + \lambda I)^{-1} Z^\top y. \]

Larger \( \lambda \) shrinks \( w \) towards zero, which reduces overfitting and limits the magnitude of the residual corrections. In the closed loop, this corresponds to a more conservative learning augmentation that is less sensitive to noise and model mismatch but may under-compensate the disturbance. From a robustness viewpoint, increasing \( \lambda \) often helps maintain stability by bounding the residual torque.

Problem 5 (Data Coverage for Safety): Suppose the operational region of the robot is specified by \( \|q\| \leq q_{\max}, \|\dot q\| \leq v_{\max} \). Describe a strategy to design the data-collection trajectories \( q_d^{(j)}(t) \) so that the training dataset sufficiently covers this region, and explain how this impacts the confidence in the learned residual when deployed online.

Solution: To cover the operational region, one can design a trajectory ensemble consisting of:

  • Multiple sinusoids with different frequencies and phases so that their superposition visits a dense set of positions and velocities in \( [-q_{\max}, q_{\max}] \times [-v_{\max}, v_{\max}] \).
  • Chirp signals that sweep through relevant frequency bands.
  • Step-like segments (with smooth transitions) to probe acceleration-rich regimes.

By covering the state space region where the controller will be used, the regression model \( \hat r_\theta \) is trained on representative data, reducing extrapolation when deployed. Extrapolation is dangerous because prediction errors can be large and poorly characterized, leading to unforeseen residual torques. Hence coverage of the operational set is critical for safety and performance.

11. Summary

In this lab, we implemented a learning-augmented robot controller based on residual dynamics modeling. Starting from a nominal computed-torque controller, we:

  • Derived the error dynamics and identified the lumped disturbance term that limits tracking performance.
  • Formulated a residual torque model \( \hat r_\theta(z) \) and showed how to construct supervised learning targets from experimental data.
  • Implemented a complete Python simulation with Gaussian Process residuals and sketched C++, Java, MATLAB/Simulink, and Mathematica implementations.
  • Discussed Lyapunov-based ISS arguments that justify bounded residual usage and posed problems emphasizing stability and data coverage.

This lab connects the abstract ideas of residual / delta-model learning from previous lessons with concrete control code, highlighting the trade-offs between tracking accuracy, data requirements, and robustness.

12. References

  1. Ljung, L. (1978). Convergence analysis of parametric identification methods. IEEE Transactions on Automatic Control, 23(5), 770–783.
  2. Narendra, K. S., & Annaswamy, A. M. (1987). A new adaptive law for robust adaptation without persistent excitation. IEEE Transactions on Automatic Control, 32(2), 134–145.
  3. Campi, M. C., & Kumar, P. R. (1998). Adaptive linear systems: identification, stability, and robustness. SIAM Journal on Control and Optimization, 36(6), 1889–1920.
  4. Schölkopf, B., & Smola, A. J. (2002). Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press.
  5. Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
  6. Deisenroth, M. P., Rasmussen, C. E., & Peters, J. (2009). Gaussian process dynamic models for control. Proceedings of the 26th International Conference on Machine Learning, 225–232.
  7. Berkenkamp, F., Schoellig, A. P., & Krause, A. (2015). Safe and robust learning control with Gaussian processes. European Control Conference, 2501–2506.
  8. Khansari-Zadeh, S. M., & Billard, A. (2011). Learning stable nonlinear dynamical systems with Gaussian mixture models. IEEE Transactions on Robotics, 27(5), 943–957.
  9. Hjalmarsson, H. (2005). From experiment design to closed-loop control. Automatica, 41(3), 393–438.
  10. Gevers, M. (2002). A personal view of the development of system identification: A 30-year journey through an exciting field. IEEE Control Systems Magazine, 22(6), 93–105.