In Chapters 5–7, the adaptive law was selected so that parameter-error terms cancel in the derivative of a Lyapunov function. That cancellation is exact only when the assumed plant structure, measured regressor, and implemented control signal are exact. A physical loop instead contains sensor noise, dynamics outside the chosen model order, sampling and computation, communication latency, actuator lag, and saturation.

The central distinction is therefore between nominal stability and robust stability. A nominal theorem proves properties for a mathematical model. A robustness result quantifies what remains true for a family of perturbed models. Plain MRAC generally has strong nominal guarantees, but its adaptive integrators can accumulate small persistent errors and may inject energy into neglected fast dynamics.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  R --> AC["Adaptive controller"]
  YM["Measured output y plus noise"] --> AC
  RM --> E["Tracking error"]
  YM --> E
  AC --> D["Computation and input delay"]
  D --> A["Actuator and neglected fast dynamics"]
  A --> P["Nominal plant states"]
  P --> S["Sensor"]
  S --> YM
  E --> UP["Parameter update integrators"]
  UP --> AC
        

The feedback path through the parameter update is nonlinear and dynamic. Consequently, a perturbation that would be harmless for a fixed-gain controller can be integrated over time by the adaptive law. This is the origin of parameter drift and several classical adaptive-control instability mechanisms.

Consider the unknown first-order plant and desired reference model

\[ \dot{x}=-a x+b u, \qquad \dot{x}_m=-a_m x_m+b_m r, \qquad a_m>0,\; b>0. \]

The direct adaptive controller uses the regressor \( \phi=[x\; r]^T \):

\[ u=\hat{\theta}_x x+\hat{\theta}_r r=\hat{\theta}^{T}\phi. \]

The matching parameters satisfy

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

With tracking error \( e=x-x_m \) and parameter error \( \tilde{\theta}=\hat{\theta}-\theta^{*} \), the nominal error equation is

\[ \dot{e}=-a_m e+b\tilde{\theta}^{T}\phi. \]

Choose the gradient update law and Lyapunov candidate

\[ \dot{\hat{\theta}}=-\Gamma\phi e, \qquad V=\frac{1}{2}e^2+\frac{b}{2}\tilde{\theta}^{T}\Gamma^{-1}\tilde{\theta}, \qquad \Gamma=\Gamma^T>0. \]

Since the ideal parameters are constant, \( \dot{\tilde{\theta}}=\dot{\hat{\theta}} \), and therefore

\[ \begin{aligned} \dot{V} &=e\left(-a_m e+b\tilde{\theta}^{T}\phi\right) +b\tilde{\theta}^{T}\Gamma^{-1}\left(-\Gamma\phi e\right) \\ &=-a_m e^2\le 0. \end{aligned} \]

The cancellation is exact. It does not create a negative term in \( \tilde{\theta} \); bounded parameter estimates follow from bounded \( V \), while parameter convergence requires excitation conditions studied later in Chapter 10. More importantly here, any mismatch in the error equation or update law leaves a residual term in \( \dot V \).

Collect implementation imperfections into a scalar perturbation \( \Delta(t) \) acting on the tracking-error dynamics:

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

If the nominal update law still uses the exact regressor, then

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

Suppose a local bound is available in the operating region, \( |\Delta(t)|\le\delta_0+\delta_1|e| \). Then

\[ \dot{V}\le-\left(a_m-\delta_1\right)e^2+\delta_0|e|. \]

When \( a_m>\delta_1 \), the derivative is negative whenever

\[ |e|>\frac{\delta_0}{a_m-\delta_1}. \]

This calculation motivates an ultimate tracking-error bound, not asymptotic convergence. It is conditional: one must also establish bounded regressors and parameter estimates. Plain MRAC does not automatically provide those bounds when the perturbation enters the adaptive law itself, which is exactly what measurement noise does.

Using Young’s inequality, \( |e\Delta|\le\frac{\varepsilon}{2}e^2+\frac{1}{2\varepsilon}\Delta^2 \), gives the alternative estimate

\[ \dot{V}\le-\left(a_m-\frac{\varepsilon}{2}\right)e^2+\frac{\Delta^2}{2\varepsilon}, \qquad 0<\varepsilon<2a_m. \]

This inequality makes the performance trade-off explicit: increasing the nominal decay rate \( a_m \) can improve disturbance attenuation in the reduced model, but demanding a faster reference model also raises controller bandwidth and may excite neglected dynamics.

Let the measured output be \( y=x+n \), where \( n \) is zero-mean, band-limited sensor noise. The implemented controller and update law use

\[ \phi_y=\begin{bmatrix}y\\r\end{bmatrix}, \qquad e_y=y-x_m=e+n, \qquad \dot{\hat{\theta}}=-\Gamma\phi_y e_y. \]

The state-feedback component of the update contains

\[ y e_y=(x+n)(e+n)=xe+xn+en+n^2. \]

If \( n \) is zero mean and approximately independent of the slowly varying signals \( x \) and \( e \), then over an averaging window

\[ \mathbb{E}[y e_y]\approx xe+\mathbb{E}[n^2]=xe+\sigma_n^2. \]

Thus zero-mean measurement noise does not imply a zero-mean parameter update. The quadratic term \( n^2 \) generates a bias. For the scalar state-feedback estimate,

\[ \mathbb{E}[\dot{\hat{\theta}}_x]\approx-\gamma_x\left(xe+\sigma_n^2\right). \]

Near good tracking, \( xe \) becomes small while \( \sigma_n^2 \) remains positive. The parameter integrator can therefore continue moving even though the true tracking error is small. This is a basic drift mechanism. In addition, large adaptation gains multiply the high-frequency stochastic component, increasing parameter variance and control chatter.

Continuous-time white noise has infinite bandwidth and is not a physically realizable signal; practical analysis should use sampled or band-limited noise. Filtering can reduce high-frequency content, but every causal filter adds dynamics and phase lag, so filtering cannot be treated as cost-free.

Suppose the controller assumes that its command reaches the plant directly, but the physical actuator has a neglected second-order mode:

\[ \ddot{q}+2\zeta_f\omega_f\dot{q}+\omega_f^2 q=\omega_f^2 u, \qquad \dot{x}=-a x+bq. \]

The nominal design assumes \( q=u \). Define the actuator mismatch \( \eta=q-u \). Then

\[ \dot{e}=-a_m e+b\tilde{\theta}^{T}\phi+b\eta. \]

The Lyapunov derivative becomes \( \dot V=-a_m e^2+b e\eta \). The term \( \eta \) is not an independent bounded disturbance: it depends on the control command, and the command depends on adapting parameters. This feedback interconnection can amplify fast content.

The neglected transfer function is

\[ H_f(s)=\frac{q(s)}{u(s)}=\frac{\omega_f^2}{s^2+2\zeta_f\omega_f s+\omega_f^2}. \]

At frequencies well below \( \omega_f \), \( H_f(j\omega)\approx1 \). Near the parasitic mode, magnitude peaking and phase lag can invalidate the reduced-order matching argument. Increasing \( \Gamma \) makes the parameter estimates react faster, but it can also move substantial adaptive activity toward the neglected mode. Therefore “faster adaptation” is not synonymous with “better tracking.”

A useful engineering separation is to require the effective adaptation and closed-loop bandwidths to remain comfortably below the first important neglected mode. This is a design heuristic, not a universal theorem; the exact admissible separation depends on damping, loop gains, excitation, and the selected robust modification.

With an input delay \( \tau \), the actual plant receives \( u(t-\tau) \):

\[ \dot{x}(t)=-a x(t)+b u(t-\tau). \]

Add and subtract \( bu(t) \) to expose the nominal model plus a delay perturbation:

\[ \dot{e}=-a_m e+b\tilde{\theta}^{T}\phi+\Delta_{\tau}, \qquad \Delta_{\tau}=b\left[u(t-\tau)-u(t)\right]. \]

If \( u \) is differentiable, the fundamental theorem of calculus gives

\[ u(t-\tau)-u(t)=-\int_{t-\tau}^{t}\dot{u}(s)\,ds, \]

and hence

\[ |\Delta_{\tau}(t)|\le |b|\tau\sup_{s\in[t-\tau,t]}|\dot{u}(s)|. \]

Delay sensitivity therefore grows with both \( \tau \) and the rate of change of the adaptive control command. Because

\[ \dot{u}=\dot{\hat{\theta}}^{T}\phi+\hat{\theta}^{T}\dot{\phi}, \]

aggressive adaptation directly increases one component of the delay perturbation. In frequency-domain terms, a pure delay contributes phase \( -\omega\tau \) radians without attenuating magnitude. A fixed-gain loop with crossover \( \omega_c \) and phase margin \( \varphi_m \) has the familiar approximate screening condition

\[ \omega_c\tau<\varphi_m. \]

For an adaptive loop, \( \omega_c \) is not fixed because the controller parameters evolve. The inequality is therefore only a diagnostic guide. Formal adaptive delay margins require a dedicated robustness analysis and usually bounded parameter sets or projection, which is introduced in Chapter 9.

Noise, neglected modes, and delay rarely act independently. Their interaction creates a feedback chain: noise perturbs the update, the parameter integrators generate a rapidly varying command, delay and actuator modes convert that command into phase-lagged plant motion, and the resulting tracking error drives still more adaptation.

flowchart TD
  N["Sensor noise or small modeling residual"] --> U["Nonzero adaptive update near good tracking"]
  U --> P["Parameter motion or drift"]
  P --> C["Larger or faster control command"]
  C --> L["Delay phase loss"]
  C --> F["Excitation of neglected fast mode"]
  L --> E["Larger tracking error"]
  F --> E
  E --> U
  G["Higher adaptation gain"] --> U
        

Three practical symptoms should be distinguished:

Parameter drift: estimates move persistently even though tracking looks acceptable. Control chatter: high-frequency command activity grows because the update reacts to noisy regressors. Bursting or instability: long quiet periods are interrupted by large transients after the adaptive parameters enter an unfavorable region or excite a parasitic mode.

These symptoms motivate the robust modifications in the next lessons: leakage terms, error-dependent leakage, dead zones, and projection. This lesson deliberately does not use those modifications, so their necessity is visible.

The implementations use the following numerical benchmark:

\[ a=1.35,\quad b=1,\quad a_m=2,\quad b_m=2,\quad \Gamma=\operatorname{diag}(7,7). \]

The command is deliberately multi-frequency:

\[ r(t)=\mathbf{1}_{t\ge1}+0.25\sin(0.7t)+0.15\sin(2.1t). \]

Five cases are compared:

Nominal: exact input channel and noiseless measurement. Noise: measurement standard deviation \( 0.035 \). Unmodeled dynamics: a neglected actuator mode with \( \omega_f=18\;\mathrm{rad/s} \) and \( \zeta_f=0.18 \). Delay: input delay \( \tau=0.08\;\mathrm{s} \). Combined: all three imperfections.

The command is clipped to \( [-25,25] \) solely to keep the numerical experiment within a finite actuator range. The update law is not informed about clipping; hence the code is still a plain MRAC experiment rather than an anti-windup or robust-adaptive design.

Compare final-window RMS tracking error, maximum command magnitude, parameter trajectories, and frequency content. A small RMS error alone is insufficient: a controller may track while its parameters drift or its command develops unacceptable high-frequency energy.

Chapter8_Lesson1.py

"""Chapter8_Lesson1.py
Plain MRAC under measurement noise, unmodeled actuator dynamics, and input delay.
The controller intentionally uses no robust modification so the practical failure
mechanisms are visible. Requires: numpy, matplotlib.
"""
from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Scenario:
    name: str
    noise_std: float = 0.0
    actuator_omega: float | None = None
    actuator_zeta: float = 0.18
    delay: float = 0.0


def reference(t: float) -> float:
    return (1.0 if t >= 1.0 else 0.0) + 0.25 * np.sin(0.7 * t) + 0.15 * np.sin(2.1 * t)


def simulate(s: Scenario, *, dt: float = 1e-3, tf: float = 25.0, seed: int = 7) -> Dict[str, np.ndarray]:
    # Actual plant: x_dot = -a*x + b*v. The controller assumes v == u.
    a, b = 1.35, 1.0
    # Reference model: xm_dot = -am*xm + bm*r.
    am, bm = 2.0, 2.0
    gamma_x, gamma_r = 7.0, 7.0

    n = int(round(tf / dt)) + 1
    t = np.linspace(0.0, tf, n)
    rng = np.random.default_rng(seed)
    delay_steps = max(0, int(round(s.delay / dt)))
    u_buffer = np.zeros(delay_steps + 1)

    x = xm = 0.0
    # q is the output of a neglected second-order actuator/flexible mode.
    q = q_dot = 0.0
    theta_x = theta_r = 0.0

    out = {key: np.zeros(n) for key in ("t", "r", "x", "xm", "e", "ym", "u", "theta_x", "theta_r")}
    out["t"] = t

    for k, tk in enumerate(t):
        r = reference(tk)
        noise = s.noise_std * rng.standard_normal()
        ym = x + noise
        e_measured = ym - xm

        # Plain direct MRAC law. Saturation is only a numerical/actuator guard;
        # the adaptive law is deliberately unaware of saturation.
        u_cmd = float(np.clip(theta_x * ym + theta_r * r, -25.0, 25.0))
        u_buffer[-1] = u_cmd
        u_delayed = float(u_buffer[0])

        if s.actuator_omega is None:
            v = u_delayed
        else:
            omega, zeta = s.actuator_omega, s.actuator_zeta
            q_ddot = omega**2 * (u_delayed - q) - 2.0 * zeta * omega * q_dot
            q_dot += dt * q_ddot
            q += dt * q_dot
            v = q

        x_dot = -a * x + b * v
        xm_dot = -am * xm + bm * r
        theta_x_dot = -gamma_x * ym * e_measured
        theta_r_dot = -gamma_r * r * e_measured

        out["r"][k], out["x"][k], out["xm"][k] = r, x, xm
        out["e"][k], out["ym"][k], out["u"][k] = x - xm, ym, u_cmd
        out["theta_x"][k], out["theta_r"][k] = theta_x, theta_r

        x += dt * x_dot
        xm += dt * xm_dot
        theta_x += dt * theta_x_dot
        theta_r += dt * theta_r_dot

        if delay_steps:
            u_buffer[:-1] = u_buffer[1:]

        if not np.isfinite([x, xm, q, q_dot, theta_x, theta_r]).all():
            raise FloatingPointError(f"Numerical divergence in scenario {s.name!r} at t={tk:.3f} s")

    return out


def rms(x: np.ndarray, start_index: int) -> float:
    return float(np.sqrt(np.mean(np.square(x[start_index:]))))


def main() -> None:
    scenarios = [
        Scenario("Nominal"),
        Scenario("Noise", noise_std=0.035),
        Scenario("Unmodeled dynamics", actuator_omega=18.0),
        Scenario("Delay", delay=0.080),
        Scenario("Combined", noise_std=0.035, actuator_omega=18.0, delay=0.080),
    ]
    results = {s.name: simulate(s) for s in scenarios}

    print("Scenario diagnostics over the final 5 seconds")
    for name, z in results.items():
        i0 = int(20.0 / (z["t"][1] - z["t"][0]))
        print(
            f"{name:20s} RMS(e)={rms(z['e'], i0):.5f}, "
            f"max|u|={np.max(np.abs(z['u'])):.3f}, "
            f"theta=({z['theta_x'][-1]:.3f}, {z['theta_r'][-1]:.3f})"
        )

    fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
    for name, z in results.items():
        axes[0].plot(z["t"], z["e"], label=name, linewidth=1.0)
        axes[1].plot(z["t"], z["theta_x"], label=name, linewidth=1.0)
        axes[2].plot(z["t"], z["u"], label=name, linewidth=1.0)
    axes[0].set_ylabel("tracking error")
    axes[1].set_ylabel("theta_x hat")
    axes[2].set_ylabel("command u")
    axes[2].set_xlabel("time (s)")
    axes[0].legend(ncol=2)
    for ax in axes:
        ax.grid(True)
    fig.tight_layout()
    output = Path("Chapter8_Lesson1_python.png")
    fig.savefig(output, dpi=180)
    print(f"Saved {output.resolve()}")
    plt.show()


if __name__ == "__main__":
    main()

Required libraries are numpy and matplotlib. The program simulates all five cases, prints final-window diagnostics, and saves a comparison plot.

Chapter8_Lesson1.cpp

// Chapter8_Lesson1.cpp
// Plain MRAC with noise, neglected actuator dynamics, and input delay.
// Build: g++ -std=c++17 -O2 Chapter8_Lesson1.cpp -o Chapter8_Lesson1
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

struct Scenario {
    std::string name;
    double noise_std{};
    double actuator_omega{}; // <= 0 means ideal actuator
    double actuator_zeta{0.18};
    double delay{};
};

double reference(double t) {
    return (t >= 1.0 ? 1.0 : 0.0) + 0.25 * std::sin(0.7 * t) + 0.15 * std::sin(2.1 * t);
}

void simulate(const Scenario& s, double dt = 1e-3, double tf = 25.0) {
    const double a = 1.35, b = 1.0, am = 2.0, bm = 2.0;
    const double gamma_x = 7.0, gamma_r = 7.0;
    const std::size_t n = static_cast<std::size_t>(std::llround(tf / dt)) + 1;
    const std::size_t delay_steps = static_cast<std::size_t>(std::max(0.0, std::round(s.delay / dt)));
    std::vector<double> delay_line(delay_steps + 1, 0.0);
    std::mt19937_64 rng(7);
    std::normal_distribution<double> normal(0.0, 1.0);

    double x = 0.0, xm = 0.0, q = 0.0, q_dot = 0.0;
    double theta_x = 0.0, theta_r = 0.0;
    double sum_e2 = 0.0;
    std::size_t rms_count = 0;

    std::ofstream csv("Chapter8_Lesson1_" + s.name + ".csv");
    if (!csv) throw std::runtime_error("Cannot open output CSV file");
    csv << "t,r,x,xm,e,ym,u,theta_x,theta_r\n";
    csv << std::setprecision(12);

    for (std::size_t k = 0; k < n; ++k) {
        const double t = k * dt;
        const double r = reference(t);
        const double ym = x + s.noise_std * normal(rng);
        const double e_measured = ym - xm;
        const double u = std::clamp(theta_x * ym + theta_r * r, -25.0, 25.0);

        delay_line.back() = u;
        const double u_delayed = delay_line.front();
        double v = u_delayed;
        if (s.actuator_omega > 0.0) {
            const double w = s.actuator_omega;
            const double q_ddot = w * w * (u_delayed - q) - 2.0 * s.actuator_zeta * w * q_dot;
            q_dot += dt * q_ddot;
            q += dt * q_dot;
            v = q;
        }

        const double x_dot = -a * x + b * v;
        const double xm_dot = -am * xm + bm * r;
        const double theta_x_dot = -gamma_x * ym * e_measured;
        const double theta_r_dot = -gamma_r * r * e_measured;

        csv << t << ',' << r << ',' << x << ',' << xm << ',' << (x - xm) << ','
            << ym << ',' << u << ',' << theta_x << ',' << theta_r << '\n';

        if (t >= 20.0) {
            sum_e2 += (x - xm) * (x - xm);
            ++rms_count;
        }

        x += dt * x_dot;
        xm += dt * xm_dot;
        theta_x += dt * theta_x_dot;
        theta_r += dt * theta_r_dot;

        if (delay_steps > 0) {
            std::rotate(delay_line.begin(), delay_line.begin() + 1, delay_line.end());
        }

        if (!(std::isfinite(x) && std::isfinite(xm) && std::isfinite(theta_x) && std::isfinite(theta_r))) {
            throw std::runtime_error("Numerical divergence in scenario " + s.name);
        }
    }

    std::cout << s.name << ": final theta=(" << theta_x << ", " << theta_r
              << "), final-5s RMS(e)=" << std::sqrt(sum_e2 / static_cast<double>(rms_count)) << '\n';
}

int main() {
    try {
        const std::vector<Scenario> scenarios{
            {"Nominal", 0.0, 0.0, 0.18, 0.0},
            {"Noise", 0.035, 0.0, 0.18, 0.0},
            {"Unmodeled", 0.0, 18.0, 0.18, 0.0},
            {"Delay", 0.0, 0.0, 0.18, 0.080},
            {"Combined", 0.035, 18.0, 0.18, 0.080}
        };
        for (const auto& s : scenarios) simulate(s);
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
    return 0;
}

This implementation uses only the C++17 standard library and writes one CSV file per scenario. Plot the files with any numerical or plotting environment.

Chapter8_Lesson1.java

// Chapter8_Lesson1.java
// Plain MRAC with noise, neglected actuator dynamics, and input delay.
// Build/run: javac Chapter8_Lesson1.java && java Chapter8_Lesson1
import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;
import java.util.Random;

public final class Chapter8_Lesson1 {
    private record Scenario(String name, double noiseStd, double actuatorOmega,
                            double actuatorZeta, double delay) {}

    private static double reference(double t) {
        return (t >= 1.0 ? 1.0 : 0.0)
                + 0.25 * Math.sin(0.7 * t) + 0.15 * Math.sin(2.1 * t);
    }

    private static double clip(double value, double low, double high) {
        return Math.max(low, Math.min(high, value));
    }

    private static void simulate(Scenario s, double dt, double tf) throws IOException {
        final double a = 1.35, b = 1.0, am = 2.0, bm = 2.0;
        final double gammaX = 7.0, gammaR = 7.0;
        final int n = (int) Math.round(tf / dt) + 1;
        final int delaySteps = Math.max(0, (int) Math.round(s.delay() / dt));
        final double[] delayLine = new double[delaySteps + 1];
        final Random random = new Random(7L);

        double x = 0.0, xm = 0.0, q = 0.0, qDot = 0.0;
        double thetaX = 0.0, thetaR = 0.0;
        double sumE2 = 0.0;
        int rmsCount = 0;

        Path output = Path.of("Chapter8_Lesson1_" + s.name() + ".csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output)) {
            writer.write("t,r,x,xm,e,ym,u,theta_x,theta_r\n");
            for (int k = 0; k < n; k++) {
                double t = k * dt;
                double r = reference(t);
                double ym = x + s.noiseStd() * random.nextGaussian();
                double eMeasured = ym - xm;
                double u = clip(thetaX * ym + thetaR * r, -25.0, 25.0);

                delayLine[delayLine.length - 1] = u;
                double uDelayed = delayLine[0];
                double v = uDelayed;
                if (s.actuatorOmega() > 0.0) {
                    double w = s.actuatorOmega();
                    double qDDot = w * w * (uDelayed - q) - 2.0 * s.actuatorZeta() * w * qDot;
                    qDot += dt * qDDot;
                    q += dt * qDot;
                    v = q;
                }

                double xDot = -a * x + b * v;
                double xmDot = -am * xm + bm * r;
                double thetaXDot = -gammaX * ym * eMeasured;
                double thetaRDot = -gammaR * r * eMeasured;

                writer.write(String.format(Locale.ROOT,
                        "%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f%n",
                        t, r, x, xm, x - xm, ym, u, thetaX, thetaR));

                if (t >= 20.0) {
                    sumE2 += (x - xm) * (x - xm);
                    rmsCount++;
                }

                x += dt * xDot;
                xm += dt * xmDot;
                thetaX += dt * thetaXDot;
                thetaR += dt * thetaRDot;

                if (delaySteps > 0) {
                    System.arraycopy(delayLine, 1, delayLine, 0, delayLine.length - 1);
                }

                if (!(Double.isFinite(x) && Double.isFinite(xm)
                        && Double.isFinite(thetaX) && Double.isFinite(thetaR))) {
                    throw new ArithmeticException("Numerical divergence in " + s.name());
                }
            }
        }

        System.out.printf(Locale.ROOT,
                "%s: final theta=(%.5f, %.5f), final-5s RMS(e)=%.6f%n",
                s.name(), thetaX, thetaR, Math.sqrt(sumE2 / rmsCount));
    }

    public static void main(String[] args) {
        Scenario[] scenarios = {
                new Scenario("Nominal", 0.0, 0.0, 0.18, 0.0),
                new Scenario("Noise", 0.035, 0.0, 0.18, 0.0),
                new Scenario("Unmodeled", 0.0, 18.0, 0.18, 0.0),
                new Scenario("Delay", 0.0, 0.0, 0.18, 0.080),
                new Scenario("Combined", 0.035, 18.0, 0.18, 0.080)
        };
        try {
            for (Scenario s : scenarios) simulate(s, 1e-3, 25.0);
        } catch (IOException | ArithmeticException ex) {
            System.err.println("Error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

The Java program uses the standard library, deterministic random seeds, and CSV output. It requires a Java version supporting records; replace the record with a conventional class on older Java versions.

Chapter8_Lesson1.m

% Chapter8_Lesson1.m
% Plain MRAC under noise, neglected actuator dynamics, and input delay.
% The script also contains an optional Simulink model builder.
clear; close all; clc;

scenarios = struct( ...
    'name', {'Nominal','Noise','Unmodeled dynamics','Delay','Combined'}, ...
    'noiseStd', {0, 0.035, 0, 0, 0.035}, ...
    'omega', {0, 0, 18, 0, 18}, ...
    'zeta', {0.18, 0.18, 0.18, 0.18, 0.18}, ...
    'delay', {0, 0, 0, 0.080, 0.080});

results = cell(size(scenarios));
for i = 1:numel(scenarios)
    results{i} = simulateMRAC(scenarios(i), 1e-3, 25, 7);
    z = results{i};
    idx = z.t >= 20;
    fprintf('%-20s RMS(e)=%.6f, max|u|=%.4f, theta=(%.4f, %.4f)\n', ...
        scenarios(i).name, rms(z.e(idx)), max(abs(z.u)), z.thetaX(end), z.thetaR(end));
end

figure('Name','Chapter 8 Lesson 1: practical MRAC issues');
tiledlayout(3,1);
nexttile; hold on; grid on;
for i = 1:numel(results), plot(results{i}.t, results{i}.e, 'DisplayName', scenarios(i).name); end
ylabel('tracking error'); legend('Location','best');
nexttile; hold on; grid on;
for i = 1:numel(results), plot(results{i}.t, results{i}.thetaX, 'DisplayName', scenarios(i).name); end
ylabel('theta_x hat');
nexttile; hold on; grid on;
for i = 1:numel(results), plot(results{i}.t, results{i}.u, 'DisplayName', scenarios(i).name); end
ylabel('command u'); xlabel('time (s)');
exportgraphics(gcf, 'Chapter8_Lesson1_matlab.png', 'Resolution', 180);

% Set true only when Simulink is installed and a block-diagram skeleton is desired.
buildSimulinkSkeleton = false;
if buildSimulinkSkeleton
    buildChapter8Lesson1SimulinkModel;
end

function z = simulateMRAC(s, dt, tf, seed)
    rng(seed);
    a = 1.35; b = 1.0; am = 2.0; bm = 2.0;
    gammaX = 7.0; gammaR = 7.0;
    t = (0:dt:tf)'; n = numel(t);
    delaySteps = max(0, round(s.delay / dt));
    delayLine = zeros(delaySteps + 1, 1);
    x = 0; xm = 0; q = 0; qDot = 0; thetaX = 0; thetaR = 0;
    z.t = t;
    fields = {'r','x','xm','e','ym','u','thetaX','thetaR'};
    for f = 1:numel(fields), z.(fields{f}) = zeros(n,1); end

    for k = 1:n
        tk = t(k);
        r = double(tk >= 1) + 0.25*sin(0.7*tk) + 0.15*sin(2.1*tk);
        ym = x + s.noiseStd*randn;
        eMeasured = ym - xm;
        u = min(25, max(-25, thetaX*ym + thetaR*r));
        delayLine(end) = u;
        uDelayed = delayLine(1);

        if s.omega > 0
            qDDot = s.omega^2*(uDelayed - q) - 2*s.zeta*s.omega*qDot;
            qDot = qDot + dt*qDDot;
            q = q + dt*qDot;
            v = q;
        else
            v = uDelayed;
        end

        xDot = -a*x + b*v;
        xmDot = -am*xm + bm*r;
        thetaXDot = -gammaX*ym*eMeasured;
        thetaRDot = -gammaR*r*eMeasured;

        z.r(k)=r; z.x(k)=x; z.xm(k)=xm; z.e(k)=x-xm; z.ym(k)=ym;
        z.u(k)=u; z.thetaX(k)=thetaX; z.thetaR(k)=thetaR;

        x = x + dt*xDot;
        xm = xm + dt*xmDot;
        thetaX = thetaX + dt*thetaXDot;
        thetaR = thetaR + dt*thetaRDot;
        if delaySteps > 0, delayLine(1:end-1) = delayLine(2:end); end
        if any(~isfinite([x xm q qDot thetaX thetaR]))
            error('Numerical divergence in scenario %s at t=%.3f', s.name, tk);
        end
    end
end

function buildChapter8Lesson1SimulinkModel
    model = 'Chapter8_Lesson1_Simulink';
    if bdIsLoaded(model), close_system(model, 0); end
    new_system(model); open_system(model);
    set_param(model, 'StopTime', '25', 'Solver', 'ode4', 'FixedStep', '0.001');

    add_block('simulink/Sources/Step', [model '/Reference'], ...
        'Time', '1', 'Before', '0', 'After', '1', 'Position', [30 55 60 85]);
    add_block('simulink/Continuous/State-Space', [model '/Reference Model'], ...
        'A', '-2', 'B', '2', 'C', '1', 'D', '0', 'Position', [120 45 230 95]);
    add_block('simulink/Math Operations/Sum', [model '/Tracking Error'], ...
        'Inputs', '+-', 'Position', [300 120 325 160]);
    add_block('simulink/Continuous/Transport Delay', [model '/Input Delay'], ...
        'DelayTime', '0.08', 'Position', [500 225 590 255]);
    add_block('simulink/Continuous/Transfer Fcn', [model '/Neglected Actuator'], ...
        'Numerator', '[18^2]', 'Denominator', '[1 2*0.18*18 18^2]', ...
        'Position', [630 215 760 265]);
    add_block('simulink/Continuous/State-Space', [model '/Plant'], ...
        'A', '-1.35', 'B', '1', 'C', '1', 'D', '0', 'Position', [810 210 920 270]);
    add_block('simulink/Sinks/Scope', [model '/Scope'], ...
        'NumInputPorts', '3', 'Position', [990 90 1030 180]);

    % The adaptive controller is intentionally represented as a subsystem so
    % students can implement theta-dot = -Gamma*[y;r]*(y-xm) with Product,
    % Gain, Integrator, and Sum blocks during the laboratory exercise.
    add_block('simulink/Ports & Subsystems/Subsystem', [model '/Plain MRAC Controller'], ...
        'Position', [350 205 455 280]);
    delete_line([model '/Plain MRAC Controller'], 'In1/1', 'Out1/1');
    delete_block([model '/Plain MRAC Controller/In1']);
    delete_block([model '/Plain MRAC Controller/Out1']);
    add_block('simulink/Ports & Subsystems/In1', [model '/Plain MRAC Controller/y'], ...
        'Port', '1', 'Position', [25 38 55 52]);
    add_block('simulink/Ports & Subsystems/In1', [model '/Plain MRAC Controller/r'], ...
        'Port', '2', 'Position', [25 88 55 102]);
    add_block('simulink/Ports & Subsystems/In1', [model '/Plain MRAC Controller/e'], ...
        'Port', '3', 'Position', [25 138 55 152]);
    add_block('simulink/Ports & Subsystems/Out1', [model '/Plain MRAC Controller/u'], ...
        'Port', '1', 'Position', [310 83 340 97]);
    add_block('simulink/Signal Routing/Mux', [model '/Plain MRAC Controller/Regressor'], ...
        'Inputs', '2', 'Position', [90 45 95 105]);
    add_block('simulink/Math Operations/Product', [model '/Plain MRAC Controller/phi_e'], ...
        'Multiplication', 'Element-wise(.*)', 'Position', [125 85 165 125]);
    add_block('simulink/Math Operations/Gain', [model '/Plain MRAC Controller/-Gamma'], ...
        'Gain', '-diag([7 7])', 'Multiplication', 'Matrix(K*u)', 'Position', [190 85 240 125]);
    add_block('simulink/Continuous/Integrator', [model '/Plain MRAC Controller/theta'], ...
        'InitialCondition', '[0;0]', 'Position', [260 85 290 115]);
    add_block('simulink/Math Operations/Dot Product', [model '/Plain MRAC Controller/thetaT_phi'], ...
        'Position', [245 30 280 60]);

    add_line([model '/Plain MRAC Controller'], 'y/1', 'Regressor/1');
    add_line([model '/Plain MRAC Controller'], 'r/1', 'Regressor/2');
    add_line([model '/Plain MRAC Controller'], 'Regressor/1', 'phi_e/1');
    add_line([model '/Plain MRAC Controller'], 'e/1', 'phi_e/2');
    add_line([model '/Plain MRAC Controller'], 'phi_e/1', '-Gamma/1');
    add_line([model '/Plain MRAC Controller'], '-Gamma/1', 'theta/1');
    add_line([model '/Plain MRAC Controller'], 'theta/1', 'thetaT_phi/1');
    add_line([model '/Plain MRAC Controller'], 'Regressor/1', 'thetaT_phi/2');
    add_line([model '/Plain MRAC Controller'], 'thetaT_phi/1', 'u/1');

    add_line(model, 'Reference/1', 'Reference Model/1');
    add_line(model, 'Reference/1', 'Plain MRAC Controller/2');
    add_line(model, 'Reference Model/1', 'Tracking Error/2');
    add_line(model, 'Plant/1', 'Tracking Error/1');
    add_line(model, 'Plant/1', 'Plain MRAC Controller/1');
    add_line(model, 'Tracking Error/1', 'Plain MRAC Controller/3');
    add_line(model, 'Plain MRAC Controller/1', 'Input Delay/1');
    add_line(model, 'Input Delay/1', 'Neglected Actuator/1');
    add_line(model, 'Neglected Actuator/1', 'Plant/1');
    add_line(model, 'Plant/1', 'Scope/1');
    add_line(model, 'Reference Model/1', 'Scope/2');
    add_line(model, 'Tracking Error/1', 'Scope/3');

    save_system(model);
    fprintf('Created %s.slx. Add a noise source at the measured-output path to compare cases.\n', model);
end

The script runs the numerical experiment directly in MATLAB and includes an optional programmatic Simulink skeleton. Set buildSimulinkSkeleton = true when Simulink is installed. The generated subsystem implements the plain update law from Product, Gain, Integrator, Mux, and Dot Product blocks; a Transport Delay and neglected actuator transfer function are placed in the input path.

Chapter8_Lesson1.nb

Notebook[{
 Cell["Chapter 8, Lesson 1 — Practical Issues in Plain MRAC", "Title"],
 Cell["This notebook simulates nominal MRAC and cases with measurement noise, neglected actuator dynamics, and input delay. No robust modification is used.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];
dt = 0.001; tf = 25.;
a = 1.35; b = 1.; am = 2.; bm = 2.; gammaX = 7.; gammaR = 7.;
reference[t_] := Boole[t >= 1.] + 0.25 Sin[0.7 t] + 0.15 Sin[2.1 t];
scenarios = {
 <|\"Name\" -> \"Nominal\", \"Noise\" -> 0., \"Omega\" -> 0., \"Zeta\" -> 0.18, \"Delay\" -> 0.|>,
 <|\"Name\" -> \"Noise\", \"Noise\" -> 0.035, \"Omega\" -> 0., \"Zeta\" -> 0.18, \"Delay\" -> 0.|>,
 <|\"Name\" -> \"Unmodeled dynamics\", \"Noise\" -> 0., \"Omega\" -> 18., \"Zeta\" -> 0.18, \"Delay\" -> 0.|>,
 <|\"Name\" -> \"Delay\", \"Noise\" -> 0., \"Omega\" -> 0., \"Zeta\" -> 0.18, \"Delay\" -> 0.080|>,
 <|\"Name\" -> \"Combined\", \"Noise\" -> 0.035, \"Omega\" -> 18., \"Zeta\" -> 0.18, \"Delay\" -> 0.080|>
};

simulate[s_Association, seed_: 7] := Module[
 {n = Round[tf/dt] + 1, delaySteps, buffer, x = 0., xm = 0., q = 0., qd = 0.,
  thetaX = 0., thetaR = 0., rows = {}, t, r, ym, eMeasured, u, ud, v, qdd},
 SeedRandom[seed];
 delaySteps = Max[0, Round[s[\"Delay\"]/dt]]; buffer = ConstantArray[0., delaySteps + 1];
 Do[
  t = (k - 1) dt; r = reference[t]; ym = x + s[\"Noise\"] RandomVariate[NormalDistribution[]];
  eMeasured = ym - xm; u = Clip[thetaX ym + thetaR r, {-25., 25.}];
  buffer[[-1]] = u; ud = First[buffer];
  If[s[\"Omega\"] > 0.,
   qdd = s[\"Omega\"]^2 (ud - q) - 2 s[\"Zeta\"] s[\"Omega\"] qd;
   qd += dt qdd; q += dt qd; v = q,
   v = ud
  ];
  AppendTo[rows, {t, r, x, xm, x - xm, ym, u, thetaX, thetaR}];
  x += dt (-a x + b v); xm += dt (-am xm + bm r);
  thetaX += dt (-gammaX ym eMeasured); thetaR += dt (-gammaR r eMeasured);
  If[delaySteps > 0, buffer = RotateLeft[buffer]],
  {k, 1, n}];
 Association[\"Name\" -> s[\"Name\"], \"Data\" -> rows]
];

results = simulate /@ scenarios;
headers = {\"t\", \"r\", \"x\", \"xm\", \"e\", \"ym\", \"u\", \"thetaX\", \"thetaR\"};
Do[Export[\"Chapter8_Lesson1_\" <> StringReplace[result[\"Name\"], \" \" -> \"_\"] <> \".csv\", Prepend[result[\"Data\"], headers]], {result, results}];

errorPlot = ListLinePlot[Table[result[\"Data\"][[All, {1, 5}]], {result, results}],
 PlotLegends -> (#[\"Name\"] & /@ results), Frame -> True,
 FrameLabel -> {\"time (s)\", \"tracking error\"}, PlotRange -> All, ImageSize -> Large];
parameterPlot = ListLinePlot[Table[result[\"Data\"][[All, {1, 8}]], {result, results}],
 PlotLegends -> (#[\"Name\"] & /@ results), Frame -> True,
 FrameLabel -> {\"time (s)\", \"theta_x hat\"}, PlotRange -> All, ImageSize -> Large];
Column[{errorPlot, parameterPlot}]"
 ], "Input"]
}, WindowTitle -> "Chapter8_Lesson1", Saveable -> True]

The notebook stores each trajectory as a list, exports scenario CSV files, and plots tracking error and parameter evolution. The fixed-step update mirrors the other languages so results can be compared directly.

Before deploying an adaptive controller, perform a robustness-oriented test matrix rather than a single nominal simulation. Vary sensor noise amplitude and bandwidth, actuator-mode frequency and damping, input delay, sampling period, adaptation gain, reference-model bandwidth, initial parameter estimates, command richness, and saturation level.

Record at least: RMS and peak tracking error; RMS and peak control effort; parameter norm \( \|\hat\theta\| \); parameter-rate norm \( \|\dot{\hat\theta}\| \); spectral energy near neglected modes; saturation occupancy; and recovery after commands or disturbances. A robust design should be assessed over uncertainty ranges and initial conditions, not only at one parameter point.

Use progressively higher-fidelity models. Begin with the nominal reduced-order model, add isolated imperfections one at a time, combine them, then perform software-in-the-loop and hardware-in-the-loop tests. Adaptation should be disabled or safely frozen when signal validity, timing, or actuator authority falls outside certified limits.

Problem 1 (Nominal cancellation): For the scalar MRAC benchmark, derive the ideal parameters and prove that the selected update law yields \( \dot V=-a_m e^2 \).

Solution: Matching the plant closed loop to the reference model requires

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

so \( \theta_x^{*}=(a-a_m)/b \) and \( \theta_r^{*}=b_m/b \). Substitution gives \( \dot e=-a_m e+b\tilde\theta^T\phi \). With

\[ V=\frac12e^2+\frac b2\tilde\theta^T\Gamma^{-1}\tilde\theta, \qquad \dot{\hat\theta}=-\Gamma\phi e, \]

the two cross terms \( be\tilde\theta^T\phi \) cancel exactly, leaving \( \dot V=-a_m e^2 \).

Problem 2 (Ultimate-bound estimate): Assume the perturbed error equation satisfies \( |\Delta|\le0.08+0.25|e| \) and \( a_m=2 \). Find the error radius outside which the Lyapunov derivative is negative.

Solution:

\[ \dot V\le-(2-0.25)e^2+0.08|e|=-1.75e^2+0.08|e|. \]

For nonzero \( e \), negativity is guaranteed when

\[ |e|>\frac{0.08}{1.75}\approx0.0457. \]

This is a conditional tracking-error estimate; it does not by itself prove bounded adaptive parameters.

Problem 3 (Noise-induced bias): Let \( y=x+n \), \( e_y=e+n \), with zero-mean noise independent of \( x \) and \( e \) and variance \( \sigma_n^2 \). Compute \( \mathbb E[y e_y] \).

Solution:

\[ \begin{aligned} \mathbb E[y e_y] &=\mathbb E[(x+n)(e+n)]\\ &=xe+x\mathbb E[n]+e\mathbb E[n]+\mathbb E[n^2]\\ &=xe+\sigma_n^2. \end{aligned} \]

Therefore the state-feedback parameter update has a nonzero mean component even when the noise itself has zero mean.

Problem 4 (Delay perturbation): Show that \( |u(t-\tau)-u(t)|\le\tau L_u \) if \( |\dot u(s)|\le L_u \) over \( [t-\tau,t] \).

Solution: By the fundamental theorem of calculus,

\[ u(t)-u(t-\tau)=\int_{t-\tau}^{t}\dot u(s)\,ds. \]

Taking absolute values and applying the derivative bound gives

\[ |u(t-\tau)-u(t)|\le\int_{t-\tau}^{t}|\dot u(s)|\,ds\le\tau L_u. \]

Hence the plant perturbation obeys \( |\Delta_\tau|\le|b|\tau L_u \).

Problem 5 (Parasitic actuator mismatch): For \( H_f(s)=\omega_f^2/(s^2+2\zeta_f\omega_fs+\omega_f^2) \), explain why increasing adaptation gain can worsen behavior near \( \omega_f \).

Solution: A larger adaptation gain increases the bandwidth and amplitude of parameter motion. Since the control command contains \( \dot{\hat\theta}^T\phi \) through its time derivative, more command energy reaches frequencies near the neglected mode. Near \( \omega_f \), the actuator has significant phase lag and may have magnitude peaking for small damping. The assumed relation \( q\approx u \) then fails, producing a phase-lagged error that drives further adaptation. Thus the adaptive loop and parasitic dynamics can form a positive feedback mechanism.

Problem 6 (Numerical experiment): Run the supplied code with \( \Gamma=\operatorname{diag}(2,2) \), \( \operatorname{diag}(7,7) \), and \( \operatorname{diag}(20,20) \). Compare tracking error, parameter variation, and control bandwidth.

Solution: The exact values depend on integration step and random realization. The expected trend is that moderate gain improves nominal transient adaptation, while high gain increases sensitivity to noise, delay, and the neglected actuator mode. A correct report must include all three metrics; selecting the gain only from nominal RMS error is insufficient.

Plain MRAC obtains nominal stability through exact Lyapunov cross-term cancellation. Measurement noise corrupts both the tracking error and regressor and can create biased parameter updates. Neglected actuator or structural modes introduce state-dependent mismatches whose phase lag can be amplified by aggressive adaptation. Input delay produces a perturbation proportional to recent control variation and removes phase margin. These effects interact, so practical evaluation must monitor parameter motion and control bandwidth in addition to tracking error. The next lessons introduce modifications that trade exact nominal adaptation for boundedness and improved robustness.