In ideal MRAC, the parameter update law interprets every nonzero tracking error as evidence of parameter mismatch. In a real loop, however, the measured error also contains sensor noise, bounded disturbances, discretization effects, and residual unmodeled dynamics. When the true tracking error is already small, continued integration of these residual signals may move the parameter estimate even though no useful model information is present. This is the parameter-drift mechanism introduced in the preceding lessons.

Two complementary design ideas are used here:

• Dead-zone modification: adapt only when the magnitude of the error exceeds a selected threshold.
• Leakage modification: add a dissipative term that prevents estimates from behaving like unconstrained integrators.

Dead-zone adaptation produces a set-valued objective: the controller is required to drive the error into a prescribed neighborhood of zero, not necessarily to the exact origin. Leakage produces an explicit restoring action in parameter space, but generally introduces steady-state bias.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> C["Adaptive controller"]
  P["Plant with uncertainty and disturbance"] --> Y["Measured output"]
  C --> P
  RM --> E["Tracking error"]
  Y --> E
  E --> DZ["Dead-zone error processor"]
  DZ --> U["Gradient update"]
  L["Leakage feedback from \nparameter estimate"] --> U
  U --> TH["Parameter estimate"]
  TH --> C
        

To isolate the robust-update mechanism, consider the normalized scalar tracking-error model obtained after applying the matching conditions from Chapters 5–7:

\[ \dot e = -a_m e + \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} + d(t), \qquad a_m > 0, \qquad |d(t)| \le \bar d. \]

Here \(e=y-y_m\), the regressor \(\boldsymbol{\phi}\in\mathbb{R}^{p}\) contains measured plant/reference signals, and \(\widetilde{\boldsymbol{\theta} }=\widehat{\boldsymbol{\theta} }- \boldsymbol{\theta}^{*}\) is the parameter-estimation error. The constant vector \(\boldsymbol{\theta}^{*}\) denotes an ideal matching parameter vector for the nominal plant.

For the first-order plant and reference model

\[ \dot y=-a y+b u+d(t), \qquad \dot y_m=-a_m y_m+b_m r, \]

with controller

\[ u=\widehat\theta_y y+\widehat\theta_r r, \]

the matching parameters are

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

After normalization by the known control-gain sign and magnitude, the nominal gradient update is

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

With the Lyapunov function

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

the parameter cross-term cancels and

\[ \dot V=-a_m e^2+e d(t). \]

The disturbance term has indefinite sign. More importantly, when \(e\) is a small noise-dominated signal, the nominal update law still integrates it indefinitely.

A hard dead-zone gate can be written as

\[ g_{\delta}(e)= \begin{cases} 0, & |e|\le\delta,\\ 1, & |e|>\delta, \end{cases} \qquad \dot{\widehat{\boldsymbol{\theta} } } =-g_{\delta}(e)\boldsymbol{\Gamma}\boldsymbol{\phi}e. \]

The hard gate is intuitive, but the update changes discontinuously at \(|e|=\delta\). A continuous alternative uses the dead-zone error

\[ e_{\delta}=\mathcal{D}_{\delta}(e)= \begin{cases} 0, & |e|\le\delta,\\ e-\delta\operatorname{sgn}(e), & |e|>\delta, \end{cases} \]

and the update law

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_{\delta} } \]

The quantity \(e_{\delta}\) is the signed distance from \(e\) to the closed interval \(\mathcal{E}_{\delta}=[-\delta,\delta]\). Therefore, \(e_{\delta}=0\) does not imply exact tracking; it means that tracking lies inside the selected robustness tube.

\[ e=e_{\delta}+\delta\operatorname{sgn}(e), \qquad e e_{\delta}=e_{\delta}^{2}+\delta|e_{\delta}| \quad \text{whenever } |e|>\delta. \]

Theorem 1 (Convergence to a dead-zone set).

Assume the normalized error model of Section 2, bounded disturbance \(|d(t)|\le\bar d\), bounded regressor, and bounded derivatives of the closed-loop signals. Apply \(\dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_{\delta}\). If

\[ \boxed{\delta>\frac{\bar d}{a_m} }, \]

then the parameter estimate remains bounded and the distance of the tracking error from \(\mathcal{E}_{\delta}\) converges to zero.

Proof. Use the set-distance Lyapunov function

\[ V_{\delta}=\frac{1}{2}e_{\delta}^{2}+ \frac{1}{2}\widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1}\widetilde{\boldsymbol{\theta} }. \]

Inside the dead zone, \(e_{\delta}=0\) and the parameter update is zero. Outside the dead zone, \(\dot e_{\delta}=\dot e\) almost everywhere. Hence

\[ \begin{aligned} \dot V_{\delta} &=e_{\delta}\left(-a_m e+ \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}+d\right) +\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \left(-\boldsymbol{\Gamma}\boldsymbol{\phi}e_{\delta}\right)\\ &=-a_m e e_{\delta}+d e_{\delta}\\ &=-a_m e_{\delta}^{2}-a_m\delta|e_{\delta}|+d e_{\delta}\\ &\le -a_m e_{\delta}^{2}-\left(a_m\delta-\bar d\right)|e_{\delta}|. \end{aligned} \]

The threshold condition makes the right-hand side nonpositive. Thus \(V_{\delta}\) and \(\widetilde{\boldsymbol{\theta} }\) are bounded, while \(e_{\delta}\in L_2\cap L_{\infty}\). Under the stated bounded-derivative assumption, Barbalat's lemma gives \(e_{\delta}(t)\to0\). Therefore the distance from \(e(t)\) to the dead-zone interval converges to zero. \(\square\)

The theorem does not claim convergence of the estimate to the ideal parameter. Once the error enters the dead zone, learning stops. The final estimate is therefore history-dependent, which is intentional: robustness is obtained by refusing to learn from low-quality residual error.

Leakage augments the gradient law with negative feedback from the parameter estimate:

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e -\sigma\boldsymbol{\Gamma}\widehat{\boldsymbol{\theta} } }, \qquad \sigma>0. \]

The first term learns from tracking error; the second makes the estimator a leaky integrator. Using the nominal Lyapunov function from Section 2,

\[ \begin{aligned} \dot V &=-a_m e^2+ed- \sigma\widetilde{\boldsymbol{\theta} }^{T} \widehat{\boldsymbol{\theta} }\\ &=-a_m e^2+ed- \sigma\|\widetilde{\boldsymbol{\theta} }\|^2- \sigma\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\theta}^{*}. \end{aligned} \]

Young's inequality gives

\[ ed\le\frac{a_m}{2}e^2+\frac{\bar d^2}{2a_m}, \qquad -\sigma\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\theta}^{*} \le\frac{\sigma}{2}\|\widetilde{\boldsymbol{\theta} }\|^2+ \frac{\sigma}{2}\|\boldsymbol{\theta}^{*}\|^2. \]

Therefore

\[ \boxed{ \dot V\le-\frac{a_m}{2}e^2- \frac{\sigma}{2}\|\widetilde{\boldsymbol{\theta} }\|^2+c_0}, \qquad c_0=\frac{\bar d^2}{2a_m}+ \frac{\sigma}{2}\|\boldsymbol{\theta}^{*}\|^2. \]

Thus the state and parameter errors are uniformly ultimately bounded: the derivative is negative outside the ellipsoid

\[ a_m e^2+\sigma\|\widetilde{\boldsymbol{\theta} }\|^2 \le 2c_0. \]

Leakage prevents unbounded drift, but the term involving \(\boldsymbol{\theta}^{*}\) reveals the bias. Even with no disturbance, leakage pulls \(\widehat{\boldsymbol{\theta} }\) toward zero rather than toward an unknown nonzero ideal vector. Increasing \(\sigma\) strengthens drift suppression while usually enlarging tracking bias.

A useful implementation combines the continuous dead-zone error with leakage that is active only outside the dead zone:

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_{\delta} -\sigma g_{\delta}(e)\widehat{\boldsymbol{\theta} } }. \]

When \(|e|\le\delta\), both terms vanish and the estimate freezes. When the error leaves the tube, gradient correction and leakage act together. The corresponding derivative outside the tube satisfies

\[ \dot V_{\delta}\le -a_m e_{\delta}^{2}-(a_m\delta-\bar d)|e_{\delta}| -\sigma\widetilde{\boldsymbol{\theta} }^{T} \widehat{\boldsymbol{\theta} }. \]

The final leakage term can be bounded exactly as in Section 5, yielding an ultimate-bound result. Compared with always-on leakage, gating avoids continued decay of a satisfactory parameter estimate after the tracking error has entered the accepted tube.

6.1 Tuning the Dead-Zone Width

The theorem suggests \(\delta>\bar d/a_m\) for the normalized error model. In practice, the effective error floor includes sensor noise and unmodeled dynamics, so a data-based rule is

\[ \delta=\kappa e_{\mathrm{floor} }, \qquad 1.2\le\kappa\le2.0, \]

where \(e_{\mathrm{floor} }\) is estimated from a nonadaptive or low-gain closed-loop experiment. A threshold that is too small fails to stop noise-driven adaptation; one that is too large permits unnecessarily large tracking error.

6.2 Tuning Leakage

The leakage time scale is approximately \(1/\sigma\) when the gradient term is weak. Choose it slower than the desired tracking transient but fast enough to arrest long-term drift. Start with a small value, inspect the parameter norm and steady tracking error, and increase only when drift remains unacceptable.

All five implementations simulate the same plant, reference model, deterministic disturbance, and deterministic measurement noise:

\[ \dot y=-1.4y+u+d(t), \qquad \dot y_m=-2y_m+2r, \qquad \boldsymbol{\phi}=\begin{bmatrix}y_{\mathrm{meas} }&r\end{bmatrix}^{T}. \]

\[ d(t)=0.08\sin(7t)+0.04\operatorname{sgn}(\sin(0.7t)), \]

\[ n(t)=0.02\sin(37t)+0.01\sin(73t), \qquad y_{\mathrm{meas} }=y+n. \]

The numerical values are \(\gamma=8\), \(\delta=0.08\), \(\sigma=0.08\), and \(\Delta t=10^{-3}\,\mathrm{s}\). The reference is multisine for the first ten seconds, constant for the next ten seconds, and zero afterward. This final interval makes low-level noise-driven adaptation visible.

In the validated run, Python, C++, and Java produced identical metrics to displayed precision. The plain law achieved the smallest RMS error in this particular benchmark, while leakage strongly reduced the final parameter norm. That outcome is expected: robust modifications are not designed to dominate nominal adaptation in every finite simulation; they trade some tracking accuracy for stronger drift control and boundedness margins.

The Python version uses NumPy for vectorized storage and Matplotlib for comparison plots. It exports one CSV file per update law.

Chapter8_Lesson3.py

"""Chapter 8, Lesson 3: Dead-Zone and Leakage Modifications.

Compares four normalized first-order MRAC update laws:
plain, dead-zone, leakage, and dead-zone + gated leakage.
Requires: numpy, matplotlib.
"""
from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 1.0e-3
    tf: float = 30.0
    a: float = 1.4
    b: float = 1.0
    am: float = 2.0
    bm: float = 2.0
    gamma: float = 8.0
    dead_zone: float = 0.08
    leakage: float = 0.08


def reference(t: float) -> float:
    if t < 10.0:
        return 0.8 * np.sin(0.7 * t) + 0.3 * np.sin(1.9 * t)
    if t < 20.0:
        return 0.5
    return 0.0


def disturbance(t: float) -> float:
    square = 1.0 if np.sin(0.7 * t) >= 0.0 else -1.0
    return 0.08 * np.sin(7.0 * t) + 0.04 * square


def measurement_noise(t: float) -> float:
    return 0.02 * np.sin(37.0 * t) + 0.01 * np.sin(73.0 * t)


def dead_zone_error(e: float, delta: float) -> float:
    if abs(e) <= delta:
        return 0.0
    return e - np.sign(e) * delta


def simulate(mode: str, cfg: Config) -> dict[str, np.ndarray]:
    valid_modes = {"plain", "deadzone", "leakage", "combined"}
    if mode not in valid_modes:
        raise ValueError(f"mode must be one of {sorted(valid_modes)}")

    n = int(round(cfg.tf / cfg.dt)) + 1
    t = np.linspace(0.0, cfg.tf, n)
    y = np.zeros(n)
    ym = np.zeros(n)
    u = np.zeros(n)
    r_hist = np.zeros(n)
    e_hist = np.zeros(n)
    theta = np.zeros((n, 2))

    for k in range(n - 1):
        tk = t[k]
        r = reference(tk)
        y_meas = y[k] + measurement_noise(tk)
        e = y_meas - ym[k]
        phi = np.array([y_meas, r])

        if mode in {"deadzone", "combined"}:
            adaptive_error = dead_zone_error(e, cfg.dead_zone)
            gate = 1.0 if abs(e) > cfg.dead_zone else 0.0
        else:
            adaptive_error = e
            gate = 1.0

        leak = cfg.gamma * cfg.leakage * theta[k] if mode in {"leakage", "combined"} else 0.0
        if mode == "combined":
            leak = gate * cfg.gamma * cfg.leakage * theta[k]

        theta_dot = -cfg.gamma * phi * adaptive_error - leak
        control = theta[k, 0] * y_meas + theta[k, 1] * r
        y_dot = -cfg.a * y[k] + cfg.b * control + disturbance(tk)
        ym_dot = -cfg.am * ym[k] + cfg.bm * r

        y[k + 1] = y[k] + cfg.dt * y_dot
        ym[k + 1] = ym[k] + cfg.dt * ym_dot
        theta[k + 1] = theta[k] + cfg.dt * theta_dot
        u[k] = control
        r_hist[k] = r
        e_hist[k] = e

    r_hist[-1] = reference(t[-1])
    e_hist[-1] = y[-1] + measurement_noise(t[-1]) - ym[-1]
    u[-1] = theta[-1, 0] * (y[-1] + measurement_noise(t[-1])) + theta[-1, 1] * r_hist[-1]

    return {
        "t": t,
        "r": r_hist,
        "y": y,
        "ym": ym,
        "e": e_hist,
        "u": u,
        "theta_y": theta[:, 0],
        "theta_r": theta[:, 1],
    }


def metrics(result: dict[str, np.ndarray], cfg: Config) -> tuple[float, float, float]:
    e = result["e"]
    tail = result["t"] >= 20.0
    rms = float(np.sqrt(np.mean(e**2)))
    tail_rms = float(np.sqrt(np.mean(e[tail] ** 2)))
    theta_norm = float(np.hypot(result["theta_y"][-1], result["theta_r"][-1]))
    return rms, tail_rms, theta_norm


def main() -> None:
    cfg = Config()
    modes = ["plain", "deadzone", "leakage", "combined"]
    results = {mode: simulate(mode, cfg) for mode in modes}

    output_dir = Path(__file__).resolve().parent
    print("mode,rms_error,tail_rms_error,final_parameter_norm")
    for mode, result in results.items():
        rms, tail_rms, theta_norm = metrics(result, cfg)
        print(f"{mode},{rms:.6f},{tail_rms:.6f},{theta_norm:.6f}")
        matrix = np.column_stack([result[key] for key in ("t", "r", "y", "ym", "e", "u", "theta_y", "theta_r")])
        np.savetxt(
            output_dir / f"Chapter8_Lesson3_{mode}.csv",
            matrix,
            delimiter=",",
            header="t,r,y,ym,e,u,theta_y,theta_r",
            comments="",
        )

    plt.figure(figsize=(10, 5))
    for mode, result in results.items():
        plt.plot(result["t"], result["e"], label=mode)
    plt.axhline(cfg.dead_zone, linestyle="--", linewidth=1)
    plt.axhline(-cfg.dead_zone, linestyle="--", linewidth=1)
    plt.xlabel("Time (s)")
    plt.ylabel("Measured tracking error")
    plt.title("MRAC tracking error: robust update-law comparison")
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter8_Lesson3_error_comparison.png", dpi=180)

    plt.figure(figsize=(10, 5))
    for mode, result in results.items():
        theta_norm = np.hypot(result["theta_y"], result["theta_r"])
        plt.plot(result["t"], theta_norm, label=mode)
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter-estimate norm")
    plt.title("Parameter drift and leakage/dead-zone suppression")
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter8_Lesson3_parameter_norm.png", dpi=180)

    plt.show()


if __name__ == "__main__":
    main()

The C++17 program uses only the standard library, writes CSV outputs, and reports RMS and parameter-norm metrics. A typical build command is g++ -std=c++17 -O2 Chapter8_Lesson3.cpp -o Chapter8_Lesson3.

Chapter8_Lesson3.cpp

// Chapter 8, Lesson 3: Dead-Zone and Leakage Modifications.
// Standard-library-only normalized first-order MRAC simulation.
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

struct Config {
    double dt = 1.0e-3;
    double tf = 30.0;
    double a = 1.4;
    double b = 1.0;
    double am = 2.0;
    double bm = 2.0;
    double gamma = 8.0;
    double deadZone = 0.08;
    double leakage = 0.08;
};

struct Sample {
    double t{}, r{}, y{}, ym{}, e{}, u{}, thetaY{}, thetaR{};
};

double reference(double t) {
    if (t < 10.0) return 0.8 * std::sin(0.7 * t) + 0.3 * std::sin(1.9 * t);
    if (t < 20.0) return 0.5;
    return 0.0;
}

double disturbance(double t) {
    const double square = std::sin(0.7 * t) >= 0.0 ? 1.0 : -1.0;
    return 0.08 * std::sin(7.0 * t) + 0.04 * square;
}

double measurementNoise(double t) {
    return 0.02 * std::sin(37.0 * t) + 0.01 * std::sin(73.0 * t);
}

double deadZoneError(double e, double delta) {
    if (std::abs(e) <= delta) return 0.0;
    return e - std::copysign(delta, e);
}

std::vector<Sample> simulate(const std::string& mode, const Config& cfg) {
    if (mode != "plain" && mode != "deadzone" && mode != "leakage" && mode != "combined") {
        throw std::invalid_argument("Unknown mode: " + mode);
    }

    const std::size_t n = static_cast<std::size_t>(std::llround(cfg.tf / cfg.dt)) + 1U;
    std::vector<Sample> out(n);
    double y = 0.0;
    double ym = 0.0;
    std::array<double, 2> theta{0.0, 0.0};

    for (std::size_t k = 0; k < n; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        const double r = reference(t);
        const double yMeas = y + measurementNoise(t);
        const double e = yMeas - ym;
        const std::array<double, 2> phi{yMeas, r};
        const bool hasDeadZone = mode == "deadzone" || mode == "combined";
        const bool hasLeakage = mode == "leakage" || mode == "combined";
        const double adaptiveError = hasDeadZone ? deadZoneError(e, cfg.deadZone) : e;
        const double gate = (!hasDeadZone || std::abs(e) > cfg.deadZone) ? 1.0 : 0.0;
        const double control = theta[0] * yMeas + theta[1] * r;

        out[k] = Sample{t, r, y, ym, e, control, theta[0], theta[1]};
        if (k + 1U == n) break;

        std::array<double, 2> thetaDot{};
        for (std::size_t i = 0; i < 2U; ++i) {
            double leak = hasLeakage ? cfg.gamma * cfg.leakage * theta[i] : 0.0;
            if (mode == "combined") leak *= gate;
            thetaDot[i] = -cfg.gamma * phi[i] * adaptiveError - leak;
        }

        const double yDot = -cfg.a * y + cfg.b * control + disturbance(t);
        const double ymDot = -cfg.am * ym + cfg.bm * r;
        y += cfg.dt * yDot;
        ym += cfg.dt * ymDot;
        theta[0] += cfg.dt * thetaDot[0];
        theta[1] += cfg.dt * thetaDot[1];
    }
    return out;
}

void writeCsv(const std::string& filename, const std::vector<Sample>& data) {
    std::ofstream file(filename);
    if (!file) throw std::runtime_error("Cannot open output file: " + filename);
    file << "t,r,y,ym,e,u,theta_y,theta_r\n" << std::setprecision(12);
    for (const auto& s : data) {
        file << s.t << ',' << s.r << ',' << s.y << ',' << s.ym << ',' << s.e << ','
             << s.u << ',' << s.thetaY << ',' << s.thetaR << '\n';
    }
}

int main() {
    try {
        const Config cfg;
        const std::array<std::string, 4> modes{"plain", "deadzone", "leakage", "combined"};
        std::cout << "mode,rms_error,tail_rms_error,final_parameter_norm\n";
        for (const auto& mode : modes) {
            const auto data = simulate(mode, cfg);
            double sumSq = 0.0;
            double tailSumSq = 0.0;
            std::size_t tailCount = 0U;
            for (const auto& s : data) {
                sumSq += s.e * s.e;
                if (s.t >= 20.0) {
                    tailSumSq += s.e * s.e;
                    ++tailCount;
                }
            }
            const double rms = std::sqrt(sumSq / static_cast<double>(data.size()));
            const double tailRms = std::sqrt(tailSumSq / static_cast<double>(tailCount));
            const auto& last = data.back();
            const double thetaNorm = std::hypot(last.thetaY, last.thetaR);
            std::cout << mode << ',' << rms << ',' << tailRms << ',' << thetaNorm << '\n';
            writeCsv("Chapter8_Lesson3_" + mode + ".csv", data);
        }
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses records for simulation samples and standard NIO file output. Compile with javac Chapter8_Lesson3.java and run with java Chapter8_Lesson3.

Chapter8_Lesson3.java

// Chapter 8, Lesson 3: Dead-Zone and Leakage Modifications.
// Standard-library-only normalized first-order MRAC simulation.
import java.io.BufferedWriter;
import java.io.IOException;
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 Chapter8_Lesson3 {
    private static final class Config {
        final double dt = 1.0e-3;
        final double tf = 30.0;
        final double a = 1.4;
        final double b = 1.0;
        final double am = 2.0;
        final double bm = 2.0;
        final double gamma = 8.0;
        final double deadZone = 0.08;
        final double leakage = 0.08;
    }

    private record Sample(double t, double r, double y, double ym, double e,
                          double u, double thetaY, double thetaR) {}

    private static double reference(double t) {
        if (t < 10.0) return 0.8 * Math.sin(0.7 * t) + 0.3 * Math.sin(1.9 * t);
        if (t < 20.0) return 0.5;
        return 0.0;
    }

    private static double disturbance(double t) {
        double square = Math.sin(0.7 * t) >= 0.0 ? 1.0 : -1.0;
        return 0.08 * Math.sin(7.0 * t) + 0.04 * square;
    }

    private static double measurementNoise(double t) {
        return 0.02 * Math.sin(37.0 * t) + 0.01 * Math.sin(73.0 * t);
    }

    private static double deadZoneError(double e, double delta) {
        if (Math.abs(e) <= delta) return 0.0;
        return e - Math.copySign(delta, e);
    }

    private static List<Sample> simulate(String mode, Config cfg) {
        if (!List.of("plain", "deadzone", "leakage", "combined").contains(mode)) {
            throw new IllegalArgumentException("Unknown mode: " + mode);
        }

        int n = (int) Math.round(cfg.tf / cfg.dt) + 1;
        List<Sample> output = new ArrayList<>(n);
        double y = 0.0;
        double ym = 0.0;
        double[] theta = {0.0, 0.0};

        for (int k = 0; k < n; k++) {
            double t = k * cfg.dt;
            double r = reference(t);
            double yMeas = y + measurementNoise(t);
            double e = yMeas - ym;
            double[] phi = {yMeas, r};
            boolean hasDeadZone = mode.equals("deadzone") || mode.equals("combined");
            boolean hasLeakage = mode.equals("leakage") || mode.equals("combined");
            double adaptiveError = hasDeadZone ? deadZoneError(e, cfg.deadZone) : e;
            double gate = (!hasDeadZone || Math.abs(e) > cfg.deadZone) ? 1.0 : 0.0;
            double control = theta[0] * yMeas + theta[1] * r;

            output.add(new Sample(t, r, y, ym, e, control, theta[0], theta[1]));
            if (k == n - 1) break;

            double[] thetaDot = new double[2];
            for (int i = 0; i < 2; i++) {
                double leak = hasLeakage ? cfg.gamma * cfg.leakage * theta[i] : 0.0;
                if (mode.equals("combined")) leak *= gate;
                thetaDot[i] = -cfg.gamma * phi[i] * adaptiveError - leak;
            }

            double yDot = -cfg.a * y + cfg.b * control + disturbance(t);
            double ymDot = -cfg.am * ym + cfg.bm * r;
            y += cfg.dt * yDot;
            ym += cfg.dt * ymDot;
            theta[0] += cfg.dt * thetaDot[0];
            theta[1] += cfg.dt * thetaDot[1];
        }
        return output;
    }

    private static void writeCsv(Path filename, List<Sample> data) throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(filename)) {
            writer.write("t,r,y,ym,e,u,theta_y,theta_r\n");
            for (Sample s : data) {
                writer.write(String.format(Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    s.t(), s.r(), s.y(), s.ym(), s.e(), s.u(), s.thetaY(), s.thetaR()));
            }
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Config cfg = new Config();
        String[] modes = {"plain", "deadzone", "leakage", "combined"};
        System.out.println("mode,rms_error,tail_rms_error,final_parameter_norm");

        for (String mode : modes) {
            List<Sample> data = simulate(mode, cfg);
            double sumSq = 0.0;
            double tailSumSq = 0.0;
            int tailCount = 0;
            for (Sample s : data) {
                sumSq += s.e() * s.e();
                if (s.t() >= 20.0) {
                    tailSumSq += s.e() * s.e();
                    tailCount++;
                }
            }
            double rms = Math.sqrt(sumSq / data.size());
            double tailRms = Math.sqrt(tailSumSq / tailCount);
            Sample last = data.get(data.size() - 1);
            double thetaNorm = Math.hypot(last.thetaY(), last.thetaR());
            System.out.printf(Locale.US, "%s,%.6f,%.6f,%.6f%n", mode, rms, tailRms, thetaNorm);
            try {
                writeCsv(Path.of("Chapter8_Lesson3_" + mode + ".csv"), data);
            } catch (IOException ex) {
                System.err.println("Could not write CSV: " + ex.getMessage());
                System.exit(1);
            }
        }
    }
}

The MATLAB file is toolbox-free and includes a direct block-level mapping for Simulink. In Simulink, use separate Integrator blocks for the plant, reference model, and two parameter estimates; implement the dead-zone function in a MATLAB Function block; and form the controller with Product and Sum blocks.

Chapter8_Lesson3.m

% Chapter 8, Lesson 3: Dead-Zone and Leakage Modifications
% Toolbox-free normalized first-order MRAC simulation.
clear; clc; close all;

cfg.dt = 1.0e-3;
cfg.tf = 30.0;
cfg.a = 1.4;
cfg.b = 1.0;
cfg.am = 2.0;
cfg.bm = 2.0;
cfg.gamma = 8.0;
cfg.deadZone = 0.08;
cfg.leakage = 0.08;

modes = {'plain', 'deadzone', 'leakage', 'combined'};
results = struct();
fprintf('mode,rms_error,tail_rms_error,final_parameter_norm\n');

for i = 1:numel(modes)
    mode = modes{i};
    result = simulateMRAC(mode, cfg);
    results.(mode) = result;

    rmsError = sqrt(mean(result.e.^2));
    tail = result.t >= 20.0;
    tailRms = sqrt(mean(result.e(tail).^2));
    thetaNorm = hypot(result.thetaY(end), result.thetaR(end));
    fprintf('%s,%.6f,%.6f,%.6f\n', mode, rmsError, tailRms, thetaNorm);

    T = table(result.t, result.r, result.y, result.ym, result.e, result.u, ...
        result.thetaY, result.thetaR, 'VariableNames', ...
        {'t','r','y','ym','e','u','theta_y','theta_r'});
    writetable(T, ['Chapter8_Lesson3_' mode '.csv']);
end

figure('Name', 'Tracking-error comparison'); hold on;
for i = 1:numel(modes)
    mode = modes{i};
    plot(results.(mode).t, results.(mode).e, 'DisplayName', mode);
end
yline(cfg.deadZone, '--'); yline(-cfg.deadZone, '--');
xlabel('Time (s)'); ylabel('Measured tracking error');
title('MRAC tracking error: robust update-law comparison');
legend('Location','best'); grid on;

figure('Name', 'Parameter-norm comparison'); hold on;
for i = 1:numel(modes)
    mode = modes{i};
    thetaNorm = hypot(results.(mode).thetaY, results.(mode).thetaR);
    plot(results.(mode).t, thetaNorm, 'DisplayName', mode);
end
xlabel('Time (s)'); ylabel('Parameter-estimate norm');
title('Parameter drift and leakage/dead-zone suppression');
legend('Location','best'); grid on;

% Simulink mapping:
% 1. Plant Integrator: y_dot = -a*y + b*u + d.
% 2. Reference Integrator: ym_dot = -am*ym + bm*r.
% 3. Error block: e = y_measured - ym.
% 4. MATLAB Function block: edz = sign(e)*max(abs(e)-delta,0).
% 5. Product/Gain blocks: theta_dot = -gamma*(phi*eadapt + sigma*theta) when leakage is active.
% 6. Two Integrator blocks hold thetaY and thetaR; u = [y_measured r]*theta.

function result = simulateMRAC(mode, cfg)
    validModes = {'plain', 'deadzone', 'leakage', 'combined'};
    assert(any(strcmp(mode, validModes)), 'Unknown mode.');

    t = (0:cfg.dt:cfg.tf)';
    n = numel(t);
    y = zeros(n,1); ym = zeros(n,1); e = zeros(n,1); u = zeros(n,1);
    rHist = zeros(n,1); theta = zeros(n,2);

    for k = 1:n-1
        tk = t(k);
        r = referenceSignal(tk);
        yMeasured = y(k) + measurementNoise(tk);
        err = yMeasured - ym(k);
        phi = [yMeasured; r];

        hasDeadZone = strcmp(mode,'deadzone') || strcmp(mode,'combined');
        hasLeakage = strcmp(mode,'leakage') || strcmp(mode,'combined');
        if hasDeadZone
            adaptiveError = deadZoneError(err, cfg.deadZone);
            gate = double(abs(err) > cfg.deadZone);
        else
            adaptiveError = err;
            gate = 1.0;
        end

        leak = zeros(2,1);
        if hasLeakage
            leak = cfg.gamma * cfg.leakage * theta(k,:)';
        end
        if strcmp(mode,'combined')
            leak = gate * leak;
        end

        thetaDot = -cfg.gamma * phi * adaptiveError - leak;
        control = theta(k,1) * yMeasured + theta(k,2) * r;
        yDot = -cfg.a * y(k) + cfg.b * control + disturbanceSignal(tk);
        ymDot = -cfg.am * ym(k) + cfg.bm * r;

        y(k+1) = y(k) + cfg.dt * yDot;
        ym(k+1) = ym(k) + cfg.dt * ymDot;
        theta(k+1,:) = theta(k,:) + cfg.dt * thetaDot';
        e(k) = err; u(k) = control; rHist(k) = r;
    end

    rHist(end) = referenceSignal(t(end));
    e(end) = y(end) + measurementNoise(t(end)) - ym(end);
    u(end) = theta(end,1) * (y(end) + measurementNoise(t(end))) + theta(end,2) * rHist(end);

    result = struct('t',t,'r',rHist,'y',y,'ym',ym,'e',e,'u',u, ...
        'thetaY',theta(:,1),'thetaR',theta(:,2));
end

function r = referenceSignal(t)
    if t < 10.0
        r = 0.8*sin(0.7*t) + 0.3*sin(1.9*t);
    elseif t < 20.0
        r = 0.5;
    else
        r = 0.0;
    end
end

function d = disturbanceSignal(t)
    square = 1.0;
    if sin(0.7*t) < 0.0, square = -1.0; end
    d = 0.08*sin(7.0*t) + 0.04*square;
end

function n = measurementNoise(t)
    n = 0.02*sin(37.0*t) + 0.01*sin(73.0*t);
end

function edz = deadZoneError(e, delta)
    if abs(e) <= delta
        edz = 0.0;
    else
        edz = e - sign(e)*delta;
    end
end

The notebook uses built-in Wolfram Language list operations, Association, Export, and ListLinePlot. No external package is required.

Chapter8_Lesson3.nb

Notebook[{
  Cell["Chapter 8, Lesson 3: Dead-Zone and Leakage Modifications", "Title"],
  Cell["Normalized first-order MRAC comparison implemented with built-in Wolfram Language functions.", "Text"],
  Cell["ClearAll[reference, disturbance, measurementNoise, deadZoneError, simulateMRAC, metrics];
reference[t_] := Piecewise[
  { {0.8 Sin[0.7 t] + 0.3 Sin[1.9 t], t < 10}, {0.5, t < 20} }, 0.0];
disturbance[t_] :=
  0.08 Sin[7.0 t] + 0.04 If[Sin[0.7 t] >= 0, 1.0, -1.0];
measurementNoise[t_] := 0.02 Sin[37.0 t] + 0.01 Sin[73.0 t];
deadZoneError[e_, delta_] :=
  If[Abs[e] <= delta, 0.0, e - Sign[e] delta];", "Input"],
  Cell["simulateMRAC[mode_String, dt_: 0.001, tf_: 30.0] := Module[
  {a = 1.4, b = 1.0, am = 2.0, bm = 2.0, gamma = 8.0,
   delta = 0.08, sigma = 0.08, n, t, y, ym, theta, rows, k, r,
   yMeasured, e, phi, hasDeadZone, hasLeakage, adaptiveError,
   gate, leak, thetaDot, u, yDot, ymDot},
  n = Round[tf/dt] + 1; t = N@Range[0, tf, dt];
  y = ConstantArray[0.0, n]; ym = ConstantArray[0.0, n];
  theta = ConstantArray[{0.0, 0.0}, n];
  rows = ConstantArray[{}, n];
  Do[
    r = reference[t[[k]]];
    yMeasured = y[[k]] + measurementNoise[t[[k]]];
    e = yMeasured - ym[[k]]; phi = {yMeasured, r};
    hasDeadZone = MemberQ[{\"deadzone\", \"combined\"}, mode];
    hasLeakage = MemberQ[{\"leakage\", \"combined\"}, mode];
    adaptiveError = If[hasDeadZone, deadZoneError[e, delta], e];
    gate = If[! hasDeadZone || Abs[e] > delta, 1.0, 0.0];
    leak = If[hasLeakage, gamma sigma theta[[k]], {0.0, 0.0}];
    If[mode == \"combined\", leak = gate leak];
    thetaDot = -gamma phi adaptiveError - leak;
    u = theta[[k]].phi;
    rows[[k]] = {t[[k]], r, y[[k]], ym[[k]], e, u,
                 theta[[k, 1]], theta[[k, 2]]};
    If[k < n,
      yDot = -a y[[k]] + b u + disturbance[t[[k]]];
      ymDot = -am ym[[k]] + bm r;
      y[[k + 1]] = y[[k]] + dt yDot;
      ym[[k + 1]] = ym[[k]] + dt ymDot;
      theta[[k + 1]] = theta[[k]] + dt thetaDot;],
    {k, 1, n}];
  AssociationThread[
    {\"Columns\", \"Data\"},
    { {\"t\", \"r\", \"y\", \"ym\", \"e\", \"u\", \"thetaY\", \"thetaR\"}, rows}]
];", "Input"],
  Cell["modes = {\"plain\", \"deadzone\", \"leakage\", \"combined\"};
results = AssociationMap[simulateMRAC, modes];
Do[
  Export[\"Chapter8_Lesson3_\" <> mode <> \".csv\",
    Prepend[results[mode][\"Data\"], results[mode][\"Columns\"]]],
  {mode, modes}];
errorPlot = ListLinePlot[
  AssociationMap[results[#][\"Data\"][[All, {1, 5}]] &, modes],
  PlotLegends -> modes, Frame -> True,
  FrameLabel -> {\"Time (s)\", \"Measured tracking error\"},
  PlotLabel -> \"MRAC update-law comparison\"];
parameterPlot = ListLinePlot[
  AssociationMap[
    ({ #[[1]], Sqrt[#[[7]]^2 + #[[8]]^2]} & /@ results[#][\"Data\"]) &,
    modes],
  PlotLegends -> modes, Frame -> True,
  FrameLabel -> {\"Time (s)\", \"Parameter-estimate norm\"},
  PlotLabel -> \"Parameter drift suppression\"];
Column[{errorPlot, parameterPlot}]", "Input"]
},
WindowTitle -> "Chapter8_Lesson3",
FrontEndVersion -> "13.0 or later"]
      

14.1 Use the Same Error in Analysis and Code

If the theoretical update uses the true tracking error but the code uses a filtered or noisy error, the cancellation in the Lyapunov derivative is no longer exact. Model the measurement path explicitly and interpret the difference as an additional bounded perturbation.

14.2 Avoid Chattering at the Dead-Zone Boundary

The continuous operator \(e_{\delta}\) is preferable to a hard on/off update because the update magnitude approaches zero continuously at the boundary. In sampled implementations, hysteresis may also be used: activate adaptation at \(|e|\ge\delta_{\mathrm{on} }\) and deactivate it at \(|e|\le\delta_{\mathrm{off} }\), with \(\delta_{\mathrm{on} }>\delta_{\mathrm{off} }\).

14.3 Scale Regressors Before Selecting One Adaptation Gain

If components of \(\boldsymbol{\phi}\) differ greatly in magnitude, a scalar adaptation gain produces very different effective update rates. Use a diagonal \(\boldsymbol{\Gamma}\) or scale the regressor signals.

14.4 Leakage Is Not a Substitute for Physical Parameter Bounds

Leakage discourages large estimates but does not impose hard limits. Parameter projection, introduced in Chapter 9, is the appropriate tool when estimates must remain in a known admissible set.

14.5 Discretization Changes the Continuous-Time Guarantee

Forward Euler is used here for transparent cross-language comparison. Real implementations should verify that the sampling period is small relative to plant, reference-model, and adaptation time scales. A large product of adaptation gain, regressor energy, and sample period can cause oscillatory or unstable numerical parameter updates.

Problem 1 (Matching parameters): For \(\dot y=-1.4y+u+d\) and \(\dot y_m=-2y_m+2r\), determine the ideal controller parameters for \(u=\theta_y y+\theta_r r\) when \(d=0\).

Solution: Matching coefficients gives

\[ -1.4+\theta_y^{*}=-2, \qquad \theta_r^{*}=2. \]

Therefore

\[ \boxed{\theta_y^{*}=-0.6, \qquad \theta_r^{*}=2}. \]

Problem 2 (Minimum dead-zone width): Suppose the normalized error dynamics have \(a_m=2\) and \(|d(t)|\le0.12\). Determine the theoretical lower bound on \(\delta\) from Theorem 1. Is \(\delta=0.08\) sufficient?

Solution:

\[ \frac{\bar d}{a_m}=\frac{0.12}{2}=0.06. \]

The theorem requires \(\delta>0.06\). Hence \(\delta=0.08\) is sufficient and leaves a margin of \(0.02\) in the normalized bound.

Problem 3 (Derive the dead-zone Lyapunov inequality): Starting with \(V_{\delta}\) from Section 4, show that outside the dead zone

\[ \dot V_{\delta}\le-a_m e_{\delta}^{2}- (a_m\delta-\bar d)|e_{\delta}|. \]

Solution: Outside the tube,

\[ e=e_{\delta}+\delta\operatorname{sgn}(e), \qquad \dot e_{\delta}=\dot e. \]

Substitution of the error dynamics and adaptive law cancels the parameter cross-term:

\[ \dot V_{\delta}=-a_m e e_{\delta}+d e_{\delta}. \]

Now use

\[ e e_{\delta}=e_{\delta}^{2}+\delta|e_{\delta}|, \qquad d e_{\delta}\le\bar d|e_{\delta}|, \]

which yields the required result.

Problem 4 (Leakage ultimate set): Let \(a_m=2\), \(\bar d=0.1\), \(\sigma=0.05\), and \(\|\boldsymbol{\theta}^{*}\|=2\). Compute \(c_0\) and write the ellipsoid outside which the derived Lyapunov bound is strictly decreasing.

Solution:

\[ c_0=\frac{0.1^2}{2(2)}+\frac{0.05}{2}(2^2) =0.0025+0.1=0.1025. \]

The non-increasing ultimate set is

\[ 2e^2+0.05\|\widetilde{\boldsymbol{\theta} }\|^2 \le 0.205. \]

Problem 5 (Behavior at zero tracking error): Assume \(e=0\) and \(\widehat{\boldsymbol{\theta} }\ne\mathbf{0}\). Compare the instantaneous parameter updates under the plain, dead-zone, always-on leakage, and gated combined laws.

Solution:

• Plain gradient: \(\dot{\widehat{\boldsymbol{\theta} } }=0\).
• Dead-zone: \(e_{\delta}=0\), so the update is zero.
• Always-on leakage: \(\dot{\widehat{\boldsymbol{\theta} } }=-\sigma \boldsymbol{\Gamma}\widehat{\boldsymbol{\theta} }\); the estimate decays toward zero.
• Gated combined law: \(g_{\delta}(0)=0\), so the update is zero and the estimate freezes.

This comparison exposes the main distinction between always-on and gated leakage: only the former changes an estimate after satisfactory tracking has been reached.

Problem 6 (Sampled update): Write the forward-Euler form of the combined update law with sample period \(T_s\).

Solution:

\[ \widehat{\boldsymbol{\theta} }_{k+1}= \widehat{\boldsymbol{\theta} }_{k}-T_s \boldsymbol{\Gamma}\boldsymbol{\phi}_{k}e_{\delta,k} -T_s\sigma g_{\delta}(e_k) \widehat{\boldsymbol{\theta} }_{k}. \]

Equivalently, the leakage-only part multiplies the estimate by \(1-T_s\sigma\) when the gate is active. A basic numerical requirement is \(0<T_s\sigma<2\) for this isolated scalar decay factor to remain within the unit circle; the complete adaptive loop generally requires a much smaller sample period.

Dead-zone modification replaces exact-error convergence with convergence to an explicitly chosen tracking tube and prevents adaptation from integrating low-level residual error. Its continuous form uses the signed distance \(e_{\delta}\), enabling a direct Lyapunov proof when the threshold dominates the normalized disturbance bound. Leakage adds parameter-space dissipation and yields an ultimate-boundedness result, but introduces bias toward its leakage center. A gated combination acts only when the error is outside the accepted tube, balancing drift suppression against unnecessary parameter decay. The next lesson studies the broader robustness–parameter-convergence trade-off created by these modifications.