After completing this lesson, students should be able to:

• derive how a scalar adaptation gain changes the poles of a local adaptive error model;
• derive transfer functions from measurement noise to tracking error and parameter estimates;
• compute stationary covariance caused by white measurement noise;
• obtain a sampling-dependent upper bound for an explicitly integrated adaptive law;
• select a matrix adaptation gain that accounts for regressor scaling and directionality;
• tune gains experimentally using tracking, parameter-activity, and control-activity metrics.

The preceding lesson established that noise and unmodeled dynamics enter an adaptive loop through both the feedback controller and the update law. Chapters 8 and 9 introduced leakage, dead zones, projection, and normalization. Here those devices are retained, but the focus is the remaining free design quantity: the adaptation gain \( \Gamma \) or its scalar counterpart \( \gamma \).

Consider a stable first-order tracking-error model with a linearly parameterized uncertainty:

\[ \dot e = -a e + c^{\mathsf T}\widetilde{\theta} + d, \qquad a > 0, \qquad \widetilde{\theta}=\widehat{\theta}-\theta^{\ast}. \]

In a complete MRAC derivation, \(c\) is formed from the plant high-frequency gain and the regressor. For local gain analysis, freeze the slowly varying regressor over a short interval and use \(c\) as a constant vector. A normalized gradient law with measured error \(e_m=e+n\) is

\[ \dot{\widehat{\theta}} =-\frac{\Gamma c}{m^2}(e+n) -\sigma\Gamma(\widehat{\theta}-\theta_0), \qquad m^2=1+\|c\|^2. \]

The matrix \(\Gamma=\Gamma^{\mathsf T}\succ0\) sets the learning rate in each parameter direction. Normalization limits the update caused by a large regressor, leakage limits drift, and projection may impose hard parameter bounds. None of these eliminates the need to tune \(\Gamma\).

For the scalar derivations below, absorb the normalizer and known sign into a scalar coupling coefficient \(c\neq0\) and write

\[ \dot e=-a e+c\widetilde\theta+d, \qquad \dot{\widetilde\theta}=-\gamma c(e+n). \]

With \(d=n=0\), consider

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

Differentiation gives

\[ \begin{aligned} \dot V &=e(-ae+c\widetilde\theta) +\frac{1}{\gamma}\widetilde\theta(-\gamma c e) \\ &=-a e^2. \end{aligned} \]

Therefore the ideal reduced-order system is stable for every \(\gamma>0\). This result is essential but incomplete for implementation. The expression \(\dot V=-ae^2\) does not reveal damping, oscillatory transient behavior, noise variance, actuator activity, sampling restrictions, or sensitivity to neglected dynamics. Those properties require additional analysis.

The homogeneous frozen-regressor dynamics are

\[ \frac{d}{dt} \begin{bmatrix}e\\ \widetilde\theta\end{bmatrix} = \begin{bmatrix}-a & c\\ -\gamma c & 0\end{bmatrix} \begin{bmatrix}e\\ \widetilde\theta\end{bmatrix}. \]

The characteristic polynomial is

\[ p(s)=s^2+a s+\gamma c^2. \]

Comparison with \(s^2+2\zeta\omega_n s+\omega_n^2\) gives

\[ \omega_n=|c|\sqrt{\gamma}, \qquad \zeta=\frac{a}{2|c|\sqrt{\gamma}}. \]

Thus the learning bandwidth increases as \(\sqrt{\gamma}\), but the damping ratio decreases as \(1/\sqrt{\gamma}\). The critical-damping gain is

\[ \gamma_{\mathrm{crit}}=\frac{a^2}{4c^2}. \]

A gain below this value produces two real modes; a gain above it produces a complex pair. A high gain can therefore make the parameter-error loop fast while simultaneously creating resonant tracking and control transients.

flowchart TD
  G["Increase adaptation gain gamma"] --> B["Learning bandwidth rises \napproximately as sqrt(gamma)"]
  B --> F["Faster initial tracking \nand parameter motion"]
  G --> D["Damping ratio falls \napproximately as \n1/sqrt(gamma)"]
  D --> R["More overshoot \nand resonance risk"]
  G --> N["Noise enters update \nwith larger multiplier"]
  N --> P["Parameter and control \nactivity increase"]
  P --> U["Greater exposure to sampling \nand unmodeled dynamics"]
        

Let the measurement noise enter only the update law. Taking Laplace transforms with zero initial conditions gives

\[ sE=-aE+c\widetilde\Theta, \qquad s\widetilde\Theta=-\gamma c(E+N). \]

Eliminating either state produces

\[ \frac{E(s)}{N(s)} =-\frac{\gamma c^2}{s^2+a s+\gamma c^2}, \]

\[ \frac{\widetilde\Theta(s)}{N(s)} =-\frac{\gamma c(s+a)}{s^2+a s+\gamma c^2}. \]

For sinusoidal noise at angular frequency \(\omega\), the parameter-noise gain is

\[ \left|H_{\theta n}(j\omega)\right|^2 =\frac{\gamma^2c^2(\omega^2+a^2)} {(\gamma c^2-\omega^2)^2+a^2\omega^2}. \]

Two mechanisms are visible. First, away from resonance and at high frequency, the magnitude is approximately \(\gamma|c|/\omega\); therefore a larger gain directly increases parameter response to high-frequency noise. Second, the denominator becomes small near \(\omega\approx|c|\sqrt{\gamma}\) when damping is low. Increasing the gain moves this resonant region upward in frequency and may place it near sensor noise, structural modes, actuator dynamics, or the Nyquist region.

Model the measurement noise as ideal white noise of intensity \(q\). The local stochastic system is

\[ d x=A x\,dt+B\sqrt{q}\,dW, \quad x=\begin{bmatrix}e\\ \widetilde\theta\end{bmatrix}, \quad A=\begin{bmatrix}-a&c\\-\gamma c&0\end{bmatrix}, \quad B=\begin{bmatrix}0\\-\gamma c\end{bmatrix}. \]

The stationary covariance \(P=\mathbb E[xx^{\mathsf T}]\) satisfies the continuous Lyapunov equation

\[ AP+PA^{\mathsf T}+qBB^{\mathsf T}=0. \]

Let

\[ P=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}. \]

Equating matrix entries gives the exact solution

\[ p_{11}=\frac{c^2\gamma q}{2a}, \qquad p_{12}=\frac{c\gamma q}{2}, \qquad p_{22}=\frac{q}{2a} \left(a^2\gamma+c^2\gamma^2\right). \]

Hence

\[ \operatorname{Var}(e)=\frac{c^2\gamma q}{2a}, \qquad \operatorname{Var}(\widetilde\theta) =\frac{q}{2a}\left(a^2\gamma+c^2\gamma^2\right). \]

This is a precise speed-versus-noise statement. The tracking-error variance grows linearly with the adaptation gain. The parameter variance contains both linear and quadratic terms and is asymptotically proportional to \(\gamma^2\). An arbitrarily large gain is therefore unacceptable even though the noise-free Lyapunov proof remains valid.

Implement the local dynamics with forward Euler and sample period \(T_s\):

\[ x_{k+1}=(I+T_sA)x_k. \]

The discrete characteristic polynomial is

\[ z^2-(2-aT_s)z+ \left(1-aT_s+\gamma c^2T_s^2\right)=0. \]

Applying the second-order Jury conditions, and assuming the ordinary implementation regime \(0<aT_s<2\), gives the key restriction

\[ 0<\gamma<\frac{a}{T_sc^2}. \]

This bound is not a universal adaptive-control theorem; it is the exact frozen-model bound for this Euler implementation. It demonstrates the general principle that the admissible adaptation gain shrinks as the sample period grows. Normalization changes the effective coupling to \(c/m\), but does not remove the need for a sampling-dependent limit.

A conservative digital design normally keeps the adaptation-loop natural frequency well below Nyquist:

\[ |c|\sqrt{\gamma} \leq \eta\frac{\pi}{T_s}, \qquad 0<\eta<1, \]

where \(\eta\) is chosen substantially below one when sensor filtering, computation delay, zero-order hold, or neglected dynamics are significant.

With a vector regressor, a single scalar gain can make one parameter direction excessively fast and another excessively slow. Define the normalized regressor covariance

\[ R_{\phi}=\mathbb E\left[ \frac{\phi\phi^{\mathsf T}}{m^2}\right], \qquad m^2=1+\phi^{\mathsf T}\phi. \]

Under an averaged prediction-error model, the mean parameter error obeys

\[ \frac{d}{dt}\mathbb E[\widetilde\theta] \approx-\Gamma R_{\phi}\mathbb E[\widetilde\theta]. \]

The learning rates are the eigenvalues of \(\Gamma R_{\phi}\). A diagonal gain may be selected as

\[ \Gamma=\operatorname{diag}(\gamma_1,\ldots,\gamma_p), \qquad \gamma_i=\frac{\bar\gamma}{\widehat R_{\phi,ii}+\varepsilon}, \]

where \(\varepsilon>0\) prevents division by a small excitation estimate. A fuller preconditioner is

\[ \Gamma=\bar\gamma (\widehat R_{\phi}+\varepsilon I)^{-1}. \]

This approximately equalizes deterministic convergence directions, but it may magnify poorly excited, noisy directions. Therefore eigenvalues of \(\Gamma\) should be capped and projection should be used when physical parameter bounds are known.

Let \(S_n=\mathbb E[\phi\phi^{\mathsf T}n^2/m^4]\). A local parameter-covariance approximation has the form

\[ \Gamma R_{\phi}P_{\theta} +P_{\theta}R_{\phi}\Gamma \approx\Gamma S_n\Gamma. \]

The right-hand side is quadratic in the gain matrix. Preconditioning must therefore consider both excitation and noise covariance, not only regressor magnitude.

flowchart TD
  A["Specify sample time, sensor bandwidth, actuator limits, and parameter bounds"] --> B["Scale states and regressors; enable normalization and projection"]
  B --> C["Estimate regressor covariance \nand measurement-noise level"]
  C --> D["Choose a small initial \nscalar or diagonal gain"]
  D --> E["Run noiseless command \nand uncertainty tests"]
  E --> F["Increase gain until tracking \ntransient is acceptable"]
  F --> G["Inject measured or \nworst-case sensor noise"]
  G --> H["Evaluate tracking RMS, \nparameter-rate RMS, control-rate RMS, \nand saturation"]
  H --> I{"All limits satisfied?"}
  I -->|"no"| J["Reduce gain or add \nfiltering, leakage, dead zone, \nor directional caps"]
  J --> G
  I -->|"yes"| K["Validate across operating \npoints, delays, and \nunmodeled-dynamics cases"]
  K --> L["Freeze documented gain bounds \nand monitoring thresholds"]
        

9.1 Recommended sequence

1. Scale first. Put states, commands, and regressors near comparable numerical magnitudes.
2. Set safety mechanisms. Configure control saturation, parameter projection, and finite-value checks before increasing the gain.
3. Find a speed floor. In a low-noise nominal model, increase the gain until the slowest acceptable transient is reached.
4. Measure noise costs. Use recorded sensor noise whenever available; synthetic white noise alone is insufficient for narrowband interference and quantization.
5. Check digital margins. Verify the gain against the integration step, computation delay, and anti-alias filtering.
6. Stress neglected dynamics. Include actuator lags, flexible modes, transport delays, and high-frequency plant uncertainty.
7. Prefer the smallest passing gain. Once all performance requirements are met, additional gain generally buys less tracking improvement while increasing robustness cost.

Use several metrics rather than tracking error alone:

\[ J_e=\sqrt{\frac{1}{T_2-T_1} \int_{T_1}^{T_2}e^2(t)\,dt}, \]

\[ J_{\dot\theta}=\sqrt{\frac{1}{T} \int_0^T\|\dot{\widehat\theta}(t)\|^2\,dt}, \qquad J_{\mathrm{TV}}=\int_0^T\|d\widehat\theta(t)\|, \]

\[ J_{\dot u}=\sqrt{\frac{1}{T} \int_0^T\dot u^2(t)\,dt}. \]

The supplied implementations compare \(\gamma\in\{0.5,5,40\}\) for a normalized first-order direct MRAC with measurement noise, mild leakage, projection, and control saturation. A representative Python run gives:

The largest gain has the smallest tracking RMSE in this reduced model, but its parameter total variation is nearly ten times that of the smallest gain. In a real system, the final selection must also pass actuator, sampling, structural-mode, and uncertainty tests. The table is therefore an illustration of the trade-off, not a universal choice of \(\gamma=40\).

This implementation uses NumPy for simulation and Matplotlib for plots. It writes one CSV file per gain and reports all four tuning metrics.

Chapter24_Lesson2.py

"""Chapter 24, Lesson 2: adaptation-gain tuning experiment.

Compares slow, moderate, and fast normalized-gradient MRAC gains under
measurement noise.  Requires NumPy and Matplotlib.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Iterable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class SimulationResult:
    gamma: float
    time: np.ndarray
    output: np.ndarray
    model_output: np.ndarray
    control: np.ndarray
    parameters: np.ndarray
    tracking_rmse: float
    parameter_total_variation: float
    parameter_rate_rms: float
    control_rate_rms: float


def reference_signal(t: float) -> float:
    """A bounded, switching command that excites both adaptive parameters."""
    return 1.0 if int(t / 4.0) % 2 == 0 else -0.6


def simulate(
    gamma: float,
    *,
    dt: float = 1.0e-3,
    final_time: float = 20.0,
    seed: int = 7,
) -> SimulationResult:
    """Simulate a first-order plant with normalized-gradient direct MRAC."""
    if gamma <= 0.0:
        raise ValueError("gamma must be positive")
    if dt <= 0.0 or final_time <= dt:
        raise ValueError("invalid integration interval")

    rng = np.random.default_rng(seed)
    time = np.arange(0.0, final_time + dt, dt)

    # Plant: y_dot = -a_p y + b_p u.
    a_p, b_p = 0.7, 1.0
    # Reference model: y_m_dot = -a_m y_m + b_m r.
    a_m, b_m = 1.5, 1.5

    # Sigma modification and projection are deliberately mild: the experiment
    # is about gain/noise trade-offs, not aggressive robustification.
    sigma = 2.0e-3
    theta_center = np.zeros(2)
    theta_limit = 5.0
    control_limit = 8.0

    y = 0.0
    y_m = 0.0
    theta = np.zeros(2)  # [theta_r, theta_y]

    output = np.empty_like(time)
    model_output = np.empty_like(time)
    control = np.empty_like(time)
    parameters = np.empty((time.size, 2))

    for k, t in enumerate(time):
        r = reference_signal(t)
        noise = (
            0.020 * rng.standard_normal()
            + 0.015 * np.sin(2.0 * np.pi * 35.0 * t)
        )

        y_m += dt * (-a_m * y_m + b_m * r)
        y_measured = y + noise
        e_measured = y_measured - y_m

        # u = theta_r r - theta_y y_measured = theta^T phi.
        phi = np.array([r, -y_measured])
        u = float(np.clip(theta @ phi, -control_limit, control_limit))
        y += dt * (-a_p * y + b_p * u)

        normalizer = 1.0 + float(phi @ phi)
        theta_dot = (
            -gamma * phi * e_measured / normalizer
            - sigma * gamma * (theta - theta_center)
        )
        theta += dt * theta_dot
        theta = np.clip(theta, -theta_limit, theta_limit)

        output[k] = y
        model_output[k] = y_m
        control[k] = u
        parameters[k] = theta

    evaluation = time >= 5.0
    tracking_rmse = float(
        np.sqrt(np.mean((output[evaluation] - model_output[evaluation]) ** 2))
    )
    parameter_increment = np.diff(parameters, axis=0)
    parameter_total_variation = float(
        np.sum(np.linalg.norm(parameter_increment, axis=1))
    )
    parameter_rate_rms = float(
        np.sqrt(np.mean(np.sum((parameter_increment / dt) ** 2, axis=1)))
    )
    control_rate_rms = float(np.sqrt(np.mean((np.diff(control) / dt) ** 2)))

    return SimulationResult(
        gamma=gamma,
        time=time,
        output=output,
        model_output=model_output,
        control=control,
        parameters=parameters,
        tracking_rmse=tracking_rmse,
        parameter_total_variation=parameter_total_variation,
        parameter_rate_rms=parameter_rate_rms,
        control_rate_rms=control_rate_rms,
    )


def save_csv(result: SimulationResult, directory: Path) -> Path:
    safe_gain = str(result.gamma).replace(".", "p")
    destination = directory / f"Chapter24_Lesson2_gamma_{safe_gain}.csv"
    data = np.column_stack(
        (
            result.time,
            result.output,
            result.model_output,
            result.control,
            result.parameters[:, 0],
            result.parameters[:, 1],
        )
    )
    np.savetxt(
        destination,
        data,
        delimiter=",",
        header="time,y,y_model,u,theta_r,theta_y",
        comments="",
    )
    return destination


def plot_results(results: Iterable[SimulationResult], directory: Path) -> None:
    results = list(results)

    plt.figure(figsize=(10, 6))
    plt.plot(results[0].time, results[0].model_output, "k--", label="reference model")
    for result in results:
        plt.plot(result.time, result.output, label=f"gamma={result.gamma:g}")
    plt.xlabel("time [s]")
    plt.ylabel("output")
    plt.title("Tracking speed versus adaptation gain")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(directory / "Chapter24_Lesson2_tracking.png", dpi=180)
    plt.close()

    plt.figure(figsize=(10, 6))
    for result in results:
        plt.plot(
            result.time,
            result.parameters[:, 0],
            label=f"theta_r, gamma={result.gamma:g}",
        )
    plt.xlabel("time [s]")
    plt.ylabel("adaptive parameter")
    plt.title("Parameter activity under measurement noise")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(directory / "Chapter24_Lesson2_parameters.png", dpi=180)
    plt.close()


def main() -> None:
    output_directory = Path(__file__).resolve().parent
    gains = (0.5, 5.0, 40.0)
    results = [simulate(gain) for gain in gains]

    print("gamma, tracking_RMSE, parameter_TV, parameter_rate_RMS, control_rate_RMS")
    for result in results:
        print(
            f"{result.gamma:6.2f}, "
            f"{result.tracking_rmse:13.6f}, "
            f"{result.parameter_total_variation:12.6f}, "
            f"{result.parameter_rate_rms:18.6f}, "
            f"{result.control_rate_rms:16.6f}"
        )
        save_csv(result, output_directory)

    plot_results(results, output_directory)
    print(f"Results written to: {output_directory}")


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library and exports CSV data for external plotting or comparison.

Chapter24_Lesson2.cpp

// Chapter 24, Lesson 2: adaptation-gain tuning experiment.
// Build: g++ -std=c++17 -O2 Chapter24_Lesson2.cpp -o Chapter24_Lesson2

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

struct Metrics {
    double tracking_rmse{};
    double parameter_total_variation{};
    double parameter_rate_rms{};
    double control_rate_rms{};
};

static double reference_signal(double t) {
    return (static_cast<int>(t / 4.0) % 2 == 0) ? 1.0 : -0.6;
}

static double clamp(double value, double lower, double upper) {
    return std::max(lower, std::min(value, upper));
}

static std::string gain_tag(double gamma) {
    std::ostringstream stream;
    stream << gamma;
    std::string tag = stream.str();
    std::replace(tag.begin(), tag.end(), '.', 'p');
    return tag;
}

static Metrics simulate(double gamma, double dt = 1.0e-3, double final_time = 20.0) {
    if (gamma <= 0.0 || dt <= 0.0 || final_time <= dt) {
        throw std::invalid_argument("invalid simulation parameter");
    }

    constexpr double pi = 3.14159265358979323846;
    constexpr double a_p = 0.7;
    constexpr double b_p = 1.0;
    constexpr double a_m = 1.5;
    constexpr double b_m = 1.5;
    constexpr double sigma = 2.0e-3;
    constexpr double theta_limit = 5.0;
    constexpr double control_limit = 8.0;

    std::mt19937 generator(7U);
    std::normal_distribution<double> normal(0.0, 1.0);

    double y = 0.0;
    double y_model = 0.0;
    std::array<double, 2> theta{0.0, 0.0};
    std::array<double, 2> previous_theta = theta;
    double previous_control = 0.0;

    double squared_error_sum = 0.0;
    std::size_t error_samples = 0U;
    double total_variation = 0.0;
    double parameter_rate_squared_sum = 0.0;
    double control_rate_squared_sum = 0.0;
    std::size_t rate_samples = 0U;

    const std::string filename =
        "Chapter24_Lesson2_gamma_" + gain_tag(gamma) + ".csv";
    std::ofstream csv(filename);
    if (!csv) {
        throw std::runtime_error("cannot open output CSV file");
    }
    csv << "time,y,y_model,u,theta_r,theta_y\n";
    csv << std::setprecision(12);

    const std::size_t steps = static_cast<std::size_t>(std::llround(final_time / dt));
    for (std::size_t k = 0; k <= steps; ++k) {
        const double t = static_cast<double>(k) * dt;
        const double r = reference_signal(t);
        const double noise =
            0.020 * normal(generator) + 0.015 * std::sin(2.0 * pi * 35.0 * t);

        y_model += dt * (-a_m * y_model + b_m * r);
        const double y_measured = y + noise;
        const double e_measured = y_measured - y_model;
        const std::array<double, 2> phi{r, -y_measured};

        const double raw_control = theta[0] * phi[0] + theta[1] * phi[1];
        const double control = clamp(raw_control, -control_limit, control_limit);
        y += dt * (-a_p * y + b_p * control);

        const double normalizer = 1.0 + phi[0] * phi[0] + phi[1] * phi[1];
        for (std::size_t i = 0; i < theta.size(); ++i) {
            const double theta_dot =
                -gamma * phi[i] * e_measured / normalizer - sigma * gamma * theta[i];
            theta[i] = clamp(theta[i] + dt * theta_dot, -theta_limit, theta_limit);
        }

        if (t >= 5.0) {
            const double error = y - y_model;
            squared_error_sum += error * error;
            ++error_samples;
        }

        if (k > 0U) {
            const double d0 = theta[0] - previous_theta[0];
            const double d1 = theta[1] - previous_theta[1];
            total_variation += std::hypot(d0, d1);
            parameter_rate_squared_sum += (d0 * d0 + d1 * d1) / (dt * dt);
            const double control_rate = (control - previous_control) / dt;
            control_rate_squared_sum += control_rate * control_rate;
            ++rate_samples;
        }

        previous_theta = theta;
        previous_control = control;
        csv << t << ',' << y << ',' << y_model << ',' << control << ','
            << theta[0] << ',' << theta[1] << '\n';
    }

    Metrics metrics;
    metrics.tracking_rmse = std::sqrt(squared_error_sum / static_cast<double>(error_samples));
    metrics.parameter_total_variation = total_variation;
    metrics.parameter_rate_rms =
        std::sqrt(parameter_rate_squared_sum / static_cast<double>(rate_samples));
    metrics.control_rate_rms =
        std::sqrt(control_rate_squared_sum / static_cast<double>(rate_samples));
    return metrics;
}

int main() {
    try {
        const std::vector<double> gains{0.5, 5.0, 40.0};
        std::cout << "gamma, tracking_RMSE, parameter_TV, parameter_rate_RMS, control_rate_RMS\n";
        for (double gamma : gains) {
            const Metrics metrics = simulate(gamma);
            std::cout << std::fixed << std::setprecision(6)
                      << gamma << ", " << metrics.tracking_rmse << ", "
                      << metrics.parameter_total_variation << ", "
                      << metrics.parameter_rate_rms << ", "
                      << metrics.control_rate_rms << '\n';
        }
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}

The Java version uses a record for summary metrics, deterministic random seeding, bounded updates, and standard CSV output.

Chapter24_Lesson2.java

// Chapter 24, Lesson 2: adaptation-gain tuning experiment.
// Build: javac Chapter24_Lesson2.java
// Run:   java Chapter24_Lesson2

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

public final class Chapter24_Lesson2 {
    private Chapter24_Lesson2() {
    }

    private record Metrics(
            double trackingRmse,
            double parameterTotalVariation,
            double parameterRateRms,
            double controlRateRms) {
    }

    private static double referenceSignal(double time) {
        return (((int) (time / 4.0)) % 2 == 0) ? 1.0 : -0.6;
    }

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

    private static String gainTag(double gamma) {
        return Double.toString(gamma).replace('.', 'p');
    }

    private static Metrics simulate(double gamma, double dt, double finalTime)
            throws IOException {
        if (gamma <= 0.0 || dt <= 0.0 || finalTime <= dt) {
            throw new IllegalArgumentException("invalid simulation parameter");
        }

        final double plantA = 0.7;
        final double plantB = 1.0;
        final double modelA = 1.5;
        final double modelB = 1.5;
        final double sigma = 2.0e-3;
        final double thetaLimit = 5.0;
        final double controlLimit = 8.0;

        final Random random = new Random(7L);
        double y = 0.0;
        double yModel = 0.0;
        final double[] theta = {0.0, 0.0};
        final double[] previousTheta = {0.0, 0.0};
        double previousControl = 0.0;

        double squaredErrorSum = 0.0;
        long errorSamples = 0L;
        double totalVariation = 0.0;
        double parameterRateSquaredSum = 0.0;
        double controlRateSquaredSum = 0.0;
        long rateSamples = 0L;

        final Path csvPath = Path.of(
                "Chapter24_Lesson2_gamma_" + gainTag(gamma) + ".csv");
        try (BufferedWriter csv = Files.newBufferedWriter(
                csvPath, StandardCharsets.UTF_8)) {
            csv.write("time,y,y_model,u,theta_r,theta_y\n");
            final int steps = (int) Math.round(finalTime / dt);
            for (int k = 0; k <= steps; ++k) {
                final double time = k * dt;
                final double reference = referenceSignal(time);
                final double noise = 0.020 * random.nextGaussian()
                        + 0.015 * Math.sin(2.0 * Math.PI * 35.0 * time);

                yModel += dt * (-modelA * yModel + modelB * reference);
                final double measuredY = y + noise;
                final double measuredError = measuredY - yModel;
                final double[] phi = {reference, -measuredY};

                final double rawControl = theta[0] * phi[0] + theta[1] * phi[1];
                final double control = clamp(rawControl, -controlLimit, controlLimit);
                y += dt * (-plantA * y + plantB * control);

                final double normalizer = 1.0 + phi[0] * phi[0] + phi[1] * phi[1];
                for (int i = 0; i < theta.length; ++i) {
                    final double thetaDot = -gamma * phi[i] * measuredError / normalizer
                            - sigma * gamma * theta[i];
                    theta[i] = clamp(
                            theta[i] + dt * thetaDot,
                            -thetaLimit,
                            thetaLimit);
                }

                if (time >= 5.0) {
                    final double error = y - yModel;
                    squaredErrorSum += error * error;
                    ++errorSamples;
                }

                if (k > 0) {
                    final double d0 = theta[0] - previousTheta[0];
                    final double d1 = theta[1] - previousTheta[1];
                    totalVariation += Math.hypot(d0, d1);
                    parameterRateSquaredSum += (d0 * d0 + d1 * d1) / (dt * dt);
                    final double controlRate = (control - previousControl) / dt;
                    controlRateSquaredSum += controlRate * controlRate;
                    ++rateSamples;
                }

                previousTheta[0] = theta[0];
                previousTheta[1] = theta[1];
                previousControl = control;
                csv.write(String.format(
                        Locale.US,
                        "%.9f,%.9f,%.9f,%.9f,%.9f,%.9f%n",
                        time,
                        y,
                        yModel,
                        control,
                        theta[0],
                        theta[1]));
            }
        }

        return new Metrics(
                Math.sqrt(squaredErrorSum / errorSamples),
                totalVariation,
                Math.sqrt(parameterRateSquaredSum / rateSamples),
                Math.sqrt(controlRateSquaredSum / rateSamples));
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        final double[] gains = {0.5, 5.0, 40.0};
        System.out.println(
                "gamma, tracking_RMSE, parameter_TV, parameter_rate_RMS, control_rate_RMS");
        try {
            for (double gamma : gains) {
                final Metrics metrics = simulate(gamma, 1.0e-3, 20.0);
                System.out.printf(
                        Locale.US,
                        "%.2f, %.6f, %.6f, %.6f, %.6f%n",
                        gamma,
                        metrics.trackingRmse(),
                        metrics.parameterTotalVariation(),
                        metrics.parameterRateRms(),
                        metrics.controlRateRms());
            }
        } catch (IOException | IllegalArgumentException error) {
            System.err.println("Error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB file performs the numerical experiment and documents a direct Simulink block mapping. For deployment-oriented studies, replace the forward-Euler loop with the intended fixed-step solver and include the real sensor and actuator blocks.

Chapter24_Lesson2.m

%% Chapter 24, Lesson 2: adaptation-gain tuning experiment
% Compares slow, moderate, and fast normalized-gradient MRAC gains under
% measurement noise. Uses only base MATLAB functions.

clear; clc; close all;
rng(7, 'twister');

gains = [0.5, 5.0, 40.0];
results = cell(size(gains));

fprintf('gamma, tracking_RMSE, parameter_TV, parameter_rate_RMS, control_rate_RMS\n');
for i = 1:numel(gains)
    results{i} = simulateGain(gains(i), 1.0e-3, 20.0, 7);
    result = results{i};
    fprintf('%.2f, %.6f, %.6f, %.6f, %.6f\n', ...
        gains(i), result.trackingRMSE, result.parameterTV, ...
        result.parameterRateRMS, result.controlRateRMS);

    tag = strrep(num2str(gains(i)), '.', 'p');
    outputTable = table(result.time, result.y, result.yModel, result.u, ...
        result.theta(:, 1), result.theta(:, 2), ...
        'VariableNames', {'time', 'y', 'y_model', 'u', 'theta_r', 'theta_y'});
    writetable(outputTable, ['Chapter24_Lesson2_gamma_', tag, '.csv']);
end

figure;
plot(results{1}.time, results{1}.yModel, 'k--', 'LineWidth', 1.5);
hold on;
for i = 1:numel(gains)
    plot(results{i}.time, results{i}.y, 'LineWidth', 1.0);
end
hold off; grid on;
xlabel('time [s]'); ylabel('output');
title('Tracking speed versus adaptation gain');
legend('reference model', 'gamma=0.5', 'gamma=5', 'gamma=40', ...
    'Location', 'best');
exportgraphics(gcf, 'Chapter24_Lesson2_tracking_matlab.png', 'Resolution', 180);

figure;
hold on;
for i = 1:numel(gains)
    plot(results{i}.time, results{i}.theta(:, 1), 'LineWidth', 1.0);
end
hold off; grid on;
xlabel('time [s]'); ylabel('theta_r');
title('Parameter activity under measurement noise');
legend('gamma=0.5', 'gamma=5', 'gamma=40', 'Location', 'best');
exportgraphics(gcf, 'Chapter24_Lesson2_parameters_matlab.png', 'Resolution', 180);

%% Optional Simulink mapping
% Implement the same experiment with these blocks:
% 1. Plant State-Space block: A=-0.7, B=1, C=1, D=0.
% 2. Reference Model State-Space block: A=-1.5, B=1.5, C=1, D=0.
% 3. MATLAB Function block implementing the normalized adaptive law.
% 4. Band-Limited White Noise plus a 35-Hz Sine Wave before the measured y.
% 5. Saturation blocks on u and theta. Use fixed-step ode4 with Ts=1e-3.

function result = simulateGain(gamma, dt, finalTime, seed)
    arguments
        gamma (1,1) double {mustBePositive}
        dt (1,1) double {mustBePositive}
        finalTime (1,1) double {mustBePositive}
        seed (1,1) double {mustBeInteger}
    end

    rng(seed, 'twister');
    time = (0:dt:finalTime).';
    sampleCount = numel(time);

    plantA = 0.7; plantB = 1.0;
    modelA = 1.5; modelB = 1.5;
    sigma = 2.0e-3;
    thetaLimit = 5.0;
    controlLimit = 8.0;

    y = zeros(sampleCount, 1);
    yModel = zeros(sampleCount, 1);
    u = zeros(sampleCount, 1);
    theta = zeros(sampleCount, 2);

    yState = 0.0;
    modelState = 0.0;
    thetaState = [0.0; 0.0];

    for k = 1:sampleCount
        t = time(k);
        if mod(floor(t / 4.0), 2) == 0
            reference = 1.0;
        else
            reference = -0.6;
        end

        noise = 0.020 * randn + 0.015 * sin(2 * pi * 35.0 * t);
        modelState = modelState + dt * (-modelA * modelState + modelB * reference);
        measuredY = yState + noise;
        measuredError = measuredY - modelState;
        phi = [reference; -measuredY];

        control = min(max(thetaState.' * phi, -controlLimit), controlLimit);
        yState = yState + dt * (-plantA * yState + plantB * control);

        normalizer = 1.0 + phi.' * phi;
        thetaDot = -gamma * phi * measuredError / normalizer ...
            - sigma * gamma * thetaState;
        thetaState = thetaState + dt * thetaDot;
        thetaState = min(max(thetaState, -thetaLimit), thetaLimit);

        y(k) = yState;
        yModel(k) = modelState;
        u(k) = control;
        theta(k, :) = thetaState.';
    end

    evaluation = time >= 5.0;
    result.time = time;
    result.y = y;
    result.yModel = yModel;
    result.u = u;
    result.theta = theta;
    result.trackingRMSE = sqrt(mean((y(evaluation) - yModel(evaluation)).^2));
    parameterIncrement = diff(theta, 1, 1);
    result.parameterTV = sum(vecnorm(parameterIncrement, 2, 2));
    result.parameterRateRMS = sqrt(mean(sum((parameterIncrement / dt).^2, 2)));
    result.controlRateRMS = sqrt(mean((diff(u) / dt).^2));
end

The notebook expression contains a complete input cell, summary dataset, plots, and CSV export. It uses built-in numerical and statistical functions.

Chapter24_Lesson2.nb

Notebook[{
 Cell["Chapter 24, Lesson 2: Tuning of Adaptation Gains", "Title"],
 Cell["Evaluate the input cell to compare normalized-gradient MRAC gains under measurement noise.", "Text"],
 Cell[BoxData[RowBox[{"ToExpression", "[",
"\"ClearAll[\\\"Global`*\\\"];
referenceSignal[t_] := If[EvenQ[Floor[t/4.0]], 1.0, -0.6];
clamp[x_, lo_, hi_] := Min[Max[x, lo], hi];
simulateGain[gamma_?Positive, dt_: 0.001, finalTime_: 20.0, seed_: 7] :=
 Module[
  {time, plantA = 0.7, plantB = 1.0, modelA = 1.5, modelB = 1.5,
   sigma = 0.002, thetaLimit = 5.0, controlLimit = 8.0, y = 0.0,
   yModel = 0.0, theta = {0.0, 0.0}, rows = {}, r, noise, measuredY,
   measuredError, phi, u, normalizer, thetaDot, data, eval, parameterIncrement,
   trackingRMSE, parameterTV, parameterRateRMS, controlRateRMS},
  SeedRandom[seed]; time = Range[0.0, finalTime, dt];
  Do[
   r = referenceSignal[t];
   noise = 0.020 RandomVariate[NormalDistribution[0, 1]] + 0.015 Sin[2 Pi 35.0 t];
   yModel = yModel + dt (-modelA yModel + modelB r);
   measuredY = y + noise; measuredError = measuredY - yModel;
   phi = {r, -measuredY}; u = clamp[theta.phi, -controlLimit, controlLimit];
   y = y + dt (-plantA y + plantB u); normalizer = 1.0 + phi.phi;
   thetaDot = -gamma phi measuredError/normalizer - sigma gamma theta;
   theta = Map[clamp[#, -thetaLimit, thetaLimit] &, theta + dt thetaDot];
   AppendTo[rows, {t, y, yModel, u, theta[[1]], theta[[2]]}],
   {t, time}];
  data = N[rows]; eval = Select[data, #[[1]] >= 5.0 &];
  trackingRMSE = Sqrt[Mean[(eval[[All, 2]] - eval[[All, 3]])^2]];
  parameterIncrement = Differences[data[[All, 5 ;; 6]]];
  parameterTV = Total[Norm /@ parameterIncrement];
  parameterRateRMS = Sqrt[Mean[Total[(parameterIncrement/dt)^2, {2}]]];
  controlRateRMS = Sqrt[Mean[(Differences[data[[All, 4]]]/dt)^2]];
  <|\\\"Gamma\\\" -> gamma, \\\"Data\\\" -> data,
    \\\"TrackingRMSE\\\" -> trackingRMSE, \\\"ParameterTV\\\" -> parameterTV,
    \\\"ParameterRateRMS\\\" -> parameterRateRMS,
    \\\"ControlRateRMS\\\" -> controlRateRMS|>
 ];
gains = {0.5, 5.0, 40.0}; results = simulateGain /@ gains;
summary = Dataset[Map[KeyDrop[#, {\\\"Data\\\"}] &, results]]; Print[summary];
trackingPlot = Show[
 ListLinePlot[results[[1, \\\"Data\\\"]][[All, {1, 3}]],
  PlotStyle -> {Black, Dashed}, PlotLegends -> {\\\"reference model\\\"}],
 ListLinePlot[Table[results[[i, \\\"Data\\\"]][[All, {1, 2}]], {i, Length[gains]}],
  PlotLegends -> (\\\"gamma=\\\" <> ToString[#] & /@ gains)],
 Frame -> True, FrameLabel -> {\\\"time [s]\\\", \\\"output\\\"},
 PlotLabel -> \\\"Tracking speed versus adaptation gain\\\"];
Print[trackingPlot];
parameterPlot = ListLinePlot[
 Table[results[[i, \\\"Data\\\"]][[All, {1, 5}]], {i, Length[gains]}],
 Frame -> True, FrameLabel -> {\\\"time [s]\\\", \\\"theta_r\\\"},
 PlotLegends -> (\\\"gamma=\\\" <> ToString[#] & /@ gains),
 PlotLabel -> \\\"Parameter activity under measurement noise\\\"];
Print[parameterPlot];
Do[
 Export[\\\"Chapter24_Lesson2_gamma_\\\" <>
   StringReplace[ToString[gains[[i]]], \\\".\\\" -> \\\"p\\\"] <> \\\".csv\\\",
  Prepend[results[[i, \\\"Data\\\"]],
   {\\\"time\\\", \\\"y\\\", \\\"y_model\\\", \\\"u\\\", \\\"theta_r\\\", \\\"theta_y\\\"}]],
 {i, Length[gains]}];\"",
 "]"}]], "Input"]}, WindowTitle -> "Chapter24_Lesson2"]        

A gain is too small when tracking remains unacceptably slow under verified excitation, parameter estimates move smoothly but fail to compensate uncertainty within the required time, and no actuator or noise limitation is active.

A gain is too large when one or more of the following appears:

• parameter estimates show high-frequency chatter or repeated projection-bound impacts;
• control-rate RMS increases strongly while tracking improvement becomes marginal;
• error or control spectra develop a peak near the adaptive-loop natural frequency;
• performance changes sharply when the sample period or solver is changed slightly;
• small sensor bias causes large parameter drift before leakage or projection intervenes;
• actuator saturation becomes frequent, invalidating the nominal error model;
• including a realistic actuator lag or flexible mode reverses the apparent benefit of the high gain.

Filtering may reduce measurement-noise energy, but it also adds phase lag. Therefore “increase the gain and then filter harder” is not a generally valid remedy. Gain, filter bandwidth, sample rate, leakage, dead-zone threshold, projection limits, and actuator constraints must be validated as one coupled implementation.

Problem 1: Pole Placement by Adaptation Gain

For \(a=4\) and \(c=2\), determine the adaptation gain that gives critical damping. Then determine the damping ratio for \(\gamma=4\).

Solution:

\[ \gamma_{\mathrm{crit}}=\frac{a^2}{4c^2} =\frac{16}{16}=1. \]

For \(\gamma=4\),

\[ \omega_n=|c|\sqrt\gamma=2(2)=4, \qquad \zeta=\frac{a}{2|c|\sqrt\gamma} =\frac{4}{8}=0.5. \]

The higher gain doubles the natural frequency relative to the critically damped case, but reduces the damping ratio to 0.5.

Problem 2: Parameter Variance Growth

Let \(a=2\), \(c=1\), and \(q=0.01\). Compute the stationary parameter variance for \(\gamma=1\) and \(\gamma=10\).

Solution:

\[ \operatorname{Var}(\widetilde\theta) =\frac{0.01}{4}(4\gamma+\gamma^2). \]

At \(\gamma=1\),

\[ \operatorname{Var}(\widetilde\theta)=0.0025(5)=0.0125. \]

At \(\gamma=10\),

\[ \operatorname{Var}(\widetilde\theta)=0.0025(140)=0.35. \]

Increasing the gain by a factor of ten increases this variance by a factor of twenty-eight because the quadratic term is already dominant.

Problem 3: Euler Sampling Bound

A frozen adaptive error model has \(a=5\), \(c=2\), and sample period \(T_s=0.01\) s. Find the Euler stability upper bound on \(\gamma\). Would \(\gamma=100\) be admissible under this local model?

Solution:

\[ \gamma<\frac{a}{T_sc^2} =\frac{5}{0.01(4)}=125. \]

The value 100 is inside the local Euler stability interval. It is not automatically acceptable because the bound excludes delays, noise, actuator dynamics, and unmodeled plant modes. A substantial margin below 125 is normally required.

Problem 4: Directional Gain Scaling

Suppose the estimated normalized regressor covariance is \(\widehat R_\phi=\operatorname{diag}(100,1)\). Compare the averaged learning rates for \(\Gamma=I\) and for \(\Gamma=5(\widehat R_\phi+0.1I)^{-1}\).

Solution:

With \(\Gamma=I\), the learning-rate eigenvalues are 100 and 1, so the first direction is one hundred times faster. With the preconditioned gain,

\[ \Gamma R_\phi =5\operatorname{diag} \left(\frac{100}{100.1},\frac{1}{1.1}\right) \approx\operatorname{diag}(4.995,4.545). \]

The deterministic rates are nearly equalized. Before using this gain, the noise level in the weakly excited second direction must be checked, because inverse-covariance scaling can amplify that direction.

Problem 5: Selecting a Gain from Experimental Metrics

Three candidates produce the following results: gain A has tracking RMSE 0.08 and parameter-rate RMS 0.20; gain B has RMSE 0.04 and parameter-rate RMS 0.55; gain C has RMSE 0.035 and parameter-rate RMS 2.4. The tracking requirement is RMSE no greater than 0.05, and the parameter-rate limit is 0.8. Which gain should be selected if all other tests pass?

Solution:

Gain A fails the tracking requirement. Gains B and C pass tracking, but C violates the parameter-rate limit. Gain B is the only feasible choice. It is also consistent with the “smallest passing gain” rule: moving from B to C gives little tracking improvement while greatly increasing adaptive activity.

The adaptation gain is a bandwidth-setting parameter, not merely a convergence-speed knob. For the frozen scalar error model, the natural frequency grows as \(\sqrt\gamma\) while the damping ratio decreases as \(1/\sqrt\gamma\). Measurement noise is multiplied directly by the gain in the update law, and exact covariance analysis shows linear growth of tracking-error variance and linear-plus-quadratic growth of parameter variance. Digital implementation adds a sample-period-dependent upper bound. Consequently, gains should be scaled by regressor statistics, increased only until tracking requirements are met, and then validated using parameter and control activity, sampling, delay, saturation, and unmodeled-dynamics tests.