A continuous-time adaptive controller is usually analyzed as an interconnection of differential equations. A processor, however, measures signals only at sampling instants, computes for a nonzero time, updates parameters in finite increments, and holds the actuator command between updates. The implemented loop is therefore a hybrid sampled-data system, not merely the original equations evaluated by a computer.

Let the continuous plant be \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u},\boldsymbol{\theta}) \), with adaptive controller and update law

\[ \mathbf{u}(t)=\boldsymbol{\kappa}(\mathbf{x}(t),\hat{\boldsymbol{\theta}}(t),r(t)), \qquad \dot{\hat{\boldsymbol{\theta}}}(t)= \mathbf{g}(\mathbf{x}(t),e(t),r(t)). \]

In a periodic implementation with sample times \( t_k \) and intervals \( T_k=t_{k+1}-t_k \), the actuator commonly applies

\[ \mathbf{u}(t)=\mathbf{u}_k, \qquad t\in[t_k,t_{k+1}), \qquad \mathbf{u}_k= \boldsymbol{\kappa}(\mathbf{x}_{k-d_x}, \hat{\boldsymbol{\theta}}_{k-d_\theta},r_{k-d_r}). \]

The integers \(d_x,d_\theta,d_r\) represent measurement, scheduling, communication, or pipeline delays. Even when all are zero, the held control creates a sampling error because the continuous feedback value changes inside the interval while \(\mathbf{u}_k\) does not.

\[ \boldsymbol{\varepsilon}_s(t)= \boldsymbol{\kappa}(\mathbf{x}(t_k),\hat{\boldsymbol{\theta}}(t_k),r(t_k))- \boldsymbol{\kappa}(\mathbf{x}(t),\hat{\boldsymbol{\theta}}(t),r(t)). \]

Thus the digital controller introduces a structured perturbation whose magnitude depends on the sample interval, closed-loop bandwidth, adaptation gain, signal derivatives, and scheduling jitter.

Consider a continuous linear plant during one interval with held input \(\mathbf{u}(t)=\mathbf{u}_k\):

\[ \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}_k, \qquad t\in[t_k,t_k+T_k). \]

Variation of constants gives the exact sampled transition

\[ \mathbf{x}_{k+1}=\mathbf{A}_{d,k}\mathbf{x}_k+ \mathbf{B}_{d,k}\mathbf{u}_k, \quad \mathbf{A}_{d,k}=e^{\mathbf{A}T_k}, \quad \mathbf{B}_{d,k}=\int_0^{T_k}e^{\mathbf{A}\tau}\mathbf{B}\,d\tau. \]

For the scalar plant \(\dot{x}=-a x+b u\), \(a>0\),

\[ x_{k+1}=e^{-aT_k}x_k+ \frac{b}{a}\left(1-e^{-aT_k}\right)u_k. \]

Forward Euler instead uses \(x_{k+1}^{E}=(1-aT_k)x_k+bT_k u_k\). Expanding the exact model,

\[ e^{-aT_k}=1-aT_k+\frac{a^2T_k^2}{2}+\mathcal{O}(T_k^3), \qquad \frac{1-e^{-aT_k}}{a}=T_k-\frac{aT_k^2}{2}+\mathcal{O}(T_k^3). \]

Therefore the one-step Euler error is

\[ x_{k+1}-x_{k+1}^{E}= \frac{T_k^2}{2}\left(a^2x_k-ab u_k\right)+\mathcal{O}(T_k^3). \]

The local truncation error is second order, but the accumulated global error over a fixed horizon is first order in the maximum step. In an adaptive loop, this is not merely a plotting error: the distorted state and error signals feed the parameter update and can alter the controller itself.

A minimal scalar stability check illustrates the danger. Euler discretization of \(\dot{x}=-a x\) gives \(x_{k+1}=(1-aT_s)x_k\). Discrete asymptotic stability requires

\[ |1-aT_s|<1 \quad\Longleftrightarrow\quad 0<T_s<\frac{2}{a}. \]

The continuous system is stable for every \(a>0\), whereas its Euler implementation becomes oscillatory and then unstable when the numerical step is too large.

To make the implementation effects explicit, consider the scalar uncertain plant and reference model

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

with direct adaptive control

\[ u=\hat{\theta}_x x+\hat{\theta}_r r. \]

For known positive control direction \(b>0\), the ideal matching parameters are

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

Defining \(e=x-x_m\), \(\tilde{\boldsymbol{\theta}}= [\hat{\theta}_x-\theta_x^\star,\hat{\theta}_r-\theta_r^\star]^T\), and \(\boldsymbol{\phi}=[x,r]^T\), the continuous error model is

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

Under zero-order hold with constant \(u_k\) and \(r_k\) during the interval, exact integration gives

\[ e_{k+1}=A_m e_k+B_m b\, \tilde{\boldsymbol{\theta}}_k^T\boldsymbol{\phi}_k+\eta_k, \quad A_m=e^{-a_mT_s}, \quad B_m=\frac{1-e^{-a_mT_s}}{a_m}, \]

where \(\eta_k\) collects intersample regressor variation, discretization mismatch, quantization, disturbance, saturation, and delay effects. Continuous-time proofs effectively set this perturbation to zero. Real code must either bound it or monitor the conditions that keep it small.

A continuous gradient update for the scalar MRAC is

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

A direct Euler implementation is

\[ \hat{\boldsymbol{\theta}}_{k+1}= \hat{\boldsymbol{\theta}}_k-T_s\boldsymbol{\Gamma} \boldsymbol{\phi}_k e_k. \]

The factor \(T_s\) is essential. Omitting it changes the effective continuous adaptation gain from \(\boldsymbol{\Gamma}\) to roughly \(\boldsymbol{\Gamma}/T_s\); changing the sampling rate would then change the controller aggressiveness by the inverse ratio.

A normalized and projected implementation is usually safer:

\[ m_k^2=1+\nu\boldsymbol{\phi}_k^T\boldsymbol{\phi}_k, \qquad \hat{\boldsymbol{\theta}}_{k+1}= \operatorname{Proj}_{\Omega}\left[ \hat{\boldsymbol{\theta}}_k- T_s\boldsymbol{\Gamma}\frac{\boldsymbol{\phi}_k e_k}{m_k^2} \right]. \]

Normalization limits the update produced by large regressors; projection enforces known parameter or controller-gain bounds. Neither device compensates for an arbitrarily large sampling interval or delay, but both prevent common numerical runaway mechanisms.

4.1 Why the Continuous Lyapunov Cancellation Is Not Exact

In continuous time, choose

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

Substitution of the error and update equations cancels the mixed term and yields \(\dot{V}=-a_m e^2\). For a digital step, write generically

\[ e_{k+1}=e_k+T_s f_e(e_k,\tilde{\boldsymbol{\theta}}_k,\boldsymbol{\phi}_k) +\mathcal{O}(T_s^2), \] \[ \tilde{\boldsymbol{\theta}}_{k+1}= \tilde{\boldsymbol{\theta}}_k- T_s\boldsymbol{\Gamma}\boldsymbol{\phi}_k e_k. \]

Expanding the Lyapunov difference gives

\[ \Delta V_k=V_{k+1}-V_k =-T_s a_m e_k^2+T_s^2\Psi_k, \]

where \(\Psi_k\) contains squared state and parameter increments, intersample variation, and delay terms. If, on an operating set, \(|\Psi_k|\le C\|\mathbf{z}_k\|^2\) and \(a_m e_k^2\ge c\|\mathbf{z}_k\|^2\) for the component being dissipated, then

\[ \Delta V_k\le -T_s(c-T_s C)\|\mathbf{z}_k\|^2. \]

A sufficient small-step condition is therefore \(0<T_s<c/C\). This argument explains a general engineering fact: a continuous adaptive proof often becomes a semiglobal practical sampled-data result whose admissible sampling period depends on bounds for states, regressors, gains, and unmodeled dynamics.

Suppose the sampled scalar plant is \(x_{k+1}=A_d x_k+B_d u_k\), while the controller uses a delayed state:

\[ u_k=K x_{k-d}. \]

For one-sample delay, define \(\boldsymbol{\xi}_k=[x_k,x_{k-1}]^T\). Then

\[ \boldsymbol{\xi}_{k+1}= \begin{bmatrix} A_d & B_d K\\ 1 & 0 \end{bmatrix} \boldsymbol{\xi}_k. \]

Stability is determined by the eigenvalues of this augmented matrix, not by \(A_d+B_dK\). For an arbitrary integer delay \(d\), the characteristic polynomial becomes

\[ \lambda^{d+1}-A_d\lambda^d-B_dK=0. \]

Increasing delay adds roots and can move them outside the unit disk. In frequency terms, a pure delay \(\tau\) contributes phase \(-\omega\tau\). A rough fixed-controller margin check is

\[ \omega_c\tau<\phi_m, \]

where \(\omega_c\) is crossover frequency and \(\phi_m\) is available phase margin in radians. Adaptive gains may increase effective bandwidth, so a delay that was harmless initially can become critical after parameter transients.

A physical delay rarely equals an integer number of samples. Write \(\tau=dT_s+\delta\), with \(0\le\delta<T_s\). Rounding \(\tau/T_s\) discards the fractional delay. Higher-fidelity implementations use timestamped data, interpolation, state prediction, a delay-line model, or a controller designed directly for sampled-data delay dynamics.

flowchart TD
  S["Sensor acquisition at t_k"] --> Q["ADC and signal conditioning"]
  Q --> E["State and error calculation"]
  E --> A["Adaptive parameter update"]
  A --> C["Control calculation and limits"]
  C --> D["DAC or communication delay"]
  D --> H["Zero-order hold on actuator"]
  H --> P["Physical plant during sample interval"]
  P --> S
        

Define the actual interval by \(T_k=T_s+\delta_k\), where \(\delta_k\) is sampling jitter. If \(|\delta_k|\le\bar{\delta}\), then \(T_k\in[T_{\min},T_{\max}]\). A variable-step Lyapunov estimate often has the form

\[ V_{k+1}-V_k\le -\alpha T_k\|\mathbf{z}_k\|^2+ \beta T_k^2\|\mathbf{z}_k\|^2. \]

Uniform decrease follows from the sufficient condition

\[ T_{\max}<\frac{\alpha}{\beta}. \]

Consequently, testing only the average execution period is inadequate. The relevant quantities are worst-case execution time, worst-case sensor-to-actuator latency, maximum interval, deadline-miss policy, and the age of each signal used by the update law.

When a deadline is missed, do not silently perform multiple large parameter updates with stale data. A defined policy should either hold the previous command and freeze adaptation, execute a bounded update using the measured elapsed time, or enter a certified fallback controller.

7.1 Separate Controller Sampling from Simulation Integration

In simulation, the controller sample time and numerical integration step are different design quantities. Let \(T_s\) be the controller period and \(h=T_s/M\) the internal plant-integration step. The control must remain held across the \(M\) substeps. Recomputing the controller at each solver substep accidentally simulates a faster controller and hides digital implementation problems.

7.2 Euler, RK4, and Exact Linear Propagation

For a smooth ODE, forward Euler has global error \(\mathcal{O}(h)\), while classical RK4 has global error \(\mathcal{O}(h^4)\). For a linear model under zero-order hold, matrix-exponential propagation is exact up to floating-point and matrix-exponential computation errors. RK4 is often a practical reference for nonlinear simulation, but it does not remove sampling or delay errors.

7.3 Stiffness Created by Fast Adaptation

Large adaptation gains create a fast parameter time scale next to slower plant dynamics. A simplified two-time-scale model is

\[ \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\hat{\boldsymbol{\theta}}), \qquad \epsilon\dot{\hat{\boldsymbol{\theta}}}=\mathbf{g}(\mathbf{x},e), \qquad 0<\epsilon\ll1. \]

Explicit solvers then require a step small enough for the fast mode, even when the visible plant trajectory is slow. Symptoms include alternating parameter estimates, NaN/Inf values, solver-dependent conclusions, and apparently improved performance when the simulation step is reduced. In offline analysis, compare multiple solver tolerances or use a stiff solver when justified. In embedded code, reduce adaptation gain, normalize updates, apply projection, and bound the parameter increment per sample.

7.4 Parameter Increment Limiting

A practical update may additionally enforce

\[ \Delta\hat{\boldsymbol{\theta}}_k^{raw}= -T_s\boldsymbol{\Gamma}\frac{\boldsymbol{\phi}_ke_k}{m_k^2}, \qquad \Delta\hat{\boldsymbol{\theta}}_k= \operatorname{sat}_{\Delta\theta_{\max}} \left(\Delta\hat{\boldsymbol{\theta}}_k^{raw}\right). \]

Increment limiting is a software safety device, not a substitute for stability analysis. It must be represented in validation models because it changes the adaptation dynamics.

With uniform quantization step \(q\), the measured signal is \(y_k^q=y_k+n_k^q\), where ideally \(|n_k^q|\le q/2\). The adaptive increment then contains a persistent perturbation:

\[ \Delta\hat{\boldsymbol{\theta}}_k= -T_s\boldsymbol{\Gamma}\boldsymbol{\phi}_k(e_k+n_k^q). \]

Even when \(e_k\) is small, quantization can keep parameters moving. A dead zone

\[ \mathcal{D}(e_k)= \begin{cases} 0, & |e_k|\le e_0,\\ e_k-e_0\operatorname{sgn}(e_k), & |e_k|>e_0, \end{cases} \]

with \(e_0\) chosen above the combined noise and quantization floor, can prevent noise-driven drift. The choice must remain consistent with the robust modification concepts introduced earlier in the course.

Floating-point implementations must also guard against overflow, loss of symmetry, loss of positive definiteness, and subtraction cancellation. For RLS-based adaptive controllers, the covariance matrix should be symmetrized numerically,

\[ \mathbf{P}_k\leftarrow\frac{1}{2} \left(\mathbf{P}_k+\mathbf{P}_k^T\right), \]

while square-root RLS or QR-based updates are preferable when conditioning is poor. Testing should log the minimum eigenvalue and condition number of \(\mathbf{P}_k\), not only the tracking error.

A frequent implementation error is adapting as though the requested input was applied when the actuator actually saturated. Let

\[ u_k^{cmd}=\hat{\boldsymbol{\theta}}_k^T\boldsymbol{\phi}_k, \qquad u_k^{act}=\operatorname{sat}_{u_{\max}}(u_k^{cmd}). \]

The mismatch \(u_k^{act}-u_k^{cmd}\) violates the nominal matching model. If the update law interprets the resulting tracking error as parametric uncertainty, it may drive gains to their limits. Common policies are adaptation freeze during saturation, an error modification based on the saturation mismatch, or a reference governor that prevents persistent infeasible commands.

The per-sample computation order should be explicit and deterministic: timestamp and validate measurements; reconstruct delayed signals; compute reference model and tracking error; evaluate the adaptation law; apply projection and increment limits; compute the command; apply actuator limits; publish the command; then record diagnostics. Changing this order changes which sample of each variable participates in the closed-loop map.

The following five implementations use the same uncertain plant, reference model, normalized gradient update, projection bounds, saturation-aware adaptation, delayed measurement queue, and RK4 substepping. Two cases are compared:

Nominal: \(T_s=0.01\,\text{s}\), zero delay, five RK4 substeps. Stressed: \(T_s=0.05\,\text{s}\), two-sample delay, one RK4 step, and measurement quantization. The experiment is pedagogical: it demonstrates the direction of degradation rather than certifying a universal maximum sample period.

\[ \dot{x}=-1.2x+u, \qquad \dot{x}_m=-2x_m+2r, \] \[ u_k=\operatorname{sat}_{8} \left(\hat{\theta}_{x,k}x_{k-d}+\hat{\theta}_{r,k}r_{k-d}\right). \]

Chapter24_Lesson3.py

Libraries: numpy for array calculations and matplotlib for plots. The adaptive algorithm and RK4 integrator are implemented directly rather than hidden in a control library.

"""Chapter24_Lesson3.py
Sampled-data MRAC experiment with delay, projection, normalization, saturation,
and RK4 plant integration.
"""
from dataclasses import dataclass
from collections import deque
import math
import numpy as np
import matplotlib.pyplot as plt


@dataclass
class Config:
    name: str
    sample_time: float
    delay_steps: int
    integration_substeps: int
    gamma_x: float = 4.0
    gamma_r: float = 4.0
    normalization: float = 0.25
    theta_limit: float = 8.0
    u_limit: float = 8.0
    quantization: float = 0.0
    final_time: float = 18.0


def reference(t: float) -> float:
    if t < 4.0:
        return 1.0
    if t < 8.0:
        return -0.5
    if t < 13.0:
        return 0.8
    return 0.2


def quantize(value: float, quantum: float) -> float:
    if quantum <= 0.0:
        return value
    return quantum * round(value / quantum)


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


def rk4_scalar(x: float, dt: float, rhs) -> float:
    k1 = rhs(x)
    k2 = rhs(x + 0.5 * dt * k1)
    k3 = rhs(x + 0.5 * dt * k2)
    k4 = rhs(x + dt * k3)
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(cfg: Config):
    # Unknown plant: x_dot = -a*x + b*u
    a, b = 1.2, 1.0
    # Reference model: xm_dot = -am*xm + bm*r
    am, bm = 2.0, 2.0

    steps = int(round(cfg.final_time / cfg.sample_time))
    t = np.arange(steps + 1, dtype=float) * cfg.sample_time
    x = np.zeros(steps + 1)
    xm = np.zeros(steps + 1)
    u = np.zeros(steps + 1)
    theta_x = np.zeros(steps + 1)
    theta_r = np.zeros(steps + 1)
    r_log = np.zeros(steps + 1)

    history = deque(maxlen=cfg.delay_steps + 1)
    for _ in range(cfg.delay_steps + 1):
        history.append((0.0, 0.0, reference(0.0)))

    saturated_samples = 0

    for k in range(steps):
        r_now = reference(t[k])
        r_log[k] = r_now
        x_meas = quantize(x[k], cfg.quantization)
        history.append((x_meas, xm[k], r_now))
        x_d, xm_d, r_d = history[0]

        e_d = x_d - xm_d
        phi_norm_sq = x_d * x_d + r_d * r_d
        normalizer = 1.0 + cfg.normalization * phi_norm_sq

        u_raw = theta_x[k] * x_d + theta_r[k] * r_d
        u[k] = min(cfg.u_limit, max(-cfg.u_limit, u_raw))
        saturated = abs(u_raw - u[k]) > 1.0e-12
        saturated_samples += int(saturated)

        # Freeze adaptation during saturation; otherwise use a normalized gradient step.
        if saturated:
            tx_next, tr_next = theta_x[k], theta_r[k]
        else:
            tx_next = theta_x[k] - cfg.sample_time * cfg.gamma_x * e_d * x_d / normalizer
            tr_next = theta_r[k] - cfg.sample_time * cfg.gamma_r * e_d * r_d / normalizer

        theta_x[k + 1] = project(tx_next, cfg.theta_limit)
        theta_r[k + 1] = project(tr_next, cfg.theta_limit)

        # The controller output is held constant over the complete sample interval.
        dt = cfg.sample_time / cfg.integration_substeps
        x_local, xm_local = x[k], xm[k]
        for j in range(cfg.integration_substeps):
            tj = t[k] + j * dt
            r_sub = reference(tj)
            x_local = rk4_scalar(x_local, dt, lambda z: -a * z + b * u[k])
            xm_local = rk4_scalar(xm_local, dt, lambda z: -am * z + bm * r_sub)
        x[k + 1] = x_local
        xm[k + 1] = xm_local

        if not all(math.isfinite(v) for v in (x[k + 1], xm[k + 1], theta_x[k + 1], theta_r[k + 1])):
            raise FloatingPointError(f"Non-finite value at sample {k}")

    r_log[-1] = reference(t[-1])
    u[-1] = u[-2]
    error = x - xm
    metrics = {
        "rms_error": float(np.sqrt(np.mean(error * error))),
        "max_error": float(np.max(np.abs(error))),
        "theta_x_final": float(theta_x[-1]),
        "theta_r_final": float(theta_r[-1]),
        "saturation_fraction": saturated_samples / max(1, steps),
    }
    return {
        "t": t,
        "x": x,
        "xm": xm,
        "u": u,
        "theta_x": theta_x,
        "theta_r": theta_r,
        "r": r_log,
        "error": error,
        "metrics": metrics,
    }


def main() -> None:
    cases = [
        Config("nominal", sample_time=0.01, delay_steps=0, integration_substeps=5),
        Config("stressed", sample_time=0.05, delay_steps=2, integration_substeps=1, quantization=0.002),
    ]
    results = {cfg.name: simulate(cfg) for cfg in cases}

    for name, result in results.items():
        print(name, result["metrics"])

    plt.figure(figsize=(10, 5))
    for name, result in results.items():
        plt.plot(result["t"], result["x"], label=f"x ({name})")
    plt.plot(results["nominal"]["t"], results["nominal"]["xm"], "--", label="reference model")
    plt.xlabel("time [s]")
    plt.ylabel("state")
    plt.title("Sampled-data MRAC: sampling and delay sensitivity")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(10, 5))
    for name, result in results.items():
        plt.plot(result["t"], result["error"], label=f"tracking error ({name})")
    plt.xlabel("time [s]")
    plt.ylabel("e = x - xm")
    plt.title("Tracking-error comparison")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Chapter24_Lesson3.cpp

This C++17 version uses only the standard library. It is suitable for studying the deterministic sample loop before porting the logic to an RTOS task or hardware abstraction layer.

// Chapter24_Lesson3.cpp
// Sampled-data MRAC with integer-sample delay, normalization, projection,
// saturation-aware adaptation, and RK4 integration. No external libraries.
#include <algorithm>
#include <cmath>
#include <deque>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <tuple>
#include <vector>

struct Config {
    std::string name;
    double sample_time;
    int delay_steps;
    int integration_substeps;
    double gamma_x = 4.0;
    double gamma_r = 4.0;
    double normalization = 0.25;
    double theta_limit = 8.0;
    double u_limit = 8.0;
    double quantization = 0.0;
    double final_time = 18.0;
};

struct Metrics {
    double rms_error{};
    double max_error{};
    double theta_x_final{};
    double theta_r_final{};
    double saturation_fraction{};
};

struct Result {
    std::vector<double> t, x, xm, u, theta_x, theta_r, error;
    Metrics metrics;
};

double reference(double t) {
    if (t < 4.0) return 1.0;
    if (t < 8.0) return -0.5;
    if (t < 13.0) return 0.8;
    return 0.2;
}

double quantize(double value, double quantum) {
    if (quantum <= 0.0) return value;
    return quantum * std::round(value / quantum);
}

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

template <typename F>
double rk4Scalar(double x, double dt, F rhs) {
    const double k1 = rhs(x);
    const double k2 = rhs(x + 0.5 * dt * k1);
    const double k3 = rhs(x + 0.5 * dt * k2);
    const double k4 = rhs(x + dt * k3);
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

Result simulate(const Config& cfg) {
    const double a = 1.2, b = 1.0;
    const double am = 2.0, bm = 2.0;
    const int steps = static_cast<int>(std::lround(cfg.final_time / cfg.sample_time));

    Result out;
    out.t.resize(steps + 1);
    out.x.assign(steps + 1, 0.0);
    out.xm.assign(steps + 1, 0.0);
    out.u.assign(steps + 1, 0.0);
    out.theta_x.assign(steps + 1, 0.0);
    out.theta_r.assign(steps + 1, 0.0);
    out.error.assign(steps + 1, 0.0);
    for (int k = 0; k <= steps; ++k) out.t[k] = k * cfg.sample_time;

    std::deque<std::tuple<double, double, double>> history;
    for (int i = 0; i <= cfg.delay_steps; ++i) {
        history.emplace_back(0.0, 0.0, reference(0.0));
    }

    int saturated_samples = 0;
    for (int k = 0; k < steps; ++k) {
        const double r_now = reference(out.t[k]);
        history.emplace_back(quantize(out.x[k], cfg.quantization), out.xm[k], r_now);
        history.pop_front();
        const auto [x_d, xm_d, r_d] = history.front();

        const double e_d = x_d - xm_d;
        const double normalizer = 1.0 + cfg.normalization * (x_d * x_d + r_d * r_d);
        const double u_raw = out.theta_x[k] * x_d + out.theta_r[k] * r_d;
        out.u[k] = std::clamp(u_raw, -cfg.u_limit, cfg.u_limit);
        const bool saturated = std::abs(u_raw - out.u[k]) > 1.0e-12;
        saturated_samples += saturated ? 1 : 0;

        double tx_next = out.theta_x[k];
        double tr_next = out.theta_r[k];
        if (!saturated) {
            tx_next -= cfg.sample_time * cfg.gamma_x * e_d * x_d / normalizer;
            tr_next -= cfg.sample_time * cfg.gamma_r * e_d * r_d / normalizer;
        }
        out.theta_x[k + 1] = project(tx_next, cfg.theta_limit);
        out.theta_r[k + 1] = project(tr_next, cfg.theta_limit);

        const double dt = cfg.sample_time / cfg.integration_substeps;
        double x_local = out.x[k];
        double xm_local = out.xm[k];
        for (int j = 0; j < cfg.integration_substeps; ++j) {
            const double tj = out.t[k] + j * dt;
            const double r_sub = reference(tj);
            x_local = rk4Scalar(x_local, dt, [&](double z) { return -a * z + b * out.u[k]; });
            xm_local = rk4Scalar(xm_local, dt, [&](double z) { return -am * z + bm * r_sub; });
        }
        out.x[k + 1] = x_local;
        out.xm[k + 1] = xm_local;

        if (!std::isfinite(out.x[k + 1]) || !std::isfinite(out.xm[k + 1]) ||
            !std::isfinite(out.theta_x[k + 1]) || !std::isfinite(out.theta_r[k + 1])) {
            throw std::runtime_error("Non-finite value at sample " + std::to_string(k));
        }
    }

    out.u.back() = out.u[steps - 1];
    double sum_sq = 0.0;
    double max_abs = 0.0;
    for (int k = 0; k <= steps; ++k) {
        out.error[k] = out.x[k] - out.xm[k];
        sum_sq += out.error[k] * out.error[k];
        max_abs = std::max(max_abs, std::abs(out.error[k]));
    }
    out.metrics.rms_error = std::sqrt(sum_sq / (steps + 1));
    out.metrics.max_error = max_abs;
    out.metrics.theta_x_final = out.theta_x.back();
    out.metrics.theta_r_final = out.theta_r.back();
    out.metrics.saturation_fraction = static_cast<double>(saturated_samples) / std::max(1, steps);
    return out;
}

void printResult(const Config& cfg, const Result& result) {
    const auto& m = result.metrics;
    std::cout << cfg.name << ": RMS=" << m.rms_error
              << ", max|e|=" << m.max_error
              << ", theta_x=" << m.theta_x_final
              << ", theta_r=" << m.theta_r_final
              << ", saturation=" << m.saturation_fraction << '\n';
}

int main() {
    try {
        const Config nominal{"nominal", 0.01, 0, 5};
        Config stressed{"stressed", 0.05, 2, 1};
        stressed.quantization = 0.002;

        const Result r1 = simulate(nominal);
        const Result r2 = simulate(stressed);
        std::cout << std::fixed << std::setprecision(6);
        printResult(nominal, r1);
        printResult(stressed, r2);

        std::cout << "\nCSV preview: t,x_nominal,x_stressed,xm\n";
        const std::size_t stride = std::max<std::size_t>(1, r1.t.size() / 10);
        for (std::size_t i = 0; i < r1.t.size(); i += stride) {
            const std::size_t j = std::min<std::size_t>(
                static_cast<std::size_t>(std::lround(r1.t[i] / stressed.sample_time)),
                r2.t.size() - 1);
            std::cout << r1.t[i] << ',' << r1.x[i] << ',' << r2.x[j] << ',' << r1.xm[i] << '\n';
        }
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter24_Lesson3.java

The Java version uses records for immutable configuration and result containers. A production real-time implementation must additionally control garbage creation, thread priority, timer behavior, and worst-case execution time.

// Chapter24_Lesson3.java
// Sampled-data MRAC with delayed measurements, normalization, projection,
// saturation-aware adaptation, and RK4 plant integration.
import java.util.ArrayDeque;
import java.util.Deque;
import java.util.Locale;

public final class Chapter24_Lesson3 {
    private record Config(
            String name,
            double sampleTime,
            int delaySteps,
            int integrationSubsteps,
            double gammaX,
            double gammaR,
            double normalization,
            double thetaLimit,
            double uLimit,
            double quantization,
            double finalTime) {
        static Config nominal() {
            return new Config("nominal", 0.01, 0, 5, 4.0, 4.0, 0.25, 8.0, 8.0, 0.0, 18.0);
        }

        static Config stressed() {
            return new Config("stressed", 0.05, 2, 1, 4.0, 4.0, 0.25, 8.0, 8.0, 0.002, 18.0);
        }
    }

    private record Sample(double x, double xm, double r) {}

    private record Metrics(
            double rmsError,
            double maxError,
            double thetaXFinal,
            double thetaRFinal,
            double saturationFraction) {}

    private record Result(
            double[] t,
            double[] x,
            double[] xm,
            double[] u,
            double[] thetaX,
            double[] thetaR,
            double[] error,
            Metrics metrics) {}

    @FunctionalInterface
    private interface ScalarRhs {
        double value(double x);
    }

    private static double reference(double t) {
        if (t < 4.0) return 1.0;
        if (t < 8.0) return -0.5;
        if (t < 13.0) return 0.8;
        return 0.2;
    }

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

    private static double quantize(double value, double quantum) {
        if (quantum <= 0.0) return value;
        return quantum * Math.rint(value / quantum);
    }

    private static double rk4Scalar(double x, double dt, ScalarRhs rhs) {
        double k1 = rhs.value(x);
        double k2 = rhs.value(x + 0.5 * dt * k1);
        double k3 = rhs.value(x + 0.5 * dt * k2);
        double k4 = rhs.value(x + dt * k3);
        return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    private static Result simulate(Config cfg) {
        final double a = 1.2, b = 1.0;
        final double am = 2.0, bm = 2.0;
        final int steps = (int) Math.round(cfg.finalTime() / cfg.sampleTime());

        double[] t = new double[steps + 1];
        double[] x = new double[steps + 1];
        double[] xm = new double[steps + 1];
        double[] u = new double[steps + 1];
        double[] thetaX = new double[steps + 1];
        double[] thetaR = new double[steps + 1];
        double[] error = new double[steps + 1];
        for (int k = 0; k <= steps; k++) t[k] = k * cfg.sampleTime();

        Deque<Sample> history = new ArrayDeque<>();
        for (int i = 0; i <= cfg.delaySteps(); i++) {
            history.addLast(new Sample(0.0, 0.0, reference(0.0)));
        }

        int saturatedSamples = 0;
        for (int k = 0; k < steps; k++) {
            double rNow = reference(t[k]);
            history.addLast(new Sample(quantize(x[k], cfg.quantization()), xm[k], rNow));
            history.removeFirst();
            Sample delayed = history.getFirst();

            double eDelayed = delayed.x() - delayed.xm();
            double normalizer = 1.0 + cfg.normalization()
                    * (delayed.x() * delayed.x() + delayed.r() * delayed.r());
            double uRaw = thetaX[k] * delayed.x() + thetaR[k] * delayed.r();
            u[k] = clamp(uRaw, -cfg.uLimit(), cfg.uLimit());
            boolean saturated = Math.abs(uRaw - u[k]) > 1.0e-12;
            if (saturated) saturatedSamples++;

            double txNext = thetaX[k];
            double trNext = thetaR[k];
            if (!saturated) {
                txNext -= cfg.sampleTime() * cfg.gammaX() * eDelayed * delayed.x() / normalizer;
                trNext -= cfg.sampleTime() * cfg.gammaR() * eDelayed * delayed.r() / normalizer;
            }
            thetaX[k + 1] = clamp(txNext, -cfg.thetaLimit(), cfg.thetaLimit());
            thetaR[k + 1] = clamp(trNext, -cfg.thetaLimit(), cfg.thetaLimit());

            double dt = cfg.sampleTime() / cfg.integrationSubsteps();
            double xLocal = x[k];
            double xmLocal = xm[k];
            for (int j = 0; j < cfg.integrationSubsteps(); j++) {
                double tj = t[k] + j * dt;
                double rSub = reference(tj);
                final double heldControl = u[k];
                xLocal = rk4Scalar(xLocal, dt, z -> -a * z + b * heldControl);
                final double heldReference = rSub;
                xmLocal = rk4Scalar(xmLocal, dt, z -> -am * z + bm * heldReference);
            }
            x[k + 1] = xLocal;
            xm[k + 1] = xmLocal;

            if (!Double.isFinite(x[k + 1]) || !Double.isFinite(xm[k + 1])
                    || !Double.isFinite(thetaX[k + 1]) || !Double.isFinite(thetaR[k + 1])) {
                throw new ArithmeticException("Non-finite value at sample " + k);
            }
        }

        u[steps] = u[steps - 1];
        double sumSquared = 0.0;
        double maxAbsolute = 0.0;
        for (int k = 0; k <= steps; k++) {
            error[k] = x[k] - xm[k];
            sumSquared += error[k] * error[k];
            maxAbsolute = Math.max(maxAbsolute, Math.abs(error[k]));
        }
        Metrics metrics = new Metrics(
                Math.sqrt(sumSquared / (steps + 1)),
                maxAbsolute,
                thetaX[steps],
                thetaR[steps],
                (double) saturatedSamples / Math.max(1, steps));
        return new Result(t, x, xm, u, thetaX, thetaR, error, metrics);
    }

    private static void printResult(Config cfg, Result result) {
        Metrics m = result.metrics();
        System.out.printf(Locale.US,
                "%s: RMS=%.6f, max|e|=%.6f, theta_x=%.6f, theta_r=%.6f, saturation=%.6f%n",
                cfg.name(), m.rmsError(), m.maxError(), m.thetaXFinal(),
                m.thetaRFinal(), m.saturationFraction());
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Config nominal = Config.nominal();
        Config stressed = Config.stressed();
        Result r1 = simulate(nominal);
        Result r2 = simulate(stressed);
        printResult(nominal, r1);
        printResult(stressed, r2);

        System.out.println("\nCSV preview: t,x_nominal,x_stressed,xm");
        int stride = Math.max(1, r1.t().length / 10);
        for (int i = 0; i < r1.t().length; i += stride) {
            int j = Math.min((int) Math.round(r1.t()[i] / stressed.sampleTime()), r2.t().length - 1);
            System.out.printf(Locale.US, "%.6f,%.6f,%.6f,%.6f%n",
                    r1.t()[i], r1.x()[i], r2.x()[j], r1.xm()[i]);
        }
    }
}

Chapter24_Lesson3.m

MATLAB uses local functions in one script. In Simulink, the same architecture can be assembled with a Zero-Order Hold, Integer Delay or Transport Delay, Discrete-Time Integrator for the parameters, Saturation, and Rate Transition blocks. Use a fixed-step solver for code-generation studies and ensure that the adaptive subsystem sample time is explicit rather than inherited accidentally.

% Chapter24_Lesson3.m
% Sampled-data MRAC with delay, projection, normalization, saturation-aware
% adaptation, and RK4 integration. This script compares nominal and stressed cases.
clear; clc; close all;

nominal = makeConfig("nominal", 0.01, 0, 5, 0.0);
stressed = makeConfig("stressed", 0.05, 2, 1, 0.002);

R1 = simulateMRAC(nominal);
R2 = simulateMRAC(stressed);

disp(nominal.name); disp(R1.metrics);
disp(stressed.name); disp(R2.metrics);

figure;
plot(R1.t, R1.x, 'DisplayName', 'x (nominal)'); hold on;
plot(R2.t, R2.x, 'DisplayName', 'x (stressed)');
plot(R1.t, R1.xm, '--', 'DisplayName', 'reference model');
xlabel('time [s]'); ylabel('state');
title('Sampled-data MRAC: sampling and delay sensitivity');
grid on; legend('Location', 'best');

figure;
plot(R1.t, R1.error, 'DisplayName', 'error (nominal)'); hold on;
plot(R2.t, R2.error, 'DisplayName', 'error (stressed)');
xlabel('time [s]'); ylabel('e = x - x_m');
title('Tracking-error comparison');
grid on; legend('Location', 'best');

function cfg = makeConfig(name, Ts, delaySteps, substeps, quantum)
    cfg.name = name;
    cfg.Ts = Ts;
    cfg.delaySteps = delaySteps;
    cfg.substeps = substeps;
    cfg.gammaX = 4.0;
    cfg.gammaR = 4.0;
    cfg.normalization = 0.25;
    cfg.thetaLimit = 8.0;
    cfg.uLimit = 8.0;
    cfg.quantization = quantum;
    cfg.finalTime = 18.0;
end

function r = commandReference(t)
    r = zeros(size(t));
    r(t < 4.0) = 1.0;
    r(t >= 4.0 & t < 8.0) = -0.5;
    r(t >= 8.0 & t < 13.0) = 0.8;
    r(t >= 13.0) = 0.2;
end

function value = quantizeSignal(value, quantum)
    if quantum > 0.0
        value = quantum * round(value / quantum);
    end
end

function xNext = rk4Scalar(x, dt, rhs)
    k1 = rhs(x);
    k2 = rhs(x + 0.5 * dt * k1);
    k3 = rhs(x + 0.5 * dt * k2);
    k4 = rhs(x + dt * k3);
    xNext = x + (dt / 6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
end

function out = simulateMRAC(cfg)
    a = 1.2; b = 1.0;
    am = 2.0; bm = 2.0;
    steps = round(cfg.finalTime / cfg.Ts);
    t = (0:steps)' * cfg.Ts;

    x = zeros(steps + 1, 1);
    xm = zeros(steps + 1, 1);
    u = zeros(steps + 1, 1);
    thetaX = zeros(steps + 1, 1);
    thetaR = zeros(steps + 1, 1);

    historyX = zeros(cfg.delaySteps + 1, 1);
    historyXm = zeros(cfg.delaySteps + 1, 1);
    historyR = commandReference(zeros(cfg.delaySteps + 1, 1));
    saturatedSamples = 0;

    for k = 1:steps
        rNow = commandReference(t(k));
        historyX = [historyX(2:end); quantizeSignal(x(k), cfg.quantization)]; %#ok<AGROW>
        historyXm = [historyXm(2:end); xm(k)]; %#ok<AGROW>
        historyR = [historyR(2:end); rNow]; %#ok<AGROW>
        xDelayed = historyX(1);
        xmDelayed = historyXm(1);
        rDelayed = historyR(1);

        eDelayed = xDelayed - xmDelayed;
        normalizer = 1.0 + cfg.normalization * (xDelayed^2 + rDelayed^2);
        uRaw = thetaX(k) * xDelayed + thetaR(k) * rDelayed;
        u(k) = min(cfg.uLimit, max(-cfg.uLimit, uRaw));
        saturated = abs(uRaw - u(k)) > 1.0e-12;
        saturatedSamples = saturatedSamples + saturated;

        if saturated
            txNext = thetaX(k);
            trNext = thetaR(k);
        else
            txNext = thetaX(k) - cfg.Ts * cfg.gammaX * eDelayed * xDelayed / normalizer;
            trNext = thetaR(k) - cfg.Ts * cfg.gammaR * eDelayed * rDelayed / normalizer;
        end
        thetaX(k + 1) = min(cfg.thetaLimit, max(-cfg.thetaLimit, txNext));
        thetaR(k + 1) = min(cfg.thetaLimit, max(-cfg.thetaLimit, trNext));

        dt = cfg.Ts / cfg.substeps;
        xLocal = x(k);
        xmLocal = xm(k);
        for j = 0:(cfg.substeps - 1)
            tj = t(k) + j * dt;
            rSub = commandReference(tj);
            heldU = u(k);
            xLocal = rk4Scalar(xLocal, dt, @(z) -a*z + b*heldU);
            xmLocal = rk4Scalar(xmLocal, dt, @(z) -am*z + bm*rSub);
        end
        x(k + 1) = xLocal;
        xm(k + 1) = xmLocal;

        if any(~isfinite([x(k + 1), xm(k + 1), thetaX(k + 1), thetaR(k + 1)]))
            error('Non-finite value at sample %d.', k);
        end
    end

    u(end) = u(end - 1);
    trackingError = x - xm;
    metrics.rmsError = sqrt(mean(trackingError.^2));
    metrics.maxError = max(abs(trackingError));
    metrics.thetaXFinal = thetaX(end);
    metrics.thetaRFinal = thetaR(end);
    metrics.saturationFraction = saturatedSamples / max(1, steps);

    out = struct('t', t, 'x', x, 'xm', xm, 'u', u, ...
        'thetaX', thetaX, 'thetaR', thetaR, 'error', trackingError, ...
        'metrics', metrics);
end

Chapter24_Lesson3.nb

This Wolfram Language source implements the same sampled loop with associations and lists. It can be pasted into a Mathematica notebook or opened as textual Wolfram Language source for evaluation.

(* Chapter24_Lesson3.nb
   Wolfram Language source that can be pasted into a Mathematica notebook.
   Sampled-data MRAC with delay, normalization, projection, saturation-aware
   adaptation, and RK4 integration. *)

ClearAll[reference, quantize, clip, rk4Scalar, simulateMRAC];

reference[t_?NumericQ] := Piecewise[{
   {1.0, t < 4.0},
   {-0.5, t < 8.0},
   {0.8, t < 13.0}}, 0.2];

quantize[value_?NumericQ, quantum_?NumericQ] :=
  If[quantum <= 0.0, value, quantum Round[value/quantum]];

clip[value_?NumericQ, limit_?NumericQ] := Clip[value, {-limit, limit}];

rk4Scalar[x_?NumericQ, dt_?NumericQ, rhs_] := Module[{k1, k2, k3, k4},
  k1 = rhs[x];
  k2 = rhs[x + 0.5 dt k1];
  k3 = rhs[x + 0.5 dt k2];
  k4 = rhs[x + dt k3];
  x + (dt/6.0) (k1 + 2.0 k2 + 2.0 k3 + k4)
];

simulateMRAC[cfg_Association] := Module[
  {a = 1.2, b = 1.0, am = 2.0, bm = 2.0, steps, t, x, xm, u,
   thetaX, thetaR, history, saturatedSamples = 0, k, rNow, delayed,
   xDelayed, xmDelayed, rDelayed, eDelayed, normalizer, uRaw,
   saturated, txNext, trNext, dt, xLocal, xmLocal, j, tj, rSub,
   trackingError, metrics},

  steps = Round[cfg["FinalTime"]/cfg["SampleTime"]];
  t = N[Range[0, steps] cfg["SampleTime"]];
  x = ConstantArray[0.0, steps + 1];
  xm = ConstantArray[0.0, steps + 1];
  u = ConstantArray[0.0, steps + 1];
  thetaX = ConstantArray[0.0, steps + 1];
  thetaR = ConstantArray[0.0, steps + 1];
  history = ConstantArray[{0.0, 0.0, reference[0.0]}, cfg["DelaySteps"] + 1];

  Do[
    rNow = reference[t[[k]]];
    history = Append[Rest[history],
      {quantize[x[[k]], cfg["Quantization"]], xm[[k]], rNow}];
    delayed = First[history];
    {xDelayed, xmDelayed, rDelayed} = delayed;
    eDelayed = xDelayed - xmDelayed;
    normalizer = 1.0 + cfg["Normalization"] (xDelayed^2 + rDelayed^2);

    uRaw = thetaX[[k]] xDelayed + thetaR[[k]] rDelayed;
    u[[k]] = clip[uRaw, cfg["ULimit"]];
    saturated = Abs[uRaw - u[[k]]] > 10^-12;
    If[saturated, saturatedSamples++];

    If[saturated,
      txNext = thetaX[[k]]; trNext = thetaR[[k]],
      txNext = thetaX[[k]] - cfg["SampleTime"] cfg["GammaX"] eDelayed xDelayed/normalizer;
      trNext = thetaR[[k]] - cfg["SampleTime"] cfg["GammaR"] eDelayed rDelayed/normalizer;
    ];
    thetaX[[k + 1]] = clip[txNext, cfg["ThetaLimit"]];
    thetaR[[k + 1]] = clip[trNext, cfg["ThetaLimit"]];

    dt = cfg["SampleTime"]/cfg["IntegrationSubsteps"];
    xLocal = x[[k]]; xmLocal = xm[[k]];
    Do[
      tj = t[[k]] + j dt;
      rSub = reference[tj];
      xLocal = rk4Scalar[xLocal, dt, Function[z, -a z + b u[[k]]]];
      xmLocal = rk4Scalar[xmLocal, dt, Function[z, -am z + bm rSub]],
      {j, 0, cfg["IntegrationSubsteps"] - 1}
    ];
    x[[k + 1]] = xLocal;
    xm[[k + 1]] = xmLocal;

    If[!And @@ (FiniteQ /@ {x[[k + 1]], xm[[k + 1]], thetaX[[k + 1]], thetaR[[k + 1]]}),
      Print["Non-finite value at sample ", k]; Abort[]],
    {k, 1, steps}
  ];

  u[[-1]] = u[[-2]];
  trackingError = x - xm;
  metrics = <|
    "RMSError" -> Sqrt[Mean[trackingError^2]],
    "MaxError" -> Max[Abs[trackingError]],
    "ThetaXFinal" -> Last[thetaX],
    "ThetaRFinal" -> Last[thetaR],
    "SaturationFraction" -> N[saturatedSamples/Max[1, steps]]
  |>;
  <|"t" -> t, "x" -> x, "xm" -> xm, "u" -> u,
    "thetaX" -> thetaX, "thetaR" -> thetaR,
    "error" -> trackingError, "metrics" -> metrics|>
];

nominal = <|"Name" -> "nominal", "SampleTime" -> 0.01,
  "DelaySteps" -> 0, "IntegrationSubsteps" -> 5,
  "GammaX" -> 4.0, "GammaR" -> 4.0, "Normalization" -> 0.25,
  "ThetaLimit" -> 8.0, "ULimit" -> 8.0, "Quantization" -> 0.0,
  "FinalTime" -> 18.0|>;

stressed = <|"Name" -> "stressed", "SampleTime" -> 0.05,
  "DelaySteps" -> 2, "IntegrationSubsteps" -> 1,
  "GammaX" -> 4.0, "GammaR" -> 4.0, "Normalization" -> 0.25,
  "ThetaLimit" -> 8.0, "ULimit" -> 8.0, "Quantization" -> 0.002,
  "FinalTime" -> 18.0|>;

r1 = simulateMRAC[nominal];
r2 = simulateMRAC[stressed];
Print["Nominal metrics: ", r1["metrics"]];
Print["Stressed metrics: ", r2["metrics"]];

ListLinePlot[
  {Transpose[{r1["t"], r1["x"]}],
   Transpose[{r2["t"], r2["x"]}],
   Transpose[{r1["t"], r1["xm"]}]},
  PlotLegends -> {"x (nominal)", "x (stressed)", "reference model"},
  AxesLabel -> {"time [s]", "state"},
  PlotLabel -> "Sampled-data MRAC: sampling and delay sensitivity",
  GridLines -> Automatic,
  ImageSize -> Large
]

ListLinePlot[
  {Transpose[{r1["t"], r1["error"]}],
   Transpose[{r2["t"], r2["error"]}]},
  PlotLegends -> {"error (nominal)", "error (stressed)"},
  AxesLabel -> {"time [s]", "e = x - xm"},
  PlotLabel -> "Tracking-error comparison",
  GridLines -> Automatic,
  ImageSize -> Large
]

flowchart TD
  A["Start from continuous adaptive design"] --> B["Derive sampled and held closed-loop model"]
  B --> C["Include sensor, computation, network, and actuator delays"]
  C --> D["Select sample time from bandwidth and stability analysis"]
  D --> E["Choose discrete adaptation, normalization, and projection"]
  E --> F["Run solver and step-size convergence tests"]
  F --> G["Inject jitter, quantization, saturation, and deadline misses"]
  G --> H["Software-in-the-loop and processor-in-the-loop tests"]
  H --> I["Hardware test with safety monitor and fallback controller"]
  I --> J["Record timing, parameters, limits, and fault events"]
        

The implementation should be considered acceptable only when conclusions are robust to smaller numerical steps, realistic timing variation, finite sensor resolution, command limits, parameter initialization, and bounded model mismatch. A successful nominal plot is not sufficient evidence.

Sampling and timing: Specify nominal, minimum, and maximum sample intervals; measure worst-case execution time; timestamp every sensor packet; define deadline-miss behavior; log sensor-to-actuator latency.

Adaptive update: Include the elapsed-time factor; normalize regressors; project parameters; limit per-sample increments; freeze or modify adaptation under saturation and invalid data.

Numerics: Separate control sampling from plant integration; perform step-size convergence tests; monitor NaN/Inf; preserve covariance symmetry and positive definiteness; scale states and regressors to comparable magnitudes.

Delay: Model integer and fractional delays; verify augmented closed-loop eigenvalues or a sampled-data Lyapunov condition; do not rely only on a continuous phase-margin calculation.

Safety: Enforce parameter and command bounds independently of the adaptive code; retain a fixed-gain fallback controller; monitor tracking error, parameter rate, timing overrun, saturation duty cycle, and excitation quality.

Reproducibility: Store software version, compiler options, numeric precision, controller period, solver settings, random seeds, hardware timestamps, and all active limits with each experiment.

Problem 1 (Euler Stability Limit): A reference-model state obeys \(\dot{x}_m=-8x_m+8r\). If forward Euler is used for the model, determine the interval of sample times that makes the homogeneous numerical model asymptotically stable. What happens at \(T_s=0.30\,\text{s}\)?

Solution: The homogeneous Euler recursion is

\[ x_{m,k+1}=(1-8T_s)x_{m,k}. \]

Stability requires

\[ |1-8T_s|<1 \quad\Longleftrightarrow\quad 0<T_s<0.25\,\text{s}. \]

At \(T_s=0.30\), the pole is \(1-8(0.30)=-1.4\), whose magnitude exceeds one. The numerical reference model diverges with alternating sign even though the continuous reference model is stable.

Problem 2 (One-Sample Delay): Consider \(x_{k+1}=0.8x_k+0.2u_k\) with \(u_k=-2x_{k-1}\). Form the augmented matrix and determine whether the delayed closed loop is asymptotically stable.

Solution: With \(\boldsymbol{\xi}_k=[x_k,x_{k-1}]^T\),

\[ \boldsymbol{\xi}_{k+1}= \begin{bmatrix} 0.8 & -0.4\\ 1 & 0 \end{bmatrix}\boldsymbol{\xi}_k. \]

The characteristic equation is

\[ \lambda^2-0.8\lambda+0.4=0. \]

The roots are a complex-conjugate pair whose product is 0.4; therefore each has magnitude \(\sqrt{0.4}\approx0.6325<1\). This delayed loop is stable. The calculation also shows why delay must be represented by augmented dynamics rather than ignored.

Problem 3 (Adaptation Gain and Sampling Rate): An implementation uses \(\hat{\theta}_{k+1}=\hat{\theta}_k-0.02\,\phi_ke_k\) at \(T_s=0.01\,\text{s}\). Estimate the equivalent continuous gradient gain. If the sample time is reduced to 0.002 s but the coefficient 0.02 is left unchanged, what equivalent gain results?

Solution: Comparing with

\[ \hat{\theta}_{k+1}=\hat{\theta}_k-T_s\gamma\phi_ke_k, \]

gives \(\gamma=0.02/0.01=2\). After changing the period without changing the coefficient, \(\gamma=0.02/0.002=10\). The effective adaptation becomes five times faster. Correct sample-time-independent code stores \(\gamma\) and explicitly multiplies by measured \(T_k\).

Problem 4 (Quantization Dead Zone): A position sensor has resolution \(q=0.004\) units, and other bounded measurement noise is at most 0.003 units. Propose a conservative error dead-zone threshold that prevents adaptation due solely to these effects.

Solution: Uniform quantization contributes at most \(q/2=0.002\). The worst-case combined magnitude is bounded by

\[ |n_k|\le0.002+0.003=0.005. \]

A threshold slightly above this bound, for example \(e_0=0.006\), prevents updates caused only by the modeled noise floor. The price is practical rather than exact tracking inside the dead zone.

Problem 5 (Jitter Bound): A sampled-data Lyapunov analysis yields

\[ \Delta V_k\le-3T_k\|z_k\|^2+40T_k^2\|z_k\|^2. \]

Find a sufficient upper bound on \(T_k\). If the nominal period is 50 ms, what positive jitter can be tolerated by this bound?

Solution: Require

\[ -3T_k+40T_k^2<0 \quad\Longleftrightarrow\quad 0<T_k<\frac{3}{40}=0.075\,\text{s}. \]

With \(T_s=0.050\,\text{s}\), the sufficient positive jitter margin is \(0.075-0.050=0.025\,\text{s}\), or 25 ms. This is a theoretical bound under the assumptions used to obtain the inequality, not an allowance to ignore execution-time verification.

Problem 6 (Solver-Convergence Test): A nonlinear adaptive simulation gives RMS tracking errors 0.180, 0.151, and 0.149 for integration steps 10 ms, 5 ms, and 2.5 ms, while the controller period remains 20 ms. Interpret the result.

Solution: The change from 10 ms to 5 ms is large, so the 10 ms integration is not converged. The change from 5 ms to 2.5 ms is small, suggesting convergence near 0.149, although another refinement is prudent. The controller period must remain 20 ms in all runs; otherwise the experiment confounds solver error with a different sampled controller.

Digital adaptive control is a sampled-data, delayed, finite-precision hybrid system. Exact zero-order-hold models reveal the difference between physical sampling and numerical approximation. Delays add state and alter characteristic roots. Euler parameter updates introduce second-order Lyapunov terms, making the sampling period and adaptation increment part of the stability problem. Jitter, quantization, saturation, solver error, stale data, and finite-precision covariance updates can all drive parameter drift or instability. Safe implementation therefore combines sampled analysis with normalization, projection, increment limits, explicit timing policies, solver-convergence studies, realistic fault injection, and an independent fallback controller.