After completing this lesson, students should be able to:

• model amplitude saturation as a static nonlinearity between the controller command and the plant input;
• derive the tracking-error dynamics of a first-order MRAC system in the presence of actuator saturation;
• identify the additional term that destroys the nominal Lyapunov cancellation;
• test whether a reference-model trajectory is physically feasible under a specified actuator limit;
• explain why persistent tracking error can produce adaptive-parameter drift even when the actual plant input remains bounded;
• use simulation diagnostics to separate transient saturation, feasible sustained saturation, and infeasible commands.

The present lesson diagnoses the problem. Detailed modifications of the update law, reference model, and anti-windup structure are reserved for Lessons 2 and 3 of this chapter.

Let \(u_c(t)\) denote the signal calculated by the adaptive controller and let \(u_a(t)\) denote the signal actually delivered by an amplitude-limited actuator. For a symmetric actuator limit \(u_{\max} > 0\),

\[ u_a = \operatorname{sat}_{u_{\max}}(u_c) = \begin{cases} -u_{\max}, & u_c < -u_{\max},\\ u_c, & |u_c| \leq u_{\max},\\ u_{\max}, & u_c > u_{\max}. \end{cases} \]

The difference between the realized and requested controls is the control-deficiency signal \(\Delta u\):

\[ \Delta u \triangleq u_a-u_c. \]

Thus \(\Delta u=0\) when the actuator is linear and \(\Delta u\neq 0\) during saturation. The plant always receives

\[ u_a=u_c+\Delta u, \]

so the deficiency behaves as an additional, controller-dependent input disturbance. Unlike an independent bounded disturbance, \(\Delta u\) depends on the adaptive parameters through \(u_c\).

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  R --> AC["Adaptive controller"]
  X["Plant state x"] --> AC
  AC --> UC["Command u_c"]
  UC --> SAT["Amplitude saturation"]
  SAT --> UA["Actual input u_a"]
  UA --> P["Plant"]
  P --> X
  RM --> E["Tracking error e = x - x_m"]
  X --> E
  E --> LAW["Parameter update law"]
  LAW --> AC
  UC --> DU["Control deficiency Delta_u = u_a - u_c"]
  UA --> DU
        

Other actuator restrictions include rate limits, dead zones, backlash, quantization, and asymmetric bounds. The amplitude-limited model is the simplest setting in which the fundamental adaptive-control difficulty is already visible.

Consider the scalar plant

\[ \dot{x}=a x+b u_a, \qquad b\neq 0, \]

and the stable reference model

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

We use the direct adaptive control command

\[ u_c=\hat{k}_x x+\hat{k}_r r =\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}, \qquad \hat{\boldsymbol{\theta}} = \begin{bmatrix} \hat{k}_x\\ \hat{k}_r \end{bmatrix}, \qquad \boldsymbol{\phi} = \begin{bmatrix} x\\ r \end{bmatrix}. \]

If the actuator were unconstrained, exact model matching would be possible whenever constant ideal gains \(k_x^{*}\) and \(k_r^{*}\) satisfy

\[ a+b k_x^{*}=a_m, \qquad b k_r^{*}=b_m. \]

Therefore,

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

Define the tracking error and parameter errors as

\[ e=x-x_m, \qquad \tilde{k}_x=\hat{k}_x-k_x^{*}, \qquad \tilde{k}_r=\hat{k}_r-k_r^{*}, \qquad \tilde{\boldsymbol{\theta}} =\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}. \]

Substitute \(u_a=u_c+\Delta u\) into the plant:

\[ \dot{x} = a x+b\hat{k}_x x+b\hat{k}_r r+b\Delta u. \]

Add and subtract the ideal gains and use the matching conditions:

\[ \begin{aligned} \dot{x} &= (a+b k_x^{*})x+b k_r^{*}r +b\tilde{k}_x x+b\tilde{k}_r r+b\Delta u\\ &= a_m x+b_m r +b\tilde{k}_x x+b\tilde{k}_r r+b\Delta u. \end{aligned} \]

Since \(\dot{x}_m=a_m x_m+b_m r\), subtraction gives

\[ \boxed{ \dot{e} = a_m e +b\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +b\Delta u }. \]

The unconstrained MRAC error model contains only the first two terms. Saturation contributes the additional term \(b\Delta u\). This term is not a small modeling detail: it is the exact mathematical location at which the nominal adaptive proof changes.

Assume the sign of \(b\) is known. Let \(\gamma>0\) and use the standard scalar update law

\[ \dot{\hat{\boldsymbol{\theta}}} = -\gamma\,\operatorname{sgn}(b)\, \boldsymbol{\phi}e. \]

Consider the Lyapunov candidate

\[ V = \frac{1}{2}e^2 + \frac{|b|}{2\gamma} \tilde{\boldsymbol{\theta}}^{T} \tilde{\boldsymbol{\theta}}. \]

Because the ideal parameters are constant, \(\dot{\tilde{\boldsymbol{\theta}}} =\dot{\hat{\boldsymbol{\theta}}}\). Hence

\[ \begin{aligned} \dot{V} &= e\dot{e} + \frac{|b|}{\gamma} \tilde{\boldsymbol{\theta}}^{T} \dot{\hat{\boldsymbol{\theta}}}\\ &= a_m e^2 +b e\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +b e\Delta u - |b|\operatorname{sgn}(b) \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}e\\ &= a_m e^2+b e\Delta u. \end{aligned} \]

Writing \(\lambda=-a_m>0\),

\[ \boxed{ \dot{V} = -\lambda e^2+b e\Delta u }. \]

In the unsaturated region, \(\Delta u=0\) and therefore \(\dot{V}=-\lambda e^2\leq 0\). During saturation, however, the product \(b e\Delta u\) can be positive. The standard update law cancels the parameter-error cross term, but it cannot cancel a control action that the actuator did not deliver.

5.1 A useful bound

Young's inequality gives, for any \(\varepsilon>0\),

\[ |b e\Delta u| \leq \frac{\varepsilon}{2}e^2 + \frac{b^2}{2\varepsilon}(\Delta u)^2. \]

Selecting \(\varepsilon=\lambda\) yields

\[ \dot{V} \leq -\frac{\lambda}{2}e^2 + \frac{b^2}{2\lambda}(\Delta u)^2. \]

This resembles an input-to-state bound with \(\Delta u\) as an input. Nevertheless, one must not silently assume that \(\Delta u\) is uniformly small: while \(u_a\) is bounded by \(u_{\max}\), the requested signal \(u_c=\hat{\boldsymbol{\theta}}^T\boldsymbol{\phi}\) can grow with the adaptive parameters. Consequently, \(|\Delta u|=|u_a-u_c|\) can also grow.

No adaptive law can force a physical actuator to produce more authority than its limit. Before studying parameter convergence, the desired reference motion must therefore be checked for actuator feasibility.

6.1 Pointwise ideal-control condition

Along an exactly matched motion \(x=x_m\), the required ideal command is

\[ u^{*}(t) = k_x^{*}x_m(t)+k_r^{*}r(t). \]

A necessary pointwise condition for exact model following is

\[ |u^{*}(t)|\leq u_{\max} \qquad \text{for all relevant }t. \]

6.2 Constant-reference steady-state condition

For a constant reference \(r(t)=r_0\), the stable reference model has equilibrium

\[ x_{m,\mathrm{ss}} = -\frac{b_m}{a_m}r_0. \]

If exact steady-state tracking is achieved, the plant equilibrium must satisfy

\[ 0=a x_{m,\mathrm{ss}}+b u_{\mathrm{ss}}^{*}, \]

and therefore

\[ u_{\mathrm{ss}}^{*} = -\frac{a}{b}x_{m,\mathrm{ss}} = \frac{a b_m}{a_m b}r_0. \]

A necessary steady-state feasibility condition is

\[ \boxed{ \left| \frac{a b_m}{a_m b}r_0 \right| \leq u_{\max} }. \]

This condition is necessary but not sufficient for the complete transient: even when the final equilibrium is feasible, a fast reference model may require a transient input larger than the available limit.

The standard update law treats tracking error as evidence that the controller parameters are incorrect. Under saturation, however, part or all of the error may be caused by missing actuator authority rather than parameter mismatch. If the command is infeasible, the tracking error cannot vanish; the update law can then continue integrating a persistent product of regressor and error.

\[ \dot{\hat{\boldsymbol{\theta}}} = -\gamma\,\operatorname{sgn}(b)\, \boldsymbol{\phi}e. \]

If, over a time interval, a component \(\phi_i e\) has a nonzero mean, then

\[ \hat{\theta}_i(t) = \hat{\theta}_i(t_0) - \gamma\operatorname{sgn}(b) \int_{t_0}^{t}\phi_i(\tau)e(\tau)\,d\tau. \]

A nonzero average integrand produces continuing parameter motion. The actual input remains clipped, but the internal command can become very large. This phenomenon is commonly described as adaptive parameter windup or parameter drift induced by saturation.

flowchart TD
  A["Reference model requests rapid or large motion"] --> B["Controller command exceeds actuator limit"]
  B --> C["Actual input differs from commanded input"]
  C --> D["Tracking error persists"]
  D --> E["Unmodified adaptive law \nkeeps changing gains"]
  E --> F["Internal command magnitude grows"]
  F --> B
  D --> G["Nominal Lyapunov cancellation \nis incomplete"]
        

7.1 Loss of parameter sensitivity

Away from the two saturation corners, the local slope of the actuator map is

\[ \frac{\partial u_a}{\partial u_c} = \begin{cases} 1, & |u_c|<u_{\max},\\ 0, & |u_c|>u_{\max}. \end{cases} \]

Since \(u_c=\hat{\boldsymbol{\theta}}^T\boldsymbol{\phi}\),

\[ \frac{\partial u_a} {\partial \hat{\boldsymbol{\theta}}} = \frac{\partial u_a}{\partial u_c} \boldsymbol{\phi}. \]

Deep in saturation this local sensitivity is zero: changing the adaptive parameters does not change the plant input. An update law derived from the unsaturated sensitivity therefore continues moving parameters in a direction that has no immediate effect on the plant.

8.1 Three regimes

1. No saturation: \(|u_c|\leq u_{\max}\) and \(\Delta u=0\). The nominal adaptive proof applies.
2. Transient or feasible saturation: the command clips temporarily, but the desired steady motion is attainable. Tracking can recover, although the final adaptive parameters may be far from ideal because the transient update was distorted.
3. Infeasible sustained saturation: the desired motion requires more actuator authority than available. Persistent error and unbounded internal commands can result under an unmodified law.

8.2 Saturation indicator

Define

\[ \chi_{\mathrm{sat}}(t) = \begin{cases} 1, & |u_c(t)|>u_{\max},\\ 0, & |u_c(t)|\leq u_{\max}. \end{cases} \]

Over a simulation interval \([0,T]\), the saturation duty fraction is

\[ \rho_{\mathrm{sat}} = \frac{1}{T} \int_0^T \chi_{\mathrm{sat}}(t)\,dt. \]

8.3 Command amplification ratio

A simple measure of hidden controller windup is

\[ \kappa_u = \frac{\max_{t\in[0,T]}|u_c(t)|} {u_{\max}}. \]

Values much larger than one show that the bounded actuator signal hides a much larger internal command.

8.4 Parameter-energy diagnostic

For \(\hat{\boldsymbol{\theta}} =[\hat{k}_x,\hat{k}_r]^T\), monitor

\[ E_{\theta}(t) = \|\hat{\boldsymbol{\theta}}(t)\|_2^2 = \hat{k}_x^2(t)+\hat{k}_r^2(t). \]

A growing \(E_\theta\) while \(u_a\) remains fixed at its limit is a clear numerical signature of adaptive windup.

The implementations use

\[ a=-0.3,\quad b=1,\quad a_m=-1.2,\quad b_m=1.2,\quad \gamma=8. \]

The ideal unconstrained gains are

\[ k_x^{*} = \frac{-1.2-(-0.3)}{1} =-0.9, \qquad k_r^{*} = \frac{1.2}{1} =1.2. \]

The reference is a step applied at \(t=1\) second. Three cases are compared:

1. Practically unsaturated: \(u_{\max}=10\), \(r_0=1\).
2. Feasible but saturated: \(u_{\max}=0.45\), \(r_0=1\). The required steady input is \(u_{\mathrm{ss}}^{*}=0.3\), so the final equilibrium is feasible even though transient saturation occurs.
3. Infeasible under saturation: \(u_{\max}=0.45\), \(r_0=2\). The required steady input is \(u_{\mathrm{ss}}^{*}=0.6\), which exceeds the available actuator magnitude.

All programs implement the same fixed-step fourth-order Runge-Kutta method. They save time histories of \(x,x_m,e,u_c,u_a,\hat{k}_x,\hat{k}_r\). The expected qualitative result is that the feasible saturated case eventually tracks but can finish with nonideal gains, whereas the infeasible case retains error and exhibits rapidly increasing parameter and command magnitudes.

The Python implementation uses NumPy for numerical arrays and Matplotlib for plots. SciPy.solve_ivp can replace the included RK4 routine, and the python-control package is useful for complementary linear-system analysis; neither is required by this file.

Chapter19_Lesson1.py

#!/usr/bin/env python3
"""
Chapter19_Lesson1.py
Effects of actuator saturation on an unmodified first-order MRAC law.

Three cases are simulated:
1. practically unsaturated actuator;
2. saturated but steady-state-feasible command;
3. saturated and steady-state-infeasible command.

Dependencies:
    pip install numpy matplotlib
"""

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 Case:
    name: str
    file_tag: str
    u_max: float
    r_amplitude: float


@dataclass
class Result:
    time: np.ndarray
    x: np.ndarray
    x_m: np.ndarray
    error: np.ndarray
    u_command: np.ndarray
    u_actual: np.ndarray
    k_x: np.ndarray
    k_r: np.ndarray
    saturation_fraction: float


A_PLANT = -0.3
B_PLANT = 1.0
A_MODEL = -1.2
B_MODEL = 1.2
GAMMA = 8.0
STEP_TIME = 1.0
DT = 0.002
T_FINAL = 25.0


def reference(t: float, amplitude: float) -> float:
    """Unit-step family used by all experiments."""
    return amplitude if t >= STEP_TIME else 0.0


def saturate(value: float, limit: float) -> float:
    return float(np.clip(value, -limit, limit))


def rhs(t: float, state: np.ndarray, case: Case) -> np.ndarray:
    """Closed-loop differential equations with the *unmodified* adaptive law."""
    x, x_m, k_x, k_r = state
    r = reference(t, case.r_amplitude)
    u_c = k_x * x + k_r * r
    u_a = saturate(u_c, case.u_max)
    e = x - x_m

    x_dot = A_PLANT * x + B_PLANT * u_a
    x_m_dot = A_MODEL * x_m + B_MODEL * r

    # Standard scalar MRAC update, derived under the assumption u_a == u_c.
    k_x_dot = -GAMMA * x * e
    k_r_dot = -GAMMA * r * e
    return np.array([x_dot, x_m_dot, k_x_dot, k_r_dot], dtype=float)


def rk4_step(
    function: Callable[[float, np.ndarray, Case], np.ndarray],
    t: float,
    state: np.ndarray,
    dt: float,
    case: Case,
) -> np.ndarray:
    k1 = function(t, state, case)
    k2 = function(t + 0.5 * dt, state + 0.5 * dt * k1, case)
    k3 = function(t + 0.5 * dt, state + 0.5 * dt * k2, case)
    k4 = function(t + dt, state + dt * k3, case)
    return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(case: Case) -> Result:
    time = np.arange(0.0, T_FINAL + DT, DT)
    state_history = np.zeros((time.size, 4), dtype=float)
    u_command = np.zeros(time.size, dtype=float)
    u_actual = np.zeros(time.size, dtype=float)

    state = np.zeros(4, dtype=float)
    for index, t in enumerate(time):
        state_history[index] = state
        x, _, k_x, k_r = state
        r = reference(t, case.r_amplitude)
        u_command[index] = k_x * x + k_r * r
        u_actual[index] = saturate(u_command[index], case.u_max)

        if index + 1 < time.size:
            state = rk4_step(rhs, t, state, DT, case)
            if not np.all(np.isfinite(state)):
                raise FloatingPointError(f"Non-finite state in case {case.name!r}.")

    x = state_history[:, 0]
    x_m = state_history[:, 1]
    error = x - x_m
    k_x = state_history[:, 2]
    k_r = state_history[:, 3]
    saturated = np.abs(u_command - u_actual) > 1.0e-10

    return Result(
        time=time,
        x=x,
        x_m=x_m,
        error=error,
        u_command=u_command,
        u_actual=u_actual,
        k_x=k_x,
        k_r=k_r,
        saturation_fraction=float(np.mean(saturated)),
    )


def save_csv(result: Result, case: Case, output_dir: Path) -> None:
    data = np.column_stack(
        (
            result.time,
            result.x,
            result.x_m,
            result.error,
            result.u_command,
            result.u_actual,
            result.k_x,
            result.k_r,
        )
    )
    header = "time,x,x_model,error,u_command,u_actual,k_x,k_r"
    np.savetxt(
        output_dir / f"Chapter19_Lesson1_{case.file_tag}.csv",
        data,
        delimiter=",",
        header=header,
        comments="",
    )


def print_summary(result: Result, case: Case) -> None:
    active = result.time >= STEP_TIME
    parameter_norm = np.hypot(result.k_x, result.k_r)
    required_steady_input = (
        -A_PLANT
        * (-B_MODEL / A_MODEL)
        * case.r_amplitude
        / B_PLANT
    )
    print(f"\nCase: {case.name}")
    print(f"  actuator limit             = {case.u_max: .4f}")
    print(f"  required steady input      = {required_steady_input: .4f}")
    print(f"  final tracking error       = {result.error[-1]: .6f}")
    print(f"  peak |tracking error|      = {np.max(np.abs(result.error[active])): .6f}")
    print(f"  saturation time fraction   = {result.saturation_fraction: .4f}")
    print(f"  peak |commanded input|     = {np.max(np.abs(result.u_command)): .6f}")
    print(f"  final gains [k_x, k_r]     = [{result.k_x[-1]: .6f}, {result.k_r[-1]: .6f}]")
    print(f"  peak parameter norm        = {np.max(parameter_norm): .6f}")


def make_plots(results: list[tuple[Case, Result]], output_dir: Path) -> None:
    figure, axes = plt.subplots(3, 1, figsize=(11, 12), sharex=True)

    for case, result in results:
        axes[0].plot(result.time, result.error, label=case.name)
        axes[1].plot(result.time, result.u_command, linestyle="--", label=f"{case.name}: u_c")
        axes[1].plot(result.time, result.u_actual, label=f"{case.name}: u_a")
        axes[2].plot(
            result.time,
            np.hypot(result.k_x, result.k_r),
            label=case.name,
        )

    axes[0].set_ylabel("e = x - x_m")
    axes[0].set_title("Tracking error")
    axes[0].grid(True)
    axes[0].legend()

    axes[1].set_ylabel("control input")
    axes[1].set_title("Commanded and actual actuator signals")
    axes[1].grid(True)
    axes[1].legend(ncol=2)

    axes[2].set_xlabel("time [s]")
    axes[2].set_ylabel("||theta_hat||_2")
    axes[2].set_title("Adaptive-parameter norm")
    axes[2].grid(True)
    axes[2].legend()

    figure.tight_layout()
    figure.savefig(output_dir / "Chapter19_Lesson1_results.png", dpi=180)
    plt.close(figure)


def main() -> None:
    output_dir = Path("Chapter19_Lesson1_output")
    output_dir.mkdir(parents=True, exist_ok=True)

    cases = [
        Case("unsaturated reference", "unsaturated", 10.0, 1.0),
        Case("feasible but saturated", "feasible_saturation", 0.45, 1.0),
        Case("infeasible under saturation", "infeasible_saturation", 0.45, 2.0),
    ]

    results: list[tuple[Case, Result]] = []
    for case in cases:
        result = simulate(case)
        results.append((case, result))
        save_csv(result, case, output_dir)
        print_summary(result, case)

    make_plots(results, output_dir)
    print(f"\nWrote CSV files and plot to: {output_dir.resolve()}")


if __name__ == "__main__":
    main()

This implementation uses only the C++17 standard library. For larger adaptive-control simulations, Eigen can provide matrix operations and Boost.Odeint can provide adaptive-step integration.

Chapter19_Lesson1.cpp

// Chapter19_Lesson1.cpp
// Effects of actuator saturation on an unmodified first-order MRAC law.
//
// Build:
//   g++ -std=c++17 -O2 -Wall -Wextra -pedantic Chapter19_Lesson1.cpp -o Chapter19_Lesson1
// Run:
//   ./Chapter19_Lesson1

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

namespace {

constexpr double aPlant = -0.3;
constexpr double bPlant = 1.0;
constexpr double aModel = -1.2;
constexpr double bModel = 1.2;
constexpr double gammaAdapt = 8.0;
constexpr double stepTime = 1.0;
constexpr double dt = 0.002;
constexpr double finalTime = 25.0;

struct Case {
    std::string name;
    std::string fileTag;
    double uMax;
    double referenceAmplitude;
};

struct Sample {
    double time{};
    double x{};
    double xModel{};
    double error{};
    double uCommand{};
    double uActual{};
    double kx{};
    double kr{};
};

using State = std::array<double, 4>;

double reference(const double time, const double amplitude) {
    return time >= stepTime ? amplitude : 0.0;
}

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

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

State rhs(const double time, const State& state, const Case& experiment) {
    const double x = state[0];
    const double xModel = state[1];
    const double kx = state[2];
    const double kr = state[3];

    const double r = reference(time, experiment.referenceAmplitude);
    const double uCommand = kx * x + kr * r;
    const double uActual = saturate(uCommand, experiment.uMax);
    const double error = x - xModel;

    return {
        aPlant * x + bPlant * uActual,
        aModel * xModel + bModel * r,
        -gammaAdapt * x * error,
        -gammaAdapt * r * error
    };
}

State rk4Step(const double time, const State& state, const Case& experiment) {
    const State k1 = rhs(time, state, experiment);
    const State k2 = rhs(time + 0.5 * dt, addScaled(state, k1, 0.5 * dt), experiment);
    const State k3 = rhs(time + 0.5 * dt, addScaled(state, k2, 0.5 * dt), experiment);
    const State k4 = rhs(time + dt, addScaled(state, k3, dt), experiment);

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

std::vector<Sample> simulate(const Case& experiment) {
    const std::size_t sampleCount =
        static_cast<std::size_t>(std::llround(finalTime / dt)) + 1U;

    std::vector<Sample> samples;
    samples.reserve(sampleCount);
    State state{0.0, 0.0, 0.0, 0.0};

    for (std::size_t index = 0; index < sampleCount; ++index) {
        const double time = static_cast<double>(index) * dt;
        const double r = reference(time, experiment.referenceAmplitude);
        const double uCommand = state[2] * state[0] + state[3] * r;
        const double uActual = saturate(uCommand, experiment.uMax);

        samples.push_back({
            time,
            state[0],
            state[1],
            state[0] - state[1],
            uCommand,
            uActual,
            state[2],
            state[3]
        });

        if (index + 1U < sampleCount) {
            state = rk4Step(time, state, experiment);
        }
    }
    return samples;
}

void saveCsv(
    const std::vector<Sample>& samples,
    const Case& experiment,
    const std::filesystem::path& outputDirectory
) {
    const auto path =
        outputDirectory / ("Chapter19_Lesson1_" + experiment.fileTag + ".csv");
    std::ofstream stream(path);
    if (!stream) {
        throw std::runtime_error("Cannot open output file: " + path.string());
    }

    stream << "time,x,x_model,error,u_command,u_actual,k_x,k_r\n";
    stream << std::setprecision(12);
    for (const auto& sample : samples) {
        stream
            << sample.time << ','
            << sample.x << ','
            << sample.xModel << ','
            << sample.error << ','
            << sample.uCommand << ','
            << sample.uActual << ','
            << sample.kx << ','
            << sample.kr << '\n';
    }
}

void printSummary(const std::vector<Sample>& samples, const Case& experiment) {
    double peakError = 0.0;
    double peakCommand = 0.0;
    double peakParameterNorm = 0.0;
    std::size_t saturatedCount = 0U;

    for (const auto& sample : samples) {
        if (sample.time >= stepTime) {
            peakError = std::max(peakError, std::abs(sample.error));
        }
        peakCommand = std::max(peakCommand, std::abs(sample.uCommand));
        peakParameterNorm = std::max(
            peakParameterNorm,
            std::hypot(sample.kx, sample.kr)
        );
        if (std::abs(sample.uCommand - sample.uActual) > 1.0e-10) {
            ++saturatedCount;
        }
    }

    const auto& last = samples.back();
    const double requiredSteadyInput =
        -aPlant * (-bModel / aModel) * experiment.referenceAmplitude / bPlant;
    const double saturationFraction =
        static_cast<double>(saturatedCount) / static_cast<double>(samples.size());

    std::cout << "\nCase: " << experiment.name << '\n'
              << "  actuator limit             = " << experiment.uMax << '\n'
              << "  required steady input      = " << requiredSteadyInput << '\n'
              << "  final tracking error       = " << last.error << '\n'
              << "  peak |tracking error|      = " << peakError << '\n'
              << "  saturation time fraction   = " << saturationFraction << '\n'
              << "  peak |commanded input|     = " << peakCommand << '\n'
              << "  final gains [k_x, k_r]     = ["
              << last.kx << ", " << last.kr << "]\n"
              << "  peak parameter norm        = " << peakParameterNorm << '\n';
}

}  // namespace

int main() {
    try {
        const std::filesystem::path outputDirectory{"Chapter19_Lesson1_output_cpp"};
        std::filesystem::create_directories(outputDirectory);

        const std::vector<Case> cases{
            {"unsaturated reference", "unsaturated", 10.0, 1.0},
            {"feasible but saturated", "feasible_saturation", 0.45, 1.0},
            {"infeasible under saturation", "infeasible_saturation", 0.45, 2.0}
        };

        for (const auto& experiment : cases) {
            const auto samples = simulate(experiment);
            saveCsv(samples, experiment, outputDirectory);
            printSummary(samples, experiment);
        }

        std::cout << "\nWrote CSV files to: "
                  << std::filesystem::absolute(outputDirectory) << '\n';
        return 0;
    } catch (const std::exception& exception) {
        std::cerr << "Error: " << exception.what() << '\n';
        return 1;
    }
}

The Java version uses the standard library and Java records. Apache Commons Math is a suitable alternative when reusable ODE solvers and matrix decompositions are required.

Chapter19_Lesson1.java

// Chapter19_Lesson1.java
// Effects of actuator saturation on an unmodified first-order MRAC law.
//
// Build:
//   javac Chapter19_Lesson1.java
// Run:
//   java Chapter19_Lesson1

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.nio.file.Paths;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter19_Lesson1 {
    private static final double A_PLANT = -0.3;
    private static final double B_PLANT = 1.0;
    private static final double A_MODEL = -1.2;
    private static final double B_MODEL = 1.2;
    private static final double GAMMA = 8.0;
    private static final double STEP_TIME = 1.0;
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 25.0;

    private Chapter19_Lesson1() {
    }

    private record ExperimentCase(
        String name,
        String fileTag,
        double uMax,
        double referenceAmplitude
    ) {
    }

    private record Sample(
        double time,
        double x,
        double xModel,
        double error,
        double uCommand,
        double uActual,
        double kx,
        double kr
    ) {
    }

    private static double reference(double time, double amplitude) {
        return time >= STEP_TIME ? amplitude : 0.0;
    }

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

    private static double[] rhs(
        double time,
        double[] state,
        ExperimentCase experiment
    ) {
        double x = state[0];
        double xModel = state[1];
        double kx = state[2];
        double kr = state[3];

        double r = reference(time, experiment.referenceAmplitude());
        double uCommand = kx * x + kr * r;
        double uActual = saturate(uCommand, experiment.uMax());
        double error = x - xModel;

        return new double[] {
            A_PLANT * x + B_PLANT * uActual,
            A_MODEL * xModel + B_MODEL * r,
            -GAMMA * x * error,
            -GAMMA * r * error
        };
    }

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

    private static double[] rk4Step(
        double time,
        double[] state,
        ExperimentCase experiment
    ) {
        double[] k1 = rhs(time, state, experiment);
        double[] k2 = rhs(
            time + 0.5 * DT,
            addScaled(state, k1, 0.5 * DT),
            experiment
        );
        double[] k3 = rhs(
            time + 0.5 * DT,
            addScaled(state, k2, 0.5 * DT),
            experiment
        );
        double[] k4 = rhs(
            time + DT,
            addScaled(state, k3, DT),
            experiment
        );

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

    private static List<Sample> simulate(ExperimentCase experiment) {
        int sampleCount = (int) Math.round(FINAL_TIME / DT) + 1;
        List<Sample> samples = new ArrayList<>(sampleCount);
        double[] state = {0.0, 0.0, 0.0, 0.0};

        for (int index = 0; index < sampleCount; ++index) {
            double time = index * DT;
            double r = reference(time, experiment.referenceAmplitude());
            double uCommand = state[2] * state[0] + state[3] * r;
            double uActual = saturate(uCommand, experiment.uMax());

            samples.add(new Sample(
                time,
                state[0],
                state[1],
                state[0] - state[1],
                uCommand,
                uActual,
                state[2],
                state[3]
            ));

            if (index + 1 < sampleCount) {
                state = rk4Step(time, state, experiment);
            }
        }
        return samples;
    }

    private static void saveCsv(
        List<Sample> samples,
        ExperimentCase experiment,
        Path outputDirectory
    ) throws IOException {
        Path path = outputDirectory.resolve(
            "Chapter19_Lesson1_" + experiment.fileTag() + ".csv"
        );

        try (BufferedWriter writer = Files.newBufferedWriter(
            path,
            StandardCharsets.UTF_8
        )) {
            writer.write("time,x,x_model,error,u_command,u_actual,k_x,k_r");
            writer.newLine();
            for (Sample sample : samples) {
                writer.write(String.format(
                    Locale.ROOT,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    sample.time(),
                    sample.x(),
                    sample.xModel(),
                    sample.error(),
                    sample.uCommand(),
                    sample.uActual(),
                    sample.kx(),
                    sample.kr()
                ));
            }
        }
    }

    private static void printSummary(
        List<Sample> samples,
        ExperimentCase experiment
    ) {
        double peakError = 0.0;
        double peakCommand = 0.0;
        double peakParameterNorm = 0.0;
        int saturatedCount = 0;

        for (Sample sample : samples) {
            if (sample.time() >= STEP_TIME) {
                peakError = Math.max(peakError, Math.abs(sample.error()));
            }
            peakCommand = Math.max(peakCommand, Math.abs(sample.uCommand()));
            peakParameterNorm = Math.max(
                peakParameterNorm,
                Math.hypot(sample.kx(), sample.kr())
            );
            if (Math.abs(sample.uCommand() - sample.uActual()) > 1.0e-10) {
                ++saturatedCount;
            }
        }

        Sample last = samples.get(samples.size() - 1);
        double requiredSteadyInput =
            -A_PLANT * (-B_MODEL / A_MODEL)
                * experiment.referenceAmplitude() / B_PLANT;
        double saturationFraction =
            (double) saturatedCount / (double) samples.size();

        System.out.printf(Locale.ROOT, "%nCase: %s%n", experiment.name());
        System.out.printf(
            Locale.ROOT,
            "  actuator limit             = % .6f%n",
            experiment.uMax()
        );
        System.out.printf(
            Locale.ROOT,
            "  required steady input      = % .6f%n",
            requiredSteadyInput
        );
        System.out.printf(
            Locale.ROOT,
            "  final tracking error       = % .6f%n",
            last.error()
        );
        System.out.printf(
            Locale.ROOT,
            "  peak |tracking error|      = % .6f%n",
            peakError
        );
        System.out.printf(
            Locale.ROOT,
            "  saturation time fraction   = % .6f%n",
            saturationFraction
        );
        System.out.printf(
            Locale.ROOT,
            "  peak |commanded input|     = % .6f%n",
            peakCommand
        );
        System.out.printf(
            Locale.ROOT,
            "  final gains [k_x, k_r]     = [% .6f, % .6f]%n",
            last.kx(),
            last.kr()
        );
        System.out.printf(
            Locale.ROOT,
            "  peak parameter norm        = % .6f%n",
            peakParameterNorm
        );
    }

    public static void main(String[] args) {
        Path outputDirectory = Paths.get("Chapter19_Lesson1_output_java");

        List<ExperimentCase> cases = List.of(
            new ExperimentCase(
                "unsaturated reference",
                "unsaturated",
                10.0,
                1.0
            ),
            new ExperimentCase(
                "feasible but saturated",
                "feasible_saturation",
                0.45,
                1.0
            ),
            new ExperimentCase(
                "infeasible under saturation",
                "infeasible_saturation",
                0.45,
                2.0
            )
        );

        try {
            Files.createDirectories(outputDirectory);
            for (ExperimentCase experiment : cases) {
                List<Sample> samples = simulate(experiment);
                saveCsv(samples, experiment, outputDirectory);
                printSummary(samples, experiment);
            }
            System.out.println(
                "\nWrote CSV files to: "
                    + outputDirectory.toAbsolutePath().normalize()
            );
        } catch (IOException | ArithmeticException exception) {
            System.err.println("Error: " + exception.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script uses a custom RK4 integrator, tables, and tiledlayout. It requires no specialized toolbox. A Simulink realization uses four Integrator blocks for \(x,x_m,\hat{k}_x,\hat{k}_r\), a Saturation block between \(u_c\) and the plant, Product blocks for \(xe\) and \(re\), and gains \(-\gamma\) in the adaptation channels.

Chapter19_Lesson1.m

% Chapter19_Lesson1.m
% Effects of actuator saturation on an unmodified first-order MRAC law.
%
% The script compares:
%   1) a practically unsaturated actuator,
%   2) a saturated but steady-state-feasible command,
%   3) a saturated and steady-state-infeasible command.
%
% No specialized toolbox is required. The integration is performed with a
% fixed-step fourth-order Runge-Kutta method so that the implementation
% parallels the Python, C++, and Java versions.

clear;
clc;
close all;

parameters.aPlant = -0.3;
parameters.bPlant = 1.0;
parameters.aModel = -1.2;
parameters.bModel = 1.2;
parameters.gamma = 8.0;
parameters.stepTime = 1.0;
parameters.dt = 0.002;
parameters.finalTime = 25.0;

cases = struct( ...
    'name', { ...
        'unsaturated reference', ...
        'feasible but saturated', ...
        'infeasible under saturation' ...
    }, ...
    'fileTag', { ...
        'unsaturated', ...
        'feasible_saturation', ...
        'infeasible_saturation' ...
    }, ...
    'uMax', {10.0, 0.45, 0.45}, ...
    'referenceAmplitude', {1.0, 1.0, 2.0} ...
);

outputDirectory = 'Chapter19_Lesson1_output_matlab';
if ~exist(outputDirectory, 'dir')
    mkdir(outputDirectory);
end

results = cell(size(cases));

for caseIndex = 1:numel(cases)
    results{caseIndex} = simulateCase(cases(caseIndex), parameters);
    saveCaseCsv(results{caseIndex}, cases(caseIndex), outputDirectory);
    printSummary(results{caseIndex}, cases(caseIndex), parameters);
end

figureHandle = figure('Name', 'Chapter 19 Lesson 1', 'Color', 'w');
layout = tiledlayout(3, 1, 'TileSpacing', 'compact', 'Padding', 'compact');

nexttile;
hold on;
grid on;
for caseIndex = 1:numel(cases)
    plot( ...
        results{caseIndex}.time, ...
        results{caseIndex}.error, ...
        'DisplayName', ...
        cases(caseIndex).name ...
    );
end
ylabel('e = x - x_m');
title('Tracking error');
legend('Location', 'best');

nexttile;
hold on;
grid on;
for caseIndex = 1:numel(cases)
    plot( ...
        results{caseIndex}.time, ...
        results{caseIndex}.uCommand, ...
        '--', ...
        'DisplayName', ...
        [cases(caseIndex).name, ': u_c'] ...
    );
    plot( ...
        results{caseIndex}.time, ...
        results{caseIndex}.uActual, ...
        'DisplayName', ...
        [cases(caseIndex).name, ': u_a'] ...
    );
end
ylabel('control input');
title('Commanded and actual actuator signals');
legend('Location', 'best');

nexttile;
hold on;
grid on;
for caseIndex = 1:numel(cases)
    parameterNorm = hypot( ...
        results{caseIndex}.kx, ...
        results{caseIndex}.kr ...
    );
    plot( ...
        results{caseIndex}.time, ...
        parameterNorm, ...
        'DisplayName', ...
        cases(caseIndex).name ...
    );
end
xlabel('time [s]');
ylabel('||thetaHat||_2');
title('Adaptive-parameter norm');
legend('Location', 'best');

title(layout, 'Effects of actuator saturation on unmodified MRAC');
exportgraphics( ...
    figureHandle, ...
    fullfile(outputDirectory, 'Chapter19_Lesson1_results.png'), ...
    'Resolution', ...
    180 ...
);

fprintf('\nWrote CSV files and plot to: %s\n', ...
    fullfile(pwd, outputDirectory));

% Simulink construction note:
% A corresponding block diagram uses four continuous Integrator blocks for
% x, x_m, k_x, and k_r; a Saturation block between u_c and the plant; Product
% blocks for x*e and r*e; and gains -gamma for the update laws. Logging both
% u_c and u_a reveals the control-deficiency signal Delta_u = u_a - u_c.

function result = simulateCase(experiment, parameters)
    time = (0:parameters.dt:parameters.finalTime).';
    sampleCount = numel(time);

    stateHistory = zeros(sampleCount, 4);
    uCommand = zeros(sampleCount, 1);
    uActual = zeros(sampleCount, 1);
    state = zeros(4, 1);

    for index = 1:sampleCount
        currentTime = time(index);
        stateHistory(index, :) = state.';

        referenceValue = referenceSignal( ...
            currentTime, ...
            experiment.referenceAmplitude, ...
            parameters.stepTime ...
        );
        uCommand(index) = state(3) * state(1) ...
            + state(4) * referenceValue;
        uActual(index) = saturate( ...
            uCommand(index), ...
            experiment.uMax ...
        );

        if index < sampleCount
            state = rk4Step( ...
                currentTime, ...
                state, ...
                experiment, ...
                parameters ...
            );
            if any(~isfinite(state))
                error('Non-finite state in case "%s".', experiment.name);
            end
        end
    end

    result.time = time;
    result.x = stateHistory(:, 1);
    result.xModel = stateHistory(:, 2);
    result.error = result.x - result.xModel;
    result.uCommand = uCommand;
    result.uActual = uActual;
    result.kx = stateHistory(:, 3);
    result.kr = stateHistory(:, 4);
    result.saturationFraction = mean( ...
        abs(uCommand - uActual) > 1.0e-10 ...
    );
end

function nextState = rk4Step(time, state, experiment, parameters)
    dt = parameters.dt;

    k1 = closedLoopRhs(time, state, experiment, parameters);
    k2 = closedLoopRhs( ...
        time + 0.5 * dt, ...
        state + 0.5 * dt * k1, ...
        experiment, ...
        parameters ...
    );
    k3 = closedLoopRhs( ...
        time + 0.5 * dt, ...
        state + 0.5 * dt * k2, ...
        experiment, ...
        parameters ...
    );
    k4 = closedLoopRhs( ...
        time + dt, ...
        state + dt * k3, ...
        experiment, ...
        parameters ...
    );

    nextState = state + (dt / 6.0) ...
        * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
end

function derivative = closedLoopRhs( ...
    time, ...
    state, ...
    experiment, ...
    parameters ...
)
    x = state(1);
    xModel = state(2);
    kx = state(3);
    kr = state(4);

    referenceValue = referenceSignal( ...
        time, ...
        experiment.referenceAmplitude, ...
        parameters.stepTime ...
    );
    uCommand = kx * x + kr * referenceValue;
    uActual = saturate(uCommand, experiment.uMax);
    trackingError = x - xModel;

    derivative = [ ...
        parameters.aPlant * x + parameters.bPlant * uActual; ...
        parameters.aModel * xModel + parameters.bModel * referenceValue; ...
        -parameters.gamma * x * trackingError; ...
        -parameters.gamma * referenceValue * trackingError ...
    ];
end

function value = referenceSignal(time, amplitude, stepTime)
    if time >= stepTime
        value = amplitude;
    else
        value = 0.0;
    end
end

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

function saveCaseCsv(result, experiment, outputDirectory)
    tableData = table( ...
        result.time, ...
        result.x, ...
        result.xModel, ...
        result.error, ...
        result.uCommand, ...
        result.uActual, ...
        result.kx, ...
        result.kr, ...
        'VariableNames', { ...
            'time', ...
            'x', ...
            'x_model', ...
            'error', ...
            'u_command', ...
            'u_actual', ...
            'k_x', ...
            'k_r' ...
        } ...
    );

    writetable( ...
        tableData, ...
        fullfile( ...
            outputDirectory, ...
            ['Chapter19_Lesson1_', experiment.fileTag, '.csv'] ...
        ) ...
    );
end

function printSummary(result, experiment, parameters)
    activeSamples = result.time >= parameters.stepTime;
    parameterNorm = hypot(result.kx, result.kr);
    requiredSteadyInput = ...
        -parameters.aPlant ...
        * (-parameters.bModel / parameters.aModel) ...
        * experiment.referenceAmplitude ...
        / parameters.bPlant;

    fprintf('\nCase: %s\n', experiment.name);
    fprintf('  actuator limit             = % .6f\n', experiment.uMax);
    fprintf('  required steady input      = % .6f\n', requiredSteadyInput);
    fprintf('  final tracking error       = % .6f\n', result.error(end));
    fprintf( ...
        '  peak |tracking error|      = % .6f\n', ...
        max(abs(result.error(activeSamples))) ...
    );
    fprintf( ...
        '  saturation time fraction   = % .6f\n', ...
        result.saturationFraction ...
    );
    fprintf( ...
        '  peak |commanded input|     = % .6f\n', ...
        max(abs(result.uCommand)) ...
    );
    fprintf( ...
        '  final gains [k_x, k_r]     = [% .6f, % .6f]\n', ...
        result.kx(end), ...
        result.kr(end) ...
    );
    fprintf( ...
        '  peak parameter norm        = % .6f\n', ...
        max(parameterNorm) ...
    );
end

The notebook expression implements the same RK4 experiment using native lists, associations, CSV export, and ListLinePlot.

Chapter19_Lesson1.nb

Notebook[{
 Cell[
  "Chapter 19, Lesson 1: Effects of Actuator Saturation on Adaptive Laws",
  "Title"
 ],
 Cell[
  BoxData[
   ToBoxes[
    Defer[
     ClearAll["Global`*"];

     parameters = <|
       "aPlant" -> -0.3,
       "bPlant" -> 1.0,
       "aModel" -> -1.2,
       "bModel" -> 1.2,
       "gamma" -> 8.0,
       "stepTime" -> 1.0,
       "dt" -> 0.002,
       "finalTime" -> 25.0
     |>;

     cases = {
       <|
         "name" -> "unsaturated reference",
         "fileTag" -> "unsaturated",
         "uMax" -> 10.0,
         "referenceAmplitude" -> 1.0
       |>,
       <|
         "name" -> "feasible but saturated",
         "fileTag" -> "feasible_saturation",
         "uMax" -> 0.45,
         "referenceAmplitude" -> 1.0
       |>,
       <|
         "name" -> "infeasible under saturation",
         "fileTag" -> "infeasible_saturation",
         "uMax" -> 0.45,
         "referenceAmplitude" -> 2.0
       |>
     };

     outputDirectory = "Chapter19_Lesson1_output_mathematica";
     If[!DirectoryQ[outputDirectory], CreateDirectory[outputDirectory]];

     referenceSignal[t_, amplitude_] :=
       If[t >= parameters["stepTime"], amplitude, 0.0];

     saturate[value_, limit_] := Clip[value, {-limit, limit}];

     closedLoopRhs[t_, state_List, experiment_Association] :=
       Module[
        {x, xModel, kx, kr, r, uCommand, uActual, error},
        {x, xModel, kx, kr} = state;
        r = referenceSignal[t, experiment["referenceAmplitude"]];
        uCommand = kx*x + kr*r;
        uActual = saturate[uCommand, experiment["uMax"]];
        error = x - xModel;
        {
          parameters["aPlant"]*x + parameters["bPlant"]*uActual,
          parameters["aModel"]*xModel + parameters["bModel"]*r,
          -parameters["gamma"]*x*error,
          -parameters["gamma"]*r*error
        }
       ];

     rk4Step[t_, state_List, experiment_Association] :=
       Module[{dt, k1, k2, k3, k4},
        dt = parameters["dt"];
        k1 = closedLoopRhs[t, state, experiment];
        k2 = closedLoopRhs[t + dt/2.0, state + dt*k1/2.0, experiment];
        k3 = closedLoopRhs[t + dt/2.0, state + dt*k2/2.0, experiment];
        k4 = closedLoopRhs[t + dt, state + dt*k3, experiment];
        state + dt*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0
       ];

     simulateCase[experiment_Association] :=
       Module[
        {
          time, sampleCount, states, uCommand, uActual, state,
          index, t, r, error, saturated
        },
        time = Range[
          0.0,
          parameters["finalTime"],
          parameters["dt"]
        ];
        sampleCount = Length[time];
        states = ConstantArray[0.0, {sampleCount, 4}];
        uCommand = ConstantArray[0.0, sampleCount];
        uActual = ConstantArray[0.0, sampleCount];
        state = {0.0, 0.0, 0.0, 0.0};

        Do[
         t = time[[index]];
         states[[index]] = state;
         r = referenceSignal[t, experiment["referenceAmplitude"]];
         uCommand[[index]] = state[[3]]*state[[1]] + state[[4]]*r;
         uActual[[index]] = saturate[
           uCommand[[index]],
           experiment["uMax"]
         ];
         If[
          index < sampleCount,
          state = rk4Step[t, state, experiment];
          If[
           !And @@ (FiniteQ /@ N[state]),
           Throw["Non-finite state in " <> experiment["name"]]
          ]
         ],
         {index, 1, sampleCount}
        ];

        error = states[[All, 1]] - states[[All, 2]];
        saturated = MapThread[
          Abs[#1 - #2] > 1.0*^-10 &,
          {uCommand, uActual}
        ];

        <|
          "time" -> time,
          "x" -> states[[All, 1]],
          "xModel" -> states[[All, 2]],
          "error" -> error,
          "uCommand" -> uCommand,
          "uActual" -> uActual,
          "kx" -> states[[All, 3]],
          "kr" -> states[[All, 4]],
          "saturationFraction" -> N[Count[saturated, True]/sampleCount]
        |>
       ];

     saveCaseCsv[result_Association, experiment_Association] :=
       Module[{header, rows, path},
        header = {
          "time", "x", "x_model", "error",
          "u_command", "u_actual", "k_x", "k_r"
        };
        rows = Transpose[{
          result["time"],
          result["x"],
          result["xModel"],
          result["error"],
          result["uCommand"],
          result["uActual"],
          result["kx"],
          result["kr"]
        }];
        path = FileNameJoin[{
          outputDirectory,
          "Chapter19_Lesson1_" <> experiment["fileTag"] <> ".csv"
        }];
        Export[path, Prepend[rows, header], "CSV"]
       ];

     printSummary[result_Association, experiment_Association] :=
       Module[
        {
          activeIndices, parameterNorm, requiredSteadyInput,
          finalIndex
        },
        activeIndices = Flatten @ Position[
          Map[# >= parameters["stepTime"] &, result["time"]],
          True
        ];
        parameterNorm = MapThread[
          Sqrt[#1^2 + #2^2] &,
          {result["kx"], result["kr"]}
        ];
        requiredSteadyInput =
          -parameters["aPlant"]
          *(-parameters["bModel"]/parameters["aModel"])
          *experiment["referenceAmplitude"]
          /parameters["bPlant"];
        finalIndex = Length[result["time"]];

        Print["\nCase: ", experiment["name"]];
        Print["  actuator limit             = ", N[experiment["uMax"], 8]];
        Print["  required steady input      = ", N[requiredSteadyInput, 8]];
        Print["  final tracking error       = ", N[result["error"][[finalIndex]], 8]];
        Print[
          "  peak |tracking error|      = ",
          N[Max[Abs[result["error"][[activeIndices]]]], 8]
        ];
        Print[
          "  saturation time fraction   = ",
          N[result["saturationFraction"], 8]
        ];
        Print[
          "  peak |commanded input|     = ",
          N[Max[Abs[result["uCommand"]]], 8]
        ];
        Print[
          "  final gains [k_x, k_r]     = ",
          N[{result["kx"][[-1]], result["kr"][[-1]]}, 8]
        ];
        Print[
          "  peak parameter norm        = ",
          N[Max[parameterNorm], 8]
        ]
       ];

     results = Association @ Map[
       Function[
        experiment,
        experiment["fileTag"] -> simulateCase[experiment]
       ],
       cases
     ];

     Scan[
       Function[
        experiment,
        saveCaseCsv[results[experiment["fileTag"]], experiment];
        printSummary[results[experiment["fileTag"]], experiment]
       ],
       cases
     ];

     trackingErrorPlot = ListLinePlot[
       Table[
        Transpose[{
          results[experiment["fileTag"]]["time"],
          results[experiment["fileTag"]]["error"]
        }],
        {experiment, cases}
       ],
       PlotLegends -> (Lookup[cases, "name"]),
       Frame -> True,
       FrameLabel -> {"time [s]", "e = x - x_m"},
       PlotLabel -> "Tracking error",
       ImageSize -> Large
     ];

     controlPlot = ListLinePlot[
       Flatten[
        Table[
         {
           Transpose[{
             results[experiment["fileTag"]]["time"],
             results[experiment["fileTag"]]["uCommand"]
           }],
           Transpose[{
             results[experiment["fileTag"]]["time"],
             results[experiment["fileTag"]]["uActual"]
           }]
         },
         {experiment, cases}
        ],
        1
       ],
       PlotLegends -> Flatten[
         Table[
          {
            experiment["name"] <> ": u_c",
            experiment["name"] <> ": u_a"
          },
          {experiment, cases}
         ]
       ],
       Frame -> True,
       FrameLabel -> {"time [s]", "control input"},
       PlotLabel -> "Commanded and actual actuator signals",
       ImageSize -> Large
     ];

     parameterPlot = ListLinePlot[
       Table[
        Transpose[{
          results[experiment["fileTag"]]["time"],
          MapThread[
           Sqrt[#1^2 + #2^2] &,
           {
             results[experiment["fileTag"]]["kx"],
             results[experiment["fileTag"]]["kr"]
           }
          ]
        }],
        {experiment, cases}
       ],
       PlotLegends -> (Lookup[cases, "name"]),
       Frame -> True,
       FrameLabel -> {"time [s]", "||thetaHat||_2"},
       PlotLabel -> "Adaptive-parameter norm",
       ImageSize -> Large
     ];

     combinedPlot = GraphicsColumn[
       {trackingErrorPlot, controlPlot, parameterPlot},
       Spacings -> 0.5
     ];

     Export[
       FileNameJoin[{
         outputDirectory,
         "Chapter19_Lesson1_results.png"
       }],
       combinedPlot,
       ImageResolution -> 180
     ];

     Print[
       "\nWrote CSV files and plot to: ",
       FileNameJoin[{Directory[], outputDirectory}]
     ];

     combinedPlot
    ]
   ]
  ],
  "Input"
 ]
},
WindowTitle -> "Chapter19_Lesson1",
Saveable -> True
]

15.1 Unsaturated case

The actuator reproduces the requested command, so \(\Delta u=0\). The nominal error dynamics and Lyapunov derivative are recovered. Tracking error tends to zero, although the adaptive gains need not converge exactly to \(\boldsymbol{\theta}^{*}\) without persistent excitation.

15.2 Feasible but saturated case

The initial adaptive command exceeds the actuator limit. During this interval, the Lyapunov derivative contains \(b e\Delta u\). The error can later vanish because the desired equilibrium is physically reachable. Nevertheless, the parameter values reached after the saturated transient can differ substantially from the ideal gains. Zero tracking error does not imply parameter correctness.

15.3 Infeasible case

The reference model asks for a steady input larger than \(u_{\max}\). The actual plant state cannot reach the reference-model equilibrium. The persistent error drives the standard update law, while the actuator remains clipped. Thus \(|u_c|\) and \(\|\hat{\boldsymbol{\theta}}\|\) can grow even though \(|u_a|\leq u_{\max}\).

15.4 Practical logging requirements

Logging only the physical actuator signal can conceal the problem. A useful implementation should record at least:

• the unsaturated command \(u_c\);
• the realized input \(u_a\);
• the deficiency \(\Delta u=u_a-u_c\);
• the tracking error \(e\);
• the parameter estimates and their norms;
• the saturation duty fraction.

Problem 1 (Error dynamics): For \(\dot{x}=a x+b\operatorname{sat}_{u_{\max}}(u_c)\), \(\dot{x}_m=a_mx_m+b_mr\), and \(u_c=\hat{k}_xx+\hat{k}_rr\), derive the tracking-error equation in terms of the parameter errors and \(\Delta u=u_a-u_c\).

Solution:

Write \(u_a=u_c+\Delta u\). Then

\[ \dot{x} = a x+b\hat{k}_x x+b\hat{k}_r r+b\Delta u. \]

Insert \(\hat{k}_x=k_x^{*}+\tilde{k}_x\) and \(\hat{k}_r=k_r^{*}+\tilde{k}_r\), and use \(a+b k_x^{*}=a_m\) and \(b k_r^{*}=b_m\):

\[ \dot{x} = a_m x+b_m r +b\tilde{k}_x x+b\tilde{k}_r r+b\Delta u. \]

Subtracting the reference-model equation gives

\[ \dot{e} = a_m e+b\tilde{k}_x x+b\tilde{k}_r r+b\Delta u. \]

Problem 2 (Feasibility test): Let \(a=-0.3\), \(b=1\), \(a_m=-1.2\), \(b_m=1.2\), and \(u_{\max}=0.45\). Determine whether constant commands \(r_0=1\) and \(r_0=2\) are steady-state feasible.

Solution:

The model equilibrium is

\[ x_{m,\mathrm{ss}} = -\frac{b_m}{a_m}r_0 = r_0. \]

The plant input required to hold that state is

\[ u_{\mathrm{ss}}^{*} = -\frac{a}{b}x_{m,\mathrm{ss}} = 0.3r_0. \]

For \(r_0=1\), \(|u_{\mathrm{ss}}^{*}|=0.3<0.45\); the equilibrium is feasible. For \(r_0=2\), \(|u_{\mathrm{ss}}^{*}|=0.6>0.45\); the equilibrium is infeasible.

Problem 3 (When can the Lyapunov function increase?): Suppose \(b>0\) and \(\dot{V}=-\lambda e^2+b e\Delta u\). Find a condition under which \(\dot{V}>0\).

Solution:

We require

\[ b e\Delta u>\lambda e^2. \]

For \(e\neq0\), this is equivalent to requiring \(e\Delta u>0\) and

\[ b|\Delta u|>\lambda|e|. \]

Thus the deficiency must have the same sign as the tracking error and a sufficiently large magnitude. The existence of such regions proves that the nominal Lyapunov function is not guaranteed to decrease under the unmodified law.

Problem 4 (Sensitivity collapse): Show why changing an adaptive parameter has no first-order effect on the actual actuator input when the actuator is strictly inside its saturated region.

Solution:

By the chain rule,

\[ \frac{\partial u_a}{\partial\hat{\theta}_i} = \frac{\partial u_a}{\partial u_c} \frac{\partial u_c}{\partial\hat{\theta}_i} = \frac{\partial u_a}{\partial u_c}\phi_i. \]

For \(|u_c|>u_{\max}\), the saturation map is locally constant, so \(\partial u_a/\partial u_c=0\). Hence

\[ \frac{\partial u_a}{\partial\hat{\theta}_i}=0. \]

Parameter updates can therefore enlarge the internal command without changing the physical input.

Problem 5 (Parameter drift under a constant residual error): Assume \(b>0\), a constant reference \(r_0\neq0\), and a persistent asymptotic tracking error \(e_{\infty}\neq0\). Under the standard law \(\dot{\hat{k}}_r=-\gamma r e\), determine the asymptotic growth of \(\hat{k}_r\).

Solution:

If \(r(t)=r_0\) and \(e(t)\) approaches the nonzero constant \(e_\infty\), then for large time

\[ \dot{\hat{k}}_r \approx -\gamma r_0 e_\infty. \]

Integration gives approximately

\[ \hat{k}_r(t) \approx \hat{k}_r(t_0) - \gamma r_0 e_\infty(t-t_0). \]

Therefore \(\hat{k}_r\) drifts linearly unless the average product \(r e\) becomes zero or the update law is modified.

Actuator saturation separates the controller command from the plant input. In first-order MRAC this creates the exact error-model term \(b\Delta u\) and changes the nominal Lyapunov derivative from \(-\lambda e^2\) to \(-\lambda e^2+b e\Delta u\). The additional term is not sign definite. When the reference motion is infeasible, persistent error can drive unmodified adaptive gains indefinitely even though the actual actuator input is bounded. A credible adaptive implementation must therefore check reference feasibility and log the commanded input, realized input, deficiency, error, and parameter norms. The next lesson develops update-law modifications intended to prevent or limit these effects.