After completing this lesson, students should be able to:

• distinguish tracking-error convergence from parameter convergence;
• explain why bounded disturbances and modeling errors may cause parameter drift in an unmodified adaptive law;
• prove uniform ultimate boundedness for a sigma-modified MRAC error system;
• calculate the steady-state tracking and parameter bias introduced by leakage;
• compare sigma-modification, e-modification, and dead-zone adaptation;
• tune a robust modification from noise bounds and acceptable bias.

The lesson assumes the Lyapunov-based MRAC constructions of Chapters 5–7 and the robust modifications introduced in Lessons 2 and 3 of this chapter. A formal excitation condition for exact parameter convergence is intentionally deferred to Chapter 10.

Consider the scalar tracking-error model

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

where \(e=y-y_m\) is the tracking error, \(\boldsymbol{\phi}(t)\) is the measurable regressor, \(b\) has known positive sign, and

\[ \widetilde{\boldsymbol{\theta}} = \widehat{\boldsymbol{\theta}} - \boldsymbol{\theta}^{*} \]

is the parameter estimation error. The constant vector \(\boldsymbol{\theta}^{*}\) denotes ideal controller parameters under the matching assumptions developed in Chapter 2. The signal \(d(t)\) collects bounded disturbances and the part of the unmodeled dynamics represented at the error-model level.

The plain gradient/Lyapunov update law is

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

In the ideal case \(d(t)=0\), use

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

Differentiation and cancellation of the cross term give

\[ \begin{aligned} \dot V &= e\dot e + \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1} \dot{\widehat{\boldsymbol{\theta}}} \\ &= -a_m e^2 + b e \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\phi} - b \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\phi}e \\ &= -a_m e^2 \le 0. \end{aligned} \]

This proves boundedness and supports convergence of the tracking error under the usual regularity assumptions. It does not, by itself, prove \(\widehat{\boldsymbol{\theta}} \longrightarrow \boldsymbol{\theta}^{*}\) . Exact parameter recovery requires additional information in the regressor. Chapter 10 formalizes that requirement.

An adaptive controller uses nonzero error as information. Unfortunately, the error may be caused by two fundamentally different mechanisms:

1. a genuine controller-parameter mismatch that adaptation should remove;
2. noise, delay, disturbances, or unmodeled dynamics that cannot be removed by changing the nominal parameters.

A plain adaptive law cannot perfectly distinguish these mechanisms. Therefore it may continue moving the parameters in response to residual errors that are not informative about \(\boldsymbol{\theta}^{*}\). Robust modifications deliberately reduce this sensitivity.

flowchart TD
  A["Measured tracking error"] --> B{"What generated \nthe error?"}
  B --> C["True parameter mismatch"]
  B --> D["Noise, delay, disturbance, \nor unmodeled dynamics"]
  C --> E["Fast adaptation is useful"]
  D --> F["Fast adaptation may create drift \nor high-gain behavior"]
  E --> G["Favor convergence and \nlow nominal bias"]
  F --> H["Add leakage, error scaling, \nor a dead zone"]
  G --> I["Design compromise"]
  H --> I
  I --> J["Bounded parameters with \nacceptable tracking residual"]
        

The central design conflict is therefore:

\[ \text{stronger suppression of parameter motion} \Longleftrightarrow \text{greater bias or slower parameter convergence}. \]

In practical robust adaptive control, the goal is often not exact identification. The goal is bounded internal signals and acceptable tracking despite uncertainty.

The sigma-modified law adds constant leakage:

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

Assume the ideal parameters are constant and \(|d(t)|\le\bar d\). Using the same Lyapunov function,

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

Young's inequality gives

\[ e d \le \frac{a_m}{2}e^2 + \frac{\bar d^2}{2a_m}, \]

\[ -\sigma \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\theta}^{*} \le \frac{\sigma}{2} \left\| \widetilde{\boldsymbol{\theta}} \right\|^2 + \frac{\sigma}{2} \left\| \boldsymbol{\theta}^{*} \right\|^2. \]

Consequently,

\[ \dot V \le -\frac{a_m}{2}e^2 - \frac{\sigma}{2} \left\| \widetilde{\boldsymbol{\theta}} \right\|^2 + \frac{\bar d^2}{2a_m} + \frac{\sigma}{2} \left\| \boldsymbol{\theta}^{*} \right\|^2. \]

Define

\[ \alpha = \min \left\{ a_m,\, \sigma\lambda_{\min}(\boldsymbol{\Gamma}) \right\}, \qquad C = \frac{\bar d^2}{2a_m} + \frac{\sigma}{2} \left\| \boldsymbol{\theta}^{*} \right\|^2. \]

Since \(\|\widetilde{\boldsymbol{\theta}}\|^2 \ge \lambda_{\min}(\boldsymbol{\Gamma}) \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1} \widetilde{\boldsymbol{\theta}}\) , we obtain

\[ \dot V \le -\alpha V+C. \]

The comparison lemma yields

\[ V(t) \le \exp(-\alpha t)V(0) + \frac{C}{\alpha} \left( 1-\exp(-\alpha t) \right), \]

\[ \limsup_{t\longrightarrow\infty}V(t) \le \frac{C}{\alpha}. \]

Thus sigma-modification changes an ideal asymptotic statement into a robust ultimate-boundedness statement. The leakage supplies direct damping in parameter space, but the term involving \(\boldsymbol{\theta}^{*}\) prevents the Lyapunov derivative from being negative definite at the ideal parameter vector.

The trade-off is especially clear for a scalar parameter, constant regressor \(\phi_0\ne0\), and zero disturbance:

\[ \dot e = -a_m e + b\phi_0 \left( \widehat\theta-\theta^{*} \right), \]

\[ \dot{\widehat\theta} = -\gamma b\phi_0 e - \gamma\sigma\widehat\theta. \]

At equilibrium,

\[ 0 = -a_m e_{\mathrm{eq}} + b\phi_0 \left( \widehat\theta_{\mathrm{eq}}-\theta^{*} \right), \]

\[ 0 = -b\phi_0 e_{\mathrm{eq}} - \sigma\widehat\theta_{\mathrm{eq}}. \]

Solving these two equations gives

\[ \boxed{ \widehat\theta_{\mathrm{eq}} = \frac{b^2\phi_0^2} {b^2\phi_0^2+a_m\sigma} \theta^{*} } \]

\[ \boxed{ e_{\mathrm{eq}} = - \frac{\sigma b\phi_0} {b^2\phi_0^2+a_m\sigma} \theta^{*} }. \]

Several important conclusions follow:

• the adaptation gain \(\gamma\) affects transient speed but not this equilibrium;
• any \(\sigma>0\) generally prevents exact convergence to \(\theta^{*}\);
• increasing \(\sigma\) increases parameter damping but also increases the nominal tracking residual;
• stronger regressor magnitude reduces the relative leakage bias.

For positive \(b,\phi_0,\theta^{*}\),

\[ \frac{\partial} {\partial\sigma} \left( \frac{\widehat\theta_{\mathrm{eq}}}{\theta^{*}} \right) = - \frac{ a_m b^2\phi_0^2 }{ \left( b^2\phi_0^2+a_m\sigma \right)^2 } < 0, \]

\[ \frac{\partial |e_{\mathrm{eq}}|} {\partial\sigma} = \frac{ b^3\phi_0^3|\theta^{*}| }{ \left( b^2\phi_0^2+a_m\sigma \right)^2 } > 0. \]

6.1 E-Modification

The e-modification scales leakage by the measured error magnitude:

\[ \dot{\widehat{\boldsymbol{\theta}}} = -\boldsymbol{\Gamma}b\boldsymbol{\phi}e_m - \sigma_e |e_m| \boldsymbol{\Gamma} \widehat{\boldsymbol{\theta}}, \qquad e_m=e+n. \]

When the tracking error is large, strong damping protects the adaptive loop. When the tracking error approaches zero, the leakage approaches zero and the nominal parameter bias can be much smaller than with constant sigma-modification. The cost is that the parameter damping is weakest near the origin, exactly where small persistent disturbances or measurement errors may continue to act.

6.2 Dead-Zone Adaptation

Define the switching factor

\[ D_{\delta}(e_m) = \begin{cases} 0, & |e_m|\le\delta, \\ 1, & |e_m|>\delta. \end{cases} \]

The dead-zone update is

\[ \dot{\widehat{\boldsymbol{\theta}}} = - D_{\delta}(e_m) \boldsymbol{\Gamma}b \boldsymbol{\phi}e_m. \]

Once the error lies inside a region consistent with measurement noise and unavoidable model mismatch, adaptation stops. This prevents noise-driven parameter motion and introduces no continuous pull toward zero. However, the final parameter estimate depends on the point at which adaptation enters the dead zone, so exact parameter convergence is generally lost.

7.1 Tracking Convergence Is Not Identification

A controller can achieve small tracking error while the parameter estimate remains different from the physical or ideal parameter. Closed-loop feedback may compensate for several parameter combinations, disturbances may be absorbed into the estimate, and robust terms may intentionally bias the estimate. Therefore,

\[ e(t)\longrightarrow0 \quad \nRightarrow \quad \widetilde{\boldsymbol{\theta}}(t) \longrightarrow\boldsymbol{0}. \]

7.2 A Practical Tuning Sequence

flowchart TD
  A["Start from the nominal MRAC design"] --> B["Estimate sensor-noise and modeling-error floor"]
  B --> C{"Is a reliable error \nbound available?"}
  C -->|yes| D["Set a dead zone slightly \nabove the bound"]
  C -->|no| E{"Does the estimate drift \nunder stress tests?"}
  E -->|no| F["Keep the plain law or \nuse very mild modification"]
  E -->|yes| G{"Is nominal parameter \nbias critical?"}
  G -->|yes| H["Use e-modification \nor small sigma"]
  G -->|no| I["Use sigma-modification \nfor stronger damping"]
  D --> J["Test ideal and disturbed cases"]
  F --> J
  H --> J
  I --> J
  J --> K["Check tracking RMS, parameter range, and control effort"]
  K --> L["Retune gamma and robustification separately"]
        

Adaptation gain and robustification gain should not be treated as the same tuning knob. The gain \(\boldsymbol{\Gamma}\) determines how rapidly error information moves the estimate. The leakage or dead-zone parameter determines how much of that motion is suppressed.

The downloadable implementations simulate

\[ \dot e = -2e + \phi(t) \left( \widehat\theta-1.5 \right) + d(t), \qquad \phi(t) = 1+0.8\sin(0.7t), \]

with \(\gamma=8\), \(\sigma=0.15\), \(\sigma_e=0.8\), and dead-zone radius \(\delta=0.08\). The disturbed case uses

\[ d(t) = 0.12\sin(5t) + 0.06\operatorname{sgn} \left( \sin(0.35t) \right), \]

\[ n(t)=0.02\sin(37t), \qquad e_m=e+n. \]

The following values are produced by the supplied deterministic implementation using a midpoint integration step of \(10^{-3}\,\mathrm{s}\).

These numerical values are not a universal ranking. They demonstrate that a robustification coefficient selected for strong damping can produce a visible nominal bias. A fair design comparison must test a range of disturbance levels, adaptation gains, and robustification parameters.

8.1 Software Libraries Used in the Implementations

• Python: NumPy stores trajectories and computes metrics; Matplotlib plots tracking and parameter histories. scipy.integrate.solve_ivp is a suitable alternative when adaptive dynamics are embedded in larger models.
• C++: the supplied scalar implementation uses only the C++17 standard library. For vector or MIMO MRAC, Eigen is commonly used for matrix and vector operations.
• Java: the supplied program uses the Java standard library. EJML or Apache Commons Math can support higher-dimensional matrix computations.
• MATLAB/Simulink: base MATLAB performs numerical integration and plotting. Simulink provides block-diagram simulation; Control System Toolbox can be used to construct and analyze the nominal reference model.
• Wolfram Mathematica: NDSolveValue integrates the coupled error and adaptation equations, while symbolic functions can verify equilibrium and Lyapunov calculations.

Chapter8_Lesson4.py

# Chapter8_Lesson4.py
# Trade-off between robustness and parameter convergence in scalar MRAC.
# Requires: numpy, matplotlib

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
import csv
import math

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    a_m: float = 2.0
    b: float = 1.0
    gamma: float = 8.0
    theta_star: float = 1.5
    sigma: float = 0.15
    sigma_e: float = 0.80
    dead_zone: float = 0.08
    dt: float = 0.001
    final_time: float = 40.0


def regressor(t: float) -> float:
    """Bounded time-varying regressor used by every controller."""
    return 1.0 + 0.8 * math.sin(0.7 * t)


def disturbance(t: float, enabled: bool) -> float:
    if not enabled:
        return 0.0
    square = 1.0 if math.sin(0.35 * t) >= 0.0 else -1.0
    return 0.12 * math.sin(5.0 * t) + 0.06 * square


def sensor_noise(t: float, enabled: bool) -> float:
    return 0.02 * math.sin(37.0 * t) if enabled else 0.0


def derivatives(
    t: float,
    error: float,
    theta_hat: float,
    mode: str,
    disturbed: bool,
    cfg: Config,
) -> tuple[float, float]:
    phi = regressor(t)
    measured_error = error + sensor_noise(t, disturbed)

    # Baseline gradient term.
    gradient = -cfg.gamma * cfg.b * phi * measured_error

    if mode == "plain":
        leakage = 0.0
    elif mode == "sigma":
        leakage = cfg.gamma * cfg.sigma * theta_hat
    elif mode == "e":
        leakage = cfg.gamma * cfg.sigma_e * abs(measured_error) * theta_hat
    elif mode == "deadzone":
        leakage = 0.0
        if abs(measured_error) <= cfg.dead_zone:
            gradient = 0.0
    else:
        raise ValueError(f"Unknown mode: {mode}")

    error_dot = (
        -cfg.a_m * error
        + cfg.b * phi * (theta_hat - cfg.theta_star)
        + disturbance(t, disturbed)
    )
    theta_dot = gradient - leakage
    return error_dot, theta_dot


def simulate(mode: str, disturbed: bool, cfg: Config) -> dict[str, np.ndarray]:
    steps = int(round(cfg.final_time / cfg.dt)) + 1
    time = np.linspace(0.0, cfg.final_time, steps)
    error = np.zeros(steps)
    theta_hat = np.zeros(steps)
    error[0] = 1.0
    theta_hat[0] = 0.0

    # Explicit midpoint (RK2): simple, deterministic, and more accurate than Euler.
    for k in range(steps - 1):
        t = float(time[k])
        e = float(error[k])
        th = float(theta_hat[k])

        k1_e, k1_th = derivatives(t, e, th, mode, disturbed, cfg)
        mid_t = t + 0.5 * cfg.dt
        mid_e = e + 0.5 * cfg.dt * k1_e
        mid_th = th + 0.5 * cfg.dt * k1_th
        k2_e, k2_th = derivatives(mid_t, mid_e, mid_th, mode, disturbed, cfg)

        error[k + 1] = e + cfg.dt * k2_e
        theta_hat[k + 1] = th + cfg.dt * k2_th

    return {"time": time, "error": error, "theta_hat": theta_hat}


def metrics(result: dict[str, np.ndarray], cfg: Config) -> dict[str, float]:
    time = result["time"]
    error = result["error"]
    theta_hat = result["theta_hat"]
    tail = time >= (cfg.final_time - 10.0)
    return {
        "rms_error_last_10s": float(np.sqrt(np.mean(error[tail] ** 2))),
        "max_error_last_10s": float(np.max(np.abs(error[tail]))),
        "final_theta_hat": float(theta_hat[-1]),
        "parameter_error": float(theta_hat[-1] - cfg.theta_star),
    }


def save_csv(
    path: Path,
    results: dict[str, dict[str, np.ndarray]],
) -> None:
    modes = list(results)
    time = results[modes[0]]["time"]
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.writer(handle)
        header = ["time"]
        for mode in modes:
            header.extend([f"{mode}_error", f"{mode}_theta_hat"])
        writer.writerow(header)

        for k, t in enumerate(time):
            row: list[float] = [float(t)]
            for mode in modes:
                row.extend(
                    [
                        float(results[mode]["error"][k]),
                        float(results[mode]["theta_hat"][k]),
                    ]
                )
            writer.writerow(row)


def plot_results(
    path: Path,
    results: dict[str, dict[str, np.ndarray]],
    cfg: Config,
    title_suffix: str,
) -> None:
    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

    for mode, result in results.items():
        axes[0].plot(result["time"], result["error"], label=mode)
        axes[1].plot(result["time"], result["theta_hat"], label=mode)

    axes[0].set_ylabel("tracking error e")
    axes[0].set_title(f"MRAC robustness/convergence comparison: {title_suffix}")
    axes[0].grid(True)
    axes[0].legend()

    axes[1].axhline(cfg.theta_star, linestyle="--", label="ideal parameter")
    axes[1].set_xlabel("time (s)")
    axes[1].set_ylabel("parameter estimate")
    axes[1].grid(True)
    axes[1].legend()

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


def run_case(disturbed: bool, cfg: Config) -> None:
    modes = ("plain", "sigma", "e", "deadzone")
    results = {mode: simulate(mode, disturbed, cfg) for mode in modes}
    label = "disturbed" if disturbed else "ideal"

    print(f"\nCase: {label}")
    for mode in modes:
        print(f"{mode:8s}: {metrics(results[mode], cfg)}")

    save_csv(Path(f"Chapter8_Lesson4_{label}.csv"), results)
    plot_results(
        Path(f"Chapter8_Lesson4_{label}.png"),
        results,
        cfg,
        title_suffix=label,
    )


def main() -> None:
    cfg = Config()
    run_case(disturbed=False, cfg=cfg)
    run_case(disturbed=True, cfg=cfg)


if __name__ == "__main__":
    main()

Chapter8_Lesson4.cpp

// Chapter8_Lesson4.cpp
// Trade-off between robustness and parameter convergence in scalar MRAC.
// Build: g++ -std=c++17 -O2 Chapter8_Lesson4.cpp -o Chapter8_Lesson4

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

struct Config {
    double a_m = 2.0;
    double b = 1.0;
    double gamma = 8.0;
    double theta_star = 1.5;
    double sigma = 0.15;
    double sigma_e = 0.80;
    double dead_zone = 0.08;
    double dt = 0.001;
    double final_time = 40.0;
};

struct Result {
    std::vector<double> time;
    std::vector<double> error;
    std::vector<double> theta_hat;
};

double regressor(double t) {
    return 1.0 + 0.8 * std::sin(0.7 * t);
}

double disturbance(double t, bool enabled) {
    if (!enabled) {
        return 0.0;
    }
    const double square = std::sin(0.35 * t) >= 0.0 ? 1.0 : -1.0;
    return 0.12 * std::sin(5.0 * t) + 0.06 * square;
}

double sensor_noise(double t, bool enabled) {
    return enabled ? 0.02 * std::sin(37.0 * t) : 0.0;
}

std::pair<double, double> derivatives(
    double t,
    double error,
    double theta_hat,
    const std::string& mode,
    bool disturbed,
    const Config& cfg
) {
    const double phi = regressor(t);
    const double measured_error = error + sensor_noise(t, disturbed);
    double gradient = -cfg.gamma * cfg.b * phi * measured_error;
    double leakage = 0.0;

    if (mode == "plain") {
        leakage = 0.0;
    } else if (mode == "sigma") {
        leakage = cfg.gamma * cfg.sigma * theta_hat;
    } else if (mode == "e") {
        leakage = cfg.gamma * cfg.sigma_e * std::abs(measured_error) * theta_hat;
    } else if (mode == "deadzone") {
        if (std::abs(measured_error) <= cfg.dead_zone) {
            gradient = 0.0;
        }
    } else {
        throw std::invalid_argument("Unknown mode: " + mode);
    }

    const double error_dot =
        -cfg.a_m * error
        + cfg.b * phi * (theta_hat - cfg.theta_star)
        + disturbance(t, disturbed);
    const double theta_dot = gradient - leakage;
    return {error_dot, theta_dot};
}

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

    Result result;
    result.time.resize(steps);
    result.error.assign(steps, 0.0);
    result.theta_hat.assign(steps, 0.0);
    result.error[0] = 1.0;

    for (std::size_t k = 0; k + 1U < steps; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        result.time[k] = t;
        const double e = result.error[k];
        const double th = result.theta_hat[k];

        const auto k1 = derivatives(t, e, th, mode, disturbed, cfg);
        const double mid_t = t + 0.5 * cfg.dt;
        const double mid_e = e + 0.5 * cfg.dt * k1.first;
        const double mid_th = th + 0.5 * cfg.dt * k1.second;
        const auto k2 = derivatives(mid_t, mid_e, mid_th, mode, disturbed, cfg);

        result.error[k + 1U] = e + cfg.dt * k2.first;
        result.theta_hat[k + 1U] = th + cfg.dt * k2.second;
    }
    result.time.back() = cfg.final_time;
    return result;
}

void print_metrics(const std::string& mode, const Result& result, const Config& cfg) {
    double sum_sq = 0.0;
    double max_abs = 0.0;
    std::size_t count = 0U;

    for (std::size_t k = 0; k < result.time.size(); ++k) {
        if (result.time[k] >= cfg.final_time - 10.0) {
            const double abs_error = std::abs(result.error[k]);
            sum_sq += result.error[k] * result.error[k];
            max_abs = std::max(max_abs, abs_error);
            ++count;
        }
    }

    const double rms = std::sqrt(sum_sq / static_cast<double>(count));
    const double final_theta = result.theta_hat.back();
    std::cout << std::setw(8) << mode
              << "  RMS(last 10 s)=" << rms
              << "  max|e|=" << max_abs
              << "  theta_hat(T)=" << final_theta
              << "  parameter error=" << final_theta - cfg.theta_star
              << '\n';
}

void save_csv(
    const std::string& filename,
    const std::vector<std::string>& modes,
    const std::map<std::string, Result>& results
) {
    std::ofstream file(filename);
    if (!file) {
        throw std::runtime_error("Cannot open " + filename);
    }

    file << "time";
    for (const auto& mode : modes) {
        file << ',' << mode << "_error," << mode << "_theta_hat";
    }
    file << '\n';

    const auto& time = results.at(modes.front()).time;
    file << std::setprecision(12);
    for (std::size_t k = 0; k < time.size(); ++k) {
        file << time[k];
        for (const auto& mode : modes) {
            file << ',' << results.at(mode).error[k]
                 << ',' << results.at(mode).theta_hat[k];
        }
        file << '\n';
    }
}

void run_case(bool disturbed, const Config& cfg) {
    const std::vector<std::string> modes = {
        "plain", "sigma", "e", "deadzone"
    };
    std::map<std::string, Result> results;

    for (const auto& mode : modes) {
        results.emplace(mode, simulate(mode, disturbed, cfg));
    }

    const std::string label = disturbed ? "disturbed" : "ideal";
    std::cout << "\nCase: " << label << '\n';
    for (const auto& mode : modes) {
        print_metrics(mode, results.at(mode), cfg);
    }

    save_csv("Chapter8_Lesson4_" + label + "_cpp.csv", modes, results);
}

int main() {
    try {
        const Config cfg;
        run_case(false, cfg);
        run_case(true, cfg);
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
    return 0;
}

Chapter8_Lesson4.java

// Chapter8_Lesson4.java
// Trade-off between robustness and parameter convergence in scalar MRAC.
// Build and run:
//   javac Chapter8_Lesson4.java
//   java Chapter8_Lesson4

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.LinkedHashMap;
import java.util.Locale;
import java.util.Map;

public final class Chapter8_Lesson4 {
    private Chapter8_Lesson4() {}

    private static final class Config {
        final double aM = 2.0;
        final double b = 1.0;
        final double gamma = 8.0;
        final double thetaStar = 1.5;
        final double sigma = 0.15;
        final double sigmaE = 0.80;
        final double deadZone = 0.08;
        final double dt = 0.001;
        final double finalTime = 40.0;
    }

    private static final class Result {
        final double[] time;
        final double[] error;
        final double[] thetaHat;

        Result(int steps) {
            time = new double[steps];
            error = new double[steps];
            thetaHat = new double[steps];
        }
    }

    private static double regressor(double t) {
        return 1.0 + 0.8 * Math.sin(0.7 * t);
    }

    private static double disturbance(double t, boolean enabled) {
        if (!enabled) {
            return 0.0;
        }
        double square = Math.sin(0.35 * t) >= 0.0 ? 1.0 : -1.0;
        return 0.12 * Math.sin(5.0 * t) + 0.06 * square;
    }

    private static double sensorNoise(double t, boolean enabled) {
        return enabled ? 0.02 * Math.sin(37.0 * t) : 0.0;
    }

    private static double[] derivatives(
        double t,
        double error,
        double thetaHat,
        String mode,
        boolean disturbed,
        Config cfg
    ) {
        double phi = regressor(t);
        double measuredError = error + sensorNoise(t, disturbed);
        double gradient = -cfg.gamma * cfg.b * phi * measuredError;
        double leakage = 0.0;

        switch (mode) {
            case "plain":
                break;
            case "sigma":
                leakage = cfg.gamma * cfg.sigma * thetaHat;
                break;
            case "e":
                leakage =
                    cfg.gamma * cfg.sigmaE * Math.abs(measuredError) * thetaHat;
                break;
            case "deadzone":
                if (Math.abs(measuredError) <= cfg.deadZone) {
                    gradient = 0.0;
                }
                break;
            default:
                throw new IllegalArgumentException("Unknown mode: " + mode);
        }

        double errorDot =
            -cfg.aM * error
            + cfg.b * phi * (thetaHat - cfg.thetaStar)
            + disturbance(t, disturbed);
        double thetaDot = gradient - leakage;
        return new double[] {errorDot, thetaDot};
    }

    private static Result simulate(
        String mode,
        boolean disturbed,
        Config cfg
    ) {
        int steps = (int) Math.round(cfg.finalTime / cfg.dt) + 1;
        Result result = new Result(steps);
        result.error[0] = 1.0;

        for (int k = 0; k < steps - 1; ++k) {
            double t = k * cfg.dt;
            result.time[k] = t;
            double e = result.error[k];
            double theta = result.thetaHat[k];

            double[] k1 = derivatives(
                t, e, theta, mode, disturbed, cfg
            );
            double midT = t + 0.5 * cfg.dt;
            double midE = e + 0.5 * cfg.dt * k1[0];
            double midTheta = theta + 0.5 * cfg.dt * k1[1];
            double[] k2 = derivatives(
                midT, midE, midTheta, mode, disturbed, cfg
            );

            result.error[k + 1] = e + cfg.dt * k2[0];
            result.thetaHat[k + 1] = theta + cfg.dt * k2[1];
        }
        result.time[steps - 1] = cfg.finalTime;
        return result;
    }

    private static void printMetrics(
        String mode,
        Result result,
        Config cfg
    ) {
        double sumSquared = 0.0;
        double maxAbsolute = 0.0;
        int count = 0;

        for (int k = 0; k < result.time.length; ++k) {
            if (result.time[k] >= cfg.finalTime - 10.0) {
                double absolute = Math.abs(result.error[k]);
                sumSquared += result.error[k] * result.error[k];
                maxAbsolute = Math.max(maxAbsolute, absolute);
                ++count;
            }
        }

        double rms = Math.sqrt(sumSquared / count);
        double finalTheta = result.thetaHat[result.thetaHat.length - 1];
        System.out.printf(
            Locale.US,
            "%8s  RMS(last 10 s)=%.8f  max|e|=%.8f"
                + "  theta_hat(T)=%.8f  parameter error=%.8f%n",
            mode,
            rms,
            maxAbsolute,
            finalTheta,
            finalTheta - cfg.thetaStar
        );
    }

    private static void saveCsv(
        Path path,
        String[] modes,
        Map<String, Result> results
    ) throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(path)) {
            writer.write("time");
            for (String mode : modes) {
                writer.write("," + mode + "_error," + mode + "_theta_hat");
            }
            writer.newLine();

            Result first = results.get(modes[0]);
            for (int k = 0; k < first.time.length; ++k) {
                writer.write(String.format(Locale.US, "%.12f", first.time[k]));
                for (String mode : modes) {
                    Result result = results.get(mode);
                    writer.write(
                        String.format(
                            Locale.US,
                            ",%.12f,%.12f",
                            result.error[k],
                            result.thetaHat[k]
                        )
                    );
                }
                writer.newLine();
            }
        }
    }

    private static void runCase(
        boolean disturbed,
        Config cfg
    ) throws IOException {
        String[] modes = {"plain", "sigma", "e", "deadzone"};
        Map<String, Result> results = new LinkedHashMap<>();

        for (String mode : modes) {
            results.put(mode, simulate(mode, disturbed, cfg));
        }

        String label = disturbed ? "disturbed" : "ideal";
        System.out.println("\nCase: " + label);
        for (String mode : modes) {
            printMetrics(mode, results.get(mode), cfg);
        }

        saveCsv(
            Path.of("Chapter8_Lesson4_" + label + "_java.csv"),
            modes,
            results
        );
    }

    public static void main(String[] args) {
        Config cfg = new Config();
        try {
            runCase(false, cfg);
            runCase(true, cfg);
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

Chapter8_Lesson4.m

% Chapter8_Lesson4.m
% Trade-off between robustness and parameter convergence in scalar MRAC.
% This script compares plain, sigma-modified, e-modified, and dead-zone laws.

clear; clc; close all;

cfg.a_m = 2.0;
cfg.b = 1.0;
cfg.gamma = 8.0;
cfg.theta_star = 1.5;
cfg.sigma = 0.15;
cfg.sigma_e = 0.80;
cfg.dead_zone = 0.08;
cfg.dt = 0.001;
cfg.final_time = 40.0;

modes = ["plain", "sigma", "e", "deadzone"];
run_case(false, cfg, modes);
run_case(true, cfg, modes);

function run_case(disturbed, cfg, modes)
    results = cell(size(modes));
    for i = 1:numel(modes)
        results{i} = simulate_mrac(modes(i), disturbed, cfg);
    end

    if disturbed
        label = "disturbed";
    else
        label = "ideal";
    end

    fprintf('\nCase: %s\n', label);
    for i = 1:numel(modes)
        print_metrics(modes(i), results{i}, cfg);
    end

    save_case_csv("Chapter8_Lesson4_" + label + "_matlab.csv", ...
                  modes, results);

    figure('Name', "Chapter08 Lesson4 " + label);
    tiledlayout(2, 1);

    nexttile;
    hold on;
    for i = 1:numel(modes)
        plot(results{i}.time, results{i}.error, ...
             'DisplayName', modes(i));
    end
    grid on;
    ylabel('tracking error e');
    title("MRAC robustness/convergence comparison: " + label);
    legend('Location', 'best');

    nexttile;
    hold on;
    for i = 1:numel(modes)
        plot(results{i}.time, results{i}.theta_hat, ...
             'DisplayName', modes(i));
    end
    yline(cfg.theta_star, '--', 'ideal parameter');
    grid on;
    xlabel('time (s)');
    ylabel('parameter estimate');
    legend('Location', 'best');

    exportgraphics(gcf, ...
        "Chapter8_Lesson4_" + label + "_matlab.png", ...
        'Resolution', 180);
end

function result = simulate_mrac(mode, disturbed, cfg)
    steps = round(cfg.final_time / cfg.dt) + 1;
    time = linspace(0.0, cfg.final_time, steps)';
    error = zeros(steps, 1);
    theta_hat = zeros(steps, 1);
    error(1) = 1.0;

    for k = 1:(steps - 1)
        t = time(k);
        e = error(k);
        theta = theta_hat(k);

        [k1_e, k1_theta] = derivatives( ...
            t, e, theta, mode, disturbed, cfg);

        mid_t = t + 0.5 * cfg.dt;
        mid_e = e + 0.5 * cfg.dt * k1_e;
        mid_theta = theta + 0.5 * cfg.dt * k1_theta;

        [k2_e, k2_theta] = derivatives( ...
            mid_t, mid_e, mid_theta, mode, disturbed, cfg);

        error(k + 1) = e + cfg.dt * k2_e;
        theta_hat(k + 1) = theta + cfg.dt * k2_theta;
    end

    result.time = time;
    result.error = error;
    result.theta_hat = theta_hat;
end

function [error_dot, theta_dot] = derivatives( ...
    t, error, theta_hat, mode, disturbed, cfg)

    phi = 1.0 + 0.8 * sin(0.7 * t);

    if disturbed
        square_wave = 2.0 * (sin(0.35 * t) >= 0.0) - 1.0;
        d = 0.12 * sin(5.0 * t) + 0.06 * square_wave;
        noise = 0.02 * sin(37.0 * t);
    else
        d = 0.0;
        noise = 0.0;
    end

    measured_error = error + noise;
    gradient = -cfg.gamma * cfg.b * phi * measured_error;
    leakage = 0.0;

    switch mode
        case "plain"
            % No robust modification.
        case "sigma"
            leakage = cfg.gamma * cfg.sigma * theta_hat;
        case "e"
            leakage = cfg.gamma * cfg.sigma_e ...
                * abs(measured_error) * theta_hat;
        case "deadzone"
            if abs(measured_error) <= cfg.dead_zone
                gradient = 0.0;
            end
        otherwise
            error('Unknown mode: %s', mode);
    end

    error_dot = -cfg.a_m * error ...
        + cfg.b * phi * (theta_hat - cfg.theta_star) + d;
    theta_dot = gradient - leakage;
end

function print_metrics(mode, result, cfg)
    tail = result.time >= (cfg.final_time - 10.0);
    rms_error = sqrt(mean(result.error(tail).^2));
    max_error = max(abs(result.error(tail)));
    final_theta = result.theta_hat(end);

    fprintf(['%8s  RMS(last 10 s)=%.8f  max|e|=%.8f  ' ...
             'theta_hat(T)=%.8f  parameter error=%.8f\n'], ...
        mode, rms_error, max_error, final_theta, ...
        final_theta - cfg.theta_star);
end

function save_case_csv(filename, modes, results)
    output = table(results{1}.time, 'VariableNames', {'time'});
    for i = 1:numel(modes)
        output.(modes(i) + "_error") = results{i}.error;
        output.(modes(i) + "_theta_hat") = results{i}.theta_hat;
    end
    writetable(output, filename);
end

The following script creates a continuous-time Simulink model for the constant-regressor sigma-modified equations used in Section 5. It logs the tracking error and parameter estimate and compares the simulated final values with the analytical equilibrium.

Chapter8_Lesson4_Simulink.m

% Chapter8_Lesson4_Simulink.m
% Programmatically builds a Simulink model for sigma-modified scalar MRAC.
% The model uses constant regressor phi = 1 and zero disturbance so that the
% leakage-induced equilibrium bias can be compared directly with the formula
% derived in the lesson.

clear; clc;

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

new_system(model);
open_system(model);

% Parameters placed in the base workspace.
assignin('base', 'a_m', 2.0);
assignin('base', 'b', 1.0);
assignin('base', 'gamma', 8.0);
assignin('base', 'sigma', 0.15);
assignin('base', 'theta_star', 1.5);

add_block('simulink/Continuous/Integrator', ...
    model + "/Error state e", ...
    'InitialCondition', '1.0', ...
    'Position', [550 100 580 130]);

add_block('simulink/Continuous/Integrator', ...
    model + "/Parameter estimate theta_hat", ...
    'InitialCondition', '0.0', ...
    'Position', [550 260 580 290]);

add_block('simulink/Math Operations/Gain', ...
    model + "/-a_m", ...
    'Gain', '-a_m', ...
    'Position', [250 85 330 115]);

add_block('simulink/Sources/Constant', ...
    model + "/theta_star", ...
    'Value', 'theta_star', ...
    'Position', [90 185 150 215]);

add_block('simulink/Math Operations/Sum', ...
    model + "/theta_tilde", ...
    'Inputs', '+-', ...
    'Position', [245 180 275 220]);

add_block('simulink/Math Operations/Gain', ...
    model + "/b_phi", ...
    'Gain', 'b', ...
    'Position', [330 180 400 220]);

add_block('simulink/Math Operations/Sum', ...
    model + "/error_dot", ...
    'Inputs', '++', ...
    'Position', [465 105 495 145]);

add_block('simulink/Math Operations/Gain', ...
    model + "/-gamma_b", ...
    'Gain', '-gamma*b', ...
    'Position', [250 285 340 315]);

add_block('simulink/Math Operations/Gain', ...
    model + "/-gamma_sigma", ...
    'Gain', '-gamma*sigma', ...
    'Position', [250 350 350 380]);

add_block('simulink/Math Operations/Sum', ...
    model + "/theta_dot", ...
    'Inputs', '++', ...
    'Position', [465 285 495 325]);

add_block('simulink/Sinks/Scope', ...
    model + "/Scope", ...
    'NumInputPorts', '2', ...
    'Position', [690 155 730 235]);

add_block('simulink/Sinks/To Workspace', ...
    model + "/e_out", ...
    'VariableName', 'e_out', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [665 75 765 105]);

add_block('simulink/Sinks/To Workspace', ...
    model + "/theta_out", ...
    'VariableName', 'theta_out', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [665 275 765 305]);

% Error dynamics: e_dot = -a_m e + b(theta_hat - theta_star).
add_line(model, ...
    'Error state e/1', '-a_m/1', ...
    'autorouting', 'on');
add_line(model, ...
    '-a_m/1', 'error_dot/1', ...
    'autorouting', 'on');
add_line(model, ...
    'Parameter estimate theta_hat/1', 'theta_tilde/1', ...
    'autorouting', 'on');
add_line(model, ...
    'theta_star/1', 'theta_tilde/2', ...
    'autorouting', 'on');
add_line(model, ...
    'theta_tilde/1', 'b_phi/1', ...
    'autorouting', 'on');
add_line(model, ...
    'b_phi/1', 'error_dot/2', ...
    'autorouting', 'on');
add_line(model, ...
    'error_dot/1', 'Error state e/1', ...
    'autorouting', 'on');

% Adaptation: theta_dot = -gamma*b*e - gamma*sigma*theta_hat.
add_line(model, ...
    'Error state e/1', '-gamma_b/1', ...
    'autorouting', 'on');
add_line(model, ...
    '-gamma_b/1', 'theta_dot/1', ...
    'autorouting', 'on');
add_line(model, ...
    'Parameter estimate theta_hat/1', '-gamma_sigma/1', ...
    'autorouting', 'on');
add_line(model, ...
    '-gamma_sigma/1', 'theta_dot/2', ...
    'autorouting', 'on');
add_line(model, ...
    'theta_dot/1', 'Parameter estimate theta_hat/1', ...
    'autorouting', 'on');

% Logging and visualization.
add_line(model, ...
    'Error state e/1', 'Scope/1', ...
    'autorouting', 'on');
add_line(model, ...
    'Parameter estimate theta_hat/1', 'Scope/2', ...
    'autorouting', 'on');
add_line(model, ...
    'Error state e/1', 'e_out/1', ...
    'autorouting', 'on');
add_line(model, ...
    'Parameter estimate theta_hat/1', 'theta_out/1', ...
    'autorouting', 'on');

set_param(model, ...
    'StopTime', '20', ...
    'Solver', 'ode45', ...
    'MaxStep', '0.01');

save_system(model);
sim(model);

b_value = evalin('base', 'b');
a_value = evalin('base', 'a_m');
sigma_value = evalin('base', 'sigma');
theta_value = evalin('base', 'theta_star');

theta_equilibrium = ...
    (b_value^2 / (b_value^2 + a_value * sigma_value)) * theta_value;
error_equilibrium = ...
    -(sigma_value * b_value / ...
    (b_value^2 + a_value * sigma_value)) * theta_value;

fprintf('Predicted theta equilibrium: %.8f\n', theta_equilibrium);
fprintf('Predicted error equilibrium: %.8f\n', error_equilibrium);
fprintf('Simulated final theta: %.8f\n', theta_out.signals.values(end));
fprintf('Simulated final error: %.8f\n', e_out.signals.values(end));

open_system(model);

Chapter8_Lesson4.nb

Notebook[{
 Cell["Chapter 8, Lesson 4: Robustness vs. Parameter Convergence", "Title"],
 Cell["This notebook compares plain, sigma-modified, e-modified, and dead-zone scalar MRAC laws in ideal and disturbed cases.", "Text"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"cfg", "=", 
    RowBox[{"<|", 
     RowBox[{
      RowBox[{"\"am\"", "->", "2.0"}], ",",
      RowBox[{"\"b\"", "->", "1.0"}], ",",
      RowBox[{"\"gamma\"", "->", "8.0"}], ",",
      RowBox[{"\"thetaStar\"", "->", "1.5"}], ",",
      RowBox[{"\"sigma\"", "->", "0.15"}], ",",
      RowBox[{"\"sigmaE\"", "->", "0.80"}], ",",
      RowBox[{"\"deadZone\"", "->", "0.08"}], ",",
      RowBox[{"\"finalTime\"", "->", "40.0"}]
     }], "|>"}]}], ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"phi", "[", "t_", "]"}], ":=", 
   RowBox[{"1.0", "+", 
    RowBox[{"0.8", " ", RowBox[{"Sin", "[", RowBox[{"0.7", " ", "t"}], "]"}]}]}]}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"disturbance", "[", RowBox[{"t_", ",", "enabled_"}], "]"}], ":=",
   RowBox[{"If", "[",
    RowBox[{"enabled", ",",
     RowBox[{
      RowBox[{"0.12", " ", RowBox[{"Sin", "[", RowBox[{"5.0", " ", "t"}], "]"}]}], "+",
      RowBox[{"0.06", " ", RowBox[{"Sign", "[", RowBox[{"Sin", "[", RowBox[{"0.35", " ", "t"}], "]"}], "]"}]}]}], ",",
     "0.0"}], "]"}]}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"sensorNoise", "[", RowBox[{"t_", ",", "enabled_"}], "]"}], ":=",
   RowBox[{"If", "[",
    RowBox[{"enabled", ",",
     RowBox[{"0.02", " ", RowBox[{"Sin", "[", RowBox[{"37.0", " ", "t"}], "]"}]}], ",",
     "0.0"}], "]"}]}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"simulate", "[", RowBox[{"mode_String", ",", "disturbed_"}], "]"}], ":=",
   RowBox[{"Module", "[",
    RowBox[{
     RowBox[{"{", RowBox[{"am", ",", "b", ",", "gamma", ",", "thetaStar", ",",
       "sigma", ",", "sigmaE", ",", "deadZone", ",", "law", ",", "sol"}], "}"}], ",",
     RowBox[{
      RowBox[{"am", "=", RowBox[{"cfg", "[", "\"am\"", "]"}]}], ";",
      RowBox[{"b", "=", RowBox[{"cfg", "[", "\"b\"", "]"}]}], ";",
      RowBox[{"gamma", "=", RowBox[{"cfg", "[", "\"gamma\"", "]"}]}], ";",
      RowBox[{"thetaStar", "=", RowBox[{"cfg", "[", "\"thetaStar\"", "]"}]}], ";",
      RowBox[{"sigma", "=", RowBox[{"cfg", "[", "\"sigma\"", "]"}]}], ";",
      RowBox[{"sigmaE", "=", RowBox[{"cfg", "[", "\"sigmaE\"", "]"}]}], ";",
      RowBox[{"deadZone", "=", RowBox[{"cfg", "[", "\"deadZone\"", "]"}]}], ";",
      RowBox[{
       RowBox[{"law", "[", RowBox[{"t_", ",", "ee_", ",", "th_"}], "]"}], ":=",
       RowBox[{"Module", "[",
        RowBox[{
         RowBox[{"{", RowBox[{"em", ",", "gradient", ",", "leakage"}], "}"}], ",",
         RowBox[{
          RowBox[{"em", "=", RowBox[{"ee", "+", RowBox[{"sensorNoise", "[", RowBox[{"t", ",", "disturbed"}], "]"}]}]}], ";",
          RowBox[{"gradient", "=", RowBox[{"-", "gamma", " ", "b", " ", RowBox[{"phi", "[", "t", "]"}], " ", "em"}]}], ";",
          RowBox[{"leakage", "=", "0"}], ";",
          RowBox[{"Switch", "[",
           RowBox[{"mode", ",",
            "\"plain\"", ",", "Null", ",",
            "\"sigma\"", ",", RowBox[{"leakage", "=", RowBox[{"gamma", " ", "sigma", " ", "th"}]}], ",",
            "\"e\"", ",", RowBox[{"leakage", "=", RowBox[{"gamma", " ", "sigmaE", " ", RowBox[{"Abs", "[", "em", "]"}], " ", "th"}]}], ",",
            "\"deadzone\"", ",", RowBox[{"If", "[", RowBox[{RowBox[{RowBox[{"Abs", "[", "em", "]"}], "<=", "deadZone"}], ",", RowBox[{"gradient", "=", "0"}]}], "]"}]
           }], "]"}], ";",
          RowBox[{"gradient", "-", "leakage"}]
         }]}], "]"}]}], ";",
      RowBox[{"sol", "=", RowBox[{"NDSolveValue", "[",
       RowBox[{
        RowBox[{"{",
         RowBox[{
          RowBox[{
           RowBox[{"e", "'", "[", "t", "]"}], "==",
           RowBox[{
            RowBox[{"-", "am", " ", RowBox[{"e", "[", "t", "]"}]}], "+",
            RowBox[{"b", " ", RowBox[{"phi", "[", "t", "]"}], " ",
             RowBox[{"(", RowBox[{RowBox[{"th", "[", "t", "]"}], "-", "thetaStar"}], ")"}]}], "+",
            RowBox[{"disturbance", "[", RowBox[{"t", ",", "disturbed"}], "]"}]
           }]}], ",",
          RowBox[{
           RowBox[{"th", "'", "[", "t", "]"}], "==",
           RowBox[{"law", "[", RowBox[{"t", ",", RowBox[{"e", "[", "t", "]"}], ",", RowBox[{"th", "[", "t", "]"}]}], "]"}]}], ",",
          RowBox[{RowBox[{"e", "[", "0", "]"}], "==", "1.0"}], ",",
          RowBox[{RowBox[{"th", "[", "0", "]"}], "==", "0.0"}]
         }], "}"}], ",",
        RowBox[{"{", RowBox[{"e", ",", "th"}], "}"}], ",",
        RowBox[{"{", RowBox[{"t", ",", "0", ",", RowBox[{"cfg", "[", "\"finalTime\"", "]"}]}], "}"}],
        ",", RowBox[{"MaxStepFraction", "->", RowBox[{"1", "/", "5000"}]}]
       }], "]"}]}], ";",
      "sol"
     }]}], "]"}]}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"modes", "=", RowBox[{"{", RowBox[{"\"plain\"", ",", "\"sigma\"", ",", "\"e\"", ",", "\"deadzone\""}], "}"}]}], ";",
   RowBox[{"ideal", "=", RowBox[{"AssociationMap", "[", RowBox[{RowBox[{RowBox[{"simulate", "[", RowBox[{"#", ",", "False"}], "]"}], "&"}], ",", "modes"}], "]"}]}], ";",
   RowBox[{"disturbed", "=", RowBox[{"AssociationMap", "[", RowBox[{RowBox[{RowBox[{"simulate", "[", RowBox[{"#", ",", "True"}], "]"}], "&"}], ",", "modes"}], "]"}]}], ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{"GraphicsGrid", "[",
   RowBox[{"{",
    RowBox[{
     RowBox[{"{",
      RowBox[{
       RowBox[{"Plot", "[",
        RowBox[{
         RowBox[{"Evaluate", "[", RowBox[{"Table", "[", RowBox[{RowBox[{RowBox[{"ideal", "[", "m", "]"}], "[", RowBox[{"[", "1", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"m", ",", "modes"}], "}"}]}], "]"}], "]"}], ",",
         RowBox[{"{", RowBox[{"t", ",", "0", ",", "40"}], "}"}], ",",
         RowBox[{"PlotLegends", "->", "modes"}], ",",
         RowBox[{"PlotLabel", "->", "\"Ideal: tracking error\""}]
        }], "]"}], ",",
       RowBox[{"Plot", "[",
        RowBox[{
         RowBox[{"Evaluate", "[", RowBox[{"Table", "[", RowBox[{RowBox[{RowBox[{"ideal", "[", "m", "]"}], "[", RowBox[{"[", "2", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"m", ",", "modes"}], "}"}]}], "]"}], "]"}], ",",
         RowBox[{"{", RowBox[{"t", ",", "0", ",", "40"}], "}"}], ",",
         RowBox[{"PlotLegends", "->", "modes"}], ",",
         RowBox[{"Epilog", "->", RowBox[{"{", RowBox[{"Dashed", ",", RowBox[{"InfiniteLine", "[", RowBox[{"{", RowBox[{"{", RowBox[{"0", ",", "1.5"}], "}"}], ",", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], "]"}]}], "}"}]}], ",",
         RowBox[{"PlotLabel", "->", "\"Ideal: parameter estimate\""}]
        }], "]"}]
      }], "}"}], ",",
     RowBox[{"{",
      RowBox[{
       RowBox[{"Plot", "[",
        RowBox[{
         RowBox[{"Evaluate", "[", RowBox[{"Table", "[", RowBox[{RowBox[{RowBox[{"disturbed", "[", "m", "]"}], "[", RowBox[{"[", "1", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"m", ",", "modes"}], "}"}]}], "]"}], "]"}], ",",
         RowBox[{"{", RowBox[{"t", ",", "0", ",", "40"}], "}"}], ",",
         RowBox[{"PlotLegends", "->", "modes"}], ",",
         RowBox[{"PlotLabel", "->", "\"Disturbed: tracking error\""}]
        }], "]"}], ",",
       RowBox[{"Plot", "[",
        RowBox[{
         RowBox[{"Evaluate", "[", RowBox[{"Table", "[", RowBox[{RowBox[{RowBox[{"disturbed", "[", "m", "]"}], "[", RowBox[{"[", "2", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"m", ",", "modes"}], "}"}]}], "]"}], "]"}], ",",
         RowBox[{"{", RowBox[{"t", ",", "0", ",", "40"}], "}"}], ",",
         RowBox[{"PlotLegends", "->", "modes"}], ",",
         RowBox[{"Epilog", "->", RowBox[{"{", RowBox[{"Dashed", ",", RowBox[{"InfiniteLine", "[", RowBox[{"{", RowBox[{"{", RowBox[{"0", ",", "1.5"}], "}"}], ",", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], "]"}]}], "}"}]}], ",",
         RowBox[{"PlotLabel", "->", "\"Disturbed: parameter estimate\""}]
        }], "]"}]
      }], "}"}]
    }], "}"}], "]"}]], "Input"]
},
WindowSize->{1200, 800},
StyleDefinitions->"Default.nb"
]

Problem 1: Ultimate Bound with Sigma-Modification

For

\[ \dot e = -a_m e + b \widetilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\phi} + d, \qquad |d|\le\bar d, \]

and the sigma-modified law of Section 4, prove that \(\dot V\le-\alpha V+C\) and state an ultimate bound for \(V\).

Solution:

Substitute the error and adaptive dynamics into the derivative of \(V\). The gradient cross terms cancel, leaving

\[ \dot V = -a_m e^2 + ed - \sigma \widetilde{\boldsymbol{\theta}}^{\mathsf T} \widehat{\boldsymbol{\theta}}. \]

Write \(\widehat{\boldsymbol{\theta}} = \widetilde{\boldsymbol{\theta}} + \boldsymbol{\theta}^{*}\) and apply Young's inequality to the disturbance and ideal-parameter cross terms. This gives

\[ \dot V \le -\frac{a_m}{2}e^2 - \frac{\sigma}{2} \left\| \widetilde{\boldsymbol{\theta}} \right\|^2 + C, \]

\[ C = \frac{\bar d^2}{2a_m} + \frac{\sigma}{2} \left\| \boldsymbol{\theta}^{*} \right\|^2. \]

With \(\alpha=\min\{a_m, \sigma\lambda_{\min}(\boldsymbol{\Gamma})\}\) ,

\[ \dot V\le-\alpha V+C, \qquad \limsup_{t\longrightarrow\infty}V(t) \le \frac{C}{\alpha}. \]

Problem 2: Numerical Leakage Bias

Let \(a_m=2,\ b=1,\ \phi_0=1,\ \theta^{*}=1.5,\ \sigma=0.15\) . Calculate the equilibrium parameter estimate and tracking error.

Solution:

\[ \widehat\theta_{\mathrm{eq}} = \frac{1}{1+2(0.15)}(1.5) = \frac{1.5}{1.3} \approx 1.15385, \]

\[ e_{\mathrm{eq}} = - \frac{0.15}{1+2(0.15)}(1.5) \approx -0.17308. \]

The parameter estimate is attenuated toward zero and the tracking error must remain nonzero to balance the leakage term.

Problem 3: Maximum Sigma for a Bias Constraint

In the scalar constant-regressor model, require

\[ \frac{ \widehat\theta_{\mathrm{eq}} }{ \theta^{*} } \ge 0.9. \]

Find an upper bound on \(\sigma\).

Solution:

\[ \frac{b^2\phi_0^2} {b^2\phi_0^2+a_m\sigma} \ge 0.9. \]

Rearranging,

\[ b^2\phi_0^2 \ge 0.9b^2\phi_0^2 + 0.9a_m\sigma, \]

\[ \boxed{ \sigma \le \frac{ b^2\phi_0^2 }{ 9a_m } }. \]

This is only a nominal equilibrium-bias constraint. A larger sigma may still be required to satisfy robustness requirements under disturbances.

Problem 4: Dead-Zone Selection

A measured error contains at most \(0.02\) units of sensor noise. Nominal-model validation shows an unavoidable tracking residual of at most \(0.04\). Select a dead-zone radius with an additional \(0.02\) safety margin.

Solution:

\[ \delta = 0.02+0.04+0.02 = 0.08. \]

Adaptation is disabled for \(|e_m|\le0.08\). A smaller dead zone may react to noise; an unnecessarily large dead zone can freeze adaptation before useful correction is complete.

Problem 5: Why E-Modification Usually Has Less Nominal Bias

Compare the leakage magnitudes in sigma- and e-modification as \(|e_m|\) approaches zero.

Solution:

Sigma-modification applies the parameter-space damping

\[ \sigma \boldsymbol{\Gamma} \widehat{\boldsymbol{\theta}} \]

even when the tracking error is arbitrarily small. E-modification applies

\[ \sigma_e|e_m| \boldsymbol{\Gamma} \widehat{\boldsymbol{\theta}}, \]

which approaches zero with \(|e_m|\). Therefore the e-modified law can approach the plain adaptive law near successful tracking and often produces less nominal bias. The same property means that its robust damping also becomes weaker near the origin.

• Tracking-error convergence and parameter convergence are distinct closed-loop properties.
• Plain adaptation has the best nominal opportunity for exact parameter recovery but can react incorrectly to nonparametric residual errors.
• Sigma-modification supplies constant parameter damping and supports a direct uniform-ultimate-boundedness proof.
• Constant leakage generally creates nonzero parameter and tracking bias; the scalar equilibrium quantifies this effect exactly.
• E-modification reduces leakage near good tracking, whereas a dead zone freezes adaptation inside a known uncertainty band.
• Robustification should be selected from disturbance/noise evidence and an explicit acceptable-bias requirement, not from adaptation speed alone.