We recall the standard rigid robot dynamics in joint space (known from the robotics kinematics/dynamics course):

\[ \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{q} \in \mathbb{R}^n \) are joint positions, \( \mathbf{M}(\mathbf{q}) \) is the inertia matrix, \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) collects Coriolis/centrifugal terms, and \( \mathbf{g}(\mathbf{q}) \) is the gravity vector. Under mild assumptions, the dynamics are linearly parameterizable:

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

where \( \mathbf{Y}(\cdot) \) is the regressor matrix and \( \boldsymbol{\theta} \) is the (unknown but constant) parameter vector containing masses, inertias, COM locations, etc.

Following Lesson 2, define tracking error and filtered (sliding) error for a desired trajectory \( \mathbf{q}_d(t) \):

\[ \mathbf{e} = \mathbf{q} - \mathbf{q}_d, \qquad \dot{\mathbf{e}} = \dot{\mathbf{q}} - \dot{\mathbf{q}}_d \]

\[ \mathbf{s} = \dot{\mathbf{e}} + \boldsymbol{\Lambda}\,\mathbf{e} \]

with \( \boldsymbol{\Lambda} \) a positive definite design matrix. The reference motion is

\[ \dot{\mathbf{q}}_r = \dot{\mathbf{q}}_d - \boldsymbol{\Lambda}\mathbf{e}, \qquad \ddot{\mathbf{q}}_r = \ddot{\mathbf{q}}_d - \boldsymbol{\Lambda}\dot{\mathbf{e}}, \qquad \mathbf{s} = \dot{\mathbf{q}} - \dot{\mathbf{q}}_r. \]

The adaptive computed-torque controller has the structure

\[ \boldsymbol{\tau} = \mathbf{Y}(\mathbf{q},\dot{\mathbf{q}},\dot{\mathbf{q}}_r,\ddot{\mathbf{q}}_r)\,\hat{\boldsymbol{\theta}} - \mathbf{K}_d\,\mathbf{s}, \]

where \( \hat{\boldsymbol{\theta}} \) is the online parameter estimate, \( \mathbf{K}_d \) is a positive definite gain matrix, and \( \mathbf{s} \) is the error signal to be driven to zero. The design problem of this lesson is: choose a parameter update law for \( \hat{\boldsymbol{\theta}} \) that guarantees stability and tracking.

flowchart TD
  Qd["Desired trajectory q_d(t)"] --> E["Tracking error e = q - q_d"]
  E --> S["Filtered error s = dot_e + Lambda e"]
  S --> YREG["Compute regressor \nY(q, dot_q, q_r, dot_q_r, ddot_q_r)"]
  YREG --> UPD["Parameter update: \ndot_hat_theta = -Gamma * Y^T * s"]
  S --> CTRL["Control torque: \ntau = Y * hat_theta - Kd * s"]
  CTRL --> PLANT["Robot joints and dynamics"]
  PLANT -->|measured q, dot_q| E
        

Substituting the control law into the robot dynamics and using the regressor representation yields the error dynamics in terms of \( \mathbf{s} \). Define the parameter error

\[ \tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}. \]

It can be shown (a standard derivation in adaptive control for robots) that

\[ \mathbf{M}(\mathbf{q})\,\dot{\mathbf{s}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\,\mathbf{s} + \mathbf{K}_d\,\mathbf{s} = \mathbf{Y}(\mathbf{q},\dot{\mathbf{q}},\dot{\mathbf{q}}_r,\ddot{\mathbf{q}}_r)\, \tilde{\boldsymbol{\theta}}. \]

This relation is crucial: it shows how the filtered error dynamics depend linearly on the parameter error \( \tilde{\boldsymbol{\theta}} \).

Two standard properties of robot dynamics (derived in robot dynamics textbooks) are used in Lyapunov analysis:

1. Positive definiteness of inertia: \( \mathbf{M}(\mathbf{q}) \) is symmetric and uniformly positive definite: there exist constants \( m_{\min}, m_{\max} > 0 \) such that \( m_{\min}\|\mathbf{x}\|^2 \leq \mathbf{x}^\top \mathbf{M}(\mathbf{q})\mathbf{x} \leq m_{\max}\|\mathbf{x}\|^2 \) for all \( \mathbf{q},\mathbf{x} \).
2. Skew-symmetry property: the matrix \( \dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) is skew-symmetric, so for all \( \mathbf{x} \):

   \[ \mathbf{x}^\top\big(\dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\big)\mathbf{x} = 0. \]

These properties allow us to construct a Lyapunov function that behaves nicely along the trajectories of the adaptive closed-loop system.

We choose a composite Lyapunov candidate incorporating both the tracking error \( \mathbf{s} \) and the parameter error \( \tilde{\boldsymbol{\theta}} \):

\[ V(\mathbf{s},\tilde{\boldsymbol{\theta}}) = \tfrac{1}{2}\,\mathbf{s}^\top \mathbf{M}(\mathbf{q})\,\mathbf{s} + \tfrac{1}{2}\,\tilde{\boldsymbol{\theta}}^\top \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}, \]

where \( \boldsymbol{\Gamma} \) is a symmetric positive definite adaptation gain matrix. Using the skew-symmetry property, its time derivative along trajectories satisfies

\[ \dot{V} = \mathbf{s}^\top\big(\mathbf{M}(\mathbf{q})\dot{\mathbf{s}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\mathbf{s}\big) + \tilde{\boldsymbol{\theta}}^\top \boldsymbol{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta}}}. \]

Substituting the error dynamics \( \mathbf{M}(\mathbf{q})\dot{\mathbf{s}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\mathbf{s} + \mathbf{K}_d\mathbf{s} = \mathbf{Y}\tilde{\boldsymbol{\theta}} \), we obtain

\[ \dot{V} = \mathbf{s}^\top\big(\mathbf{Y}\tilde{\boldsymbol{\theta}} - \mathbf{K}_d\,\mathbf{s}\big) + \tilde{\boldsymbol{\theta}}^\top \boldsymbol{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta}}} = -\mathbf{s}^\top \mathbf{K}_d\,\mathbf{s} + \mathbf{s}^\top \mathbf{Y}\tilde{\boldsymbol{\theta}} + \tilde{\boldsymbol{\theta}}^\top \boldsymbol{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta}}}. \]

To ensure that \( \dot{V} \) is negative semi-definite, we choose the parameter update law so that the cross term \( \mathbf{s}^\top \mathbf{Y}\tilde{\boldsymbol{\theta}} + \tilde{\boldsymbol{\theta}}^\top\boldsymbol{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta}}} \) cancels. Since the true parameters are constant (\( \dot{\boldsymbol{\theta}} = 0 \)), \( \dot{\tilde{\boldsymbol{\theta}}} = \dot{\hat{\boldsymbol{\theta}}} \). Let

\[ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\,\mathbf{Y}^\top(\mathbf{q},\dot{\mathbf{q}},\dot{\mathbf{q}}_r,\ddot{\mathbf{q}}_r)\,\mathbf{s}. \]

Then

\[ \tilde{\boldsymbol{\theta}}^\top\boldsymbol{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta}}} = \tilde{\boldsymbol{\theta}}^\top\boldsymbol{\Gamma}^{-1} \big(-\boldsymbol{\Gamma}\mathbf{Y}^\top\mathbf{s}\big) = -\tilde{\boldsymbol{\theta}}^\top\mathbf{Y}^\top\mathbf{s} = -\mathbf{s}^\top\mathbf{Y}\tilde{\boldsymbol{\theta}}. \]

Hence the cross terms cancel and we obtain

\[ \dot{V} = -\mathbf{s}^\top\mathbf{K}_d\,\mathbf{s} \leq 0. \]

Stability sketch.

• Since \( \mathbf{M}(\mathbf{q}) \) and \( \boldsymbol{\Gamma}^{-1} \) are positive definite, \( V \) is positive definite and radially unbounded in \( (\mathbf{s},\tilde{\boldsymbol{\theta}}) \).
• The derivative \( \dot{V} = -\mathbf{s}^\top\mathbf{K}_d\mathbf{s} \) is negative semi-definite, implying that \( \mathbf{s} \in L_2 \cap L_\infty \), \( \tilde{\boldsymbol{\theta}} \in L_\infty \), and all closed-loop signals are bounded.
• Using the relation between \( \mathbf{s} \) and \( \mathbf{e} \) from Lesson 2 and standard tools (e.g., Barbalat's lemma), we can show \( \mathbf{s}(t) \to \mathbf{0} \) and \( \mathbf{e}(t) \to \mathbf{0} \).
• Convergence of the parameter error \( \tilde{\boldsymbol{\theta}} \) to zero requires an additional excitation condition on the regressor (discussed in Lesson 4).

Thus, the Lyapunov-based update law \( \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\mathbf{Y}^\top\mathbf{s} \) ensures global stability and tracking for the adaptive computed-torque controller under the usual rigid-robot assumptions.

The same update law can be interpreted as gradient descent on an instantaneous performance index. Consider the quadratic cost

\[ J(\hat{\boldsymbol{\theta}}) = \tfrac{1}{2}\,\mathbf{s}^\top\mathbf{s}. \]

Using the error dynamics, one can show that the sensitivity of \( \mathbf{s} \) with respect to \( \hat{\boldsymbol{\theta}} \) is proportional to \( \mathbf{Y} \), and the gradient satisfies

\[ \frac{\partial J}{\partial \hat{\boldsymbol{\theta}}} = \mathbf{Y}^\top\mathbf{s}. \]

The gradient-descent update

\[ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\,\frac{\partial J}{\partial \hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\,\mathbf{Y}^\top\mathbf{s} \]

coincides with the Lyapunov-based law derived earlier. The Lyapunov analysis adds a global stability guarantee, which a naive gradient interpretation does not automatically provide.

Two simple practical modifications are often used:

1. Normalization: to reduce sensitivity to large regressors, a normalized update may be used:

   \[ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\, \frac{\mathbf{Y}^\top\mathbf{s}} {1 + \|\mathbf{Y}\|^2}. \]

   This preserves the direction of adaptation but scales its magnitude.
2. Leakage (or "σ-modification"): to avoid slow drift of parameters in the presence of small unmodeled dynamics, a leakage term can be added:

   \[ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\,\mathbf{Y}^\top\mathbf{s} - \sigma \hat{\boldsymbol{\theta}}, \qquad \sigma > 0. \]

   This drives the estimates toward zero when the error signal \( \mathbf{s} \) is small for a long time, leading to bounded parameter estimates and an ultimately bounded tracking error.

Full Lyapunov proofs for these variants use slightly modified Lyapunov functions that add terms depending on \( \hat{\boldsymbol{\theta}} \) itself, but the basic structure remains: choose update laws that cancel cross terms and introduce negative terms in \( V \).

Real robot controllers are digital, with sampling period \( T_s \). A straightforward discretization of the continuous-time update law at sampling instants \( t_k = kT_s \) is:

\[ \hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k - T_s\,\boldsymbol{\Gamma}\, \mathbf{Y}^\top_k\,\mathbf{s}_k, \]

where \( \mathbf{Y}_k = \mathbf{Y}(\mathbf{q}_k,\dot{\mathbf{q}}_k,\dot{\mathbf{q}}_{r,k},\ddot{\mathbf{q}}_{r,k}) \) and \( \mathbf{s}_k \) are evaluated at time \( t_k \). Writing \( \boldsymbol{\Gamma}_d = T_s\boldsymbol{\Gamma} \), we get

\[ \hat{\boldsymbol{\theta}}_{k+1} = \hat{\boldsymbol{\theta}}_k - \boldsymbol{\Gamma}_d \mathbf{Y}^\top_k\mathbf{s}_k. \]

The matrix \( \boldsymbol{\Gamma}_d \) plays the role of a discrete-time learning rate. Too large values can induce oscillations or numerical instability; too small values yield very slow adaptation.

A practical implementation often includes:

• Bounds on \( \hat{\boldsymbol{\theta}}_k \) to avoid unrealistically large parameters.
• Dead zones where adaptation is frozen when \( \|\mathbf{s}_k\| \) is below a threshold.
• Optional normalization, replacing \( \mathbf{Y}^\top_k\mathbf{s}_k \) by a normalized version.

flowchart TD
  SAMP["At each sample k"] --> MEAS["Read q_k, dot_q_k"]
  MEAS --> ERRK["Compute e_k, s_k"]
  ERRK --> YK["Compute regressor Y_k"]
  YK --> UPDK["Update hat_theta_{k+1} = hat_theta_k - Gamma_d * Y_k^T * s_k"]
  UPDK --> CTRLK["Compute tau_k = Y_k * hat_theta_k - Kd * s_k"]
  CTRLK --> ACTK["Send tau_k to robot"]
  ACTK --> SAMP
        

We now provide skeleton implementations of the discrete-time update law for a generic joint-space adaptive computed-torque controller. In all examples, we assume that a function computing the regressor \( \mathbf{Y} \) is available from the robot model (e.g., from a robotics toolbox or implemented from known dynamics).

6.1 Python (NumPy, typical robotics toolbox)

import numpy as np

class AdaptiveCTControllerPython:
    def __init__(self, n_dof, Kd, Lambda, Gamma_d, robot_model):
        """
        n_dof      : number of joints
        Kd         : (n_dof x n_dof) positive definite damping matrix
        Lambda     : (n_dof x n_dof) positive definite filtering matrix
        Gamma_d    : (p x p) discrete adaptation gain matrix
        robot_model: object providing regressor Y(q, dq, qr, dqr, ddqr)
        """
        self.n = n_dof
        self.Kd = Kd
        self.Lambda = Lambda
        self.Gamma_d = Gamma_d
        # p = number of base parameters of the robot
        self.robot_model = robot_model
        self.theta_hat = np.zeros(robot_model.num_params)

    def step(self, q, dq, qd, dqd, ddqd, dt):
        # tracking error
        e = q - qd
        de = dq - dqd
        # filtered error
        s = de + self.Lambda @ e
        # reference motion
        dqr = dqd - self.Lambda @ e
        ddqr = ddqd - self.Lambda @ de
        # regressor matrix
        Y = self.robot_model.regressor(q, dq, dqr, ddqr)  # shape (n, p)
        # control torque (use current theta_hat)
        tau = Y @ self.theta_hat - self.Kd @ s
        # parameter update (discrete form)
        grad = Y.T @ s
        self.theta_hat = self.theta_hat - self.Gamma_d @ grad

        return tau, self.theta_hat
      

6.2 C++ (Eigen-based implementation, joint-space controller)

#include "Eigen/Dense"

class RobotModel {
public:
    int numParams() const;
    Eigen::MatrixXd regressor(const Eigen::VectorXd& q,
                              const Eigen::VectorXd& dq,
                              const Eigen::VectorXd& qr,
                              const Eigen::VectorXd& dqr,
                              const Eigen::VectorXd& ddqr) const;
};

class AdaptiveCTControllerCpp {
public:
    AdaptiveCTControllerCpp(int n_dof,
                            const Eigen::MatrixXd& Kd,
                            const Eigen::MatrixXd& Lambda,
                            const Eigen::MatrixXd& Gamma_d,
                            const RobotModel& model)
        : n(n_dof),
          Kd_(Kd),
          Lambda_(Lambda),
          Gamma_d_(Gamma_d),
          robot_(model)
    {
        theta_hat_.setZero(robot_.numParams());
    }

    void step(const Eigen::VectorXd& q,
              const Eigen::VectorXd& dq,
              const Eigen::VectorXd& qd,
              const Eigen::VectorXd& dqd,
              const Eigen::VectorXd& ddqd,
              Eigen::VectorXd& tau_out)
    {
        Eigen::VectorXd e  = q - qd;
        Eigen::VectorXd de = dq - dqd;
        Eigen::VectorXd s  = de + Lambda_ * e;

        Eigen::VectorXd dqr  = dqd - Lambda_ * e;
        Eigen::VectorXd ddqr = ddqd - Lambda_ * de;

        Eigen::MatrixXd Y = robot_.regressor(q, dq, dqr, ddqr);

        tau_out = Y * theta_hat_ - Kd_ * s;

        Eigen::VectorXd grad = Y.transpose() * s;
        theta_hat_ = theta_hat_ - Gamma_d_ * grad;
    }

    const Eigen::VectorXd& thetaHat() const { return theta_hat_; }

private:
    int n;
    Eigen::MatrixXd Kd_, Lambda_, Gamma_d_;
    Eigen::VectorXd theta_hat_;
    const RobotModel& robot_;
};
      

6.3 Java (simple array-based controller skeleton)

public interface RobotRegressor {
    // returns Y(q, dq, qr, dqr, ddqr) as n x p matrix
    double[][] regressor(double[] q, double[] dq,
                         double[] qr, double[] dqr, double[] ddqr);
    int numParams();
}

public class AdaptiveCTControllerJava {
    private final int n;
    private final int p;
    private final double[][] Kd;
    private final double[][] Lambda;
    private final double[][] Gamma_d;
    private final RobotRegressor model;
    private final double[] thetaHat;

    public AdaptiveCTControllerJava(int nDof,
                                    double[][] Kd,
                                    double[][] Lambda,
                                    double[][] Gamma_d,
                                    RobotRegressor model) {
        this.n = nDof;
        this.Kd = Kd;
        this.Lambda = Lambda;
        this.Gamma_d = Gamma_d;
        this.model = model;
        this.p = model.numParams();
        this.thetaHat = new double[p];
    }

    public double[] step(double[] q, double[] dq,
                         double[] qd, double[] dqd, double[] ddqd) {
        double[] e  = sub(q, qd);
        double[] de = sub(dq, dqd);

        double[] s = add(de, matVec(Lambda, e));

        double[] dqr  = sub(dqd, matVec(Lambda, e));
        double[] ddqr = sub(ddqd, matVec(Lambda, de));

        double[][] Y = model.regressor(q, dq, dqr, dqr, ddqr);

        double[] Ytheta = matVec(Y, thetaHat);
        double[] Kds    = matVec(Kd, s);

        double[] tau = sub(Ytheta, Kds);

        double[] grad = matTVec(Y, s);
        double[] GammaGrad = matVec(Gamma_d, grad);
        for (int i = 0; i != p; ++i) {
            thetaHat[i] -= GammaGrad[i];
        }
        return tau;
    }

    // helper linear algebra (sub, add, matVec, matTVec) to be implemented
}
      

6.4 MATLAB / Simulink

In MATLAB, a similar discrete-time controller can be written as:

function [tau, theta_hat_next] = adaptive_ct_step(q, dq, qd, dqd, ddqd, ...
                                                  theta_hat, Kd, Lambda, Gamma_d, robot_model)
    % tracking error
    e  = q - qd;
    de = dq - dqd;

    % filtered error
    s = de + Lambda * e;

    % reference motion
    dqr  = dqd - Lambda * e;
    ddqr = ddqd - Lambda * de;

    % regressor matrix (n x p)
    Y = robot_model.regressor(q, dq, dqr, ddqr);

    % control torque
    tau = Y * theta_hat - Kd * s;

    % parameter update
    grad = Y.' * s;
    theta_hat_next = theta_hat - Gamma_d * grad;
end
      

In Simulink, this logic can be placed inside a "MATLAB Function" block implementing the parameter update and torque computation. The state theta_hat is stored using a Unit Delay or Memory block, feeding the previous estimate into the function block at each sampling instant.

6.5 Wolfram Mathematica (continuous-time simulation)

Mathematica can be used to simulate the coupled differential equations for \( \mathbf{q}(t) \) and \( \hat{\boldsymbol{\theta}}(t) \). For illustration, consider a simple one-DOF joint with dynamics \( m\ddot{q} = \tau \) and regressor \( Y(q,\dot{q},\ddot{q}_r) = \ddot{q}_r \).

(* parameters *)
Lambda = 5.0;
Kd     = 10.0;
Gamma  = 2.0;

qd[t_]   := Sin[t];
dqd[t_]  := Cos[t];
ddqd[t_] := -Sin[t];

(* state variables: q[t], dq[t], thetaHat[t] *)
eqs = {
  (* reference terms *)
  qr[t]   == qd[t] - Lambda (q[t] - qd[t]),
  dqr[t]  == dqd[t] - Lambda (q'[t] - dqd[t]),
  ddqr[t] == ddqd[t] - Lambda (q''[t] - ddqd[t]),

  (* sliding error *)
  s[t] == q'[t] - dqr[t],

  (* regressor and torque *)
  Y[t]   == ddqr[t],
  tau[t] == Y[t] thetaHat[t] - Kd s[t],

  (* plant dynamics: m ddq = tau with unknown m *)
  q''[t] == tau[t],  (* true m = 1 *)

  (* parameter update law *)
  thetaHat'[t] == -Gamma Y[t] s[t]
};

ics = {q[0] == 0.0, q'[0] == 0.0, thetaHat[0] == 0.5};

sol = NDSolve[Flatten[{eqs, ics}], {q, thetaHat, s}, {t, 0, 10}];

Plot[{Evaluate[q[t] /. sol], qd[t]}, {t, 0, 10},
     PlotLegends -> {"q(t)", "qd(t)"}]

Plot[Evaluate[thetaHat[t] /. sol], {t, 0, 10},
     PlotLegends -> {"thetaHat(t)"}]
      

This code integrates the true system and the adaptive law simultaneously, illustrating how \( q(t) \) tracks \( q_d(t) \) and how the estimate \( \hat{\theta}(t) \) evolves.

In this lesson we derived parameter update laws for adaptive computed-torque control of robot manipulators. Exploiting the linear parameterization of robot dynamics and the skew-symmetry property of inertia and Coriolis terms, we constructed a composite Lyapunov function in the filtered tracking error and the parameter error. By selecting the update law \( \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}\mathbf{Y}^\top\mathbf{s} \), we cancelled cross terms and obtained a negative semi-definite Lyapunov derivative, yielding boundedness and asymptotic tracking. We showed how this update can be interpreted as gradient descent, extended it to discrete time, and provided implementation templates in Python, C++, Java, MATLAB/Simulink, and Mathematica. The next lesson will analyze conditions under which the parameters themselves converge, emphasizing the role of excitation.