After completing this lesson, students should be able to:

• distinguish parametric uncertainty from unmodeled dynamics and disturbances;
• derive the nominal MRAC error model and its Lyapunov stability proof;
• identify the perturbation term that invalidates the nominal cancellation;
• explain parameter drift as a consequence of an undamped adaptive integrator;
• relate adaptation gain, closed-loop bandwidth, and neglected high-frequency modes;
• simulate sensitivity to a neglected stable actuator pole in five programming environments.

This lesson concerns diagnosis. Robust modifications such as leakage, projection, dead zones, and the sigma modification were introduced earlier and will be analyzed systematically in Lesson 2 of this chapter.

A linearly parameterized adaptive design assumes that all relevant uncertainty can be represented by an unknown constant vector \( \theta^* \) multiplying a known regressor \( \phi(t) \). The actual plant is more accurately written as

\[ y = G_0(s,\theta^*)u + W_\Delta(s)u + d, \]

where \(G_0\) is the modeled family, \(W_\Delta\) represents neglected stable modes, delays, sensor filters, flexible dynamics, or actuator dynamics, and \(d\) is an exogenous disturbance. The crucial distinction is that changing \(\theta\) cannot generally reproduce \(W_\Delta(s)\) over all frequencies.

flowchart TD
  R["Reference r"] --> C["Adaptive controller"]
  C --> U["Command u"]
  U --> A["Unmodeled actuator / fast modes"]
  A --> P["Modeled plant dynamics"]
  P --> Y["Measured output y"]
  Y --> E["Tracking error e"]
  E --> AD["Parameter update integrator"]
  AD --> C
  D["Disturbance and noise"] --> Y
  A -. "not represented by theta" .-> E
        

Even when the neglected block is internally stable, its phase lag and gain at frequencies reached by the adaptive loop may change negative feedback into weak damping or positive feedback.

Let the nominal transfer function be \(G_0(s)\). Two common uncertainty descriptions are

\[ G_p(s)=G_0(s)+W_a(s)\Delta_a(s), \qquad \|\Delta_a\|_\infty \leq 1, \]

\[ G_p(s)=G_0(s)\left[1+W_m(s)\Delta_m(s)\right], \qquad \|\Delta_m\|_\infty \leq 1. \]

In robust-control language, \(W_a\) and \(W_m\) describe the frequency distribution of model error. Here they are not used to synthesize an H-infinity controller; they are used to reveal why a low-frequency adaptive model may be accurate while the same model becomes unreliable at higher frequencies.

Consider a modeled first-order plant and an actual plant with an actuator lag:

\[ G_0(s)=\frac{b}{s+a}, \qquad G_p(s)=\frac{b}{(s+a)(\tau s+1)}. \]

The induced multiplicative error is

\[ \Delta_m(s)=\frac{G_p(s)}{G_0(s)}-1 =-\frac{\tau s}{\tau s+1}. \]

Thus \(|\Delta_m(j\omega)|\approx \tau\omega\) for \(\omega\tau \ll 1\), but \(|\Delta_m(j\omega)| → 1\) as frequency grows. A mode that appears negligible in the intended control band is therefore not negligible to an adaptive law that creates substantial high-frequency activity.

Consider the ideal plant

\[ \dot{x}=-a x+b u, \qquad a > 0,\quad b > 0, \]

and reference model

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

Choose the direct adaptive control law

\[ u=\theta_x x+\theta_r r. \]

Ideal matching parameters exist when the plant is exactly first order:

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

Define tracking and parameter errors \(e=x-x_m\) and \(\widetilde{\theta}=\theta-\theta^*\), with \(\phi=[x\; r]^T\). Then

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

Use the gradient law

\[ \dot{\theta}=-\gamma b\phi e, \qquad \gamma > 0. \]

For the Lyapunov candidate

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

its derivative is

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

The proof depends on exact matching and exact cancellation. It establishes boundedness of \(e\) and \(\widetilde{\theta}\), but it does not provide a robustness margin for dynamics that are absent from the error equation.

Suppose neglected dynamics and disturbances generate an additional signal \(\delta(t)\) in the tracking-error equation:

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

With the same adaptive law and Lyapunov function,

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

Young's inequality, \(pq\leq \frac{\varepsilon}{2}p^2+ \frac{1}{2\varepsilon}q^2\), gives

\[ \dot{V} \leq -\frac{a_m}{2}e^2+\frac{1}{2a_m}\delta^2(t). \]

This inequality can bound the tracking error under additional assumptions, but it contains no negative term in \(\|\widetilde{\theta}\|^2\). Therefore, a small persistent error can continue integrating into the parameters. In particular,

\[ \theta(t)=\theta(0)-\gamma b\int_0^t \phi(\sigma)e(\sigma)\,d\sigma. \]

If the integral has a nonzero average component, the parameter vector can drift even when \(e(t)\) remains numerically small for a long time. This is the central robustness defect of an unmodified gradient law.

flowchart TD
  A["Small model error or sensor disturbance"] --> B["Persistent tracking error component"]
  B --> C["Adaptive integrator changes parameters"]
  C --> D["Controller gain and bandwidth increase"]
  D --> E["Neglected pole adds phase lag"]
  E --> F["More oscillatory control and output"]
  F --> B
  C --> G["Parameter drift"]
  D --> H["Control saturation or numerical stiffness"]
        

Several mechanisms can coexist:

1. Parameter-drift mechanism: disturbance-induced error is integrated by the adaptation law without parameter damping.
2. Fast-adaptation mechanism: increasing \(\gamma\) accelerates parameter motion and can inject energy near neglected poles.
3. High-frequency mechanism: differentiation, sensor noise, actuator lag, or flexible modes introduce phase loss where the adaptive loop has high effective gain.
4. Throughput mechanism: an exogenous disturbance propagates through the plant, tracking error, and parameter estimator repeatedly.
5. Linearized instability: around a slowly varying operating point, the frozen adaptive closed loop can possess poorly damped or unstable local modes despite nominal nonlinear stability claims for the ideal model.

The phrase “stable neglected dynamics” is therefore insufficient. Robustness is a property of the interconnection, not of the neglected subsystem alone.

The actual plant is

\[ \dot{x}=-a x+b z+d(t), \qquad \tau\dot{z}=-z+u, \qquad \tau > 0. \]

The controller was designed under the approximation \(z\approx u\). Let the actuator mismatch be \(\eta=z-u\). The tracking-error equation becomes

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

Hence the perturbation in Section 5 is

\[ \delta(t)=b\eta(t)+d(t). \]

From \(\tau\dot{z}=-z+u\), the mismatch satisfies

\[ \tau\dot{\eta}=-\eta-\tau\dot{u}. \]

Therefore rapid controller motion makes \(\eta\) large. Since \(u=\theta^T\phi\),

\[ \dot{u}=\dot{\theta}^{T}\phi+\theta^{T}\dot{\phi} =-\gamma b e\,\phi^{T}\phi+\theta^{T}\dot{\phi}. \]

The neglected-mode forcing therefore grows directly with the adaptation gain, tracking error, and regressor energy. This equation gives a concrete reason why “larger adaptation gain” is not equivalent to “better tracking.”

Singular-perturbation intuition suggests that the actuator approximation is reasonable only when the command varies slowly relative to \(1/\tau\). From the stable filter relation,

\[ \eta(t)=-\int_0^t e^{-(t-\sigma)/\tau}\dot{u}(\sigma)\,d\sigma +e^{-t/\tau}\eta(0), \]

so that

\[ |\eta(t)| \leq e^{-t/\tau}|\eta(0)| +\tau\sup_{0\leq \sigma\leq t}|\dot{u}(\sigma)|. \]

A useful engineering condition is consequently

\[ \tau\,\|\dot{u}\|_\infty \ll \|u\|_\infty. \]

Because \(\dot{u}\) depends on \(\gamma\), this condition couples model bandwidth and adaptation speed. The limit \(\tau → 0\) is not uniform over arbitrarily large adaptation gains.

A credible robustness study should include at least the following checks:

• sweep the adaptation gain over at least one decade;
• insert realistic actuator and sensor poles rather than ideal algebraic blocks;
• test delays, quantization, saturation, and sample-and-hold effects;
• inject low-amplitude broadband or high-frequency disturbances;
• monitor parameter norms, control-rate magnitude, and saturation duty cycle;
• repeat simulations with smaller integration steps to separate physical instability from numerical instability;
• compare commanded input \(u\) with delivered input \(z\).

Tracking error alone is not a sufficient diagnostic. A system can track acceptably while parameters drift toward unsafe values or while control effort approaches actuator limits.

The program uses a fixed-step fourth-order Runge–Kutta integrator. It compares the ideal first-order plant with the same plant followed by an unmodeled actuator pole of time constant \(\tau=0.08\) s. A small high-frequency disturbance is activated after 12 s.

Chapter22_Lesson1.py

"""Chapter22_Lesson1.py

Sensitivity of a first-order MRAC loop to an unmodeled actuator pole.
The same gradient adaptive law is simulated for a nominal plant and for a
plant containing a neglected stable first-order actuator dynamic.
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    a: float = 1.0
    b: float = 1.0
    a_m: float = 2.0
    b_m: float = 2.0
    gamma: float = 30.0
    tau: float = 0.0
    dt: float = 1.0e-3
    t_final: float = 30.0
    u_limit: float = 50.0


def reference(t: float) -> float:
    return 1.0 + 0.35 * np.sin(0.6 * t)


def disturbance(t: float) -> float:
    # A small high-frequency output disturbance is introduced after 12 s.
    return 0.0 if t < 12.0 else 0.05 * np.sin(18.0 * t)


def rhs(t: float, y: np.ndarray, cfg: Config) -> np.ndarray:
    x, z, x_m, theta_x, theta_r = y
    r = reference(t)
    e = x - x_m

    u = float(np.clip(theta_x * x + theta_r * r, -cfg.u_limit, cfg.u_limit))

    if cfg.tau <= 0.0:
        plant_input = u
        z_dot = 0.0
    else:
        plant_input = z
        z_dot = (-z + u) / cfg.tau

    x_dot = -cfg.a * x + cfg.b * plant_input + disturbance(t)
    x_m_dot = -cfg.a_m * x_m + cfg.b_m * r

    # Standard gradient/MRAC update for known positive high-frequency gain.
    theta_x_dot = -cfg.gamma * x * e
    theta_r_dot = -cfg.gamma * r * e

    return np.array([x_dot, z_dot, x_m_dot, theta_x_dot, theta_r_dot], dtype=float)


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


def simulate(cfg: Config) -> dict[str, np.ndarray]:
    count = int(round(cfg.t_final / cfg.dt)) + 1
    t = np.linspace(0.0, cfg.t_final, count)
    y = np.zeros((count, 5), dtype=float)
    u = np.zeros(count, dtype=float)
    r = np.zeros(count, dtype=float)

    for k in range(count - 1):
        x, _, _, theta_x, theta_r = y[k]
        r[k] = reference(t[k])
        u[k] = np.clip(theta_x * x + theta_r * r[k], -cfg.u_limit, cfg.u_limit)
        y[k + 1] = rk4_step(rhs, t[k], y[k], cfg.dt, cfg)

        if not np.all(np.isfinite(y[k + 1])):
            raise FloatingPointError(f"Non-finite state at t={t[k + 1]:.6f} s")

    r[-1] = reference(t[-1])
    u[-1] = np.clip(y[-1, 3] * y[-1, 0] + y[-1, 4] * r[-1], -cfg.u_limit, cfg.u_limit)

    return {
        "t": t,
        "x": y[:, 0],
        "z": y[:, 1],
        "x_m": y[:, 2],
        "theta_x": y[:, 3],
        "theta_r": y[:, 4],
        "u": u,
        "r": r,
        "e": y[:, 0] - y[:, 2],
    }


def summarize(name: str, data: dict[str, np.ndarray]) -> None:
    rms_error = float(np.sqrt(np.mean(data["e"] ** 2)))
    max_error = float(np.max(np.abs(data["e"])))
    max_parameter = float(
        max(np.max(np.abs(data["theta_x"])), np.max(np.abs(data["theta_r"])))
    )
    print(
        f"{name:20s} RMS(e)={rms_error:8.4f}, "
        f"max|e|={max_error:8.4f}, max|theta|={max_parameter:8.4f}"
    )


def main() -> None:
    nominal = simulate(Config(tau=0.0))
    unmodeled = simulate(Config(tau=0.08))

    summarize("Nominal model", nominal)
    summarize("Unmodeled pole", unmodeled)

    output = Path("Chapter22_Lesson1_results.csv")
    np.savetxt(
        output,
        np.column_stack(
            [
                nominal["t"],
                nominal["x"],
                nominal["x_m"],
                nominal["e"],
                nominal["theta_x"],
                nominal["theta_r"],
                unmodeled["x"],
                unmodeled["x_m"],
                unmodeled["e"],
                unmodeled["theta_x"],
                unmodeled["theta_r"],
            ]
        ),
        delimiter=",",
        header=(
            "t,x_nominal,xm_nominal,e_nominal,theta_x_nominal,theta_r_nominal,"
            "x_unmodeled,xm_unmodeled,e_unmodeled,theta_x_unmodeled,theta_r_unmodeled"
        ),
        comments="",
    )

    plt.figure(figsize=(10, 6))
    plt.plot(nominal["t"], nominal["e"], label="nominal model")
    plt.plot(unmodeled["t"], unmodeled["e"], label="unmodeled actuator pole")
    plt.xlabel("time [s]")
    plt.ylabel("tracking error e")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter22_Lesson1_tracking_error.png", dpi=180)

    plt.figure(figsize=(10, 6))
    plt.plot(unmodeled["t"], unmodeled["theta_x"], label="theta_x")
    plt.plot(unmodeled["t"], unmodeled["theta_r"], label="theta_r")
    plt.xlabel("time [s]")
    plt.ylabel("adaptive parameters")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter22_Lesson1_parameters.png", dpi=180)

    print(f"Wrote {output}")


if __name__ == "__main__":
    main()

Required packages are numpy and matplotlib. The script writes a CSV file and two figures so that tracking error and parameter motion can be inspected separately.

Chapter22_Lesson1.cpp

// Chapter22_Lesson1.cpp
// Sensitivity of MRAC to an unmodeled actuator pole.

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

struct Config {
    double a = 1.0;
    double b = 1.0;
    double am = 2.0;
    double bm = 2.0;
    double gamma = 30.0;
    double tau = 0.0;
    double dt = 1.0e-3;
    double tf = 30.0;
    double uLimit = 50.0;
};

using State = std::array<double, 5>; // x, z, xm, theta_x, theta_r

struct Metrics {
    double sumErrorSquared = 0.0;
    double maxAbsError = 0.0;
    double maxAbsParameter = 0.0;
    std::size_t samples = 0;
};

double reference(double t) {
    return 1.0 + 0.35 * std::sin(0.6 * t);
}

double disturbance(double t) {
    return t < 12.0 ? 0.0 : 0.05 * std::sin(18.0 * t);
}

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

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

State rhs(double t, const State& y, const Config& cfg) {
    const double x = y[0];
    const double z = y[1];
    const double xm = y[2];
    const double thetaX = y[3];
    const double thetaR = y[4];

    const double r = reference(t);
    const double e = x - xm;
    const double u = clamp(thetaX * x + thetaR * r, cfg.uLimit);

    const double plantInput = cfg.tau <= 0.0 ? u : z;
    const double zDot = cfg.tau <= 0.0 ? 0.0 : (-z + u) / cfg.tau;

    return State{
        -cfg.a * x + cfg.b * plantInput + disturbance(t),
        zDot,
        -cfg.am * xm + cfg.bm * r,
        -cfg.gamma * x * e,
        -cfg.gamma * r * e
    };
}

State rk4(double t, const State& y, const Config& cfg) {
    const double h = cfg.dt;
    const State k1 = rhs(t, y, cfg);
    const State k2 = rhs(t + 0.5 * h, addScaled(y, k1, 0.5 * h), cfg);
    const State k3 = rhs(t + 0.5 * h, addScaled(y, k2, 0.5 * h), cfg);
    const State k4 = rhs(t + h, addScaled(y, k3, h), cfg);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = y[i] + h * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Non-finite state encountered");
        }
    }
    return next;
}

Metrics simulate(const Config& cfg, const std::string& filename) {
    std::ofstream out(filename);
    if (!out) {
        throw std::runtime_error("Cannot open output file: " + filename);
    }
    out << "t,x,z,xm,e,theta_x,theta_r,u\n";
    out << std::setprecision(12);

    State y{};
    Metrics metrics{};
    const std::size_t steps = static_cast<std::size_t>(std::llround(cfg.tf / cfg.dt));

    for (std::size_t k = 0; k <= steps; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        const double r = reference(t);
        const double e = y[0] - y[2];
        const double u = clamp(y[3] * y[0] + y[4] * r, cfg.uLimit);

        out << t << ',' << y[0] << ',' << y[1] << ',' << y[2] << ',' << e << ','
            << y[3] << ',' << y[4] << ',' << u << '\n';

        metrics.sumErrorSquared += e * e;
        metrics.maxAbsError = std::max(metrics.maxAbsError, std::abs(e));
        metrics.maxAbsParameter = std::max(
            metrics.maxAbsParameter,
            std::max(std::abs(y[3]), std::abs(y[4]))
        );
        ++metrics.samples;

        if (k < steps) {
            y = rk4(t, y, cfg);
        }
    }
    return metrics;
}

void printMetrics(const std::string& label, const Metrics& m) {
    const double rms = std::sqrt(m.sumErrorSquared / static_cast<double>(m.samples));
    std::cout << std::left << std::setw(20) << label
              << " RMS(e)=" << std::setw(10) << rms
              << " max|e|=" << std::setw(10) << m.maxAbsError
              << " max|theta|=" << m.maxAbsParameter << '\n';
}

int main() {
    try {
        Config nominal;
        nominal.tau = 0.0;
        Config unmodeled = nominal;
        unmodeled.tau = 0.08;

        const Metrics m1 = simulate(nominal, "Chapter22_Lesson1_nominal.csv");
        const Metrics m2 = simulate(unmodeled, "Chapter22_Lesson1_unmodeled.csv");
        printMetrics("Nominal model", m1);
        printMetrics("Unmodeled pole", m2);
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

The C++17 program has no external dependency. It writes one CSV file for the nominal case and one for the unmodeled-pole case.

Chapter22_Lesson1.java

// Chapter22_Lesson1.java
// Sensitivity of MRAC to an unmodeled actuator pole.

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

public final class Chapter22_Lesson1 {
    private static final class Config {
        double a = 1.0;
        double b = 1.0;
        double am = 2.0;
        double bm = 2.0;
        double gamma = 30.0;
        double tau = 0.0;
        double dt = 1.0e-3;
        double tf = 30.0;
        double uLimit = 50.0;
    }

    private static final class Metrics {
        double sumErrorSquared;
        double maxAbsError;
        double maxAbsParameter;
        long samples;
    }

    private static double reference(double t) {
        return 1.0 + 0.35 * Math.sin(0.6 * t);
    }

    private static double disturbance(double t) {
        return t < 12.0 ? 0.0 : 0.05 * Math.sin(18.0 * t);
    }

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

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

    private static double[] rhs(double t, double[] y, Config cfg) {
        double x = y[0];
        double z = y[1];
        double xm = y[2];
        double thetaX = y[3];
        double thetaR = y[4];

        double r = reference(t);
        double e = x - xm;
        double u = clamp(thetaX * x + thetaR * r, cfg.uLimit);
        double plantInput = cfg.tau <= 0.0 ? u : z;
        double zDot = cfg.tau <= 0.0 ? 0.0 : (-z + u) / cfg.tau;

        return new double[] {
            -cfg.a * x + cfg.b * plantInput + disturbance(t),
            zDot,
            -cfg.am * xm + cfg.bm * r,
            -cfg.gamma * x * e,
            -cfg.gamma * r * e
        };
    }

    private static double[] rk4(double t, double[] y, Config cfg) {
        double h = cfg.dt;
        double[] k1 = rhs(t, y, cfg);
        double[] k2 = rhs(t + 0.5 * h, addScaled(y, k1, 0.5 * h), cfg);
        double[] k3 = rhs(t + 0.5 * h, addScaled(y, k2, 0.5 * h), cfg);
        double[] k4 = rhs(t + h, addScaled(y, k3, h), cfg);

        double[] next = new double[y.length];
        for (int i = 0; i < y.length; i++) {
            next[i] = y[i] + h * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
            if (!Double.isFinite(next[i])) {
                throw new IllegalStateException("Non-finite state encountered");
            }
        }
        return next;
    }

    private static Metrics simulate(Config cfg, String filename) throws IOException {
        Metrics metrics = new Metrics();
        double[] y = new double[5];
        long steps = Math.round(cfg.tf / cfg.dt);

        try (BufferedWriter writer = Files.newBufferedWriter(Path.of(filename))) {
            writer.write("t,x,z,xm,e,theta_x,theta_r,u\n");
            for (long k = 0; k <= steps; k++) {
                double t = k * cfg.dt;
                double r = reference(t);
                double e = y[0] - y[2];
                double u = clamp(y[3] * y[0] + y[4] * r, cfg.uLimit);

                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    t, y[0], y[1], y[2], e, y[3], y[4], u
                ));

                metrics.sumErrorSquared += e * e;
                metrics.maxAbsError = Math.max(metrics.maxAbsError, Math.abs(e));
                metrics.maxAbsParameter = Math.max(
                    metrics.maxAbsParameter,
                    Math.max(Math.abs(y[3]), Math.abs(y[4]))
                );
                metrics.samples++;

                if (k < steps) {
                    y = rk4(t, y, cfg);
                }
            }
        }
        return metrics;
    }

    private static void printMetrics(String label, Metrics m) {
        double rms = Math.sqrt(m.sumErrorSquared / m.samples);
        System.out.printf(
            Locale.US,
            "%-20s RMS(e)=%8.4f, max|e|=%8.4f, max|theta|=%8.4f%n",
            label, rms, m.maxAbsError, m.maxAbsParameter
        );
    }

    public static void main(String[] args) {
        try {
            Config nominal = new Config();
            nominal.tau = 0.0;
            Config unmodeled = new Config();
            unmodeled.tau = 0.08;

            Metrics m1 = simulate(nominal, "Chapter22_Lesson1_nominal_java.csv");
            Metrics m2 = simulate(unmodeled, "Chapter22_Lesson1_unmodeled_java.csv");
            printMetrics("Nominal model", m1);
            printMetrics("Unmodeled pole", m2);
        } catch (IOException | RuntimeException ex) {
            System.err.println("Error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

The Java implementation uses only the standard library and writes reproducible CSV output for external plotting.

Chapter22_Lesson1.m

% Chapter22_Lesson1.m
% Sensitivity of MRAC to an unmodeled actuator pole.

clear; clc; close all;

cfg.a = 1.0;
cfg.b = 1.0;
cfg.am = 2.0;
cfg.bm = 2.0;
cfg.gamma = 30.0;
cfg.dt = 1.0e-3;
cfg.tf = 30.0;
cfg.uLimit = 50.0;

nominal = simulateCase(cfg, 0.0);
unmodeled = simulateCase(cfg, 0.08);

printMetrics('Nominal model', nominal);
printMetrics('Unmodeled pole', unmodeled);

T = table(nominal.t, nominal.x, nominal.xm, nominal.e, ...
    nominal.thetaX, nominal.thetaR, unmodeled.x, unmodeled.xm, ...
    unmodeled.e, unmodeled.thetaX, unmodeled.thetaR, ...
    'VariableNames', {'t','x_nominal','xm_nominal','e_nominal', ...
    'theta_x_nominal','theta_r_nominal','x_unmodeled','xm_unmodeled', ...
    'e_unmodeled','theta_x_unmodeled','theta_r_unmodeled'});
writetable(T, 'Chapter22_Lesson1_results_matlab.csv');

figure;
plot(nominal.t, nominal.e, 'LineWidth', 1.2); hold on;
plot(unmodeled.t, unmodeled.e, 'LineWidth', 1.2);
grid on;
xlabel('time [s]'); ylabel('tracking error e');
legend('nominal model', 'unmodeled actuator pole', 'Location', 'best');
title('Sensitivity of the adaptive loop to neglected dynamics');

figure;
plot(unmodeled.t, unmodeled.thetaX, 'LineWidth', 1.2); hold on;
plot(unmodeled.t, unmodeled.thetaR, 'LineWidth', 1.2);
grid on;
xlabel('time [s]'); ylabel('adaptive parameters');
legend('\theta_x', '\theta_r', 'Location', 'best');
title('Parameter motion with an unmodeled actuator pole');

function data = simulateCase(cfg, tau)
    n = round(cfg.tf / cfg.dt) + 1;
    t = linspace(0, cfg.tf, n)';
    y = zeros(n, 5); % x, z, xm, theta_x, theta_r
    u = zeros(n, 1);

    for k = 1:n-1
        r = referenceSignal(t(k));
        u(k) = saturate(y(k,4)*y(k,1) + y(k,5)*r, cfg.uLimit);
        y(k+1,:) = rk4Step(t(k), y(k,:)', cfg.dt, cfg, tau)';
        if any(~isfinite(y(k+1,:)))
            error('Non-finite state at t = %.6f s', t(k+1));
        end
    end
    r = referenceSignal(t(end));
    u(end) = saturate(y(end,4)*y(end,1) + y(end,5)*r, cfg.uLimit);

    data.t = t;
    data.x = y(:,1);
    data.z = y(:,2);
    data.xm = y(:,3);
    data.thetaX = y(:,4);
    data.thetaR = y(:,5);
    data.u = u;
    data.e = data.x - data.xm;
end

function next = rk4Step(t, y, h, cfg, tau)
    k1 = dynamics(t, y, cfg, tau);
    k2 = dynamics(t + h/2, y + h*k1/2, cfg, tau);
    k3 = dynamics(t + h/2, y + h*k2/2, cfg, tau);
    k4 = dynamics(t + h, y + h*k3, cfg, tau);
    next = y + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end

function dy = dynamics(t, y, cfg, tau)
    x = y(1); z = y(2); xm = y(3);
    thetaX = y(4); thetaR = y(5);
    r = referenceSignal(t);
    e = x - xm;
    u = saturate(thetaX*x + thetaR*r, cfg.uLimit);

    if tau <= 0
        plantInput = u;
        zDot = 0;
    else
        plantInput = z;
        zDot = (-z + u)/tau;
    end

    dy = [
        -cfg.a*x + cfg.b*plantInput + disturbanceSignal(t);
        zDot;
        -cfg.am*xm + cfg.bm*r;
        -cfg.gamma*x*e;
        -cfg.gamma*r*e
    ];
end

function r = referenceSignal(t)
    r = 1.0 + 0.35*sin(0.6*t);
end

function d = disturbanceSignal(t)
    if t < 12.0
        d = 0.0;
    else
        d = 0.05*sin(18.0*t);
    end
end

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

function printMetrics(label, data)
    rmsError = sqrt(mean(data.e.^2));
    maxError = max(abs(data.e));
    maxParameter = max([abs(data.thetaX); abs(data.thetaR)]);
    fprintf('%-20s RMS(e)=%8.4f, max|e|=%8.4f, max|theta|=%8.4f\n', ...
        label, rmsError, maxError, maxParameter);
end

The same equations can be mapped directly into Simulink using five Integrator blocks for x, z, x_m, theta_x, and theta_r; Product blocks for theta_x*x, theta_r*r, x*e, and r*e; a Saturation block for u; and a Switch block that selects either the ideal input path or actuator state z. The fixed solver step should be materially smaller than the neglected time constant.

Chapter22_Lesson1.nb

Notebook[{
 Cell["Chapter 22, Lesson 1: Sensitivity of Adaptive Systems to Unmodeled Dynamics", "Title"],
 Cell["This notebook simulates the same gradient adaptive law for a nominal first-order plant and for a plant with a neglected actuator pole.", "Text"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"a", "=", "1.0"}], ";",
   RowBox[{"b", "=", "1.0"}], ";",
   RowBox[{"am", "=", "2.0"}], ";",
   RowBox[{"bm", "=", "2.0"}], ";",
   RowBox[{"gamma", "=", "30.0"}], ";",
   RowBox[{"uLimit", "=", "50.0"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"r", "[", "t_", "]", ":=", RowBox[{"1.0", "+", RowBox[{"0.35", " ", RowBox[{"Sin", "[", RowBox[{"0.6", " t"}], "]"}]}]}]}], ";",
   RowBox[{"d", "[", "t_", "]", ":=", RowBox[{"If", "[", RowBox[{RowBox[{"t", "<", "12.0"}], ",", "0.0", ",", RowBox[{"0.05", " ", RowBox[{"Sin", "[", RowBox[{"18.0", " t"}], "]"}]}]}], "]"}]}], ";",
   RowBox[{"sat", "[", "v_", "]", ":=", RowBox[{"Clip", "[", RowBox[{"v", ",", RowBox[{"{", RowBox[{RowBox[{"-", "uLimit"}], ",", "uLimit"}], "}"}]}], "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"simulate", "[", "tau_", "]", ":=", RowBox[{"Module", "[", RowBox[{
    RowBox[{"{", RowBox[{"eqns", ",", "ics", ",", "vars", ",", "sol", ",", "plantInput"}], "}"}], ",",
    RowBox[{
     RowBox[{"plantInput", "=", RowBox[{"If", "[", RowBox[{RowBox[{"tau", "==", "0"}], ",", RowBox[{"sat", "[", RowBox[{RowBox[{"thetaX", "[", "t", "]"}], " ", RowBox[{"x", "[", "t", "]"}], "+", RowBox[{RowBox[{"thetaR", "[", "t", "]"}], " ", RowBox[{"r", "[", "t", "]"}]}], "]"}], ",", RowBox[{"z", "[", "t", "]"}]}], "]"}]}], ";",
     RowBox[{"eqns", "=", RowBox[{"{",
      RowBox[{
       RowBox[{RowBox[{"x", "'", "[", "t", "]"}], "==", RowBox[{RowBox[{"-", "a"}], " ", RowBox[{"x", "[", "t", "]"}], "+", RowBox[{"b", " ", "plantInput"}], "+", RowBox[{"d", "[", "t", "]"}]}]}], ",",
       RowBox[{RowBox[{"xm", "'", "[", "t", "]"}], "==", RowBox[{RowBox[{"-", "am"}], " ", RowBox[{"xm", "[", "t", "]"}], "+", RowBox[{"bm", " ", RowBox[{"r", "[", "t", "]"}]}]}]}], ",",
       RowBox[{RowBox[{"thetaX", "'", "[", "t", "]"}], "==", RowBox[{RowBox[{"-", "gamma"}], " ", RowBox[{"x", "[", "t", "]"}], " ", RowBox[{"(", RowBox[{RowBox[{"x", "[", "t", "]"}], "-", RowBox[{"xm", "[", "t", "]"}]}], ")"}]}]}], ",",
       RowBox[{RowBox[{"thetaR", "'", "[", "t", "]"}], "==", RowBox[{RowBox[{"-", "gamma"}], " ", RowBox[{"r", "[", "t", "]"}], " ", RowBox[{"(", RowBox[{RowBox[{"x", "[", "t", "]"}], "-", RowBox[{"xm", "[", "t", "]"}]}], ")"}]}]}], ",",
       RowBox[{"If", "[", RowBox[{RowBox[{"tau", "==", "0"}], ",", RowBox[{RowBox[{"z", "'", "[", "t", "]"}], "==", "0"}], ",", RowBox[{RowBox[{"z", "'", "[", "t", "]"}], "==", RowBox[{RowBox[{"(", RowBox[{RowBox[{"-", RowBox[{"z", "[", "t", "]"}]}], "+", RowBox[{"sat", "[", RowBox[{RowBox[{RowBox[{"thetaX", "[", "t", "]"}], " ", RowBox[{"x", "[", "t", "]"}]}], "+", RowBox[{RowBox[{"thetaR", "[", "t", "]"}], " ", RowBox[{"r", "[", "t", "]"}]}]}], "]"}]}], ")"}], "/", "tau"}]}]}], "]"}]
      }], "}"}]}], ";",
     RowBox[{"ics", "=", RowBox[{"{", RowBox[{RowBox[{RowBox[{"x", "[", "0", "]"}], "==", "0"}], ",", RowBox[{RowBox[{"z", "[", "0", "]"}], "==", "0"}], ",", RowBox[{RowBox[{"xm", "[", "0", "]"}], "==", "0"}], ",", RowBox[{RowBox[{"thetaX", "[", "0", "]"}], "==", "0"}], ",", RowBox[{RowBox[{"thetaR", "[", "0", "]"}], "==", "0"}]}], "}"}]}], ";",
     RowBox[{"vars", "=", RowBox[{"{", RowBox[{"x", ",", "z", ",", "xm", ",", "thetaX", ",", "thetaR"}], "}"}]}], ";",
     RowBox[{"sol", "=", RowBox[{"NDSolveValue", "[", RowBox[{RowBox[{"Join", "[", RowBox[{"eqns", ",", "ics"}], "]"}], ",", "vars", ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}], ",", RowBox[{"Method", "->", RowBox[{"{", RowBox[{"\"TimeIntegration\"", ",", RowBox[{"\"ExplicitRungeKutta\"", "->", RowBox[{"{", RowBox[{"\"DifferenceOrder\"", "->", "4"}], "}"}]}]}], "}"}]}], ",", RowBox[{"MaxStepFraction", "->", RowBox[{"1", "/", "30000"}]}]}], "]"}]}], ";",
     "sol"
    }]}], "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"nominal", "=", RowBox[{"simulate", "[", "0", "]"}]}], ";",
   RowBox[{"unmodeled", "=", RowBox[{"simulate", "[", "0.08", "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[", RowBox[{
   RowBox[{"{", RowBox[{RowBox[{RowBox[{"nominal", "[", RowBox[{"[", "1", "]"}], "]"}], "[", "t", "]"}], "-", RowBox[{RowBox[{"nominal", "[", RowBox[{"[", "3", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{RowBox[{"unmodeled", "[", RowBox[{"[", "1", "]"}], "]"}], "[", "t", "]"}], "-", RowBox[{RowBox[{"unmodeled", "[", RowBox[{"[", "3", "]"}], "]"}], "[", "t", "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}], ",", RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"nominal\"", ",", "\"unmodeled pole\""}], "}"}]}], ",", RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t\"", ",", "\"e(t)\""}], "}"}]}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[", RowBox[{
   RowBox[{"{", RowBox[{RowBox[{RowBox[{"unmodeled", "[", RowBox[{"[", "4", "]"}], "]"}], "[", "t", "]"}], ",", RowBox[{RowBox[{"unmodeled", "[", RowBox[{"[", "5", "]"}], "]"}], "[", "t", "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}], ",", RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"theta_x\"", ",", "\"theta_r\""}], "}"}]}], ",", RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"t\"", ",", "\"parameters\""}], "}"}]}]}], "]"}]], "Input"]
}, WindowTitle -> "Chapter22_Lesson1", StyleDefinitions -> "Default.nb"]

The notebook uses NDSolveValue and plots the nominal and perturbed tracking errors together with the adaptive parameters.

With the nominal plant, the adaptive parameters move toward values that generate the desired first-order model, and the tracking error remains comparatively small. After the neglected pole is introduced, the command is no longer delivered instantaneously. The resulting mismatch enters the same channel as a disturbance, and the gradient law reacts by moving the parameters more aggressively.

The experiment should not be interpreted as a universal instability threshold. Its purpose is to demonstrate three reproducible trends:

1. the maximum parameter magnitude grows as the neglected time constant grows;
2. the RMS tracking error can worsen when adaptation is made faster;
3. small high-frequency disturbances become more consequential after phase lag is added.

Students should repeat the run for \(\tau\in\{0,0.02,0.05,0.08,0.12\}\) and \(\gamma\in\{2,10,30,60\}\), then construct a table of RMS error, maximum parameter norm, and maximum control effort.

Problem 1 — Nominal matching: For \(\dot{x}=-a x+b u\), reference model \(\dot{x}_m=-a_mx_m+b_mr\), and control \(u=\theta_xx+\theta_rr\), derive the ideal parameters.

Solution: Substituting the controller gives

\[ \dot{x}=(-a+b\theta_x)x+b\theta_r r. \]

Matching coefficients with the reference model requires

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

Therefore

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

Problem 2 — Why bounded disturbance does not automatically bound the parameters: Given \(\dot{V}=-a_me^2+e\delta\) with bounded \(\delta\), explain why the standard Lyapunov function does not prove bounded parameter error.

Solution: Young's inequality yields

\[ \dot{V}\leq -\frac{a_m}{2}e^2+\frac{\delta^2}{2a_m}. \]

The positive term prevents global negative semidefiniteness, and no term of the form \(-c\|\widetilde{\theta}\|^2\) appears. Since the update law is an integrator, a small biased product \(\phi e\) can accumulate indefinitely. A bound on the instantaneous tracking error is not a bound on the integral that defines the parameter estimate.

Problem 3 — Actuator mismatch bound: Starting from \(\tau\dot{z}=-z+u\), define \(\eta=z-u\) and derive an upper bound in terms of \(\dot{u}\).

Solution: Since \(z=\eta+u\),

\[ \tau(\dot{\eta}+\dot{u})=-\eta, \qquad \dot{\eta}=-\frac{1}{\tau}\eta-\dot{u}. \]

Variation of constants gives

\[ \eta(t)=e^{-t/\tau}\eta(0) -\int_0^t e^{-(t-\sigma)/\tau}\dot{u}(\sigma)\,d\sigma. \]

Taking absolute values and evaluating the exponential integral gives

\[ |\eta(t)|\leq e^{-t/\tau}|\eta(0)| +\tau\left(1-e^{-t/\tau}\right) \sup_{0\leq\sigma\leq t}|\dot{u}(\sigma)|. \]

Hence the steady mismatch is of order \(\tau\|\dot{u}\|_\infty\), not merely of order \(\tau\) independently of the controller.

Problem 4 — Frequency-domain interpretation: For \(\Delta_m(s)=-\tau s/(\tau s+1)\), find its low- and high-frequency limits.

Solution: At low frequency,

\[ |\Delta_m(j\omega)| =\frac{\tau\omega}{\sqrt{1+(\tau\omega)^2}} \approx \tau\omega, \qquad \tau\omega\ll 1. \]

At high frequency,

\[ \lim_{\omega → \infty}|\Delta_m(j\omega)|=1, \qquad \lim_{\omega → \infty}\angle\Delta_m(j\omega)=-180^\circ. \]

Thus the nominal approximation is accurate only below the neglected-pole frequency. An adaptive loop that pushes control activity upward in frequency can encounter nearly full relative model error and severe additional phase lag.

Problem 5 — Adaptation-gain scaling: Assume \(\|\phi\|\leq \bar{\phi}\) and \(|e|\leq \bar{e}\). Bound the contribution of parameter motion to \(|\dot{u}|\) for \(u=\theta^T\phi\) and \(\dot{\theta}=-\gamma b\phi e\).

Solution:

\[ |\dot{\theta}^{T}\phi| \leq \|\dot{\theta}\|\,\|\phi\| =\gamma b|e|\,\|\phi\|^2 \leq \gamma b\bar{e}\bar{\phi}^{2}. \]

Therefore the adaptive contribution to command rate scales linearly with \(\gamma\). Through the actuator mismatch bound, the neglected dynamics can contribute a perturbation of approximate scale \(\tau\gamma b\bar{e}\bar{\phi}^{2}\). This explains why a sufficiently large adaptation gain can invalidate a time-scale separation that appeared acceptable for a fixed-gain controller.

Problem 6 — Numerical investigation: Run the supplied program for all combinations of \(\tau\in\{0,0.02,0.05,0.08,0.12\}\) and \(\gamma\in\{2,10,30,60\}\). Identify a region in which increasing \(\gamma\) decreases nominal error but increases perturbed error or parameter magnitude.

Solution:

There is no single analytical table because the result depends on integration step, saturation, reference, and disturbance choices. A correct solution must report RMS error, maximum absolute error, maximum parameter magnitude, and maximum command magnitude for every pair. The expected qualitative result is a non-monotonic trade-off: adaptation initially improves low-frequency tracking, but beyond a problem-dependent range it excites the neglected pole and increases parameter motion. Numerical convergence must be checked by repeating selected cases with half the time step.

Nominal adaptive-control proofs rely on a model structure that permits exact cancellation between tracking-error and parameter-error cross terms. Unmodeled dynamics introduce an additional signal into the error equation, leaving an undamped parameter integrator exposed to persistent residual errors. A neglected stable actuator pole also links command rate to modeling error, so adaptation gain, regressor energy, and neglected bandwidth cannot be selected independently. Robust adaptive design therefore requires explicit mechanisms that bound parameters, suppress adaptation in noise-dominated regions, or limit the high-frequency authority of the adaptive channel.