In the ideal MRAC derivations studied earlier, the plant belongs exactly to the assumed parameterized model class. The update law cancels the state–parameter cross term in the Lyapunov derivative, and the tracking error converges to zero. Chapter 22, Lesson 1 showed why this conclusion can fail when the true plant contains disturbances or unmodeled dynamics. A small persistent error may continue to drive the parameter integrator even when further adaptation no longer improves tracking. The result is parameter drift.

Robust adaptive modifications do not attempt to identify every omitted phenomenon. Instead, they alter the parameter update mechanism so that bounded nonideal signals cannot produce unbounded estimates. The three mechanisms emphasized here are:

• Leakage: continuously pulls estimates toward a nominal vector.
• σ-modification: a constant isotropic leakage term used in adaptive control.
• Projection: prevents estimates from leaving a known admissible set.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> C["Adaptive controller"]
  P["Plant with bounded nonideal effects"] --> Y["Measured output"]
  C --> P
  RM --> E["Tracking error"]
  Y --> E
  E --> G["Gradient update direction"]
  G --> L["Add leakage toward nominal parameters"]
  L --> PR["Project onto admissible parameter set"]
  PR --> C
        

Consider the uncertain scalar plant

\[ \dot{x}=a x+b u+d(t), \qquad |d(t)|\leq \bar d, \qquad b\neq 0, \]

where the sign of \( b \) is known. The desired response is generated by the stable reference model

\[ \dot{x}_m=a_m x_m+b_m r, \qquad a_m<0. \]

Use the directly parameterized controller

\[ u=\hat\theta_x x+\hat\theta_r r =\hat{\boldsymbol\theta}^{T}\boldsymbol\phi, \qquad \boldsymbol\phi=\begin{bmatrix}x & r\end{bmatrix}^{T}. \]

If the ideal parameters exist, they satisfy the matching equations

\[ a+b\theta_x^{*}=a_m, \qquad b\theta_r^{*}=b_m, \]

\[ \theta_x^{*}=\frac{a_m-a}{b}, \qquad \theta_r^{*}=\frac{b_m}{b}. \]

Define tracking and parameter errors as \( e=x-x_m \) and \( \tilde{\boldsymbol\theta}=\hat{\boldsymbol\theta}- \boldsymbol\theta^{*} \). Their coupled dynamics are

\[ \dot e=a_m e+b\tilde{\boldsymbol\theta}^{T}\boldsymbol\phi+d(t). \]

Let \( Q>0 \) and choose \( P>0 \) from the scalar Lyapunov equation

\[ 2a_mP=-Q, \qquad P=-\frac{Q}{2a_m}. \]

The standard Lyapunov candidate is

\[ V=Pe^2+\tilde{\boldsymbol\theta}^{T}\Gamma^{-1} \tilde{\boldsymbol\theta}, \qquad \Gamma=\Gamma^{T}>0. \]

Differentiating gives

\[ \dot V=-Qe^2+2Peb\tilde{\boldsymbol\theta}^{T}\boldsymbol\phi +2\tilde{\boldsymbol\theta}^{T}\Gamma^{-1} \dot{\hat{\boldsymbol\theta} }+2Pe\,d(t). \]

The ideal gradient law

\[ \dot{\hat{\boldsymbol\theta} }_{\mathrm{grad} } =-\Gamma\boldsymbol\phi bPe \]

cancels the cross term, leaving

\[ \dot V=-Qe^2+2Pe\,d(t). \]

When \( d(t)=0 \), this is negative semidefinite. When a bounded disturbance persists, however, the second term has indefinite sign. Moreover, the update law still integrates \( \boldsymbol\phi e \). Thus even a small noise-induced error can accumulate over a long interval and move the estimates far from physically plausible values. Increasing \( \Gamma \) accelerates both desired adaptation and unwanted integration of nonideal error.

4.1 General leakage toward a nominal parameter vector

Introduce a nominal estimate \( \boldsymbol\theta_0 \) and write

\[ \dot{\hat{\boldsymbol\theta} } =-\Gamma\boldsymbol\phi bPe -\Gamma\Lambda(\hat{\boldsymbol\theta}-\boldsymbol\theta_0), \qquad \Lambda=\Lambda^{T}>0. \]

The first term learns from tracking error. The second term is a stable first-order forgetting mechanism. With no informative error signal, estimates exponentially approach \( \boldsymbol\theta_0 \) rather than freezing at an arbitrary value.

4.2 σ-modification

The classical σ-modification is obtained with \( \Lambda=\sigma I \), where \( \sigma>0 \):

\[ \dot{\hat{\boldsymbol\theta} } =-\Gamma\boldsymbol\phi bPe -\sigma\Gamma(\hat{\boldsymbol\theta}-\boldsymbol\theta_0). \]

Many texts set \( \boldsymbol\theta_0=\mathbf 0 \) and write \( -\sigma\Gamma\hat{\boldsymbol\theta} \). A nonzero nominal vector is often preferable when a credible nominal model is available. It reduces the steady bias introduced by leakage.

4.3 Relation to e-modification

For comparison, an e-modification scales leakage with tracking-error magnitude:

\[ \dot{\hat{\boldsymbol\theta} } =-\Gamma\boldsymbol\phi bPe -\sigma_e |e|\Gamma (\hat{\boldsymbol\theta}-\boldsymbol\theta_0). \]

Constant σ-modification supplies dissipation even for small errors; e-modification weakens the bias as tracking improves. This lesson uses constant leakage because it yields a direct uniform-ultimate-boundedness calculation.

Suppose prior physical knowledge provides a closed convex set \( \Omega \) containing the ideal parameter vector:

\[ \boldsymbol\theta^{*}\in\Omega. \]

Given an unconstrained update velocity \( \boldsymbol v \), projection replaces it by the closest velocity that lies in the tangent cone of \( \Omega \) at the current estimate:

\[ \operatorname{Proj}_{\Omega} (\hat{\boldsymbol\theta},\boldsymbol v) =\arg\min_{\boldsymbol z\in T_{\Omega}(\hat{\boldsymbol\theta})} (\boldsymbol z-\boldsymbol v)^{T}\Gamma^{-1} (\boldsymbol z-\boldsymbol v). \]

Its central Lyapunov property is

\[ \tilde{\boldsymbol\theta}^{T}\Gamma^{-1} \left[\operatorname{Proj}_{\Omega} (\hat{\boldsymbol\theta},\boldsymbol v)-\boldsymbol v\right] \leq 0. \]

Thus projection never makes the parameter-error contribution to the Lyapunov derivative worse, provided the ideal vector lies inside the admissible set.

5.1 Componentwise box projection

For practical software, take

\[ \Omega=\left\{\boldsymbol\theta: \underline\theta_i\leq\theta_i\leq\overline\theta_i, \;i=1,\ldots,p\right\}. \]

The projected update is

\[ \operatorname{Proj}_i(\hat\theta_i,v_i)= \begin{cases} 0, & \hat\theta_i=\underline\theta_i \text{ and } v_i<0,\\ 0, & \hat\theta_i=\overline\theta_i \text{ and } v_i>0,\\ v_i, & \text{otherwise}. \end{cases} \]

flowchart TD
  V["Candidate update velocity v"] --> I["Estimate strictly inside bounds?"]
  I -->|yes| A["Apply v unchanged"]
  I -->|no| O["Does v point outward at an active bound?"]
  O -->|yes| Z["Set outward component to zero"]
  O -->|no| A
  Z --> U["Integrate estimate and clip finite-step overshoot"]
  A --> U
        

Continuous-time projection makes \( \Omega \) forward invariant. In sampled software, a final numerical clip is also applied because a finite Euler or Runge–Kutta step can overshoot a boundary.

Define the unconstrained robust velocity

\[ \boldsymbol v= -\Gamma\boldsymbol\phi bPe -\sigma\Gamma(\hat{\boldsymbol\theta}-\boldsymbol\theta_0). \]

The complete update law is

\[ \boxed{\dot{\hat{\boldsymbol\theta} } =\operatorname{Proj}_{\Omega} (\hat{\boldsymbol\theta},\boldsymbol v)}. \]

Leakage and projection solve different problems. Leakage supplies dissipation throughout the interior and prevents slow random-walk behavior. Projection enforces hard physical limits but does not pull an interior estimate toward a preferred value. Their combination therefore gives both an ultimate-bound argument and explicit parameter safety.

Let \( \boldsymbol\delta_0=\boldsymbol\theta^{*}- \boldsymbol\theta_0 \). With the combined law and the projection inequality, the derivative of \( V=Pe^2+\tilde{\boldsymbol\theta}^{T} \Gamma^{-1}\tilde{\boldsymbol\theta} \) satisfies

\[ \dot V\leq -Qe^2+2Pe\,d(t) -2\sigma\tilde{\boldsymbol\theta}^{T} (\hat{\boldsymbol\theta}-\boldsymbol\theta_0). \]

Since

\[ \hat{\boldsymbol\theta}-\boldsymbol\theta_0 =\tilde{\boldsymbol\theta}+\boldsymbol\delta_0, \]

the leakage term becomes

\[ -2\sigma\|\tilde{\boldsymbol\theta}\|^2 -2\sigma\tilde{\boldsymbol\theta}^{T}\boldsymbol\delta_0. \]

Young’s inequality gives

\[ -2\sigma\tilde{\boldsymbol\theta}^{T}\boldsymbol\delta_0 \leq \sigma\|\tilde{\boldsymbol\theta}\|^2 +\sigma\|\boldsymbol\delta_0\|^2. \]

Similarly, because \( |d(t)|\leq\bar d \),

\[ 2P|e|\bar d \leq \frac{Q}{2}e^2+\frac{2P^2}{Q}\bar d^2. \]

Therefore,

\[ \dot V\leq -\frac{Q}{2}e^2 -\sigma\|\tilde{\boldsymbol\theta}\|^2 +\sigma\|\boldsymbol\delta_0\|^2 +\frac{2P^2}{Q}\bar d^2. \]

Define

\[ \alpha=\min\left\{\frac{Q}{2P}, \sigma\lambda_{\min}(\Gamma)\right\}, \qquad c_0=\sigma\|\boldsymbol\delta_0\|^2 +\frac{2P^2}{Q}\bar d^2. \]

Using

\[ \|\tilde{\boldsymbol\theta}\|^2 \geq\lambda_{\min}(\Gamma) \tilde{\boldsymbol\theta}^{T}\Gamma^{-1} \tilde{\boldsymbol\theta}, \]

we obtain the scalar comparison inequality

\[ \boxed{\dot V\leq-\alpha V+c_0}. \]

Hence

\[ V(t)\leq e^{-\alpha t}V(0) +\frac{c_0}{\alpha}\left(1-e^{-\alpha t}\right), \]

\[ \limsup_{t\to\infty}V(t)\leq\frac{c_0}{\alpha}, \qquad \limsup_{t\to\infty}|e(t)| \leq\sqrt{\frac{c_0}{\alpha P} }. \]

The conclusion is uniform ultimate boundedness, not exact parameter convergence. The ultimate set shrinks when the disturbance bound and nominal-parameter mismatch shrink. Increasing \( \sigma \) strengthens dissipation but can also increase bias through the term \( \sigma\|\boldsymbol\delta_0\|^2 \).

With \( \sigma=0 \), projection guarantees that \( \hat{\boldsymbol\theta}(t)\in\Omega \) and does not worsen the Lyapunov derivative. It does not, by itself, create a negative parameter-error term. One obtains

\[ \dot V\leq -Qe^2+2Pe\,d(t). \]

The estimate is bounded because it is constrained, and the tracking error can be shown to enter a disturbance-dependent neighborhood under additional boundedness arguments. However, the clean inequality \( \dot V\leq-\alpha V+c_0 \) generally requires a dissipative modification such as leakage. This is why projection and σ-modification are complementary rather than interchangeable.

1. Choose the admissible set from physics. Bounds should contain every credible ideal parameter with margin. Incorrect bounds can exclude the matching vector and create avoidable tracking bias.
2. Use a meaningful nominal vector. Setting \( \boldsymbol\theta_0 \) from a nominal plant model reduces the leakage mismatch \( \boldsymbol\delta_0 \).
3. Tune adaptation gain before leakage. Select \( \Gamma \) for acceptable transient speed without excessive noise amplification; then increase \( \sigma \) only enough to suppress drift.
4. Account for sampling. Use projected velocity plus a post-step clip. Log boundary activity because frequent clipping often indicates poor bounds, too much adaptation gain, or an inadequate plant model.
5. Do not interpret bounded estimates as identified physical parameters. Without persistent excitation, estimates may settle anywhere in a tracking-equivalent neighborhood.

Useful software libraries include NumPy/SciPy and Matplotlib in Python, Eigen or Armadillo in C++, Apache Commons Math in Java, MATLAB Control System Toolbox and Simulink, and NDSolve or list-based numerical iteration in Wolfram Language. The implementations below use minimal dependencies so the update logic remains explicit.

The common simulation uses

\[ a=-0.4,\quad b=1,\quad a_m=-1.5,\quad b_m=1.5, \quad Q=1,\quad P=\frac{1}{3}. \]

The ideal controller parameters are

\[ \boldsymbol\theta^{*}= \begin{bmatrix}-1.1 & 1.5\end{bmatrix}^{T}. \]

The robust settings are

\[ \Gamma=\operatorname{diag}(12,12),\quad \sigma=0.03,\quad \boldsymbol\theta_0= \begin{bmatrix}-0.8 & 1.2\end{bmatrix}^{T}, \]

\[ -2.5\leq\hat\theta_x\leq0.5,\qquad 0.2\leq\hat\theta_r\leq2.5. \]

The reference is sinusoidal for the first 20 seconds and then becomes zero. The loss of excitation, combined with deterministic measurement noise and a bounded disturbance, is deliberately chosen to reveal parameter drift in the baseline update. Four modes are compared: baseline, σ-modification, projection, and the combined law.

Chapter22_Lesson2.py

This version uses NumPy for array operations and Matplotlib for separate tracking-error and parameter-history figures. The integration is an explicit fixed-step Euler method so every adaptive update is visible.

"""Chapter22_Lesson2.py

Comparison of four adaptive laws for a first-order MRAC example:
1. Baseline gradient law
2. Sigma modification (constant leakage toward a nominal vector)
3. Box projection
4. Sigma modification plus box projection

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Literal

import matplotlib.pyplot as plt
import numpy as np

Mode = Literal["baseline", "sigma", "projection", "combined"]


@dataclass(frozen=True)
class Config:
    dt: float = 0.002
    final_time: float = 120.0
    a: float = -0.4
    b: float = 1.0
    a_m: float = -1.5
    b_m: float = 1.5
    q: float = 1.0
    gamma_x: float = 12.0
    gamma_r: float = 12.0
    sigma: float = 0.03
    theta0_x: float = -0.8
    theta0_r: float = 1.2
    lower_x: float = -2.5
    lower_r: float = 0.2
    upper_x: float = 0.5
    upper_r: float = 2.5

    @property
    def p(self) -> float:
        # Scalar Lyapunov equation: 2*a_m*p = -q.
        return -self.q / (2.0 * self.a_m)


def reference(t: float) -> float:
    """Exciting command first, then regulation with no persistent excitation."""
    return np.sin(0.8 * t) if t < 20.0 else 0.0


def deterministic_measurement_noise(t: float) -> float:
    return 0.06 * np.sin(37.0 * t) + 0.03 * np.sin(71.0 * t)


def disturbance(t: float) -> float:
    """Bounded disturbance representing bias and neglected dynamics."""
    return 0.25 * np.sin(12.0 * t) + 0.12 * np.sin(0.17 * t) + 0.12


def project_box(
    theta: np.ndarray,
    velocity: np.ndarray,
    lower: np.ndarray,
    upper: np.ndarray,
) -> np.ndarray:
    """Project an update velocity onto the tangent cone of a box."""
    projected = velocity.copy()
    for i in range(theta.size):
        if theta[i] <= lower[i] and projected[i] < 0.0:
            projected[i] = 0.0
        elif theta[i] >= upper[i] and projected[i] > 0.0:
            projected[i] = 0.0
    return projected


def simulate(mode: Mode, cfg: Config) -> dict[str, np.ndarray | float]:
    steps = int(round(cfg.final_time / cfg.dt)) + 1
    time = np.linspace(0.0, cfg.final_time, steps)

    x = np.zeros(steps)
    x_m = np.zeros(steps)
    u = np.zeros(steps)
    theta = np.zeros((steps, 2))

    gamma = np.array([cfg.gamma_x, cfg.gamma_r], dtype=float)
    theta0 = np.array([cfg.theta0_x, cfg.theta0_r], dtype=float)
    lower = np.array([cfg.lower_x, cfg.lower_r], dtype=float)
    upper = np.array([cfg.upper_x, cfg.upper_r], dtype=float)

    for k in range(steps - 1):
        t = time[k]
        r = reference(t)
        x_measured = x[k] + deterministic_measurement_noise(t)
        e_measured = x_measured - x_m[k]
        phi = np.array([x_measured, r], dtype=float)

        u[k] = float(theta[k] @ phi)

        x_dot = cfg.a * x[k] + cfg.b * u[k] + disturbance(t)
        x_m_dot = cfg.a_m * x_m[k] + cfg.b_m * r

        velocity = -gamma * phi * cfg.b * cfg.p * e_measured

        if mode in ("sigma", "combined"):
            velocity -= cfg.sigma * gamma * (theta[k] - theta0)

        if mode in ("projection", "combined"):
            velocity = project_box(theta[k], velocity, lower, upper)

        x[k + 1] = x[k] + cfg.dt * x_dot
        x_m[k + 1] = x_m[k] + cfg.dt * x_m_dot
        theta[k + 1] = theta[k] + cfg.dt * velocity

        # Clipping removes finite-step overshoot at a projection boundary.
        if mode in ("projection", "combined"):
            theta[k + 1] = np.clip(theta[k + 1], lower, upper)

    r_final = reference(time[-1])
    phi_final = np.array(
        [x[-1] + deterministic_measurement_noise(time[-1]), r_final]
    )
    u[-1] = float(theta[-1] @ phi_final)

    error = x - x_m
    rms_error = float(np.sqrt(np.mean(error**2)))
    max_error = float(np.max(np.abs(error)))

    return {
        "time": time,
        "x": x,
        "x_m": x_m,
        "u": u,
        "theta": theta,
        "error": error,
        "rms_error": rms_error,
        "max_error": max_error,
    }


def save_summary(results: dict[Mode, dict[str, np.ndarray | float]]) -> None:
    lines = ["mode,rms_error,max_abs_error,theta_x_final,theta_r_final"]
    for mode, data in results.items():
        theta = np.asarray(data["theta"])
        lines.append(
            f"{mode},{data['rms_error']:.9f},{data['max_error']:.9f},"
            f"{theta[-1, 0]:.9f},{theta[-1, 1]:.9f}"
        )
    Path("Chapter22_Lesson2_summary.csv").write_text(
        "\n".join(lines) + "\n", encoding="utf-8"
    )


def make_plots(results: dict[Mode, dict[str, np.ndarray | float]]) -> None:
    plt.figure(figsize=(10, 5))
    for mode, data in results.items():
        plt.plot(data["time"], data["error"], label=mode)
    plt.xlabel("Time (s)")
    plt.ylabel("Tracking error")
    plt.title("Tracking Error Comparison")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter22_Lesson2_tracking_error.png", dpi=180)

    plt.figure(figsize=(10, 5))
    for mode, data in results.items():
        theta = np.asarray(data["theta"])
        plt.plot(data["time"], theta[:, 0], label=f"{mode}: theta_x")
    plt.axhline(-2.5, linestyle="--", label="lower projection bound")
    plt.axhline(0.5, linestyle="--", label="upper projection bound")
    plt.xlabel("Time (s)")
    plt.ylabel("theta_x estimate")
    plt.title("Parameter Drift and Robustification")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter22_Lesson2_theta_x.png", dpi=180)


if __name__ == "__main__":
    configuration = Config()
    modes: tuple[Mode, ...] = ("baseline", "sigma", "projection", "combined")
    simulation_results = {
        mode: simulate(mode, configuration) for mode in modes
    }

    for mode, data in simulation_results.items():
        final_theta = np.asarray(data["theta"])[-1]
        print(
            f"{mode:10s} | RMS={data['rms_error']:.6f} | "
            f"max|e|={data['max_error']:.6f} | "
            f"theta_final={final_theta}"
        )

    save_summary(simulation_results)
    make_plots(simulation_results)
    plt.show()

Chapter22_Lesson2.cpp

The C++17 implementation uses only the standard library and exports a CSV file. For larger parameter vectors, replace std::array with Eigen or Armadillo vectors and matrices.

// Chapter22_Lesson2.cpp
// Build: g++ -std=c++17 -O2 Chapter22_Lesson2.cpp -o Chapter22_Lesson2
// Run:   ./Chapter22_Lesson2

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

enum class Mode { Baseline, Sigma, Projection, Combined };

struct Config {
    double dt = 0.002;
    double finalTime = 120.0;
    double a = -0.4;
    double b = 1.0;
    double am = -1.5;
    double bm = 1.5;
    double q = 1.0;
    std::array<double, 2> gamma{12.0, 12.0};
    double sigma = 0.03;
    std::array<double, 2> theta0{-0.8, 1.2};
    std::array<double, 2> lower{-2.5, 0.2};
    std::array<double, 2> upper{0.5, 2.5};

    double p() const {
        return -q / (2.0 * am);
    }
};

struct Sample {
    double time{};
    double x{};
    double xm{};
    double error{};
    double control{};
    std::array<double, 2> theta{};
};

struct Result {
    std::vector<Sample> samples;
    double rmsError{};
    double maxError{};
};

double reference(double t) {
    return (t < 20.0) ? std::sin(0.8 * t) : 0.0;
}

double measurementNoise(double t) {
    return 0.06 * std::sin(37.0 * t) + 0.03 * std::sin(71.0 * t);
}

double disturbance(double t) {
    return 0.25 * std::sin(12.0 * t) + 0.12 * std::sin(0.17 * t) + 0.12;
}

std::string modeName(Mode mode) {
    switch (mode) {
        case Mode::Baseline: return "baseline";
        case Mode::Sigma: return "sigma";
        case Mode::Projection: return "projection";
        case Mode::Combined: return "combined";
    }
    throw std::logic_error("Unknown mode");
}

std::array<double, 2> projectBox(
    const std::array<double, 2>& theta,
    std::array<double, 2> velocity,
    const Config& cfg
) {
    for (std::size_t i = 0; i < theta.size(); ++i) {
        if (theta[i] <= cfg.lower[i] && velocity[i] < 0.0) {
            velocity[i] = 0.0;
        } else if (theta[i] >= cfg.upper[i] && velocity[i] > 0.0) {
            velocity[i] = 0.0;
        }
    }
    return velocity;
}

Result simulate(Mode mode, const Config& cfg) {
    const std::size_t steps =
        static_cast<std::size_t>(std::llround(cfg.finalTime / cfg.dt)) + 1U;

    Result result;
    result.samples.reserve(steps);

    double x = 0.0;
    double xm = 0.0;
    std::array<double, 2> theta{0.0, 0.0};
    double squaredErrorSum = 0.0;
    double maxError = 0.0;

    for (std::size_t k = 0; k < steps; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        const double r = reference(t);
        const double xMeasured = x + measurementNoise(t);
        const double measuredError = xMeasured - xm;
        const std::array<double, 2> phi{xMeasured, r};
        const double control = theta[0] * phi[0] + theta[1] * phi[1];
        const double trueError = x - xm;

        result.samples.push_back({t, x, xm, trueError, control, theta});
        squaredErrorSum += trueError * trueError;
        maxError = std::max(maxError, std::abs(trueError));

        if (k + 1U == steps) {
            break;
        }

        const double xDot = cfg.a * x + cfg.b * control + disturbance(t);
        const double xmDot = cfg.am * xm + cfg.bm * r;

        std::array<double, 2> velocity{};
        for (std::size_t i = 0; i < velocity.size(); ++i) {
            velocity[i] = -cfg.gamma[i] * phi[i] * cfg.b * cfg.p() * measuredError;
        }

        if (mode == Mode::Sigma || mode == Mode::Combined) {
            for (std::size_t i = 0; i < velocity.size(); ++i) {
                velocity[i] -= cfg.sigma * cfg.gamma[i] * (theta[i] - cfg.theta0[i]);
            }
        }

        if (mode == Mode::Projection || mode == Mode::Combined) {
            velocity = projectBox(theta, velocity, cfg);
        }

        x += cfg.dt * xDot;
        xm += cfg.dt * xmDot;
        for (std::size_t i = 0; i < theta.size(); ++i) {
            theta[i] += cfg.dt * velocity[i];
        }

        if (mode == Mode::Projection || mode == Mode::Combined) {
            for (std::size_t i = 0; i < theta.size(); ++i) {
                theta[i] = std::clamp(theta[i], cfg.lower[i], cfg.upper[i]);
            }
        }
    }

    result.rmsError = std::sqrt(squaredErrorSum / static_cast<double>(steps));
    result.maxError = maxError;
    return result;
}

int main() {
    const Config cfg;
    const std::array<Mode, 4> modes{
        Mode::Baseline, Mode::Sigma, Mode::Projection, Mode::Combined
    };

    std::ofstream csv("Chapter22_Lesson2_results_cpp.csv");
    if (!csv) {
        std::cerr << "Cannot create output CSV file.\n";
        return 1;
    }
    csv << "mode,time,x,xm,error,control,theta_x,theta_r\n";
    csv << std::setprecision(10);

    for (const Mode mode : modes) {
        const Result result = simulate(mode, cfg);
        const Sample& last = result.samples.back();

        std::cout << std::setw(10) << modeName(mode)
                  << " | RMS=" << result.rmsError
                  << " | max|e|=" << result.maxError
                  << " | theta_final=[" << last.theta[0]
                  << ", " << last.theta[1] << "]\n";

        for (const Sample& sample : result.samples) {
            csv << modeName(mode) << ','
                << sample.time << ','
                << sample.x << ','
                << sample.xm << ','
                << sample.error << ','
                << sample.control << ','
                << sample.theta[0] << ','
                << sample.theta[1] << '\n';
        }
    }

    return 0;
}

Chapter22_Lesson2.java

The Java implementation uses records for immutable samples and writes results to CSV. Apache Commons Math can replace the primitive arrays in multivariable designs.

// Chapter22_Lesson2.java
// Build: javac Chapter22_Lesson2.java
// Run:   java Chapter22_Lesson2

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 Chapter22_Lesson2 {
    private enum Mode {
        BASELINE("baseline"),
        SIGMA("sigma"),
        PROJECTION("projection"),
        COMBINED("combined");

        private final String label;

        Mode(String label) {
            this.label = label;
        }
    }

    private static final class Config {
        final double dt = 0.002;
        final double finalTime = 120.0;
        final double a = -0.4;
        final double b = 1.0;
        final double am = -1.5;
        final double bm = 1.5;
        final double q = 1.0;
        final double[] gamma = {12.0, 12.0};
        final double sigma = 0.03;
        final double[] theta0 = {-0.8, 1.2};
        final double[] lower = {-2.5, 0.2};
        final double[] upper = {0.5, 2.5};

        double p() {
            return -q / (2.0 * am);
        }
    }

    private record Sample(
        double time,
        double x,
        double xm,
        double error,
        double control,
        double thetaX,
        double thetaR
    ) {}

    private record Result(List<Sample> samples, double rmsError, double maxError) {}

    private static double reference(double t) {
        return t < 20.0 ? Math.sin(0.8 * t) : 0.0;
    }

    private static double measurementNoise(double t) {
        return 0.06 * Math.sin(37.0 * t) + 0.03 * Math.sin(71.0 * t);
    }

    private static double disturbance(double t) {
        return 0.25 * Math.sin(12.0 * t)
            + 0.12 * Math.sin(0.17 * t)
            + 0.12;
    }

    private static double[] projectBox(
        double[] theta,
        double[] velocity,
        Config cfg
    ) {
        double[] projected = velocity.clone();
        for (int i = 0; i < theta.length; i++) {
            if (theta[i] <= cfg.lower[i] && projected[i] < 0.0) {
                projected[i] = 0.0;
            } else if (theta[i] >= cfg.upper[i] && projected[i] > 0.0) {
                projected[i] = 0.0;
            }
        }
        return projected;
    }

    private static double clamp(double value, double lower, double upper) {
        return Math.max(lower, Math.min(upper, value));
    }

    private static Result simulate(Mode mode, Config cfg) {
        int steps = (int) Math.round(cfg.finalTime / cfg.dt) + 1;
        List<Sample> samples = new ArrayList<>(steps);

        double x = 0.0;
        double xm = 0.0;
        double[] theta = {0.0, 0.0};
        double squaredErrorSum = 0.0;
        double maxError = 0.0;

        for (int k = 0; k < steps; k++) {
            double t = k * cfg.dt;
            double r = reference(t);
            double xMeasured = x + measurementNoise(t);
            double measuredError = xMeasured - xm;
            double[] phi = {xMeasured, r};
            double control = theta[0] * phi[0] + theta[1] * phi[1];
            double trueError = x - xm;

            samples.add(new Sample(
                t, x, xm, trueError, control, theta[0], theta[1]
            ));
            squaredErrorSum += trueError * trueError;
            maxError = Math.max(maxError, Math.abs(trueError));

            if (k + 1 == steps) {
                break;
            }

            double xDot = cfg.a * x + cfg.b * control + disturbance(t);
            double xmDot = cfg.am * xm + cfg.bm * r;

            double[] velocity = new double[2];
            for (int i = 0; i < velocity.length; i++) {
                velocity[i] = -cfg.gamma[i] * phi[i] * cfg.b
                    * cfg.p() * measuredError;
            }

            if (mode == Mode.SIGMA || mode == Mode.COMBINED) {
                for (int i = 0; i < velocity.length; i++) {
                    velocity[i] -= cfg.sigma * cfg.gamma[i]
                        * (theta[i] - cfg.theta0[i]);
                }
            }

            if (mode == Mode.PROJECTION || mode == Mode.COMBINED) {
                velocity = projectBox(theta, velocity, cfg);
            }

            x += cfg.dt * xDot;
            xm += cfg.dt * xmDot;
            for (int i = 0; i < theta.length; i++) {
                theta[i] += cfg.dt * velocity[i];
            }

            if (mode == Mode.PROJECTION || mode == Mode.COMBINED) {
                for (int i = 0; i < theta.length; i++) {
                    theta[i] = clamp(theta[i], cfg.lower[i], cfg.upper[i]);
                }
            }
        }

        double rmsError = Math.sqrt(squaredErrorSum / steps);
        return new Result(samples, rmsError, maxError);
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        Config cfg = new Config();
        Path output = Path.of("Chapter22_Lesson2_results_java.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(
            output, StandardCharsets.UTF_8
        )) {
            writer.write("mode,time,x,xm,error,control,theta_x,theta_r\n");

            for (Mode mode : Mode.values()) {
                Result result = simulate(mode, cfg);
                Sample last = result.samples().get(result.samples().size() - 1);

                System.out.printf(
                    "%10s | RMS=%.6f | max|e|=%.6f | theta_final=[%.6f, %.6f]%n",
                    mode.label,
                    result.rmsError(),
                    result.maxError(),
                    last.thetaX(),
                    last.thetaR()
                );

                for (Sample sample : result.samples()) {
                    writer.write(String.format(
                        Locale.US,
                        "%s,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f%n",
                        mode.label,
                        sample.time(),
                        sample.x(),
                        sample.xm(),
                        sample.error(),
                        sample.control(),
                        sample.thetaX(),
                        sample.thetaR()
                    ));
                }
            }
        }
    }
}

Chapter22_Lesson2.m

The script requires base MATLAB. In Simulink, represent the plant and reference model with Integrator and Gain blocks, compute the gradient and leakage in Product/Sum blocks, and implement tangent-cone projection in a MATLAB Function block. A Saturation block after the parameter integrator provides the finite-step safety clip.

% Chapter22_Lesson2.m
% Robust adaptive-law comparison for a first-order MRAC example.
% The script requires only base MATLAB. For Simulink, implement the same
% equations with Integrator, Sum, Gain, Product, Saturation, and MATLAB
% Function blocks; place the projection logic in the MATLAB Function block.

clear; clc; close all;

cfg.dt = 0.002;
cfg.finalTime = 120.0;
cfg.a = -0.4;
cfg.b = 1.0;
cfg.am = -1.5;
cfg.bm = 1.5;
cfg.q = 1.0;
cfg.p = -cfg.q / (2.0 * cfg.am);
cfg.gamma = [12.0; 12.0];
cfg.sigma = 0.03;
cfg.theta0 = [-0.8; 1.2];
cfg.lower = [-2.5; 0.2];
cfg.upper = [0.5; 2.5];

modes = {'baseline', 'sigma', 'projection', 'combined'};
results = struct();

for i = 1:numel(modes)
    mode = modes{i};
    results.(mode) = simulateAdaptiveLaw(mode, cfg);
    finalTheta = results.(mode).theta(:, end);
    fprintf('%10s | RMS=%.6f | max|e|=%.6f | theta_final=[%.6f, %.6f]\n', ...
        mode, results.(mode).rmsError, results.(mode).maxError, ...
        finalTheta(1), finalTheta(2));
end

figure('Name', 'Tracking error comparison');
hold on;
for i = 1:numel(modes)
    mode = modes{i};
    plot(results.(mode).time, results.(mode).error, ...
        'DisplayName', mode);
end
xlabel('Time (s)');
ylabel('Tracking error');
title('Tracking Error Comparison');
grid on; legend('Location', 'best');

figure('Name', 'Parameter drift and robustification');
hold on;
for i = 1:numel(modes)
    mode = modes{i};
    plot(results.(mode).time, results.(mode).theta(1, :), ...
        'DisplayName', [mode, ': theta_x']);
end
yline(cfg.lower(1), '--', 'Lower projection bound');
yline(cfg.upper(1), '--', 'Upper projection bound');
xlabel('Time (s)');
ylabel('theta_x estimate');
title('Parameter Drift and Robustification');
grid on; legend('Location', 'best');

summaryTable = table( ...
    string(modes(:)), ...
    cellfun(@(m) results.(m).rmsError, modes(:)), ...
    cellfun(@(m) results.(m).maxError, modes(:)), ...
    cellfun(@(m) results.(m).theta(1, end), modes(:)), ...
    cellfun(@(m) results.(m).theta(2, end), modes(:)), ...
    'VariableNames', {'Mode', 'RMSError', 'MaxAbsError', ...
    'ThetaXFinal', 'ThetaRFinal'});
writetable(summaryTable, 'Chapter22_Lesson2_summary_matlab.csv');

function result = simulateAdaptiveLaw(mode, cfg)
    steps = round(cfg.finalTime / cfg.dt) + 1;
    time = linspace(0.0, cfg.finalTime, steps);
    x = zeros(1, steps);
    xm = zeros(1, steps);
    u = zeros(1, steps);
    theta = zeros(2, steps);

    for k = 1:(steps - 1)
        t = time(k);
        r = referenceSignal(t);
        xMeasured = x(k) + measurementNoise(t);
        measuredError = xMeasured - xm(k);
        phi = [xMeasured; r];
        u(k) = theta(:, k).' * phi;

        xDot = cfg.a * x(k) + cfg.b * u(k) + disturbanceSignal(t);
        xmDot = cfg.am * xm(k) + cfg.bm * r;

        velocity = -cfg.gamma .* phi * cfg.b * cfg.p * measuredError;

        if strcmp(mode, 'sigma') || strcmp(mode, 'combined')
            velocity = velocity ...
                - cfg.sigma * cfg.gamma .* (theta(:, k) - cfg.theta0);
        end

        if strcmp(mode, 'projection') || strcmp(mode, 'combined')
            velocity = projectBox(theta(:, k), velocity, ...
                cfg.lower, cfg.upper);
        end

        x(k + 1) = x(k) + cfg.dt * xDot;
        xm(k + 1) = xm(k) + cfg.dt * xmDot;
        theta(:, k + 1) = theta(:, k) + cfg.dt * velocity;

        if strcmp(mode, 'projection') || strcmp(mode, 'combined')
            theta(:, k + 1) = min(max(theta(:, k + 1), ...
                cfg.lower), cfg.upper);
        end
    end

    finalPhi = [x(end) + measurementNoise(time(end)); ...
        referenceSignal(time(end))];
    u(end) = theta(:, end).' * finalPhi;
    error = x - xm;

    result.time = time;
    result.x = x;
    result.xm = xm;
    result.u = u;
    result.theta = theta;
    result.error = error;
    result.rmsError = sqrt(mean(error .^ 2));
    result.maxError = max(abs(error));
end

function projected = projectBox(theta, velocity, lower, upper)
    projected = velocity;
    for i = 1:numel(theta)
        if theta(i) <= lower(i) && projected(i) < 0.0
            projected(i) = 0.0;
        elseif theta(i) >= upper(i) && projected(i) > 0.0
            projected(i) = 0.0;
        end
    end
end

function r = referenceSignal(t)
    if t < 20.0
        r = sin(0.8 * t);
    else
        r = 0.0;
    end
end

function n = measurementNoise(t)
    n = 0.06 * sin(37.0 * t) + 0.03 * sin(71.0 * t);
end

function d = disturbanceSignal(t)
    d = 0.25 * sin(12.0 * t) + 0.12 * sin(0.17 * t) + 0.12;
end

Chapter22_Lesson2.nb

The notebook expression contains one executable input cell. It uses associations and NestList to reproduce the same sampled-data algorithm and exports a summary CSV file.

Notebook[{Cell["Chapter 22, Lesson 2: Robust Adaptive Laws", "Title"],
  Cell["Leakage, projection, and sigma modification", "Subtitle"],
  Cell[BoxData[
"ClearAll[\"Global`*\"];
cfg = <|\"dt\" -> 0.002, \"finalTime\" -> 120.0, \"a\" -> -0.4, \"b\" -> 1.0,
  \"am\" -> -1.5, \"bm\" -> 1.5, \"q\" -> 1.0, \"gamma\" -> {12.0, 12.0},
  \"sigma\" -> 0.03, \"theta0\" -> {-0.8, 1.2}, \"lower\" -> {-2.5, 0.2},
  \"upper\" -> {0.5, 2.5}|>;
p = -cfg[\"q\"]/(2.0 cfg[\"am\"]);
reference[t_] := If[t < 20.0, Sin[0.8 t], 0.0];
measurementNoise[t_] := 0.06 Sin[37.0 t] + 0.03 Sin[71.0 t]; disturbance[t_] := 0.25 Sin[12.0 t] + 0.12 Sin[0.17 t] + 0.12;
projectBox[theta_, velocity_] := MapThread[
  Which[#1 <= #3 && #2 < 0.0, 0.0, #1 >= #4 && #2 > 0.0, 0.0, True, #2] &,
  {theta, velocity, cfg[\"lower\"], cfg[\"upper\"]}];
simulate[mode_String] := Module[{steps, initial, trajectory, update, rows, error},
  steps = Round[cfg[\"finalTime\"]/cfg[\"dt\"]];
  initial = <|\"t\" -> 0.0, \"x\" -> 0.0, \"xm\" -> 0.0,
    \"theta\" -> {0.0, 0.0}, \"u\" -> 0.0|>;
  update[state_] := Module[
    {t, x, xm, theta, r, xMeasured, eMeasured, phi, u, xDot, xmDot, velocity, thetaNext},
    t = state[\"t\"]; x = state[\"x\"]; xm = state[\"xm\"];
    theta = state[\"theta\"]; r = reference[t];
    xMeasured = x + measurementNoise[t]; eMeasured = xMeasured - xm;
    phi = {xMeasured, r}; u = theta.phi;
    xDot = cfg[\"a\"] x + cfg[\"b\"] u + disturbance[t];
    xmDot = cfg[\"am\"] xm + cfg[\"bm\"] r;
    velocity = -cfg[\"gamma\"] phi cfg[\"b\"] p eMeasured;
    If[MemberQ[{\"sigma\", \"combined\"}, mode],
      velocity = velocity - cfg[\"sigma\"] cfg[\"gamma\"] (theta - cfg[\"theta0\"])];
    If[MemberQ[{\"projection\", \"combined\"}, mode],
      velocity = projectBox[theta, velocity]];
    thetaNext = theta + cfg[\"dt\"] velocity;
    If[MemberQ[{\"projection\", \"combined\"}, mode],
      thetaNext = MapThread[Clip[#1, { #2, #3 }] &,
        {thetaNext, cfg[\"lower\"], cfg[\"upper\"]}]];
    <|\"t\" -> t + cfg[\"dt\"], \"x\" -> x + cfg[\"dt\"] xDot,
      \"xm\" -> xm + cfg[\"dt\"] xmDot, \"theta\" -> thetaNext, \"u\" -> u|>];
  trajectory = NestList[update, initial, steps];
  rows = ({ [\"t\"], #[\"x\"], #[\"xm\"], #[\"x\"] - #[\"xm\"],
      #[\"u\"], #[\"theta\"][[1]], #[\"theta\"][[2]] } &) /@ trajectory;
  error = rows[[All, 4]];
  <|\"rows\" -> rows, \"rmsError\" -> Sqrt[Mean[error^2]],
    \"maxError\" -> Max[Abs[error]], \"finalTheta\" -> trajectory[[-1]][\"theta\"]|>];
modes = {\"baseline\", \"sigma\", \"projection\", \"combined\"}; results = AssociationMap[simulate, modes];
summary = Table[{mode, results[mode][\"rmsError\"], results[mode][\"maxError\"],
  results[mode][\"finalTheta\"][[1]], results[mode][\"finalTheta\"][[2]]}, {mode, modes}];
Print[Grid[Prepend[summary, {\"Mode\", \"RMS Error\", \"Max |e|\",
  \"theta_x final\", \"theta_r final\"}], Frame -> All]];
Export[\"Chapter22_Lesson2_summary_wolfram.csv\", Prepend[summary,
  {\"mode\", \"rms_error\", \"max_abs_error\", \"theta_x_final\", \"theta_r_final\"}]];
trackingPlot = ListLinePlot[
  Table[results[mode][\"rows\"][[All, {1, 4}]], {mode, modes}],
  PlotLegends -> modes, Frame -> True, FrameLabel -> {\"Time (s)\", \"Tracking error\"},
  PlotLabel -> \"Tracking Error Comparison\", ImageSize -> Large];
parameterPlot = ListLinePlot[
  Table[results[mode][\"rows\"][[All, {1, 6}]], {mode, modes}],
  PlotLegends -> (StringJoin[#, \" theta_x\"] & /@ modes), Frame -> True,
  FrameLabel -> {\"Time (s)\", \"theta_x estimate\"}, PlotLabel -> \"Parameter Drift and Robustification\", ImageSize -> Large];
Column[{trackingPlot, parameterPlot}]
"], "Input"]}, WindowTitle -> "Chapter22_Lesson2"]        

With the chosen deterministic signals, the baseline controller usually achieves small tracking error but its estimate \( \hat\theta_x \) continues to move after excitation is removed. Projection clips this drift at the lower admissible bound. Leakage keeps the estimate near the nominal vector and reduces the largest error excursion, but introduces a modest tracking bias. The combined law inherits the nominal attraction of leakage and the hard safety guarantee of projection.

These observations illustrate a central robust-adaptive trade-off: exact ideal-model convergence is weakened in exchange for bounded estimates and predictable behavior under nonideal signals. The result is not that one modification is universally best. The correct design depends on the quality of prior parameter bounds, the disturbance level, the available excitation, and the acceptable steady-state bias.

Problem 1: Matching Parameters and Gradient Cancellation

For \( \dot x=ax+bu+d(t) \), reference model \( \dot x_m=a_mx_m+b_mr \), and controller \( u=\hat\theta_xx+\hat\theta_rr \), derive the ideal parameters and the baseline gradient law that cancels the state–parameter cross term.

Solution:

Matching the plant coefficients to the reference model requires

\[ a+b\theta_x^{*}=a_m,\qquad b\theta_r^{*}=b_m, \]

\[ \theta_x^{*}=\frac{a_m-a}{b},\qquad \theta_r^{*}=\frac{b_m}{b}. \]

The error dynamics are

\[ \dot e=a_me+b\tilde{\boldsymbol\theta}^{T} \boldsymbol\phi+d(t). \]

For \( V=Pe^2+\tilde{\boldsymbol\theta}^{T} \Gamma^{-1}\tilde{\boldsymbol\theta} \), the cross term is \( 2Peb\tilde{\boldsymbol\theta}^{T}\boldsymbol\phi \). Choosing

\[ \dot{\hat{\boldsymbol\theta} } =-\Gamma\boldsymbol\phi bPe \]

produces its exact negative in \( \dot V \).

Problem 2: Scalar Projection Inequality

Let the true scalar parameter satisfy \( \theta^{*}\in[\underline\theta,\overline\theta] \). Prove that the interval projection satisfies

\[ (\hat\theta-\theta^{*}) \left[\operatorname{Proj}(\hat\theta,v)-v\right]\leq0. \]

Solution:

Inside the interval, projection returns \( v \), so the product is zero. At the lower bound with \( v<0 \), the projected velocity is zero. Since \( \hat\theta-\theta^{*} =\underline\theta-\theta^{*}\leq0 \) and \( 0-v>0 \), their product is nonpositive. At the upper bound with \( v>0 \), both \( \hat\theta-\theta^{*}\geq0 \) and \( 0-v<0 \), so the product is again nonpositive. All remaining boundary directions are inward and are not modified.

Problem 3: Numerical Ultimate-Bound Estimate

Use \( Q=1 \), \( P=1/3 \), \( \Gamma=12I \), \( \sigma=0.03 \), \( \bar d=0.5 \), and \( \boldsymbol\delta_0=[-0.3,0.3]^T \). Compute the conservative tracking-error ultimate bound derived in Section 7.

Solution:

\[ \|\boldsymbol\delta_0\|^2=0.18, \qquad \alpha=\min\left\{\frac{1}{2/3},0.03(12)\right\}=0.36. \]

\[ c_0=0.03(0.18)+2\left(\frac{1}{3}\right)^2(0.5)^2 \approx0.06096. \]

\[ \limsup V\leq\frac{0.06096}{0.36}\approx0.1693, \]

\[ \limsup|e|\leq \sqrt{\frac{0.1693}{1/3} }\approx0.713. \]

This analytical number is intentionally conservative because Young’s inequality replaces time-varying signal information by worst-case constants.

Problem 4: Why a Nominal Leakage Center Matters

Compare leakage toward zero with leakage toward \( \boldsymbol\theta_0 \). Which choice gives the smaller theoretical bias term?

Solution:

Leakage toward zero yields a residual term proportional to \( \sigma\|\boldsymbol\theta^{*}\|^2 \). Leakage toward \( \boldsymbol\theta_0 \) yields \( \sigma\|\boldsymbol\theta^{*}- \boldsymbol\theta_0\|^2 \). Therefore, any nominal vector closer to the ideal vector than the origin reduces the Lyapunov ultimate-bound constant. A poor nominal model can do the opposite.

Problem 5: Sampled Projection and Boundary Overshoot

An Euler update uses \( \hat\theta_{k+1}=\hat\theta_k+h \operatorname{Proj}(\hat\theta_k,v_k) \). Explain why a final clip is still useful even though the continuous-time projected dynamics preserve the interval.

Solution:

Projection evaluates the active constraint at the beginning of the sample. If \( \hat\theta_k \) lies slightly inside a bound and \( hv_k \) is large, the finite step can cross the boundary before projection becomes active at the next sample. The operation

\[ \hat\theta_{k+1}\leftarrow \min\left(\overline\theta, \max(\underline\theta,\hat\theta_{k+1})\right) \]

restores discrete invariance. It does not replace proper gain and sampling-time selection; frequent large clips indicate poor numerical or control tuning.

The baseline adaptive law cancels the ideal state–parameter cross term but can integrate persistent nonideal error and drift. Leakage adds dissipative forgetting toward a nominal vector; σ-modification is its constant isotropic form. Projection imposes a known convex parameter constraint and has a nonpositive Lyapunov correction. Combining both produces hard parameter bounds and the comparison inequality \( \dot V\leq-\alpha V+c_0 \), establishing uniform ultimate boundedness under bounded disturbances. The price is a deliberate robustness–accuracy trade-off: parameter convergence and exact asymptotic tracking may be replaced by predictable residual bounds.