After completing this lesson, the student should be able to:

• distinguish unknown physical parameters from unmodeled mechanical dynamics;
• express translational and rotational equations in the regression form \(Y\theta\);
• construct a filtered tracking error and reference velocity;
• derive an adaptive computed-torque controller and its update law;
• prove boundedness and asymptotic tracking using the skew-symmetry property of rigid-body dynamics; and
• explain why tracking convergence does not automatically imply exact recovery of mass or inertia.

Typical sources of parametric uncertainty include an unknown payload, uncertain center of mass, imprecise link inertia, temperature-dependent viscous friction, and replacement of an end effector. These quantities may change the coefficients of a known mechanical model while leaving the model structure intact. Adaptive control is especially suitable for this case.

In contrast, gearbox backlash, flexible modes, Coulomb friction near zero velocity, actuator dead zones, and time-varying impacts may not be representable by the selected constant parameter vector. Such effects are unmodeled dynamics or disturbances rather than ordinary parameter uncertainty.

Consider a nonlinear mechanical system with generalized coordinate \(q\). Its equation may contain nonlinear functions such as \(\sin q\), products of velocities, and configuration-dependent inertia. Nevertheless, the unknown physical constants can frequently be collected in a vector \(\theta\) so that the dynamics are affine in those constants.

\[ M(q,\theta)\ddot q+C(q,\dot q,\theta)\dot q+g(q,\theta) =Y(q,\dot q,v,a)\theta \]

The equality above is interpreted for selected auxiliary velocity and acceleration signals \(v\) and \(a\). The regressor \(Y\) is known because it depends on measured or commanded signals; the parameter vector \(\theta\) is unknown.

This distinction is fundamental. The map from state to torque is nonlinear, but the map from the unknown physical coefficients to torque is linear. Adaptive laws operate in this parameter space.

2.1 Translational example

A mass-spring-damper system with unknown mass \(m\), damping \(c\), and stiffness \(k\) satisfies

\[ m\ddot q+c\dot q+kq=u. \]

For arbitrary known signals \(a_r\), measured \(\dot q\), and \(q\), define

\[ Y=\begin{bmatrix}a_r & \dot q & q\end{bmatrix},\qquad \theta=\begin{bmatrix}m & c & k\end{bmatrix}^{T}. \]

Then \(Y\theta=ma_r+c\dot q+kq\).

2.2 Rotational example

A single rotary link with unknown inertia \(J\), viscous coefficient \(b\), payload mass \(m_p\), and payload moment arm \(\ell\) can be modeled as

\[ J\ddot q+b\dot q+m_p g\ell\sin q=\tau. \]

Let \(h=m_p g\ell\). For a selected reference acceleration \(\ddot q_r\),

\[ Y(q,\dot q,\ddot q_r)= \begin{bmatrix}\ddot q_r & \dot q & \sin q\end{bmatrix}, \qquad \theta=\begin{bmatrix}J & b & h\end{bmatrix}^{T}. \]

The unknown payload mass appears through the gravity coefficient \(h\), and the unknown rotational inertia appears directly as \(J\).

Let \(q_d(t)\) be a twice differentiable desired trajectory. Define the position and velocity errors as

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

Choose \(\lambda>0\) and introduce the filtered error

\[ s=\dot e+\lambda e. \]

Equivalently, define the reference velocity and acceleration

\[ \dot q_r=\dot q_d-\lambda e,\qquad \ddot q_r=\ddot q_d-\lambda\dot e, \]

so that

\[ s=\dot q-\dot q_r,\qquad \dot s=\ddot q-\ddot q_r. \]

If \(s(t)→0\), then the stable first-order error equation \(\dot e+\lambda e=s\) implies \(e(t)→0\). Therefore the adaptive design may focus on driving \(s\) to zero.

flowchart TD
  QD["Desired motion qd, qd_dot, qd_ddot"] --> REF["Reference filter"]
  Q["Measured q and q_dot"] --> ERR["Compute e and s"]
  QD --> ERR
  ERR --> REF
  REF --> REG["Mechanical regressor Y"]
  Q --> REG
  EST["Parameter estimates"] --> CTRL["Adaptive feedforward \nplus damping"]
  REG --> CTRL
  ERR --> CTRL
  CTRL --> PLANT["Unknown mechanical plant"]
  PLANT --> Q
  REG --> UPDATE["Parameter update law"]
  ERR --> UPDATE
  UPDATE --> EST
        

For the rotary system, define

\[ Y=\begin{bmatrix}\ddot q_r & \dot q & \sin q\end{bmatrix}, \qquad \hat\theta= \begin{bmatrix}\hat J & \hat b & \hat h\end{bmatrix}^{T}. \]

Choose the certainty-equivalent adaptive control law

\[ \tau=Y\hat\theta-k_s s,\qquad k_s>0, \]

and the gradient update law

\[ \dot{\hat\theta}=-\Gamma Y^{T}s, \qquad \Gamma=\Gamma^{T}>0. \]

Define the estimation error \(\tilde\theta=\hat\theta-\theta\). Since the true parameters are constant, \(\dot{\tilde\theta}=\dot{\hat\theta}\).

4.1 Closed-loop filtered-error equation

Using \(\dot s=\ddot q-\ddot q_r\),

\[ \begin{aligned} J\dot s &=J\ddot q-J\ddot q_r\\ &=\tau-b\dot q-h\sin q-J\ddot q_r\\ &=Y\hat\theta-k_s s-Y\theta\\ &=Y\tilde\theta-k_s s. \end{aligned} \]

Hence

\[ J\dot s+k_s s=Y\tilde\theta. \]

4.2 Lyapunov proof

Consider

\[ V=\frac{1}{2}Js^2+\frac{1}{2}\tilde\theta^{T} \Gamma^{-1}\tilde\theta. \]

Because \(J>0\) and \(\Gamma>0\), this function is positive definite in \((s,\tilde\theta)\). Its derivative is

\[ \begin{aligned} \dot V &=Js\dot s+\tilde\theta^{T}\Gamma^{-1} \dot{\tilde\theta}\\ &=s\left(Y\tilde\theta-k_s s\right) +\tilde\theta^{T}\Gamma^{-1} \left(-\Gamma Y^{T}s\right)\\ &=-k_s s^2. \end{aligned} \]

Therefore \(V(t)\leq V(0)\), so \(s\) and \(\tilde\theta\) are bounded, and \(s\in L_2\). Under bounded desired signals and the stated mechanical model, \(\dot s\) is bounded. Barbalat's lemma then gives \(s(t)→0\), and the stable filter gives \(e(t)→0\).

The same construction applies to \(m\ddot q+c\dot q+kq=u\). Define

\[ Y=\begin{bmatrix}\ddot q_r & \dot q & q\end{bmatrix},\qquad \hat\theta=\begin{bmatrix}\hat m & \hat c & \hat k\end{bmatrix}^{T}. \]

The controller and adaptation law are

\[ u=Y\hat\theta-k_s s, \qquad \dot{\hat\theta}=-\Gamma Y^{T}s. \]

The resulting error dynamics are

\[ m\dot s+k_s s=Y\tilde\theta. \]

With \(V=\frac{1}{2}ms^2+\frac{1}{2}\tilde\theta^{T} \Gamma^{-1}\tilde\theta\), the same cancellation yields \(\dot V=-k_s s^2\).

An unknown attached payload changes the effective translational mass. If the spring and damping coefficients are known, the parameter vector may be reduced to \(\theta=m\) and the scalar update law becomes

\[ \dot{\hat m}=-\gamma_m\ddot q_r s. \]

This estimate should not be interpreted as an accurate physical scale measurement unless the regressor is persistently exciting. Its primary role is to generate the feedforward force required for tracking.

For a fully actuated rigid mechanical system,

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

the classical structural properties are:

1. \(M(q)=M^{T}(q)>0\) for every admissible configuration;
2. \(\dot M(q)-2C(q,\dot q)\) is skew-symmetric for a consistent choice of \(C\); and
3. the dynamics can be written as a known regressor times a constant physical parameter vector:

\[ M(q)\ddot q_r+C(q,\dot q)\dot q_r+g(q) =Y(q,\dot q,\dot q_r,\ddot q_r)\theta. \]

For vector tracking error, choose

\[ e=q-q_d,\qquad \dot q_r=\dot q_d-\Lambda e,\qquad s=\dot q-\dot q_r=\dot e+\Lambda e, \]

where \(\Lambda=\Lambda^{T}>0\). The adaptive controller is

\[ \tau=Y\hat\theta-K_Ds, \qquad \dot{\hat\theta}=-\Gamma Y^{T}s, \]

with \(K_D=K_D^{T}>0\) and \(\Gamma=\Gamma^{T}>0\).

6.1 Vector error dynamics

Substitute \(\dot q=s+\dot q_r\) and \(\ddot q=\dot s+\ddot q_r\) into the plant:

\[ M\dot s+Cs+M\ddot q_r+C\dot q_r+g=\tau. \]

By the regression identity and the controller,

\[ M\dot s+Cs+K_Ds=Y\tilde\theta. \]

6.2 Lyapunov proof using mechanical skew symmetry

Take

\[ V=\frac{1}{2}s^{T}M(q)s+\frac{1}{2} \tilde\theta^{T}\Gamma^{-1}\tilde\theta. \]

Its derivative is

\[ \begin{aligned} \dot V &=s^{T}M\dot s+\frac{1}{2}s^{T}\dot M s +\tilde\theta^{T}\Gamma^{-1}\dot{\tilde\theta}\\ &=s^{T}\left(Y\tilde\theta-Cs-K_Ds\right) +\frac{1}{2}s^{T}\dot M s -\tilde\theta^{T}Y^{T}s\\ &=-s^{T}K_Ds+\frac{1}{2}s^{T} \left(\dot M-2C\right)s\\ &=-s^{T}K_Ds. \end{aligned} \]

The final equality follows because \(x^{T}Ax=0\) for every skew-symmetric matrix \(A\). This is the key mechanical cancellation that extends the scalar proof to nonlinear coupled manipulators.

7.1 Guaranteed conclusions

Under the ideal model assumptions, bounded desired trajectory and derivatives, constant true parameters, and positive-definite design matrices, the Lyapunov argument establishes:

\[ s,\tilde\theta\in L_{\infty},\qquad s\in L_2,\qquad s(t)→0,\qquad e(t)→0. \]

7.2 Parameter convergence requires excitation

The derivative \(\dot V=-s^{T}K_Ds\) is negative semidefinite, not negative definite in \(\tilde\theta\). Consequently, the controller can achieve excellent tracking while the estimates settle at values different from the true mass or inertia.

A sufficient excitation condition is that there exist constants \(T>0\) and \(\alpha>0\) such that

\[ \int_{t}^{t+T}Y^{T}(\tau)Y(\tau)\,d\tau \geq \alpha I \quad\text{for all }t\geq0. \]

If the commanded motion is nearly constant, several columns of the mechanical regressor may be weak or dependent. In that case, multiple parameter vectors can produce nearly identical torques along the observed trajectory.

7.3 Physical estimates and controller parameters

The vector \(\theta\) normally contains base inertial combinations rather than every geometric parameter separately. For example, a gravity term may identify the product \(m_p\ell\) rather than the payload mass and center-of-mass distance individually. Therefore observability and identifiability must be interpreted at the level of the chosen regressor.

The downloadable programs use the rotary plant

\[ 2.5\ddot q+0.8\dot q+5\sin q=\tau, \]

but initialize the controller with

\[ \hat\theta(0)= \begin{bmatrix}1.2 & 0.2 & 2.0\end{bmatrix}^{T}. \]

The desired motion is

\[ q_d(t)=0.6\sin(0.5t)+0.2\sin(1.2t), \]

with design constants

\[ \lambda=3,\qquad k_s=8,\qquad \Gamma=\operatorname{diag}(2,1,1.5). \]

A fixed-step fourth-order Runge-Kutta method integrates the combined plant and adaptation dynamics. Each implementation writes a CSV file containing tracking, control, filtered error, and parameter-estimate histories.

flowchart TD
  INIT["Set true plant and incorrect initial estimates"] --> DES["Evaluate desired position velocity acceleration"]
  DES --> FILT["Compute tracking error reference motion and s"]
  FILT --> Y["Build inertia friction gravity regressor"]
  Y --> TORQUE["Compute adaptive torque"]
  TORQUE --> DYN["Evaluate plant and estimator derivatives"]
  DYN --> RK["Advance one RK4 step"]
  RK --> LOG["Store tracking torque and estimates"]
  LOG --> CHECK["Final time reached"]
  CHECK -->|no| DES
  CHECK -->|yes| OUT["Write CSV and report RMSE"]
        

The Python implementation uses NumPy for vector operations and Matplotlib for tracking and parameter plots. The controller itself is implemented directly rather than through a control toolbox so that every adaptive equation remains visible.

Chapter18_Lesson4.py

"""Chapter 18, Lesson 4: adaptive control of a rotary mechanical system.

The plant is
    J*q_ddot + b*q_dot + h*sin(q) = tau,
where J, b, and h = m*g*l are unknown to the controller.

The controller uses the filtered error s and a linearly parameterized regressor:
    tau = Y @ theta_hat - k_s*s
    theta_hat_dot = -Gamma*Y*s
"""

from __future__ import annotations

import math
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


TRUE_THETA = np.array([2.5, 0.8, 5.0], dtype=float)  # [J, b, h]
LAMBDA = 3.0
K_S = 8.0
GAMMA = np.diag([2.0, 1.0, 1.5])
DT = 0.001
T_FINAL = 20.0


def desired_motion(t: float) -> tuple[float, float, float]:
    """Return q_d, q_d_dot, and q_d_ddot."""
    q_d = 0.6 * math.sin(0.5 * t) + 0.2 * math.sin(1.2 * t)
    q_d_dot = 0.3 * math.cos(0.5 * t) + 0.24 * math.cos(1.2 * t)
    q_d_ddot = -0.15 * math.sin(0.5 * t) - 0.288 * math.sin(1.2 * t)
    return q_d, q_d_dot, q_d_ddot


def closed_loop_rhs(t: float, x: np.ndarray) -> np.ndarray:
    """Combined plant and parameter-estimator dynamics."""
    q, q_dot = x[0], x[1]
    theta_hat = x[2:5]
    j_true, b_true, h_true = TRUE_THETA

    q_d, q_d_dot, q_d_ddot = desired_motion(t)
    e = q - q_d
    e_dot = q_dot - q_d_dot
    s = e_dot + LAMBDA * e

    q_r_dot = q_d_dot - LAMBDA * e
    q_r_ddot = q_d_ddot - LAMBDA * e_dot

    # J*q_r_ddot + b*q_dot + h*sin(q) = Y*theta.
    regressor = np.array([q_r_ddot, q_dot, math.sin(q)], dtype=float)
    torque = float(regressor @ theta_hat - K_S * s)

    q_ddot = (torque - b_true * q_dot - h_true * math.sin(q)) / j_true
    theta_hat_dot = -GAMMA @ regressor * s

    return np.concatenate(([q_dot, q_ddot], theta_hat_dot))


def rk4_step(t: float, x: np.ndarray, dt: float) -> np.ndarray:
    """Advance one fixed-step fourth-order Runge-Kutta step."""
    k1 = closed_loop_rhs(t, x)
    k2 = closed_loop_rhs(t + 0.5 * dt, x + 0.5 * dt * k1)
    k3 = closed_loop_rhs(t + 0.5 * dt, x + 0.5 * dt * k2)
    k4 = closed_loop_rhs(t + dt, x + dt * k3)
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate() -> dict[str, np.ndarray]:
    steps = int(round(T_FINAL / DT)) + 1
    time = np.linspace(0.0, T_FINAL, steps)
    state = np.zeros((steps, 5), dtype=float)
    state[0] = np.array([0.0, 0.0, 1.2, 0.2, 2.0], dtype=float)

    q_d = np.zeros(steps)
    s_history = np.zeros(steps)
    torque = np.zeros(steps)

    for index, t in enumerate(time):
        q, q_dot = state[index, 0], state[index, 1]
        theta_hat = state[index, 2:5]
        q_d[index], q_d_dot, q_d_ddot = desired_motion(float(t))
        e = q - q_d[index]
        e_dot = q_dot - q_d_dot
        s_history[index] = e_dot + LAMBDA * e
        q_r_ddot = q_d_ddot - LAMBDA * e_dot
        regressor = np.array([q_r_ddot, q_dot, math.sin(q)])
        torque[index] = regressor @ theta_hat - K_S * s_history[index]

        if index + 1 < steps:
            state[index + 1] = rk4_step(float(t), state[index], DT)

    return {
        "time": time,
        "state": state,
        "q_d": q_d,
        "s": s_history,
        "torque": torque,
    }


def main() -> None:
    result = simulate()
    time = result["time"]
    state = result["state"]
    q_d = result["q_d"]
    tracking_error = state[:, 0] - q_d

    rmse = float(np.sqrt(np.mean(tracking_error**2)))
    print(f"Tracking RMSE: {rmse:.6f} rad")
    print("True [J, b, h]:", TRUE_THETA)
    print("Final estimates:", state[-1, 2:5])
    print("Parameter convergence is not guaranteed without persistent excitation.")

    output = np.column_stack(
        (time, state[:, 0], q_d, tracking_error, result["s"], result["torque"], state[:, 2:5])
    )
    header = "time,q,q_d,error,s,torque,J_hat,b_hat,h_hat"
    np.savetxt("Chapter18_Lesson4_results.csv", output, delimiter=",", header=header, comments="")

    plt.figure(figsize=(9, 5))
    plt.plot(time, q_d, label="desired")
    plt.plot(time, state[:, 0], "--", label="actual")
    plt.xlabel("time [s]")
    plt.ylabel("angle [rad]")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter18_Lesson4_tracking.png", dpi=180)

    plt.figure(figsize=(9, 5))
    labels = ["J_hat", "b_hat", "h_hat"]
    for column, label in zip(range(2, 5), labels):
        plt.plot(time, state[:, column], label=label)
    plt.xlabel("time [s]")
    plt.ylabel("parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter18_Lesson4_parameters.png", dpi=180)

    if Path("Chapter18_Lesson4_tracking.png").exists():
        print("Saved CSV and PNG result files in the current directory.")


if __name__ == "__main__":
    main()

The C++17 program uses only the standard library. Fixed-size std::array objects represent the state and parameter vectors, which makes the implementation suitable as a starting point for an embedded or real-time version.

Chapter18_Lesson4.cpp

// Chapter 18, Lesson 4: adaptive control of a rotary mechanical system.
// Build: g++ -std=c++17 -O2 Chapter18_Lesson4.cpp -o Chapter18_Lesson4

#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>

namespace {
constexpr double kJTrue = 2.5;
constexpr double kBTrue = 0.8;
constexpr double kHTrue = 5.0;
constexpr double kLambda = 3.0;
constexpr double kKs = 8.0;
constexpr double kDt = 0.001;
constexpr double kFinalTime = 20.0;
constexpr std::array<double, 3> kGamma{2.0, 1.0, 1.5};

using State = std::array<double, 5>;  // q, q_dot, J_hat, b_hat, h_hat

struct Desired {
    double q;
    double q_dot;
    double q_ddot;
};

Desired desiredMotion(double t) {
    return {
        0.6 * std::sin(0.5 * t) + 0.2 * std::sin(1.2 * t),
        0.3 * std::cos(0.5 * t) + 0.24 * std::cos(1.2 * t),
        -0.15 * std::sin(0.5 * t) - 0.288 * std::sin(1.2 * t)
    };
}

State addScaled(const State& x, const State& dx, double scale) {
    State result{};
    for (std::size_t i = 0; i < x.size(); ++i) {
        result[i] = x[i] + scale * dx[i];
    }
    return result;
}

State rhs(double t, const State& x) {
    const double q = x[0];
    const double qDot = x[1];
    const Desired d = desiredMotion(t);

    const double e = q - d.q;
    const double eDot = qDot - d.q_dot;
    const double s = eDot + kLambda * e;
    const double qRddot = d.q_ddot - kLambda * eDot;

    const std::array<double, 3> y{qRddot, qDot, std::sin(q)};
    const double torque = y[0] * x[2] + y[1] * x[3] + y[2] * x[4] - kKs * s;
    const double qDdot = (torque - kBTrue * qDot - kHTrue * std::sin(q)) / kJTrue;

    State dx{};
    dx[0] = qDot;
    dx[1] = qDdot;
    for (std::size_t i = 0; i < 3; ++i) {
        dx[i + 2] = -kGamma[i] * y[i] * s;
    }
    return dx;
}

State rk4Step(double t, const State& x) {
    const State k1 = rhs(t, x);
    const State k2 = rhs(t + 0.5 * kDt, addScaled(x, k1, 0.5 * kDt));
    const State k3 = rhs(t + 0.5 * kDt, addScaled(x, k2, 0.5 * kDt));
    const State k4 = rhs(t + kDt, addScaled(x, k3, kDt));

    State next{};
    for (std::size_t i = 0; i < x.size(); ++i) {
        next[i] = x[i] + (kDt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return next;
}
}  // namespace

int main() {
    State x{0.0, 0.0, 1.2, 0.2, 2.0};
    std::ofstream csv("Chapter18_Lesson4_results_cpp.csv");
    if (!csv) {
        throw std::runtime_error("Cannot open output CSV file.");
    }
    csv << "time,q,q_d,error,s,torque,J_hat,b_hat,h_hat\n";
    csv << std::setprecision(12);

    double squaredErrorSum = 0.0;
    std::size_t sampleCount = 0;
    const int steps = static_cast<int>(std::round(kFinalTime / kDt));

    for (int step = 0; step <= steps; ++step) {
        const double t = step * kDt;
        const Desired d = desiredMotion(t);
        const double e = x[0] - d.q;
        const double eDot = x[1] - d.q_dot;
        const double s = eDot + kLambda * e;
        const double qRddot = d.q_ddot - kLambda * eDot;
        const std::array<double, 3> y{qRddot, x[1], std::sin(x[0])};
        const double torque = y[0] * x[2] + y[1] * x[3] + y[2] * x[4] - kKs * s;

        csv << t << ',' << x[0] << ',' << d.q << ',' << e << ',' << s << ',' << torque
            << ',' << x[2] << ',' << x[3] << ',' << x[4] << '\n';
        squaredErrorSum += e * e;
        ++sampleCount;

        if (step < steps) {
            x = rk4Step(t, x);
        }
    }

    const double rmse = std::sqrt(squaredErrorSum / static_cast<double>(sampleCount));
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Tracking RMSE: " << rmse << " rad\n";
    std::cout << "Final estimates [J, b, h]: " << x[2] << ", " << x[3] << ", " << x[4] << '\n';
    std::cout << "Parameter convergence is not guaranteed without persistent excitation.\n";
    return 0;
}

The Java implementation uses primitive arrays and standard file I/O. It follows the same RK4 equations and produces a CSV file for plotting in a spreadsheet, Python, MATLAB, or another analysis tool.

Chapter18_Lesson4.java

// Chapter 18, Lesson 4: adaptive control of a rotary mechanical system.
// Build: javac Chapter18_Lesson4.java
// Run:   java Chapter18_Lesson4

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;

public final class Chapter18_Lesson4 {
    private static final double J_TRUE = 2.5;
    private static final double B_TRUE = 0.8;
    private static final double H_TRUE = 5.0;
    private static final double LAMBDA = 3.0;
    private static final double K_S = 8.0;
    private static final double DT = 0.001;
    private static final double FINAL_TIME = 20.0;
    private static final double[] GAMMA = {2.0, 1.0, 1.5};

    private Chapter18_Lesson4() {}

    private static double[] desiredMotion(double t) {
        return new double[] {
            0.6 * Math.sin(0.5 * t) + 0.2 * Math.sin(1.2 * t),
            0.3 * Math.cos(0.5 * t) + 0.24 * Math.cos(1.2 * t),
            -0.15 * Math.sin(0.5 * t) - 0.288 * Math.sin(1.2 * t)
        };
    }

    private static double[] addScaled(double[] x, double[] dx, double scale) {
        double[] result = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            result[i] = x[i] + scale * dx[i];
        }
        return result;
    }

    private static double[] rhs(double t, double[] x) {
        double q = x[0];
        double qDot = x[1];
        double[] desired = desiredMotion(t);

        double e = q - desired[0];
        double eDot = qDot - desired[1];
        double s = eDot + LAMBDA * e;
        double qRddot = desired[2] - LAMBDA * eDot;

        double[] y = {qRddot, qDot, Math.sin(q)};
        double torque = y[0] * x[2] + y[1] * x[3] + y[2] * x[4] - K_S * s;
        double qDdot = (torque - B_TRUE * qDot - H_TRUE * Math.sin(q)) / J_TRUE;

        double[] dx = new double[5];
        dx[0] = qDot;
        dx[1] = qDdot;
        for (int i = 0; i < 3; i++) {
            dx[i + 2] = -GAMMA[i] * y[i] * s;
        }
        return dx;
    }

    private static double[] rk4Step(double t, double[] x) {
        double[] k1 = rhs(t, x);
        double[] k2 = rhs(t + 0.5 * DT, addScaled(x, k1, 0.5 * DT));
        double[] k3 = rhs(t + 0.5 * DT, addScaled(x, k2, 0.5 * DT));
        double[] k4 = rhs(t + DT, addScaled(x, k3, DT));

        double[] next = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            next[i] = x[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        return next;
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        double[] x = {0.0, 0.0, 1.2, 0.2, 2.0};
        int steps = (int) Math.round(FINAL_TIME / DT);
        double squaredErrorSum = 0.0;
        int sampleCount = 0;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter18_Lesson4_results_java.csv"))) {
            writer.write("time,q,q_d,error,s,torque,J_hat,b_hat,h_hat\n");

            for (int step = 0; step <= steps; step++) {
                double t = step * DT;
                double[] desired = desiredMotion(t);
                double e = x[0] - desired[0];
                double eDot = x[1] - desired[1];
                double s = eDot + LAMBDA * e;
                double qRddot = desired[2] - LAMBDA * eDot;
                double[] y = {qRddot, x[1], Math.sin(x[0])};
                double torque = y[0] * x[2] + y[1] * x[3] + y[2] * x[4] - K_S * s;

                writer.write(String.format(Locale.US,
                    "%.9f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    t, x[0], desired[0], e, s, torque, x[2], x[3], x[4]));
                squaredErrorSum += e * e;
                sampleCount++;

                if (step < steps) {
                    x = rk4Step(t, x);
                }
            }
        }

        double rmse = Math.sqrt(squaredErrorSum / sampleCount);
        System.out.printf(Locale.US, "Tracking RMSE: %.6f rad%n", rmse);
        System.out.printf(Locale.US, "Final estimates [J, b, h]: %.6f, %.6f, %.6f%n",
            x[2], x[3], x[4]);
        System.out.println("Parameter convergence is not guaranteed without persistent excitation.");
    }
}

The MATLAB script performs the numerical simulation and plots the main signals. When Simulink is installed and licensed, the final function also constructs a model containing a state integrator and a MATLAB Function block for the adaptive plant-controller dynamics.

Chapter18_Lesson4.m

%% Chapter18_Lesson4.m
% Adaptive control of a rotary mechanical system with unknown inertia,
% viscous friction, and gravity coefficient.
clear; clc; close all;

J_true = 2.5;
b_true = 0.8;
h_true = 5.0;              % h = m*g*l
lambda = 3.0;
k_s = 8.0;
Gamma = diag([2.0, 1.0, 1.5]);
dt = 0.001;
t_final = 20.0;

time = (0:dt:t_final).';
state = zeros(numel(time), 5);  % q, q_dot, J_hat, b_hat, h_hat
state(1, :) = [0.0, 0.0, 1.2, 0.2, 2.0];
q_d_history = zeros(size(time));
s_history = zeros(size(time));
torque_history = zeros(size(time));

for k = 1:numel(time)
    t = time(k);
    x = state(k, :).';
    [q_d, q_d_dot, q_d_ddot] = desired_motion(t);
    q_d_history(k) = q_d;

    e = x(1) - q_d;
    e_dot = x(2) - q_d_dot;
    s = e_dot + lambda * e;
    q_r_ddot = q_d_ddot - lambda * e_dot;
    Y = [q_r_ddot, x(2), sin(x(1))];
    torque = Y * x(3:5) - k_s * s;

    s_history(k) = s;
    torque_history(k) = torque;

    if k < numel(time)
        k1 = closed_loop_rhs(t, x, J_true, b_true, h_true, lambda, k_s, Gamma);
        k2 = closed_loop_rhs(t + dt/2, x + dt*k1/2, J_true, b_true, h_true, lambda, k_s, Gamma);
        k3 = closed_loop_rhs(t + dt/2, x + dt*k2/2, J_true, b_true, h_true, lambda, k_s, Gamma);
        k4 = closed_loop_rhs(t + dt, x + dt*k3, J_true, b_true, h_true, lambda, k_s, Gamma);
        state(k + 1, :) = (x + dt*(k1 + 2*k2 + 2*k3 + k4)/6).';
    end
end

tracking_error = state(:, 1) - q_d_history;
rmse = sqrt(mean(tracking_error.^2));
fprintf('Tracking RMSE: %.6f rad\n', rmse);
fprintf('Final estimates [J, b, h]: %.6f, %.6f, %.6f\n', state(end, 3:5));
fprintf('Parameter convergence is not guaranteed without persistent excitation.\n');

results = table(time, state(:,1), q_d_history, tracking_error, s_history, ...
    torque_history, state(:,3), state(:,4), state(:,5), ...
    'VariableNames', {'time','q','q_d','error','s','torque','J_hat','b_hat','h_hat'});
writetable(results, 'Chapter18_Lesson4_results_matlab.csv');

figure;
plot(time, q_d_history, 'LineWidth', 1.4); hold on;
plot(time, state(:,1), '--', 'LineWidth', 1.4);
grid on; xlabel('time [s]'); ylabel('angle [rad]');
legend('desired', 'actual', 'Location', 'best');
title('Adaptive tracking with unknown inertia and mass-related gravity');

figure;
plot(time, state(:,3:5), 'LineWidth', 1.2);
grid on; xlabel('time [s]'); ylabel('parameter estimate');
legend('J hat', 'b hat', 'h hat', 'Location', 'best');
title('Online parameter estimates');

% Optional Simulink construction. The model is generated only when Simulink
% is licensed. It provides a reusable shell with a MATLAB Function block that
% evaluates the adaptive controller and estimator derivatives.
if license('test', 'Simulink') && exist('simulink', 'file') == 2
    build_simulink_model();
end

function dx = closed_loop_rhs(t, x, J_true, b_true, h_true, lambda, k_s, Gamma)
    [q_d, q_d_dot, q_d_ddot] = desired_motion(t);
    e = x(1) - q_d;
    e_dot = x(2) - q_d_dot;
    s = e_dot + lambda * e;
    q_r_ddot = q_d_ddot - lambda * e_dot;
    Y = [q_r_ddot, x(2), sin(x(1))];
    torque = Y * x(3:5) - k_s * s;
    q_ddot = (torque - b_true*x(2) - h_true*sin(x(1))) / J_true;
    theta_hat_dot = -Gamma * Y.' * s;
    dx = [x(2); q_ddot; theta_hat_dot];
end

function [q_d, q_d_dot, q_d_ddot] = desired_motion(t)
    q_d = 0.6*sin(0.5*t) + 0.2*sin(1.2*t);
    q_d_dot = 0.3*cos(0.5*t) + 0.24*cos(1.2*t);
    q_d_ddot = -0.15*sin(0.5*t) - 0.288*sin(1.2*t);
end

function build_simulink_model()
    model = 'Chapter18_Lesson4_Simulink';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Clock', [model '/Clock'], ...
        'Position', [40 55 70 85]);
    add_block('simulink/Continuous/Integrator', [model '/State Integrator'], ...
        'InitialCondition', '[0;0;1.2;0.2;2.0]', ...
        'Position', [310 40 350 110]);
    add_block('simulink/User-Defined Functions/MATLAB Function', ...
        [model '/Adaptive Dynamics'], 'Position', [135 35 245 115]);
    add_block('simulink/Sinks/To Workspace', [model '/State Output'], ...
        'VariableName', 'adaptive_state', 'SaveFormat', 'Structure With Time', ...
        'Position', [430 50 520 90]);

    block = [model '/Adaptive Dynamics'];
    script = sprintf([ ...
        'function dx = fcn(t,x)\n' ...
        'J=2.5; b=0.8; h=5.0; lambda=3.0; ks=8.0;\n' ...
        'Gamma=diag([2.0 1.0 1.5]);\n' ...
        'qd=0.6*sin(0.5*t)+0.2*sin(1.2*t);\n' ...
        'qdd=0.3*cos(0.5*t)+0.24*cos(1.2*t);\n' ...
        'qddd=-0.15*sin(0.5*t)-0.288*sin(1.2*t);\n' ...
        'e=x(1)-qd; ed=x(2)-qdd; s=ed+lambda*e;\n' ...
        'qrdd=qddd-lambda*ed; Y=[qrdd x(2) sin(x(1))];\n' ...
        'tau=Y*x(3:5)-ks*s;\n' ...
        'acc=(tau-b*x(2)-h*sin(x(1)))/J;\n' ...
        'thd=-Gamma*Y''*s; dx=[x(2);acc;thd];\n' ...
        'end\n']);
    root = sfroot;
    chart = root.find('-isa', 'Stateflow.EMChart', 'Path', block);
    if isempty(chart)
        error('Unable to locate the MATLAB Function chart.');
    end
    chart.Script = script;

    add_line(model, 'Clock/1', 'Adaptive Dynamics/1');
    add_line(model, 'State Integrator/1', 'Adaptive Dynamics/2');
    add_line(model, 'Adaptive Dynamics/1', 'State Integrator/1');
    add_line(model, 'State Integrator/1', 'State Output/1');
    set_param(model, 'StopTime', '20', 'Solver', 'ode4', 'FixedStep', '0.001');
    save_system(model);
    fprintf('Created %s.slx\n', model);
end

The Mathematica notebook expresses the coupled differential equations symbolically and solves them with NDSolve. This form is useful for checking derivative definitions and experimenting with alternative desired trajectories.

Chapter18_Lesson4.nb

Notebook[{
  Cell[
    "Chapter 18, Lesson 4: Adaptive Control of a Mechanical System",
    "Title"
  ],

  Cell[
    StringJoin[
      "The rotary plant has unknown inertia J, viscous coefficient b, ",
      "and gravity coefficient h = m g l. The adaptive controller uses ",
      "a filtered tracking error and a linear regressor."
    ],
    "Text"
  ],

  Cell[
    BoxData[
      StringJoin[
        "ClearAll[\"Global`*\"];\n\n",
        "jTrue = 2.5;\nbTrue = 0.8;\nhTrue = 5.0;\n",
        "lambda = 3.0;\nks = 8.0;\n",
        "gamma = DiagonalMatrix[{2.0, 1.0, 1.5}];\n\n",
        "qd[t_] := 0.6 Sin[0.5 t] + 0.2 Sin[1.2 t];\n",
        "qd1[t_] := 0.3 Cos[0.5 t] + 0.24 Cos[1.2 t];\n",
        "qd2[t_] := -0.15 Sin[0.5 t] - 0.288 Sin[1.2 t];\n\n",
        "e[t_] := q[t] - qd[t];\n",
        "ed[t_] := q'[t] - qd1[t];\n",
        "s[t_] := ed[t] + lambda e[t];\n",
        "qrdd[t_] := qd2[t] - lambda ed[t];\n",
        "y[t_] := {qrdd[t], q'[t], Sin[q[t]]};\n",
        "thetaHat[t_] := {jHat[t], bHat[t], hHat[t]};\n",
        "tau[t_] := y[t].thetaHat[t] - ks s[t];\n\n",
        "sol = First@NDSolve[\n",
        "  {\n",
        "    jTrue q''[t] + bTrue q'[t] + hTrue Sin[q[t]] == tau[t],\n",
        "    jHat'[t] == -gamma[[1, 1]] y[t][[1]] s[t],\n",
        "    bHat'[t] == -gamma[[2, 2]] y[t][[2]] s[t],\n",
        "    hHat'[t] == -gamma[[3, 3]] y[t][[3]] s[t],\n",
        "    q[0] == 0,\n    q'[0] == 0,\n",
        "    jHat[0] == 1.2,\n    bHat[0] == 0.2,\n",
        "    hHat[0] == 2.0\n",
        "  },\n",
        "  {q, jHat, bHat, hHat},\n",
        "  {t, 0, 20},\n",
        "  MaxStepSize -> 0.001\n",
        "];\n\n",
        "trackingPlot = Plot[\n",
        "  Evaluate[{q[t], qd[t]} /. sol],\n",
        "  {t, 0, 20},\n",
        "  PlotLegends -> {\"actual\", \"desired\"},\n",
        "  AxesLabel -> {\"time [s]\", \"angle [rad]\"},\n",
        "  PlotRange -> All\n",
        "];\n\n",
        "parameterPlot = Plot[\n",
        "  Evaluate[{jHat[t], bHat[t], hHat[t]} /. sol],\n",
        "  {t, 0, 20},\n",
        "  PlotLegends -> {\"J hat\", \"b hat\", \"h hat\"},\n",
        "  AxesLabel -> {\"time [s]\", \"estimate\"},\n",
        "  PlotRange -> All\n",
        "];\n\n",
        "Column[{trackingPlot, parameterPlot}]\n"
      ]
    ],
    "Input"
  ]
},
WindowSize -> {1200, 800},
StyleDefinitions -> "Default.nb"
]        

14.1 Tracking gain

Increasing \(k_s\) strengthens damping in the filtered error equation. It usually reduces tracking error but demands greater torque and can amplify measurement noise in practical implementations.

14.2 Adaptation gain

Larger entries of \(\Gamma\) produce faster parameter motion. They do not directly strengthen the Lyapunov dissipation term, which remains \(-s^{T}K_Ds\). Excessively large adaptation gains can create rapid estimate variation under sampled, noisy, or unmodeled dynamics.

14.3 Numerical integration

The continuous-time proof assumes exact differential equations. A digital implementation approximates both plant and adaptation dynamics. The sampling interval must resolve the fastest closed-loop mode, and actuator saturation must be handled explicitly because the ideal cancellation no longer holds while the commanded torque is clipped.

14.4 Acceleration measurement is not required

The selected regressor uses the known reference acceleration \(\ddot q_r\), not the measured joint acceleration \(\ddot q\). This is an important practical advantage of the filtered-error construction.

14.5 Structural knowledge remains essential

The adaptive law does not discover arbitrary mechanics. The designer must know which basis functions belong in \(Y\). Omitting a significant gravity, friction, or coupling term changes the problem from exact parametric uncertainty to model mismatch.

Problem 1 (Unknown Payload Mass): A vertical actuator is modeled by \(m\ddot q+c\dot q+mg=u\), where \(m\) is unknown and \(c\) is known. Construct a scalar regressor, adaptive controller, and update law.

Solution:

Define \(s=\dot e+\lambda e\) and \(\ddot q_r=\ddot q_d-\lambda\dot e\). Then

\[ m\ddot q_r+mg=m(\ddot q_r+g)=Y_m m, \qquad Y_m=\ddot q_r+g. \]

Use

\[ u=\hat m(\ddot q_r+g)+c\dot q-k_s s, \qquad \dot{\hat m}=-\gamma_m(\ddot q_r+g)s. \]

The closed-loop equation is \(m\dot s+k_s s=Y_m\tilde m\). With \(V=\frac{1}{2}ms^2+\frac{1}{2\gamma_m}\tilde m^2\), one obtains \(\dot V=-k_s s^2\).

Problem 2 (Sign Error in the Update Law): Suppose the estimation error is defined as \(\tilde\theta=\hat\theta-\theta\), but the update law is chosen as \(\dot{\hat\theta}=+\Gamma Y^{T}s\). Determine the Lyapunov derivative.

Solution:

Starting from \(J\dot s=Y\tilde\theta-k_s s\),

\[ \dot V=sY\tilde\theta-k_s s^2+ \tilde\theta^{T}Y^{T}s =2sY\tilde\theta-k_s s^2. \]

The cross term doubles instead of cancelling, so negative semidefiniteness is lost. The sign of the update law must be consistent with the definition of the estimation error.

Problem 3 (Skew-Symmetry Cancellation): Prove that if \(A^{T}=-A\), then \(x^{T}Ax=0\) for every real vector \(x\).

Solution:

The scalar \(x^{T}Ax\) equals its transpose:

\[ x^{T}Ax=(x^{T}Ax)^{T}=x^{T}A^{T}x=-x^{T}Ax. \]

Therefore \(2x^{T}Ax=0\), which implies \(x^{T}Ax=0\). Taking \(A=\dot M-2C\) removes the kinetic cross term in the robot Lyapunov derivative.

Problem 4 (Tracking Without Identification): A rotary system tracks a constant setpoint after a transient. Explain why the estimates of \(J\) and \(b\) may stop far from their true values.

Solution:

At a constant setpoint, the reference acceleration and velocity approach zero. Consequently the inertia and damping columns of \(Y=[\ddot q_r,\dot q,\sin q]\) approach zero. The update signals for \(\hat J\) and \(\hat b\) vanish, so those estimates freeze. The feedback term can still drive the tracking error to zero. The trajectory therefore supplies insufficient persistent excitation for parameter identification.

Problem 5 (n-DOF Lyapunov Derivative): For \(M\dot s+Cs+K_Ds=Y\tilde\theta\), derive \(\dot V\) for

\[ V=\frac{1}{2}s^{T}Ms+\frac{1}{2} \tilde\theta^{T}\Gamma^{-1}\tilde\theta \]

under \(\dot{\hat\theta}=-\Gamma Y^{T}s\).

Solution:

\[ \begin{aligned} \dot V &=s^{T}(Y\tilde\theta-Cs-K_Ds) +\frac{1}{2}s^{T}\dot M s -\tilde\theta^{T}Y^{T}s\\ &=-s^{T}K_Ds+\frac{1}{2}s^{T}(\dot M-2C)s\\ &=-s^{T}K_Ds. \end{aligned} \]

The parameter cross terms cancel by the adaptive law, and the mechanical cross term vanishes by skew symmetry.

Problem 6 (Effect of an Incorrect Regressor): The true rotary plant contains Coulomb friction \(f_c\operatorname{sgn}(\dot q)\), but the controller regressor contains only inertia, viscous friction, and gravity terms. Write the modified filtered-error dynamics and explain the implication.

Solution:

The true plant is

\[ J\ddot q+b\dot q+h\sin q+f_c\operatorname{sgn}(\dot q)=\tau. \]

Therefore

\[ J\dot s+k_s s=Y\tilde\theta-f_c\operatorname{sgn}(\dot q). \]

The Lyapunov derivative becomes

\[ \dot V=-k_s s^2-f_c s\operatorname{sgn}(\dot q), \]

whose second term is not guaranteed to be nonpositive. The ideal proof no longer establishes asymptotic tracking. The missing friction term should be included in the regressor when its structure is credible, or the controller should be robustified using methods introduced earlier in the course.

Mechanical systems may be nonlinear in their states while remaining linear in unknown mass, inertia, damping, gravity, and other physical coefficients. By introducing a filtered tracking error, expressing the reference dynamics as \(Y\theta\), and combining certainty-equivalent feedforward with damping feedback, one obtains an adaptive controller whose parameter cross term cancels exactly in a Lyapunov derivative. For multi-DOF rigid mechanics, the additional cancellation follows from the skew-symmetry of \(\dot M-2C\). The resulting ideal controller guarantees bounded estimates and asymptotic tracking, but exact mass or inertia convergence requires sufficiently rich excitation and a correct mechanical regressor.