In the preceding lessons, the plant, estimator, and control update were written at ideal sampling instants. A real implementation inserts an analog-to-digital converter, finite-word arithmetic, a task scheduler, an actuator interface, and a zero-order hold. The controller therefore operates on quantized measurements and its command reaches the plant only after a finite computation and communication time.

For a scalar plant, the implementation model used throughout this lesson is

\[ y_{k+1}=a y_k+b u_k+w_k,\qquad u_k=c_{k-d},\qquad y_k^{q}=Q_{\Delta_y}(y_k+v_k), \]

where \( c_k \) is the newly calculated command, \( d \) is an integer delay in samples, \( u_k \) is the command actually applied during the current interval, and \( Q_{\Delta_y} \) is a quantizer. The estimator must use the applied input \( u_k \), not merely the command \( c_k \); otherwise the regressor is time-misaligned.

flowchart TD
  R["Reference r(k)"] --> C["Adaptive controller"]
  YQ["Quantized measurement yq(k)"] --> C
  C --> CQ["Command quantizer and saturation"]
  CQ --> D["Computation and transport delay"]
  D --> H["Zero-order hold"]
  H --> P["Physical plant"]
  P --> S["Sensor and sampler"]
  S --> YQ
  YQ --> E["Online estimator"]
  H --> E
  E --> C
        

The three implementation effects are coupled. A larger sampling period changes the discrete plant coefficients, converts a fixed computation time into a larger or smaller fractional delay, and changes the physical meaning of one quantization count per sample. Consequently, sampling, quantization, and delay should be analyzed as one closed-loop design problem rather than as independent software details.

Consider the continuous-time linear plant \( \dot{x}(t)=A_c x(t)+B_c u(t) \). If the digital controller holds \( u(t)=u_k \) for \( t\in[kT_s,(k+1)T_s) \), exact integration gives

\[ x_{k+1}=A_d x_k+B_d u_k, \qquad A_d=e^{A_cT_s}, \qquad B_d=\int_0^{T_s}e^{A_c\tau}B_c\,d\tau. \]

This is not an approximation: it is the exact sample-to-sample model under a zero-order hold. For the first-order process

\[ \tau_p\dot{y}(t)+y(t)=K_pu(t), \]

the exact sampled coefficients are

\[ a(T_s)=e^{-T_s/\tau_p}, \qquad b(T_s)=K_p\left(1-e^{-T_s/\tau_p}\right). \]

Hence changing \( T_s \) changes the parameter vector being estimated. If a controller changes its sampling rate online but retains old estimates of \( a \) and \( b \), the estimator experiences an artificial parameter jump even though the physical plant has not changed.

2.1 Euler approximation and its error

The forward-Euler approximation is \( A_E=I+A_cT_s \) and \( B_E=B_cT_s \). Expanding the matrix exponential,

\[ e^{A_cT_s} =I+A_cT_s+\frac{A_c^2T_s^2}{2!} +\frac{A_c^3T_s^3}{3!}+\cdots . \]

Therefore the local matrix error is second order in the sampling period. Using a compatible matrix norm,

\[ \left\|A_d-A_E\right\| \le e^{\|A_c\|T_s}-1-\|A_c\|T_s =\mathcal{O}(T_s^2). \]

Adaptive control can magnify this modeling error because the estimator may interpret discretization error as a parameter change. Exact zero-order-hold discretization should therefore be used whenever the continuous model is known and the matrix exponential is available offline.

The sampling period must satisfy several competing requirements: adequate closed-loop phase margin, sufficiently small intersample motion, enough processor time to finish the adaptive update, and enough excitation per observation for reliable estimation. Sampling faster is not automatically better: very small changes in output may become comparable with sensor quantization, while covariance updates occur more frequently and can amplify numerical roundoff.

3.1 Bandwidth and delay-phase constraint

A pure delay of \( dT_s \) contributes the frequency response factor \( e^{-j\omega dT_s} \) and phase lag

\[ \phi_d(\omega)=-\omega dT_s. \]

If the delay is allowed to consume at most \( \phi_{\mathrm{allow} } \) radians at the desired closed-loop bandwidth \( \omega_b \), then

\[ T_s\le \frac{\phi_{\mathrm{allow} } }{d\omega_b}. \]

This relation is a delay budget, not a complete stability theorem. It is useful for rejecting sampling choices that necessarily consume too much phase margin before the adaptive dynamics are considered.

3.2 Sampling jitter

Let the actual interval be \( T_k=T_s+\delta T_k \), with bounded jitter \( |\delta T_k|\le J \). For the first-order process,

\[ a_k=e^{-T_k/\tau_p}, \qquad |a_k-a(T_s)| \le \frac{J}{\tau_p} e^{-(T_s-J)/\tau_p}, \qquad J<T_s. \]

The inequality follows from the mean-value theorem because the maximum magnitude of \( da/dT \) on \( [T_s-J,T_s+J] \) occurs at the smallest interval. Thus scheduler jitter appears to the estimator as bounded time-varying parameter uncertainty.

3.3 Deadline condition

A deterministic implementation should satisfy

\[ t_{\mathrm{ADC} }+ t_{\mathrm{filter} }+ t_{\mathrm{estimator} }+ t_{\mathrm{controller} }+ t_{\mathrm{DAC} } \le T_s-t_{\mathrm{margin} }. \]

If this condition is occasionally violated, the effective delay switches between different integers. Such deadline misses are more damaging than a known constant delay because the controller and regressor alignment become time varying.

A uniform mid-tread quantizer with resolution \( \Delta \) is

\[ Q_\Delta(z) =\Delta\left\lfloor\frac{z}{\Delta}+\frac{1}{2}\right\rfloor. \]

When the quantizer is not saturated, define \( \eta_q=Q_\Delta(z)-z \). Then

\[ |\eta_q|\le\frac{\Delta}{2}. \]

This deterministic bound is preferable to automatically modeling quantization as white noise. The white-noise approximation requires additional assumptions, such as sufficiently rich signals and weak correlation between the signal and the quantization error. Near a setpoint, those assumptions often fail and deterministic limit cycles may appear.

4.1 Quantized regression equation

Let both output and applied input be quantized:

\[ y_k^q=y_k+\eta_k^y, \qquad u_k^q=u_k+\eta_k^u. \]

Substitution into the scalar plant gives the measured regression

\[ y_{k+1}^q = \begin{bmatrix}y_k^q & u_k^q\end{bmatrix} \begin{bmatrix}a\\b\end{bmatrix} +w_k+d_k^q, \]

where the quantization-induced disturbance is

\[ d_k^q = \eta_{k+1}^y-a\eta_k^y-b\eta_k^u. \]

Hence, before saturation,

\[ |d_k^q| \le \frac{1+|a|}{2}\Delta_y +\frac{|b|}{2}\Delta_u. \]

This bound explains why reducing only sensor quantization may not improve identification if the actuator command remains coarse. Both channels enter the prediction error.

4.2 Finite range and saturation

A real converter has both resolution and range. If \( |z| \) exceeds the representable range, the quantization error is no longer bounded by \( \Delta/2 \). Saturation must therefore be treated as a separate nonlinearity. Parameter projection, command saturation, and anti-drift logic should be active before hardware limits are reached.

Suppose command \( c_k \) is applied after \( d \) complete sample intervals. The plant equation can be written as

\[ y_{k+1}=a y_k+b c_{k-d}+w_k. \]

A delay can be represented by an augmented state containing the command queue:

\[ \xi_k= \begin{bmatrix} y_k & c_{k-d} & c_{k-d+1} & \cdots & c_{k-1} \end{bmatrix}^T. \]

The augmented model is delay free but has higher order. For the simple certainty-equivalent controller in this lesson, it is more transparent to predict through the queue directly. Let \( q_i(k) \) denote the inputs already scheduled over the delay horizon, with \( q_0(k)=u_k \) equal to the input applied in the current interval. Starting from the quantized measurement,

\[ \widehat y_{k|k}=y_k^q, \qquad \widehat y_{k+i+1|k} = \widehat a_k\widehat y_{k+i|k} +\widehat b_k q_i(k), \quad i=0,\ldots,d-1. \]

The newly calculated command affects \( y_{k+d+1} \). With desired scalar closed-loop pole \( a_c \), choose the target

\[ y_{k+d+1}^\star = a_c\widehat y_{k+d|k} +(1-a_c)r_k, \qquad |a_c|<1, \]

and compute

\[ c_k= \operatorname{sat}_{u_{\max} } \left( \frac{ y_{k+d+1}^\star- \widehat a_k\widehat y_{k+d|k} }{\widehat b_k} \right). \]

Projection must keep \( \widehat b_k \) away from zero; otherwise the certainty-equivalent division can produce an unbounded command. The implementation example assumes the sign of \( b \) is known, consistent with the sign-definite case studied earlier in the course.

Define the regressor from the quantized output and the input actually applied during the interval:

\[ \varphi_k= \begin{bmatrix}y_k^q\\u_k\end{bmatrix}, \qquad \widehat y_{k+1|k}=\widehat\theta_k^T\varphi_k, \qquad \widehat\theta_k= \begin{bmatrix}\widehat a_k\\\widehat b_k\end{bmatrix}. \]

The exponentially weighted RLS recursion is

\[ L_k= \frac{P_k\varphi_k} {\lambda+\varphi_k^TP_k\varphi_k}, \qquad \varepsilon_{k+1} =y_{k+1}^q-\widehat\theta_k^T\varphi_k, \]

\[ \widehat\theta_{k+1} = \operatorname{Proj}_\Theta \left( \widehat\theta_k+L_k\varepsilon_{k+1} \right), \]

\[ P_{k+1} = \frac{1}{\lambda} \left( P_k-L_k\varphi_k^TP_k \right), \qquad 0<\lambda\le1. \]

In exact arithmetic, the covariance matrix remains symmetric. In finite precision it should be explicitly symmetrized:

\[ P_{k+1} ← \frac{1}{2} \left(P_{k+1}+P_{k+1}^T\right) +\epsilon_P I, \qquad \epsilon_P>0. \]

The small diagonal regularizer is a numerical safeguard, not a substitute for excitation. Quantization may cause \( \varphi_k \) to remain constant over many samples, so the estimator can temporarily receive almost no new information.

Under ideal parameters and exact delay prediction, the scalar tracking error satisfies \( e_{k+1}=a_c e_k \). Sampling mismatch, parameter error, quantization, disturbances, and imperfect delay prediction combine into a perturbation \( \zeta_k \):

\[ e_{k+1}=a_c e_k+\zeta_k, \qquad |a_c|<1, \qquad |\zeta_k|\le\overline\zeta. \]

7.1 Direct geometric bound

Repeated substitution yields

\[ e_k=a_c^k e_0+ \sum_{i=0}^{k-1}a_c^{k-1-i}\zeta_i. \]

Taking absolute values and using the geometric series,

\[ |e_k| \le |a_c|^k|e_0| + \frac{1-|a_c|^k}{1-|a_c|}\overline\zeta. \]

Consequently,

\[ \limsup_{k→\infty}|e_k| \le \frac{\overline\zeta}{1-|a_c|}. \]

Thus fixed quantization and bounded implementation errors generally imply practical convergence to a neighborhood, not exact asymptotic tracking. The neighborhood shrinks when converter resolution improves and when the desired pole is not chosen too close to one.

7.2 Lyapunov-difference proof

Let \( V_k=e_k^2 \). Then

\[ \Delta V_k = (a_c^2-1)e_k^2 +2a_c e_k\zeta_k+\zeta_k^2. \]

Young's inequality gives, for \( 0<\epsilon<1-a_c^2 \),

\[ 2|a_c||e_k||\zeta_k| \le \epsilon e_k^2+ \frac{a_c^2}{\epsilon}\zeta_k^2. \]

Therefore

\[ \Delta V_k \le -\alpha e_k^2+\beta\overline\zeta^2, \]

with

\[ \alpha=1-a_c^2-\epsilon>0, \qquad \beta=1+\frac{a_c^2}{\epsilon}. \]

Outside the set \( e_k^2\le(\beta/\alpha)\overline\zeta^2 \), the Lyapunov difference is negative. This establishes uniform ultimate boundedness of the tracking error under the stated bounded-perturbation model.

Converter quantization is only one source of finite precision. Fixed-point or floating-point arithmetic also affects the parameter update, covariance matrix, predictor, and command. For a signed fixed-point representation with \( F \) fractional bits, the numerical step is

\[ \Delta_{\mathrm{fxp} }=2^{-F}. \]

The word length must cover worst-case transients in \( P_k \), not only nominal plant signals. The following protections are recommended:

Covariance symmetry: replace \( P \) by \( (P+P^T)/2 \) after each update.

Positive diagonal floor: add a very small \( \epsilon_PI \) in floating-point simulation, or use a square-root RLS implementation when numerical conditioning is critical.

Parameter projection: enforce a physically meaningful compact set \( \Theta \), including a lower bound on \( |\widehat b| \).

Innovation gating: suspend or reduce adaptation when the regressor energy is below a threshold or when the actuator is saturated for a prolonged interval. Otherwise the estimator may learn the saturation nonlinearity rather than the plant.

Use of applied input: read back the held actuator command whenever possible. A command calculated by software is not necessarily the command delivered by a rate limiter, safety interlock, or lower-level drive.

A deterministic sample should preserve a fixed order so that the model indices are unambiguous. The estimator update uses the new measurement and the input applied over the preceding interval. The new control command is then calculated, quantized, queued, and committed before the deadline.

flowchart TD
  A["Timer interrupt at sample k"] --> B["Latch sensor and actuator readback"]
  B --> C["Quantize or convert to engineering units"]
  C --> D["Update estimator with y(k) and applied u(k-1)"]
  D --> E["Project parameters and condition covariance"]
  E --> F["Predict through queued delayed commands"]
  F --> G["Compute command c(k)"]
  G --> H["Apply saturation and actuator quantization"]
  H --> I["Append command to delay queue"]
  I --> J["Log timing, estimates, saturation, and deadline status"]
  J --> K["Commit output before next sample"]
        

The software log should include the nominal and actual sampling interval, task execution time, deadline misses, raw and quantized measurements, commanded and applied inputs, parameter estimates, covariance diagnostics, and saturation flags. Without these signals, a delay or quantization fault can be misdiagnosed as poor adaptive-law tuning.

This implementation uses numpy for matrix operations and optionally matplotlib for plots. It performs exact zero-order-hold discretization of a first-order continuous plant, models sensor and actuator quantization, maintains an integer command queue, and identifies the sampled coefficients by RLS.

Chapter14_Lesson4.py

# Chapter14_Lesson4.py
"""
Discrete-time adaptive control with sampling, quantization, and integer-sample delay.

Dependencies:
    numpy
    matplotlib (optional, only for plotting)

The estimator identifies y[k+1] = a*y[k] + b*u_applied[k].
The controller predicts through the queued delayed commands and computes a
certainty-equivalent command with projection, saturation, and quantization.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import List, Tuple
import math
import numpy as np


@dataclass
class Config:
    sample_time: float = 0.05
    duration: float = 30.0
    plant_time_constant: float = 0.8
    plant_dc_gain: float = 1.25
    delay_samples: int = 2
    sensor_step: float = 0.01
    actuator_step: float = 0.01
    actuator_limit: float = 3.0
    forgetting_factor: float = 0.995
    desired_pole: float = 0.65
    seed: int = 7


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


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


def project_parameters(theta: np.ndarray) -> np.ndarray:
    # Stable first-order model and known-positive input direction.
    projected = theta.copy()
    projected[0] = np.clip(projected[0], -0.98, 0.999)
    projected[1] = np.clip(projected[1], 0.05, 2.50)
    return projected


def reference_signal(time_value: float) -> float:
    if time_value < 3.0:
        return 0.0
    if time_value < 12.0:
        return 1.0
    if time_value < 20.0:
        return -0.6
    return 0.8


def simulate(config: Config) -> Tuple[np.ndarray, ...]:
    rng = np.random.default_rng(config.seed)
    steps = int(round(config.duration / config.sample_time))

    # Exact zero-order-hold discretization of K/(tau*s + 1).
    a_true = math.exp(-config.sample_time / config.plant_time_constant)
    b_true = config.plant_dc_gain * (1.0 - a_true)

    time = np.arange(steps + 1, dtype=float) * config.sample_time
    y = np.zeros(steps + 1)
    y_measured = np.zeros(steps + 1)
    reference = np.array([reference_signal(t) for t in time])
    applied_input = np.zeros(steps)
    command = np.zeros(steps)
    theta_history = np.zeros((steps + 1, 2))

    theta = np.array([0.70, 0.30], dtype=float)
    covariance = 100.0 * np.eye(2)
    theta_history[0] = theta

    # queue[0] is the command applied now; a new command is appended at the end.
    queue: List[float] = [0.0] * max(config.delay_samples, 1)

    y_measured[0] = quantize(y[0], config.sensor_step)

    for k in range(steps):
        u_now = queue.pop(0) if config.delay_samples > 0 else 0.0
        applied_input[k] = u_now

        # Predict the state to the instant when the new command will be applied.
        predicted = y_measured[k]
        if config.delay_samples > 0:
            predicted = theta[0] * predicted + theta[1] * u_now
            for queued_input in queue:
                predicted = theta[0] * predicted + theta[1] * queued_input

        target_next = (
            config.desired_pole * predicted
            + (1.0 - config.desired_pole) * reference[k]
        )
        raw_command = (target_next - theta[0] * predicted) / theta[1]
        new_command = quantize(
            saturate(raw_command, config.actuator_limit),
            config.actuator_step,
        )
        command[k] = new_command

        if config.delay_samples > 0:
            queue.append(new_command)
        else:
            # Zero-delay implementation: the command is applied in this interval.
            applied_input[k] = new_command
            u_now = new_command

        disturbance = 0.004 * rng.standard_normal()
        y[k + 1] = a_true * y[k] + b_true * u_now + disturbance
        y_measured[k + 1] = quantize(y[k + 1], config.sensor_step)

        # RLS uses the input that actually acted on the plant, not the command.
        phi = np.array([y_measured[k], u_now])
        prediction_error = y_measured[k + 1] - float(phi @ theta)
        denominator = config.forgetting_factor + float(phi @ covariance @ phi)
        gain = (covariance @ phi) / denominator

        theta = project_parameters(theta + gain * prediction_error)
        covariance = (
            covariance - np.outer(gain, phi) @ covariance
        ) / config.forgetting_factor
        covariance = 0.5 * (covariance + covariance.T)
        covariance += 1.0e-10 * np.eye(2)
        theta_history[k + 1] = theta

    return (
        time,
        reference,
        y,
        y_measured,
        applied_input,
        command,
        theta_history,
        np.array([a_true, b_true]),
    )


def main() -> None:
    config = Config()
    (
        time,
        reference,
        y,
        y_measured,
        applied_input,
        command,
        theta_history,
        true_parameters,
    ) = simulate(config)

    start = int(5.0 / config.sample_time)
    rmse = math.sqrt(np.mean((reference[start:] - y[start:]) ** 2))

    print(f"Sampling period: {config.sample_time:.4f} s")
    print(f"Delay: {config.delay_samples} samples")
    print(f"Tracking RMSE after 5 s: {rmse:.6f}")
    print(f"True [a, b]: {true_parameters}")
    print(f"Estimated [a, b]: {theta_history[-1]}")

    data = np.column_stack(
        (
            time,
            reference,
            y,
            y_measured,
            np.r_[applied_input, np.nan],
            np.r_[command, np.nan],
            theta_history[:, 0],
            theta_history[:, 1],
        )
    )
    np.savetxt(
        "Chapter14_Lesson4_results.csv",
        data,
        delimiter=",",
        header="time,reference,y,y_quantized,u_applied,u_command,a_hat,b_hat",
        comments="",
    )

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(time, reference, "--", label="reference")
        plt.plot(time, y, label="plant output")
        plt.step(time, y_measured, where="post", label="quantized measurement")
        plt.xlabel("time (s)")
        plt.ylabel("output")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.savefig("Chapter14_Lesson4_tracking.png", dpi=180)

        plt.figure()
        plt.plot(time, theta_history[:, 0], label="a_hat")
        plt.plot(time, theta_history[:, 1], label="b_hat")
        plt.axhline(true_parameters[0], linestyle="--", label="a true")
        plt.axhline(true_parameters[1], linestyle="--", label="b true")
        plt.xlabel("time (s)")
        plt.ylabel("parameter")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.savefig("Chapter14_Lesson4_parameters.png", dpi=180)
        plt.show()
    except ImportError:
        print("matplotlib is not installed; CSV output was still generated.")


if __name__ == "__main__":
    main()
      

The C++ version uses only the standard library. Fixed-size arrays are used for the two-parameter RLS update, while std::deque represents the delayed command queue. The code writes a CSV file suitable for plotting or hardware-in-the-loop comparison.

Chapter14_Lesson4.cpp

// Chapter14_Lesson4.cpp
// Discrete-time adaptive control with sampling, quantization, and delay.
// Build: g++ -std=c++17 -O2 Chapter14_Lesson4.cpp -o Chapter14_Lesson4

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

struct Config {
    double sampleTime = 0.05;
    double duration = 30.0;
    double plantTimeConstant = 0.8;
    double plantDcGain = 1.25;
    int delaySamples = 2;
    double sensorStep = 0.01;
    double actuatorStep = 0.01;
    double actuatorLimit = 3.0;
    double forgettingFactor = 0.995;
    double desiredPole = 0.65;
    unsigned int seed = 7;
};

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

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

double referenceSignal(double timeValue) {
    if (timeValue < 3.0) return 0.0;
    if (timeValue < 12.0) return 1.0;
    if (timeValue < 20.0) return -0.6;
    return 0.8;
}

std::array<double, 2> projectParameters(std::array<double, 2> theta) {
    theta[0] = std::clamp(theta[0], -0.98, 0.999);
    theta[1] = std::clamp(theta[1], 0.05, 2.50);
    return theta;
}

int main() {
    const Config cfg;
    if (cfg.delaySamples < 0) {
        throw std::invalid_argument("delaySamples must be nonnegative");
    }

    const int steps = static_cast<int>(std::llround(cfg.duration / cfg.sampleTime));
    const double aTrue = std::exp(-cfg.sampleTime / cfg.plantTimeConstant);
    const double bTrue = cfg.plantDcGain * (1.0 - aTrue);

    std::vector<double> time(steps + 1, 0.0);
    std::vector<double> reference(steps + 1, 0.0);
    std::vector<double> y(steps + 1, 0.0);
    std::vector<double> yMeasured(steps + 1, 0.0);
    std::vector<double> uApplied(steps, 0.0);
    std::vector<double> command(steps, 0.0);
    std::vector<std::array<double, 2>> thetaHistory(steps + 1);

    for (int k = 0; k <= steps; ++k) {
        time[k] = k * cfg.sampleTime;
        reference[k] = referenceSignal(time[k]);
    }

    std::array<double, 2> theta{0.70, 0.30};
    std::array<std::array<double, 2>, 2> P{ {
        { {100.0, 0.0} },
        { {0.0, 100.0} }
    } };
    thetaHistory[0] = theta;
    yMeasured[0] = quantize(y[0], cfg.sensorStep);

    std::deque<double> queue;
    for (int i = 0; i < cfg.delaySamples; ++i) {
        queue.push_back(0.0);
    }

    std::mt19937 generator(cfg.seed);
    std::normal_distribution<double> normal(0.0, 1.0);

    for (int k = 0; k < steps; ++k) {
        double uNow = 0.0;
        if (cfg.delaySamples > 0) {
            uNow = queue.front();
            queue.pop_front();
        }

        double predicted = yMeasured[k];
        if (cfg.delaySamples > 0) {
            predicted = theta[0] * predicted + theta[1] * uNow;
            for (double queuedInput : queue) {
                predicted = theta[0] * predicted + theta[1] * queuedInput;
            }
        }

        const double targetNext =
            cfg.desiredPole * predicted
            + (1.0 - cfg.desiredPole) * reference[k];
        const double rawCommand =
            (targetNext - theta[0] * predicted) / theta[1];
        const double newCommand = quantize(
            saturate(rawCommand, cfg.actuatorLimit),
            cfg.actuatorStep
        );
        command[k] = newCommand;

        if (cfg.delaySamples > 0) {
            queue.push_back(newCommand);
        } else {
            uNow = newCommand;
        }
        uApplied[k] = uNow;

        const double disturbance = 0.004 * normal(generator);
        y[k + 1] = aTrue * y[k] + bTrue * uNow + disturbance;
        yMeasured[k + 1] = quantize(y[k + 1], cfg.sensorStep);

        const std::array<double, 2> phi{yMeasured[k], uNow};
        const double predictedOutput = phi[0] * theta[0] + phi[1] * theta[1];
        const double error = yMeasured[k + 1] - predictedOutput;

        const std::array<double, 2> Pphi{
            P[0][0] * phi[0] + P[0][1] * phi[1],
            P[1][0] * phi[0] + P[1][1] * phi[1]
        };
        const double denominator =
            cfg.forgettingFactor + phi[0] * Pphi[0] + phi[1] * Pphi[1];
        const std::array<double, 2> gain{
            Pphi[0] / denominator,
            Pphi[1] / denominator
        };

        theta[0] += gain[0] * error;
        theta[1] += gain[1] * error;
        theta = projectParameters(theta);

        std::array<std::array<double, 2>, 2> updatedP{};
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                const double phiTPj = phi[0] * P[0][j] + phi[1] * P[1][j];
                updatedP[i][j] =
                    (P[i][j] - gain[i] * phiTPj) / cfg.forgettingFactor;
            }
        }
        const double offDiagonal =
            0.5 * (updatedP[0][1] + updatedP[1][0]);
        updatedP[0][1] = offDiagonal;
        updatedP[1][0] = offDiagonal;
        updatedP[0][0] += 1.0e-10;
        updatedP[1][1] += 1.0e-10;
        P = updatedP;
        thetaHistory[k + 1] = theta;
    }

    const int start = static_cast<int>(5.0 / cfg.sampleTime);
    double squaredErrorSum = 0.0;
    int count = 0;
    for (int k = start; k <= steps; ++k) {
        const double error = reference[k] - y[k];
        squaredErrorSum += error * error;
        ++count;
    }
    const double rmse = std::sqrt(squaredErrorSum / count);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Sampling period: " << cfg.sampleTime << " s\n";
    std::cout << "Delay: " << cfg.delaySamples << " samples\n";
    std::cout << "Tracking RMSE after 5 s: " << rmse << "\n";
    std::cout << "True [a, b]: [" << aTrue << ", " << bTrue << "]\n";
    std::cout << "Estimated [a, b]: ["
              << thetaHistory.back()[0] << ", "
              << thetaHistory.back()[1] << "]\n";

    std::ofstream csv("Chapter14_Lesson4_results_cpp.csv");
    csv << "time,reference,y,y_quantized,u_applied,u_command,a_hat,b_hat\n";
    csv << std::setprecision(12);
    for (int k = 0; k <= steps; ++k) {
        csv << time[k] << ',' << reference[k] << ',' << y[k] << ','
            << yMeasured[k] << ',';
        if (k < steps) {
            csv << uApplied[k] << ',' << command[k] << ',';
        } else {
            csv << ",,";
        }
        csv << thetaHistory[k][0] << ',' << thetaHistory[k][1] << '\n';
    }

    return 0;
}
      

The Java version uses Java SE classes only. It mirrors the C++ timing and indexing so numerical differences arise primarily from the random-number generator rather than from a different control algorithm.

Chapter14_Lesson4.java

// Chapter14_Lesson4.java
// Discrete-time adaptive control with sampling, quantization, and delay.
// Build and run:
//   javac Chapter14_Lesson4.java
//   java Chapter14_Lesson4

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayDeque;
import java.util.ArrayList;
import java.util.Deque;
import java.util.List;
import java.util.Locale;
import java.util.Random;

public final class Chapter14_Lesson4 {
    private static final class Config {
        final double sampleTime = 0.05;
        final double duration = 30.0;
        final double plantTimeConstant = 0.8;
        final double plantDcGain = 1.25;
        final int delaySamples = 2;
        final double sensorStep = 0.01;
        final double actuatorStep = 0.01;
        final double actuatorLimit = 3.0;
        final double forgettingFactor = 0.995;
        final double desiredPole = 0.65;
        final long seed = 7L;
    }

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

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

    private static double referenceSignal(double timeValue) {
        if (timeValue < 3.0) {
            return 0.0;
        }
        if (timeValue < 12.0) {
            return 1.0;
        }
        if (timeValue < 20.0) {
            return -0.6;
        }
        return 0.8;
    }

    private static void projectParameters(double[] theta) {
        theta[0] = Math.max(-0.98, Math.min(0.999, theta[0]));
        theta[1] = Math.max(0.05, Math.min(2.50, theta[1]));
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        Config cfg = new Config();
        if (cfg.delaySamples < 0) {
            throw new IllegalArgumentException("delaySamples must be nonnegative");
        }

        int steps = (int) Math.round(cfg.duration / cfg.sampleTime);
        double aTrue = Math.exp(-cfg.sampleTime / cfg.plantTimeConstant);
        double bTrue = cfg.plantDcGain * (1.0 - aTrue);

        double[] time = new double[steps + 1];
        double[] reference = new double[steps + 1];
        double[] y = new double[steps + 1];
        double[] yMeasured = new double[steps + 1];
        double[] uApplied = new double[steps];
        double[] command = new double[steps];
        double[][] thetaHistory = new double[steps + 1][2];

        for (int k = 0; k <= steps; ++k) {
            time[k] = k * cfg.sampleTime;
            reference[k] = referenceSignal(time[k]);
        }

        double[] theta = {0.70, 0.30};
        double[][] covariance = {
            {100.0, 0.0},
            {0.0, 100.0}
        };
        thetaHistory[0][0] = theta[0];
        thetaHistory[0][1] = theta[1];
        yMeasured[0] = quantize(y[0], cfg.sensorStep);

        Deque<Double> queue = new ArrayDeque<>();
        for (int i = 0; i < cfg.delaySamples; ++i) {
            queue.addLast(0.0);
        }

        Random random = new Random(cfg.seed);

        for (int k = 0; k < steps; ++k) {
            double uNow = 0.0;
            if (cfg.delaySamples > 0) {
                uNow = queue.removeFirst();
            }

            double predicted = yMeasured[k];
            if (cfg.delaySamples > 0) {
                predicted = theta[0] * predicted + theta[1] * uNow;
                List<Double> queuedValues = new ArrayList<>(queue);
                for (double queuedInput : queuedValues) {
                    predicted = theta[0] * predicted + theta[1] * queuedInput;
                }
            }

            double targetNext =
                cfg.desiredPole * predicted
                + (1.0 - cfg.desiredPole) * reference[k];
            double rawCommand =
                (targetNext - theta[0] * predicted) / theta[1];
            double newCommand = quantize(
                saturate(rawCommand, cfg.actuatorLimit),
                cfg.actuatorStep
            );
            command[k] = newCommand;

            if (cfg.delaySamples > 0) {
                queue.addLast(newCommand);
            } else {
                uNow = newCommand;
            }
            uApplied[k] = uNow;

            double disturbance = 0.004 * random.nextGaussian();
            y[k + 1] = aTrue * y[k] + bTrue * uNow + disturbance;
            yMeasured[k + 1] = quantize(y[k + 1], cfg.sensorStep);

            double[] phi = {yMeasured[k], uNow};
            double prediction = phi[0] * theta[0] + phi[1] * theta[1];
            double error = yMeasured[k + 1] - prediction;

            double[] pPhi = {
                covariance[0][0] * phi[0] + covariance[0][1] * phi[1],
                covariance[1][0] * phi[0] + covariance[1][1] * phi[1]
            };
            double denominator =
                cfg.forgettingFactor + phi[0] * pPhi[0] + phi[1] * pPhi[1];
            double[] gain = {
                pPhi[0] / denominator,
                pPhi[1] / denominator
            };

            theta[0] += gain[0] * error;
            theta[1] += gain[1] * error;
            projectParameters(theta);

            double[][] updated = new double[2][2];
            for (int i = 0; i < 2; ++i) {
                for (int j = 0; j < 2; ++j) {
                    double phiTPj =
                        phi[0] * covariance[0][j]
                        + phi[1] * covariance[1][j];
                    updated[i][j] =
                        (covariance[i][j] - gain[i] * phiTPj)
                        / cfg.forgettingFactor;
                }
            }
            double offDiagonal = 0.5 * (updated[0][1] + updated[1][0]);
            updated[0][1] = offDiagonal;
            updated[1][0] = offDiagonal;
            updated[0][0] += 1.0e-10;
            updated[1][1] += 1.0e-10;
            covariance = updated;

            thetaHistory[k + 1][0] = theta[0];
            thetaHistory[k + 1][1] = theta[1];
        }

        int start = (int) (5.0 / cfg.sampleTime);
        double sumSquaredError = 0.0;
        int count = 0;
        for (int k = start; k <= steps; ++k) {
            double error = reference[k] - y[k];
            sumSquaredError += error * error;
            ++count;
        }
        double rmse = Math.sqrt(sumSquaredError / count);

        System.out.printf("Sampling period: %.6f s%n", cfg.sampleTime);
        System.out.printf("Delay: %d samples%n", cfg.delaySamples);
        System.out.printf("Tracking RMSE after 5 s: %.6f%n", rmse);
        System.out.printf("True [a, b]: [%.6f, %.6f]%n", aTrue, bTrue);
        System.out.printf(
            "Estimated [a, b]: [%.6f, %.6f]%n",
            thetaHistory[steps][0],
            thetaHistory[steps][1]
        );

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter14_Lesson4_results_java.csv"))) {
            writer.write(
                "time,reference,y,y_quantized,u_applied,u_command,a_hat,b_hat\n"
            );
            for (int k = 0; k <= steps; ++k) {
                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,",
                    time[k],
                    reference[k],
                    y[k],
                    yMeasured[k]
                ));
                if (k < steps) {
                    writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12f,",
                        uApplied[k],
                        command[k]
                    ));
                } else {
                    writer.write(",,");
                }
                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f%n",
                    thetaHistory[k][0],
                    thetaHistory[k][1]
                ));
            }
        }
    }
}
      

The MATLAB portion is a complete numerical simulation using base MATLAB. Setting buildSimulinkModel = true programmatically creates a fixed-step Simulink model whose MATLAB Function block contains the plant, RLS estimator, quantizers, delay queue, projection, and controller. The model requires Simulink; fixed-point conversion can later be investigated with Fixed-Point Designer when that product is available.

Chapter14_Lesson4.m

% Chapter14_Lesson4.m
% Discrete-time adaptive control with sampling, quantization, and delay.
% The first part is a base-MATLAB simulation. If Simulink is installed,
% set buildSimulinkModel = true to create Chapter14_Lesson4_Simulink.slx.

clear; clc; close all;

cfg.sampleTime = 0.05;
cfg.duration = 30.0;
cfg.plantTimeConstant = 0.8;
cfg.plantDcGain = 1.25;
cfg.delaySamples = 2;
cfg.sensorStep = 0.01;
cfg.actuatorStep = 0.01;
cfg.actuatorLimit = 3.0;
cfg.forgettingFactor = 0.995;
cfg.desiredPole = 0.65;
cfg.seed = 7;
buildSimulinkModel = false;

rng(cfg.seed);
steps = round(cfg.duration / cfg.sampleTime);
aTrue = exp(-cfg.sampleTime / cfg.plantTimeConstant);
bTrue = cfg.plantDcGain * (1 - aTrue);

time = (0:steps)' * cfg.sampleTime;
reference = arrayfun(@referenceSignal, time);
y = zeros(steps + 1, 1);
yMeasured = zeros(steps + 1, 1);
uApplied = zeros(steps, 1);
command = zeros(steps, 1);
thetaHistory = zeros(steps + 1, 2);

theta = [0.70; 0.30];
P = 100 * eye(2);
thetaHistory(1, :) = theta.';
yMeasured(1) = quantizeValue(y(1), cfg.sensorStep);
queue = zeros(cfg.delaySamples, 1);

for k = 1:steps
    if cfg.delaySamples > 0
        uNow = queue(1);
        queue(1) = [];
    else
        uNow = 0;
    end

    predicted = yMeasured(k);
    if cfg.delaySamples > 0
        predicted = theta(1) * predicted + theta(2) * uNow;
        for j = 1:numel(queue)
            predicted = theta(1) * predicted + theta(2) * queue(j);
        end
    end

    targetNext = cfg.desiredPole * predicted ...
        + (1 - cfg.desiredPole) * reference(k);
    rawCommand = (targetNext - theta(1) * predicted) / theta(2);
    newCommand = quantizeValue( ...
        saturateValue(rawCommand, cfg.actuatorLimit), ...
        cfg.actuatorStep);
    command(k) = newCommand;

    if cfg.delaySamples > 0
        queue(end + 1, 1) = newCommand;
    else
        uNow = newCommand;
    end
    uApplied(k) = uNow;

    disturbance = 0.004 * randn();
    y(k + 1) = aTrue * y(k) + bTrue * uNow + disturbance;
    yMeasured(k + 1) = quantizeValue(y(k + 1), cfg.sensorStep);

    phi = [yMeasured(k); uNow];
    predictionError = yMeasured(k + 1) - phi.' * theta;
    denominator = cfg.forgettingFactor + phi.' * P * phi;
    gain = P * phi / denominator;

    theta = projectParameters(theta + gain * predictionError);
    P = (P - gain * phi.' * P) / cfg.forgettingFactor;
    P = 0.5 * (P + P.') + 1.0e-10 * eye(2);
    thetaHistory(k + 1, :) = theta.';
end

startIndex = floor(5.0 / cfg.sampleTime) + 1;
rmse = sqrt(mean((reference(startIndex:end) - y(startIndex:end)).^2));

fprintf('Sampling period: %.6f s\n', cfg.sampleTime);
fprintf('Delay: %d samples\n', cfg.delaySamples);
fprintf('Tracking RMSE after 5 s: %.6f\n', rmse);
fprintf('True [a, b]: [%.6f, %.6f]\n', aTrue, bTrue);
fprintf('Estimated [a, b]: [%.6f, %.6f]\n', ...
    thetaHistory(end, 1), thetaHistory(end, 2));

resultTable = table( ...
    time, reference, y, yMeasured, ...
    [uApplied; NaN], [command; NaN], ...
    thetaHistory(:, 1), thetaHistory(:, 2), ...
    'VariableNames', {'time', 'reference', 'y', 'yQuantized', ...
    'uApplied', 'uCommand', 'aHat', 'bHat'});
writetable(resultTable, 'Chapter14_Lesson4_results_matlab.csv');

figure;
plot(time, reference, '--', 'DisplayName', 'reference'); hold on;
plot(time, y, 'DisplayName', 'plant output');
stairs(time, yMeasured, 'DisplayName', 'quantized measurement');
grid on; xlabel('time (s)'); ylabel('output'); legend('Location', 'best');
title('Sampling, quantization, and delayed adaptive control');

figure;
plot(time, thetaHistory(:, 1), 'DisplayName', 'a hat'); hold on;
plot(time, thetaHistory(:, 2), 'DisplayName', 'b hat');
yline(aTrue, '--', 'DisplayName', 'a true');
yline(bTrue, '--', 'DisplayName', 'b true');
grid on; xlabel('time (s)'); ylabel('parameter'); legend('Location', 'best');
title('RLS parameter estimates');

if buildSimulinkModel
    if license('test', 'Simulink')
        buildAdaptiveSimulinkModel(cfg);
    else
        warning('Simulink is not licensed; the MATLAB simulation completed.');
    end
end

function value = quantizeValue(value, step)
    if step > 0
        value = step * round(value / step);
    end
end

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

function theta = projectParameters(theta)
    theta(1) = max(-0.98, min(0.999, theta(1)));
    theta(2) = max(0.05, min(2.50, theta(2)));
end

function r = referenceSignal(t)
    if t < 3.0
        r = 0.0;
    elseif t < 12.0
        r = 1.0;
    elseif t < 20.0
        r = -0.6;
    else
        r = 0.8;
    end
end

function buildAdaptiveSimulinkModel(cfg)
    model = 'Chapter14_Lesson4_Simulink';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Step', [model '/Reference'], ...
        'Time', '3', 'Before', '0', 'After', '1', ...
        'Position', [40 90 80 120]);

    adaptiveBlock = add_block( ...
        'simulink/User-Defined Functions/MATLAB Function', ...
        [model '/AdaptiveLoop'], ...
        'Position', [160 55 360 155]);

    add_block('simulink/Sinks/Scope', [model '/OutputScope'], ...
        'Position', [455 40 485 70]);
    add_block('simulink/Sinks/Scope', [model '/InputScope'], ...
        'Position', [455 85 485 115]);
    add_block('simulink/Sinks/Scope', [model '/ParameterScope'], ...
        'NumInputPorts', '2', ...
        'Position', [455 130 485 190]);

    add_line(model, 'Reference/1', 'AdaptiveLoop/1');
    add_line(model, 'AdaptiveLoop/1', 'OutputScope/1');
    add_line(model, 'AdaptiveLoop/2', 'InputScope/1');
    add_line(model, 'AdaptiveLoop/3', 'ParameterScope/1');
    add_line(model, 'AdaptiveLoop/4', 'ParameterScope/2');

    chart = find(sfroot, '-isa', 'Stateflow.EMChart', ...
        'Path', [model '/AdaptiveLoop']);

    nl = newline;
    script = [ ...
        'function [yOut,uOut,aHat,bHat] = fcn(r)' nl ...
        '%#codegen' nl ...
        'persistent yState theta P queue initialized' nl ...
        'Ts = ' num2str(cfg.sampleTime, 17) ';' nl ...
        'tau = ' num2str(cfg.plantTimeConstant, 17) ';' nl ...
        'Kp = ' num2str(cfg.plantDcGain, 17) ';' nl ...
        'd = ' num2str(cfg.delaySamples) ';' nl ...
        'dy = ' num2str(cfg.sensorStep, 17) ';' nl ...
        'du = ' num2str(cfg.actuatorStep, 17) ';' nl ...
        'uMax = ' num2str(cfg.actuatorLimit, 17) ';' nl ...
        'lambda = ' num2str(cfg.forgettingFactor, 17) ';' nl ...
        'ac = ' num2str(cfg.desiredPole, 17) ';' nl ...
        'aTrue = exp(-Ts/tau); bTrue = Kp*(1-aTrue);' nl ...
        'if isempty(initialized)' nl ...
        '  yState = 0; theta = [0.70;0.30]; P = 100*eye(2);' nl ...
        '  queue = zeros(max(d,1),1); initialized = true;' nl ...
        'end' nl ...
        'yq = dy*round(yState/dy);' nl ...
        'if d > 0' nl ...
        '  uNow = queue(1); remaining = queue(2:end);' nl ...
        'else' nl ...
        '  uNow = 0; remaining = zeros(0,1);' nl ...
        'end' nl ...
        'predicted = yq;' nl ...
        'if d > 0' nl ...
        '  predicted = theta(1)*predicted + theta(2)*uNow;' nl ...
        'end' nl ...
        'for i = 1:numel(remaining)' nl ...
        '  predicted = theta(1)*predicted + theta(2)*remaining(i);' nl ...
        'end' nl ...
        'target = ac*predicted + (1-ac)*r;' nl ...
        'cmd = (target-theta(1)*predicted)/theta(2);' nl ...
        'cmd = max(-uMax,min(uMax,cmd)); cmd = du*round(cmd/du);' nl ...
        'if d > 0' nl ...
        '  queue = [remaining;cmd];' nl ...
        'else' nl ...
        '  uNow = cmd;' nl ...
        'end' nl ...
        'yNext = aTrue*yState + bTrue*uNow;' nl ...
        'yqNext = dy*round(yNext/dy); phi = [yq;uNow];' nl ...
        'err = yqNext-phi''*theta;' nl ...
        'gain = P*phi/(lambda+phi''*P*phi);' nl ...
        'theta = theta+gain*err;' nl ...
        'theta(1)=max(-0.98,min(0.999,theta(1)));' nl ...
        'theta(2)=max(0.05,min(2.50,theta(2)));' nl ...
        'P=(P-gain*phi''*P)/lambda; P=0.5*(P+P'')+1e-10*eye(2);' nl ...
        'yState=yNext; yOut=yNext; uOut=uNow;' nl ...
        'aHat=theta(1); bHat=theta(2);' nl ...
        'end' ...
    ];
    chart.Script = script;

    set_param(model, ...
        'Solver', 'FixedStepDiscrete', ...
        'FixedStep', num2str(cfg.sampleTime, 17), ...
        'StopTime', num2str(cfg.duration, 17));
    save_system(model, [model '.slx']);
    open_system(model);
end
      

The notebook uses built-in Wolfram Language matrix, probability, plotting, and export functions. The downloadable file is a textual Notebook[...] expression containing one executable input cell.

Chapter14_Lesson4.nb

Notebook[{
  Cell["Chapter 14, Lesson 4: Sampling, Quantization, and Delay", "Title"],
  Cell[BoxData["(* Chapter14_Lesson4.nb *)
ClearAll[\"Global`*\"];

config = <|
   \"SampleTime\" -> 0.05,
   \"Duration\" -> 30.0,
   \"PlantTimeConstant\" -> 0.8,
   \"PlantDcGain\" -> 1.25,
   \"DelaySamples\" -> 2,
   \"SensorStep\" -> 0.01,
   \"ActuatorStep\" -> 0.01,
   \"ActuatorLimit\" -> 3.0,
   \"ForgettingFactor\" -> 0.995,
   \"DesiredPole\" -> 0.65,
   \"Seed\" -> 7
|>;

quantizeValue[value_, step_] :=
  If[step <= 0.0, value, step Round[value/step]];

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

projectParameters[theta_] := {
   Clip[theta[[1]], {-0.98, 0.999}],
   Clip[theta[[2]], {0.05, 2.50}]
};

referenceSignal[t_] := Piecewise[{
    {0.0, t < 3.0},
    {1.0, t < 12.0},
    {-0.6, t < 20.0}
  }, 0.8];

SeedRandom[config[\"Seed\"]];

steps = Round[config[\"Duration\"]/config[\"SampleTime\"]];
aTrue = Exp[-config[\"SampleTime\"]/config[\"PlantTimeConstant\"]];
bTrue = config[\"PlantDcGain\"] (1.0 - aTrue);

time = N[Range[0, steps] config[\"SampleTime\"]];
reference = referenceSignal /@ time;
y = ConstantArray[0.0, steps + 1];
yMeasured = ConstantArray[0.0, steps + 1];
uApplied = ConstantArray[0.0, steps];
command = ConstantArray[0.0, steps];
thetaHistory = ConstantArray[{0.0, 0.0}, steps + 1];

theta = {0.70, 0.30};
covariance = 100.0 IdentityMatrix[2];
thetaHistory[[1]] = theta;
yMeasured[[1]] = quantizeValue[y[[1]], config[\"SensorStep\"]];
queue = ConstantArray[0.0, config[\"DelaySamples\"]];

Do[
  If[config[\"DelaySamples\"] > 0,
    uNow = First[queue];
    queue = Rest[queue],
    uNow = 0.0
  ];

  predicted = yMeasured[[k]];
  If[config[\"DelaySamples\"] > 0,
    predicted = theta[[1]] predicted + theta[[2]] uNow;
    Do[
      predicted = theta[[1]] predicted + theta[[2]] queuedInput,
      {queuedInput, queue}
    ];
  ];

  targetNext =
    config[\"DesiredPole\"] predicted +
    (1.0 - config[\"DesiredPole\"]) reference[[k]];
  rawCommand = (targetNext - theta[[1]] predicted)/theta[[2]];
  newCommand = quantizeValue[
    saturateValue[rawCommand, config[\"ActuatorLimit\"]],
    config[\"ActuatorStep\"]
  ];
  command[[k]] = newCommand;

  If[config[\"DelaySamples\"] > 0,
    queue = Append[queue, newCommand],
    uNow = newCommand
  ];
  uApplied[[k]] = uNow;

  disturbance = 0.004 RandomVariate[NormalDistribution[0.0, 1.0]];
  y[[k + 1]] = aTrue y[[k]] + bTrue uNow + disturbance;
  yMeasured[[k + 1]] =
    quantizeValue[y[[k + 1]], config[\"SensorStep\"]];

  phi = {yMeasured[[k]], uNow};
  predictionError = yMeasured[[k + 1]] - phi.theta;
  denominator =
    config[\"ForgettingFactor\"] + phi.covariance.phi;
  gain = covariance.phi/denominator;

  theta = projectParameters[theta + gain predictionError];
  covariance =
    (covariance - Outer[Times, gain, phi].covariance)/
      config[\"ForgettingFactor\"];
  covariance =
    0.5 (covariance + Transpose[covariance]) +
    1.0*^-10 IdentityMatrix[2];
  thetaHistory[[k + 1]] = theta,
  {k, 1, steps}
];

startIndex = Floor[5.0/config[\"SampleTime\"]] + 1;
rmse = Sqrt[
   Mean[(reference[[startIndex ;;]] - y[[startIndex ;;]])^2]
];

Print[\"Sampling period: \", NumberForm[config[\"SampleTime\"], {6, 4}], \" s\"];
Print[\"Delay: \", config[\"DelaySamples\"], \" samples\"];
Print[\"Tracking RMSE after 5 s: \", NumberForm[rmse, {8, 6}]];
Print[\"True [a, b]: \", N[{aTrue, bTrue}, 8]];
Print[\"Estimated [a, b]: \", N[Last[thetaHistory], 8]];

resultRows = Prepend[
   Table[{
     time[[k]],
     reference[[k]],
     y[[k]],
     yMeasured[[k]],
     If[k <= steps, uApplied[[k]], Missing[\"NotAvailable\"]],
     If[k <= steps, command[[k]], Missing[\"NotAvailable\"]],
     thetaHistory[[k, 1]],
     thetaHistory[[k, 2]]
   }, {k, 1, steps + 1}],
   {\"time\", \"reference\", \"y\", \"y_quantized\", \"u_applied\",
    \"u_command\", \"a_hat\", \"b_hat\"}
];
Export[\"Chapter14_Lesson4_results_mathematica.csv\", resultRows];

trackingPlot = ListLinePlot[
   {
     Transpose[{time, reference}],
     Transpose[{time, y}],
     Transpose[{time, yMeasured}]
   },
   PlotLegends -> {\"reference\", \"plant output\", \"quantized measurement\"},
   Frame -> True,
   FrameLabel -> {\"time (s)\", \"output\"},
   PlotRange -> All,
   ImageSize -> Large
];

parameterPlot = ListLinePlot[
   {
     Transpose[{time, thetaHistory[[All, 1]]}],
     Transpose[{time, thetaHistory[[All, 2]]}],
     Transpose[{time, ConstantArray[aTrue, steps + 1]}],
     Transpose[{time, ConstantArray[bTrue, steps + 1]}]
   },
   PlotLegends -> {\"a hat\", \"b hat\", \"a true\", \"b true\"},
   Frame -> True,
   FrameLabel -> {\"time (s)\", \"parameter\"},
   PlotRange -> All,
   ImageSize -> Large
];

Column[{trackingPlot, parameterPlot}]
"], "Input"]
}, WindowTitle -> "Chapter14_Lesson4"]
      

Problem 1 (Exact sampled parameters): A process satisfies \( 0.5\dot y+y=2u \). Find the exact discrete-time coefficients for zero-order-hold sampling at \( T_s=0.1 \) s.

Solution:

\[ a=e^{-T_s/\tau_p}=e^{-0.1/0.5}=e^{-0.2}\approx0.81873, \]

\[ b=K_p(1-a)=2(1-0.81873)\approx0.36254. \]

The exact model is \( y_{k+1}=0.81873y_k+0.36254u_k \). Using the continuous coefficients directly in an RLS regression would estimate the wrong parameterization.

Problem 2 (Quantization disturbance): For \( a=0.9 \), \( b=0.2 \), \( \Delta_y=0.01 \), and \( \Delta_u=0.02 \), compute an upper bound on \( |d_k^q| \).

Solution:

\[ |d_k^q| \le \frac{1+0.9}{2}(0.01) +\frac{0.2}{2}(0.02) =0.0095+0.002=0.0115. \]

The result is a deterministic worst-case bound before converter saturation. It should not be interpreted as a standard deviation.

Problem 3 (Delay phase budget): A controller has desired bandwidth \( \omega_b=12 \) rad/s and a two-sample total delay. Restrict delay phase to 20 degrees. Find the maximum sampling period.

Solution: Convert 20 degrees to radians:

\[ \phi_{\mathrm{allow} }=20\frac{\pi}{180} \approx0.3491. \]

\[ T_s\le \frac{0.3491}{2(12)} \approx0.01455\;\text{s}. \]

Thus a nominal choice such as 10 ms leaves additional margin for model uncertainty and adaptive transients.

Problem 4 (Ultimate tracking bound): Suppose \( a_c=0.6 \) and the combined implementation perturbation satisfies \( |\zeta_k|\le0.02 \). Find the asymptotic geometric bound.

Solution:

\[ \limsup_{k→\infty}|e_k| \le \frac{0.02}{1-0.6} =0.05. \]

Exact convergence cannot be guaranteed from this bounded-disturbance model; the result is convergence to a neighborhood of radius 0.05.

Problem 5 (Regressor alignment): A command computed at sample \( k \) is applied two samples later. A programmer forms \( \varphi_k=[y_k^q,c_k]^T \) for estimating \( y_{k+1} \). Explain the error and give the correct regressor.

Solution:

During interval \( [kT_s,(k+1)T_s) \), the plant receives \( u_k=c_{k-2} \), not \( c_k \). The correct regression is

\[ y_{k+1}^q = \begin{bmatrix}y_k^q & u_k\end{bmatrix} \begin{bmatrix}a\\b\end{bmatrix} +\text{disturbance}, \qquad u_k=c_{k-2}. \]

Using \( c_k \) introduces a deterministic input-output misalignment. RLS then attempts to fit delay error as though it were a change in \( a \) or \( b \).

Problem 6 (Jitter-induced coefficient variation): For \( \tau_p=0.4 \) s, \( T_s=0.02 \) s, and \( J=0.001 \) s, bound the variation in \( a_k=e^{-T_k/\tau_p} \).

Solution:

\[ |a_k-a(T_s)| \le \frac{0.001}{0.4} e^{-(0.02-0.001)/0.4} \approx0.00238. \]

Even though the physical time constant is constant, the sampled coefficient can fluctuate by approximately 0.00238 under this jitter bound.

Sampling determines the plant parameterization seen by the adaptive law, and exact zero-order-hold discretization separates physical uncertainty from numerical approximation. Quantization enters the regression as a bounded, signal-dependent disturbance and commonly changes asymptotic convergence into practical convergence. Computation and transport delay require command-queue modeling, use of the actually applied input in the estimator, and prediction to the command-application instant. Projection, covariance conditioning, saturation handling, deadline monitoring, and deterministic indexing are essential parts of the adaptive controller, not optional software refinements.