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

• express a simple mechanical plant in a linear-in-the-parameters form;
• derive a tracking-error-only adaptive controller and identify what its Lyapunov proof does and does not guarantee;
• construct a prediction regression without directly measuring acceleration;
• form an information matrix and a measurable prediction residual;
• derive a composite adaptive law using tracking and prediction information;
• prove boundedness and asymptotic tracking, and state an excitation condition for parameter convergence; and
• reproduce the numerical comparison in five programming environments.

The example deliberately uses a one-degree-of-freedom mechanical model. This keeps the algebra visible while preserving the same structural ingredients used in multi-degree-of-freedom Euler-Lagrange composite adaptive controllers.

Consider a translating mechanical body with unknown mass \( m \) and unknown viscous damping \( b \). Its measured velocity is \( v(t) \), and the applied force is \( u(t) \):

\[ m\dot v + bv = u, \qquad m > 0, \qquad b \ge 0. \]

Define the constant parameter vector and the plant regression:

\[ \boldsymbol{\theta} = \begin{bmatrix}m\\b\end{bmatrix}, \qquad u = \underbrace{\begin{bmatrix}\dot v & v\end{bmatrix}}_ {\mathbf{Y}(v,\dot v)} \boldsymbol{\theta}. \]

Let the desired velocity \( v_d(t) \) and desired acceleration \( \dot v_d(t) \) be bounded and available. The tracking error is

\[ e = v-v_d. \]

For the numerical experiment, the reference contains two frequencies:

\[ v_d(t) = 0.8\sin(0.7t)+0.35\sin(1.9t), \] \[ \dot v_d(t) = 0.56\cos(0.7t)+0.665\cos(1.9t). \]

The two-frequency signal is selected to provide richer parameter information than a constant set point. A constant or single-mode trajectory can still yield good tracking while leaving some parameter directions weakly observable.

Let \( \hat{\boldsymbol{\theta}} = [\hat m,\hat b]^T \) denote the parameter estimate. Choose the control law

\[ u = \hat m\dot v_d + \hat b v - ke, \qquad k > 0. \]

Introduce the tracking regressor and parameter error:

\[ \boldsymbol{\phi} = \begin{bmatrix}\dot v_d\\v\end{bmatrix}, \qquad \tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}} - \boldsymbol{\theta}. \]

Substituting the controller into the plant and using \( \dot e=\dot v-\dot v_d \) gives

\[ \begin{aligned} m\dot e &= u-bv-m\dot v_d\\ &= -ke + (\hat m-m)\dot v_d + (\hat b-b)v\\ &= -ke + \boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}}. \end{aligned} \]

Thus the tracking subsystem is stable except for a term generated by parameter error.

A conventional Lyapunov-gradient update uses only the instantaneous tracking error:

\[ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma} \boldsymbol{\phi}e, \qquad \boldsymbol{\Gamma} = \boldsymbol{\Gamma}^T >0. \]

Consider the Lyapunov function

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

Since the true parameters are constant, \( \dot{\tilde{\boldsymbol{\theta}}} = \dot{\hat{\boldsymbol{\theta}}} \) . Therefore,

\[ \begin{aligned} \dot V &= me\dot e + \tilde{\boldsymbol{\theta}}^T \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}\\ &= -ke^2 + e\boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}} - \tilde{\boldsymbol{\theta}}^T \boldsymbol{\phi}e\\ &= -ke^2. \end{aligned} \]

This proves boundedness and, under standard boundedness conditions, asymptotic tracking. It does not by itself force \( \tilde{\boldsymbol{\theta}} \) to zero. Once \( e \) becomes small, the parameter update may lose its driving signal. Composite adaptation addresses this limitation by adding a second source of parameter information.

Direct use of \( u=m\dot v+bv \) in a predictor would require measured acceleration. Instead, apply the stable first-order filter

\[ H(s)=\frac{\lambda}{s+\lambda}, \qquad \lambda >0. \]

Generate filtered velocity and filtered input through

\[ \dot v_f=-\lambda v_f+\lambda v, \qquad \dot u_f=-\lambda u_f+\lambda u. \]

Define the filtered derivative signal

\[ p_f = \lambda(v-v_f). \]

If the filter is initialized consistently, for example \( v_f(0)=v(0) \), then \( p_f=H(s)[\dot v] \). Filtering the plant equation yields the measurable regression

\[ u_f = mp_f+bv_f = \mathbf{Y}_f\boldsymbol{\theta}, \qquad \mathbf{Y}_f = \begin{bmatrix}p_f & v_f\end{bmatrix}. \]

With inconsistent initial filter states, an exponentially decaying filter-transient term must be included. The simulations use \( v(0)=v_f(0)=0 \), so the exact filtered regression is obtained.

Accumulate filtered regression information:

\[ \mathbf{R}(t) = \int_0^t \mathbf{Y}_f^T(\tau) \mathbf{Y}_f(\tau) \,d\tau, \] \[ \mathbf{q}(t) = \int_0^t \mathbf{Y}_f^T(\tau) u_f(\tau) \,d\tau. \]

Since \( u_f=\mathbf{Y}_f\boldsymbol{\theta} \),

\[ \mathbf{q} = \int_0^t \mathbf{Y}_f^T \mathbf{Y}_f \boldsymbol{\theta} \,d\tau = \mathbf{R}\boldsymbol{\theta}. \]

To avoid a residual whose numerical magnitude grows with the amount of stored data, introduce the normalized quantities

\[ \bar{\mathbf{R}} = \frac{\mathbf{R}}{\rho_0+\operatorname{tr}(\mathbf{R})}, \qquad \bar{\mathbf{q}} = \frac{\mathbf{q}}{\rho_0+\operatorname{tr}(\mathbf{R})}, \qquad \rho_0 >0. \]

The composite prediction residual is

\[ \boldsymbol{\zeta} = \bar{\mathbf{R}} \hat{\boldsymbol{\theta}} - \bar{\mathbf{q}} = \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}}. \]

Unlike the unavailable parameter error, this residual is computable from the estimated parameters and filtered input-output data.

flowchart TD
  RD["Reference: vd and dvd"] --> CT["Controller: u = mhat*dvd + bhat*v - k*e"]
  CT --> PL["Plant: m*vdot + b*v = u"]
  PL --> VE["Measured velocity v"]
  VE --> ER["Tracking error e = v - vd"]
  ER --> AD["Composite parameter update"]
  VE --> VF["Stable velocity filter"]
  CT --> UF["Stable input filter"]
  VF --> RG["Filtered regressor Yf"]
  UF --> IR["Information vector q"]
  RG --> IM["Information matrix R"]
  IM --> PR["Prediction residual zeta"]
  IR --> PR
  PR --> AD
  AD --> CT
        

Combine the tracking term and prediction term:

\[ \boxed{ \dot{\hat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma} \left( \boldsymbol{\phi}e + k_c\boldsymbol{\zeta} \right) } \]

where \( k_c\ge0 \) is the composite learning gain. In expanded form,

\[ \begin{bmatrix} \dot{\hat m}\\ \dot{\hat b} \end{bmatrix} = - \begin{bmatrix} \gamma_m & 0\\ 0 & \gamma_b \end{bmatrix} \left[ \begin{bmatrix} \dot v_d\\ v \end{bmatrix}e + k_c \left( \bar{\mathbf{R}} \begin{bmatrix} \hat m\\ \hat b \end{bmatrix} - \bar{\mathbf{q}} \right) \right]. \]

Setting \( k_c=0 \) recovers the tracking-error-only adaptive controller. The controller structure itself is unchanged; only the update law receives an additional identification channel.

Use the same Lyapunov candidate:

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

Substituting the error dynamics and composite update law gives

\[ \begin{aligned} \dot V &= e\left( -ke + \boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}} \right) + \tilde{\boldsymbol{\theta}}^T \boldsymbol{\Gamma}^{-1} \left[ -\boldsymbol{\Gamma} \left( \boldsymbol{\phi}e + k_c\bar{\mathbf{R}} \tilde{\boldsymbol{\theta}} \right) \right]\\ &= -ke^2 + e\boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}} - \tilde{\boldsymbol{\theta}}^T \boldsymbol{\phi}e - k_c \tilde{\boldsymbol{\theta}}^T \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}}\\ &= \boxed{ -ke^2 - k_c \tilde{\boldsymbol{\theta}}^T \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}} } \le0. \end{aligned} \]

The information matrix is positive semidefinite because

\[ \mathbf{x}^T\mathbf{R}\mathbf{x} = \int_0^t \left( \mathbf{Y}_f\mathbf{x} \right)^2d\tau \ge0 \qquad \text{for all }\mathbf{x}. \]

Therefore all signals appearing in \( V \) remain bounded. If the reference and its derivative are bounded, then the closed-loop equations imply bounded \( \dot e \). Since \( e\in L_2\cap L_\infty \), Barbalat's lemma gives

\[ \lim_{t\rightarrow\infty}e(t)=0. \]

Tracking convergence does not automatically imply parameter convergence. Suppose that after an information-collection interval there exist \( t_e \) and \( \alpha_R >0 \) such that

\[ \bar{\mathbf{R}}(t) \ge \alpha_R\mathbf{I}, \qquad t\ge t_e. \]

Then

\[ \dot V \le -ke^2 - k_c\alpha_R \left\| \tilde{\boldsymbol{\theta}} \right\|^2. \]

Because \( V \) is quadratically equivalent to \( e^2+ \|\tilde{\boldsymbol{\theta}}\|^2 \) , there exists \( c >0 \) such that

\[ \dot V\le-cV, \qquad t\ge t_e, \] \[ V(t) \le V(t_e)e^{-c(t-t_e)}. \]

Hence both tracking and parameter errors converge exponentially once the accumulated data matrix is uniformly positive definite. For the two-parameter mass-damper model, this requires the filtered regressor directions to span two independent directions over time. Merely making the input nonzero is not enough; the data must be informative with respect to both \( m \) and \( b \).

The condition above is closely related to finite or interval excitation. It differs from requiring persistent excitation forever: with a non-forgetting memory matrix, sufficiently rich information collected during a finite interval remains available to the update law.

At every integration or control update, execute the following sequence. In a digital controller, the differential equations are replaced by a numerically stable discretization.

flowchart TD
  S["Read time t and velocity v"] --> R["Compute vd and dvd"]
  R --> E["Compute tracking error e"]
  E --> U["Compute control force u"]
  U --> P["Advance physical plant or apply u to hardware"]
  P --> F["Update stable filters vf and uf"]
  F --> Y["Compute filtered derivative p and regressor Yf"]
  Y --> M["Update information matrix R and vector q"]
  M --> Z["Compute normalized prediction residual zeta"]
  Z --> A["Update mhat and bhat"]
  A --> S
        

1. Compute \( e=v-v_d \), \( \boldsymbol{\phi}=[\dot v_d,v]^T \).
2. Apply \( u=\hat m\dot v_d+\hat bv-ke \).
3. Update \( v_f \) and \( u_f \).
4. Compute \( p_f=\lambda(v-v_f) \) and \( \mathbf{Y}_f=[p_f,v_f] \).
5. Update \( \dot{\mathbf{R}} =\mathbf{Y}_f^T\mathbf{Y}_f \) and \( \dot{\mathbf{q}} =\mathbf{Y}_f^Tu_f \).
6. Form \( \boldsymbol{\zeta} = (\mathbf{R}\hat{\boldsymbol{\theta}}-\mathbf{q}) /(\rho_0+\operatorname{tr}\mathbf{R}) \).
7. Update the estimate with the composite law.

The simulation uses

\[ m=2.5, \qquad b=1.2, \qquad k=6, \] \[ \boldsymbol{\Gamma} = \operatorname{diag}(4,4), \qquad k_c=2, \qquad \lambda=8, \] \[ \hat m(0)=1, \qquad \hat b(0)=0.2. \]

The software runs both controllers with the same fourth-order Runge-Kutta integrator. The tracking-only case is obtained by setting \( k_c=0 \); the composite case uses \( k_c=2 \).

With the supplied parameters and a step size of \( 10^{-3} \) s, the reference implementation produces approximately:

• tracking-only RMSE: \( 5.82\times10^{-2} \) m/s;
• composite RMSE: \( 3.81\times10^{-2} \) m/s;
• tracking-only terminal mass and damping errors of roughly \( 1.76\times10^{-2} \) each; and
• composite terminal parameter errors close to numerical precision for this ideal noise-free experiment.

These numbers are not universal performance guarantees. They describe this model, reference, gain set, initialization, and numerical method. The theoretical conclusion is the sign of \( \dot V \) and the role of information-matrix rank, not a fixed percentage improvement.

Chapter21_Lesson5.py

Libraries: numpy for vector and matrix calculations and matplotlib for plots. The implementation uses an explicit RK4 integrator so the adaptive equations remain visible.

"""
Chapter21_Lesson5.py
Composite adaptive control of an uncertain mass-damper velocity servo.

The program compares:
1. Tracking-error-only adaptive control.
2. Composite adaptive control using tracking error plus a filtered-data
   prediction residual.

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Configuration:
    mass_true: float = 2.5
    damping_true: float = 1.2
    feedback_gain: float = 6.0
    gamma_mass: float = 4.0
    gamma_damping: float = 4.0
    composite_gain: float = 2.0
    filter_rate: float = 8.0
    dt: float = 0.001
    final_time: float = 20.0
    sample_stride: int = 10


CFG = Configuration()


def reference(t: float) -> tuple[float, float]:
    """Return desired velocity and its derivative."""
    velocity = 0.8 * np.sin(0.7 * t) + 0.35 * np.sin(1.9 * t)
    acceleration = 0.56 * np.cos(0.7 * t) + 0.665 * np.cos(1.9 * t)
    return float(velocity), float(acceleration)


def derivative(
    t: float,
    state: np.ndarray,
    composite_gain: float,
    cfg: Configuration = CFG,
) -> np.ndarray:
    """
    State order:
      [v, v_f, u_f, m_hat, b_hat, R11, R12, R22, q1, q2]
    """
    (
        velocity,
        velocity_filtered,
        input_filtered,
        mass_hat,
        damping_hat,
        r11,
        r12,
        r22,
        q1,
        q2,
    ) = state

    desired_velocity, desired_acceleration = reference(t)
    tracking_error = velocity - desired_velocity

    control_input = (
        mass_hat * desired_acceleration
        + damping_hat * velocity
        - cfg.feedback_gain * tracking_error
    )
    velocity_dot = (
        control_input - cfg.damping_true * velocity
    ) / cfg.mass_true

    velocity_filtered_dot = (
        -cfg.filter_rate * velocity_filtered
        + cfg.filter_rate * velocity
    )
    input_filtered_dot = (
        -cfg.filter_rate * input_filtered
        + cfg.filter_rate * control_input
    )

    filtered_acceleration = cfg.filter_rate * (
        velocity - velocity_filtered
    )
    tracking_regressor = np.array(
        [desired_acceleration, velocity], dtype=float
    )
    prediction_regressor = np.array(
        [filtered_acceleration, velocity_filtered], dtype=float
    )

    information_matrix = np.array(
        [[r11, r12], [r12, r22]], dtype=float
    )
    information_vector = np.array([q1, q2], dtype=float)
    parameter_estimate = np.array(
        [mass_hat, damping_hat], dtype=float
    )

    prediction_residual = (
        information_matrix @ parameter_estimate - information_vector
    )
    normalization = 1.0 + float(np.trace(information_matrix))
    gamma = np.diag([cfg.gamma_mass, cfg.gamma_damping])

    estimate_dot = -gamma @ (
        tracking_regressor * tracking_error
        + composite_gain * prediction_residual / normalization
    )

    information_matrix_dot = np.outer(
        prediction_regressor, prediction_regressor
    )
    information_vector_dot = (
        prediction_regressor * input_filtered
    )

    return np.array(
        [
            velocity_dot,
            velocity_filtered_dot,
            input_filtered_dot,
            estimate_dot[0],
            estimate_dot[1],
            information_matrix_dot[0, 0],
            information_matrix_dot[0, 1],
            information_matrix_dot[1, 1],
            information_vector_dot[0],
            information_vector_dot[1],
        ],
        dtype=float,
    )


def rk4_step(
    rhs: Callable[[float, np.ndarray, float], np.ndarray],
    t: float,
    state: np.ndarray,
    dt: float,
    composite_gain: float,
) -> np.ndarray:
    """One classical fourth-order Runge-Kutta step."""
    k1 = rhs(t, state, composite_gain)
    k2 = rhs(t + 0.5 * dt, state + 0.5 * dt * k1, composite_gain)
    k3 = rhs(t + 0.5 * dt, state + 0.5 * dt * k2, composite_gain)
    k4 = rhs(t + dt, state + dt * k3, composite_gain)
    return state + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0


def simulate(
    composite: bool,
    cfg: Configuration = CFG,
) -> dict[str, np.ndarray]:
    """Simulate one adaptive controller."""
    gain = cfg.composite_gain if composite else 0.0
    state = np.array(
        [0.0, 0.0, 0.0, 1.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0],
        dtype=float,
    )
    steps = int(round(cfg.final_time / cfg.dt))

    records: list[list[float]] = []
    for step in range(steps + 1):
        t = step * cfg.dt
        if step % cfg.sample_stride == 0:
            desired_velocity, desired_acceleration = reference(t)
            velocity = state[0]
            error = velocity - desired_velocity
            mass_hat = state[3]
            damping_hat = state[4]
            control_input = (
                mass_hat * desired_acceleration
                + damping_hat * velocity
                - cfg.feedback_gain * error
            )
            information_matrix = np.array(
                [[state[5], state[6]], [state[6], state[7]]],
                dtype=float,
            )
            minimum_eigenvalue = float(
                np.linalg.eigvalsh(information_matrix)[0]
            )
            records.append(
                [
                    t,
                    desired_velocity,
                    velocity,
                    error,
                    control_input,
                    mass_hat,
                    damping_hat,
                    minimum_eigenvalue,
                ]
            )

        if step < steps:
            state = rk4_step(
                derivative, t, state, cfg.dt, gain
            )

    data = np.asarray(records, dtype=float)
    keys = (
        "time",
        "desired_velocity",
        "velocity",
        "tracking_error",
        "control_input",
        "mass_hat",
        "damping_hat",
        "information_min_eigenvalue",
    )
    return {key: data[:, index] for index, key in enumerate(keys)}


def metrics(
    data: dict[str, np.ndarray],
    cfg: Configuration = CFG,
) -> dict[str, float]:
    """Compute tracking and terminal parameter metrics."""
    return {
        "tracking_rmse": float(
            np.sqrt(np.mean(data["tracking_error"] ** 2))
        ),
        "terminal_mass_error": float(
            abs(data["mass_hat"][-1] - cfg.mass_true)
        ),
        "terminal_damping_error": float(
            abs(data["damping_hat"][-1] - cfg.damping_true)
        ),
    }


def save_csv(
    direct: dict[str, np.ndarray],
    composite: dict[str, np.ndarray],
    destination: Path,
) -> None:
    """Save aligned comparison data."""
    matrix = np.column_stack(
        [
            direct["time"],
            direct["desired_velocity"],
            direct["velocity"],
            composite["velocity"],
            direct["tracking_error"],
            composite["tracking_error"],
            direct["mass_hat"],
            composite["mass_hat"],
            direct["damping_hat"],
            composite["damping_hat"],
            composite["information_min_eigenvalue"],
        ]
    )
    header = (
        "time,desired_velocity,direct_velocity,composite_velocity,"
        "direct_error,composite_error,direct_mass_hat,"
        "composite_mass_hat,direct_damping_hat,"
        "composite_damping_hat,information_min_eigenvalue"
    )
    np.savetxt(
        destination,
        matrix,
        delimiter=",",
        header=header,
        comments="",
    )


def create_plots(
    direct: dict[str, np.ndarray],
    composite: dict[str, np.ndarray],
    destination: Path,
) -> None:
    """Create a four-panel comparison figure."""
    time = direct["time"]
    figure, axes = plt.subplots(2, 2, figsize=(12, 8))

    axes[0, 0].plot(
        time, direct["desired_velocity"], label="reference"
    )
    axes[0, 0].plot(
        time, direct["velocity"], label="tracking-only"
    )
    axes[0, 0].plot(
        time, composite["velocity"], label="composite"
    )
    axes[0, 0].set_title("Velocity tracking")
    axes[0, 0].set_xlabel("Time [s]")
    axes[0, 0].set_ylabel("Velocity [m/s]")
    axes[0, 0].grid(True)
    axes[0, 0].legend()

    axes[0, 1].plot(
        time, direct["tracking_error"], label="tracking-only"
    )
    axes[0, 1].plot(
        time, composite["tracking_error"], label="composite"
    )
    axes[0, 1].set_title("Tracking error")
    axes[0, 1].set_xlabel("Time [s]")
    axes[0, 1].set_ylabel("Error [m/s]")
    axes[0, 1].grid(True)
    axes[0, 1].legend()

    axes[1, 0].axhline(
        CFG.mass_true, linestyle="--", label="true mass"
    )
    axes[1, 0].plot(
        time, direct["mass_hat"], label="tracking-only estimate"
    )
    axes[1, 0].plot(
        time, composite["mass_hat"], label="composite estimate"
    )
    axes[1, 0].set_title("Mass estimate")
    axes[1, 0].set_xlabel("Time [s]")
    axes[1, 0].set_ylabel("Mass [kg]")
    axes[1, 0].grid(True)
    axes[1, 0].legend()

    axes[1, 1].axhline(
        CFG.damping_true, linestyle="--", label="true damping"
    )
    axes[1, 1].plot(
        time,
        direct["damping_hat"],
        label="tracking-only estimate",
    )
    axes[1, 1].plot(
        time,
        composite["damping_hat"],
        label="composite estimate",
    )
    axes[1, 1].set_title("Damping estimate")
    axes[1, 1].set_xlabel("Time [s]")
    axes[1, 1].set_ylabel("Damping [N s/m]")
    axes[1, 1].grid(True)
    axes[1, 1].legend()

    figure.tight_layout()
    figure.savefig(destination, dpi=180)
    plt.close(figure)


def main() -> None:
    output_directory = Path(__file__).resolve().parent
    direct = simulate(composite=False)
    composite = simulate(composite=True)

    direct_metrics = metrics(direct)
    composite_metrics = metrics(composite)

    print("Tracking-error-only controller:")
    for name, value in direct_metrics.items():
        print(f"  {name}: {value:.6f}")

    print("Composite controller:")
    for name, value in composite_metrics.items():
        print(f"  {name}: {value:.6f}")

    csv_path = output_directory / "Chapter21_Lesson5_results_python.csv"
    plot_path = output_directory / "Chapter21_Lesson5_python.png"
    save_csv(direct, composite, csv_path)
    create_plots(direct, composite, plot_path)
    print(f"Saved {csv_path.name}")
    print(f"Saved {plot_path.name}")


if __name__ == "__main__":
    main()

Chapter21_Lesson5.cpp

The C++17 version uses only the standard library. It writes a CSV file for plotting in a separate tool.

/*
Chapter21_Lesson5.cpp
Composite adaptive control of an uncertain mass-damper velocity servo.

Build:
    g++ -std=c++17 -O2 -Wall -Wextra -pedantic Chapter21_Lesson5.cpp -o Chapter21_Lesson5

Run:
    ./Chapter21_Lesson5
*/

#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

namespace {

constexpr double MASS_TRUE = 2.5;
constexpr double DAMPING_TRUE = 1.2;
constexpr double FEEDBACK_GAIN = 6.0;
constexpr double GAMMA_MASS = 4.0;
constexpr double GAMMA_DAMPING = 4.0;
constexpr double COMPOSITE_GAIN = 2.0;
constexpr double FILTER_RATE = 8.0;
constexpr double DT = 0.001;
constexpr double FINAL_TIME = 20.0;
constexpr int SAMPLE_STRIDE = 10;

using State = std::array<double, 10>;

struct Sample {
    double time{};
    double desired_velocity{};
    double velocity{};
    double tracking_error{};
    double control_input{};
    double mass_hat{};
    double damping_hat{};
    double information_min_eigenvalue{};
};

struct Metrics {
    double tracking_rmse{};
    double terminal_mass_error{};
    double terminal_damping_error{};
};

std::pair<double, double> reference(const double t) {
    const double velocity =
        0.8 * std::sin(0.7 * t) + 0.35 * std::sin(1.9 * t);
    const double acceleration =
        0.56 * std::cos(0.7 * t) + 0.665 * std::cos(1.9 * t);
    return {velocity, acceleration};
}

double minimumEigenvalueSymmetric2x2(
    const double a,
    const double b,
    const double d
) {
    const double trace = a + d;
    const double discriminant =
        std::sqrt((a - d) * (a - d) + 4.0 * b * b);
    return 0.5 * (trace - discriminant);
}

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

State derivative(
    const double t,
    const State& state,
    const double composite_gain
) {
    const double velocity = state[0];
    const double velocity_filtered = state[1];
    const double input_filtered = state[2];
    const double mass_hat = state[3];
    const double damping_hat = state[4];
    const double r11 = state[5];
    const double r12 = state[6];
    const double r22 = state[7];
    const double q1 = state[8];
    const double q2 = state[9];

    const auto [desired_velocity, desired_acceleration] = reference(t);
    const double tracking_error = velocity - desired_velocity;

    const double control_input =
        mass_hat * desired_acceleration
        + damping_hat * velocity
        - FEEDBACK_GAIN * tracking_error;

    const double velocity_dot =
        (control_input - DAMPING_TRUE * velocity) / MASS_TRUE;
    const double velocity_filtered_dot =
        -FILTER_RATE * velocity_filtered + FILTER_RATE * velocity;
    const double input_filtered_dot =
        -FILTER_RATE * input_filtered + FILTER_RATE * control_input;

    const double filtered_acceleration =
        FILTER_RATE * (velocity - velocity_filtered);
    const double tracking_phi_1 = desired_acceleration;
    const double tracking_phi_2 = velocity;
    const double prediction_y_1 = filtered_acceleration;
    const double prediction_y_2 = velocity_filtered;

    const double residual_1 =
        r11 * mass_hat + r12 * damping_hat - q1;
    const double residual_2 =
        r12 * mass_hat + r22 * damping_hat - q2;
    const double normalization = 1.0 + r11 + r22;

    const double mass_hat_dot =
        -GAMMA_MASS
        * (
            tracking_phi_1 * tracking_error
            + composite_gain * residual_1 / normalization
        );
    const double damping_hat_dot =
        -GAMMA_DAMPING
        * (
            tracking_phi_2 * tracking_error
            + composite_gain * residual_2 / normalization
        );

    const double r11_dot = prediction_y_1 * prediction_y_1;
    const double r12_dot = prediction_y_1 * prediction_y_2;
    const double r22_dot = prediction_y_2 * prediction_y_2;
    const double q1_dot = prediction_y_1 * input_filtered;
    const double q2_dot = prediction_y_2 * input_filtered;

    return {
        velocity_dot,
        velocity_filtered_dot,
        input_filtered_dot,
        mass_hat_dot,
        damping_hat_dot,
        r11_dot,
        r12_dot,
        r22_dot,
        q1_dot,
        q2_dot
    };
}

State rk4Step(
    const double t,
    const State& state,
    const double composite_gain
) {
    const State k1 = derivative(t, state, composite_gain);
    const State k2 = derivative(
        t + 0.5 * DT,
        addScaled(state, k1, 0.5 * DT),
        composite_gain
    );
    const State k3 = derivative(
        t + 0.5 * DT,
        addScaled(state, k2, 0.5 * DT),
        composite_gain
    );
    const State k4 = derivative(
        t + DT,
        addScaled(state, k3, DT),
        composite_gain
    );

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

std::vector<Sample> simulate(const bool composite) {
    const double composite_gain =
        composite ? COMPOSITE_GAIN : 0.0;
    State state{
        0.0, 0.0, 0.0, 1.0, 0.2,
        0.0, 0.0, 0.0, 0.0, 0.0
    };

    const int steps = static_cast<int>(
        std::llround(FINAL_TIME / DT)
    );
    std::vector<Sample> samples;
    samples.reserve(
        static_cast<std::size_t>(steps / SAMPLE_STRIDE + 1)
    );

    for (int step = 0; step <= steps; ++step) {
        const double t = static_cast<double>(step) * DT;

        if (step % SAMPLE_STRIDE == 0) {
            const auto [desired_velocity, desired_acceleration] =
                reference(t);
            const double tracking_error =
                state[0] - desired_velocity;
            const double control_input =
                state[3] * desired_acceleration
                + state[4] * state[0]
                - FEEDBACK_GAIN * tracking_error;

            samples.push_back(
                Sample{
                    t,
                    desired_velocity,
                    state[0],
                    tracking_error,
                    control_input,
                    state[3],
                    state[4],
                    minimumEigenvalueSymmetric2x2(
                        state[5], state[6], state[7]
                    )
                }
            );
        }

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

    return samples;
}

Metrics computeMetrics(const std::vector<Sample>& samples) {
    if (samples.empty()) {
        throw std::runtime_error("No simulation samples.");
    }

    double squared_error_sum = 0.0;
    for (const Sample& sample : samples) {
        squared_error_sum +=
            sample.tracking_error * sample.tracking_error;
    }

    const Sample& terminal = samples.back();
    return Metrics{
        std::sqrt(
            squared_error_sum
            / static_cast<double>(samples.size())
        ),
        std::abs(terminal.mass_hat - MASS_TRUE),
        std::abs(terminal.damping_hat - DAMPING_TRUE)
    };
}

void writeCsv(
    const std::vector<Sample>& direct,
    const std::vector<Sample>& composite,
    const std::string& file_name
) {
    if (direct.size() != composite.size()) {
        throw std::runtime_error(
            "Direct and composite data sizes differ."
        );
    }

    std::ofstream file(file_name);
    if (!file) {
        throw std::runtime_error(
            "Could not open output CSV file."
        );
    }

    file
        << "time,desired_velocity,direct_velocity,"
        << "composite_velocity,direct_error,composite_error,"
        << "direct_mass_hat,composite_mass_hat,"
        << "direct_damping_hat,composite_damping_hat,"
        << "information_min_eigenvalue\n";
    file << std::setprecision(12);

    for (std::size_t index = 0; index < direct.size(); ++index) {
        file
            << direct[index].time << ','
            << direct[index].desired_velocity << ','
            << direct[index].velocity << ','
            << composite[index].velocity << ','
            << direct[index].tracking_error << ','
            << composite[index].tracking_error << ','
            << direct[index].mass_hat << ','
            << composite[index].mass_hat << ','
            << direct[index].damping_hat << ','
            << composite[index].damping_hat << ','
            << composite[index].information_min_eigenvalue
            << '\n';
    }
}

void printMetrics(
    const std::string& label,
    const Metrics& metrics
) {
    std::cout << label << '\n';
    std::cout
        << "  tracking_rmse: "
        << metrics.tracking_rmse << '\n';
    std::cout
        << "  terminal_mass_error: "
        << metrics.terminal_mass_error << '\n';
    std::cout
        << "  terminal_damping_error: "
        << metrics.terminal_damping_error << '\n';
}

}  // namespace

int main() {
    try {
        const std::vector<Sample> direct = simulate(false);
        const std::vector<Sample> composite = simulate(true);

        std::cout << std::fixed << std::setprecision(6);
        printMetrics(
            "Tracking-error-only controller:",
            computeMetrics(direct)
        );
        printMetrics(
            "Composite controller:",
            computeMetrics(composite)
        );

        const std::string output_file =
            "Chapter21_Lesson5_results_cpp.csv";
        writeCsv(direct, composite, output_file);
        std::cout << "Saved " << output_file << '\n';
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }

    return 0;
}

Chapter21_Lesson5.java

The Java implementation uses records for structured samples and writes the same comparison variables to CSV.

/*
Chapter21_Lesson5.java
Composite adaptive control of an uncertain mass-damper velocity servo.

Build:
    javac Chapter21_Lesson5.java

Run:
    java Chapter21_Lesson5
*/

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter21_Lesson5 {
    private static final double MASS_TRUE = 2.5;
    private static final double DAMPING_TRUE = 1.2;
    private static final double FEEDBACK_GAIN = 6.0;
    private static final double GAMMA_MASS = 4.0;
    private static final double GAMMA_DAMPING = 4.0;
    private static final double COMPOSITE_GAIN = 2.0;
    private static final double FILTER_RATE = 8.0;
    private static final double DT = 0.001;
    private static final double FINAL_TIME = 20.0;
    private static final int SAMPLE_STRIDE = 10;

    private Chapter21_Lesson5() {
        // Utility class.
    }

    private record Reference(
        double velocity,
        double acceleration
    ) {}

    private record Sample(
        double time,
        double desiredVelocity,
        double velocity,
        double trackingError,
        double controlInput,
        double massHat,
        double dampingHat,
        double informationMinEigenvalue
    ) {}

    private record Metrics(
        double trackingRmse,
        double terminalMassError,
        double terminalDampingError
    ) {}

    private static Reference reference(double t) {
        double velocity =
            0.8 * Math.sin(0.7 * t)
            + 0.35 * Math.sin(1.9 * t);
        double acceleration =
            0.56 * Math.cos(0.7 * t)
            + 0.665 * Math.cos(1.9 * t);
        return new Reference(velocity, acceleration);
    }

    private static double minimumEigenvalueSymmetric2x2(
        double a,
        double b,
        double d
    ) {
        double trace = a + d;
        double discriminant = Math.sqrt(
            (a - d) * (a - d) + 4.0 * b * b
        );
        return 0.5 * (trace - discriminant);
    }

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

    private static double[] derivative(
        double t,
        double[] state,
        double compositeGain
    ) {
        double velocity = state[0];
        double velocityFiltered = state[1];
        double inputFiltered = state[2];
        double massHat = state[3];
        double dampingHat = state[4];
        double r11 = state[5];
        double r12 = state[6];
        double r22 = state[7];
        double q1 = state[8];
        double q2 = state[9];

        Reference desired = reference(t);
        double trackingError =
            velocity - desired.velocity();

        double controlInput =
            massHat * desired.acceleration()
            + dampingHat * velocity
            - FEEDBACK_GAIN * trackingError;

        double velocityDot =
            (
                controlInput
                - DAMPING_TRUE * velocity
            )
            / MASS_TRUE;
        double velocityFilteredDot =
            -FILTER_RATE * velocityFiltered
            + FILTER_RATE * velocity;
        double inputFilteredDot =
            -FILTER_RATE * inputFiltered
            + FILTER_RATE * controlInput;

        double filteredAcceleration =
            FILTER_RATE * (velocity - velocityFiltered);
        double trackingPhi1 = desired.acceleration();
        double trackingPhi2 = velocity;
        double predictionY1 = filteredAcceleration;
        double predictionY2 = velocityFiltered;

        double residual1 =
            r11 * massHat + r12 * dampingHat - q1;
        double residual2 =
            r12 * massHat + r22 * dampingHat - q2;
        double normalization = 1.0 + r11 + r22;

        double massHatDot =
            -GAMMA_MASS
            * (
                trackingPhi1 * trackingError
                + compositeGain
                    * residual1
                    / normalization
            );
        double dampingHatDot =
            -GAMMA_DAMPING
            * (
                trackingPhi2 * trackingError
                + compositeGain
                    * residual2
                    / normalization
            );

        double r11Dot = predictionY1 * predictionY1;
        double r12Dot = predictionY1 * predictionY2;
        double r22Dot = predictionY2 * predictionY2;
        double q1Dot = predictionY1 * inputFiltered;
        double q2Dot = predictionY2 * inputFiltered;

        return new double[] {
            velocityDot,
            velocityFilteredDot,
            inputFilteredDot,
            massHatDot,
            dampingHatDot,
            r11Dot,
            r12Dot,
            r22Dot,
            q1Dot,
            q2Dot
        };
    }

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

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

    private static List<Sample> simulate(boolean composite) {
        double compositeGain =
            composite ? COMPOSITE_GAIN : 0.0;
        double[] state = {
            0.0, 0.0, 0.0, 1.0, 0.2,
            0.0, 0.0, 0.0, 0.0, 0.0
        };

        int steps = (int) Math.round(FINAL_TIME / DT);
        List<Sample> samples = new ArrayList<>(
            steps / SAMPLE_STRIDE + 1
        );

        for (int step = 0; step <= steps; step++) {
            double t = step * DT;

            if (step % SAMPLE_STRIDE == 0) {
                Reference desired = reference(t);
                double trackingError =
                    state[0] - desired.velocity();
                double controlInput =
                    state[3] * desired.acceleration()
                    + state[4] * state[0]
                    - FEEDBACK_GAIN * trackingError;

                samples.add(
                    new Sample(
                        t,
                        desired.velocity(),
                        state[0],
                        trackingError,
                        controlInput,
                        state[3],
                        state[4],
                        minimumEigenvalueSymmetric2x2(
                            state[5],
                            state[6],
                            state[7]
                        )
                    )
                );
            }

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

        return samples;
    }

    private static Metrics computeMetrics(
        List<Sample> samples
    ) {
        if (samples.isEmpty()) {
            throw new IllegalArgumentException(
                "No simulation samples."
            );
        }

        double squaredErrorSum = 0.0;
        for (Sample sample : samples) {
            squaredErrorSum +=
                sample.trackingError()
                * sample.trackingError();
        }

        Sample terminal = samples.get(samples.size() - 1);
        return new Metrics(
            Math.sqrt(
                squaredErrorSum / samples.size()
            ),
            Math.abs(
                terminal.massHat() - MASS_TRUE
            ),
            Math.abs(
                terminal.dampingHat() - DAMPING_TRUE
            )
        );
    }

    private static void writeCsv(
        List<Sample> direct,
        List<Sample> composite,
        Path destination
    ) throws IOException {
        if (direct.size() != composite.size()) {
            throw new IllegalArgumentException(
                "Direct and composite data sizes differ."
            );
        }

        try (
            BufferedWriter writer = Files.newBufferedWriter(
                destination,
                StandardCharsets.UTF_8
            )
        ) {
            writer.write(
                "time,desired_velocity,direct_velocity,"
                + "composite_velocity,direct_error,"
                + "composite_error,direct_mass_hat,"
                + "composite_mass_hat,direct_damping_hat,"
                + "composite_damping_hat,"
                + "information_min_eigenvalue"
            );
            writer.newLine();

            for (int index = 0; index < direct.size(); index++) {
                Sample d = direct.get(index);
                Sample c = composite.get(index);
                writer.write(
                    String.format(
                        Locale.US,
                        "%.12g,%.12g,%.12g,%.12g,"
                            + "%.12g,%.12g,%.12g,%.12g,"
                            + "%.12g,%.12g,%.12g",
                        d.time(),
                        d.desiredVelocity(),
                        d.velocity(),
                        c.velocity(),
                        d.trackingError(),
                        c.trackingError(),
                        d.massHat(),
                        c.massHat(),
                        d.dampingHat(),
                        c.dampingHat(),
                        c.informationMinEigenvalue()
                    )
                );
                writer.newLine();
            }
        }
    }

    private static void printMetrics(
        String label,
        Metrics metrics
    ) {
        System.out.println(label);
        System.out.printf(
            Locale.US,
            "  tracking_rmse: %.6f%n",
            metrics.trackingRmse()
        );
        System.out.printf(
            Locale.US,
            "  terminal_mass_error: %.6f%n",
            metrics.terminalMassError()
        );
        System.out.printf(
            Locale.US,
            "  terminal_damping_error: %.6f%n",
            metrics.terminalDampingError()
        );
    }

    public static void main(String[] args) {
        try {
            List<Sample> direct = simulate(false);
            List<Sample> composite = simulate(true);

            printMetrics(
                "Tracking-error-only controller:",
                computeMetrics(direct)
            );
            printMetrics(
                "Composite controller:",
                computeMetrics(composite)
            );

            Path output = Path.of(
                "Chapter21_Lesson5_results_java.csv"
            );
            writeCsv(direct, composite, output);
            System.out.println(
                "Saved " + output
            );
        } catch (IOException | RuntimeException error) {
            System.err.println(
                "Error: " + error.getMessage()
            );
            System.exit(1);
        }
    }
}

Chapter21_Lesson5.m

The MATLAB script simulates both controllers, exports results and a figure, and contains an optional function that constructs a Simulink model. Set buildSimulinkModel = true when Simulink is installed.

%% Chapter21_Lesson5.m
% Composite adaptive control of an uncertain mass-damper velocity servo.
%
% This script compares:
%   1. Tracking-error-only adaptive control.
%   2. Composite adaptive control using tracking error and a filtered-data
%      prediction residual.
%
% Required MATLAB products:
%   - MATLAB
% Optional:
%   - Simulink, for automatic construction of a continuous-time model.
%
% Run:
%   Chapter21_Lesson5

clear;
close all;
clc;

configuration = struct( ...
    'massTrue', 2.5, ...
    'dampingTrue', 1.2, ...
    'feedbackGain', 6.0, ...
    'gammaMass', 4.0, ...
    'gammaDamping', 4.0, ...
    'compositeGain', 2.0, ...
    'filterRate', 8.0, ...
    'dt', 0.001, ...
    'finalTime', 20.0, ...
    'sampleStride', 10);

direct = simulateController(false, configuration);
composite = simulateController(true, configuration);

directMetrics = calculateMetrics(direct, configuration);
compositeMetrics = calculateMetrics(composite, configuration);

fprintf('Tracking-error-only controller:\n');
printMetrics(directMetrics);
fprintf('Composite controller:\n');
printMetrics(compositeMetrics);

results = table( ...
    direct.time, ...
    direct.desiredVelocity, ...
    direct.velocity, ...
    composite.velocity, ...
    direct.trackingError, ...
    composite.trackingError, ...
    direct.massHat, ...
    composite.massHat, ...
    direct.dampingHat, ...
    composite.dampingHat, ...
    composite.informationMinEigenvalue, ...
    'VariableNames', { ...
        'time', ...
        'desired_velocity', ...
        'direct_velocity', ...
        'composite_velocity', ...
        'direct_error', ...
        'composite_error', ...
        'direct_mass_hat', ...
        'composite_mass_hat', ...
        'direct_damping_hat', ...
        'composite_damping_hat', ...
        'information_min_eigenvalue'});

writetable(results, 'Chapter21_Lesson5_results_matlab.csv');

figure('Name', 'Chapter 21 Lesson 5');

subplot(2, 2, 1);
plot(direct.time, direct.desiredVelocity, 'LineWidth', 1.2);
hold on;
plot(direct.time, direct.velocity, 'LineWidth', 1.0);
plot(composite.time, composite.velocity, 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Velocity [m/s]');
title('Velocity tracking');
legend('Reference', 'Tracking-only', 'Composite', ...
    'Location', 'best');

subplot(2, 2, 2);
plot(direct.time, direct.trackingError, 'LineWidth', 1.0);
hold on;
plot(composite.time, composite.trackingError, 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Error [m/s]');
title('Tracking error');
legend('Tracking-only', 'Composite', 'Location', 'best');

subplot(2, 2, 3);
yline(configuration.massTrue, '--', 'True mass');
hold on;
plot(direct.time, direct.massHat, 'LineWidth', 1.0);
plot(composite.time, composite.massHat, 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Mass [kg]');
title('Mass estimate');
legend('True', 'Tracking-only', 'Composite', ...
    'Location', 'best');

subplot(2, 2, 4);
yline(configuration.dampingTrue, '--', 'True damping');
hold on;
plot(direct.time, direct.dampingHat, 'LineWidth', 1.0);
plot(composite.time, composite.dampingHat, 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Damping [N s/m]');
title('Damping estimate');
legend('True', 'Tracking-only', 'Composite', ...
    'Location', 'best');

exportgraphics(gcf, 'Chapter21_Lesson5_matlab.png', ...
    'Resolution', 180);

fprintf('Saved Chapter21_Lesson5_results_matlab.csv\n');
fprintf('Saved Chapter21_Lesson5_matlab.png\n');

% Set this flag to true to create Chapter21_Lesson5_Simulink.slx.
buildSimulinkModel = false;
if buildSimulinkModel
    buildCompositeAdaptiveSimulinkModel(configuration);
end

%% Local functions

function data = simulateController(composite, configuration)
    if composite
        compositeGain = configuration.compositeGain;
    else
        compositeGain = 0.0;
    end

    % State order:
    % [v, v_f, u_f, m_hat, b_hat, R11, R12, R22, q1, q2]
    state = [ ...
        0.0; 0.0; 0.0; 1.0; 0.2; ...
        0.0; 0.0; 0.0; 0.0; 0.0];

    steps = round(configuration.finalTime / configuration.dt);
    sampleCount = floor(steps / configuration.sampleStride) + 1;
    records = zeros(sampleCount, 8);
    recordIndex = 1;

    for step = 0:steps
        time = step * configuration.dt;

        if mod(step, configuration.sampleStride) == 0
            [desiredVelocity, desiredAcceleration] = ...
                referenceSignal(time);
            trackingError = state(1) - desiredVelocity;
            controlInput = ...
                state(4) * desiredAcceleration ...
                + state(5) * state(1) ...
                - configuration.feedbackGain * trackingError;

            informationMatrix = [ ...
                state(6), state(7); ...
                state(7), state(8)];
            eigenvalues = eig(informationMatrix);

            records(recordIndex, :) = [ ...
                time, ...
                desiredVelocity, ...
                state(1), ...
                trackingError, ...
                controlInput, ...
                state(4), ...
                state(5), ...
                min(eigenvalues)];
            recordIndex = recordIndex + 1;
        end

        if step < steps
            state = rk4Step( ...
                time, ...
                state, ...
                compositeGain, ...
                configuration);
        end
    end

    data = struct( ...
        'time', records(:, 1), ...
        'desiredVelocity', records(:, 2), ...
        'velocity', records(:, 3), ...
        'trackingError', records(:, 4), ...
        'controlInput', records(:, 5), ...
        'massHat', records(:, 6), ...
        'dampingHat', records(:, 7), ...
        'informationMinEigenvalue', records(:, 8));
end

function nextState = rk4Step( ...
    time, state, compositeGain, configuration)

    dt = configuration.dt;
    k1 = stateDerivative( ...
        time, state, compositeGain, configuration);
    k2 = stateDerivative( ...
        time + 0.5 * dt, ...
        state + 0.5 * dt * k1, ...
        compositeGain, ...
        configuration);
    k3 = stateDerivative( ...
        time + 0.5 * dt, ...
        state + 0.5 * dt * k2, ...
        compositeGain, ...
        configuration);
    k4 = stateDerivative( ...
        time + dt, ...
        state + dt * k3, ...
        compositeGain, ...
        configuration);

    nextState = state ...
        + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
end

function derivative = stateDerivative( ...
    time, state, compositeGain, configuration)

    velocity = state(1);
    velocityFiltered = state(2);
    inputFiltered = state(3);
    massHat = state(4);
    dampingHat = state(5);
    r11 = state(6);
    r12 = state(7);
    r22 = state(8);
    q1 = state(9);
    q2 = state(10);

    [desiredVelocity, desiredAcceleration] = ...
        referenceSignal(time);
    trackingError = velocity - desiredVelocity;

    controlInput = ...
        massHat * desiredAcceleration ...
        + dampingHat * velocity ...
        - configuration.feedbackGain * trackingError;

    velocityDot = ...
        (controlInput ...
        - configuration.dampingTrue * velocity) ...
        / configuration.massTrue;
    velocityFilteredDot = ...
        -configuration.filterRate * velocityFiltered ...
        + configuration.filterRate * velocity;
    inputFilteredDot = ...
        -configuration.filterRate * inputFiltered ...
        + configuration.filterRate * controlInput;

    filteredAcceleration = ...
        configuration.filterRate ...
        * (velocity - velocityFiltered);
    trackingRegressor = [desiredAcceleration; velocity];
    predictionRegressor = [ ...
        filteredAcceleration; velocityFiltered];

    informationMatrix = [r11, r12; r12, r22];
    informationVector = [q1; q2];
    parameterEstimate = [massHat; dampingHat];

    predictionResidual = ...
        informationMatrix * parameterEstimate ...
        - informationVector;
    normalization = 1.0 + trace(informationMatrix);
    adaptationGain = diag([ ...
        configuration.gammaMass, ...
        configuration.gammaDamping]);

    estimateDot = ...
        -adaptationGain ...
        * ( ...
            trackingRegressor * trackingError ...
            + compositeGain ...
                * predictionResidual ...
                / normalization);

    informationMatrixDot = ...
        predictionRegressor * predictionRegressor.';
    informationVectorDot = ...
        predictionRegressor * inputFiltered;

    derivative = [ ...
        velocityDot; ...
        velocityFilteredDot; ...
        inputFilteredDot; ...
        estimateDot(1); ...
        estimateDot(2); ...
        informationMatrixDot(1, 1); ...
        informationMatrixDot(1, 2); ...
        informationMatrixDot(2, 2); ...
        informationVectorDot(1); ...
        informationVectorDot(2)];
end

function [velocity, acceleration] = referenceSignal(time)
    velocity = ...
        0.8 * sin(0.7 * time) ...
        + 0.35 * sin(1.9 * time);
    acceleration = ...
        0.56 * cos(0.7 * time) ...
        + 0.665 * cos(1.9 * time);
end

function values = calculateMetrics(data, configuration)
    values = struct( ...
        'trackingRMSE', ...
            sqrt(mean(data.trackingError .^ 2)), ...
        'terminalMassError', ...
            abs(data.massHat(end) - configuration.massTrue), ...
        'terminalDampingError', ...
            abs( ...
                data.dampingHat(end) ...
                - configuration.dampingTrue));
end

function printMetrics(values)
    fprintf('  tracking_rmse: %.6f\n', values.trackingRMSE);
    fprintf( ...
        '  terminal_mass_error: %.6f\n', ...
        values.terminalMassError);
    fprintf( ...
        '  terminal_damping_error: %.6f\n', ...
        values.terminalDampingError);
end

function buildCompositeAdaptiveSimulinkModel(configuration)
    % This function builds a compact Simulink realization whose single
    % vector Integrator stores all ten continuous states. A MATLAB Function
    % block computes the plant, filters, information matrix, and composite
    % adaptive law.

    modelName = 'Chapter21_Lesson5_Simulink';

    if ~license('test', 'Simulink')
        error('A Simulink license is required.');
    end

    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end

    modelFile = [modelName, '.slx'];
    if isfile(modelFile)
        delete(modelFile);
    end

    new_system(modelName);
    open_system(modelName);

    add_block( ...
        'simulink/Sources/Clock', ...
        [modelName, '/Clock'], ...
        'Position', [40, 90, 70, 110]);

    add_block( ...
        'simulink/Continuous/Integrator', ...
        [modelName, '/State Integrator'], ...
        'InitialCondition', ...
            '[0;0;0;1;0.2;0;0;0;0;0]', ...
        'Position', [420, 80, 455, 120]);

    add_block( ...
        'simulink/User-Defined Functions/MATLAB Function', ...
        [modelName, '/Composite Dynamics'], ...
        'Position', [170, 55, 340, 145]);

    add_block( ...
        'simulink/Sinks/To Workspace', ...
        [modelName, '/Logged Signals'], ...
        'VariableName', 'simulinkCompositeLog', ...
        'SaveFormat', 'Array', ...
        'Position', [420, 175, 515, 205]);

    chart = find( ...
        sfroot, ...
        '-isa', 'Stateflow.EMChart', ...
        'Path', [modelName, '/Composite Dynamics']);
    if isempty(chart)
        error('Could not access the MATLAB Function block.');
    end

    chart.Script = sprintf([ ...
        'function [dx, logData] = fcn(t, x)\n' ...
        '%%#codegen\n' ...
        'massTrue = %.17g;\n' ...
        'dampingTrue = %.17g;\n' ...
        'feedbackGain = %.17g;\n' ...
        'gammaMass = %.17g;\n' ...
        'gammaDamping = %.17g;\n' ...
        'compositeGain = %.17g;\n' ...
        'filterRate = %.17g;\n' ...
        'v = x(1); vf = x(2); uf = x(3);\n' ...
        'mHat = x(4); bHat = x(5);\n' ...
        'R11 = x(6); R12 = x(7); R22 = x(8);\n' ...
        'q1 = x(9); q2 = x(10);\n' ...
        'vd = 0.8*sin(0.7*t) + 0.35*sin(1.9*t);\n' ...
        'dvd = 0.56*cos(0.7*t) + 0.665*cos(1.9*t);\n' ...
        'e = v - vd;\n' ...
        'u = mHat*dvd + bHat*v - feedbackGain*e;\n' ...
        'vDot = (u - dampingTrue*v)/massTrue;\n' ...
        'vfDot = -filterRate*vf + filterRate*v;\n' ...
        'ufDot = -filterRate*uf + filterRate*u;\n' ...
        'p = filterRate*(v-vf);\n' ...
        'res1 = R11*mHat + R12*bHat - q1;\n' ...
        'res2 = R12*mHat + R22*bHat - q2;\n' ...
        'normalization = 1 + R11 + R22;\n' ...
        'mHatDot = -gammaMass*(dvd*e + compositeGain*res1/normalization);\n' ...
        'bHatDot = -gammaDamping*(v*e + compositeGain*res2/normalization);\n' ...
        'dx = [vDot;vfDot;ufDot;mHatDot;bHatDot;p*p;p*vf;vf*vf;p*uf;vf*uf];\n' ...
        'logData = [t;vd;v;e;u;mHat;bHat];\n' ...
        'end\n'], ...
        configuration.massTrue, ...
        configuration.dampingTrue, ...
        configuration.feedbackGain, ...
        configuration.gammaMass, ...
        configuration.gammaDamping, ...
        configuration.compositeGain, ...
        configuration.filterRate);

    set_param(modelName, ...
        'StopTime', num2str(configuration.finalTime), ...
        'Solver', 'ode4', ...
        'FixedStep', num2str(configuration.dt));

    set_param(modelName, 'SimulationCommand', 'update');

    add_line( ...
        modelName, ...
        'Clock/1', ...
        'Composite Dynamics/1', ...
        'autorouting', 'on');
    add_line( ...
        modelName, ...
        'State Integrator/1', ...
        'Composite Dynamics/2', ...
        'autorouting', 'on');
    add_line( ...
        modelName, ...
        'Composite Dynamics/1', ...
        'State Integrator/1', ...
        'autorouting', 'on');
    add_line( ...
        modelName, ...
        'Composite Dynamics/2', ...
        'Logged Signals/1', ...
        'autorouting', 'on');

    save_system(modelName, modelFile);
    fprintf('Created %s\n', modelFile);
end

Chapter21_Lesson5.nb

The notebook expression contains a complete input cell based on NDSolveValue, exports comparison data, and creates four plots.

(* Chapter21_Lesson5.nb *)
Notebook[{
  Cell["Chapter 21, Lesson 5: Composite Adaptive Control", "Title"], Cell["Comparison of tracking-error-only and composite adaptation for an uncertain mass-damper velocity servo.", "Text"],
  Cell[BoxData["(* Composite adaptive control simulation for Chapter 21, Lesson 5. *)
ClearAll[\"Global`*\"];
massTrue = 2.5; dampingTrue = 1.2; feedbackGain = 6.0; gammaMass = 4.0; gammaDamping = 4.0; compositeGain = 2.0; filterRate = 8.0; finalTime = 20.0; samplePeriod = 0.01;
desiredVelocity[t_] := 0.8 Sin[0.7 t] + 0.35 Sin[1.9 t]; desiredAcceleration[t_] := 0.56 Cos[0.7 t] + 0.665 Cos[1.9 t];
simulate[gain_] := Module[
  {v, vf, uf, mHat, bHat, r11, r12, r22, q1, q2, error, control, filteredAcceleration, residual1, residual2, normalization, solution},
  error[t_] := v[t] - desiredVelocity[t];
  control[t_] := mHat[t] desiredAcceleration[t] + bHat[t] v[t] - feedbackGain error[t];
  filteredAcceleration[t_] := filterRate (v[t] - vf[t]);
  residual1[t_] := r11[t] mHat[t] + r12[t] bHat[t] - q1[t]; residual2[t_] := r12[t] mHat[t] + r22[t] bHat[t] - q2[t];
  normalization[t_] := 1 + r11[t] + r22[t];
  solution = NDSolveValue[
    {
      v'[t] == (control[t] - dampingTrue v[t])/massTrue, vf'[t] == -filterRate vf[t] + filterRate v[t],
      uf'[t] == -filterRate uf[t] + filterRate control[t],
      mHat'[t] == -gammaMass (desiredAcceleration[t] error[t] + gain residual1[t]/normalization[t]),
      bHat'[t] == -gammaDamping (v[t] error[t] + gain residual2[t]/normalization[t]),
      r11'[t] == filteredAcceleration[t]^2, r12'[t] == filteredAcceleration[t] vf[t], r22'[t] == vf[t]^2,
      q1'[t] == filteredAcceleration[t] uf[t], q2'[t] == vf[t] uf[t], v[0] == 0, vf[0] == 0, uf[0] == 0,
      mHat[0] == 1.0, bHat[0] == 0.2, r11[0] == 0, r12[0] == 0, r22[0] == 0, q1[0] == 0, q2[0] == 0
    },
    {v, vf, uf, mHat, bHat, r11, r12, r22, q1, q2}, {t, 0, finalTime},
    Method -> {\"EquationSimplification\" -> \"Residual\"}
  ];
  <|\"Velocity\" -> solution[[1]], \"VelocityFiltered\" -> solution[[2]], \"InputFiltered\" -> solution[[3]], \"MassHat\" -> solution[[4]], \"DampingHat\" -> solution[[5]],
    \"R11\" -> solution[[6]], \"R12\" -> solution[[7]], \"R22\" -> solution[[8]], \"Q1\" -> solution[[9]], \"Q2\" -> solution[[10]]|>
];
direct = simulate[0.0]; composite = simulate[compositeGain]; timeGrid = Range[0, finalTime, samplePeriod];
trackingError[data_, time_] := data[\"Velocity\"][time] - desiredVelocity[time]; trackingRMSE[data_] := Sqrt[Mean[(trackingError[data, #]^2) & /@ timeGrid]];
terminalMassError[data_] := Abs[data[\"MassHat\"][finalTime] - massTrue]; terminalDampingError[data_] := Abs[data[\"DampingHat\"][finalTime] - dampingTrue];
Print[\"Tracking-error-only controller:\"]; Print[\"  tracking_rmse: \", NumberForm[trackingRMSE[direct], {8, 6}]];
Print[\"  terminal_mass_error: \", NumberForm[terminalMassError[direct], {8, 6}]]; Print[\"  terminal_damping_error: \", NumberForm[terminalDampingError[direct], {8, 6}]];
Print[\"Composite controller:\"]; Print[\"  tracking_rmse: \", NumberForm[trackingRMSE[composite], {8, 6}]];
Print[\"  terminal_mass_error: \", NumberForm[terminalMassError[composite], {8, 6}]]; Print[\"  terminal_damping_error: \", NumberForm[terminalDampingError[composite], {8, 6}]];
minimumInformationEigenvalue[data_, time_] := Min[Eigenvalues[{{data[\"R11\"][time], data[\"R12\"][time]}, {data[\"R12\"][time], data[\"R22\"][time]}}]];
csvData = Prepend[Table[{time, desiredVelocity[time], direct[\"Velocity\"][time], composite[\"Velocity\"][time],
  trackingError[direct, time], trackingError[composite, time], direct[\"MassHat\"][time],
  composite[\"MassHat\"][time], direct[\"DampingHat\"][time], composite[\"DampingHat\"][time],
  minimumInformationEigenvalue[composite, time]}, {time, timeGrid}],
  {\"time\", \"desired_velocity\", \"direct_velocity\", \"composite_velocity\", \"direct_error\", \"composite_error\", \"direct_mass_hat\", \"composite_mass_hat\",
   \"direct_damping_hat\", \"composite_damping_hat\", \"information_min_eigenvalue\"}];
Export[\"Chapter21_Lesson5_results_mathematica.csv\", csvData];
trackingPlot = Plot[{desiredVelocity[t], direct[\"Velocity\"][t], composite[\"Velocity\"][t]}, {t, 0, finalTime},
  PlotLegends -> {\"Reference\", \"Tracking-only\", \"Composite\"}, Frame -> True, FrameLabel -> {\"Time [s]\", \"Velocity [m/s]\"}, PlotLabel -> \"Velocity Tracking\", ImageSize -> Large];
errorPlot = Plot[{trackingError[direct, t], trackingError[composite, t]}, {t, 0, finalTime},
  PlotLegends -> {\"Tracking-only\", \"Composite\"}, Frame -> True, FrameLabel -> {\"Time [s]\", \"Error [m/s]\"}, PlotLabel -> \"Tracking Error\", ImageSize -> Large];
massPlot = Plot[{massTrue, direct[\"MassHat\"][t], composite[\"MassHat\"][t]}, {t, 0, finalTime},
  PlotLegends -> {\"True mass\", \"Tracking-only estimate\", \"Composite estimate\"}, Frame -> True, FrameLabel -> {\"Time [s]\", \"Mass [kg]\"}, PlotLabel -> \"Mass Estimate\", ImageSize -> Large];
dampingPlot = Plot[{dampingTrue, direct[\"DampingHat\"][t], composite[\"DampingHat\"][t]}, {t, 0, finalTime},
  PlotLegends -> {\"True damping\", \"Tracking-only estimate\", \"Composite estimate\"}, Frame -> True, FrameLabel -> {\"Time [s]\", \"Damping [N s/m]\"}, PlotLabel -> \"Damping Estimate\", ImageSize -> Large];
comparisonFigure = GraphicsGrid[{{trackingPlot, errorPlot}, {massPlot, dampingPlot}}, ImageSize -> 1200];
Export[\"Chapter21_Lesson5_mathematica.png\", comparisonFigure, ImageResolution -> 180]; comparisonFigure
"], "Input"]
},
WindowSize -> {1280, 800}, WindowMargins -> {{Automatic, 0}, {Automatic, 0}}, StyleDefinitions -> "Default.nb"]                

17.1 Adaptation Gains

Larger entries of \( \boldsymbol{\Gamma} \) accelerate parameter motion but also increase sensitivity to noise, discretization, and unmodeled dynamics. The composite gain \( k_c \) controls how strongly stored prediction information influences adaptation.

17.2 Filter Bandwidth

A large \( \lambda \) makes \( p_f \) follow acceleration more closely, but it also transmits more high-frequency measurement noise. A small \( \lambda \) reduces noise but increases phase lag and slows the information channel.

17.3 Information-Matrix Conditioning

Monitor the minimum eigenvalue or condition number of \( \mathbf{R} \). A large trace with a very small minimum eigenvalue means that much data have been collected, but the data remain nearly collinear and do not identify all parameters.

17.4 Noise and Model Mismatch

In real hardware, the identity \( \mathbf{q}=\mathbf{R}\boldsymbol{\theta} \) is only approximate because of sensor noise, Coulomb friction, actuator dynamics, and numerical errors. Projection, leakage, dead zones, bounded-memory updates, or covariance weighting may be added using the robust modifications introduced earlier in the course.

17.5 Positivity of the Mass Estimate

The nominal proof does not require dividing by \( \hat m \), so temporary estimate error does not create an algebraic singularity in this example. Nevertheless, projection onto a physically meaningful set such as \( \hat m\ge m_{\min} >0 \) is advisable in safety-critical software.

17.6 Sampling

The filter and adaptation dynamics may be substantially faster than the mechanical plant. The sample period must resolve the fastest selected pole and update gain. Reducing the step size should produce convergent numerical results; if it does not, the implementation is numerically unreliable.

Problem 1 (Tracking Error Dynamics): Starting from \( m\dot v+bv=u \) and \( u=\hat m\dot v_d+\hat bv-ke \), derive the closed-loop tracking error equation.

Solution:

\[ \begin{aligned} m\dot e &= m(\dot v-\dot v_d)\\ &= u-bv-m\dot v_d\\ &= \hat m\dot v_d + \hat bv - ke - bv - m\dot v_d\\ &= -ke + \tilde m\dot v_d + \tilde bv\\ &= -ke + \boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}}. \end{aligned} \]

Problem 2 (Filtered Derivative Identity): Show that \( p_f=\lambda(v-v_f) \), with \( \dot v_f=-\lambda v_f+\lambda v \), represents a stable filtered version of \( \dot v \).

Solution:

Taking Laplace transforms with consistent initial conditions gives

\[ V_f(s) = \frac{\lambda}{s+\lambda}V(s). \]

Therefore,

\[ P_f(s) = \lambda\left[V(s)-V_f(s)\right] = \lambda \left( 1-\frac{\lambda}{s+\lambda} \right)V(s) = \frac{\lambda s}{s+\lambda}V(s) = H(s)\,sV(s). \]

Hence \( p_f=H(s)[\dot v] \). An unmatched initial condition introduces a decaying term proportional to \( e^{-\lambda t} \).

Problem 3 (Measurable Parameter-Error Signal): Prove that \( \boldsymbol{\zeta} = \bar{\mathbf{R}} \hat{\boldsymbol{\theta}} - \bar{\mathbf{q}} \) equals \( \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}} \).

Solution:

\[ \mathbf{q} = \int_0^t \mathbf{Y}_f^Tu_f\,d\tau = \int_0^t \mathbf{Y}_f^T \mathbf{Y}_f \boldsymbol{\theta}\,d\tau = \mathbf{R}\boldsymbol{\theta}. \]

Dividing both sides by \( \rho_0+\operatorname{tr}\mathbf{R} \) gives \( \bar{\mathbf{q}} = \bar{\mathbf{R}}\boldsymbol{\theta} \). Consequently,

\[ \boldsymbol{\zeta} = \bar{\mathbf{R}} \left( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta} \right) = \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}}. \]

Problem 4 (Composite Lyapunov Derivative): Derive \( \dot V \) under the composite law and explain the cancellation mechanism.

Solution:

The tracking dynamics create the mixed term \( e\boldsymbol{\phi}^T \tilde{\boldsymbol{\theta}} \). The tracking component of the adaptive law creates its negative:

\[ \tilde{\boldsymbol{\theta}}^T \boldsymbol{\Gamma}^{-1} \left( -\boldsymbol{\Gamma} \boldsymbol{\phi}e \right) = - \tilde{\boldsymbol{\theta}}^T \boldsymbol{\phi}e. \]

The two scalar mixed terms are equal and opposite. The remaining prediction component is nonpositive:

\[ \dot V = -ke^2 - k_c \tilde{\boldsymbol{\theta}}^T \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}} \le0. \]

Problem 5 (Loss of Identifiability): Suppose the desired velocity is a constant and the system reaches exact steady state. What happens to the filtered regressor rank?

Solution:

At steady state, \( \dot v_d=0 \), \( \dot v=0 \), and \( p_f=0 \). The filtered regressor becomes \( \mathbf{Y}_f=[0,v_d] \). Therefore,

\[ \mathbf{Y}_f^T\mathbf{Y}_f = \begin{bmatrix} 0 & 0\\ 0 & v_d^2 \end{bmatrix}, \]

which has rank one when \( v_d\ne0 \). The damping can be informed by steady-state data, but the mass cannot be identified because no acceleration is present. Tracking may still be excellent.

Problem 6 (Initial Control Calculation): For the numerical values in Section 11, compute the initial control input.

Solution:

At \( t=0 \), \( v_d(0)=0 \), \( \dot v_d(0)=0.56+0.665=1.225 \), \( v(0)=0 \), and \( e(0)=0 \). Hence

\[ u(0) = \hat m(0)\dot v_d(0) + \hat b(0)v(0) - ke(0) = 1.0(1.225) = 1.225\;\text{N}. \]

The mass-damper example demonstrated the complete composite adaptive control construction. A certainty-equivalent controller produced an error model linear in parameter error. A conventional tracking-error update guaranteed asymptotic tracking but did not directly penalize parameter error. Stable filtering created a measurable linear regression without acceleration measurement. Integrating that regression produced an information matrix and a computable prediction residual. Adding the residual to the adaptive law yielded

\[ \dot V = -ke^2 - k_c \tilde{\boldsymbol{\theta}}^T \bar{\mathbf{R}} \tilde{\boldsymbol{\theta}} \le0. \]

Thus the tracking channel and prediction channel contribute distinct negative terms. Tracking converges under bounded closed-loop signals, while parameter convergence requires sufficiently informative filtered data. This distinction between control success and identification success is central to the correct interpretation of adaptive systems.