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

• explain why a nominal MRAC parameter law can drift under bounded disturbances and unmodeled dynamics;
• derive constant-leakage and error-dependent-leakage update laws;
• construct a Lyapunov proof that establishes bounded tracking and parameter errors;
• identify the robustness-versus-bias trade-off created by leakage;
• select and tune \( \sigma \) and \( \sigma_e \) for a specified MRAC structure; and
• implement the same design in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The required background is the direct MRAC structure, tracking-error dynamics, quadratic Lyapunov functions, and the Lyapunov equation from Chapters 3–7. No robust-control machinery outside the adaptive loop is required.

Consider a matched MRAC error model in which the reference-model matrix \( A_m \) is Hurwitz:

\[ \dot{e}=A_m e+B\widetilde{\theta}^{T}\phi+w, \qquad \widetilde{\theta}=\widehat{\theta}-\theta^{\ast}. \]

Here \( e\in\mathbb{R}^{n} \) is the tracking error, \( \phi\in\mathbb{R}^{p} \) is the known regressor, \( \widehat{\theta} \) is the adaptive parameter vector, \( \theta^{\ast} \) is the unknown ideal parameter vector, and \( w \) collects bounded disturbance and reduced-order modeling error. Assume

\[ \lVert w(t)\rVert\leq \overline{w}. \]

Let \( P=P^{T}\succ 0 \) solve the Lyapunov equation

\[ A_m^{T}P+PA_m=-Q, \qquad Q=Q^{T}\succ 0. \]

The nominal gradient law is

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

With no modeling error, this law cancels the state–parameter cross term in the Lyapunov derivative. When \( w\neq 0 \), however, a small persistent error can continuously integrate through the adaptation law. The parameter estimate may then move far from the set of useful controller parameters even though the disturbance itself is bounded. This phenomenon is called parameter drift.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> C["Adaptive controller"]
  P["Plant plus neglected dynamics"] --> Y["Measured output y"]
  C --> P
  RM --> E["Tracking error e = y - ym"]
  Y --> E
  E --> A["Nominal parameter update"]
  Y --> A
  R --> A
  A --> C
  W["Bounded disturbance \nand modeling error"] --> P
  W --> E
  E --> D["Small persistent error \ncan accumulate into \nparameter drift"]
        

Robust modifications do not remove \( w \). Instead, they alter the parameter law so that the estimator contains a restoring mechanism. The price is that the ideal parameter vector is generally no longer an exact equilibrium of the adaptive law.

The \( \sigma \)-modification adds a linear leakage term directed toward the origin of parameter space:

\[ \boxed{\dot{\widehat{\theta} } =-\Gamma\left(\phi e^{T}PB+\sigma\widehat{\theta}\right)}, \qquad \sigma>0. \]

The leakage is active at every instant, including intervals during which the tracking error is small. In an implementation with diagonal adaptation gains,

\[ \Gamma=\operatorname{diag}(\gamma_1,\ldots,\gamma_p), \]

the effective decay rate of parameter component \( \widehat{\theta}_i \) is \( \gamma_i\sigma \). Consequently, changing \( \Gamma \) changes both adaptation speed and leakage speed under this parameterization.

At exact tracking, the gradient term vanishes, but the leakage remains:

\[ e=0 \quad \Longrightarrow \quad \dot{\widehat{\theta} }=-\Gamma\sigma\widehat{\theta}. \]

Therefore, unless the ideal parameter vector is the origin, constant leakage tends to move the estimate away from the ideal matching value. The resulting tracking error activates adaptation again. This mechanism explains why the method can suppress drift while producing a nonzero steady-state error neighborhood.

Choose the composite Lyapunov function

\[ V=e^{T}Pe+\widetilde{\theta}^{T}\Gamma^{-1}\widetilde{\theta}. \]

Since \( \theta^{\ast} \) is constant, \( \dot{\widetilde{\theta} }=\dot{\widehat{\theta} } \). Differentiating and substituting the error model gives

\[ \begin{aligned} \dot V &=e^{T}(A_m^{T}P+PA_m)e +2e^{T}PB\widetilde{\theta}^{T}\phi +2e^{T}Pw +2\widetilde{\theta}^{T}\Gamma^{-1} \dot{\widehat{\theta} } \\ &=-e^{T}Qe+2e^{T}Pw -2\sigma\widetilde{\theta}^{T}\widehat{\theta}. \end{aligned} \]

The identity

\[ 2\widetilde{\theta}^{T}\widehat{\theta} =\lVert\widetilde{\theta}\rVert^{2} +\lVert\widehat{\theta}\rVert^{2} -\lVert\theta^{\ast}\rVert^{2} \]

implies

\[ -2\sigma\widetilde{\theta}^{T}\widehat{\theta} \leq -\sigma\lVert\widetilde{\theta}\rVert^{2} +\sigma\lVert\theta^{\ast}\rVert^{2}. \]

Let \( \lambda_Q=\lambda_{\min}(Q) \). Using the disturbance bound and Young's inequality,

\[ 2e^{T}Pw \leq 2\lVert P\rVert\overline{w}\lVert e\rVert \leq \frac{\lambda_Q}{2}\lVert e\rVert^{2} +\frac{2\lVert P\rVert^{2}\overline{w}^{2} }{\lambda_Q}. \]

Hence

\[ \dot V \leq -\frac{\lambda_Q}{2}\lVert e\rVert^{2} -\sigma\lVert\widetilde{\theta}\rVert^{2} +C_{\sigma}, \]

\[ C_{\sigma}=\sigma\lVert\theta^{\ast}\rVert^{2} +\frac{2\lVert P\rVert^{2}\overline{w}^{2} }{\lambda_Q}. \]

Define

\[ \alpha_{\sigma}=\min\left\{ \frac{\lambda_Q}{2\lambda_{\max}(P)}, \sigma\lambda_{\min}(\Gamma) \right\}. \]

Then

\[ \dot V\leq -\alpha_{\sigma}V+C_{\sigma}. \]

By the comparison principle,

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

and therefore

\[ \limsup_{t\rightarrow\infty}V(t) \leq \frac{C_{\sigma} }{\alpha_{\sigma} }. \]

Thus the tracking and parameter errors are uniformly ultimately bounded. Increasing \( \sigma \) strengthens direct parameter damping, but also enlarges the term involving \( \lVert\theta^{\ast}\rVert \). The bound therefore displays the central design trade-off rather than implying that a larger leakage coefficient is always better.

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

\[ \boxed{\dot{\widehat{\theta} } =-\Gamma\left(\phi e^{T}PB +\sigma_e\lVert e\rVert\widehat{\theta}\right)}, \qquad \sigma_e>0. \]

For a scalar tracking error, \( \lVert e\rVert=|e| \). The leakage is strong during a large tracking transient and weak near accurate tracking. In particular,

\[ e=0 \quad \Longrightarrow \quad \dot{\widehat{\theta} }=0. \]

Unlike constant leakage, e-modification does not force a useful nonzero estimate to decay when exact tracking has already been reached. This commonly reduces leakage-induced steady-state bias. The limitation is equally important: when the measured error is very small, the restoring action against parameter drift is also very small.

In noisy measurements, replacing \( |e| \) by a filtered error magnitude can avoid rapidly varying leakage. That filter must be included in the implemented dynamics and in any rigorous proof; it cannot be added only in software without changing the mathematical system.

Using the same Lyapunov function and the e-modification update law, cross-term cancellation yields

\[ \dot V=-e^{T}Qe+2e^{T}Pw -2\sigma_e\lVert e\rVert \widetilde{\theta}^{T}\widehat{\theta}. \]

Applying the parameter identity gives

\[ \dot V\leq -\lambda_Q\lVert e\rVert^{2} +2\lVert P\rVert\overline{w}\lVert e\rVert -\sigma_e\lVert e\rVert \lVert\widetilde{\theta}\rVert^{2} +\sigma_e\lVert e\rVert \lVert\theta^{\ast}\rVert^{2}. \]

Define

\[ c_e=2\lVert P\rVert\overline{w} +\sigma_e\lVert\theta^{\ast}\rVert^{2}. \]

Then

\[ \dot V\leq -\lVert e\rVert\left( \lambda_Q\lVert e\rVert +\sigma_e\lVert\widetilde{\theta}\rVert^{2} -c_e\right). \]

Therefore, away from \( e=0 \), the Lyapunov function decreases whenever

\[ \lambda_Q\lVert e\rVert +\sigma_e\lVert\widetilde{\theta}\rVert^{2}>c_e. \]

A conservative residual region is consequently described by

\[ \mathcal{R}_e=\left\{(e,\widetilde{\theta}): \lambda_Q\lVert e\rVert +\sigma_e\lVert\widetilde{\theta}\rVert^{2} \leq c_e\right\}. \]

This proof highlights a structural difference from constant leakage. The damping of \( \widetilde{\theta} \) is multiplied by \( \lVert e\rVert \); the method acts most strongly precisely when poor tracking signals that the controller is far from satisfactory. The displayed residual set is a sufficient bound, not an exact prediction of the observed steady-state error.

flowchart TD
  S["Start from a verified \nnominal MRAC design"] --> M["Add representative bounded disturbance \nand neglected dynamics"]
  M --> D["Does the plain parameter \nestimate drift or burst?"]
  D -->|"No"| K["Keep nominal law and \ndocument tested uncertainty range"]
  D -->|"Yes"| B["Is continuous parameter \ndamping acceptable?"]
  B -->|"Yes"| SG["Select sigma-modification"]
  B -->|"No"| EM["Select e-modification"]
  SG --> T1["Increase sigma gradually while \nchecking tracking bias and control effort"]
  EM --> T2["Increase sigma_e gradually while c\nhecking small-error drift"]
  T1 --> V["Validate error, parameter norm, \nand control signal over worst tested cases"]
  T2 --> V
        

Practical tuning sequence

1. Retain the nominal adaptation gain \( \Gamma \) that gives acceptable ideal-model transients.
2. Introduce realistic bounded disturbances and the neglected dynamics identified in Lesson 1.
3. Begin with a small leakage coefficient and inspect tracking error, parameter norm, and control activity simultaneously.
4. Increase leakage only until parameter drift is acceptably suppressed.
5. Repeat the test over initial conditions, command changes, disturbance phases, and uncertain plant parameters.

A leakage coefficient cannot be selected independently of signal and parameter scaling. If one parameter component is numerically much larger than another, an isotropic term \( \sigma\widehat{\theta} \) may damp the components unevenly in physical terms. Parameter scaling should therefore be fixed before final tuning.

To compare the three laws under identical conditions, consider the actual plant

\[ \dot y=-a_p y+b_p z+d(t), \qquad \tau_a\dot z=-z+u. \]

The controller is designed from the reduced model that neglects the actuator state \( z \). The reference model and control law are

\[ \dot y_m=-a_m y_m+b_m r, \qquad u=\widehat{\theta}_y y+\widehat{\theta}_r r. \]

If \( \tau_a=0 \) and \( d=0 \), the matching parameters are

\[ \theta_y^{\ast}=\frac{a_p-a_m}{b_p}, \qquad \theta_r^{\ast}=\frac{b_m}{b_p}. \]

For the numerical values

\[ a_p=0.8,\quad b_p=1.2,\quad a_m=b_m=1.5,\quad \tau_a=0.08, \]

the reduced-model ideal parameters are

\[ \theta_y^{\ast}=-0.5833, \qquad \theta_r^{\ast}=1.25. \]

The regressor and scalar error are \( \phi=[y\;r]^{T} \) and \( e=y-y_m \). The implementations use

\[ \begin{aligned} \dot{\widehat{\theta} }_{\mathrm{plain} } &=-\gamma\phi e,\\ \dot{\widehat{\theta} }_{\sigma} &=-\gamma\left(\phi e+\sigma\widehat{\theta}_{\sigma}\right),\\ \dot{\widehat{\theta} }_{e} &=-\gamma\left(\phi e+\sigma_e|e|\widehat{\theta}_{e}\right), \end{aligned} \]

with

\[ \gamma=6,\qquad \sigma=0.08,\qquad \sigma_e=0.20. \]

The bounded disturbance is

\[ d(t)=0.12\sin(2.3t)+0.05\sin(17t) +0.10\,\mathbf{1}_{[12,\infty)}(t). \]

The actuator lag supplies a controlled example of unmodeled dynamics, while the sinusoidal and bias terms supply persistent bounded mismatch. Because the same fixed-step fourth-order Runge–Kutta solver and parameters are used in Python, C++, Java, and MATLAB, the reported metrics should agree up to floating-point and output-format differences.

For the included parameter set, the plain law generally obtains the smallest RMS tracking error but permits the largest parameter norm. The \( \sigma \)-modified law gives the strongest continuous pull toward the origin and therefore the smallest parameter excursion, at the cost of the largest tracking bias among the three tested laws. The e-modified law typically lies between those behaviors.

These numerical observations do not rank the methods universally. A different command, disturbance, ideal parameter magnitude, actuator lag, or adaptation gain can change the ordering. The mathematically important conclusions are structural:

• constant leakage provides damping even when the error is small;
• error-dependent leakage reduces its own bias near accurate tracking;
• both modifications replace exact nominal convergence claims with boundedness or residual-set guarantees under mismatch; and
• tracking performance alone is insufficient for evaluation because parameter and control trajectories can reveal hidden fragility.

The control commands in the examples include a large numerical safety limit. It is not used as part of the stability proof; it only prevents a failed numerical experiment from producing unbounded machine values.

The Python program uses NumPy for array operations and a from-scratch RK4 integrator. Matplotlib is optional and is used only for plotting. The CSV output is generated with the standard library.

Chapter8_Lesson2.py

"""Chapter 8, Lesson 2: sigma- and e-modification for scalar MRAC.

The script compares a plain gradient MRAC law with sigma-modification and
error-dependent e-modification.  It writes a CSV file and, when matplotlib is
available, a PNG figure.
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path
from typing import Callable, Dict, List

import numpy as np


@dataclass(frozen=True)
class Config:
    final_time: float = 30.0
    dt: float = 0.002
    a_plant: float = 0.8
    b_plant: float = 1.2
    a_model: float = 1.5
    b_model: float = 1.5
    actuator_tau: float = 0.08
    gamma: float = 6.0
    sigma: float = 0.08
    sigma_e: float = 0.20
    control_limit: float = 15.0


def reference(t: float) -> float:
    if t < 10.0:
        return 1.0
    if t < 20.0:
        return -0.7
    return 0.5


def disturbance(t: float) -> float:
    bias = 0.10 if t >= 12.0 else 0.0
    return 0.12 * np.sin(2.3 * t) + 0.05 * np.sin(17.0 * t) + bias


def derivative(t: float, state: np.ndarray, mode: str, cfg: Config) -> np.ndarray:
    y, y_m, actuator_state, theta_y, theta_r = state
    r = reference(t)
    e = y - y_m
    theta = np.array([theta_y, theta_r], dtype=float)
    phi = np.array([y, r], dtype=float)

    commanded_u = float(theta @ phi)
    u = float(np.clip(commanded_u, -cfg.control_limit, cfg.control_limit))

    y_dot = -cfg.a_plant * y + cfg.b_plant * actuator_state + disturbance(t)
    y_m_dot = -cfg.a_model * y_m + cfg.b_model * r
    actuator_dot = (-actuator_state + u) / cfg.actuator_tau

    gradient = phi * e
    if mode == "plain":
        theta_dot = -cfg.gamma * gradient
    elif mode == "sigma":
        theta_dot = -cfg.gamma * (gradient + cfg.sigma * theta)
    elif mode == "e":
        theta_dot = -cfg.gamma * (gradient + cfg.sigma_e * abs(e) * theta)
    else:
        raise ValueError(f"Unknown mode: {mode}")

    return np.array(
        [y_dot, y_m_dot, actuator_dot, theta_dot[0], theta_dot[1]],
        dtype=float,
    )


def rk4_step(
    f: Callable[[float, np.ndarray, str, Config], np.ndarray],
    t: float,
    state: np.ndarray,
    mode: str,
    cfg: Config,
) -> np.ndarray:
    dt = cfg.dt
    k1 = f(t, state, mode, cfg)
    k2 = f(t + 0.5 * dt, state + 0.5 * dt * k1, mode, cfg)
    k3 = f(t + 0.5 * dt, state + 0.5 * dt * k2, mode, cfg)
    k4 = f(t + dt, state + dt * k3, mode, cfg)
    return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(mode: str, cfg: Config) -> List[Dict[str, float | str]]:
    state = np.zeros(5, dtype=float)
    rows: List[Dict[str, float | str]] = []
    number_of_steps = int(round(cfg.final_time / cfg.dt))

    for step in range(number_of_steps + 1):
        t = step * cfg.dt
        y, y_m, actuator_state, theta_y, theta_r = state
        r = reference(t)
        e = y - y_m
        u_command = theta_y * y + theta_r * r
        u = float(np.clip(u_command, -cfg.control_limit, cfg.control_limit))

        rows.append(
            {
                "mode": mode,
                "t": t,
                "y": y,
                "y_m": y_m,
                "e": e,
                "theta_y": theta_y,
                "theta_r": theta_r,
                "u": u,
                "actuator_state": actuator_state,
                "disturbance": disturbance(t),
            }
        )

        if step < number_of_steps:
            state = rk4_step(derivative, t, state, mode, cfg)

    return rows


def write_csv(rows: List[Dict[str, float | str]], path: Path) -> None:
    fieldnames = list(rows[0].keys())
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.DictWriter(stream, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(rows)


def report(rows: List[Dict[str, float | str]]) -> None:
    for mode in ("plain", "sigma", "e"):
        selected = [row for row in rows if row["mode"] == mode]
        errors = np.array([float(row["e"]) for row in selected])
        parameter_norm = np.array(
            [np.hypot(float(row["theta_y"]), float(row["theta_r"])) for row in selected]
        )
        print(
            f"{mode:>5s}: RMS(e)={np.sqrt(np.mean(errors**2)):.5f}, "
            f"max|e|={np.max(np.abs(errors)):.5f}, "
            f"max||theta||={np.max(parameter_norm):.5f}"
        )


def make_plot(rows: List[Dict[str, float | str]], path: Path) -> None:
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print("matplotlib is unavailable; CSV output was still generated.")
        return

    figure, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
    for mode in ("plain", "sigma", "e"):
        selected = [row for row in rows if row["mode"] == mode]
        t = np.array([float(row["t"]) for row in selected])
        y = np.array([float(row["y"]) for row in selected])
        y_m = np.array([float(row["y_m"]) for row in selected])
        e = np.array([float(row["e"]) for row in selected])
        theta_norm = np.array(
            [np.hypot(float(row["theta_y"]), float(row["theta_r"])) for row in selected]
        )
        axes[0].plot(t, y, label=f"y ({mode})")
        axes[1].plot(t, e, label=mode)
        axes[2].plot(t, theta_norm, label=mode)

    # The reference trajectory is identical in all runs.
    selected = [row for row in rows if row["mode"] == "plain"]
    axes[0].plot(
        [float(row["t"]) for row in selected],
        [float(row["y_m"]) for row in selected],
        "k--",
        linewidth=1.5,
        label="reference model",
    )
    axes[0].set_ylabel("output")
    axes[1].set_ylabel("tracking error")
    axes[2].set_ylabel("parameter norm")
    axes[2].set_xlabel("time (s)")
    for axis in axes:
        axis.grid(True)
        axis.legend()
    figure.tight_layout()
    figure.savefig(path, dpi=180)
    print(f"Saved plot: {path}")


def main() -> None:
    cfg = Config()
    all_rows: List[Dict[str, float | str]] = []
    for mode in ("plain", "sigma", "e"):
        all_rows.extend(simulate(mode, cfg))

    output_dir = Path(__file__).resolve().parent
    csv_path = output_dir / "Chapter8_Lesson2_results.csv"
    plot_path = output_dir / "Chapter8_Lesson2_results.png"
    write_csv(all_rows, csv_path)
    report(all_rows)
    make_plot(all_rows, plot_path)
    print(f"Saved data: {csv_path}")


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library. It mirrors the Python state ordering, RK4 step, reference command, disturbance, and adaptation laws.

Chapter8_Lesson2.cpp

// Chapter 8, Lesson 2: sigma- and e-modification for scalar MRAC.
// Build: g++ -std=c++17 -O2 Chapter8_Lesson2.cpp -o Chapter8_Lesson2

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

struct Config {
    double finalTime = 30.0;
    double dt = 0.002;
    double aPlant = 0.8;
    double bPlant = 1.2;
    double aModel = 1.5;
    double bModel = 1.5;
    double actuatorTau = 0.08;
    double gamma = 6.0;
    double sigma = 0.08;
    double sigmaE = 0.20;
    double controlLimit = 15.0;
};

using State = std::array<double, 5>;  // y, y_m, actuator, theta_y, theta_r

struct Sample {
    std::string mode;
    double t;
    double y;
    double yModel;
    double error;
    double thetaY;
    double thetaR;
    double control;
    double actuator;
    double disturbance;
};

double reference(double t) {
    if (t < 10.0) return 1.0;
    if (t < 20.0) return -0.7;
    return 0.5;
}

double disturbance(double t) {
    const double bias = (t >= 12.0) ? 0.10 : 0.0;
    return 0.12 * std::sin(2.3 * t) + 0.05 * std::sin(17.0 * t) + bias;
}

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

State derivative(double t, const State& x, const std::string& mode, const Config& cfg) {
    const double y = x[0];
    const double yModel = x[1];
    const double actuator = x[2];
    const double thetaY = x[3];
    const double thetaR = x[4];
    const double r = reference(t);
    const double error = y - yModel;
    const double commandedControl = thetaY * y + thetaR * r;
    const double control = std::clamp(commandedControl, -cfg.controlLimit, cfg.controlLimit);

    State dx{};
    dx[0] = -cfg.aPlant * y + cfg.bPlant * actuator + disturbance(t);
    dx[1] = -cfg.aModel * yModel + cfg.bModel * r;
    dx[2] = (-actuator + control) / cfg.actuatorTau;

    const double gradientY = y * error;
    const double gradientR = r * error;
    if (mode == "plain") {
        dx[3] = -cfg.gamma * gradientY;
        dx[4] = -cfg.gamma * gradientR;
    } else if (mode == "sigma") {
        dx[3] = -cfg.gamma * (gradientY + cfg.sigma * thetaY);
        dx[4] = -cfg.gamma * (gradientR + cfg.sigma * thetaR);
    } else if (mode == "e") {
        const double leakage = cfg.sigmaE * std::abs(error);
        dx[3] = -cfg.gamma * (gradientY + leakage * thetaY);
        dx[4] = -cfg.gamma * (gradientR + leakage * thetaR);
    } else {
        throw std::invalid_argument("Unknown mode: " + mode);
    }
    return dx;
}

State rk4Step(double t, const State& x, const std::string& mode, const Config& cfg) {
    const State k1 = derivative(t, x, mode, cfg);
    const State k2 = derivative(t + 0.5 * cfg.dt, addScaled(x, k1, 0.5 * cfg.dt), mode, cfg);
    const State k3 = derivative(t + 0.5 * cfg.dt, addScaled(x, k2, 0.5 * cfg.dt), mode, cfg);
    const State k4 = derivative(t + cfg.dt, addScaled(x, k3, cfg.dt), mode, cfg);

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

std::vector<Sample> simulate(const std::string& mode, const Config& cfg) {
    State state{0.0, 0.0, 0.0, 0.0, 0.0};
    const int steps = static_cast<int>(std::llround(cfg.finalTime / cfg.dt));
    std::vector<Sample> samples;
    samples.reserve(static_cast<std::size_t>(steps + 1));

    for (int step = 0; step <= steps; ++step) {
        const double t = step * cfg.dt;
        const double r = reference(t);
        const double error = state[0] - state[1];
        const double control = std::clamp(
            state[3] * state[0] + state[4] * r,
            -cfg.controlLimit,
            cfg.controlLimit
        );
        samples.push_back({
            mode, t, state[0], state[1], error, state[3], state[4],
            control, state[2], disturbance(t)
        });
        if (step < steps) state = rk4Step(t, state, mode, cfg);
    }
    return samples;
}

void writeCsv(const std::vector<Sample>& samples, const std::string& filename) {
    std::ofstream stream(filename);
    if (!stream) throw std::runtime_error("Cannot open output file: " + filename);
    stream << "mode,t,y,y_m,e,theta_y,theta_r,u,actuator_state,disturbance\n";
    stream << std::setprecision(12);
    for (const auto& s : samples) {
        stream << s.mode << ',' << s.t << ',' << s.y << ',' << s.yModel << ','
               << s.error << ',' << s.thetaY << ',' << s.thetaR << ','
               << s.control << ',' << s.actuator << ',' << s.disturbance << '\n';
    }
}

void report(const std::vector<Sample>& samples, const std::string& mode) {
    double errorSquareSum = 0.0;
    double maxAbsError = 0.0;
    double maxParameterNorm = 0.0;
    std::size_t count = 0;
    for (const auto& s : samples) {
        if (s.mode != mode) continue;
        errorSquareSum += s.error * s.error;
        maxAbsError = std::max(maxAbsError, std::abs(s.error));
        maxParameterNorm = std::max(maxParameterNorm, std::hypot(s.thetaY, s.thetaR));
        ++count;
    }
    std::cout << std::setw(5) << mode
              << ": RMS(e)=" << std::sqrt(errorSquareSum / static_cast<double>(count))
              << ", max|e|=" << maxAbsError
              << ", max||theta||=" << maxParameterNorm << '\n';
}

int main() {
    try {
        const Config cfg;
        std::vector<Sample> allSamples;
        for (const std::string mode : {"plain", "sigma", "e"}) {
            auto samples = simulate(mode, cfg);
            allSamples.insert(allSamples.end(), samples.begin(), samples.end());
        }
        writeCsv(allSamples, "Chapter8_Lesson2_results_cpp.csv");
        report(allSamples, "plain");
        report(allSamples, "sigma");
        report(allSamples, "e");
        std::cout << "Saved Chapter8_Lesson2_results_cpp.csv\n";
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses records for output samples and no external numerical package. It is compatible with Java 21 and can be adapted to earlier Java versions by replacing the record with a conventional class.

Chapter8_Lesson2.java

// Chapter 8, Lesson 2: sigma- and e-modification for scalar MRAC.
// Build: javac Chapter8_Lesson2.java
// Run:   java Chapter8_Lesson2

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_Lesson2 {
    private static final class Config {
        final double finalTime = 30.0;
        final double dt = 0.002;
        final double aPlant = 0.8;
        final double bPlant = 1.2;
        final double aModel = 1.5;
        final double bModel = 1.5;
        final double actuatorTau = 0.08;
        final double gamma = 6.0;
        final double sigma = 0.08;
        final double sigmaE = 0.20;
        final double controlLimit = 15.0;
    }

    private record Sample(
        String mode,
        double t,
        double y,
        double yModel,
        double error,
        double thetaY,
        double thetaR,
        double control,
        double actuator,
        double disturbance
    ) {}

    private static double reference(double t) {
        if (t < 10.0) return 1.0;
        if (t < 20.0) return -0.7;
        return 0.5;
    }

    private static double disturbance(double t) {
        double bias = t >= 12.0 ? 0.10 : 0.0;
        return 0.12 * Math.sin(2.3 * t) + 0.05 * Math.sin(17.0 * t) + bias;
    }

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

    private static double[] derivative(double t, double[] x, String mode, Config cfg) {
        double y = x[0];
        double yModel = x[1];
        double actuator = x[2];
        double thetaY = x[3];
        double thetaR = x[4];
        double r = reference(t);
        double error = y - yModel;
        double control = clamp(thetaY * y + thetaR * r, -cfg.controlLimit, cfg.controlLimit);

        double[] dx = new double[5];
        dx[0] = -cfg.aPlant * y + cfg.bPlant * actuator + disturbance(t);
        dx[1] = -cfg.aModel * yModel + cfg.bModel * r;
        dx[2] = (-actuator + control) / cfg.actuatorTau;

        double gradientY = y * error;
        double gradientR = r * error;
        switch (mode) {
            case "plain" -> {
                dx[3] = -cfg.gamma * gradientY;
                dx[4] = -cfg.gamma * gradientR;
            }
            case "sigma" -> {
                dx[3] = -cfg.gamma * (gradientY + cfg.sigma * thetaY);
                dx[4] = -cfg.gamma * (gradientR + cfg.sigma * thetaR);
            }
            case "e" -> {
                double leakage = cfg.sigmaE * Math.abs(error);
                dx[3] = -cfg.gamma * (gradientY + leakage * thetaY);
                dx[4] = -cfg.gamma * (gradientR + leakage * thetaR);
            }
            default -> throw new IllegalArgumentException("Unknown mode: " + mode);
        }
        return dx;
    }

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

    private static double[] rk4Step(double t, double[] x, String mode, Config cfg) {
        double[] k1 = derivative(t, x, mode, cfg);
        double[] k2 = derivative(t + 0.5 * cfg.dt, addScaled(x, k1, 0.5 * cfg.dt), mode, cfg);
        double[] k3 = derivative(t + 0.5 * cfg.dt, addScaled(x, k2, 0.5 * cfg.dt), mode, cfg);
        double[] k4 = derivative(t + cfg.dt, addScaled(x, k3, cfg.dt), mode, cfg);
        double[] next = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            next[i] = x[i] + (cfg.dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        return next;
    }

    private static List<Sample> simulate(String mode, Config cfg) {
        double[] state = {0.0, 0.0, 0.0, 0.0, 0.0};
        int steps = (int) Math.round(cfg.finalTime / cfg.dt);
        List<Sample> samples = new ArrayList<>(steps + 1);

        for (int step = 0; step <= steps; step++) {
            double t = step * cfg.dt;
            double r = reference(t);
            double error = state[0] - state[1];
            double control = clamp(
                state[3] * state[0] + state[4] * r,
                -cfg.controlLimit,
                cfg.controlLimit
            );
            samples.add(new Sample(
                mode, t, state[0], state[1], error, state[3], state[4],
                control, state[2], disturbance(t)
            ));
            if (step < steps) state = rk4Step(t, state, mode, cfg);
        }
        return samples;
    }

    private static void writeCsv(List<Sample> samples, Path path) throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(path)) {
            writer.write("mode,t,y,y_m,e,theta_y,theta_r,u,actuator_state,disturbance\n");
            for (Sample s : samples) {
                writer.write(String.format(
                    Locale.ROOT,
                    "%s,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    s.mode(), s.t(), s.y(), s.yModel(), s.error(), s.thetaY(),
                    s.thetaR(), s.control(), s.actuator(), s.disturbance()
                ));
            }
        }
    }

    private static void report(List<Sample> samples, String mode) {
        double errorSquareSum = 0.0;
        double maxAbsError = 0.0;
        double maxParameterNorm = 0.0;
        int count = 0;
        for (Sample s : samples) {
            if (!s.mode().equals(mode)) continue;
            errorSquareSum += s.error() * s.error();
            maxAbsError = Math.max(maxAbsError, Math.abs(s.error()));
            maxParameterNorm = Math.max(maxParameterNorm, Math.hypot(s.thetaY(), s.thetaR()));
            count++;
        }
        System.out.printf(
            Locale.ROOT,
            "%5s: RMS(e)=%.5f, max|e|=%.5f, max||theta||=%.5f%n",
            mode, Math.sqrt(errorSquareSum / count), maxAbsError, maxParameterNorm
        );
    }

    public static void main(String[] args) {
        try {
            Config cfg = new Config();
            List<Sample> allSamples = new ArrayList<>();
            for (String mode : List.of("plain", "sigma", "e")) {
                allSamples.addAll(simulate(mode, cfg));
            }
            Path output = Path.of("Chapter8_Lesson2_results_java.csv");
            writeCsv(allSamples, output);
            report(allSamples, "plain");
            report(allSamples, "sigma");
            report(allSamples, "e");
            System.out.println("Saved " + output.toAbsolutePath());
        } catch (IOException | RuntimeException error) {
            System.err.println("Error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The first MATLAB file performs the same RK4 comparison and produces plots and a CSV table. The second file programmatically creates an educational fixed-step Simulink model containing a MATLAB Function block. Set the integer mode inside that block to select plain, \( \sigma \)-modified, or e-modified adaptation.

Chapter8_Lesson2.m

% Chapter8_Lesson2.m
% Sigma- and e-modification for scalar MRAC using a fixed-step RK4 solver.

clear; clc; close all;

cfg.finalTime = 30.0;
cfg.dt = 0.002;
cfg.aPlant = 0.8;
cfg.bPlant = 1.2;
cfg.aModel = 1.5;
cfg.bModel = 1.5;
cfg.actuatorTau = 0.08;
cfg.gamma = 6.0;
cfg.sigma = 0.08;
cfg.sigmaE = 0.20;
cfg.controlLimit = 15.0;

modes = {'plain', 'sigma', 'e'};
results = struct();
for k = 1:numel(modes)
    mode = modes{k};
    results.(mode) = simulateMode(mode, cfg);
end

figure('Name', 'Chapter 8 Lesson 2: Robust MRAC Modifications');
tiledlayout(3, 1);

nexttile; hold on; grid on;
plot(results.plain.t, results.plain.yModel, 'k--', 'LineWidth', 1.5, ...
    'DisplayName', 'reference model');
for k = 1:numel(modes)
    mode = modes{k};
    plot(results.(mode).t, results.(mode).y, 'DisplayName', mode);
end
ylabel('output'); legend('Location', 'best');

nexttile; hold on; grid on;
for k = 1:numel(modes)
    mode = modes{k};
    plot(results.(mode).t, results.(mode).error, 'DisplayName', mode);
end
ylabel('tracking error'); legend('Location', 'best');

nexttile; hold on; grid on;
for k = 1:numel(modes)
    mode = modes{k};
    thetaNorm = hypot(results.(mode).thetaY, results.(mode).thetaR);
    plot(results.(mode).t, thetaNorm, 'DisplayName', mode);
end
ylabel('parameter norm'); xlabel('time (s)'); legend('Location', 'best');

allRows = table();
for k = 1:numel(modes)
    mode = modes{k};
    data = results.(mode);
    modeColumn = repmat(string(mode), numel(data.t), 1);
    modeTable = table(modeColumn, data.t, data.y, data.yModel, data.error, ...
        data.thetaY, data.thetaR, data.control, data.actuator, data.disturbance, ...
        'VariableNames', {'mode','t','y','y_m','e','theta_y','theta_r', ...
        'u','actuator_state','disturbance'});
    allRows = [allRows; modeTable]; %#ok<AGROW>

    rmsError = sqrt(mean(data.error.^2));
    maxError = max(abs(data.error));
    maxTheta = max(hypot(data.thetaY, data.thetaR));
    fprintf('%5s: RMS(e)=%.5f, max|e|=%.5f, max||theta||=%.5f\n', ...
        mode, rmsError, maxError, maxTheta);
end
writetable(allRows, 'Chapter8_Lesson2_results_matlab.csv');

function data = simulateMode(mode, cfg)
    steps = round(cfg.finalTime / cfg.dt);
    t = (0:steps)' * cfg.dt;
    x = zeros(steps + 1, 5); % y, y_m, actuator, theta_y, theta_r

    for n = 1:steps
        tn = t(n);
        xn = x(n, :)';
        k1 = dynamics(tn, xn, mode, cfg);
        k2 = dynamics(tn + cfg.dt/2, xn + cfg.dt*k1/2, mode, cfg);
        k3 = dynamics(tn + cfg.dt/2, xn + cfg.dt*k2/2, mode, cfg);
        k4 = dynamics(tn + cfg.dt, xn + cfg.dt*k3, mode, cfg);
        x(n + 1, :) = (xn + cfg.dt*(k1 + 2*k2 + 2*k3 + k4)/6)';
    end

    r = arrayfun(@referenceSignal, t);
    d = arrayfun(@disturbanceSignal, t);
    uCommand = x(:,4).*x(:,1) + x(:,5).*r;
    u = min(max(uCommand, -cfg.controlLimit), cfg.controlLimit);

    data.t = t;
    data.y = x(:,1);
    data.yModel = x(:,2);
    data.actuator = x(:,3);
    data.thetaY = x(:,4);
    data.thetaR = x(:,5);
    data.error = x(:,1) - x(:,2);
    data.control = u;
    data.disturbance = d;
end

function dx = dynamics(t, x, mode, cfg)
    y = x(1);
    yModel = x(2);
    actuator = x(3);
    theta = x(4:5);
    r = referenceSignal(t);
    error = y - yModel;
    phi = [y; r];
    uCommand = theta' * phi;
    u = min(max(uCommand, -cfg.controlLimit), cfg.controlLimit);

    dx = zeros(5,1);
    dx(1) = -cfg.aPlant*y + cfg.bPlant*actuator + disturbanceSignal(t);
    dx(2) = -cfg.aModel*yModel + cfg.bModel*r;
    dx(3) = (-actuator + u) / cfg.actuatorTau;

    gradient = phi * error;
    switch mode
        case 'plain'
            thetaDot = -cfg.gamma * gradient;
        case 'sigma'
            thetaDot = -cfg.gamma * (gradient + cfg.sigma*theta);
        case 'e'
            thetaDot = -cfg.gamma * (gradient + cfg.sigmaE*abs(error)*theta);
        otherwise
            error('Unknown adaptation mode: %s', mode);
    end
    dx(4:5) = thetaDot;
end

function r = referenceSignal(t)
    if t < 10.0
        r = 1.0;
    elseif t < 20.0
        r = -0.7;
    else
        r = 0.5;
    end
end

function d = disturbanceSignal(t)
    if t >= 12.0
        bias = 0.10;
    else
        bias = 0.0;
    end
    d = 0.12*sin(2.3*t) + 0.05*sin(17.0*t) + bias;
end

Chapter8_Lesson2_Simulink.m

% Chapter8_Lesson2_Simulink.m
% Programmatically creates a discrete educational Simulink model for one
% selectable MRAC law: 1 = plain, 2 = sigma-modification, 3 = e-modification.
% This script requires Simulink and Stateflow (for the MATLAB Function block).

clear; clc;

modelName = 'Chapter8_Lesson2_Simulink_Model';
if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
new_system(modelName);
open_system(modelName);

set_param(modelName, ...
    'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '0.002', ...
    'StopTime', '30');

add_block('simulink/Sources/Clock', [modelName '/Clock'], ...
    'Position', [40 90 70 120]);
add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [modelName '/MRAC Simulation'], 'Position', [140 45 360 175]);
add_block('simulink/Signal Routing/Mux', [modelName '/Mux'], ...
    'Inputs', '6', 'Position', [420 45 425 185]);
add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
    'Position', [500 85 535 120]);
add_block('simulink/Sinks/To Workspace', [modelName '/To Workspace'], ...
    'VariableName', 'chapter8Lesson2Simulink', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [500 145 610 175]);

add_line(modelName, 'Clock/1', 'MRAC Simulation/1');
for port = 1:6
    add_line(modelName, sprintf('MRAC Simulation/%d', port), ...
        sprintf('Mux/%d', port));
end
add_line(modelName, 'Mux/1', 'Scope/1');
add_line(modelName, 'Mux/1', 'To Workspace/1');

rootObject = sfroot;
chart = find(rootObject, '-isa', 'Stateflow.EMChart', ...
    'Path', [modelName '/MRAC Simulation']);
if isempty(chart)
    error('The MATLAB Function block chart could not be located.');
end

chart.Script = sprintf([ ...
'function [y, yModel, trackingError, thetaY, thetaR, control] = fcn(t)\n' ...
'%%#codegen\n' ...
'persistent x previousTime\n' ...
'if isempty(x)\n' ...
'    x = zeros(5,1);\n' ...
'    previousTime = t;\n' ...
'end\n' ...
'dt = max(t - previousTime, 0.002);\n' ...
'previousTime = t;\n' ...
'mode = 2; %% 1 plain, 2 sigma, 3 e-modification\n' ...
'aPlant = 0.8; bPlant = 1.2; aModel = 1.5; bModel = 1.5;\n' ...
'actuatorTau = 0.08; gamma = 6.0; sigma = 0.08; sigmaE = 0.20;\n' ...
'if t < 10.0\n' ...
'    r = 1.0;\n' ...
'elseif t < 20.0\n' ...
'    r = -0.7;\n' ...
'else\n' ...
'    r = 0.5;\n' ...
'end\n' ...
'bias = 0.0;\n' ...
'if t >= 12.0\n' ...
'    bias = 0.10;\n' ...
'end\n' ...
'd = 0.12*sin(2.3*t) + 0.05*sin(17.0*t) + bias;\n' ...
'y = x(1); yModel = x(2); actuator = x(3); thetaY = x(4); thetaR = x(5);\n' ...
'trackingError = y - yModel;\n' ...
'control = min(max(thetaY*y + thetaR*r, -15.0), 15.0);\n' ...
'phi = [y; r];\n' ...
'gradient = phi*trackingError;\n' ...
'if mode == 1\n' ...
'    thetaDot = -gamma*gradient;\n' ...
'elseif mode == 2\n' ...
'    thetaDot = -gamma*(gradient + sigma*[thetaY; thetaR]);\n' ...
'else\n' ...
'    thetaDot = -gamma*(gradient + sigmaE*abs(trackingError)*[thetaY; thetaR]);\n' ...
'end\n' ...
'dx = [-aPlant*y + bPlant*actuator + d; ...\n' ...
'      -aModel*yModel + bModel*r; ...\n' ...
'      (-actuator + control)/actuatorTau; ...\n' ...
'      thetaDot];\n' ...
'x = x + dt*dx;\n' ...
'y = x(1); yModel = x(2); trackingError = y-yModel; thetaY=x(4); thetaR=x(5);\n' ...
'control = min(max(thetaY*y + thetaR*r, -15.0), 15.0);\n' ...
'end\n']);

save_system(modelName);
set_param(modelName, 'SimulationCommand', 'update');
open_system([modelName '/Scope']);
fprintf('Created %s.slx. Change mode inside the MATLAB Function block as needed.\n', modelName);

The notebook file loads the adjacent Wolfram Language source. The source uses NDSolveValue, produces tracking, error, and parameter plots, and exports sampled results to CSV.

Chapter8_Lesson2.nb

Notebook[{
  Cell["Chapter 8, Lesson 2: Sigma- and e-Modification", "Title"],
  Cell["This notebook loads and evaluates the adjacent Wolfram Language implementation file.", "Text"],
  Cell[BoxData[
    RowBox[{"Get", "[",
      RowBox[{"FileNameJoin", "[",
        RowBox[{"{",
          RowBox[{RowBox[{"NotebookDirectory", "[", "]"}], ",", "\"Chapter8_Lesson2.wl\""}],
        "}"}],
      "]"}],
    "]"}]], "Input"]
}, WindowTitle -> "Chapter8_Lesson2", StyleDefinitions -> "Default.nb"]

Chapter8_Lesson2.wl

(* Chapter8_Lesson2.wl
   Sigma- and e-modification for scalar MRAC. *)

ClearAll[reference, disturbance, simulate, modes, config];

config = <|
   "FinalTime" -> 30.0,
   "APlant" -> 0.8,
   "BPlant" -> 1.2,
   "AModel" -> 1.5,
   "BModel" -> 1.5,
   "ActuatorTau" -> 0.08,
   "Gamma" -> 6.0,
   "Sigma" -> 0.08,
   "SigmaE" -> 0.20,
   "ControlLimit" -> 15.0
|>;

reference[t_?NumericQ] := Piecewise[{ {1.0, t < 10.0}, {-0.7, t < 20.0} }, 0.5];
disturbance[t_?NumericQ] :=
  0.12 Sin[2.3 t] + 0.05 Sin[17.0 t] + Piecewise[{ {0.10, t >= 12.0} }, 0.0];

simulate[mode_String] := Module[
  {y, ym, z, thetaY, thetaR, e, r, phi, theta, u, leakage, equations, solution},
  equations = {
    e[t] == y[t] - ym[t],
    r[t] == reference[t],
    phi[t] == {y[t], r[t]},
    theta[t] == {thetaY[t], thetaR[t]},
    u[t] == Clip[theta[t].phi[t], {-config["ControlLimit"], config["ControlLimit"]}],
    y'[t] == -config["APlant"] y[t] + config["BPlant"] z[t] + disturbance[t],
    ym'[t] == -config["AModel"] ym[t] + config["BModel"] r[t],
    z'[t] == (-z[t] + u[t])/config["ActuatorTau"],
    leakage[t] == Switch[mode,
      "plain", 0.0,
      "sigma", config["Sigma"],
      "e", config["SigmaE"] Abs[e[t]],
      _, Message[simulate::mode, mode]; 0.0
    ],
    thetaY'[t] == -config["Gamma"] (y[t] e[t] + leakage[t] thetaY[t]),
    thetaR'[t] == -config["Gamma"] (r[t] e[t] + leakage[t] thetaR[t]),
    y[0] == 0.0, ym[0] == 0.0, z[0] == 0.0,
    thetaY[0] == 0.0, thetaR[0] == 0.0
  };
  solution = NDSolveValue[
    equations,
    {y, ym, z, thetaY, thetaR, e, u},
    {t, 0.0, config["FinalTime"]},
    Method -> {"EquationSimplification" -> "Residual"},
    MaxStepFraction -> 1/3000
  ];
  <|
    "y" -> solution[[1]], "ym" -> solution[[2]], "z" -> solution[[3]],
    "thetaY" -> solution[[4]], "thetaR" -> solution[[5]],
    "e" -> solution[[6]], "u" -> solution[[7]]
  |>
];

simulate::mode = "Unknown adaptation mode `1`.";
modes = {"plain", "sigma", "e"};
results = AssociationMap[simulate, modes];

trackingPlot = Plot[
  Evaluate@Join[
    Table[results[mode, "y"][t], {mode, modes}],
    {results["plain", "ym"][t]}
  ],
  {t, 0.0, config["FinalTime"]},
  PlotLegends -> Join[modes, {"reference model"}],
  AxesLabel -> {"time (s)", "output"},
  PlotRange -> All,
  ImageSize -> Large
];

errorPlot = Plot[
  Evaluate@Table[results[mode, "e"][t], {mode, modes}],
  {t, 0.0, config["FinalTime"]},
  PlotLegends -> modes,
  AxesLabel -> {"time (s)", "tracking error"},
  PlotRange -> All,
  ImageSize -> Large
];

parameterPlot = Plot[
  Evaluate@Table[
    Sqrt[results[mode, "thetaY"][t]^2 + results[mode, "thetaR"][t]^2],
    {mode, modes}
  ],
  {t, 0.0, config["FinalTime"]},
  PlotLegends -> modes,
  AxesLabel -> {"time (s)", "parameter norm"},
  PlotRange -> All,
  ImageSize -> Large
];

Print[trackingPlot];
Print[errorPlot];
Print[parameterPlot];

sampleTimes = Range[0.0, config["FinalTime"], 0.01];
rows = Flatten[
  Table[
    Table[
      {
        mode, time,
        results[mode, "y"][time],
        results[mode, "ym"][time],
        results[mode, "e"][time],
        results[mode, "thetaY"][time],
        results[mode, "thetaR"][time],
        results[mode, "u"][time],
        results[mode, "z"][time],
        disturbance[time]
      },
      {time, sampleTimes}
    ],
    {mode, modes}
  ],
  1
];
Export[
  FileNameJoin[{DirectoryName[$InputFileName], "Chapter8_Lesson2_results_mathematica.csv"}],
  Prepend[rows, {"mode", "t", "y", "y_m", "e", "theta_y", "theta_r", "u", "actuator_state", "disturbance"}]
];

Problem 1 (Cross-Term Cancellation): Starting from \( \dot e=A_m e+B\widetilde{\theta}^{T}\phi+w \) and \( V=e^{T}Pe+\widetilde{\theta}^{T}\Gamma^{-1}\widetilde{\theta} \), show explicitly that the nominal update law cancels the mixed term in \( \dot V \).

Solution: Differentiation gives

\[ \begin{aligned} \dot V &=e^{T}(A_m^{T}P+PA_m)e +2e^{T}PB\widetilde{\theta}^{T}\phi +2e^{T}Pw +2\widetilde{\theta}^{T}\Gamma^{-1} \dot{\widehat{\theta} }. \end{aligned} \]

Substituting

\[ \dot{\widehat{\theta} }=-\Gamma\phi e^{T}PB \]

produces

\[ 2\widetilde{\theta}^{T}\Gamma^{-1} \dot{\widehat{\theta} } =-2\widetilde{\theta}^{T}\phi e^{T}PB. \]

Since both mixed expressions are scalars, they are identical with opposite signs. Thus

\[ \dot V=-e^{T}Qe+2e^{T}Pw. \]

Problem 2 (Ultimate Bound for Constant Leakage): Given \( \dot V\leq-\alpha_{\sigma}V+C_{\sigma} \), derive the finite-time bound and the limiting bound on \( V \).

Solution: Multiply by the integrating factor

\[ e^{\alpha_{\sigma}t}. \]

Then

\[ \frac{d}{dt}\left(e^{\alpha_{\sigma}t}V(t)\right) \leq C_{\sigma}e^{\alpha_{\sigma}t}. \]

Integrating from \( 0 \) to \( t \) gives

\[ V(t)\leq e^{-\alpha_{\sigma}t}V(0) +\frac{C_{\sigma} }{\alpha_{\sigma} } \left(1-e^{-\alpha_{\sigma}t}\right). \]

Therefore

\[ \limsup_{t\rightarrow\infty}V(t) \leq \frac{C_{\sigma} }{\alpha_{\sigma} }. \]

Problem 3 (Why Constant Leakage Creates Bias): Suppose the plant is ideal, the tracking error is momentarily zero, and \( \widehat{\theta}=\theta^{\ast}\neq 0 \). Compare the immediate parameter derivatives under the two robust modifications.

Solution: For constant leakage,

\[ \dot{\widehat{\theta} }_{\sigma} =-\Gamma\sigma\theta^{\ast}\neq 0. \]

The estimate leaves the ideal value even though tracking is exact. For e-modification,

\[ \dot{\widehat{\theta} }_{e} =-\Gamma\sigma_e\lVert e\rVert\theta^{\ast}=0. \]

Thus e-modification preserves the ideal estimate at exact tracking, whereas constant leakage does not.

Problem 4 (Scalar Matching Parameters): For \( \dot y=-a_p y+b_p u \), \( \dot y_m=-a_m y_m+b_m r \), and \( u=\theta_y y+\theta_r r \), derive the matching parameters.

Solution: Substitution gives

\[ \dot y=(-a_p+b_p\theta_y)y+b_p\theta_r r. \]

Matching the coefficients to the reference-model dynamics requires

\[ -a_p+b_p\theta_y^{\ast}=-a_m, \qquad b_p\theta_r^{\ast}=b_m. \]

Therefore

\[ \theta_y^{\ast}=\frac{a_p-a_m}{b_p}, \qquad \theta_r^{\ast}=\frac{b_m}{b_p}. \]

Problem 5 (Leakage Scaling): Consider \( \Gamma=\operatorname{diag}(2,20) \) and \( \sigma=0.1 \). With the gradient term set to zero, find the exponential decay rates of the two parameter estimates.

Solution: The leakage-only dynamics are

\[ \dot{\widehat{\theta} } =-\Gamma\sigma\widehat{\theta}. \]

Thus

\[ \dot{\widehat{\theta} }_1=-0.2\widehat{\theta}_1, \qquad \dot{\widehat{\theta} }_2=-2\widehat{\theta}_2. \]

The second estimate decays ten times faster. This demonstrates why adaptation-gain and parameter scaling must be considered jointly.

Problem 6 (Design Choice): An adaptive controller shows slow drift during long periods of small nonzero tracking error. The ideal controller parameters are known to be far from the origin. Which modification is the more natural first test, and what failure mode must still be checked?

Solution: e-modification is a natural first test because it reduces leakage as the tracking error becomes small and therefore tends to introduce less origin-centered bias. However, its leakage also becomes weak in exactly that regime, so the long-duration small-error drift must be tested explicitly. If the drift persists, constant leakage or a later parameter-bounding method may be required.

Nominal MRAC cancels state–parameter cross terms but can integrate bounded mismatch into parameter drift. The \( \sigma \)-modification adds continuous linear leakage and yields a direct uniform-ultimate-boundedness inequality. Its constant damping is effective against drift but biases estimates toward the origin. The e-modification scales leakage by tracking-error magnitude, reducing forced parameter decay near exact tracking while providing weaker protection during small-error intervals. Both methods must be evaluated using tracking error, parameter trajectories, and control activity, and both require tuning against representative unmodeled dynamics and disturbances.