In a fixed-gain PI or PID controller, actuator saturation can cause an integral state to continue accumulating even though the commanded input cannot be delivered. In an adaptive controller there is a second memory: the parameter estimate. A persistent tracking error created by saturation can therefore drive both the integral state and the adaptive parameters in directions that do not represent the plant uncertainty.

Let \(v(t)\) denote the unconstrained controller command and let \(u(t)\) denote the input actually applied by the actuator:

\[ u(t)=\operatorname{sat}_{\bar u}(v(t)),\qquad \Delta_u(t)=u(t)-v(t). \]

The mismatch \(\Delta_u\) is measurable whenever both the command and the actuator output are available. It is zero in the linear operating region. At the upper limit, \(v>\bar u\) implies \(\Delta_u<0\); at the lower limit, \(v<-\bar u\) implies \(\Delta_u>0\).

The central design principle is therefore: do not allow an adaptive law to interpret actuator deficiency as parametric plant error. Anti-windup feeds the known discrepancy \(\Delta_u\) into controller states, error models, or update laws so that their internal dynamics remain compatible with the constrained plant.

A useful architecture separates four signals: the nominal command \(v\), the applied input \(u\), the mismatch \(\Delta_u=u-v\), and a modified error \(\varepsilon\) used for adaptation.

flowchart TD
  R["Reference r"] --> RM["Reference model"]
  RM --> E["Tracking error e"]
  X["Plant state or output"] --> E
  E --> AC["Adaptive controller \nand integral state"]
  R --> AC
  X --> AC
  AC --> V["Requested input v"]
  V --> SAT["Actuator saturation"]
  SAT --> U["Applied input u"]
  U --> P["Uncertain plant"]
  P --> X
  U --> D["Mismatch: delta_u = u - v"]
  V --> D
  D --> BC["Back-calculation \nand mismatch model"]
  BC --> AC
  BC --> ME["Modified error epsilon"]
  E --> ME
  ME --> AL["Adaptive update law"]
  AL --> AC
        

The architecture does not make an infeasible reference feasible. Its purpose is to prevent hidden controller states and parameter estimates from growing merely because the actuator is unable to realize the nominal command.

Consider a SISO-input linear plant with state \(x\in\mathbb{R}^n\):

\[ \dot x=A x+b u,\qquad u=\operatorname{sat}_{\bar u}(v). \]

Let the stable reference model and the nominal adaptive command be

\[ \dot x_m=A_m x_m+b_m r,\qquad v=\hat\theta^{T}\phi, \]

where \(\phi\in\mathbb{R}^p\) is a regressor assembled from available states and reference signals. Assume an ideal parameter vector \(\theta^\ast\) exists such that the matching conditions hold. With \(\tilde\theta=\hat\theta-\theta^\ast\) and \(e=x-x_m\), the error system becomes

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

Choose \(P=P^T>0\) from

\[ A_m^T P+P A_m=-Q,\qquad Q=Q^T>0. \]

The ordinary gradient update law is

\[ \dot{\hat\theta}=-\Gamma\phi e^T P b,\qquad \Gamma=\Gamma^T>0. \]

For the standard Lyapunov function

\[ V=e^T P e+\tilde\theta^T\Gamma^{-1}\tilde\theta, \]

differentiation and cancellation of the parameter cross term give

\[ \dot V=-e^TQe+2e^TPb\Delta_u. \]

The final term has no fixed sign. Thus, the unconstrained MRAC proof no longer establishes monotonic decrease of \(V\). More importantly, a nonzero error caused solely by actuator saturation can continue driving \(\hat\theta\). This is adaptive windup or parameter drift.

4.1 Controller-State Back-Calculation

Suppose an integral state contributes positively to the command:

\[ v=v_a+k_I\eta,\qquad k_I>0, \]

where \(v_a\) is the nonintegral adaptive command. A back-calculation law is

\[ \dot\eta=e_c+k_{aw}\Delta_u,\qquad k_{aw}>0, \]

with servo error \(e_c=r-y\). If the upper actuator limit is active, then \(\Delta_u<0\); consequently, the correction term decreases \(\eta\) and discharges the integral contribution. The characteristic unwinding time is approximately \(T_{aw}=1/(k_I k_{aw})\) when other signals vary slowly.

A bounded implementation may combine back-calculation with tangent-cone projection:

\[ \dot\eta= \Pi_{[-\bar\eta,\bar\eta]} \left(\eta,e_c+k_{aw}\Delta_u\right). \]

Conditional integration is an alternative: suspend integration only when the actuator is saturated and the current error would push the command farther into saturation.

4.2 Adaptive-Error Conditioning

Back-calculation protects an integral state, but it does not by itself restore the MRAC Lyapunov cancellation. For the adaptive parameters, a mismatch model can separate the part of the tracking error caused by actuator deficiency from the part caused by parameter error.

Define an auxiliary state driven by the known actuator discrepancy:

\[ \dot x_{aw}=A_m x_{aw}+b\Delta_u,\qquad x_{aw}(0)=0. \]

The modified error is

\[ \varepsilon=e-x_{aw}. \]

Subtracting the auxiliary dynamics from the physical tracking-error dynamics gives

\[ \begin{aligned} \dot\varepsilon &=\dot e-\dot x_{aw}\\ &=A_m e+b\tilde\theta^T\phi+b\Delta_u -A_m x_{aw}-b\Delta_u\\ &=A_m\varepsilon+b\tilde\theta^T\phi. \end{aligned} \]

The saturation mismatch has disappeared from the modified-error dynamics. Adaptation can therefore use

\[ \dot{\hat\theta} =-\Gamma\phi\varepsilon^TPb. \]

Conceptually, \(x_{aw}\) is the response that the reference-model dynamics would produce from the missing actuator authority. Subtracting it prevents the adaptive law from learning the actuator limit as if it were an unknown plant coefficient.

Assume:

1. \(A_m\) is Hurwitz.
2. The matching parameter \(\theta^\ast\) exists.
3. The sign and direction needed in \(b\) are known.
4. The regressor \(\phi\) is bounded.
5. The mismatch model uses the same stable matrix and input direction as the error model.

Consider

\[ V=\varepsilon^TP\varepsilon+ \tilde\theta^T\Gamma^{-1}\tilde\theta. \]

Since \(\dot{\tilde\theta}=\dot{\hat\theta}\), substitution yields

\[ \begin{aligned} \dot V &=\varepsilon^T(A_m^TP+PA_m)\varepsilon +2\varepsilon^TPb\tilde\theta^T\phi +2\tilde\theta^T\Gamma^{-1}\dot{\hat\theta}\\ &=-\varepsilon^TQ\varepsilon +2\varepsilon^TPb\tilde\theta^T\phi -2\tilde\theta^T\phi\varepsilon^TPb\\ &=-\varepsilon^TQ\varepsilon\le0. \end{aligned} \]

Therefore \(\varepsilon\) and \(\tilde\theta\) are bounded and \(\varepsilon\in L_2\). If \(\dot\varepsilon\) is bounded, Barbalat's lemma gives \(\varepsilon(t)\to0\).

This conclusion concerns the modified error. The physical tracking error satisfies

\[ e=\varepsilon+x_{aw}. \]

Consequently, asymptotic tracking requires the actuator mismatch to vanish sufficiently so that \(x_{aw}(t)\to0\). If the reference remains infeasible and saturation persists, anti-windup cannot force \(e\to0\); it can only preserve bounded, interpretable internal behavior.

The ideal mismatch model may be imperfect because the true plant input gain is uncertain, the actuator has unmodeled rate limits, or the measured actuator output is noisy. Three modifications are commonly combined:

7.1 Projection

\[ \dot{\hat\theta} =\operatorname{Proj}_{\Omega} \left(\hat\theta,y_\theta\right), \qquad y_\theta=-\Gamma\phi\varepsilon^TPb. \]

For a convex admissible set \(\Omega\) containing \(\theta^\ast\), the projection operator is chosen to satisfy

\[ \tilde\theta^T\Gamma^{-1} \left[ \operatorname{Proj}_{\Omega}(\hat\theta,y_\theta)-y_\theta \right]\le0. \]

Thus projection does not spoil the Lyapunov inequality at the boundary.

7.2 Leakage

\[ y_\theta= -\Gamma\phi\varepsilon^TPb-\sigma\hat\theta, \qquad \sigma>0. \]

Leakage discourages parameter drift under small persistent disturbances. It generally replaces exact asymptotic parameter properties with uniform ultimate boundedness.

7.3 Saturation-Dependent Adaptation Gate

\[ g(\Delta_u)= \frac{1}{1+\left(|\Delta_u|/\delta_g\right)^2}, \qquad 0<g\le1, \]

\[ y_\theta= -\Gamma\phi\varepsilon^TPb\,g(\Delta_u) -\sigma\hat\theta. \]

The gate leaves adaptation almost unchanged near the unsaturated region and slows learning when the discrepancy is large. Because \(g\neq1\) changes the exact cancellation used in Section 6, its stability claim must be established as a robust or ultimate-boundedness result rather than copied from the ideal proof.

Consider

\[ \dot x=-a x+b u,\qquad \dot x_m=-a_mx_m+b_mr, \]

\[ v=\hat\theta_xx+\hat\theta_rr+k_I\eta,\qquad |u|\le\bar u. \]

Ignoring the integral augmentation in the matching calculation, the ideal gains are

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

For a stable plant with constant saturated input, the largest steady-state magnitude is

\[ |x_{ss}|\le\frac{|b|}{a}\bar u. \]

In the numerical lab, \(a=1.2\), \(b=0.8\), and \(\bar u=1.2\), so the equilibrium magnitude cannot exceed \(0.8\). References such as \(r=2.0\) and \(r=-1.8\) are intentionally infeasible. They are followed by feasible intervals to test whether the controller recovers without large stored integral or parameter error.

flowchart TD
  A["Compute command v"] --> B["Apply actuator limits to obtain u"]
  B --> C["Compute delta_u = u - v"]
  C --> D["Update mismatch model x_aw"]
  D --> E["Form epsilon = e - x_aw"]
  E --> F["Apply adaptation gate and leakage"]
  F --> G["Project parameter rate at bounds"]
  C --> H["Back-calculate integral state"]
  H --> I["Project integral state at bounds"]
  G --> J["Integrate adaptive parameters"]
  I --> K["Advance plant and reference model"]
  J --> K
  K --> L["Log command, input, mismatch, errors, and limits"]
        

A defensible implementation logs at least: \(v\), \(u\), \(\Delta_u\), physical and modified errors, parameter estimates, integrator states, projection activity, and the fraction of time spent saturated. Anti-windup should be assessed by recovery after saturation, not only by tracking error while the reference is impossible.

The examples use forward Euler integration to expose every update explicitly. Production software should use the sampling and numerical safeguards developed in later implementation chapters.

Required libraries are numpy and matplotlib. The script generates two CSV files and one comparison plot.

Chapter19_Lesson3.py

# Chapter19_Lesson3.py
# Adaptive anti-windup demonstration for a saturated first-order plant.
# Dependencies: numpy, matplotlib
#
# The script compares:
#   1) a naive adaptive controller whose integrator and parameters continue
#      changing during saturation;
#   2) a saturation-aware controller with back-calculation, a mismatch model,
#      adaptation gating, leakage, and projection.
#
# Plant:           x_dot  = -a*x + b*u
# Reference model: xm_dot = -am*xm + bm*r
# Command:         v      = theta_x*x + theta_r*r + ki*eta
# Applied input:   u      = sat(v, umax)
# Mismatch:        delta  = u - v

from __future__ import annotations

from dataclasses import dataclass
import csv
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 0.002
    horizon: float = 20.0
    a: float = 1.2
    b: float = 0.8
    am: float = 1.5
    bm: float = 1.5
    umax: float = 1.2
    gamma_x: float = 1.8
    gamma_r: float = 1.0
    ki: float = 0.35
    kaw: float = 4.0
    lambda_aw: float = 1.5
    sigma: float = 0.01
    gate_width: float = 0.15
    eta_limit: float = 4.0
    theta_x_min: float = -5.0
    theta_x_max: float = 2.0
    theta_r_min: float = -1.0
    theta_r_max: float = 5.0


def reference(t: float) -> float:
    """Reference schedule containing infeasible and feasible intervals."""
    if t < 0.5:
        return 0.0
    if t < 6.0:
        return 2.0       # infeasible: the actuator cannot sustain this state
    if t < 12.0:
        return 0.5       # feasible
    if t < 16.0:
        return -1.8      # infeasible in the negative direction
    return 0.3           # feasible recovery interval


def saturate(value: float, limit: float) -> float:
    return max(-limit, min(limit, value))


def projected_rate(
    theta: np.ndarray,
    raw_rate: np.ndarray,
    lower: np.ndarray,
    upper: np.ndarray,
) -> np.ndarray:
    """Tangent-cone projection for a box-constrained parameter vector."""
    rate = raw_rate.copy()
    for i in range(theta.size):
        if theta[i] <= lower[i] and rate[i] < 0.0:
            rate[i] = 0.0
        elif theta[i] >= upper[i] and rate[i] > 0.0:
            rate[i] = 0.0
    return rate


def simulate(saturation_aware: bool, cfg: Config) -> dict[str, np.ndarray]:
    steps = int(round(cfg.horizon / cfg.dt)) + 1
    time = np.linspace(0.0, cfg.horizon, steps)

    data = {
        "time": time,
        "reference": np.zeros(steps),
        "x": np.zeros(steps),
        "xm": np.zeros(steps),
        "tracking_error": np.zeros(steps),
        "modified_error": np.zeros(steps),
        "command": np.zeros(steps),
        "input": np.zeros(steps),
        "mismatch": np.zeros(steps),
        "integrator": np.zeros(steps),
        "theta_x": np.zeros(steps),
        "theta_r": np.zeros(steps),
        "gate": np.ones(steps),
    }

    x = 0.0
    xm = 0.0
    eta = 0.0
    x_aw = 0.0
    theta = np.array([0.0, 0.5], dtype=float)

    gamma = np.array([cfg.gamma_x, cfg.gamma_r], dtype=float)
    lower = np.array([cfg.theta_x_min, cfg.theta_r_min], dtype=float)
    upper = np.array([cfg.theta_x_max, cfg.theta_r_max], dtype=float)

    for k, t in enumerate(time):
        r = reference(float(t))
        e = x - xm
        servo_error = r - x
        phi = np.array([x, r], dtype=float)

        v = float(theta @ phi + cfg.ki * eta)
        u = saturate(v, cfg.umax)
        delta = u - v

        if saturation_aware:
            # The mismatch model estimates the state component generated by
            # applying u instead of the unconstrained command v.
            epsilon = e - x_aw

            # Smoothly reduce learning when the actuator discrepancy is large.
            gate = 1.0 / (1.0 + (abs(delta) / cfg.gate_width) ** 2)

            raw_theta_rate = -gamma * phi * epsilon * gate - cfg.sigma * theta
            theta_rate = projected_rate(theta, raw_theta_rate, lower, upper)

            # Back-calculation discharges the integral state when v and u differ.
            eta_rate = servo_error + cfg.kaw * delta
            if eta >= cfg.eta_limit and eta_rate > 0.0:
                eta_rate = 0.0
            elif eta <= -cfg.eta_limit and eta_rate < 0.0:
                eta_rate = 0.0

            x_aw_rate = -cfg.lambda_aw * x_aw + cfg.b * delta
        else:
            epsilon = e
            gate = 1.0
            theta_rate = -gamma * phi * e
            eta_rate = servo_error
            x_aw_rate = 0.0

        x_rate = -cfg.a * x + cfg.b * u
        xm_rate = -cfg.am * xm + cfg.bm * r

        data["reference"][k] = r
        data["x"][k] = x
        data["xm"][k] = xm
        data["tracking_error"][k] = e
        data["modified_error"][k] = epsilon
        data["command"][k] = v
        data["input"][k] = u
        data["mismatch"][k] = delta
        data["integrator"][k] = eta
        data["theta_x"][k] = theta[0]
        data["theta_r"][k] = theta[1]
        data["gate"][k] = gate

        x += cfg.dt * x_rate
        xm += cfg.dt * xm_rate
        eta += cfg.dt * eta_rate
        x_aw += cfg.dt * x_aw_rate
        theta += cfg.dt * theta_rate

    return data


def metrics(data: dict[str, np.ndarray]) -> dict[str, float]:
    time = data["time"]
    dt = time[1] - time[0]
    return {
        "IAE": float(np.sum(np.abs(data["tracking_error"])) * dt),
        "max_abs_command": float(np.max(np.abs(data["command"]))),
        "max_abs_integrator": float(np.max(np.abs(data["integrator"]))),
        "max_abs_parameter": float(
            np.max(np.abs(np.column_stack((data["theta_x"], data["theta_r"]))))
        ),
        "saturation_fraction": float(np.mean(np.abs(data["mismatch"]) > 1.0e-9)),
        "final_abs_tracking_error": float(abs(data["tracking_error"][-1])),
    }


def write_csv(path: Path, data: dict[str, np.ndarray]) -> None:
    columns = list(data.keys())
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(columns)
        for values in zip(*(data[column] for column in columns)):
            writer.writerow([f"{float(value):.10g}" for value in values])


def make_plot(
    baseline: dict[str, np.ndarray],
    aware: dict[str, np.ndarray],
    output_path: Path,
) -> None:
    t = baseline["time"]
    fig, axes = plt.subplots(4, 1, figsize=(10, 12), sharex=True)

    axes[0].plot(t, baseline["reference"], "--", label="reference")
    axes[0].plot(t, baseline["x"], label="naive x")
    axes[0].plot(t, aware["x"], label="anti-windup x")
    axes[0].set_ylabel("state")
    axes[0].grid(True)
    axes[0].legend()

    axes[1].plot(t, baseline["command"], label="naive command v")
    axes[1].plot(t, aware["command"], label="anti-windup command v")
    axes[1].plot(t, aware["input"], "--", label="applied input u")
    axes[1].set_ylabel("control")
    axes[1].grid(True)
    axes[1].legend()

    axes[2].plot(t, baseline["theta_x"], label="naive theta_x")
    axes[2].plot(t, baseline["theta_r"], label="naive theta_r")
    axes[2].plot(t, aware["theta_x"], "--", label="AW theta_x")
    axes[2].plot(t, aware["theta_r"], "--", label="AW theta_r")
    axes[2].set_ylabel("parameters")
    axes[2].grid(True)
    axes[2].legend()

    axes[3].plot(t, baseline["integrator"], label="naive eta")
    axes[3].plot(t, aware["integrator"], label="AW eta")
    axes[3].plot(t, aware["gate"], "--", label="adaptation gate")
    axes[3].set_xlabel("time [s]")
    axes[3].set_ylabel("internal states")
    axes[3].grid(True)
    axes[3].legend()

    fig.tight_layout()
    fig.savefig(output_path, dpi=180)
    plt.close(fig)


def main() -> None:
    cfg = Config()
    output_dir = Path(__file__).resolve().parent

    baseline = simulate(False, cfg)
    aware = simulate(True, cfg)

    write_csv(output_dir / "Chapter19_Lesson3_naive.csv", baseline)
    write_csv(output_dir / "Chapter19_Lesson3_anti_windup.csv", aware)
    make_plot(baseline, aware, output_dir / "Chapter19_Lesson3_comparison.png")

    print("Naive adaptive controller")
    for key, value in metrics(baseline).items():
        print(f"  {key:28s}: {value:.6f}")

    print("\nSaturation-aware adaptive controller")
    for key, value in metrics(aware).items():
        print(f"  {key:28s}: {value:.6f}")

    print("\nGenerated CSV files and Chapter19_Lesson3_comparison.png")


if __name__ == "__main__":
    main()

This version uses only the C++ standard library and writes comparable CSV outputs.

Chapter19_Lesson3.cpp

// Chapter19_Lesson3.cpp
// Saturated adaptive control with and without anti-windup compensation.
// Build: g++ -std=c++17 -O2 Chapter19_Lesson3.cpp -o Chapter19_Lesson3
// Run:   ./Chapter19_Lesson3

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

struct Config {
    double dt = 0.002;
    double horizon = 20.0;
    double a = 1.2;
    double b = 0.8;
    double am = 1.5;
    double bm = 1.5;
    double umax = 1.2;
    double gamma_x = 1.8;
    double gamma_r = 1.0;
    double ki = 0.35;
    double kaw = 4.0;
    double lambda_aw = 1.5;
    double sigma = 0.01;
    double gate_width = 0.15;
    double eta_limit = 4.0;
    double theta_x_min = -5.0;
    double theta_x_max = 2.0;
    double theta_r_min = -1.0;
    double theta_r_max = 5.0;
};

struct Row {
    double time;
    double reference;
    double x;
    double xm;
    double tracking_error;
    double modified_error;
    double command;
    double input;
    double mismatch;
    double integrator;
    double theta_x;
    double theta_r;
    double gate;
};

struct Metrics {
    double iae = 0.0;
    double max_abs_command = 0.0;
    double max_abs_integrator = 0.0;
    double max_abs_parameter = 0.0;
    double saturation_fraction = 0.0;
    double final_abs_tracking_error = 0.0;
};

double reference(double t) {
    if (t < 0.5) {
        return 0.0;
    }
    if (t < 6.0) {
        return 2.0;
    }
    if (t < 12.0) {
        return 0.5;
    }
    if (t < 16.0) {
        return -1.8;
    }
    return 0.3;
}

double saturate(double value, double limit) {
    return std::max(-limit, std::min(limit, value));
}

double projectedRate(
    double theta,
    double raw_rate,
    double lower,
    double upper
) {
    if (theta <= lower && raw_rate < 0.0) {
        return 0.0;
    }
    if (theta >= upper && raw_rate > 0.0) {
        return 0.0;
    }
    return raw_rate;
}

std::vector<Row> simulate(bool saturation_aware, const Config& cfg) {
    const std::size_t steps =
        static_cast<std::size_t>(std::llround(cfg.horizon / cfg.dt)) + 1U;
    std::vector<Row> rows;
    rows.reserve(steps);

    double x = 0.0;
    double xm = 0.0;
    double eta = 0.0;
    double x_aw = 0.0;
    double theta_x = 0.0;
    double theta_r = 0.5;

    for (std::size_t k = 0; k < steps; ++k) {
        const double t = static_cast<double>(k) * cfg.dt;
        const double r = reference(t);
        const double e = x - xm;
        const double servo_error = r - x;

        const double v = theta_x * x + theta_r * r + cfg.ki * eta;
        const double u = saturate(v, cfg.umax);
        const double delta = u - v;

        double epsilon = e;
        double gate = 1.0;
        double theta_x_rate = -cfg.gamma_x * x * e;
        double theta_r_rate = -cfg.gamma_r * r * e;
        double eta_rate = servo_error;
        double x_aw_rate = 0.0;

        if (saturation_aware) {
            epsilon = e - x_aw;
            gate = 1.0 /
                (1.0 + std::pow(std::abs(delta) / cfg.gate_width, 2.0));

            const double raw_theta_x_rate =
                -cfg.gamma_x * x * epsilon * gate - cfg.sigma * theta_x;
            const double raw_theta_r_rate =
                -cfg.gamma_r * r * epsilon * gate - cfg.sigma * theta_r;

            theta_x_rate = projectedRate(
                theta_x,
                raw_theta_x_rate,
                cfg.theta_x_min,
                cfg.theta_x_max
            );
            theta_r_rate = projectedRate(
                theta_r,
                raw_theta_r_rate,
                cfg.theta_r_min,
                cfg.theta_r_max
            );

            eta_rate = servo_error + cfg.kaw * delta;
            if (eta >= cfg.eta_limit && eta_rate > 0.0) {
                eta_rate = 0.0;
            } else if (eta <= -cfg.eta_limit && eta_rate < 0.0) {
                eta_rate = 0.0;
            }

            x_aw_rate = -cfg.lambda_aw * x_aw + cfg.b * delta;
        }

        rows.push_back(Row{
            t,
            r,
            x,
            xm,
            e,
            epsilon,
            v,
            u,
            delta,
            eta,
            theta_x,
            theta_r,
            gate
        });

        const double x_rate = -cfg.a * x + cfg.b * u;
        const double xm_rate = -cfg.am * xm + cfg.bm * r;

        x += cfg.dt * x_rate;
        xm += cfg.dt * xm_rate;
        eta += cfg.dt * eta_rate;
        x_aw += cfg.dt * x_aw_rate;
        theta_x += cfg.dt * theta_x_rate;
        theta_r += cfg.dt * theta_r_rate;

        if (!std::isfinite(x) || !std::isfinite(theta_x) ||
            !std::isfinite(theta_r) || !std::isfinite(eta)) {
            throw std::runtime_error("Non-finite simulation state detected.");
        }
    }

    return rows;
}

Metrics computeMetrics(const std::vector<Row>& rows, const Config& cfg) {
    Metrics result;
    std::size_t saturated_samples = 0U;

    for (const Row& row : rows) {
        result.iae += std::abs(row.tracking_error) * cfg.dt;
        result.max_abs_command =
            std::max(result.max_abs_command, std::abs(row.command));
        result.max_abs_integrator =
            std::max(result.max_abs_integrator, std::abs(row.integrator));
        result.max_abs_parameter = std::max(
            result.max_abs_parameter,
            std::max(std::abs(row.theta_x), std::abs(row.theta_r))
        );
        if (std::abs(row.mismatch) > 1.0e-9) {
            ++saturated_samples;
        }
    }

    result.saturation_fraction =
        static_cast<double>(saturated_samples) /
        static_cast<double>(rows.size());
    result.final_abs_tracking_error =
        std::abs(rows.back().tracking_error);
    return result;
}

void writeCsv(const std::string& filename, const std::vector<Row>& rows) {
    std::ofstream output(filename);
    if (!output) {
        throw std::runtime_error("Cannot open output file: " + filename);
    }

    output << "time,reference,x,xm,tracking_error,modified_error,"
              "command,input,mismatch,integrator,theta_x,theta_r,gate\n";
    output << std::setprecision(10);

    for (const Row& row : rows) {
        output
            << row.time << ','
            << row.reference << ','
            << row.x << ','
            << row.xm << ','
            << row.tracking_error << ','
            << row.modified_error << ','
            << row.command << ','
            << row.input << ','
            << row.mismatch << ','
            << row.integrator << ','
            << row.theta_x << ','
            << row.theta_r << ','
            << row.gate << '\n';
    }
}

void printMetrics(const std::string& title, const Metrics& metrics) {
    std::cout << title << '\n'
              << "  IAE                         : " << metrics.iae << '\n'
              << "  max_abs_command             : "
              << metrics.max_abs_command << '\n'
              << "  max_abs_integrator          : "
              << metrics.max_abs_integrator << '\n'
              << "  max_abs_parameter           : "
              << metrics.max_abs_parameter << '\n'
              << "  saturation_fraction         : "
              << metrics.saturation_fraction << '\n'
              << "  final_abs_tracking_error    : "
              << metrics.final_abs_tracking_error << "\n\n";
}

int main() {
    try {
        const Config cfg;
        const std::vector<Row> baseline = simulate(false, cfg);
        const std::vector<Row> aware = simulate(true, cfg);

        writeCsv("Chapter19_Lesson3_naive_cpp.csv", baseline);
        writeCsv("Chapter19_Lesson3_anti_windup_cpp.csv", aware);

        std::cout << std::fixed << std::setprecision(6);
        printMetrics(
            "Naive adaptive controller",
            computeMetrics(baseline, cfg)
        );
        printMetrics(
            "Saturation-aware adaptive controller",
            computeMetrics(aware, cfg)
        );

        std::cout << "Generated two CSV files.\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

The Java implementation uses no third-party dependency and is suitable for Java 11 or newer.

Chapter19_Lesson3.java

// Chapter19_Lesson3.java
// Saturated adaptive control with and without anti-windup compensation.
// Build: javac Chapter19_Lesson3.java
// Run:   java Chapter19_Lesson3

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

public final class Chapter19_Lesson3 {
    private Chapter19_Lesson3() {
    }

    private static final class Config {
        final double dt = 0.002;
        final double horizon = 20.0;
        final double a = 1.2;
        final double b = 0.8;
        final double am = 1.5;
        final double bm = 1.5;
        final double umax = 1.2;
        final double gammaX = 1.8;
        final double gammaR = 1.0;
        final double ki = 0.35;
        final double kaw = 4.0;
        final double lambdaAw = 1.5;
        final double sigma = 0.01;
        final double gateWidth = 0.15;
        final double etaLimit = 4.0;
        final double thetaXMin = -5.0;
        final double thetaXMax = 2.0;
        final double thetaRMin = -1.0;
        final double thetaRMax = 5.0;
    }

    private static final class Row {
        final double time;
        final double reference;
        final double x;
        final double xm;
        final double trackingError;
        final double modifiedError;
        final double command;
        final double input;
        final double mismatch;
        final double integrator;
        final double thetaX;
        final double thetaR;
        final double gate;

        Row(
            double time,
            double reference,
            double x,
            double xm,
            double trackingError,
            double modifiedError,
            double command,
            double input,
            double mismatch,
            double integrator,
            double thetaX,
            double thetaR,
            double gate
        ) {
            this.time = time;
            this.reference = reference;
            this.x = x;
            this.xm = xm;
            this.trackingError = trackingError;
            this.modifiedError = modifiedError;
            this.command = command;
            this.input = input;
            this.mismatch = mismatch;
            this.integrator = integrator;
            this.thetaX = thetaX;
            this.thetaR = thetaR;
            this.gate = gate;
        }
    }

    private static final class Metrics {
        double iae;
        double maxAbsCommand;
        double maxAbsIntegrator;
        double maxAbsParameter;
        double saturationFraction;
        double finalAbsTrackingError;
    }

    private static double reference(double t) {
        if (t < 0.5) {
            return 0.0;
        }
        if (t < 6.0) {
            return 2.0;
        }
        if (t < 12.0) {
            return 0.5;
        }
        if (t < 16.0) {
            return -1.8;
        }
        return 0.3;
    }

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

    private static double projectedRate(
        double theta,
        double rawRate,
        double lower,
        double upper
    ) {
        if (theta <= lower && rawRate < 0.0) {
            return 0.0;
        }
        if (theta >= upper && rawRate > 0.0) {
            return 0.0;
        }
        return rawRate;
    }

    private static List<Row> simulate(
        boolean saturationAware,
        Config cfg
    ) {
        int steps = (int) Math.round(cfg.horizon / cfg.dt) + 1;
        List<Row> rows = new ArrayList<>(steps);

        double x = 0.0;
        double xm = 0.0;
        double eta = 0.0;
        double xAw = 0.0;
        double thetaX = 0.0;
        double thetaR = 0.5;

        for (int k = 0; k < steps; k++) {
            double t = k * cfg.dt;
            double r = reference(t);
            double e = x - xm;
            double servoError = r - x;

            double v = thetaX * x + thetaR * r + cfg.ki * eta;
            double u = saturate(v, cfg.umax);
            double delta = u - v;

            double epsilon = e;
            double gate = 1.0;
            double thetaXRate = -cfg.gammaX * x * e;
            double thetaRRate = -cfg.gammaR * r * e;
            double etaRate = servoError;
            double xAwRate = 0.0;

            if (saturationAware) {
                epsilon = e - xAw;
                gate = 1.0 /
                    (1.0 + Math.pow(Math.abs(delta) / cfg.gateWidth, 2.0));

                double rawThetaXRate =
                    -cfg.gammaX * x * epsilon * gate
                    - cfg.sigma * thetaX;
                double rawThetaRRate =
                    -cfg.gammaR * r * epsilon * gate
                    - cfg.sigma * thetaR;

                thetaXRate = projectedRate(
                    thetaX,
                    rawThetaXRate,
                    cfg.thetaXMin,
                    cfg.thetaXMax
                );
                thetaRRate = projectedRate(
                    thetaR,
                    rawThetaRRate,
                    cfg.thetaRMin,
                    cfg.thetaRMax
                );

                etaRate = servoError + cfg.kaw * delta;
                if (eta >= cfg.etaLimit && etaRate > 0.0) {
                    etaRate = 0.0;
                } else if (eta <= -cfg.etaLimit && etaRate < 0.0) {
                    etaRate = 0.0;
                }

                xAwRate = -cfg.lambdaAw * xAw + cfg.b * delta;
            }

            rows.add(new Row(
                t,
                r,
                x,
                xm,
                e,
                epsilon,
                v,
                u,
                delta,
                eta,
                thetaX,
                thetaR,
                gate
            ));

            double xRate = -cfg.a * x + cfg.b * u;
            double xmRate = -cfg.am * xm + cfg.bm * r;

            x += cfg.dt * xRate;
            xm += cfg.dt * xmRate;
            eta += cfg.dt * etaRate;
            xAw += cfg.dt * xAwRate;
            thetaX += cfg.dt * thetaXRate;
            thetaR += cfg.dt * thetaRRate;

            if (
                !Double.isFinite(x)
                || !Double.isFinite(thetaX)
                || !Double.isFinite(thetaR)
                || !Double.isFinite(eta)
            ) {
                throw new IllegalStateException(
                    "Non-finite simulation state detected."
                );
            }
        }

        return rows;
    }

    private static Metrics computeMetrics(
        List<Row> rows,
        Config cfg
    ) {
        Metrics metrics = new Metrics();
        int saturatedSamples = 0;

        for (Row row : rows) {
            metrics.iae += Math.abs(row.trackingError) * cfg.dt;
            metrics.maxAbsCommand = Math.max(
                metrics.maxAbsCommand,
                Math.abs(row.command)
            );
            metrics.maxAbsIntegrator = Math.max(
                metrics.maxAbsIntegrator,
                Math.abs(row.integrator)
            );
            metrics.maxAbsParameter = Math.max(
                metrics.maxAbsParameter,
                Math.max(Math.abs(row.thetaX), Math.abs(row.thetaR))
            );
            if (Math.abs(row.mismatch) > 1.0e-9) {
                saturatedSamples++;
            }
        }

        metrics.saturationFraction =
            (double) saturatedSamples / rows.size();
        metrics.finalAbsTrackingError =
            Math.abs(rows.get(rows.size() - 1).trackingError);
        return metrics;
    }

    private static void writeCsv(
        Path path,
        List<Row> rows
    ) throws IOException {
        try (
            BufferedWriter writer = Files.newBufferedWriter(
                path,
                StandardCharsets.UTF_8
            )
        ) {
            writer.write(
                "time,reference,x,xm,tracking_error,modified_error,"
                + "command,input,mismatch,integrator,theta_x,theta_r,gate"
            );
            writer.newLine();

            for (Row row : rows) {
                writer.write(String.format(
                    Locale.ROOT,
                    "%.10g,%.10g,%.10g,%.10g,%.10g,%.10g,"
                    + "%.10g,%.10g,%.10g,%.10g,%.10g,%.10g,%.10g%n",
                    row.time,
                    row.reference,
                    row.x,
                    row.xm,
                    row.trackingError,
                    row.modifiedError,
                    row.command,
                    row.input,
                    row.mismatch,
                    row.integrator,
                    row.thetaX,
                    row.thetaR,
                    row.gate
                ));
            }
        }
    }

    private static void printMetrics(
        String title,
        Metrics metrics
    ) {
        System.out.println(title);
        System.out.printf(
            Locale.ROOT,
            "  IAE                         : %.6f%n",
            metrics.iae
        );
        System.out.printf(
            Locale.ROOT,
            "  max_abs_command             : %.6f%n",
            metrics.maxAbsCommand
        );
        System.out.printf(
            Locale.ROOT,
            "  max_abs_integrator          : %.6f%n",
            metrics.maxAbsIntegrator
        );
        System.out.printf(
            Locale.ROOT,
            "  max_abs_parameter           : %.6f%n",
            metrics.maxAbsParameter
        );
        System.out.printf(
            Locale.ROOT,
            "  saturation_fraction         : %.6f%n",
            metrics.saturationFraction
        );
        System.out.printf(
            Locale.ROOT,
            "  final_abs_tracking_error    : %.6f%n%n",
            metrics.finalAbsTrackingError
        );
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Config cfg = new Config();

        List<Row> baseline = simulate(false, cfg);
        List<Row> aware = simulate(true, cfg);

        try {
            writeCsv(
                Path.of("Chapter19_Lesson3_naive_java.csv"),
                baseline
            );
            writeCsv(
                Path.of("Chapter19_Lesson3_anti_windup_java.csv"),
                aware
            );
        } catch (IOException error) {
            System.err.println("CSV error: " + error.getMessage());
            System.exit(1);
        }

        printMetrics(
            "Naive adaptive controller",
            computeMetrics(baseline, cfg)
        );
        printMetrics(
            "Saturation-aware adaptive controller",
            computeMetrics(aware, cfg)
        );
        System.out.println("Generated two CSV files.");
    }
}

The MATLAB script reproduces the same simulation, CSV files, metrics, and plots. In Simulink, map each differential equation to an Integrator block, use a Saturation block for \(u\), a Sum block for \(\Delta_u=u-v\), and separate enabled or projected subsystems for the integral and parameter updates.

Chapter19_Lesson3.m

% Chapter19_Lesson3.m
% Adaptive anti-windup demonstration for a saturated first-order plant.
%
% This script compares a naive adaptive controller with a saturation-aware
% controller using:
%   - integral back-calculation,
%   - a mismatch model,
%   - a modified tracking error,
%   - adaptation gating,
%   - sigma leakage,
%   - parameter projection.
%
% Run:
%   Chapter19_Lesson3

clear; clc; close all;

cfg.dt = 0.002;
cfg.horizon = 20.0;
cfg.a = 1.2;
cfg.b = 0.8;
cfg.am = 1.5;
cfg.bm = 1.5;
cfg.umax = 1.2;
cfg.gamma = [1.8; 1.0];
cfg.ki = 0.35;
cfg.kaw = 4.0;
cfg.lambda_aw = 1.5;
cfg.sigma = 0.01;
cfg.gate_width = 0.15;
cfg.eta_limit = 4.0;
cfg.theta_lower = [-5.0; -1.0];
cfg.theta_upper = [2.0; 5.0];

naive = simulate_case(false, cfg);
aware = simulate_case(true, cfg);

writetable(struct2table(naive), ...
    'Chapter19_Lesson3_naive_matlab.csv');
writetable(struct2table(aware), ...
    'Chapter19_Lesson3_anti_windup_matlab.csv');

disp('Naive adaptive controller');
print_metrics(compute_metrics(naive, cfg));

disp('Saturation-aware adaptive controller');
print_metrics(compute_metrics(aware, cfg));

figure('Name', 'Chapter 19 Lesson 3: Adaptive Anti-Windup', ...
       'Color', 'w', 'Position', [100, 100, 1000, 850]);

subplot(4, 1, 1);
plot(naive.time, naive.reference, '--', 'LineWidth', 1.2);
hold on;
plot(naive.time, naive.x, 'LineWidth', 1.2);
plot(aware.time, aware.x, 'LineWidth', 1.2);
grid on;
ylabel('state');
legend('reference', 'naive x', 'anti-windup x', ...
       'Location', 'best');

subplot(4, 1, 2);
plot(naive.time, naive.command, 'LineWidth', 1.1);
hold on;
plot(aware.time, aware.command, 'LineWidth', 1.1);
plot(aware.time, aware.input, '--', 'LineWidth', 1.1);
grid on;
ylabel('control');
legend('naive command v', 'AW command v', 'applied input u', ...
       'Location', 'best');

subplot(4, 1, 3);
plot(naive.time, naive.theta_x, 'LineWidth', 1.1);
hold on;
plot(naive.time, naive.theta_r, 'LineWidth', 1.1);
plot(aware.time, aware.theta_x, '--', 'LineWidth', 1.1);
plot(aware.time, aware.theta_r, '--', 'LineWidth', 1.1);
grid on;
ylabel('parameters');
legend('naive theta_x', 'naive theta_r', ...
       'AW theta_x', 'AW theta_r', 'Location', 'best');

subplot(4, 1, 4);
plot(naive.time, naive.integrator, 'LineWidth', 1.1);
hold on;
plot(aware.time, aware.integrator, 'LineWidth', 1.1);
plot(aware.time, aware.gate, '--', 'LineWidth', 1.1);
grid on;
xlabel('time [s]');
ylabel('internal states');
legend('naive eta', 'AW eta', 'adaptation gate', ...
       'Location', 'best');

exportgraphics(gcf, 'Chapter19_Lesson3_comparison_matlab.png', ...
               'Resolution', 180);

disp('Generated CSV files and Chapter19_Lesson3_comparison_matlab.png');


function data = simulate_case(saturation_aware, cfg)
    steps = round(cfg.horizon / cfg.dt) + 1;
    time = (0:(steps - 1))' * cfg.dt;

    data.time = time;
    data.reference = zeros(steps, 1);
    data.x = zeros(steps, 1);
    data.xm = zeros(steps, 1);
    data.tracking_error = zeros(steps, 1);
    data.modified_error = zeros(steps, 1);
    data.command = zeros(steps, 1);
    data.input = zeros(steps, 1);
    data.mismatch = zeros(steps, 1);
    data.integrator = zeros(steps, 1);
    data.theta_x = zeros(steps, 1);
    data.theta_r = zeros(steps, 1);
    data.gate = ones(steps, 1);

    x = 0.0;
    xm = 0.0;
    eta = 0.0;
    x_aw = 0.0;
    theta = [0.0; 0.5];

    for k = 1:steps
        t = time(k);
        r = reference_signal(t);
        e = x - xm;
        servo_error = r - x;
        phi = [x; r];

        v = theta' * phi + cfg.ki * eta;
        u = saturate(v, cfg.umax);
        delta = u - v;

        if saturation_aware
            epsilon = e - x_aw;
            gate = 1.0 / ...
                (1.0 + (abs(delta) / cfg.gate_width)^2);

            raw_theta_rate = ...
                -cfg.gamma .* phi * epsilon * gate ...
                -cfg.sigma * theta;
            theta_rate = projected_rate( ...
                theta, ...
                raw_theta_rate, ...
                cfg.theta_lower, ...
                cfg.theta_upper);

            eta_rate = servo_error + cfg.kaw * delta;
            if eta >= cfg.eta_limit && eta_rate > 0.0
                eta_rate = 0.0;
            elseif eta <= -cfg.eta_limit && eta_rate < 0.0
                eta_rate = 0.0;
            end

            x_aw_rate = -cfg.lambda_aw * x_aw + cfg.b * delta;
        else
            epsilon = e;
            gate = 1.0;
            theta_rate = -cfg.gamma .* phi * e;
            eta_rate = servo_error;
            x_aw_rate = 0.0;
        end

        x_rate = -cfg.a * x + cfg.b * u;
        xm_rate = -cfg.am * xm + cfg.bm * r;

        data.reference(k) = r;
        data.x(k) = x;
        data.xm(k) = xm;
        data.tracking_error(k) = e;
        data.modified_error(k) = epsilon;
        data.command(k) = v;
        data.input(k) = u;
        data.mismatch(k) = delta;
        data.integrator(k) = eta;
        data.theta_x(k) = theta(1);
        data.theta_r(k) = theta(2);
        data.gate(k) = gate;

        x = x + cfg.dt * x_rate;
        xm = xm + cfg.dt * xm_rate;
        eta = eta + cfg.dt * eta_rate;
        x_aw = x_aw + cfg.dt * x_aw_rate;
        theta = theta + cfg.dt * theta_rate;

        if any(~isfinite([x, xm, eta, x_aw, theta']))
            error('Non-finite simulation state detected.');
        end
    end
end


function r = reference_signal(t)
    if t < 0.5
        r = 0.0;
    elseif t < 6.0
        r = 2.0;
    elseif t < 12.0
        r = 0.5;
    elseif t < 16.0
        r = -1.8;
    else
        r = 0.3;
    end
end


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


function rate = projected_rate(theta, raw_rate, lower, upper)
    rate = raw_rate;
    for i = 1:length(theta)
        if theta(i) <= lower(i) && rate(i) < 0.0
            rate(i) = 0.0;
        elseif theta(i) >= upper(i) && rate(i) > 0.0
            rate(i) = 0.0;
        end
    end
end


function metrics = compute_metrics(data, cfg)
    metrics.IAE = sum(abs(data.tracking_error)) * cfg.dt;
    metrics.max_abs_command = max(abs(data.command));
    metrics.max_abs_integrator = max(abs(data.integrator));
    metrics.max_abs_parameter = max(abs( ...
        [data.theta_x; data.theta_r]));
    metrics.saturation_fraction = mean(abs(data.mismatch) > 1.0e-9);
    metrics.final_abs_tracking_error = ...
        abs(data.tracking_error(end));
end


function print_metrics(metrics)
    fprintf('  IAE                         : %.6f\n', ...
        metrics.IAE);
    fprintf('  max_abs_command             : %.6f\n', ...
        metrics.max_abs_command);
    fprintf('  max_abs_integrator          : %.6f\n', ...
        metrics.max_abs_integrator);
    fprintf('  max_abs_parameter           : %.6f\n', ...
        metrics.max_abs_parameter);
    fprintf('  saturation_fraction         : %.6f\n', ...
        metrics.saturation_fraction);
    fprintf('  final_abs_tracking_error    : %.6f\n\n', ...
        metrics.final_abs_tracking_error);
end

Simulink Block Mapping

The notebook source uses built-in Wolfram Language functions only. It exports CSV data and a multi-panel comparison graphic.

Chapter19_Lesson3.nb

Notebook[{
  Cell["Chapter 19, Lesson 3: Anti-Windup Ideas in the Presence of Adaptation", "Title"],
  Cell["Evaluate the input cell to compare naive and saturation-aware adaptive control.", "Text"],
  Cell["ClearAll[\"Global`*\"];

cfg = <|
   \"dt\" -> 0.002,
   \"horizon\" -> 20.0,
   \"a\" -> 1.2,
   \"b\" -> 0.8,
   \"am\" -> 1.5,
   \"bm\" -> 1.5,
   \"umax\" -> 1.2,
   \"gamma\" -> {1.8, 1.0},
   \"ki\" -> 0.35,
   \"kaw\" -> 4.0,
   \"lambdaAw\" -> 1.5,
   \"sigma\" -> 0.01,
   \"gateWidth\" -> 0.15,
   \"etaLimit\" -> 4.0,
   \"thetaLower\" -> {-5.0, -1.0},
   \"thetaUpper\" -> {2.0, 5.0}
|>;

headers = {
   \"time\", \"reference\", \"x\", \"xm\", \"tracking_error\",
   \"modified_error\", \"command\", \"input\", \"mismatch\",
   \"integrator\", \"theta_x\", \"theta_r\", \"gate\"
};

reference[t_] := Piecewise[{
    {0.0, t < 0.5},
    {2.0, t < 6.0},
    {0.5, t < 12.0},
    {-1.8, t < 16.0}
    }, 0.3];

saturate[value_] := Clip[value, {-cfg[\"umax\"], cfg[\"umax\"]}];

projectedRate[theta_, rawRate_, lower_, upper_] :=
  MapThread[
   Which[
     #1 <= #3 && #2 < 0.0, 0.0,
     #1 >= #4 && #2 > 0.0, 0.0,
     True, #2
     ] &,
   {theta, rawRate, lower, upper}
   ];

simulate[saturationAware_] := Module[
  {
   steps, rows, x = 0.0, xm = 0.0, eta = 0.0,
   xAw = 0.0, theta = {0.0, 0.5}, t, r, e,
   servoError, phi, v, u, delta, epsilon, gate,
   thetaRate, rawThetaRate, etaRate, xAwRate,
   xRate, xmRate, k
   },

  steps = Round[cfg[\"horizon\"]/cfg[\"dt\"]] + 1;
  rows = ConstantArray[0.0, {steps, Length[headers]}];

  Do[
   t = (k - 1) cfg[\"dt\"];
   r = reference[t];
   e = x - xm;
   servoError = r - x;
   phi = {x, r};

   v = theta.phi + cfg[\"ki\"] eta;
   u = saturate[v];
   delta = u - v;

   If[TrueQ[saturationAware],
    epsilon = e - xAw;
    gate = 1.0/(1.0 + (Abs[delta]/cfg[\"gateWidth\"])^2);

    rawThetaRate =
      -cfg[\"gamma\"] phi epsilon gate - cfg[\"sigma\"] theta;
    thetaRate = projectedRate[
      theta,
      rawThetaRate,
      cfg[\"thetaLower\"],
      cfg[\"thetaUpper\"]
      ];

    etaRate = servoError + cfg[\"kaw\"] delta;
    If[eta >= cfg[\"etaLimit\"] && etaRate > 0.0, etaRate = 0.0];
    If[eta <= -cfg[\"etaLimit\"] && etaRate < 0.0, etaRate = 0.0];

    xAwRate = -cfg[\"lambdaAw\"] xAw + cfg[\"b\"] delta;,
    epsilon = e;
    gate = 1.0;
    thetaRate = -cfg[\"gamma\"] phi e;
    etaRate = servoError;
    xAwRate = 0.0;
    ];

   xRate = -cfg[\"a\"] x + cfg[\"b\"] u;
   xmRate = -cfg[\"am\"] xm + cfg[\"bm\"] r;

   rows[[k]] = {
     t, r, x, xm, e, epsilon, v, u, delta,
     eta, theta[[1]], theta[[2]], gate
     };

   x = x + cfg[\"dt\"] xRate;
   xm = xm + cfg[\"dt\"] xmRate;
   eta = eta + cfg[\"dt\"] etaRate;
   xAw = xAw + cfg[\"dt\"] xAwRate;
   theta = theta + cfg[\"dt\"] thetaRate;
   ,
   {k, steps}
   ];

  AssociationThread[headers, Transpose[rows]]
  ];

metrics[data_] := <|
   \"IAE\" ->
    Total[Abs[data[\"tracking_error\"]]] cfg[\"dt\"],
   \"max_abs_command\" ->
    Max[Abs[data[\"command\"]]],
   \"max_abs_integrator\" ->
    Max[Abs[data[\"integrator\"]]],
   \"max_abs_parameter\" ->
    Max[Abs@Join[data[\"theta_x\"], data[\"theta_r\"]]],
   \"saturation_fraction\" ->
    Mean[(Boole[Abs[#] > 10^-9] &) /@ data[\"mismatch\"]],
   \"final_abs_tracking_error\" ->
    Abs[Last[data[\"tracking_error\"]]]
   |>;

exportCsv[path_, data_] :=
 Export[
  path,
  Prepend[Transpose[Lookup[data, headers]], headers],
  \"CSV\"
  ];

baseDirectory = Quiet@Check[NotebookDirectory[], Directory[]];
naive = simulate[False];
aware = simulate[True];

exportCsv[
 FileNameJoin[{baseDirectory, \"Chapter19_Lesson3_naive_mathematica.csv\"}],
 naive
 ];
exportCsv[
 FileNameJoin[{
   baseDirectory,
   \"Chapter19_Lesson3_anti_windup_mathematica.csv\"
   }],
 aware
 ];

Print[\"Naive adaptive controller\"];
Print[Dataset[metrics[naive]]];
Print[\"Saturation-aware adaptive controller\"];
Print[Dataset[metrics[aware]]];

statePlot = ListLinePlot[
   {
    Transpose[{naive[\"time\"], naive[\"reference\"]}],
    Transpose[{naive[\"time\"], naive[\"x\"]}],
    Transpose[{aware[\"time\"], aware[\"x\"]}]
    },
   PlotLegends -> {\"reference\", \"naive x\", \"anti-windup x\"},
   Frame -> True,
   FrameLabel -> {\"time [s]\", \"state\"},
   PlotRange -> All,
   ImageSize -> Large
   ];

controlPlot = ListLinePlot[
   {
    Transpose[{naive[\"time\"], naive[\"command\"]}],
    Transpose[{aware[\"time\"], aware[\"command\"]}],
    Transpose[{aware[\"time\"], aware[\"input\"]}]
    },
   PlotLegends -> {
     \"naive command v\", \"AW command v\", \"applied input u\"
     },
   Frame -> True,
   FrameLabel -> {\"time [s]\", \"control\"},
   PlotRange -> All,
   ImageSize -> Large
   ];

parameterPlot = ListLinePlot[
   {
    Transpose[{naive[\"time\"], naive[\"theta_x\"]}],
    Transpose[{naive[\"time\"], naive[\"theta_r\"]}],
    Transpose[{aware[\"time\"], aware[\"theta_x\"]}],
    Transpose[{aware[\"time\"], aware[\"theta_r\"]}]
    },
   PlotLegends -> {
     \"naive theta_x\", \"naive theta_r\",
     \"AW theta_x\", \"AW theta_r\"
     },
   Frame -> True,
   FrameLabel -> {\"time [s]\", \"adaptive parameters\"},
   PlotRange -> All,
   ImageSize -> Large
   ];

internalPlot = ListLinePlot[
   {
    Transpose[{naive[\"time\"], naive[\"integrator\"]}],
    Transpose[{aware[\"time\"], aware[\"integrator\"]}],
    Transpose[{aware[\"time\"], aware[\"gate\"]}]
    },
   PlotLegends -> {\"naive eta\", \"AW eta\", \"adaptation gate\"},
   Frame -> True,
   FrameLabel -> {\"time [s]\", \"internal states\"},
   PlotRange -> All,
   ImageSize -> Large
   ];

comparison = GraphicsGrid[
   { {statePlot}, {controlPlot}, {parameterPlot}, {internalPlot} },
   Spacings -> {0.5, 0.8}
   ];

Export[
 FileNameJoin[{
   baseDirectory,
   \"Chapter19_Lesson3_comparison_mathematica.png\"
   }],
 comparison,
 ImageResolution -> 180
 ];

comparison", "Input"]
},
WindowTitle -> "Chapter19_Lesson3",
Saveable -> True
]

The naive controller often obtains a smaller accumulated tracking error during an infeasible command by driving its requested input far beyond the actuator limit. That apparent improvement is misleading: the plant receives the same saturated input while the hidden command, integral state, and parameter estimates continue to grow.

With the supplied parameters, the verified Python, C++, and Java runs produce approximately:

The central result is not that saturation disappears. Rather, the saturation-aware loop stores substantially less fictitious control authority and returns to the feasible operating region with smaller internal distortion.

Problem 1 (Failure of the Nominal Lyapunov Proof): Starting from \(\dot e=A_me+b\tilde\theta^T\phi+b\Delta_u\) and the ordinary gradient law, derive \(\dot V\) and explain why saturation can cause parameter drift.

Solution: With \(V=e^TPe+\tilde\theta^T\Gamma^{-1}\tilde\theta\),

\[ \dot V=-e^TQe+2e^TPb\Delta_u. \]

The adaptive-law cross term cancels the parametric term, but no update term cancels the actuator discrepancy. Because \(e^TPb\Delta_u\) may be positive, the Lyapunov function need not decrease. A persistent saturation-generated error also remains inside \(\dot{\hat\theta}\), so the estimate may move even when the plant parameters are constant.

Problem 2 (Modified-Error Cancellation): Let \(\dot x_{aw}=A_mx_{aw}+b\Delta_u\) and \(\varepsilon=e-x_{aw}\). Derive the modified-error dynamics and prove boundedness of the ideal adaptive core.

Solution:

\[ \dot\varepsilon =A_m(e-x_{aw})+b\tilde\theta^T\phi =A_m\varepsilon+b\tilde\theta^T\phi. \]

Using \(\dot{\hat\theta}=-\Gamma\phi\varepsilon^TPb\) in \(V=\varepsilon^TP\varepsilon+ \tilde\theta^T\Gamma^{-1}\tilde\theta\) gives \(\dot V=-\varepsilon^TQ\varepsilon\le0\). Therefore \(\varepsilon\) and \(\tilde\theta\) are bounded. Under bounded-regressor assumptions, \(\varepsilon\to0\). Physical tracking additionally requires \(x_{aw}\to0\).

Problem 3 (Reference Feasibility): For \(\dot x=-1.2x+0.8u\) with \(|u|\le1.2\), determine the largest possible steady-state magnitude. Can a constant reference of 2 be tracked?

Solution: At equilibrium, \(0=-1.2x_{ss}+0.8u_{ss}\), hence

\[ |x_{ss}|=\frac{0.8}{1.2}|u_{ss}| \le\frac{0.8}{1.2}(1.2)=0.8. \]

A constant reference of 2 is infeasible. No anti-windup law can remove this physical limitation; it can only prevent internal states from winding up while the limitation is active.

Problem 4 (Back-Calculation Sign): The command is \(v=v_a+k_I\eta\) with \(k_I>0\). At the upper limit, \(v>\bar u\). Show why \(\dot\eta=e_c+k_{aw}(u-v)\) has the correct unwinding sign.

Solution: Upper saturation gives \(u-v<0\). Therefore \(k_{aw}(u-v)<0\), which decreases \(\eta\) and hence decreases \(k_I\eta\). The command is driven back toward the realizable interval.

Problem 5 (Projection Property): Suppose the true ideal parameter lies in a convex set \(\Omega\). Explain why a tangent-cone projection is preferable to simply clipping the parameter after numerical integration.

Solution: Projection modifies the derivative before integration and can be constructed to satisfy the Lyapunov-compatible inequality

\[ \tilde\theta^T\Gamma^{-1} \left[ \operatorname{Proj}_{\Omega}(\hat\theta,y)-y \right]\le0. \]

Post-integration clipping limits the stored value but does not generally preserve this derivative inequality, especially with finite sampling and large update steps.

Problem 6 (What Should Be Frozen?): Compare full adaptation freeze, smooth gating, and modified-error adaptation during saturation.

Solution: Full freeze is simple and prevents drift but also stops learning from informative transients. Smooth gating retains a reduced update and avoids discontinuous switching, but changes the exact Lyapunov cancellation. Modified-error adaptation is structurally stronger because it subtracts a modeled saturation effect rather than merely suppressing learning. In practice, a mismatch model is often combined with projection, small leakage, and a gate for model imperfections.

Actuator saturation introduces a known discrepancy between requested and applied control. In an adaptive loop, that discrepancy can wind up both classical controller states and parameter estimates. Back-calculation protects integral states, while a mismatch model produces a modified error whose dynamics recover the nominal adaptive structure. The ideal Lyapunov proof establishes boundedness and convergence of the modified error, not tracking of a permanently infeasible reference. Projection, leakage, and saturation-dependent gating provide practical safeguards but must be analyzed as robust modifications rather than assumed to preserve the exact nominal theorem.