After completing this lesson, students should be able to:

• distinguish commanded input from actuator-delivered input;
• derive the saturation mismatch term in MRAC error dynamics;
• explain why the nominal update can produce parameter windup;
• construct an auxiliary saturation-error model;
• prove boundedness for the compensated adaptive subsystem;
• implement gating, projection, and saturation-activated leakage;
• compare naive, frozen, and compensated adaptation numerically.

Students are assumed to know linear state-space control, Lyapunov stability, scalar and vector MRAC error models, robust modifications, and parameter projection from earlier chapters.

Consider the uncertain scalar plant

\[ \dot x = a x + b u_a,\qquad b\neq0. \]

The controller computes \(u_c\), whereas the actuator applies \(u_a\). For symmetric amplitude saturation,

\[ u_a=\operatorname{sat}_{u_{\max}}(u_c) = \begin{cases} -u_{\max}, & u_c < -u_{\max},\\ u_c, & |u_c|\le u_{\max},\\ u_{\max}, & u_c > u_{\max}. \end{cases} \]

Define the measurable actuator mismatch

\[ \Delta u=u_a-u_c,\qquad u_a=u_c+\Delta u. \]

Saturation is not merely an independent bounded disturbance: \(\Delta u\) depends on the adaptive command and hence on the parameter estimates. An update may increase \(|u_c|\) while \(|u_a|\) remains fixed, making the mismatch larger.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> AC["Adaptive controller"]
  X["Plant state x"] --> AC
  AC --> UC["Command u_c"]
  UC --> SAT["Actuator saturation"]
  SAT --> UA["Applied input u_a"]
  UA --> P["Uncertain plant"]
  P --> X
  UC --> M["Mismatch: delta_u = u_a - u_c"]
  UA --> M
  M --> AUX["Stable saturation-error model"]
  AUX --> EC["Compensated error e_c"]
  RM --> EC
  X --> EC
  EC --> UP["Modified parameter update"]
  UP --> AC
        

Let the desired reference model and adaptive controller be

\[ \dot x_m=a_mx_m+b_mr,\qquad a_m<0, \]

\[ u_c=\hat\theta_xx+\hat\theta_rr =\hat{\boldsymbol\theta}^{T}\boldsymbol\phi,\qquad \boldsymbol\phi=\begin{bmatrix}x&r\end{bmatrix}^{T}. \]

Assume the ideal matching conditions

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

With \(e=x-x_m\) and \(\widetilde{\boldsymbol\theta} =\hat{\boldsymbol\theta}-\boldsymbol\theta^{*}\),

\[ \dot e=a_me+ b\widetilde{\boldsymbol\theta}^{T}\boldsymbol\phi+ b\Delta u. \]

Choose \(P>0\) from

\[ 2a_mP=-Q,\qquad Q>0, \]

and use the nominal update

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

For

\[ V=Pe^2+ \widetilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1} \widetilde{\boldsymbol\theta}, \]

differentiation gives

\[ \dot V=-Qe^2+2Pbe\Delta u. \]

The term \(2Pbe\Delta u\) has no fixed sign. Therefore the nominal proof no longer establishes monotonic decrease of \(V\). The raw error combines parameter mismatch with actuator infeasibility, so the estimator may try to explain an error that no parameter value can remove while saturation persists.

Introduce a stable auxiliary state:

\[ \dot e_\Delta=a_me_\Delta+b\Delta u,\qquad e_\Delta(0)=0. \]

Define \(e_c=e-e_\Delta\). Then

\[ \dot e_c=a_me_c+ b\widetilde{\boldsymbol\theta}^{T}\boldsymbol\phi. \]

The actuator mismatch is removed from the error channel used for adaptation. The physical tracking error remains \(e=e_c+e_\Delta\); compensation does not claim that an infeasible reference is trackable.

4.1 Lyapunov theorem

Use

\[ \dot{\hat{\boldsymbol\theta}} =-\boldsymbol\Gamma\boldsymbol\phi\,bPe_c. \]

For

\[ V_c=Pe_c^2+ \widetilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1} \widetilde{\boldsymbol\theta}, \]

\[ \begin{aligned} \dot V_c &= 2Pe_c\left( a_me_c+ b\widetilde{\boldsymbol\theta}^{T}\boldsymbol\phi \right) + 2\widetilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1} \dot{\hat{\boldsymbol\theta}}\\ &= 2Pa_me_c^2+ 2Pbe_c\widetilde{\boldsymbol\theta}^{T}\boldsymbol\phi- 2Pbe_c\widetilde{\boldsymbol\theta}^{T}\boldsymbol\phi\\ &=-Qe_c^2\le0. \end{aligned} \]

Hence \(e_c\) and the parameter error are bounded. Under the usual bounded-regressor assumptions, \(e_c\in L_2\cap L_\infty\); if \(\dot e_c\) is bounded, Barbalat's lemma gives \(e_c(t)\to0\). The physical error also requires \(e_\Delta\to0\), which follows after saturation disappears and the stable auxiliary model decays.

A hard gate freezes adaptation during detected saturation:

\[ \dot{\hat{\boldsymbol\theta}} =-\rho(t)\boldsymbol\Gamma\boldsymbol\phi\,bPe, \]

\[ \rho(t)= \begin{cases} 1,&|\Delta u(t)|\le\varepsilon_u,\\ 0,&|\Delta u(t)|>\varepsilon_u. \end{cases} \]

This prevents obvious parameter windup but may suppress useful learning during long saturation intervals.

5.1 Directional gating

Let \(\boldsymbol\nu= -\boldsymbol\Gamma\boldsymbol\phi bPe\). Since the parameter-update contribution to the command derivative is

\[ \dot u_c\big|_{\text{parameter update}} =\boldsymbol\phi^T\boldsymbol\nu, \]

a practical directional gate is

\[ \rho_{\mathrm{dir}}= \begin{cases} 0,& |u_c|>u_{\max}-\varepsilon_u \text{ and } u_c\boldsymbol\phi^T\boldsymbol\nu>0,\\ 1,&\text{otherwise}. \end{cases} \]

This rule blocks only updates that push the command farther outward. A hysteresis band is recommended to prevent switching chatter.

5.2 Smooth gating

\[ \rho_s(\Delta u) = \frac{1}{ 1+\left(|\Delta u|/\delta_s\right)^p}, \qquad p\ge2. \]

Smooth gating reduces derivative discontinuities but is not by itself a Lyapunov proof; it must be analyzed as part of the complete modified dynamics.

Let \(\Omega_\theta\) be a known convex parameter set. The projected compensated update is

\[ \dot{\hat{\boldsymbol\theta}} = \operatorname{Proj}_{\Omega_\theta} \left( \hat{\boldsymbol\theta}, -\boldsymbol\Gamma\boldsymbol\phi bPe_c \right). \]

The standard projection inequality is

\[ \widetilde{\boldsymbol\theta}^T \boldsymbol\Gamma^{-1} \left[ \operatorname{Proj}( \hat{\boldsymbol\theta},\boldsymbol\nu) -\boldsymbol\nu \right]\le0, \]

provided the ideal parameter lies in the admissible set. Thus projection does not spoil the Lyapunov upper bound.

6.1 Saturation-activated leakage

\[ s_\Delta(t)= \min\left(1,\frac{|\Delta u(t)|}{\delta_\Delta}\right), \]

\[ \dot{\hat{\boldsymbol\theta}} = \operatorname{Proj}_{\Omega_\theta} \left( \hat{\boldsymbol\theta}, -\boldsymbol\Gamma\boldsymbol\phi bPe_c - \sigma_ss_\Delta\boldsymbol\Gamma (\hat{\boldsymbol\theta}-\boldsymbol\theta_0) \right). \]

Leakage limits accumulated parameter displacement during prolonged saturation. It introduces bias relative to the ideal parameter. A representative bound is

\[ \dot V_c \le -Qe_c^2 -\sigma_ss_\Delta \|\widetilde{\boldsymbol\theta}\|^2 + \sigma_ss_\Delta \|\boldsymbol\theta^{*}-\boldsymbol\theta_0\|^2. \]

The ultimate parameter-error radius therefore depends on the distance between the ideal parameter and the leakage anchor.

\[ \begin{aligned} u_c&=\hat{\boldsymbol\theta}^T\boldsymbol\phi,\\ u_a&=\operatorname{sat}_{u_{\max}}(u_c),\\ \Delta u&=u_a-u_c,\\ \dot e_\Delta&=a_me_\Delta+b\Delta u,\\ e_c&=(x-x_m)-e_\Delta,\\ \boldsymbol\nu&= -\rho\boldsymbol\Gamma\boldsymbol\phi bPe_c -\sigma_ss_\Delta\boldsymbol\Gamma (\hat{\boldsymbol\theta}-\boldsymbol\theta_0),\\ \dot{\hat{\boldsymbol\theta}} &= \operatorname{Proj}_{\Omega_\theta} (\hat{\boldsymbol\theta},\boldsymbol\nu). \end{aligned} \]

With reliable mismatch compensation one may normally set \(\rho=1\). Gating remains useful when actuator mismatch is noisy, delayed, rate-limited, or only approximately known.

flowchart TD
  A["Measure x and r"] --> B["Compute command u_c"]
  B --> C["Apply actuator limits to obtain u_a"]
  C --> D["Compute delta_u = u_a - u_c"]
  D --> E["Update stable saturation-error state"]
  E --> F["Compute e_c = x - x_m - e_delta"]
  F --> G["Form nominal adaptive direction"]
  G --> H["Add saturation leakage if required"]
  H --> I["Apply parameter projection"]
  I --> J["Integrate parameter estimates"]
  J --> B
  D --> K["Log saturation duty and recovery"]
        

7.1 Sampled implementation

\[ e_\Delta[k+1] = e_\Delta[k]+T_s \left(a_me_\Delta[k]+b\Delta u[k]\right), \]

\[ \hat{\boldsymbol\theta}[k+1] = \Pi_{\Omega_\theta} \left( \hat{\boldsymbol\theta}[k]+T_s\boldsymbol\nu[k] \right). \]

The commanded and actuated inputs must be logged separately, otherwise the mismatch cannot be reconstructed.

This program compares naive adaptation, hard gating, and compensated adaptation. It exports numerical results and a figure.

Chapter19_Lesson2.py

"""
Chapter19_Lesson2.py
Saturation-aware adaptive update laws for a scalar MRAC example.

The script compares:
  1. naive       - nominal update law driven by the raw tracking error,
  2. gated       - adaptation is frozen while the actuator is saturated,
  3. compensated - an auxiliary saturation-error model supplies a
                   compensated error; projection and saturation-activated
                   leakage are also applied.

Dependencies:
    Python 3.10+, 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["naive", "gated", "compensated"]


@dataclass(frozen=True)
class Configuration:
    a: float = -0.4
    b: float = 1.0
    a_m: float = -1.2
    b_m: float = 1.2
    u_max: float = 1.0
    gamma_x: float = 4.0
    gamma_r: float = 4.0
    sigma_sat: float = 0.15
    theta_min: float = -5.0
    theta_max: float = 5.0
    dt: float = 1.0e-3
    final_time: float = 24.0

    @property
    def p(self) -> float:
        """Scalar solution of 2*a_m*p = -1."""
        return -1.0 / (2.0 * self.a_m)


def reference(t: float) -> float:
    """Command containing an intentionally infeasible interval."""
    if t < 1.0:
        return 0.0
    if t < 8.0:
        return 4.0
    if t < 16.0:
        return 0.6
    return 0.6 + 0.25 * np.sin(0.8 * (t - 16.0))


def saturate(value: float, limit: float) -> float:
    return float(np.clip(value, -limit, limit))


def componentwise_projection(
    theta: np.ndarray,
    direction: np.ndarray,
    lower: float,
    upper: float,
) -> np.ndarray:
    """
    Tangent-cone projection for a box constraint.

    At a lower bound, outward negative motion is removed.
    At an upper bound, outward positive motion is removed.
    """
    projected = direction.copy()
    for i in range(theta.size):
        if theta[i] <= lower and direction[i] < 0.0:
            projected[i] = 0.0
        elif theta[i] >= upper and direction[i] > 0.0:
            projected[i] = 0.0
    return projected


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

    x = np.zeros(n_steps)
    x_m = np.zeros(n_steps)
    e_delta = np.zeros(n_steps)
    theta = np.zeros((n_steps, 2))

    r = np.zeros(n_steps)
    u_c = np.zeros(n_steps)
    u_a = np.zeros(n_steps)
    delta_u = np.zeros(n_steps)
    e = np.zeros(n_steps)
    e_c = np.zeros(n_steps)
    rho = np.ones(n_steps)

    gamma = np.diag([cfg.gamma_x, cfg.gamma_r])
    theta_anchor = np.zeros(2)

    for k in range(n_steps - 1):
        r[k] = reference(t[k])
        phi = np.array([x[k], r[k]])

        u_c[k] = float(theta[k] @ phi)
        u_a[k] = saturate(u_c[k], cfg.u_max)
        delta_u[k] = u_a[k] - u_c[k]

        e[k] = x[k] - x_m[k]
        e_c[k] = e[k] - e_delta[k]
        saturated = abs(delta_u[k]) > 1.0e-12

        if mode == "naive":
            adaptation_error = e[k]
            rho[k] = 1.0
            leakage = 0.0
        elif mode == "gated":
            adaptation_error = e[k]
            rho[k] = 0.0 if saturated else 1.0
            leakage = 0.0
        elif mode == "compensated":
            adaptation_error = e_c[k]
            rho[k] = 1.0
            leakage = cfg.sigma_sat if saturated else 0.0
        else:
            raise ValueError(f"Unsupported mode: {mode}")

        nominal_direction = (
            -rho[k] * (gamma @ phi) * (cfg.p * cfg.b * adaptation_error)
        )
        leakage_direction = (
            -leakage * (gamma @ (theta[k] - theta_anchor))
        )
        theta_dot_raw = nominal_direction + leakage_direction
        theta_dot = componentwise_projection(
            theta[k], theta_dot_raw, cfg.theta_min, cfg.theta_max
        )

        x_dot = cfg.a * x[k] + cfg.b * u_a[k]
        x_m_dot = cfg.a_m * x_m[k] + cfg.b_m * r[k]
        e_delta_dot = cfg.a_m * e_delta[k] + cfg.b * delta_u[k]

        x[k + 1] = x[k] + cfg.dt * x_dot
        x_m[k + 1] = x_m[k] + cfg.dt * x_m_dot
        e_delta[k + 1] = e_delta[k] + cfg.dt * e_delta_dot
        theta[k + 1] = theta[k] + cfg.dt * theta_dot
        theta[k + 1] = np.clip(theta[k + 1], cfg.theta_min, cfg.theta_max)

    # Fill final algebraic samples.
    r[-1] = reference(t[-1])
    phi_final = np.array([x[-1], r[-1]])
    u_c[-1] = float(theta[-1] @ phi_final)
    u_a[-1] = saturate(u_c[-1], cfg.u_max)
    delta_u[-1] = u_a[-1] - u_c[-1]
    e[-1] = x[-1] - x_m[-1]
    e_c[-1] = e[-1] - e_delta[-1]
    rho[-1] = rho[-2]

    recovery = (t >= 8.0) & (t <= 14.0)
    parameter_norm = np.linalg.norm(theta, axis=1)
    metrics = {
        "rms_error": float(np.sqrt(np.mean(e**2))),
        "recovery_rms_error": float(np.sqrt(np.mean(e[recovery] ** 2))),
        "maximum_parameter_norm": float(np.max(parameter_norm)),
        "saturation_fraction": float(np.mean(np.abs(delta_u) > 1.0e-9)),
    }

    return {
        "time": t,
        "reference": r,
        "x": x,
        "x_m": x_m,
        "e": e,
        "e_c": e_c,
        "e_delta": e_delta,
        "u_c": u_c,
        "u_a": u_a,
        "delta_u": delta_u,
        "theta_x": theta[:, 0],
        "theta_r": theta[:, 1],
        "rho": rho,
        **metrics,
    }


def save_csv(results: dict[str, dict[str, np.ndarray | float]], path: Path) -> None:
    columns = ["time"]
    arrays = [results["naive"]["time"]]
    for mode, result in results.items():
        for key in (
            "reference",
            "x",
            "x_m",
            "e",
            "e_c",
            "u_c",
            "u_a",
            "delta_u",
            "theta_x",
            "theta_r",
            "rho",
        ):
            columns.append(f"{mode}_{key}")
            arrays.append(result[key])
    data = np.column_stack(arrays)
    np.savetxt(
        path,
        data,
        delimiter=",",
        header=",".join(columns),
        comments="",
    )


def plot_results(
    results: dict[str, dict[str, np.ndarray | float]],
    output_path: Path,
) -> None:
    fig, axes = plt.subplots(3, 1, figsize=(11, 10), sharex=True)

    for mode, result in results.items():
        axes[0].plot(result["time"], result["x"], label=f"x: {mode}")
    axes[0].plot(
        results["compensated"]["time"],
        results["compensated"]["x_m"],
        "k--",
        linewidth=1.5,
        label="reference model",
    )
    axes[0].set_ylabel("state")
    axes[0].grid(True)
    axes[0].legend(ncol=2)

    for mode, result in results.items():
        axes[1].plot(result["time"], result["u_a"], label=f"u_a: {mode}")
    axes[1].axhline(1.0, color="k", linestyle=":")
    axes[1].axhline(-1.0, color="k", linestyle=":")
    axes[1].set_ylabel("actuated input")
    axes[1].grid(True)
    axes[1].legend(ncol=3)

    for mode, result in results.items():
        axes[2].plot(
            result["time"],
            result["theta_x"],
            label=f"theta_x: {mode}",
        )
        axes[2].plot(
            result["time"],
            result["theta_r"],
            linestyle="--",
            label=f"theta_r: {mode}",
        )
    axes[2].set_xlabel("time [s]")
    axes[2].set_ylabel("adaptive parameters")
    axes[2].grid(True)
    axes[2].legend(ncol=3)

    fig.tight_layout()
    fig.savefig(output_path, dpi=180)
    plt.close(fig)


def main() -> None:
    cfg = Configuration()
    results = {
        mode: simulate(mode, cfg)
        for mode in ("naive", "gated", "compensated")
    }

    print("Mode          RMS(e)   RMS recovery   max||theta||   saturation")
    print("-" * 70)
    for mode, result in results.items():
        print(
            f"{mode:12s}"
            f"{result['rms_error']:9.4f}"
            f"{result['recovery_rms_error']:15.4f}"
            f"{result['maximum_parameter_norm']:15.4f}"
            f"{result['saturation_fraction']:13.3f}"
        )

    output_dir = Path(__file__).resolve().parent
    save_csv(results, output_dir / "Chapter19_Lesson2_results.csv")
    plot_results(results, output_dir / "Chapter19_Lesson2_results.png")


if __name__ == "__main__":
    main()

The C++17 program uses only the standard library and exports a CSV file.

Chapter19_Lesson2.cpp

/*
Chapter19_Lesson2.cpp
Saturation-aware adaptive update laws for a scalar MRAC example.

Build:
    g++ -std=c++17 -O2 -Wall -Wextra -pedantic Chapter19_Lesson2.cpp -o Chapter19_Lesson2
Run:
    ./Chapter19_Lesson2
*/

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

enum class Mode { Naive, Gated, Compensated };

struct Configuration {
    double a{-0.4};
    double b{1.0};
    double am{-1.2};
    double bm{1.2};
    double umax{1.0};
    double gammaX{4.0};
    double gammaR{4.0};
    double sigmaSat{0.15};
    double thetaMin{-5.0};
    double thetaMax{5.0};
    double dt{1.0e-3};
    double finalTime{24.0};

    [[nodiscard]] double p() const {
        return -1.0 / (2.0 * am);
    }
};

struct Sample {
    double t{};
    double r{};
    double x{};
    double xm{};
    double e{};
    double ec{};
    double eDelta{};
    double uc{};
    double ua{};
    double deltaU{};
    double thetaX{};
    double thetaR{};
    double rho{};
};

struct Result {
    std::vector<Sample> samples;
    double rmsError{};
    double recoveryRmsError{};
    double maximumParameterNorm{};
    double saturationFraction{};
};

double referenceSignal(const double t) {
    if (t < 1.0) {
        return 0.0;
    }
    if (t < 8.0) {
        return 4.0;
    }
    if (t < 16.0) {
        return 0.6;
    }
    return 0.6 + 0.25 * std::sin(0.8 * (t - 16.0));
}

double saturate(const double value, const double limit) {
    return std::clamp(value, -limit, limit);
}

double projectDirection(
    const double theta,
    const double direction,
    const double lower,
    const double upper
) {
    if (theta <= lower && direction < 0.0) {
        return 0.0;
    }
    if (theta >= upper && direction > 0.0) {
        return 0.0;
    }
    return direction;
}

std::string modeName(const Mode mode) {
    switch (mode) {
        case Mode::Naive:
            return "naive";
        case Mode::Gated:
            return "gated";
        case Mode::Compensated:
            return "compensated";
    }
    throw std::logic_error("Unhandled mode.");
}

Result simulate(const Mode mode, const Configuration& 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;
    double eDelta = 0.0;
    double thetaX = 0.0;
    double thetaR = 0.0;

    double sumErrorSquared = 0.0;
    double recoveryErrorSquared = 0.0;
    std::size_t recoveryCount = 0U;
    std::size_t saturationCount = 0U;
    double maximumParameterNorm = 0.0;

    for (std::size_t k = 0; k < steps; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        const double r = referenceSignal(t);
        const double uc = thetaX * x + thetaR * r;
        const double ua = saturate(uc, cfg.umax);
        const double deltaU = ua - uc;
        const bool saturated = std::abs(deltaU) > 1.0e-12;

        const double e = x - xm;
        const double ec = e - eDelta;

        double adaptationError = e;
        double rho = 1.0;
        double leakage = 0.0;

        if (mode == Mode::Gated) {
            rho = saturated ? 0.0 : 1.0;
        } else if (mode == Mode::Compensated) {
            adaptationError = ec;
            leakage = saturated ? cfg.sigmaSat : 0.0;
        }

        result.samples.push_back(
            Sample{
                t, r, x, xm, e, ec, eDelta, uc, ua, deltaU,
                thetaX, thetaR, rho
            }
        );

        sumErrorSquared += e * e;
        if (t >= 8.0 && t <= 14.0) {
            recoveryErrorSquared += e * e;
            ++recoveryCount;
        }
        if (saturated) {
            ++saturationCount;
        }
        maximumParameterNorm = std::max(
            maximumParameterNorm,
            std::hypot(thetaX, thetaR)
        );

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

        const double rawThetaXDot =
            -rho * cfg.gammaX * x * cfg.p() * cfg.b * adaptationError
            -leakage * cfg.gammaX * thetaX;
        const double rawThetaRDot =
            -rho * cfg.gammaR * r * cfg.p() * cfg.b * adaptationError
            -leakage * cfg.gammaR * thetaR;

        const double thetaXDot = projectDirection(
            thetaX, rawThetaXDot, cfg.thetaMin, cfg.thetaMax
        );
        const double thetaRDot = projectDirection(
            thetaR, rawThetaRDot, cfg.thetaMin, cfg.thetaMax
        );

        const double xDot = cfg.a * x + cfg.b * ua;
        const double xmDot = cfg.am * xm + cfg.bm * r;
        const double eDeltaDot = cfg.am * eDelta + cfg.b * deltaU;

        x += cfg.dt * xDot;
        xm += cfg.dt * xmDot;
        eDelta += cfg.dt * eDeltaDot;
        thetaX = std::clamp(
            thetaX + cfg.dt * thetaXDot,
            cfg.thetaMin,
            cfg.thetaMax
        );
        thetaR = std::clamp(
            thetaR + cfg.dt * thetaRDot,
            cfg.thetaMin,
            cfg.thetaMax
        );
    }

    result.rmsError = std::sqrt(
        sumErrorSquared / static_cast<double>(result.samples.size())
    );
    result.recoveryRmsError = std::sqrt(
        recoveryErrorSquared / static_cast<double>(recoveryCount)
    );
    result.maximumParameterNorm = maximumParameterNorm;
    result.saturationFraction =
        static_cast<double>(saturationCount)
        / static_cast<double>(result.samples.size());

    return result;
}

void writeCsv(
    const std::array<Result, 3>& results,
    const std::array<Mode, 3>& modes,
    const std::string& fileName
) {
    std::ofstream output(fileName);
    if (!output) {
        throw std::runtime_error("Cannot open output CSV file.");
    }

    output << "time";
    for (const Mode mode : modes) {
        const std::string name = modeName(mode);
        output
            << ',' << name << "_reference"
            << ',' << name << "_x"
            << ',' << name << "_xm"
            << ',' << name << "_e"
            << ',' << name << "_ec"
            << ',' << name << "_uc"
            << ',' << name << "_ua"
            << ',' << name << "_delta_u"
            << ',' << name << "_theta_x"
            << ',' << name << "_theta_r"
            << ',' << name << "_rho";
    }
    output << '\n';

    const std::size_t count = results[0].samples.size();
    output << std::setprecision(12);
    for (std::size_t k = 0; k < count; ++k) {
        output << results[0].samples[k].t;
        for (const Result& result : results) {
            const Sample& s = result.samples[k];
            output
                << ',' << s.r
                << ',' << s.x
                << ',' << s.xm
                << ',' << s.e
                << ',' << s.ec
                << ',' << s.uc
                << ',' << s.ua
                << ',' << s.deltaU
                << ',' << s.thetaX
                << ',' << s.thetaR
                << ',' << s.rho;
        }
        output << '\n';
    }
}

int main() {
    try {
        const Configuration cfg{};
        const std::array<Mode, 3> modes{
            Mode::Naive,
            Mode::Gated,
            Mode::Compensated
        };

        std::array<Result, 3> results{
            simulate(modes[0], cfg),
            simulate(modes[1], cfg),
            simulate(modes[2], cfg)
        };

        std::cout
            << "Mode          RMS(e)   RMS recovery   max||theta||   saturation\n"
            << std::string(70, '-') << '\n'
            << std::fixed << std::setprecision(4);

        for (std::size_t i = 0; i < modes.size(); ++i) {
            const Result& result = results[i];
            std::cout
                << std::left << std::setw(13) << modeName(modes[i])
                << std::right << std::setw(8) << result.rmsError
                << std::setw(15) << result.recoveryRmsError
                << std::setw(15) << result.maximumParameterNorm
                << std::setw(13) << result.saturationFraction
                << '\n';
        }

        writeCsv(results, modes, "Chapter19_Lesson2_cpp_results.csv");
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

The Java implementation uses records for immutable simulation data.

Chapter19_Lesson2.java

/*
Chapter19_Lesson2.java
Saturation-aware adaptive update laws for a scalar MRAC example.

Build:
    javac Chapter19_Lesson2.java
Run:
    java Chapter19_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 Chapter19_Lesson2 {
    private enum Mode {
        NAIVE("naive"),
        GATED("gated"),
        COMPENSATED("compensated");

        private final String label;

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

        String label() {
            return label;
        }
    }

    private record Configuration(
        double a,
        double b,
        double am,
        double bm,
        double uMax,
        double gammaX,
        double gammaR,
        double sigmaSat,
        double thetaMin,
        double thetaMax,
        double dt,
        double finalTime
    ) {
        static Configuration defaults() {
            return new Configuration(
                -0.4, 1.0, -1.2, 1.2, 1.0,
                4.0, 4.0, 0.15, -5.0, 5.0,
                1.0e-3, 24.0
            );
        }

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

    private record Sample(
        double time,
        double reference,
        double x,
        double xm,
        double error,
        double compensatedError,
        double saturationErrorState,
        double commandedInput,
        double actuatedInput,
        double inputMismatch,
        double thetaX,
        double thetaR,
        double rho
    ) {}

    private record Result(
        List<Sample> samples,
        double rmsError,
        double recoveryRmsError,
        double maximumParameterNorm,
        double saturationFraction
    ) {}

    private Chapter19_Lesson2() {}

    private static double referenceSignal(double time) {
        if (time < 1.0) {
            return 0.0;
        }
        if (time < 8.0) {
            return 4.0;
        }
        if (time < 16.0) {
            return 0.6;
        }
        return 0.6 + 0.25 * Math.sin(0.8 * (time - 16.0));
    }

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

    private static double projectDirection(
        double theta,
        double direction,
        double lower,
        double upper
    ) {
        if (theta <= lower && direction < 0.0) {
            return 0.0;
        }
        if (theta >= upper && direction > 0.0) {
            return 0.0;
        }
        return direction;
    }

    private static Result simulate(Mode mode, Configuration 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 eDelta = 0.0;
        double thetaX = 0.0;
        double thetaR = 0.0;

        double sumErrorSquared = 0.0;
        double recoveryErrorSquared = 0.0;
        int recoveryCount = 0;
        int saturationCount = 0;
        double maximumParameterNorm = 0.0;

        for (int k = 0; k < steps; ++k) {
            double time = k * cfg.dt();
            double reference = referenceSignal(time);
            double commandedInput = thetaX * x + thetaR * reference;
            double actuatedInput = clamp(
                commandedInput,
                -cfg.uMax(),
                cfg.uMax()
            );
            double inputMismatch = actuatedInput - commandedInput;
            boolean saturated = Math.abs(inputMismatch) > 1.0e-12;

            double error = x - xm;
            double compensatedError = error - eDelta;

            double adaptationError = error;
            double rho = 1.0;
            double leakage = 0.0;

            switch (mode) {
                case NAIVE -> {
                    // Nominal update law.
                }
                case GATED -> rho = saturated ? 0.0 : 1.0;
                case COMPENSATED -> {
                    adaptationError = compensatedError;
                    leakage = saturated ? cfg.sigmaSat() : 0.0;
                }
            }

            samples.add(
                new Sample(
                    time,
                    reference,
                    x,
                    xm,
                    error,
                    compensatedError,
                    eDelta,
                    commandedInput,
                    actuatedInput,
                    inputMismatch,
                    thetaX,
                    thetaR,
                    rho
                )
            );

            sumErrorSquared += error * error;
            if (time >= 8.0 && time <= 14.0) {
                recoveryErrorSquared += error * error;
                recoveryCount++;
            }
            if (saturated) {
                saturationCount++;
            }
            maximumParameterNorm = Math.max(
                maximumParameterNorm,
                Math.hypot(thetaX, thetaR)
            );

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

            double rawThetaXDot =
                -rho * cfg.gammaX() * x * cfg.p() * cfg.b() * adaptationError
                -leakage * cfg.gammaX() * thetaX;
            double rawThetaRDot =
                -rho * cfg.gammaR() * reference
                * cfg.p() * cfg.b() * adaptationError
                -leakage * cfg.gammaR() * thetaR;

            double thetaXDot = projectDirection(
                thetaX,
                rawThetaXDot,
                cfg.thetaMin(),
                cfg.thetaMax()
            );
            double thetaRDot = projectDirection(
                thetaR,
                rawThetaRDot,
                cfg.thetaMin(),
                cfg.thetaMax()
            );

            double xDot = cfg.a() * x + cfg.b() * actuatedInput;
            double xmDot = cfg.am() * xm + cfg.bm() * reference;
            double eDeltaDot =
                cfg.am() * eDelta + cfg.b() * inputMismatch;

            x += cfg.dt() * xDot;
            xm += cfg.dt() * xmDot;
            eDelta += cfg.dt() * eDeltaDot;
            thetaX = clamp(
                thetaX + cfg.dt() * thetaXDot,
                cfg.thetaMin(),
                cfg.thetaMax()
            );
            thetaR = clamp(
                thetaR + cfg.dt() * thetaRDot,
                cfg.thetaMin(),
                cfg.thetaMax()
            );
        }

        return new Result(
            samples,
            Math.sqrt(sumErrorSquared / samples.size()),
            Math.sqrt(recoveryErrorSquared / recoveryCount),
            maximumParameterNorm,
            (double) saturationCount / samples.size()
        );
    }

    private static void writeCsv(
        List<Result> results,
        List<Mode> modes,
        Path outputPath
    ) throws IOException {
        try (
            BufferedWriter writer = Files.newBufferedWriter(
                outputPath,
                StandardCharsets.UTF_8
            )
        ) {
            writer.write("time");
            for (Mode mode : modes) {
                String name = mode.label();
                writer.write("," + name + "_reference");
                writer.write("," + name + "_x");
                writer.write("," + name + "_xm");
                writer.write("," + name + "_e");
                writer.write("," + name + "_ec");
                writer.write("," + name + "_uc");
                writer.write("," + name + "_ua");
                writer.write("," + name + "_delta_u");
                writer.write("," + name + "_theta_x");
                writer.write("," + name + "_theta_r");
                writer.write("," + name + "_rho");
            }
            writer.newLine();

            int count = results.get(0).samples().size();
            for (int k = 0; k < count; ++k) {
                writer.write(
                    String.format(
                        Locale.ROOT,
                        "%.12g",
                        results.get(0).samples().get(k).time()
                    )
                );
                for (Result result : results) {
                    Sample s = result.samples().get(k);
                    writer.write(
                        String.format(
                            Locale.ROOT,
                            ",%.12g,%.12g,%.12g,%.12g,%.12g,"
                                + "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g",
                            s.reference(),
                            s.x(),
                            s.xm(),
                            s.error(),
                            s.compensatedError(),
                            s.commandedInput(),
                            s.actuatedInput(),
                            s.inputMismatch(),
                            s.thetaX(),
                            s.thetaR(),
                            s.rho()
                        )
                    );
                }
                writer.newLine();
            }
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Configuration cfg = Configuration.defaults();
        List<Mode> modes = List.of(
            Mode.NAIVE,
            Mode.GATED,
            Mode.COMPENSATED
        );
        List<Result> results = modes.stream()
            .map(mode -> simulate(mode, cfg))
            .toList();

        System.out.println(
            "Mode          RMS(e)   RMS recovery   max||theta||   saturation"
        );
        System.out.println("-".repeat(70));
        for (int i = 0; i < modes.size(); ++i) {
            Result result = results.get(i);
            System.out.printf(
                Locale.ROOT,
                "%-13s%8.4f%15.4f%15.4f%13.4f%n",
                modes.get(i).label(),
                result.rmsError(),
                result.recoveryRmsError(),
                result.maximumParameterNorm(),
                result.saturationFraction()
            );
        }

        try {
            writeCsv(
                results,
                modes,
                Path.of("Chapter19_Lesson2_java_results.csv")
            );
        } catch (IOException error) {
            System.err.println("Could not write CSV: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script performs the three-mode comparison. Its final local function programmatically builds a fixed-step Simulink model containing the compensated adaptive loop.

Chapter19_Lesson2.m

%% Chapter19_Lesson2.m
% Saturation-aware adaptive update laws for a scalar MRAC example.
%
% The script compares:
%   1) naive       - nominal adaptation driven by raw tracking error,
%   2) gated       - adaptation frozen while the actuator saturates,
%   3) compensated - auxiliary saturation-error dynamics, projection,
%                    and saturation-activated leakage.
%
% Run:
%   Chapter19_Lesson2
%
% Optional Simulink model:
%   buildChapter19Lesson2Simulink

clear; close all; clc;

cfg = defaultConfiguration();
modes = ["naive", "gated", "compensated"];
results = struct();

for i = 1:numel(modes)
    mode = modes(i);
    results.(mode) = simulateMode(mode, cfg);
end

fprintf('Mode          RMS(e)   RMS recovery   max||theta||   saturation\n');
fprintf('%s\n', repmat('-', 1, 70));
for i = 1:numel(modes)
    mode = modes(i);
    item = results.(mode);
    fprintf('%-13s%8.4f%15.4f%15.4f%13.4f\n', ...
        mode, ...
        item.rmsError, ...
        item.recoveryRmsError, ...
        item.maximumParameterNorm, ...
        item.saturationFraction);
end

writeResultsCsv(results, modes, "Chapter19_Lesson2_matlab_results.csv");
plotResults(results);

%% Local functions

function cfg = defaultConfiguration()
    cfg.a = -0.4;
    cfg.b = 1.0;
    cfg.am = -1.2;
    cfg.bm = 1.2;
    cfg.uMax = 1.0;
    cfg.gamma = diag([4.0, 4.0]);
    cfg.sigmaSat = 0.15;
    cfg.thetaMin = -5.0;
    cfg.thetaMax = 5.0;
    cfg.dt = 1.0e-3;
    cfg.finalTime = 24.0;
    cfg.P = -1.0 / (2.0 * cfg.am);  % 2*am*P = -1
end

function value = referenceSignal(t)
    if t < 1.0
        value = 0.0;
    elseif t < 8.0
        value = 4.0;
    elseif t < 16.0
        value = 0.6;
    else
        value = 0.6 + 0.25 * sin(0.8 * (t - 16.0));
    end
end

function value = saturate(value, limit)
    value = min(limit, max(-limit, value));
end

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

function result = simulateMode(mode, cfg)
    t = (0.0:cfg.dt:cfg.finalTime).';
    n = numel(t);

    x = zeros(n, 1);
    xm = zeros(n, 1);
    eDelta = zeros(n, 1);
    theta = zeros(n, 2);

    r = zeros(n, 1);
    uc = zeros(n, 1);
    ua = zeros(n, 1);
    deltaU = zeros(n, 1);
    e = zeros(n, 1);
    ec = zeros(n, 1);
    rho = ones(n, 1);

    thetaAnchor = [0.0; 0.0];

    for k = 1:n-1
        r(k) = referenceSignal(t(k));
        phi = [x(k); r(k)];

        uc(k) = theta(k, :) * phi;
        ua(k) = saturate(uc(k), cfg.uMax);
        deltaU(k) = ua(k) - uc(k);

        e(k) = x(k) - xm(k);
        ec(k) = e(k) - eDelta(k);
        saturated = abs(deltaU(k)) > 1.0e-12;

        switch mode
            case "naive"
                adaptationError = e(k);
                leakage = 0.0;
            case "gated"
                adaptationError = e(k);
                rho(k) = double(~saturated);
                leakage = 0.0;
            case "compensated"
                adaptationError = ec(k);
                leakage = cfg.sigmaSat * double(saturated);
            otherwise
                error("Unsupported mode: %s", mode);
        end

        nominalDirection = ...
            -rho(k) * cfg.gamma * phi * (cfg.P * cfg.b * adaptationError);
        leakageDirection = ...
            -leakage * cfg.gamma * (theta(k, :).' - thetaAnchor);
        rawThetaDot = nominalDirection + leakageDirection;
        thetaDot = componentwiseProjection( ...
            theta(k, :).', ...
            rawThetaDot, ...
            cfg.thetaMin, ...
            cfg.thetaMax);

        xDot = cfg.a * x(k) + cfg.b * ua(k);
        xmDot = cfg.am * xm(k) + cfg.bm * r(k);
        eDeltaDot = cfg.am * eDelta(k) + cfg.b * deltaU(k);

        x(k+1) = x(k) + cfg.dt * xDot;
        xm(k+1) = xm(k) + cfg.dt * xmDot;
        eDelta(k+1) = eDelta(k) + cfg.dt * eDeltaDot;
        theta(k+1, :) = theta(k, :) + cfg.dt * thetaDot.';
        theta(k+1, :) = min( ...
            cfg.thetaMax, ...
            max(cfg.thetaMin, theta(k+1, :)));
    end

    r(end) = referenceSignal(t(end));
    phi = [x(end); r(end)];
    uc(end) = theta(end, :) * phi;
    ua(end) = saturate(uc(end), cfg.uMax);
    deltaU(end) = ua(end) - uc(end);
    e(end) = x(end) - xm(end);
    ec(end) = e(end) - eDelta(end);
    rho(end) = rho(end-1);

    recovery = t >= 8.0 & t <= 14.0;
    parameterNorm = vecnorm(theta, 2, 2);

    result.time = t;
    result.reference = r;
    result.x = x;
    result.xm = xm;
    result.e = e;
    result.ec = ec;
    result.eDelta = eDelta;
    result.uc = uc;
    result.ua = ua;
    result.deltaU = deltaU;
    result.thetaX = theta(:, 1);
    result.thetaR = theta(:, 2);
    result.rho = rho;
    result.rmsError = sqrt(mean(e.^2));
    result.recoveryRmsError = sqrt(mean(e(recovery).^2));
    result.maximumParameterNorm = max(parameterNorm);
    result.saturationFraction = mean(abs(deltaU) > 1.0e-9);
end

function writeResultsCsv(results, modes, fileName)
    t = results.(modes(1)).time;
    variableNames = "time";
    matrix = t;

    fields = [ ...
        "reference", "x", "xm", "e", "ec", ...
        "uc", "ua", "deltaU", "thetaX", "thetaR", "rho"];

    for i = 1:numel(modes)
        mode = modes(i);
        item = results.(mode);
        for j = 1:numel(fields)
            field = fields(j);
            matrix = [matrix, item.(field)]; %#ok<AGROW>
            variableNames(end+1) = mode + "_" + field; %#ok<AGROW>
        end
    end

    tableData = array2table(matrix, ...
        'VariableNames', cellstr(variableNames));
    writetable(tableData, fileName);
end

function plotResults(results)
    modes = ["naive", "gated", "compensated"];

    figure('Name', 'Chapter 19 Lesson 2', 'Color', 'w');
    tiledlayout(3, 1);

    nexttile;
    hold on;
    for mode = modes
        plot(results.(mode).time, results.(mode).x, ...
            'DisplayName', "x: " + mode);
    end
    plot(results.compensated.time, results.compensated.xm, ...
        'k--', 'LineWidth', 1.3, 'DisplayName', 'reference model');
    ylabel('state');
    grid on;
    legend('Location', 'best');

    nexttile;
    hold on;
    for mode = modes
        plot(results.(mode).time, results.(mode).ua, ...
            'DisplayName', "u_a: " + mode);
    end
    yline(1.0, 'k:');
    yline(-1.0, 'k:');
    ylabel('actuated input');
    grid on;
    legend('Location', 'best');

    nexttile;
    hold on;
    for mode = modes
        plot(results.(mode).time, results.(mode).thetaX, ...
            'DisplayName', "theta_x: " + mode);
        plot(results.(mode).time, results.(mode).thetaR, '--', ...
            'DisplayName', "theta_r: " + mode);
    end
    xlabel('time [s]');
    ylabel('adaptive parameters');
    grid on;
    legend('Location', 'best');

    exportgraphics(gcf, "Chapter19_Lesson2_matlab_results.png", ...
        'Resolution', 180);
end

function buildChapter19Lesson2Simulink()
%BUILDCHAPTER19LESSON2SIMULINK Build a discrete educational Simulink model.
%
% The MATLAB Function block contains the compensated update law. The model
% uses a fixed-step discrete solver with sample time 1 ms.

    model = "Chapter19_Lesson2_Simulink";
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    if isfile(model + ".slx")
        delete(model + ".slx");
    end

    new_system(model);
    open_system(model);

    add_block("simulink/Sources/Clock", model + "/Clock", ...
        'Position', [40 80 70 100]);
    add_block("simulink/User-Defined Functions/MATLAB Function", ...
        model + "/AdaptiveLoop", ...
        'Position', [150 45 350 145]);
    add_block("simulink/Signal Routing/Mux", model + "/Mux", ...
        'Inputs', '6', ...
        'Position', [410 42 415 168]);
    add_block("simulink/Sinks/Scope", model + "/Scope", ...
        'Position', [480 75 530 125]);
    add_block("simulink/Sinks/To Workspace", model + "/LoggedData", ...
        'VariableName', 'Chapter19Lesson2SimulinkData', ...
        'SaveFormat', 'Structure With Time', ...
        'Position', [480 145 590 175]);

    add_line(model, "Clock/1", "AdaptiveLoop/1");
    for port = 1:6
        add_line(model, "AdaptiveLoop/" + port, "Mux/" + port);
    end
    add_line(model, "Mux/1", "Scope/1");
    add_line(model, "Mux/1", "LoggedData/1");

    rt = sfroot;
    chart = rt.find('-isa', 'Stateflow.EMChart', ...
        'Path', model + "/AdaptiveLoop");

    chart.Script = sprintf([ ...
        'function [x,xm,uc,ua,thetaX,thetaR] = fcn(t)\n' ...
        '%%#codegen\n' ...
        'persistent xp xmp ed thx thr\n' ...
        'if isempty(xp)\n' ...
        '    xp=0; xmp=0; ed=0; thx=0; thr=0;\n' ...
        'end\n' ...
        'dt=1e-3; a=-0.4; b=1; am=-1.2; bm=1.2;\n' ...
        'umax=1; P=-1/(2*am); gx=4; gr=4; sigma=0.15;\n' ...
        'if t<1\n' ...
        '    r=0;\n' ...
        'elseif t<8\n' ...
        '    r=4;\n' ...
        'elseif t<16\n' ...
        '    r=0.6;\n' ...
        'else\n' ...
        '    r=0.6+0.25*sin(0.8*(t-16));\n' ...
        'end\n' ...
        'uc=thx*xp+thr*r;\n' ...
        'ua=min(umax,max(-umax,uc));\n' ...
        'du=ua-uc; e=xp-xmp; ec=e-ed;\n' ...
        'sat=double(abs(du)>1e-12);\n' ...
        'dthx=-gx*xp*P*b*ec-sigma*sat*gx*thx;\n' ...
        'dthr=-gr*r*P*b*ec-sigma*sat*gr*thr;\n' ...
        'if thx<=-5 && dthx<0, dthx=0; end\n' ...
        'if thx>= 5 && dthx>0, dthx=0; end\n' ...
        'if thr<=-5 && dthr<0, dthr=0; end\n' ...
        'if thr>= 5 && dthr>0, dthr=0; end\n' ...
        'xp=xp+dt*(a*xp+b*ua);\n' ...
        'xmp=xmp+dt*(am*xmp+bm*r);\n' ...
        'ed=ed+dt*(am*ed+b*du);\n' ...
        'thx=min(5,max(-5,thx+dt*dthx));\n' ...
        'thr=min(5,max(-5,thr+dt*dthr));\n' ...
        'x=xp; xm=xmp; thetaX=thx; thetaR=thr;\n' ...
        'end\n']);

    set_param(model, ...
        'SolverType', 'Fixed-step', ...
        'Solver', 'FixedStepDiscrete', ...
        'FixedStep', '0.001', ...
        'StopTime', '24');

    save_system(model);
    fprintf('Created %s.slx\n', model);
end

The downloadable notebook contains this Wolfram Language source in an executable input cell. A Chapter19_Lesson2.wl copy is included for command-line use.

Chapter19_Lesson2.nb

(* Chapter19_Lesson2.nb
   Saturation-aware adaptive update laws for a scalar MRAC example. *)

ClearAll["Global`*"];

cfg = <|
   "a" -> -0.4, "b" -> 1.0,
   "am" -> -1.2, "bm" -> 1.2,
   "uMax" -> 1.0,
   "gamma" -> DiagonalMatrix[{4.0, 4.0}],
   "sigmaSat" -> 0.15,
   "thetaMin" -> -5.0, "thetaMax" -> 5.0,
   "dt" -> 0.001, "finalTime" -> 24.0
   |>;
cfg["P"] = -1.0/(2.0 cfg["am"]);

referenceSignal[t_] := Piecewise[
   {
    {0.0, t < 1.0},
    {4.0, t < 8.0},
    {0.6, t < 16.0}
    },
   0.6 + 0.25 Sin[0.8 (t - 16.0)]
   ];

saturate[value_, limit_] := Clip[value, {-limit, limit}];

componentwiseProjection[theta_, direction_, lower_, upper_] :=
  MapThread[
   Function[{th, dir},
    Which[
     th <= lower && dir < 0.0, 0.0,
     th >= upper && dir > 0.0, 0.0,
     True, dir
     ]
    ],
   {theta, direction}
   ];

simulate[mode_String] := Module[
  {
   dt = cfg["dt"], steps, time, x, xm, eDelta, theta,
   reference, uc, ua, deltaU, error, compensatedError, rho,
   k, phi, saturated, adaptationError, leakage,
   nominalDirection, leakageDirection, rawThetaDot, thetaDot,
   xDot, xmDot, eDeltaDot, recoveryMask, parameterNorm
   },

  steps = Round[cfg["finalTime"]/dt] + 1;
  time = N@Range[0, steps - 1] dt;

  x = ConstantArray[0.0, steps];
  xm = ConstantArray[0.0, steps];
  eDelta = ConstantArray[0.0, steps];
  theta = ConstantArray[{0.0, 0.0}, steps];

  reference = ConstantArray[0.0, steps];
  uc = ConstantArray[0.0, steps];
  ua = ConstantArray[0.0, steps];
  deltaU = ConstantArray[0.0, steps];
  error = ConstantArray[0.0, steps];
  compensatedError = ConstantArray[0.0, steps];
  rho = ConstantArray[1.0, steps];

  Do[
   reference[[k]] = referenceSignal[time[[k]]];
   phi = {x[[k]], reference[[k]]};

   uc[[k]] = theta[[k]].phi;
   ua[[k]] = saturate[uc[[k]], cfg["uMax"]];
   deltaU[[k]] = ua[[k]] - uc[[k]];

   error[[k]] = x[[k]] - xm[[k]];
   compensatedError[[k]] = error[[k]] - eDelta[[k]];
   saturated = Abs[deltaU[[k]]] > 10^-12;

   Switch[mode,
    "naive",
    adaptationError = error[[k]]; leakage = 0.0,

    "gated",
    adaptationError = error[[k]];
    rho[[k]] = If[saturated, 0.0, 1.0];
    leakage = 0.0,

    "compensated",
    adaptationError = compensatedError[[k]];
    leakage = If[saturated, cfg["sigmaSat"], 0.0],

    _,
    Print["Unsupported mode: ", mode]; Abort[]
    ];

   nominalDirection =
    -rho[[k]] (cfg["gamma"].phi)
      (cfg["P"] cfg["b"] adaptationError);
   leakageDirection =
    -leakage cfg["gamma"].theta[[k]];
   rawThetaDot = nominalDirection + leakageDirection;

   thetaDot = componentwiseProjection[
     theta[[k]], rawThetaDot,
     cfg["thetaMin"], cfg["thetaMax"]
     ];

   xDot = cfg["a"] x[[k]] + cfg["b"] ua[[k]];
   xmDot = cfg["am"] xm[[k]] + cfg["bm"] reference[[k]];
   eDeltaDot = cfg["am"] eDelta[[k]] + cfg["b"] deltaU[[k]];

   x[[k + 1]] = x[[k]] + dt xDot;
   xm[[k + 1]] = xm[[k]] + dt xmDot;
   eDelta[[k + 1]] = eDelta[[k]] + dt eDeltaDot;
   theta[[k + 1]] = Clip[
     theta[[k]] + dt thetaDot,
     {cfg["thetaMin"], cfg["thetaMax"]}
     ],
   {k, 1, steps - 1}
   ];

  reference[[-1]] = referenceSignal[time[[-1]]];
  phi = {x[[-1]], reference[[-1]]};
  uc[[-1]] = theta[[-1]].phi;
  ua[[-1]] = saturate[uc[[-1]], cfg["uMax"]];
  deltaU[[-1]] = ua[[-1]] - uc[[-1]];
  error[[-1]] = x[[-1]] - xm[[-1]];
  compensatedError[[-1]] = error[[-1]] - eDelta[[-1]];
  rho[[-1]] = rho[[-2]];

  recoveryMask = Map[8.0 <= # <= 14.0 &, time];
  parameterNorm = Norm /@ theta;

  <|
   "time" -> time,
   "reference" -> reference,
   "x" -> x,
   "xm" -> xm,
   "e" -> error,
   "ec" -> compensatedError,
   "eDelta" -> eDelta,
   "uc" -> uc,
   "ua" -> ua,
   "deltaU" -> deltaU,
   "thetaX" -> theta[[All, 1]],
   "thetaR" -> theta[[All, 2]],
   "rho" -> rho,
   "rmsError" -> Sqrt[Mean[error^2]],
   "recoveryRmsError" ->
    Sqrt[Mean[Pick[error^2, recoveryMask]]],
   "maximumParameterNorm" -> Max[parameterNorm],
   "saturationFraction" ->
    Mean[Boole /@ Thread[Abs[deltaU] > 10^-9]]
   |>
  ];

modes = {"naive", "gated", "compensated"};
results = AssociationMap[simulate, modes];

Grid[
 Prepend[
  Table[
   {
    mode,
    NumberForm[results[mode]["rmsError"], {7, 4}],
    NumberForm[results[mode]["recoveryRmsError"], {7, 4}],
    NumberForm[results[mode]["maximumParameterNorm"], {7, 4}],
    NumberForm[results[mode]["saturationFraction"], {7, 4}]
    },
   {mode, modes}
   ],
  {"Mode", "RMS(e)", "RMS recovery", "max norm(theta)", "saturation"}
  ],
 Frame -> All
 ]

fields = {
   "reference", "x", "xm", "e", "ec",
   "uc", "ua", "deltaU", "thetaX", "thetaR", "rho"
   };

header = Join[
   {"time"},
   Flatten[
    Table[mode <> "_" <> field, {mode, modes}, {field, fields}]
    ]
   ];

rows = Transpose@Join[
    {results["naive"]["time"]},
    Flatten[
     Table[results[mode][field], {mode, modes}, {field, fields}],
     1
     ]
    ];

Export["Chapter19_Lesson2_mathematica_results.csv",
 Prepend[rows, header]];

trackingPlot = Show[
   Table[
    ListLinePlot[
     Transpose[{results[mode]["time"], results[mode]["x"]}],
     PlotLegends -> {("x: " <> mode)}
     ],
    {mode, modes}
    ],
   ListLinePlot[
    Transpose[{
      results["compensated"]["time"],
      results["compensated"]["xm"]
      }],
    PlotStyle -> Directive[Black, Dashed],
    PlotLegends -> {"reference model"}
    ],
   Frame -> True,
   FrameLabel -> {"time [s]", "state"},
   PlotRange -> All,
   ImageSize -> Large
   ];

inputPlot = ListLinePlot[
   Table[
    Transpose[{results[mode]["time"], results[mode]["ua"]}],
    {mode, modes}
    ],
   PlotLegends -> modes,
   Frame -> True,
   FrameLabel -> {"time [s]", "actuated input"},
   PlotRange -> All,
   ImageSize -> Large
   ];

parameterPlot = ListLinePlot[
   Flatten[
    Table[
     {
      Transpose[{results[mode]["time"], results[mode]["thetaX"]}],
      Transpose[{results[mode]["time"], results[mode]["thetaR"]}]
      },
     {mode, modes}
     ],
    1
    ],
   PlotLegends -> Flatten[
     Table[{"thetaX: " <> mode, "thetaR: " <> mode}, {mode, modes}]
     ],
   Frame -> True,
   FrameLabel -> {"time [s]", "adaptive parameters"},
   PlotRange -> All,
   ImageSize -> Large
   ];

GraphicsColumn[{trackingPlot, inputPlot, parameterPlot},
 ImageSize -> Large]

13.1 Recommended tuning sequence

1. Set physical amplitude and rate limits from actuator data.
2. Verify the unsaturated baseline loop before adaptation.
3. Tune \(\boldsymbol\Gamma\) without saturation.
4. Validate the measured mismatch \(\Delta u\).
5. Add the auxiliary model using the stable reference dynamics.
6. Apply physically defensible parameter bounds.
7. Add only enough leakage to limit prolonged drift.
8. Test feasible, marginal, and infeasible commands.

13.2 Required verification plots

Report \(x\) and \(x_m\), \(u_c\) and \(u_a\), \(\Delta u\), \(e\) and \(e_c\), parameter estimates, saturation duty cycle, and recovery time after the command becomes feasible.

13.3 Safety interpretation

Modified adaptation protects the estimator; it does not create actuator authority. If the reference requires unavailable effort, the supervisory layer must limit or reshape the command, switch modes, or redesign the reference model.

Problem 1: For \(\dot e=a_me+ b\widetilde{\boldsymbol\theta}^T\boldsymbol\phi+ b\Delta u\), derive \(\dot V\) under the nominal gradient law.

Solution:

\[ \begin{aligned} \dot V &= 2Pe\dot e+ 2\widetilde{\boldsymbol\theta}^T \boldsymbol\Gamma^{-1} \dot{\hat{\boldsymbol\theta}}\\ &= 2Pa_me^2+ 2Pbe\widetilde{\boldsymbol\theta}^T\boldsymbol\phi+ 2Pbe\Delta u- 2Pbe\widetilde{\boldsymbol\theta}^T\boldsymbol\phi\\ &=-Qe^2+2Pbe\Delta u. \end{aligned} \]

The product \(e\Delta u\) may be positive, so bounded mismatch does not imply \(\dot V\le0\).

Problem 2: Construct an auxiliary state that removes \(b\Delta u\) and prove the compensated result.

Solution:

Choose \(\dot e_\Delta=a_me_\Delta+b\Delta u\) and \(e_c=e-e_\Delta\). Then \(\dot e_c=a_me_c+ b\widetilde{\boldsymbol\theta}^T\boldsymbol\phi\). The update \(\dot{\hat{\boldsymbol\theta}}= -\boldsymbol\Gamma\boldsymbol\phi bPe_c\) cancels the parameter cross term and gives \(\dot V_c=-Qe_c^2\).

Problem 3: Give the componentwise tangent-cone projection for \(\hat\theta_i\in[ \underline\theta_i,\overline\theta_i]\).

Solution:

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

Problem 4: Show the residual term introduced by leakage toward \(\boldsymbol\theta_0\).

Solution:

\[ \hat{\boldsymbol\theta}-\boldsymbol\theta_0 = \widetilde{\boldsymbol\theta}+ (\boldsymbol\theta^{*}-\boldsymbol\theta_0). \]

Applying Young's inequality to the leakage contribution gives

\[ \dot V_c \le -Qe_c^2- \sigma_s\|\widetilde{\boldsymbol\theta}\|^2+ \sigma_s\|\boldsymbol\theta^{*}-\boldsymbol\theta_0\|^2. \]

Problem 5: Let \(T_s=0.01\), \(a_m=-2\), \(b=1\), \(e_\Delta[k]=0.2\), and \(\Delta u[k]=-0.5\). Compute the next auxiliary state.

Solution:

\[ \begin{aligned} e_\Delta[k+1] &= 0.2+0.01\left((-2)(0.2)+(1)(-0.5)\right)\\ &=0.2-0.009=0.191. \end{aligned} \]

Saturation inserts \(b\Delta u\) into the tracking error and destroys the nominal Lyapunov decrease. Gating can stop the estimator from chasing an infeasible command, while a stable auxiliary error model provides a stronger structural remedy: the compensated error obeys the same parameter-affine dynamics as the unsaturated design. Projection enforces parameter limits and saturation-activated leakage limits prolonged drift at the cost of bias.