Chapter 24: Implementation Issues and Practical Pitfalls

Lesson 5: Best Practices: Safety Checks, Saturation Handling, and Documentation

This lesson consolidates implementation practices required to move an adaptive controller from a mathematical design to a defensible real-time system. The central principle is that adaptation is a performance mechanism, not the final safety authority. We develop a layered architecture that validates signals, enforces state and actuator constraints, suppresses harmful learning during saturation, bounds parameter motion, provides deterministic fallback behavior, and records enough evidence to reproduce every closed-loop decision.

1. From Stable Equations to a Safe Implementation

A Lyapunov proof is conditional: it applies to the mathematical interconnection for which the proof was constructed. A deployed controller contains sensor conversions, sample-and-hold devices, finite-word arithmetic, schedulers, saturation nonlinearities, rate limits, communication delays, and reset logic. Each item can alter the error dynamics. Therefore, the implementation must explicitly separate four responsibilities:

  1. Nominal adaptive control: computes the unconstrained command that pursues tracking.
  2. Safety supervision: modifies or rejects commands that threaten a certified state envelope.
  3. Actuator management: enforces magnitude, rate, thermal, and duty-cycle limits.
  4. Assurance and evidence: validates inputs, monitors invariants, switches to fallback control, and logs decisions.
flowchart TD
  R["Reference and mode request"] --> V["Validate, timestamp, \nand limit"]
  Y["Sensors"] --> V
  V --> A["Nominal adaptive controller"]
  A --> S["Safety filter"]
  S --> L["Magnitude and \nslew-rate limits"]
  L --> P["Plant and actuator"]
  P --> Y
  L --> M["Constraint and \nmismatch monitor"]
  M --> G["Adaptation gate, \nleakage, projection"]
  G --> A
  V --> F["Watchdog and \nfallback controller"]
  F --> S
  A --> T["Telemetry and \nconfiguration record"]
  S --> T
  L --> T
  M --> T
        

The safety path should be simpler than the adaptive path, execute at least as reliably, and remain available when parameter estimates, optimization routines, or nonessential communications fail.

2. Commanded Input, Applied Input, and the Saturation Mismatch

Consider a linearly parameterized controller with regressor \( \boldsymbol{\phi}(t) \) and estimate \( \hat{\boldsymbol{\theta} }(t) \):

\[ u_c = \hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}, \qquad u_a = \operatorname{sat}_{[u_{\min},u_{\max}]}(u_c), \qquad \Delta u = u_a-u_c. \]

The plant receives \(u_a\), not \(u_c\). If the nominal matching derivation gives

\[ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} + \mathbf{B}\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}, \]

then the implemented error model contains an additional forcing term:

\[ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} + \mathbf{B}\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi} + \mathbf{B}\Delta u. \]

With \(\mathbf{A}_m^T\mathbf{P}+\mathbf{P}\mathbf{A}_m=-\mathbf{Q}\) and \(V=\mathbf{e}^T\mathbf{P}\mathbf{e}+\tilde{\boldsymbol{\theta} }^T\boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta} }\), the nominal cancellation no longer removes every cross term:

\[ \dot V = -\mathbf{e}^T\mathbf{Q}\mathbf{e} + 2\mathbf{e}^T\mathbf{P}\mathbf{B}\Delta u. \]

Young's inequality yields, for any \(\eta > 0\),

\[ 2\mathbf{e}^T\mathbf{P}\mathbf{B}\Delta u \leq \eta\lVert\mathbf{e}\rVert^2 + \frac{\lVert\mathbf{P}\mathbf{B}\rVert^2}{\eta}|\Delta u|^2. \]

Thus persistent saturation generally changes asymptotic tracking into an ultimate-boundedness statement. More importantly, a naive estimator can interpret the saturation-induced tracking error as parameter error and continue increasing controller gains, which enlarges \(|u_c|\) while the actuator remains fixed at its limit.

3. Saturation-Aware Adaptation, Leakage, Rate Limits, and Projection

A practical update law combines an intervention gate \(g_s\in[0,1]\), robust leakage, a parameter-rate limiter, and projection onto an admissible set \(\Omega_\theta\):

\[ \dot{\hat{\boldsymbol{\theta} } } = \operatorname{Proj}_{\Omega_\theta}\!\left( \hat{\boldsymbol{\theta} }, \operatorname{sat}_{\dot\theta_{\max} }\!\left[ -g_s\boldsymbol{\Gamma}\boldsymbol{\phi}\mathbf{e}^T\mathbf{P}\mathbf{B} -\sigma\boldsymbol{\Gamma}\hat{\boldsymbol{\theta} } \right]\right). \]

A conservative gate is

\[ g_s = \begin{cases} 1, & |u_a-u_c|\leq\varepsilon_u \text{ and no safety intervention is active},\\ 0, & \text{otherwise}. \end{cases} \]

Freezing the tracking-error term during intervention prevents the estimator from learning actuator deficiency as though it were plant uncertainty. The leakage term remains active and slowly returns unsupported estimates toward a chosen center. Projection enforces hard parameter bounds; the rate limiter prevents a single corrupted sample from producing a large parameter jump.

For a box

\[ \Omega_{\theta} = \left\{ \boldsymbol{\theta} \;:\; \theta_i^{-} \leq \theta_i \leq \theta_i^{+} \right\} \]

, a discrete implementation may use

\[ \hat{\theta}_i[k+1] = \min\!\left( \theta_i^{+}, \max\!\left( \theta_i^{-}, \hat{\theta}_i[k] + T_s\dot{\hat{\theta} }_i[k] \right) \right) \]

Bounds must represent physically plausible controller parameters and preserve the known control direction. Projection is not a substitute for correct modeling; an incorrectly chosen set can exclude every stabilizing parameter.

4. Safety Filters and Forward-Invariant Operating Envelopes

Let a scalar plant satisfy \(\dot x=f(x)+b(x)u\), with a desired safe set \(\mathcal{C}=\{x:h(x)\geq0\}\). For the symmetric state constraint \(|x|\leq x_{\max}\), choose

\[ h(x)=x_{\max}^2-x^2. \]

A zeroing barrier condition is

\[ \dot h(x,u)+\kappa h(x)\geq0, \qquad \kappa > 0. \]

Since \(\dot h=-2x[f(x)+b(x)u]\), the condition is affine in \(u\):

\[ -2x b(x)u \geq 2x f(x)-\kappa h(x). \]

A safety filter selects the command nearest the adaptive request while satisfying the barrier inequality and actuator bounds:

\[ u_{\mathrm{safe} } = \arg\min_u \frac12(u-u_c)^2 \quad \text{subject to} \quad \dot h(x,u)+\kappa h(x)\geq0,\; u_{\min}\leq u\leq u_{\max}. \]

For uncertain drift, replace \(f(x)\) with conservative upper or lower bounds appropriate to the sign of \(x\). If the constraints become infeasible, the supervisor must report infeasibility and apply a separately analyzed emergency action; silently returning the nominal command defeats the purpose of the filter.

5. Signal Validation, Watchdogs, and Deterministic Fallback

Every adaptive update should be preceded by explicit checks. At minimum:

  • Reject NaN, infinity, impossible ranges, invalid status bits, and stale timestamps.
  • Check sample interval bounds before numerical integration.
  • Limit references and reference derivatives to the tested operating envelope.
  • Cross-check redundant sensors or use analytical consistency residuals when available.
  • Disable adaptation during initialization, mode transitions, estimator resets, and actuator faults.

Let \(q[k]\in\{0,1\}\) denote the aggregate validity flag. A bumpless supervisory blend can be written

\[ u_c[k]=\lambda[k]u_{\mathrm{ad} }[k]+(1-\lambda[k])u_{\mathrm{fb} }[k], \qquad 0\leq\lambda[k]\leq1, \]

\[ \lambda[k+1]=\operatorname{sat}_{[0,1]}\!\left( \lambda[k]+\rho_{\uparrow}q[k]-\rho_{\downarrow}(1-q[k])\right). \]

Safety-critical implementations commonly use immediate transfer to fallback and a slower, qualified return to adaptive operation. Re-entry should require valid sensors, bounded errors, actuator headroom, acceptable parameter estimates, and a dwell time without new faults.

6. Numerical and Real-Time Guards

Continuous-time update laws are implemented numerically. For forward Euler,

\[ \hat{\boldsymbol{\theta} }[k+1]=\hat{\boldsymbol{\theta} }[k] +T_s\mathbf{f}_\theta[k]. \]

The product \(T_s\boldsymbol{\Gamma}\) controls the actual discrete adaptation step. A gain safe in a fast simulation may be excessive on slower hardware. The implementation should therefore monitor

\[ \delta_\theta[k]=\lVert\hat{\boldsymbol{\theta} }[k+1]- \hat{\boldsymbol{\theta} }[k]\rVert, \qquad J_{\mathrm{cpu} }[k]=\frac{t_{\mathrm{finish} }[k]-t_{\mathrm{start} }[k]}{T_s}. \]

Trigger a diagnostic or fallback when the parameter step, missed-deadline count, solver iteration count, covariance condition number, or timestamp jitter exceeds its qualified threshold. State variables and estimates should be checked for finite values after every integration step, not only at sensor ingress.

7. Documentation as Part of the Control Design

Documentation must make the implemented closed loop reconstructable. A useful controller data package records:

  • plant model assumptions, uncertainty structure, control direction, and excluded dynamics;
  • reference model, regressors, units, signs, coordinate frames, and sampling rates;
  • all gains, parameter bounds, saturation limits, rate limits, safety margins, and reset values;
  • mode-transition logic, fallback controller, watchdog thresholds, and re-entry conditions;
  • software version, compiler/runtime version, configuration hash, calibration revision, and test seed;
  • telemetry dictionary with units, rates, validity semantics, and retention policy;
  • verification evidence linking each requirement to analysis, simulation, hardware-in-the-loop, and test results.

A configuration vector \(\boldsymbol{\rho}\) should be treated as controlled data. Each run should store an immutable identifier such as

\[ H_{\mathrm{run} }=\operatorname{Hash}(\text{software revision}, \boldsymbol{\rho},\text{calibration},\text{test scenario}). \]

Without this information, an apparently successful plot cannot be independently reproduced or compared with a later regression.

8. Verification and Release Workflow

flowchart TD
  A["Freeze requirements and safe envelope"] --> B["Review equations, units, signs, and bounds"]
  B --> C["Nominal and Monte Carlo simulation"]
  C --> D["Fault injection: sensor, actuator, timing, and communication"]
  D --> E["Processor-in-the-loop and numerical stress tests"]
  E --> F["Hardware-in-the-loop with physical limit models"]
  F --> G["Low-energy commissioning with adaptation disabled"]
  G --> H["Enable bounded adaptation in stages"]
  H --> I["Independent review of logs and invariant monitors"]
  I --> J["Baseline configuration and regression tests"]
  J --> K["Operational monitoring and controlled updates"]
        

The order matters. Adaptation should not be the first mechanism used to compensate for sign errors, scaling mistakes, excessive latency, or an unstable fallback controller. Commission the deterministic baseline first; then introduce adaptation with small bounds and gains while retaining an immediate disable path.

9. Python Implementation

The following standard-library implementation simulates an uncertain first-order plant. It writes telemetry, injects a sensor fault and disturbance pulse, and reports every safety or actuator intervention.

Chapter24_Lesson5.py

"""Chapter24_Lesson5.py
Safety-supervised adaptive control demonstration.

The example is intentionally self-contained and uses only the Python standard
library.  It demonstrates:
  * sensor validity and finite-value checks;
  * reference limiting;
  * a scalar control-barrier safety filter;
  * actuator magnitude and slew-rate limits;
  * saturation-aware adaptation gating;
  * sigma leakage, parameter-rate limits, and projection;
  * watchdog fallback and CSV telemetry.
"""

from __future__ import annotations

from dataclasses import dataclass
import csv
import math
from pathlib import Path


@dataclass(frozen=True)
class Config:
    dt: float = 0.005
    duration: float = 20.0
    u_min: float = -2.0
    u_max: float = 2.0
    du_max: float = 8.0
    x_safe: float = 2.4
    reference_limit: float = 1.6
    cbf_rate: float = 3.0
    drift_bound: float = 1.25
    input_gain_min: float = 0.8
    gamma_x: float = 1.5
    gamma_r: float = 1.0
    sigma: float = 0.03
    theta_x_min: float = -4.0
    theta_x_max: float = -0.05
    theta_r_min: float = 0.0
    theta_r_max: float = 4.0
    theta_rate_max: float = 1.0
    fallback_gain: float = 1.8
    intervention_tolerance: float = 1.0e-7


def clamp(value: float, lower: float, upper: float) -> float:
    return max(lower, min(upper, value))


class SafetyAdaptiveController:
    def __init__(self, cfg: Config) -> None:
        self.cfg = cfg
        self.theta_x = -1.0
        self.theta_r = 1.0
        self.xm = 0.0
        self.previous_u = 0.0
        self.safety_interventions = 0
        self.fallback_events = 0
        self.saturation_events = 0

    def cbf_filter(self, x: float, nominal_u: float) -> tuple[float, bool]:
        """Project nominal_u onto a conservative one-dimensional CBF interval."""
        c = self.cfg
        lower, upper = c.u_min, c.u_max
        h = c.x_safe * c.x_safe - x * x
        eps = 1.0e-9

        if x > eps:
            cbf_upper = (c.cbf_rate * h / (2.0 * x) - c.drift_bound) / c.input_gain_min
            upper = min(upper, cbf_upper)
        elif x < -eps:
            q = -x
            cbf_lower = (-c.cbf_rate * h / (2.0 * q) + c.drift_bound) / c.input_gain_min
            lower = max(lower, cbf_lower)

        # If modeling conservatism makes the interval empty, choose the bound that
        # drives the state toward the origin instead of returning an invalid command.
        if lower > upper:
            safe_u = c.u_min if x > 0.0 else c.u_max
            return safe_u, True

        filtered = clamp(nominal_u, lower, upper)
        return filtered, abs(filtered - nominal_u) > c.intervention_tolerance

    def rate_limit(self, target_u: float) -> tuple[float, bool]:
        max_step = self.cfg.du_max * self.cfg.dt
        limited = clamp(target_u, self.previous_u - max_step, self.previous_u + max_step)
        active = abs(limited - target_u) > self.cfg.intervention_tolerance
        self.previous_u = limited
        return limited, active

    def update(
        self,
        x_measured: float,
        reference: float,
        sensor_valid: bool,
    ) -> dict[str, float | int | str]:
        c = self.cfg
        valid = sensor_valid and math.isfinite(x_measured) and math.isfinite(reference)
        r = clamp(reference if math.isfinite(reference) else 0.0,
                  -c.reference_limit, c.reference_limit)

        # Stable reference model: x_m_dot = -a_m x_m + b_m r.
        a_m, b_m = 1.2, 1.2
        self.xm += c.dt * (-a_m * self.xm + b_m * r)

        if not valid:
            self.fallback_events += 1
            nominal_u = clamp(-c.fallback_gain * (x_measured if math.isfinite(x_measured) else 0.0),
                              c.u_min, c.u_max)
            filtered_u = nominal_u
            safety_active = False
            status = "FALLBACK_INVALID_SENSOR"
            adaptation_gate = 0.0
        else:
            nominal_u = self.theta_x * x_measured + self.theta_r * r
            filtered_u, safety_active = self.cbf_filter(x_measured, nominal_u)
            status = "SAFETY_FILTER" if safety_active else "ADAPTIVE"
            adaptation_gate = 0.0 if safety_active else 1.0
            if safety_active:
                self.safety_interventions += 1

        magnitude_limited = clamp(filtered_u, c.u_min, c.u_max)
        magnitude_saturated = abs(magnitude_limited - filtered_u) > c.intervention_tolerance
        applied_u, rate_limited = self.rate_limit(magnitude_limited)
        constrained = safety_active or magnitude_saturated or rate_limited

        if magnitude_saturated or rate_limited:
            self.saturation_events += 1
            adaptation_gate = 0.0
            if status == "ADAPTIVE":
                status = "ACTUATOR_LIMIT"

        error = (x_measured - self.xm) if valid else 0.0

        # Saturation-aware robust update:
        # theta_dot = Proj_rate(-g Gamma phi e - sigma theta).
        adaptation_x = x_measured if valid else 0.0
        raw_theta_x_dot = -adaptation_gate * c.gamma_x * adaptation_x * error - c.sigma * self.theta_x
        raw_theta_r_dot = -adaptation_gate * c.gamma_r * r * error - c.sigma * self.theta_r
        theta_x_dot = clamp(raw_theta_x_dot, -c.theta_rate_max, c.theta_rate_max)
        theta_r_dot = clamp(raw_theta_r_dot, -c.theta_rate_max, c.theta_rate_max)

        self.theta_x = clamp(self.theta_x + c.dt * theta_x_dot,
                             c.theta_x_min, c.theta_x_max)
        self.theta_r = clamp(self.theta_r + c.dt * theta_r_dot,
                             c.theta_r_min, c.theta_r_max)

        return {
            "x_m": self.xm,
            "error": error,
            "reference_limited": r,
            "u_nominal": nominal_u,
            "u_applied": applied_u,
            "theta_x": self.theta_x,
            "theta_r": self.theta_r,
            "adaptation_gate": adaptation_gate,
            "safety_active": int(safety_active),
            "constrained": int(constrained),
            "sensor_valid": int(valid),
            "status": status,
        }


def reference_signal(t: float) -> float:
    if t < 2.0:
        return 0.0
    if t < 8.0:
        return 1.25
    if t < 14.0:
        return -1.1
    return 1.45 * math.sin(0.55 * (t - 14.0))


def main() -> None:
    cfg = Config()
    controller = SafetyAdaptiveController(cfg)

    # Unknown plant used only by the simulation, not by the adaptive controller.
    true_a = -0.45
    true_b = 1.0
    x = 0.0
    peak_abs_x = 0.0
    output = Path("Chapter24_Lesson5_telemetry.csv")

    fields = [
        "time", "x", "x_m", "error", "reference", "reference_limited",
        "u_nominal", "u_applied", "theta_x", "theta_r", "adaptation_gate",
        "safety_active", "constrained", "sensor_valid", "status",
    ]

    with output.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.DictWriter(stream, fieldnames=fields)
        writer.writeheader()

        steps = int(cfg.duration / cfg.dt)
        for k in range(steps + 1):
            t = k * cfg.dt
            r = reference_signal(t)

            # Deliberate sensor fault demonstrates deterministic fallback behavior.
            sensor_valid = not (12.0 <= t < 12.15)
            measured_x = x if sensor_valid else float("nan")
            data = controller.update(measured_x, r, sensor_valid)

            disturbance = 0.12 * math.sin(1.7 * t)
            if 9.0 <= t < 9.4:
                disturbance += 0.55
            x_dot = true_a * x + true_b * float(data["u_applied"]) + disturbance
            x += cfg.dt * x_dot
            peak_abs_x = max(peak_abs_x, abs(x))

            writer.writerow({
                "time": f"{t:.6f}",
                "x": f"{x:.9f}",
                "x_m": f"{float(data['x_m']):.9f}",
                "error": f"{float(data['error']):.9f}",
                "reference": f"{r:.9f}",
                "reference_limited": f"{float(data['reference_limited']):.9f}",
                "u_nominal": f"{float(data['u_nominal']):.9f}",
                "u_applied": f"{float(data['u_applied']):.9f}",
                "theta_x": f"{float(data['theta_x']):.9f}",
                "theta_r": f"{float(data['theta_r']):.9f}",
                "adaptation_gate": f"{float(data['adaptation_gate']):.1f}",
                "safety_active": data["safety_active"],
                "constrained": data["constrained"],
                "sensor_valid": data["sensor_valid"],
                "status": data["status"],
            })

    print(f"Telemetry: {output.resolve()}")
    print(f"Peak |x|: {peak_abs_x:.6f} (configured safe boundary {cfg.x_safe:.3f})")
    print(f"Final theta_x={controller.theta_x:.6f}, theta_r={controller.theta_r:.6f}")
    print(f"Safety interventions: {controller.safety_interventions}")
    print(f"Actuator/rate interventions: {controller.saturation_events}")
    print(f"Fallback samples: {controller.fallback_events}")


if __name__ == "__main__":
    main()

10. C++ Implementation

The C++17 version uses deterministic fixed-step arithmetic, explicit finite-value checks, bounded updates, and a nonzero exit code if the simulated safe boundary is violated.

Chapter24_Lesson5.cpp

// Chapter24_Lesson5.cpp
// Safety-supervised adaptive control demonstration (C++17).

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <tuple>

struct Config {
    double dt = 0.005;
    double duration = 20.0;
    double uMin = -2.0;
    double uMax = 2.0;
    double duMax = 8.0;
    double xSafe = 2.4;
    double referenceLimit = 1.6;
    double cbfRate = 3.0;
    double driftBound = 1.25;
    double inputGainMin = 0.8;
    double gammaX = 1.5;
    double gammaR = 1.0;
    double sigma = 0.03;
    double thetaXMin = -4.0;
    double thetaXMax = -0.05;
    double thetaRMin = 0.0;
    double thetaRMax = 4.0;
    double thetaRateMax = 1.0;
    double fallbackGain = 1.8;
    double tolerance = 1.0e-7;
};

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

struct StepResult {
    double xModel{};
    double error{};
    double referenceLimited{};
    double uNominal{};
    double uApplied{};
    double thetaX{};
    double thetaR{};
    double adaptationGate{};
    bool safetyActive{};
    bool constrained{};
    bool sensorValid{};
    std::string status;
};

class SafetyAdaptiveController {
public:
    explicit SafetyAdaptiveController(Config cfg) : cfg_(cfg) {}

    StepResult update(double xMeasured, double reference, bool sensorValid) {
        const bool valid = sensorValid && std::isfinite(xMeasured) && std::isfinite(reference);
        const double r = clampValue(std::isfinite(reference) ? reference : 0.0,
                                    -cfg_.referenceLimit, cfg_.referenceLimit);

        constexpr double aModel = 1.2;
        constexpr double bModel = 1.2;
        xModel_ += cfg_.dt * (-aModel * xModel_ + bModel * r);

        double nominalU = 0.0;
        double filteredU = 0.0;
        bool safetyActive = false;
        double adaptationGate = 0.0;
        std::string status;

        if (!valid) {
            ++fallbackEvents_;
            const double safeMeasurement = std::isfinite(xMeasured) ? xMeasured : 0.0;
            nominalU = clampValue(-cfg_.fallbackGain * safeMeasurement, cfg_.uMin, cfg_.uMax);
            filteredU = nominalU;
            status = "FALLBACK_INVALID_SENSOR";
        } else {
            nominalU = thetaX_ * xMeasured + thetaR_ * r;
            std::tie(filteredU, safetyActive) = cbfFilter(xMeasured, nominalU);
            adaptationGate = safetyActive ? 0.0 : 1.0;
            status = safetyActive ? "SAFETY_FILTER" : "ADAPTIVE";
            if (safetyActive) {
                ++safetyInterventions_;
            }
        }

        const double magnitudeLimited = clampValue(filteredU, cfg_.uMin, cfg_.uMax);
        const bool magnitudeSaturated = std::abs(magnitudeLimited - filteredU) > cfg_.tolerance;
        const auto [appliedU, rateLimited] = rateLimit(magnitudeLimited);
        const bool constrained = safetyActive || magnitudeSaturated || rateLimited;

        if (magnitudeSaturated || rateLimited) {
            ++saturationEvents_;
            adaptationGate = 0.0;
            if (status == "ADAPTIVE") {
                status = "ACTUATOR_LIMIT";
            }
        }

        const double error = valid ? (xMeasured - xModel_) : 0.0;
        const double adaptationX = valid ? xMeasured : 0.0;
        const double rawThetaXDot = -adaptationGate * cfg_.gammaX * adaptationX * error
                                    - cfg_.sigma * thetaX_;
        const double rawThetaRDot = -adaptationGate * cfg_.gammaR * r * error
                                    - cfg_.sigma * thetaR_;
        const double thetaXDot = clampValue(rawThetaXDot, -cfg_.thetaRateMax, cfg_.thetaRateMax);
        const double thetaRDot = clampValue(rawThetaRDot, -cfg_.thetaRateMax, cfg_.thetaRateMax);

        thetaX_ = clampValue(thetaX_ + cfg_.dt * thetaXDot,
                             cfg_.thetaXMin, cfg_.thetaXMax);
        thetaR_ = clampValue(thetaR_ + cfg_.dt * thetaRDot,
                             cfg_.thetaRMin, cfg_.thetaRMax);

        return {xModel_, error, r, nominalU, appliedU, thetaX_, thetaR_,
                adaptationGate, safetyActive, constrained, valid, status};
    }

    int safetyInterventions() const { return safetyInterventions_; }
    int fallbackEvents() const { return fallbackEvents_; }
    int saturationEvents() const { return saturationEvents_; }
    double thetaX() const { return thetaX_; }
    double thetaR() const { return thetaR_; }

private:
    std::pair<double, bool> cbfFilter(double x, double nominalU) const {
        double lower = cfg_.uMin;
        double upper = cfg_.uMax;
        const double h = cfg_.xSafe * cfg_.xSafe - x * x;
        constexpr double eps = 1.0e-9;

        if (x > eps) {
            const double cbfUpper = (cfg_.cbfRate * h / (2.0 * x) - cfg_.driftBound)
                                    / cfg_.inputGainMin;
            upper = std::min(upper, cbfUpper);
        } else if (x < -eps) {
            const double q = -x;
            const double cbfLower = (-cfg_.cbfRate * h / (2.0 * q) + cfg_.driftBound)
                                    / cfg_.inputGainMin;
            lower = std::max(lower, cbfLower);
        }

        if (lower > upper) {
            return {x > 0.0 ? cfg_.uMin : cfg_.uMax, true};
        }
        const double filtered = clampValue(nominalU, lower, upper);
        return {filtered, std::abs(filtered - nominalU) > cfg_.tolerance};
    }

    std::pair<double, bool> rateLimit(double targetU) {
        const double maxStep = cfg_.duMax * cfg_.dt;
        const double limited = clampValue(targetU, previousU_ - maxStep, previousU_ + maxStep);
        const bool active = std::abs(limited - targetU) > cfg_.tolerance;
        previousU_ = limited;
        return {limited, active};
    }

    Config cfg_;
    double thetaX_ = -1.0;
    double thetaR_ = 1.0;
    double xModel_ = 0.0;
    double previousU_ = 0.0;
    int safetyInterventions_ = 0;
    int fallbackEvents_ = 0;
    int saturationEvents_ = 0;
};

double referenceSignal(double t) {
    if (t < 2.0) return 0.0;
    if (t < 8.0) return 1.25;
    if (t < 14.0) return -1.1;
    return 1.45 * std::sin(0.55 * (t - 14.0));
}

int main() {
    try {
        const Config cfg;
        SafetyAdaptiveController controller(cfg);
        std::ofstream out("Chapter24_Lesson5_telemetry_cpp.csv");
        if (!out) {
            throw std::runtime_error("Cannot open telemetry output file.");
        }
        out << "time,x,x_m,error,reference,reference_limited,u_nominal,u_applied,"
               "theta_x,theta_r,adaptation_gate,safety_active,constrained,sensor_valid,status\n";
        out << std::fixed << std::setprecision(9);

        constexpr double trueA = -0.45;
        constexpr double trueB = 1.0;
        double x = 0.0;
        double peakAbsX = 0.0;
        const int steps = static_cast<int>(cfg.duration / cfg.dt);

        for (int k = 0; k <= steps; ++k) {
            const double t = k * cfg.dt;
            const double r = referenceSignal(t);
            const bool sensorValid = !(t >= 12.0 && t < 12.15);
            const double measuredX = sensorValid ? x : std::numeric_limits<double>::quiet_NaN();
            const StepResult data = controller.update(measuredX, r, sensorValid);

            double disturbance = 0.12 * std::sin(1.7 * t);
            if (t >= 9.0 && t < 9.4) disturbance += 0.55;
            const double xDot = trueA * x + trueB * data.uApplied + disturbance;
            x += cfg.dt * xDot;
            peakAbsX = std::max(peakAbsX, std::abs(x));

            out << t << ',' << x << ',' << data.xModel << ',' << data.error << ','
                << r << ',' << data.referenceLimited << ',' << data.uNominal << ','
                << data.uApplied << ',' << data.thetaX << ',' << data.thetaR << ','
                << data.adaptationGate << ',' << data.safetyActive << ','
                << data.constrained << ',' << data.sensorValid << ',' << data.status << '\n';
        }

        std::cout << std::fixed << std::setprecision(6)
                  << "Peak |x|: " << peakAbsX << " (safe boundary " << cfg.xSafe << ")\n"
                  << "Final theta_x=" << controller.thetaX()
                  << ", theta_r=" << controller.thetaR() << '\n'
                  << "Safety interventions: " << controller.safetyInterventions() << '\n'
                  << "Actuator/rate interventions: " << controller.saturationEvents() << '\n'
                  << "Fallback samples: " << controller.fallbackEvents() << '\n';
        return peakAbsX <= cfg.xSafe + 1.0e-6 ? 0 : 2;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

11. Java Implementation

The Java implementation mirrors the same supervisor and produces a CSV log suitable for regression comparison with the Python and C++ runs.

Chapter24_Lesson5.java

// Chapter24_Lesson5.java
// Safety-supervised adaptive control demonstration (Java 17+).

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

public final class Chapter24_Lesson5 {
    private Chapter24_Lesson5() {}

    private static final class Config {
        final double dt = 0.005;
        final double duration = 20.0;
        final double uMin = -2.0;
        final double uMax = 2.0;
        final double duMax = 8.0;
        final double xSafe = 2.4;
        final double referenceLimit = 1.6;
        final double cbfRate = 3.0;
        final double driftBound = 1.25;
        final double inputGainMin = 0.8;
        final double gammaX = 1.5;
        final double gammaR = 1.0;
        final double sigma = 0.03;
        final double thetaXMin = -4.0;
        final double thetaXMax = -0.05;
        final double thetaRMin = 0.0;
        final double thetaRMax = 4.0;
        final double thetaRateMax = 1.0;
        final double fallbackGain = 1.8;
        final double tolerance = 1.0e-7;
    }

    private record FilterResult(double value, boolean active) {}

    private record StepResult(
            double xModel,
            double error,
            double referenceLimited,
            double uNominal,
            double uApplied,
            double thetaX,
            double thetaR,
            double adaptationGate,
            boolean safetyActive,
            boolean constrained,
            boolean sensorValid,
            String status) {}

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

    private static final class SafetyAdaptiveController {
        private final Config cfg;
        private double thetaX = -1.0;
        private double thetaR = 1.0;
        private double xModel = 0.0;
        private double previousU = 0.0;
        private int safetyInterventions = 0;
        private int fallbackEvents = 0;
        private int saturationEvents = 0;

        SafetyAdaptiveController(Config cfg) {
            this.cfg = cfg;
        }

        StepResult update(double xMeasured, double reference, boolean sensorValid) {
            boolean valid = sensorValid && Double.isFinite(xMeasured) && Double.isFinite(reference);
            double r = clamp(Double.isFinite(reference) ? reference : 0.0,
                    -cfg.referenceLimit, cfg.referenceLimit);

            double aModel = 1.2;
            double bModel = 1.2;
            xModel += cfg.dt * (-aModel * xModel + bModel * r);

            double nominalU;
            double filteredU;
            boolean safetyActive = false;
            double adaptationGate = 0.0;
            String status;

            if (!valid) {
                fallbackEvents++;
                double safeMeasurement = Double.isFinite(xMeasured) ? xMeasured : 0.0;
                nominalU = clamp(-cfg.fallbackGain * safeMeasurement, cfg.uMin, cfg.uMax);
                filteredU = nominalU;
                status = "FALLBACK_INVALID_SENSOR";
            } else {
                nominalU = thetaX * xMeasured + thetaR * r;
                FilterResult filter = cbfFilter(xMeasured, nominalU);
                filteredU = filter.value();
                safetyActive = filter.active();
                adaptationGate = safetyActive ? 0.0 : 1.0;
                status = safetyActive ? "SAFETY_FILTER" : "ADAPTIVE";
                if (safetyActive) safetyInterventions++;
            }

            double magnitudeLimited = clamp(filteredU, cfg.uMin, cfg.uMax);
            boolean magnitudeSaturated = Math.abs(magnitudeLimited - filteredU) > cfg.tolerance;
            FilterResult rate = rateLimit(magnitudeLimited);
            double appliedU = rate.value();
            boolean constrained = safetyActive || magnitudeSaturated || rate.active();

            if (magnitudeSaturated || rate.active()) {
                saturationEvents++;
                adaptationGate = 0.0;
                if (status.equals("ADAPTIVE")) status = "ACTUATOR_LIMIT";
            }

            double error = valid ? xMeasured - xModel : 0.0;
            double adaptationX = valid ? xMeasured : 0.0;
            double rawThetaXDot = -adaptationGate * cfg.gammaX * adaptationX * error
                    - cfg.sigma * thetaX;
            double rawThetaRDot = -adaptationGate * cfg.gammaR * r * error
                    - cfg.sigma * thetaR;
            double thetaXDot = clamp(rawThetaXDot, -cfg.thetaRateMax, cfg.thetaRateMax);
            double thetaRDot = clamp(rawThetaRDot, -cfg.thetaRateMax, cfg.thetaRateMax);

            thetaX = clamp(thetaX + cfg.dt * thetaXDot, cfg.thetaXMin, cfg.thetaXMax);
            thetaR = clamp(thetaR + cfg.dt * thetaRDot, cfg.thetaRMin, cfg.thetaRMax);

            return new StepResult(xModel, error, r, nominalU, appliedU, thetaX, thetaR,
                    adaptationGate, safetyActive, constrained, valid, status);
        }

        private FilterResult cbfFilter(double x, double nominalU) {
            double lower = cfg.uMin;
            double upper = cfg.uMax;
            double h = cfg.xSafe * cfg.xSafe - x * x;
            double eps = 1.0e-9;

            if (x > eps) {
                double cbfUpper = (cfg.cbfRate * h / (2.0 * x) - cfg.driftBound)
                        / cfg.inputGainMin;
                upper = Math.min(upper, cbfUpper);
            } else if (x < -eps) {
                double q = -x;
                double cbfLower = (-cfg.cbfRate * h / (2.0 * q) + cfg.driftBound)
                        / cfg.inputGainMin;
                lower = Math.max(lower, cbfLower);
            }

            if (lower > upper) {
                return new FilterResult(x > 0.0 ? cfg.uMin : cfg.uMax, true);
            }
            double filtered = clamp(nominalU, lower, upper);
            return new FilterResult(filtered, Math.abs(filtered - nominalU) > cfg.tolerance);
        }

        private FilterResult rateLimit(double targetU) {
            double maxStep = cfg.duMax * cfg.dt;
            double limited = clamp(targetU, previousU - maxStep, previousU + maxStep);
            boolean active = Math.abs(limited - targetU) > cfg.tolerance;
            previousU = limited;
            return new FilterResult(limited, active);
        }
    }

    private static double referenceSignal(double t) {
        if (t < 2.0) return 0.0;
        if (t < 8.0) return 1.25;
        if (t < 14.0) return -1.1;
        return 1.45 * Math.sin(0.55 * (t - 14.0));
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Config cfg = new Config();
        SafetyAdaptiveController controller = new SafetyAdaptiveController(cfg);
        Path output = Path.of("Chapter24_Lesson5_telemetry_java.csv");

        double trueA = -0.45;
        double trueB = 1.0;
        double x = 0.0;
        double peakAbsX = 0.0;
        int steps = (int) (cfg.duration / cfg.dt);

        try (BufferedWriter out = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            out.write("time,x,x_m,error,reference,reference_limited,u_nominal,u_applied,"
                    + "theta_x,theta_r,adaptation_gate,safety_active,constrained,sensor_valid,status\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * cfg.dt;
                double r = referenceSignal(t);
                boolean sensorValid = !(t >= 12.0 && t < 12.15);
                double measuredX = sensorValid ? x : Double.NaN;
                StepResult data = controller.update(measuredX, r, sensorValid);

                double disturbance = 0.12 * Math.sin(1.7 * t);
                if (t >= 9.0 && t < 9.4) disturbance += 0.55;
                double xDot = trueA * x + trueB * data.uApplied() + disturbance;
                x += cfg.dt * xDot;
                peakAbsX = Math.max(peakAbsX, Math.abs(x));

                out.write(String.format(Locale.ROOT,
                        "%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.9f,%.1f,%b,%b,%b,%s%n",
                        t, x, data.xModel(), data.error(), r, data.referenceLimited(),
                        data.uNominal(), data.uApplied(), data.thetaX(), data.thetaR(),
                        data.adaptationGate(), data.safetyActive(), data.constrained(),
                        data.sensorValid(), data.status()));
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }

        System.out.printf(Locale.ROOT, "Peak |x|: %.6f (safe boundary %.3f)%n", peakAbsX, cfg.xSafe);
        System.out.printf(Locale.ROOT, "Final theta_x=%.6f, theta_r=%.6f%n",
                controller.thetaX, controller.thetaR);
        System.out.println("Safety interventions: " + controller.safetyInterventions);
        System.out.println("Actuator/rate interventions: " + controller.saturationEvents);
        System.out.println("Fallback samples: " + controller.fallbackEvents);

        if (peakAbsX > cfg.xSafe + 1.0e-6) System.exit(2);
    }
}

12. MATLAB and Simulink Implementation

The MATLAB script performs the complete adaptive simulation, exports telemetry, asserts the state envelope, and plots the principal signals. When Simulink is licensed, it also generates a reusable actuator safety shell containing magnitude limiting, slew-rate limiting, and command telemetry blocks.

Chapter24_Lesson5.m

% Chapter24_Lesson5.m
% Safety-supervised adaptive control demonstration.
% The script also creates a small Simulink safety-shell model when Simulink is
% available.  The complete adaptive simulation remains executable in MATLAB.

clear; clc; close all;

cfg.dt = 0.005;
cfg.duration = 20.0;
cfg.uMin = -2.0;
cfg.uMax = 2.0;
cfg.duMax = 8.0;
cfg.xSafe = 2.4;
cfg.referenceLimit = 1.6;
cfg.cbfRate = 3.0;
cfg.driftBound = 1.25;
cfg.inputGainMin = 0.8;
cfg.gammaX = 1.5;
cfg.gammaR = 1.0;
cfg.sigma = 0.03;
cfg.thetaXMin = -4.0;
cfg.thetaXMax = -0.05;
cfg.thetaRMin = 0.0;
cfg.thetaRMax = 4.0;
cfg.thetaRateMax = 1.0;
cfg.fallbackGain = 1.8;
cfg.tolerance = 1.0e-7;

N = floor(cfg.duration / cfg.dt) + 1;
time = (0:N-1)' * cfg.dt;
x = zeros(N,1);
xm = zeros(N,1);
errorSignal = zeros(N,1);
reference = zeros(N,1);
referenceLimited = zeros(N,1);
uNominal = zeros(N,1);
uApplied = zeros(N,1);
thetaX = zeros(N,1); thetaX(1) = -1.0;
thetaR = zeros(N,1); thetaR(1) = 1.0;
adaptationGate = zeros(N,1);
safetyActive = false(N,1);
constrained = false(N,1);
sensorValidLog = true(N,1);
status = strings(N,1);

trueA = -0.45;
trueB = 1.0;
previousU = 0.0;
safetyInterventions = 0;
saturationEvents = 0;
fallbackSamples = 0;

for k = 1:N-1
    t = time(k);
    reference(k) = referenceSignal(t);
    sensorValid = ~(t >= 12.0 && t < 12.15);
    sensorValidLog(k) = sensorValid;
    measuredX = x(k);
    if ~sensorValid
        measuredX = NaN;
    end

    r = clampValue(reference(k), -cfg.referenceLimit, cfg.referenceLimit);
    referenceLimited(k) = r;
    xm(k+1) = xm(k) + cfg.dt * (-1.2*xm(k) + 1.2*r);

    valid = sensorValid && isfinite(measuredX) && isfinite(reference(k));
    if ~valid
        fallbackSamples = fallbackSamples + 1;
        safeMeasurement = 0.0;
        if isfinite(measuredX), safeMeasurement = measuredX; end
        nominal = clampValue(-cfg.fallbackGain * safeMeasurement, cfg.uMin, cfg.uMax);
        filtered = nominal;
        gate = 0.0;
        active = false;
        status(k) = "FALLBACK_INVALID_SENSOR";
    else
        nominal = thetaX(k)*measuredX + thetaR(k)*r;
        [filtered, active] = cbfFilter(measuredX, nominal, cfg);
        gate = double(~active);
        if active
            safetyInterventions = safetyInterventions + 1;
            status(k) = "SAFETY_FILTER";
        else
            status(k) = "ADAPTIVE";
        end
    end

    magnitudeLimited = clampValue(filtered, cfg.uMin, cfg.uMax);
    magnitudeSaturated = abs(magnitudeLimited - filtered) > cfg.tolerance;
    maxStep = cfg.duMax * cfg.dt;
    applied = clampValue(magnitudeLimited, previousU-maxStep, previousU+maxStep);
    rateLimited = abs(applied - magnitudeLimited) > cfg.tolerance;
    previousU = applied;

    if magnitudeSaturated || rateLimited
        saturationEvents = saturationEvents + 1;
        gate = 0.0;
        if status(k) == "ADAPTIVE", status(k) = "ACTUATOR_LIMIT"; end
    end

    e = 0.0;
    if valid, e = measuredX - xm(k); end
    rawThetaXDot = -gate*cfg.gammaX*measuredX*e - cfg.sigma*thetaX(k);
    rawThetaRDot = -gate*cfg.gammaR*r*e - cfg.sigma*thetaR(k);
    if ~isfinite(rawThetaXDot), rawThetaXDot = -cfg.sigma*thetaX(k); end
    if ~isfinite(rawThetaRDot), rawThetaRDot = -cfg.sigma*thetaR(k); end
    thetaXDot = clampValue(rawThetaXDot, -cfg.thetaRateMax, cfg.thetaRateMax);
    thetaRDot = clampValue(rawThetaRDot, -cfg.thetaRateMax, cfg.thetaRateMax);

    thetaX(k+1) = clampValue(thetaX(k)+cfg.dt*thetaXDot, cfg.thetaXMin, cfg.thetaXMax);
    thetaR(k+1) = clampValue(thetaR(k)+cfg.dt*thetaRDot, cfg.thetaRMin, cfg.thetaRMax);

    disturbance = 0.12*sin(1.7*t);
    if t >= 9.0 && t < 9.4, disturbance = disturbance + 0.55; end
    xDot = trueA*x(k) + trueB*applied + disturbance;
    x(k+1) = x(k) + cfg.dt*xDot;

    errorSignal(k) = e;
    uNominal(k) = nominal;
    uApplied(k) = applied;
    adaptationGate(k) = gate;
    safetyActive(k) = active;
    constrained(k) = active || magnitudeSaturated || rateLimited;
end

reference(end) = referenceSignal(time(end));
referenceLimited(end) = clampValue(reference(end), -cfg.referenceLimit, cfg.referenceLimit);
errorSignal(end) = x(end)-xm(end);
uNominal(end) = uNominal(end-1);
uApplied(end) = uApplied(end-1);
adaptationGate(end) = adaptationGate(end-1);
safetyActive(end) = safetyActive(end-1);
constrained(end) = constrained(end-1);
sensorValidLog(end) = sensorValidLog(end-1);
status(end) = status(end-1);

telemetry = table(time,x,xm,errorSignal,reference,referenceLimited,uNominal,uApplied, ...
    thetaX,thetaR,adaptationGate,safetyActive,constrained,sensorValidLog,status);
writetable(telemetry, 'Chapter24_Lesson5_telemetry_matlab.csv');

fprintf('Peak |x|: %.6f (safe boundary %.3f)\n', max(abs(x)), cfg.xSafe);
fprintf('Final theta_x=%.6f, theta_r=%.6f\n', thetaX(end), thetaR(end));
fprintf('Safety interventions: %d\n', safetyInterventions);
fprintf('Actuator/rate interventions: %d\n', saturationEvents);
fprintf('Fallback samples: %d\n', fallbackSamples);
assert(max(abs(x)) <= cfg.xSafe + 1.0e-6, 'Safety boundary was exceeded.');

figure('Name','Chapter24 Lesson5 Safety-Supervised Adaptive Control');
plot(time, x, time, xm, '--', time, referenceLimited, ':', 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('State / reference');
legend('x','x_m','r_{limited}','Location','best');

figure('Name','Control and Adaptation');
plot(time, uApplied, time, uNominal, '--', 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('Control');
legend('u_{applied}','u_{nominal}','Location','best');

if license('test','Simulink')
    buildSimulinkSafetyShell(cfg);
end

function y = clampValue(value, lower, upper)
    y = max(lower, min(upper, value));
end

function r = referenceSignal(t)
    if t < 2.0
        r = 0.0;
    elseif t < 8.0
        r = 1.25;
    elseif t < 14.0
        r = -1.1;
    else
        r = 1.45*sin(0.55*(t-14.0));
    end
end

function [filtered, active] = cbfFilter(x, nominal, cfg)
    lower = cfg.uMin;
    upper = cfg.uMax;
    h = cfg.xSafe^2 - x^2;
    epsValue = 1.0e-9;
    if x > epsValue
        cbfUpper = (cfg.cbfRate*h/(2.0*x)-cfg.driftBound)/cfg.inputGainMin;
        upper = min(upper, cbfUpper);
    elseif x < -epsValue
        q = -x;
        cbfLower = (-cfg.cbfRate*h/(2.0*q)+cfg.driftBound)/cfg.inputGainMin;
        lower = max(lower, cbfLower);
    end
    if lower > upper
        if x > 0.0, filtered = cfg.uMin; else, filtered = cfg.uMax; end
        active = true;
        return;
    end
    filtered = clampValue(nominal, lower, upper);
    active = abs(filtered-nominal) > cfg.tolerance;
end

function buildSimulinkSafetyShell(cfg)
    model = 'Chapter24_Lesson5_SafetyShell';
    if bdIsLoaded(model), close_system(model,0); end
    if exist([model '.slx'],'file'), delete([model '.slx']); end
    new_system(model); open_system(model);

    add_block('simulink/Sources/In1',[model '/Nominal adaptive command'], ...
        'Position',[35 55 65 75]);
    add_block('simulink/Discontinuities/Saturation',[model '/Magnitude limit'], ...
        'UpperLimit',num2str(cfg.uMax),'LowerLimit',num2str(cfg.uMin), ...
        'Position',[125 45 195 85]);
    add_block('simulink/Discontinuities/Rate Limiter',[model '/Slew-rate limit'], ...
        'RisingSlewLimit',num2str(cfg.duMax),'FallingSlewLimit',num2str(-cfg.duMax), ...
        'Position',[245 45 325 85]);
    add_block('simulink/Sinks/Out1',[model '/Applied command'], ...
        'Position',[385 55 415 75]);
    add_block('simulink/Sinks/To Workspace',[model '/Command telemetry'], ...
        'VariableName','uAppliedLog','SaveFormat','Structure With Time', ...
        'Position',[365 115 455 145]);

    add_line(model,'Nominal adaptive command/1','Magnitude limit/1');
    add_line(model,'Magnitude limit/1','Slew-rate limit/1');
    add_line(model,'Slew-rate limit/1','Applied command/1');
    add_line(model,'Slew-rate limit/1','Command telemetry/1','autorouting','on');
    set_param(model,'StopTime','20');
    save_system(model);
    fprintf('Created Simulink shell: %s.slx\n', model);
end

13. Wolfram Mathematica Implementation

The notebook source contains an executable Wolfram Language simulation, CSV export, safety-envelope report, and state/reference plot.

Chapter24_Lesson5.nb


Notebook[{
 Cell["Chapter 24, Lesson 5: Safety-Supervised Adaptive Control", "Title"],
 Cell["This notebook implements finite-value checks, reference limiting, a scalar CBF safety filter, actuator magnitude and rate limits, saturation-aware adaptation, sigma leakage, projection, fallback logic, and CSV telemetry.", "Text"],
 Cell[BoxData["ClearAll[\"Global`*\"];
cfg = <|\"dt\" -> 0.005, \"duration\" -> 20.0, \"uMin\" -> -2.0, \"uMax\" -> 2.0, \"duMax\" -> 8.0, \"xSafe\" -> 2.4, \"referenceLimit\" -> 1.6,
  \"cbfRate\" -> 3.0, \"driftBound\" -> 1.25, \"inputGainMin\" -> 0.8, \"gammaX\" -> 1.5, \"gammaR\" -> 1.0, \"sigma\" -> 0.03,
  \"thetaXMin\" -> -4.0, \"thetaXMax\" -> -0.05, \"thetaRMin\" -> 0.0, \"thetaRMax\" -> 4.0, \"thetaRateMax\" -> 1.0, \"fallbackGain\" -> 1.8,
  \"tolerance\" -> 10^-7|>;
clamp[v_, lo_, hi_] := Max[lo, Min[hi, v]];
referenceSignal[t_] := Piecewise[{ {0.0, t < 2.0}, {1.25, t < 8.0}, {-1.1, t < 14.0} }, 1.45 Sin[0.55 (t - 14.0)]];
cbfFilter[x_, nominal_] := Module[{lower = cfg[\"uMin\"], upper = cfg[\"uMax\"],
   h, eps = 10^-9, cbfUpper, cbfLower, filtered},
  h = cfg[\"xSafe\"]^2 - x^2;
  If[x > eps, cbfUpper = (cfg[\"cbfRate\"] h/(2 x) - cfg[\"driftBound\"])/cfg[\"inputGainMin\"];
   upper = Min[upper, cbfUpper], If[x < -eps,
    cbfLower = (-cfg[\"cbfRate\"] h/(2 (-x)) + cfg[\"driftBound\"])/cfg[\"inputGainMin\"];
    lower = Max[lower, cbfLower]]];
  If[lower > upper, Return[{If[x > 0, cfg[\"uMin\"], cfg[\"uMax\"]], True}]];
  filtered = clamp[nominal, lower, upper];
  {filtered, Abs[filtered - nominal] > cfg[\"tolerance\"]}];
n = Floor[cfg[\"duration\"]/cfg[\"dt\"]] + 1; time = N[Range[0, n - 1] cfg[\"dt\"]];
x = ConstantArray[0.0, n]; xm = ConstantArray[0.0, n];
thetaX = ConstantArray[-1.0, n]; thetaR = ConstantArray[1.0, n];
uApplied = ConstantArray[0.0, n]; uNominal = ConstantArray[0.0, n];
reference = referenceSignal /@ time; gate = ConstantArray[0.0, n];
safety = ConstantArray[False, n]; validLog = ConstantArray[True, n];
status = ConstantArray[\"\", n]; previousU = 0.0; trueA = -0.45; trueB = 1.0;
Do[t = time[[k]];
 r = clamp[reference[[k]], -cfg[\"referenceLimit\"], cfg[\"referenceLimit\"]];
 sensorValid = Not[12.0 <= t < 12.15]; validLog[[k]] = sensorValid;
 measuredX = If[sensorValid, x[[k]], Indeterminate];
 xm[[k + 1]] = xm[[k]] + cfg[\"dt\"] (-1.2 xm[[k]] + 1.2 r);
 valid = sensorValid && NumberQ[measuredX] && Im[measuredX] == 0;
 If[! valid, nominal = clamp[-cfg[\"fallbackGain\"] If[NumberQ[measuredX], measuredX, 0.0], cfg[\"uMin\"], cfg[\"uMax\"]];
  filtered = nominal; active = False; g = 0.0; status[[k]] = \"FALLBACK_INVALID_SENSOR\",
  nominal = thetaX[[k]] measuredX + thetaR[[k]] r; {filtered, active} = cbfFilter[measuredX, nominal];
  g = If[active, 0.0, 1.0]; status[[k]] = If[active, \"SAFETY_FILTER\", \"ADAPTIVE\"]];
 magnitudeLimited = clamp[filtered, cfg[\"uMin\"], cfg[\"uMax\"]]; maxStep = cfg[\"duMax\"] cfg[\"dt\"];
 applied = clamp[magnitudeLimited, previousU - maxStep, previousU + maxStep];
 constrained = active || Abs[magnitudeLimited - filtered] > cfg[\"tolerance\"] || Abs[applied - magnitudeLimited] > cfg[\"tolerance\"];
 If[constrained, g = 0.0; If[status[[k]] == \"ADAPTIVE\", status[[k]] = \"ACTUATOR_LIMIT\"]];
 previousU = applied; e = If[valid, measuredX - xm[[k]], 0.0];
 txDot = clamp[-g cfg[\"gammaX\"] If[valid, measuredX, 0.0] e - cfg[\"sigma\"] thetaX[[k]], -cfg[\"thetaRateMax\"], cfg[\"thetaRateMax\"]];
 trDot = clamp[-g cfg[\"gammaR\"] r e - cfg[\"sigma\"] thetaR[[k]], -cfg[\"thetaRateMax\"], cfg[\"thetaRateMax\"]];
 thetaX[[k + 1]] = clamp[thetaX[[k]] + cfg[\"dt\"] txDot, cfg[\"thetaXMin\"], cfg[\"thetaXMax\"]];
 thetaR[[k + 1]] = clamp[thetaR[[k]] + cfg[\"dt\"] trDot, cfg[\"thetaRMin\"], cfg[\"thetaRMax\"]];
 disturbance = 0.12 Sin[1.7 t] + If[9.0 <= t < 9.4, 0.55, 0.0];
 x[[k + 1]] = x[[k]] + cfg[\"dt\"] (trueA x[[k]] + trueB applied + disturbance);
 uNominal[[k]] = nominal; uApplied[[k]] = applied; gate[[k]] = g; safety[[k]] = active,
 {k, 1, n - 1}];
headers = {\"time\", \"x\", \"x_m\", \"reference\", \"u_nominal\", \"u_applied\", \"theta_x\", \"theta_r\", \"adaptation_gate\", \"safety_active\", \"sensor_valid\", \"status\"};
rows = Transpose[{time, x, xm, reference, uNominal, uApplied, thetaX, thetaR, gate, Boole /@ safety, Boole /@ validLog, status}];
Export[\"Chapter24_Lesson5_telemetry_mathematica.csv\", Prepend[rows, headers]];
Print[\"Peak |x|: \", Max[Abs[x]], \" (safe boundary \", cfg[\"xSafe\"], \")\"];
Print[\"Final theta_x=\", Last[thetaX], \", theta_r=\", Last[thetaR]];
ListLinePlot[{Transpose[{time, x}], Transpose[{time, xm}], Transpose[{time, reference}]},
 PlotLegends -> {\"x\", \"x_m\", \"reference\"}, Frame -> True,
 FrameLabel -> {\"Time (s)\", \"State / reference\"}]
"], "Input"]
}, WindowTitle -> "Chapter24_Lesson5", StyleDefinitions -> "Default.nb"]        

14. Practical Acceptance Criteria

A release candidate should satisfy quantitative criteria defined before testing. Typical examples are:

\[ \max_t |x(t)|\leq x_{\max}, \qquad \max_t |u_a(t)|\leq u_{\max}, \qquad \max_k \frac{|u_a[k]-u_a[k-1]|}{T_s}\leq \dot u_{\max}. \]

\[ \hat{\boldsymbol{\theta} }(t)\in\Omega_\theta, \qquad \max_k\delta_\theta[k]\leq\delta_{\theta,\max}, \qquad N_{\mathrm{deadline} }=0 \text{ in the qualified operating region}. \]

In addition, every injected invalid sensor value must produce the documented fallback state within the allowed reaction time, every safety-filter intervention must be logged, and a replay using the stored configuration must reproduce the command history within a declared numerical tolerance.

15. Problems and Solutions

Problem 1 (Ultimate Bound Under Saturation): Suppose

\[ \dot V\leq-\lambda_Q\lVert\mathbf{e}\rVert^2 +2c_B\lVert\mathbf{e}\rVert|\Delta u|, \qquad |\Delta u|\leq\bar u_\Delta. \]

Derive a region outside which \(\dot V<0\).

Solution: Factor the right-hand side:

\[ \dot V\leq-\lVert\mathbf{e}\rVert \left(\lambda_Q\lVert\mathbf{e}\rVert-2c_B\bar u_\Delta\right). \]

Therefore \(\dot V<0\) whenever

\[ \lVert\mathbf{e}\rVert > \frac{2c_B\bar u_\Delta}{\lambda_Q}. \]

The saturation mismatch creates an ultimate tracking-error radius proportional to its worst-case magnitude.

Problem 2 (Scalar Barrier Constraint): For \(\dot x=ax+bu\), \(b>0\), and \(h=x_{\max}^2-x^2\), derive the admissible command when \(x>0\).

Solution: The condition \(\dot h+\kappa h\geq0\) gives

\[ -2x(ax+bu)+\kappa(x_{\max}^2-x^2)\geq0. \]

Since \(2xb>0\),

\[ u\leq\frac{\kappa(x_{\max}^2-x^2)-2ax^2}{2xb}. \]

The filter clips the nominal command to this upper bound and to the physical actuator interval.

Problem 3 (Why Freeze Adaptation?): A positive reference produces \(u_c>u_{\max}\) and a persistent positive tracking error. Explain why continuing the nominal gradient update may be harmful.

Solution: The error is partly caused by unavailable actuator authority, but the estimator attributes it to parameter mismatch. The update can increase the requested control magnitude, increasing \(|u_c-u_a|\) without changing the plant input. Freezing the tracking-driven term removes this false learning channel; leakage and projection continue to keep the estimate bounded.

Problem 4 (Discrete Parameter-Step Limit): Let \(T_s=0.01\) s and \(|\dot{\hat\theta}|\leq0.8\). What is the maximum change in one sample and in 100 consecutive samples if projection never activates?

Solution:

\[ |\Delta\hat\theta[k]|\leq T_s\dot\theta_{\max}=0.008. \]

Over 100 samples, the worst-case accumulated motion is \(100(0.008)=0.8\). This shows why both rate and absolute parameter bounds are needed.

Problem 5 (Traceability Record): A hardware test log stores states and commands but not the controller configuration or software revision. Can the test be considered reproducible?

Solution: No. The same measured states can produce different commands under different gains, bounds, fallback modes, calibration files, or code revisions. A complete run record must bind telemetry to immutable software, configuration, calibration, scenario, and timing identifiers.

Problem 6 (Safe Re-entry): Propose a logical condition for returning from fallback to adaptive operation.

Solution: One conservative condition is

\[ q_{\mathrm{reentry} } = q_{\mathrm{sensor} }\land q_{\mathrm{timing} }\land q_{\mathrm{actuator} }\land q_{\theta}\land q_{e}\land q_{\mathrm{dwell} }. \]

Here the flags require valid and fresh sensors, acceptable timing, actuator headroom, estimates inside a tighter re-entry set, bounded tracking error, and a fault-free dwell interval. The transition should be bumpless or ramped.

16. Summary

A deployable adaptive controller requires more than a nominal update law. The implemented system must distinguish commanded and applied inputs, prevent saturation-induced parameter drift, enforce parameter and state envelopes, validate every signal, handle timing and numerical faults, provide a deterministic fallback, and log the complete decision context. Safety filters and actuator management should remain authoritative over adaptation. Documentation, configuration control, fault injection, and regression testing are mathematical implementation requirements because they determine which closed-loop system is actually being analyzed and operated.

17. References

  1. Ioannou, P.A., & Kokotovic, P.V. (1984). Instability analysis and improvement of robustness of adaptive control. Automatica, 20(5), 583–594.
  2. Narendra, K.S., & Annaswamy, A.M. (1987). A new adaptive law for robust adaptation without persistent excitation. IEEE Transactions on Automatic Control, 32(2), 134–145.
  3. Ioannou, P.A., & Tsakalis, K.S. (1986). A robust direct adaptive controller. IEEE Transactions on Automatic Control, 31(11), 1033–1043.
  4. Grimm, G., Hatfield, J., Postlethwaite, I., Teel, A.R., Turner, M.C., & Zaccarian, L. (2003). Antiwindup for stable linear systems with input saturation: An LMI-based synthesis. IEEE Transactions on Automatic Control, 48(9), 1509–1525.
  5. Teel, A.R. (1992). Global stabilization and restricted tracking for multiple integrators with bounded controls. Systems & Control Letters, 18(3), 165–171.
  6. Ames, A.D., Xu, X., Grizzle, J.W., & Tabuada, P. (2017). Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control, 62(8), 3861–3876.
  7. Xu, X., Tabuada, P., Grizzle, J.W., & Ames, A.D. (2018). Robustness of control barrier functions for safety critical control. IFAC-PapersOnLine, 51(34), 54–61.
  8. Morse, A.S. (1996). Supervisory control of families of linear set-point controllers—Part 1: Exact matching. IEEE Transactions on Automatic Control, 41(10), 1413–1431.
  9. Hespanha, J.P., Liberzon, D., & Morse, A.S. (2003). Hysteresis-based switching algorithms for supervisory control of uncertain systems. Automatica, 39(2), 263–272.
Support CaaT Academy

Help keep these engineering tutorials free and growing

If these lessons, examples, and project pages help you, a small donation supports the continued creation and improvement of free control, robotics, software, and engineering education resources.

Created and maintained by Abolfazl Mohammadijoo.