After completing this lesson, students should be able to:

• derive an exact zero-order-hold model of a first-order velocity plant;
• construct a discrete reference model with a specified time constant;
• derive a certainty-equivalent adaptive velocity controller;
• prove monotonicity of a normalized-gradient parameter-error function;
• distinguish tracking convergence from physical-parameter convergence;
• incorporate projection, dead zones, quantization, and saturation safely; and
• implement and test the complete sampled-data loop.

We assume that an inner current or torque loop is sufficiently fast that its command \(u(t)\) may be treated as the effective motor current or normalized torque command. The outer-loop output is angular velocity \(\omega(t)\). The continuous mechanical model is

\[ J\dot{\omega}(t) = -B\omega(t) + K_t u(t) - \tau_L, \qquad J > 0,\; B > 0,\; K_t > 0. \]

The sign of the input gain is therefore known and positive, but the numerical values of \(J\), \(B\), \(K_t\), and the approximately constant load torque \(\tau_L\) are not assumed known to the controller.

During one sampling interval, let the input be held constant by a zero-order hold. Dividing the motor equation by \(J\) gives

\[ \dot{\omega}(t) = -\frac{B}{J}\omega(t) + \frac{K_t}{J}u(t) - \frac{\tau_L}{J}. \]

For a sample period \(T_s\), exact integration over \([kT_s,(k+1)T_s]\) yields

\[ \omega_{k+1} = a\omega_k + bu_k + d, \]

\[ a = \exp\!\left(-\frac{BT_s}{J}\right), \qquad b = \frac{K_t}{B}(1-a), \qquad d = -\frac{\tau_L}{B}(1-a). \]

The limit for negligible viscous friction is obtained from \(1-e^{-x}\approx x\):

\[ B \approx 0: \qquad a=1, \qquad b=\frac{K_tT_s}{J}, \qquad d=-\frac{\tau_LT_s}{J}. \]

Define the unknown plant-parameter vector and measured regressor as

\[ \boldsymbol{\theta} = \begin{bmatrix}a & b & d\end{bmatrix}^{\!T}, \qquad \boldsymbol{\psi}_k = \begin{bmatrix}y_k & u_k & 1\end{bmatrix}^{\!T}, \]

\[ y_{k+1}=\boldsymbol{\theta}^{T}\boldsymbol{\psi}_k, \]

where \(y_k\) is the sampled encoder velocity. Noise and quantization will be added later; the first derivations use the nominal regression equation.

Choose a stable first-order reference model with unit steady-state gain:

\[ y_{m,k+1}=a_m y_{m,k}+b_m r_k, \qquad 0<a_m<1, \qquad b_m=1-a_m. \]

A desired continuous time constant \(\tau_m\) can be mapped exactly to discrete time by

\[ a_m=\exp\!\left(-\frac{T_s}{\tau_m}\right), \qquad b_m=1-a_m. \]

If the plant parameters were known, one-step model matching would impose

\[ a y_k+b u_k+d = a_m y_k+b_m r_k. \]

Solving for the ideal control gives

\[ u_k^{\star}=\frac{(a_m-a)y_k+b_mr_k-d}{b}. \]

Replacing the unknown parameters by their current estimates produces the certainty-equivalent adaptive controller

\[ u_{ce,k}=\frac{(a_m-\hat a_k)y_k+b_mr_k-\hat d_k}{\hat b_k}. \]

Because division by a small or sign-incorrect estimate is unsafe, the estimate is projected so that \(0<b_{\min}\leq\hat b_k\leq b_{\max}\). The implemented actuator command is

\[ u_k=\operatorname{sat}_{[-u_{\max},u_{\max}]}(u_{ce,k}). \]

flowchart TD
  R["Velocity command r(k)"] --> RM["Reference model"]
  R --> C["Certainty-equivalent controller"]
  YM["Measured velocity y(k)"] --> C
  TH["Parameter estimates aHat, bHat, dHat"] --> C
  C --> SAT["Input saturation"]
  SAT --> P["Sampled motor and load"]
  P --> ENC["Encoder and quantizer"]
  ENC --> YM
  SAT --> ID["Normalized-gradient estimator"]
  YM --> ID
  ID --> TH
  RM --> E["Tracking comparison"]
  YM --> E
        

Using the estimate available before processing sample \(k+1\), define the one-step prediction and prediction error as

\[ \hat y_{k+1|k}=\hat{\boldsymbol{\theta}}_k^T \boldsymbol{\psi}_k, \qquad \varepsilon_{k+1}=y_{k+1}-\hat y_{k+1|k}. \]

The normalized-gradient update is

\[ m_k=\delta+\boldsymbol{\psi}_k^T\boldsymbol{\psi}_k, \qquad \delta>0, \]

\[ \hat{\boldsymbol{\theta}}_{k+1}^{\circ} =\hat{\boldsymbol{\theta}}_k+ \gamma\frac{\boldsymbol{\psi}_k\varepsilon_{k+1}}{m_k}, \qquad 0<\gamma<2. \]

The superscript \(\circ\) denotes the unprojected update. The actual estimate is

\[ \hat{\boldsymbol{\theta}}_{k+1} =\operatorname{Proj}_{\Omega} \left(\hat{\boldsymbol{\theta}}_{k+1}^{\circ}\right), \]

\[ \Omega=[a_{\min},a_{\max}]\times [b_{\min},b_{\max}]\times[d_{\min},d_{\max}]. \]

With noisy and quantized measurements, adaptation is disabled inside a dead zone:

\[ |\varepsilon_{k+1}|\leq\varepsilon_0 \quad \Longrightarrow \quad \hat{\boldsymbol{\theta}}_{k+1} =\hat{\boldsymbol{\theta}}_k. \]

Normalization prevents a large regressor from producing an excessively large parameter step. Projection prevents denominator collapse and encodes prior physical bounds. The dead zone prevents continuous random walk caused by encoder noise.

First neglect noise, quantization, saturation, and parameter variation. Define the parameter error using the true parameter minus the estimate:

\[ \widetilde{\boldsymbol{\theta}}_k =\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}_k. \]

Then the prediction error satisfies

\[ \varepsilon_{k+1} =\widetilde{\boldsymbol{\theta}}_k^T\boldsymbol{\psi}_k. \]

Before projection, the parameter-error recursion is

\[ \widetilde{\boldsymbol{\theta}}_{k+1}^{\circ} =\widetilde{\boldsymbol{\theta}}_k- \gamma\frac{\boldsymbol{\psi}_k\varepsilon_{k+1}}{m_k}. \]

Choose the discrete Lyapunov-like function

\[ V_k=\widetilde{\boldsymbol{\theta}}_k^T \widetilde{\boldsymbol{\theta}}_k. \]

Expanding one update gives

\[ \begin{aligned} V_{k+1}^{\circ}-V_k &=-2\gamma\frac{\varepsilon_{k+1}^2}{m_k} +\gamma^2\frac{\boldsymbol{\psi}_k^T\boldsymbol{\psi}_k} {m_k^2}\varepsilon_{k+1}^2 \\ &=-\gamma\frac{2\delta+(2-\gamma) \boldsymbol{\psi}_k^T\boldsymbol{\psi}_k}{m_k^2} \varepsilon_{k+1}^2. \end{aligned} \]

Since \(\delta>0\) and \(0<\gamma<2\),

\[ V_{k+1}^{\circ}-V_k\leq 0. \]

Moreover,

\[ 2\delta+(2-\gamma)\|\boldsymbol{\psi}_k\|^2 \geq (2-\gamma)m_k, \]

and therefore

\[ V_{k+1}^{\circ}-V_k \leq-\gamma(2-\gamma) \frac{\varepsilon_{k+1}^2}{m_k}. \]

Summing from \(k=0\) to \(N\) gives

\[ \gamma(2-\gamma)\sum_{k=0}^{N} \frac{\varepsilon_{k+1}^2}{m_k}\leq V_0. \]

Thus the parameter estimate is bounded and the normalized squared prediction error is summable. If the true parameter belongs to the convex projection set \(\Omega\), Euclidean projection is nonexpansive relative to that true parameter:

\[ \left\|\boldsymbol{\theta}- \operatorname{Proj}_{\Omega}(\mathbf{z})\right\| \leq\left\|\boldsymbol{\theta}-\mathbf{z}\right\|. \]

Consequently, projection does not destroy the nonincrease property of \(V_k\). For bounded regressors, the above summability implies \(\varepsilon_{k+1}\) approaches zero in the nominal constant-parameter case.

Let \(e_k=y_k-y_{m,k}\). Without saturation, substitution of the certainty-equivalent control law into the true plant gives

\[ y_{k+1}=a_my_k+b_mr_k+ \widetilde{\boldsymbol{\theta}}_k^T\boldsymbol{\psi}_k. \]

Subtracting the reference model yields

\[ e_{k+1}=a_me_k+\varepsilon_{k+1}. \]

Repeated substitution gives the stable convolution

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

Hence, if \(|a_m|<1\) and the prediction mismatch approaches zero, then the tracking error also approaches zero. With a dead zone, measurement noise, slowly varying parameters, or saturation, one generally expects convergence to a bounded neighborhood rather than exact asymptotic convergence.

This conclusion does not imply that each estimated physical parameter reaches its true value. Parameter convergence requires a persistent-excitation condition. One discrete form is: there exist an integer \(N_e\) and a scalar \(\alpha>0\) such that

\[ \sum_{i=k}^{k+N_e-1} \boldsymbol{\psi}_i\boldsymbol{\psi}_i^T \succeq\alpha\mathbf{I} \qquad \text{for every }k. \]

A constant or nearly constant velocity command often fails to excite all three directions associated with \(a\), \(b\), and \(d\). The closed loop can nevertheless track well because many parameter combinations produce almost the same input-output behavior over the observed trajectory.

Under saturation, define \(\Delta u_k=u_k-u_{ce,k}\). The tracking equation becomes

\[ e_{k+1}=a_me_k+ \widetilde{\boldsymbol{\theta}}_k^T\boldsymbol{\psi}_k +\hat b_k\Delta u_k. \]

The saturation residual is therefore an additional bounded forcing term. This is why aggressive reference models can worsen practical tracking even when the nominal adaptive derivation is correct.

7.1 Sampling period

A smaller \(T_s\) improves temporal resolution but can make consecutive regressors nearly collinear, reduce the information per sample, and amplify the influence of encoder quantization. A larger \(T_s\) increases discrete input gain but may violate the assumed first-order model or introduce excessive phase lag.

7.2 Encoder quantization and noise

For encoder resolution \(q_y\), the simulated measurement is

\[ y_k^{q}=q_y\operatorname{round} \left(\frac{y_k+\nu_k}{q_y}\right). \]

The estimator uses \(y_k^q\) and a dead zone slightly larger than the typical combined noise and quantization error.

7.3 Computation and transport delay

If the command calculated at sample \(k\) is not applied until sample \(k+1\), the correct regression is

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

The estimator must then use \(\boldsymbol{\psi}_k=[y_k,u_{k-1},1]^T\), and the one-step controller derived for an undelayed input is no longer exact. Ignoring this index shift can create biased estimates or oscillation.

7.4 Required safeguards

• project \(\hat b_k\) away from zero;
• limit the actuator command and its rate if required;
• use the actual saturated input in the estimator regressor;
• freeze or slow adaptation under sensor faults and severe saturation;
• monitor residuals and parameter-bound contact; and
• retain a verified fallback controller for hardware experiments.

flowchart TD
  A["Read encoder sample y(k)"] --> B["Quantize and validate measurement"]
  B --> C["Read projected estimates"]
  C --> D["Compute certainty-equivalent input"]
  D --> E["Apply saturation and output u(k)"]
  E --> F["Propagate or observe motor to y(k+1)"]
  F --> G["Build regressor from y(k), u(k), and constant term"]
  G --> H["Compute one-step prediction error"]
  H --> I["Apply dead-zone test"]
  I --> J["Normalized parameter update"]
  J --> K["Project estimates into physical bounds"]
  K --> L["Update reference model and log diagnostics"]
  L --> A
        

On real hardware, the sequence is split across interrupt boundaries. The parameter update associated with \(u_k\) can only be completed after \(y_{k+1}\) becomes available. The implementation must preserve that indexing exactly.

The supplied implementations use:

• \(T_s=0.02\text{ s}\) and \(\tau_m=0.20\text{ s}\);
• \(\gamma=0.55\), \(\delta=1\), and dead zone \(0.015\);
• input saturation at \(\pm5\);
• encoder quantum \(0.01\text{ rad/s}\); and
• a plant and load change at \(t=15\text{ s}\).

Before the change:

\[ J=0.020,\; B=0.080,\;K_t=0.500,\;\tau_L=0.040, \]

\[ a\approx0.923116,\qquad b\approx0.480525, \qquad d\approx-0.038442. \]

After the change:

\[ J=0.032,\; B=0.100,\;K_t=0.450,\;\tau_L=0.100, \]

\[ a\approx0.939413,\qquad b\approx0.272641, \qquad d\approx-0.060587. \]

The selected reference model has

\[ a_m=e^{-0.1}\approx0.904837, \qquad b_m\approx0.095163. \]

The command includes steps followed by sinusoidal segments. The richer final segment improves information content, but closed-loop correlation can still prevent exact convergence of every physical parameter. Students should compare the tracking error, prediction error, parameter estimates, and saturation activity rather than judging adaptation from one plot.

This implementation uses NumPy for arrays and random-number generation, the standard csv module for export, and optional Matplotlib for plots. Install plotting support with python -m pip install numpy matplotlib.

Chapter14_Lesson5.py

"""Chapter14_Lesson5.py

Discrete-time indirect adaptive control of a DC-motor velocity loop.
The plant is represented by y[k+1] = a*y[k] + b*u[k] + d.
A normalized-gradient estimator updates [a, b, d], and a
certainty-equivalent controller enforces a first-order reference model.
"""

from __future__ import annotations

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

import numpy as np


@dataclass(frozen=True)
class MotorParameters:
    inertia: float
    viscous_friction: float
    torque_constant: float
    load_torque: float


def zoh_parameters(motor: MotorParameters, sample_time: float) -> tuple[float, float, float]:
    """Return exact ZOH parameters for J*w_dot = -B*w + Kt*u - tau_L."""
    if motor.inertia <= 0.0:
        raise ValueError("Inertia must be positive.")
    if sample_time <= 0.0:
        raise ValueError("Sample time must be positive.")

    if abs(motor.viscous_friction) < 1.0e-12:
        a = 1.0
        b = motor.torque_constant * sample_time / motor.inertia
        d = -motor.load_torque * sample_time / motor.inertia
    else:
        a = math.exp(-motor.viscous_friction * sample_time / motor.inertia)
        one_minus_a = 1.0 - a
        b = (motor.torque_constant / motor.viscous_friction) * one_minus_a
        d = -(motor.load_torque / motor.viscous_friction) * one_minus_a
    return a, b, d


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


def reference_command(time_s: float) -> float:
    if time_s < 2.0:
        return 0.0
    if time_s < 8.0:
        return 6.0
    if time_s < 14.0:
        return 10.0
    if time_s < 20.0:
        return 4.0
    if time_s < 26.0:
        return 7.0 + 2.0 * math.sin(2.0 * math.pi * 0.40 * (time_s - 20.0))
    tau = time_s - 26.0
    return 7.0 + 1.5 * math.sin(2.0 * math.pi * 0.65 * tau) + 0.6 * math.sin(2.0 * math.pi * 1.35 * tau)


def main() -> None:
    sample_time = 0.02
    final_time = 30.0
    samples = int(final_time / sample_time)

    # Reference model: y_m[k+1] = a_m*y_m[k] + b_m*r[k], with unit DC gain.
    reference_time_constant = 0.20
    a_m = math.exp(-sample_time / reference_time_constant)
    b_m = 1.0 - a_m

    # Normalized-gradient estimator settings.
    gamma = 0.55                 # theorem requires 0 < gamma < 2
    normalization_offset = 1.0   # delta > 0
    error_dead_zone = 0.015

    # Projection set; known positive control direction is enforced by b_hat >= b_min.
    a_bounds = (0.0, 0.9995)
    b_bounds = (0.05, 1.50)
    d_bounds = (-1.50, 1.50)

    input_limit = 5.0
    encoder_quantum = 0.01
    measurement_noise_std = 0.004
    process_noise_std = 0.001

    motor_before = MotorParameters(0.020, 0.080, 0.500, 0.040)
    motor_after = MotorParameters(0.032, 0.100, 0.450, 0.100)
    change_time = 15.0

    rng = np.random.default_rng(14)

    time = np.arange(samples + 1, dtype=float) * sample_time
    r = np.zeros(samples)
    y = np.zeros(samples + 1)
    y_measured = np.zeros(samples + 1)
    y_model = np.zeros(samples + 1)
    u = np.zeros(samples)
    prediction_error = np.zeros(samples)
    theta_hat = np.zeros((samples + 1, 3))
    theta_true = np.zeros((samples, 3))

    theta_hat[0] = np.array([0.80, 0.35, 0.0])

    for k in range(samples):
        t = time[k]
        r[k] = reference_command(t)
        motor = motor_before if t < change_time else motor_after
        a, b, d = zoh_parameters(motor, sample_time)
        theta_true[k] = np.array([a, b, d])

        noisy_velocity = y[k] + rng.normal(0.0, measurement_noise_std)
        y_measured[k] = encoder_quantum * round(noisy_velocity / encoder_quantum)

        a_hat, b_hat, d_hat = theta_hat[k]
        b_safe = max(b_hat, b_bounds[0])

        # Certainty-equivalent one-step model matching.
        u_unsaturated = ((a_m - a_hat) * y_measured[k] + b_m * r[k] - d_hat) / b_safe
        u[k] = clip(u_unsaturated, -input_limit, input_limit)

        y[k + 1] = (
            a * y[k]
            + b * u[k]
            + d
            + rng.normal(0.0, process_noise_std)
        )
        y_model[k + 1] = a_m * y_model[k] + b_m * r[k]

        noisy_next = y[k + 1] + rng.normal(0.0, measurement_noise_std)
        y_measured[k + 1] = encoder_quantum * round(noisy_next / encoder_quantum)

        psi = np.array([y_measured[k], u[k], 1.0])
        prediction_error[k] = y_measured[k + 1] - float(theta_hat[k] @ psi)

        theta_next = theta_hat[k].copy()
        if abs(prediction_error[k]) > error_dead_zone:
            normalizer = normalization_offset + float(psi @ psi)
            theta_next += gamma * psi * prediction_error[k] / normalizer

        # Componentwise projection prevents an unsafe denominator and implausible estimates.
        theta_next[0] = clip(theta_next[0], *a_bounds)
        theta_next[1] = clip(theta_next[1], *b_bounds)
        theta_next[2] = clip(theta_next[2], *d_bounds)
        theta_hat[k + 1] = theta_next

    tracking_error = y - y_model
    rms_after_startup = math.sqrt(float(np.mean(tracking_error[int(2.0 / sample_time):] ** 2)))
    print(f"RMS tracking error after 2 s: {rms_after_startup:.4f} rad/s")
    print("Final estimate [a_hat, b_hat, d_hat]:", theta_hat[-1])
    print("Final true parameters [a, b, d]:", theta_true[-1])

    output_path = Path(__file__).with_name("Chapter14_Lesson5_results.csv")
    with output_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow([
            "time_s", "reference", "velocity", "reference_model", "control",
            "a_hat", "b_hat", "d_hat", "prediction_error"
        ])
        for k in range(samples):
            writer.writerow([
                time[k], r[k], y[k], y_model[k], u[k],
                theta_hat[k, 0], theta_hat[k, 1], theta_hat[k, 2], prediction_error[k]
            ])
        writer.writerow([
            time[-1], "", y[-1], y_model[-1], "",
            theta_hat[-1, 0], theta_hat[-1, 1], theta_hat[-1, 2], ""
        ])
    print(f"Saved: {output_path}")

    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print("Install matplotlib to display plots: python -m pip install matplotlib")
        return

    plt.figure()
    plt.plot(time[:-1], r, label="command")
    plt.plot(time, y, label="velocity")
    plt.plot(time, y_model, "--", label="reference model")
    plt.xlabel("Time (s)")
    plt.ylabel("Velocity (rad/s)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(time[:-1], u)
    plt.xlabel("Time (s)")
    plt.ylabel("Control input")
    plt.grid(True)
    plt.tight_layout()

    plt.figure()
    plt.plot(time, theta_hat[:, 0], label="a_hat")
    plt.plot(time, theta_hat[:, 1], label="b_hat")
    plt.plot(time, theta_hat[:, 2], label="d_hat")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. Compile with g++ -std=c++17 -O2 Chapter14_Lesson5.cpp -o Chapter14_Lesson5. It exports a CSV file suitable for plotting in Python, MATLAB, or a spreadsheet.

Chapter14_Lesson5.cpp

// Chapter14_Lesson5.cpp
// Discrete-time indirect adaptive velocity control using a normalized-gradient estimator.

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

struct MotorParameters {
    double inertia;
    double viscousFriction;
    double torqueConstant;
    double loadTorque;
};

struct DiscretePlant {
    double a;
    double b;
    double d;
};

DiscretePlant zohParameters(const MotorParameters& motor, double sampleTime) {
    if (motor.inertia <= 0.0 || sampleTime <= 0.0) {
        throw std::invalid_argument("Inertia and sample time must be positive.");
    }
    if (std::abs(motor.viscousFriction) < 1.0e-12) {
        return {1.0,
                motor.torqueConstant * sampleTime / motor.inertia,
                -motor.loadTorque * sampleTime / motor.inertia};
    }
    const double a = std::exp(-motor.viscousFriction * sampleTime / motor.inertia);
    const double oneMinusA = 1.0 - a;
    const double b = (motor.torqueConstant / motor.viscousFriction) * oneMinusA;
    const double d = -(motor.loadTorque / motor.viscousFriction) * oneMinusA;
    return {a, b, d};
}

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

double referenceCommand(double timeSeconds) {
    constexpr double pi = 3.14159265358979323846;
    if (timeSeconds < 2.0) return 0.0;
    if (timeSeconds < 8.0) return 6.0;
    if (timeSeconds < 14.0) return 10.0;
    if (timeSeconds < 20.0) return 4.0;
    if (timeSeconds < 26.0) {
        return 7.0 + 2.0 * std::sin(2.0 * pi * 0.40 * (timeSeconds - 20.0));
    }
    const double tau = timeSeconds - 26.0;
    return 7.0
        + 1.5 * std::sin(2.0 * pi * 0.65 * tau)
        + 0.6 * std::sin(2.0 * pi * 1.35 * tau);
}

int main() {
    try {
        const double sampleTime = 0.02;
        const double finalTime = 30.0;
        const int samples = static_cast<int>(finalTime / sampleTime);

        const double referenceTimeConstant = 0.20;
        const double aModel = std::exp(-sampleTime / referenceTimeConstant);
        const double bModel = 1.0 - aModel;

        const double gamma = 0.55;
        const double normalizationOffset = 1.0;
        const double errorDeadZone = 0.015;

        const std::array<double, 2> aBounds{0.0, 0.9995};
        const std::array<double, 2> bBounds{0.05, 1.50};
        const std::array<double, 2> dBounds{-1.50, 1.50};

        const double inputLimit = 5.0;
        const double encoderQuantum = 0.01;
        const double measurementNoiseStd = 0.004;
        const double processNoiseStd = 0.001;

        const MotorParameters motorBefore{0.020, 0.080, 0.500, 0.040};
        const MotorParameters motorAfter{0.032, 0.100, 0.450, 0.100};
        const double changeTime = 15.0;

        std::mt19937_64 generator(14);
        std::normal_distribution<double> measurementNoise(0.0, measurementNoiseStd);
        std::normal_distribution<double> processNoise(0.0, processNoiseStd);

        std::vector<double> time(samples + 1, 0.0);
        std::vector<double> reference(samples, 0.0);
        std::vector<double> velocity(samples + 1, 0.0);
        std::vector<double> measuredVelocity(samples + 1, 0.0);
        std::vector<double> modelVelocity(samples + 1, 0.0);
        std::vector<double> control(samples, 0.0);
        std::vector<double> predictionError(samples, 0.0);
        std::vector<std::array<double, 3>> thetaHat(samples + 1);
        std::vector<std::array<double, 3>> thetaTrue(samples);

        thetaHat[0] = {0.80, 0.35, 0.0};

        for (int k = 0; k < samples; ++k) {
            time[k] = static_cast<double>(k) * sampleTime;
            const double t = time[k];
            reference[k] = referenceCommand(t);

            const MotorParameters& motor = (t < changeTime) ? motorBefore : motorAfter;
            const DiscretePlant plant = zohParameters(motor, sampleTime);
            thetaTrue[k] = {plant.a, plant.b, plant.d};

            const double noisyVelocity = velocity[k] + measurementNoise(generator);
            measuredVelocity[k] = encoderQuantum * std::round(noisyVelocity / encoderQuantum);

            const double aHat = thetaHat[k][0];
            const double bHat = thetaHat[k][1];
            const double dHat = thetaHat[k][2];
            const double safeB = std::max(bHat, bBounds[0]);

            const double unsaturatedControl =
                ((aModel - aHat) * measuredVelocity[k] + bModel * reference[k] - dHat) / safeB;
            control[k] = clip(unsaturatedControl, -inputLimit, inputLimit);

            velocity[k + 1] = plant.a * velocity[k]
                            + plant.b * control[k]
                            + plant.d
                            + processNoise(generator);
            modelVelocity[k + 1] = aModel * modelVelocity[k] + bModel * reference[k];

            const double noisyNext = velocity[k + 1] + measurementNoise(generator);
            measuredVelocity[k + 1] = encoderQuantum * std::round(noisyNext / encoderQuantum);

            const std::array<double, 3> psi{measuredVelocity[k], control[k], 1.0};
            const double prediction = thetaHat[k][0] * psi[0]
                                    + thetaHat[k][1] * psi[1]
                                    + thetaHat[k][2] * psi[2];
            predictionError[k] = measuredVelocity[k + 1] - prediction;

            std::array<double, 3> next = thetaHat[k];
            if (std::abs(predictionError[k]) > errorDeadZone) {
                const double normalizer = normalizationOffset
                                        + psi[0] * psi[0]
                                        + psi[1] * psi[1]
                                        + psi[2] * psi[2];
                for (int i = 0; i < 3; ++i) {
                    next[i] += gamma * psi[i] * predictionError[k] / normalizer;
                }
            }

            next[0] = clip(next[0], aBounds[0], aBounds[1]);
            next[1] = clip(next[1], bBounds[0], bBounds[1]);
            next[2] = clip(next[2], dBounds[0], dBounds[1]);
            thetaHat[k + 1] = next;
        }
        time[samples] = finalTime;

        double squaredErrorSum = 0.0;
        const int startIndex = static_cast<int>(2.0 / sampleTime);
        for (int k = startIndex; k <= samples; ++k) {
            const double error = velocity[k] - modelVelocity[k];
            squaredErrorSum += error * error;
        }
        const double rmsError = std::sqrt(squaredErrorSum / (samples + 1 - startIndex));

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "RMS tracking error after 2 s: " << rmsError << " rad/s\n";
        std::cout << "Final estimate [a_hat, b_hat, d_hat]: ["
                  << thetaHat.back()[0] << ", " << thetaHat.back()[1] << ", "
                  << thetaHat.back()[2] << "]\n";
        std::cout << "Final true parameters [a, b, d]: ["
                  << thetaTrue.back()[0] << ", " << thetaTrue.back()[1] << ", "
                  << thetaTrue.back()[2] << "]\n";

        std::ofstream csv("Chapter14_Lesson5_results_cpp.csv");
        if (!csv) {
            throw std::runtime_error("Could not create CSV output file.");
        }
        csv << "time_s,reference,velocity,reference_model,control,a_hat,b_hat,d_hat,prediction_error\n";
        csv << std::setprecision(12);
        for (int k = 0; k < samples; ++k) {
            csv << time[k] << ',' << reference[k] << ',' << velocity[k] << ','
                << modelVelocity[k] << ',' << control[k] << ',' << thetaHat[k][0] << ','
                << thetaHat[k][1] << ',' << thetaHat[k][2] << ',' << predictionError[k] << '\n';
        }
        csv << time.back() << ",," << velocity.back() << ',' << modelVelocity.back()
            << ",," << thetaHat.back()[0] << ',' << thetaHat.back()[1] << ','
            << thetaHat.back()[2] << ",\n";
        std::cout << "Saved: Chapter14_Lesson5_results_cpp.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

The Java version uses the standard runtime only. Compile and run with javac Chapter14_Lesson5.java followed by java Chapter14_Lesson5.

Chapter14_Lesson5.java

// Chapter14_Lesson5.java
// Discrete-time indirect adaptive control of a DC-motor velocity loop.

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

public final class Chapter14_Lesson5 {
    private static final class MotorParameters {
        final double inertia;
        final double viscousFriction;
        final double torqueConstant;
        final double loadTorque;

        MotorParameters(double inertia, double viscousFriction,
                        double torqueConstant, double loadTorque) {
            this.inertia = inertia;
            this.viscousFriction = viscousFriction;
            this.torqueConstant = torqueConstant;
            this.loadTorque = loadTorque;
        }
    }

    private static final class DiscretePlant {
        final double a;
        final double b;
        final double d;

        DiscretePlant(double a, double b, double d) {
            this.a = a;
            this.b = b;
            this.d = d;
        }
    }

    private Chapter14_Lesson5() {}

    private static DiscretePlant zohParameters(MotorParameters motor, double sampleTime) {
        if (motor.inertia <= 0.0 || sampleTime <= 0.0) {
            throw new IllegalArgumentException("Inertia and sample time must be positive.");
        }
        if (Math.abs(motor.viscousFriction) < 1.0e-12) {
            return new DiscretePlant(
                    1.0,
                    motor.torqueConstant * sampleTime / motor.inertia,
                    -motor.loadTorque * sampleTime / motor.inertia);
        }
        double a = Math.exp(-motor.viscousFriction * sampleTime / motor.inertia);
        double oneMinusA = 1.0 - a;
        double b = (motor.torqueConstant / motor.viscousFriction) * oneMinusA;
        double d = -(motor.loadTorque / motor.viscousFriction) * oneMinusA;
        return new DiscretePlant(a, b, d);
    }

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

    private static double referenceCommand(double timeSeconds) {
        if (timeSeconds < 2.0) return 0.0;
        if (timeSeconds < 8.0) return 6.0;
        if (timeSeconds < 14.0) return 10.0;
        if (timeSeconds < 20.0) return 4.0;
        if (timeSeconds < 26.0) {
            return 7.0 + 2.0 * Math.sin(2.0 * Math.PI * 0.40 * (timeSeconds - 20.0));
        }
        double tau = timeSeconds - 26.0;
        return 7.0
                + 1.5 * Math.sin(2.0 * Math.PI * 0.65 * tau)
                + 0.6 * Math.sin(2.0 * Math.PI * 1.35 * tau);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);

        double sampleTime = 0.02;
        double finalTime = 30.0;
        int samples = (int) (finalTime / sampleTime);

        double referenceTimeConstant = 0.20;
        double aModel = Math.exp(-sampleTime / referenceTimeConstant);
        double bModel = 1.0 - aModel;

        double gamma = 0.55;
        double normalizationOffset = 1.0;
        double errorDeadZone = 0.015;

        double[] aBounds = {0.0, 0.9995};
        double[] bBounds = {0.05, 1.50};
        double[] dBounds = {-1.50, 1.50};

        double inputLimit = 5.0;
        double encoderQuantum = 0.01;
        double measurementNoiseStd = 0.004;
        double processNoiseStd = 0.001;

        MotorParameters motorBefore = new MotorParameters(0.020, 0.080, 0.500, 0.040);
        MotorParameters motorAfter = new MotorParameters(0.032, 0.100, 0.450, 0.100);
        double changeTime = 15.0;

        Random random = new Random(14L);

        double[] time = new double[samples + 1];
        double[] reference = new double[samples];
        double[] velocity = new double[samples + 1];
        double[] measuredVelocity = new double[samples + 1];
        double[] modelVelocity = new double[samples + 1];
        double[] control = new double[samples];
        double[] predictionError = new double[samples];
        double[][] thetaHat = new double[samples + 1][3];
        double[][] thetaTrue = new double[samples][3];

        thetaHat[0][0] = 0.80;
        thetaHat[0][1] = 0.35;
        thetaHat[0][2] = 0.0;

        for (int k = 0; k < samples; k++) {
            time[k] = k * sampleTime;
            double t = time[k];
            reference[k] = referenceCommand(t);

            MotorParameters motor = t < changeTime ? motorBefore : motorAfter;
            DiscretePlant plant = zohParameters(motor, sampleTime);
            thetaTrue[k][0] = plant.a;
            thetaTrue[k][1] = plant.b;
            thetaTrue[k][2] = plant.d;

            double noisyVelocity = velocity[k] + measurementNoiseStd * random.nextGaussian();
            measuredVelocity[k] = encoderQuantum * Math.rint(noisyVelocity / encoderQuantum);

            double aHat = thetaHat[k][0];
            double bHat = thetaHat[k][1];
            double dHat = thetaHat[k][2];
            double safeB = Math.max(bHat, bBounds[0]);

            double unsaturatedControl =
                    ((aModel - aHat) * measuredVelocity[k] + bModel * reference[k] - dHat) / safeB;
            control[k] = clip(unsaturatedControl, -inputLimit, inputLimit);

            velocity[k + 1] = plant.a * velocity[k]
                    + plant.b * control[k]
                    + plant.d
                    + processNoiseStd * random.nextGaussian();
            modelVelocity[k + 1] = aModel * modelVelocity[k] + bModel * reference[k];

            double noisyNext = velocity[k + 1] + measurementNoiseStd * random.nextGaussian();
            measuredVelocity[k + 1] = encoderQuantum * Math.rint(noisyNext / encoderQuantum);

            double[] psi = {measuredVelocity[k], control[k], 1.0};
            double prediction = thetaHat[k][0] * psi[0]
                    + thetaHat[k][1] * psi[1]
                    + thetaHat[k][2] * psi[2];
            predictionError[k] = measuredVelocity[k + 1] - prediction;

            thetaHat[k + 1][0] = thetaHat[k][0];
            thetaHat[k + 1][1] = thetaHat[k][1];
            thetaHat[k + 1][2] = thetaHat[k][2];

            if (Math.abs(predictionError[k]) > errorDeadZone) {
                double normalizer = normalizationOffset
                        + psi[0] * psi[0]
                        + psi[1] * psi[1]
                        + psi[2] * psi[2];
                for (int i = 0; i < 3; i++) {
                    thetaHat[k + 1][i] += gamma * psi[i] * predictionError[k] / normalizer;
                }
            }

            thetaHat[k + 1][0] = clip(thetaHat[k + 1][0], aBounds[0], aBounds[1]);
            thetaHat[k + 1][1] = clip(thetaHat[k + 1][1], bBounds[0], bBounds[1]);
            thetaHat[k + 1][2] = clip(thetaHat[k + 1][2], dBounds[0], dBounds[1]);
        }
        time[samples] = finalTime;

        int startIndex = (int) (2.0 / sampleTime);
        double squaredErrorSum = 0.0;
        for (int k = startIndex; k <= samples; k++) {
            double error = velocity[k] - modelVelocity[k];
            squaredErrorSum += error * error;
        }
        double rmsError = Math.sqrt(squaredErrorSum / (samples + 1 - startIndex));

        System.out.printf("RMS tracking error after 2 s: %.6f rad/s%n", rmsError);
        System.out.printf("Final estimate [a_hat, b_hat, d_hat]: [%.6f, %.6f, %.6f]%n",
                thetaHat[samples][0], thetaHat[samples][1], thetaHat[samples][2]);
        System.out.printf("Final true parameters [a, b, d]: [%.6f, %.6f, %.6f]%n",
                thetaTrue[samples - 1][0], thetaTrue[samples - 1][1], thetaTrue[samples - 1][2]);

        Path output = Path.of("Chapter14_Lesson5_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("time_s,reference,velocity,reference_model,control,a_hat,b_hat,d_hat,prediction_error\n");
            for (int k = 0; k < samples; k++) {
                writer.write(String.format(Locale.US,
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        time[k], reference[k], velocity[k], modelVelocity[k], control[k],
                        thetaHat[k][0], thetaHat[k][1], thetaHat[k][2], predictionError[k]));
            }
            writer.write(String.format(Locale.US,
                    "%.12f,,%.12f,%.12f,,%.12f,%.12f,%.12f,%n",
                    time[samples], velocity[samples], modelVelocity[samples],
                    thetaHat[samples][0], thetaHat[samples][1], thetaHat[samples][2]));
        } catch (IOException error) {
            System.err.println("Could not write CSV: " + error.getMessage());
            System.exit(1);
        }
        System.out.println("Saved: " + output.toAbsolutePath());
    }
}

The MATLAB function uses base numerical and plotting facilities and writes a table with writetable. Run it by placing the file on the MATLAB path and entering Chapter14_Lesson5.

Chapter14_Lesson5.m

function Chapter14_Lesson5
% Chapter14_Lesson5.m
% Discrete-time indirect adaptive velocity control using normalized gradient.
% Plant: y(k+1) = a*y(k) + b*u(k) + d.

rng(14, 'twister');

Ts = 0.02;
Tf = 30.0;
N = round(Tf / Ts);

tauModel = 0.20;
am = exp(-Ts / tauModel);
bm = 1 - am;

gamma = 0.55;              % theorem: 0 < gamma < 2
delta = 1.0;               % normalization offset
errorDeadZone = 0.015;

aBounds = [0.0, 0.9995];
bBounds = [0.05, 1.50];
dBounds = [-1.50, 1.50];

uLimit = 5.0;
encoderQuantum = 0.01;
measurementNoiseStd = 0.004;
processNoiseStd = 0.001;

motorBefore = struct('J', 0.020, 'B', 0.080, 'Kt', 0.500, 'tauL', 0.040);
motorAfter  = struct('J', 0.032, 'B', 0.100, 'Kt', 0.450, 'tauL', 0.100);
changeTime = 15.0;

time = (0:N)' * Ts;
r = zeros(N, 1);
y = zeros(N + 1, 1);
yMeasured = zeros(N + 1, 1);
yModel = zeros(N + 1, 1);
u = zeros(N, 1);
predictionError = zeros(N, 1);
thetaHat = zeros(N + 1, 3);
thetaTrue = zeros(N, 3);

thetaHat(1, :) = [0.80, 0.35, 0.0];

for k = 1:N
    t = time(k);
    r(k) = referenceCommand(t);

    if t < changeTime
        motor = motorBefore;
    else
        motor = motorAfter;
    end
    [a, b, d] = zohParameters(motor, Ts);
    thetaTrue(k, :) = [a, b, d];

    noisyVelocity = y(k) + measurementNoiseStd * randn;
    yMeasured(k) = encoderQuantum * round(noisyVelocity / encoderQuantum);

    aHat = thetaHat(k, 1);
    bHat = thetaHat(k, 2);
    dHat = thetaHat(k, 3);
    bSafe = max(bHat, bBounds(1));

    uUnsaturated = ((am - aHat) * yMeasured(k) + bm * r(k) - dHat) / bSafe;
    u(k) = min(max(uUnsaturated, -uLimit), uLimit);

    y(k + 1) = a * y(k) + b * u(k) + d + processNoiseStd * randn;
    yModel(k + 1) = am * yModel(k) + bm * r(k);

    noisyNext = y(k + 1) + measurementNoiseStd * randn;
    yMeasured(k + 1) = encoderQuantum * round(noisyNext / encoderQuantum);

    psi = [yMeasured(k); u(k); 1.0];
    predictionError(k) = yMeasured(k + 1) - thetaHat(k, :) * psi;

    thetaNext = thetaHat(k, :)';
    if abs(predictionError(k)) > errorDeadZone
        normalizer = delta + psi' * psi;
        thetaNext = thetaNext + gamma * psi * predictionError(k) / normalizer;
    end

    thetaNext(1) = min(max(thetaNext(1), aBounds(1)), aBounds(2));
    thetaNext(2) = min(max(thetaNext(2), bBounds(1)), bBounds(2));
    thetaNext(3) = min(max(thetaNext(3), dBounds(1)), dBounds(2));
    thetaHat(k + 1, :) = thetaNext';
end

startIndex = round(2.0 / Ts) + 1;
trackingError = y - yModel;
rmsError = sqrt(mean(trackingError(startIndex:end).^2));

fprintf('RMS tracking error after 2 s: %.6f rad/s\n', rmsError);
fprintf('Final estimate [aHat, bHat, dHat]: [%.6f, %.6f, %.6f]\n', thetaHat(end, :));
fprintf('Final true parameters [a, b, d]: [%.6f, %.6f, %.6f]\n', thetaTrue(end, :));

results = table(time(1:N), r, y(1:N), yModel(1:N), u, ...
    thetaHat(1:N, 1), thetaHat(1:N, 2), thetaHat(1:N, 3), predictionError, ...
    'VariableNames', {'time_s', 'reference', 'velocity', 'reference_model', ...
    'control', 'a_hat', 'b_hat', 'd_hat', 'prediction_error'});
writetable(results, 'Chapter14_Lesson5_results_matlab.csv');

figure('Name', 'Adaptive velocity tracking');
plot(time(1:N), r, 'LineWidth', 1.0); hold on;
plot(time, y, 'LineWidth', 1.2);
plot(time, yModel, '--', 'LineWidth', 1.2);
xlabel('Time (s)'); ylabel('Velocity (rad/s)'); grid on;
legend('Command', 'Velocity', 'Reference model', 'Location', 'best');

figure('Name', 'Control input');
plot(time(1:N), u, 'LineWidth', 1.0);
xlabel('Time (s)'); ylabel('Control input'); grid on;

figure('Name', 'Parameter estimates');
plot(time, thetaHat(:, 1), 'LineWidth', 1.0); hold on;
plot(time, thetaHat(:, 2), 'LineWidth', 1.0);
plot(time, thetaHat(:, 3), 'LineWidth', 1.0);
xlabel('Time (s)'); ylabel('Estimate'); grid on;
legend('aHat', 'bHat', 'dHat', 'Location', 'best');
end

function [a, b, d] = zohParameters(motor, Ts)
if motor.J <= 0 || Ts <= 0
    error('Inertia and sample time must be positive.');
end
if abs(motor.B) < 1e-12
    a = 1.0;
    b = motor.Kt * Ts / motor.J;
    d = -motor.tauL * Ts / motor.J;
else
    a = exp(-motor.B * Ts / motor.J);
    oneMinusA = 1 - a;
    b = (motor.Kt / motor.B) * oneMinusA;
    d = -(motor.tauL / motor.B) * oneMinusA;
end
end

function r = referenceCommand(t)
if t < 2.0
    r = 0.0;
elseif t < 8.0
    r = 6.0;
elseif t < 14.0
    r = 10.0;
elseif t < 20.0
    r = 4.0;
elseif t < 26.0
    r = 7.0 + 2.0 * sin(2*pi*0.40*(t - 20.0));
else
    tau = t - 26.0;
    r = 7.0 + 1.5 * sin(2*pi*0.65*tau) + 0.6 * sin(2*pi*1.35*tau);
end
end

The following MATLAB function creates a fixed-step discrete Simulink model programmatically. It places the sample-by-sample adaptive algorithm in a MATLAB Function block and routes velocity, command, control, and parameter estimates to workspace logging blocks. Simulink is required.

Chapter14_Lesson5_Simulink.m

function Chapter14_Lesson5_Simulink
% Chapter14_Lesson5_Simulink.m
% Programmatically builds a discrete Simulink demonstration of the adaptive
% velocity loop. Requires Simulink and a MATLAB Function block.

model = 'Chapter14_Lesson5_Model';
Ts = 0.02;
Tf = 30.0;

if ~license('test', 'Simulink')
    error('Simulink is required to build this model.');
end

if bdIsLoaded(model)
    close_system(model, 0);
end
if exist([model '.slx'], 'file')
    delete([model '.slx']);
end

new_system(model);
open_system(model);
set_param(model, 'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', 'FixedStep', num2str(Ts), ...
    'StopTime', num2str(Tf));

% Reference signal is generated in the MATLAB Function block from simulation time.
add_block('simulink/Sources/Clock', [model '/Clock'], ...
    'Position', [40 75 70 105]);
add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/AdaptiveVelocityLoop'], 'Position', [140 35 360 190]);

root = sfroot;
chart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [model '/AdaptiveVelocityLoop']);
if isempty(chart)
    error('Could not access the MATLAB Function block.');
end

chart.Script = sprintf([ ...
'function [r,y,ym,u,aHat,bHat,dHat] = fcn(t)\n' ...
'%%#codegen\n' ...
'persistent yk ymk th lastT\n' ...
'if isempty(yk) || t < lastT\n' ...
'    yk = 0; ymk = 0; th = [0.80;0.35;0.0]; lastT = -0.02;\n' ...
'end\n' ...
'Ts=0.02; am=exp(-Ts/0.20); bm=1-am; gamma=0.55; delta=1.0;\n' ...
'if t<2, r=0; elseif t<8, r=6; elseif t<14, r=10; elseif t<20, r=4;\n' ...
'elseif t<26, r=7+2*sin(2*pi*0.40*(t-20));\n' ...
'else, q=t-26; r=7+1.5*sin(2*pi*0.65*q)+0.6*sin(2*pi*1.35*q); end\n' ...
'if t<15, J=0.020; B=0.080; Kt=0.500; tauL=0.040;\n' ...
'else, J=0.032; B=0.100; Kt=0.450; tauL=0.100; end\n' ...
'a=exp(-B*Ts/J); b=(Kt/B)*(1-a); d=-(tauL/B)*(1-a);\n' ...
'bSafe=max(th(2),0.05); uRaw=((am-th(1))*yk+bm*r-th(3))/bSafe;\n' ...
'u=min(max(uRaw,-5),5); yNext=a*yk+b*u+d; ymNext=am*ymk+bm*r;\n' ...
'psi=[yk;u;1]; eps=yNext-th''*psi;\n' ...
'if abs(eps)>0.015, th=th+gamma*psi*eps/(delta+psi''*psi); end\n' ...
'th(1)=min(max(th(1),0),0.9995); th(2)=min(max(th(2),0.05),1.50);\n' ...
'th(3)=min(max(th(3),-1.50),1.50);\n' ...
'y=yk; ym=ymk; aHat=th(1); bHat=th(2); dHat=th(3);\n' ...
'yk=yNext; ymk=ymNext; lastT=t;\n' ...
'end\n']);

% Update once so the MATLAB Function block exposes all seven output ports.
set_param(model, 'SimulationCommand', 'update');

outputNames = {'reference', 'velocity', 'model_velocity', 'control', ...
    'a_hat', 'b_hat', 'd_hat'};
for i = 1:numel(outputNames)
    y0 = 20 + 45 * i;
    add_block('simulink/Sinks/To Workspace', ...
        [model '/' outputNames{i}], ...
        'VariableName', outputNames{i}, ...
        'SaveFormat', 'Structure With Time', ...
        'Position', [450 y0 555 y0 + 25]);
end

add_line(model, 'Clock/1', 'AdaptiveVelocityLoop/1', 'autorouting', 'on');
for i = 1:numel(outputNames)
    add_line(model, sprintf('AdaptiveVelocityLoop/%d', i), ...
        sprintf('%s/1', outputNames{i}), 'autorouting', 'on');
end

save_system(model);
set_param(model, 'SimulationCommand', 'update');
fprintf('Created %s.slx\n', model);
fprintf('Run: sim(''%s'')\n', model);
end

The readable Wolfram Language implementation is stored in Chapter14_Lesson5.wl. The companion notebook loads that file, keeping the notebook small while preserving a directly executable implementation.

Chapter14_Lesson5.wl

(* Chapter14_Lesson5.wl
   Discrete-time indirect adaptive velocity control using normalized gradient. *)

ClearAll["Global`*"];
SeedRandom[14];

ts = 0.02;
tf = 30.0;
n = Round[tf/ts];

tauModel = 0.20;
am = Exp[-ts/tauModel];
bm = 1 - am;

gamma = 0.55;
delta = 1.0;
errorDeadZone = 0.015;

aBounds = {0.0, 0.9995};
bBounds = {0.05, 1.50};
dBounds = {-1.50, 1.50};

uLimit = 5.0;
encoderQuantum = 0.01;
measurementNoiseStd = 0.004;
processNoiseStd = 0.001;

motorBefore = <|"J" -> 0.020, "B" -> 0.080, "Kt" -> 0.500, "TauL" -> 0.040|>;
motorAfter = <|"J" -> 0.032, "B" -> 0.100, "Kt" -> 0.450, "TauL" -> 0.100|>;
changeTime = 15.0;

zohParameters[motor_Association] := Module[{a, oneMinusA, b, d},
  If[motor["J"] <= 0 || ts <= 0, Message[zohParameters::bad]; Abort[]];
  If[Abs[motor["B"]] < 10^-12,
    {1.0, motor["Kt"] ts/motor["J"], -motor["TauL"] ts/motor["J"]},
    a = Exp[-motor["B"] ts/motor["J"]];
    oneMinusA = 1 - a;
    b = (motor["Kt"]/motor["B"]) oneMinusA;
    d = -(motor["TauL"]/motor["B"]) oneMinusA;
    {a, b, d}
  ]
];
zohParameters::bad = "Inertia and sample time must be positive.";

referenceCommand[t_] := Piecewise[{
   {0.0, t < 2.0},
   {6.0, t < 8.0},
   {10.0, t < 14.0},
   {4.0, t < 20.0},
   {7.0 + 2.0 Sin[2 Pi 0.40 (t - 20.0)], t < 26.0}
   }, 7.0 + 1.5 Sin[2 Pi 0.65 (t - 26.0)] + 0.6 Sin[2 Pi 1.35 (t - 26.0)]];

time = N[Range[0, n] ts];
reference = ConstantArray[0.0, n];
y = ConstantArray[0.0, n + 1];
yMeasured = ConstantArray[0.0, n + 1];
yModel = ConstantArray[0.0, n + 1];
u = ConstantArray[0.0, n];
predictionError = ConstantArray[0.0, n];
thetaHat = ConstantArray[{0.0, 0.0, 0.0}, n + 1];
thetaTrue = ConstantArray[{0.0, 0.0, 0.0}, n];
thetaHat[[1]] = {0.80, 0.35, 0.0};

Do[
  t = time[[k]];
  reference[[k]] = referenceCommand[t];
  motor = If[t < changeTime, motorBefore, motorAfter];
  {a, b, d} = zohParameters[motor];
  thetaTrue[[k]] = {a, b, d};

  noisyVelocity = y[[k]] + RandomVariate[NormalDistribution[0, measurementNoiseStd]];
  yMeasured[[k]] = encoderQuantum Round[noisyVelocity/encoderQuantum];

  {aHat, bHat, dHat} = thetaHat[[k]];
  bSafe = Max[bHat, bBounds[[1]]];
  uUnsaturated = ((am - aHat) yMeasured[[k]] + bm reference[[k]] - dHat)/bSafe;
  u[[k]] = Clip[uUnsaturated, {-uLimit, uLimit}];

  y[[k + 1]] = a y[[k]] + b u[[k]] + d +
    RandomVariate[NormalDistribution[0, processNoiseStd]];
  yModel[[k + 1]] = am yModel[[k]] + bm reference[[k]];

  noisyNext = y[[k + 1]] + RandomVariate[NormalDistribution[0, measurementNoiseStd]];
  yMeasured[[k + 1]] = encoderQuantum Round[noisyNext/encoderQuantum];

  psi = {yMeasured[[k]], u[[k]], 1.0};
  predictionError[[k]] = yMeasured[[k + 1]] - thetaHat[[k]].psi;

  thetaNext = thetaHat[[k]];
  If[Abs[predictionError[[k]]] > errorDeadZone,
    normalizer = delta + psi.psi;
    thetaNext = thetaNext + gamma psi predictionError[[k]]/normalizer;
  ];
  thetaHat[[k + 1]] = {
    Clip[thetaNext[[1]], aBounds],
    Clip[thetaNext[[2]], bBounds],
    Clip[thetaNext[[3]], dBounds]
  };
, {k, 1, n}];

startIndex = Round[2.0/ts] + 1;
trackingError = y - yModel;
rmsError = Sqrt[Mean[trackingError[[startIndex ;;]]^2]];

Print["RMS tracking error after 2 s: ", NumberForm[rmsError, {8, 6}], " rad/s"];
Print["Final estimate [aHat, bHat, dHat]: ", thetaHat[[-1]]];
Print["Final true parameters [a, b, d]: ", thetaTrue[[-1]]];

rows = Prepend[
   Table[{time[[k]], reference[[k]], y[[k]], yModel[[k]], u[[k]],
     thetaHat[[k, 1]], thetaHat[[k, 2]], thetaHat[[k, 3]], predictionError[[k]]},
    {k, 1, n}],
   {"time_s", "reference", "velocity", "reference_model", "control",
    "a_hat", "b_hat", "d_hat", "prediction_error"}
];
Export["Chapter14_Lesson5_results_wolfram.csv", rows];

trackingPlot = ListLinePlot[
  {Transpose[{time[[1 ;; n]], reference}], Transpose[{time, y}], Transpose[{time, yModel}]},
  PlotLegends -> {"Command", "Velocity", "Reference model"},
  Frame -> True, FrameLabel -> {"Time (s)", "Velocity (rad/s)"},
  PlotRange -> All, ImageSize -> Large
];

controlPlot = ListLinePlot[Transpose[{time[[1 ;; n]], u}],
  Frame -> True, FrameLabel -> {"Time (s)", "Control input"},
  PlotRange -> All, ImageSize -> Large];

parameterPlot = ListLinePlot[
  Table[Transpose[{time, thetaHat[[All, i]]}], {i, 1, 3}],
  PlotLegends -> {"aHat", "bHat", "dHat"},
  Frame -> True, FrameLabel -> {"Time (s)", "Estimate"},
  PlotRange -> All, ImageSize -> Large
];

Column[{trackingPlot, controlPlot, parameterPlot}]

Chapter14_Lesson5.nb

Notebook[{
  Cell["Chapter 14, Lesson 5: Discrete-Time Adaptive Control for a Velocity Loop", "Title"],
  Cell["Keep Chapter14_Lesson5.wl in the same folder, then evaluate the input cell.", "Text"],
  Cell[BoxData[RowBox[{"Get", "[", "\"Chapter14_Lesson5.wl\"", "]"}]], "Input"]
}, WindowTitle -> "Chapter14_Lesson5", Saveable -> True]

Problem 1 (Exact discretization): A velocity plant has \(J=0.025\), \(B=0.10\), \(K_t=0.40\), \(\tau_L=0.05\), and \(T_s=0.01\text{ s}\). Derive the exact discrete model.

Solution:

\[ a=\exp\!\left(-\frac{0.10(0.01)}{0.025}\right) =e^{-0.04}\approx0.960789. \]

\[ b=\frac{0.40}{0.10}(1-a) =4(0.039211)\approx0.156842. \]

\[ d=-\frac{0.05}{0.10}(1-a) =-0.5(0.039211)\approx-0.019605. \]

Thus

\[ y_{k+1}=0.960789y_k+0.156842u_k-0.019605. \]

Problem 2 (Ideal model matching): For the plant in Problem 1, choose \(a_m=0.90\) and \(b_m=0.10\). Derive the ideal unsaturated controller.

Solution:

\[ u_k^{\star}=\frac{(0.90-0.960789)y_k+0.10r_k-(-0.019605)}{0.156842}. \]

\[ u_k^{\star}\approx-0.38759y_k+0.63758r_k+0.12500. \]

Substitution into the plant produces \(y_{k+1}=0.90y_k+0.10r_k\), so the tracking error obeys \(e_{k+1}=0.90e_k\) when the reference model starts from the corresponding initial condition.

Problem 3 (Normalized-gradient bound): Starting from the unprojected update, prove that

\[ \sum_{k=0}^{\infty} \frac{\varepsilon_{k+1}^2}{m_k}<\infty \]

for \(0<\gamma<2\).

Solution: Section 5 established

\[ V_{k+1}-V_k\leq-\gamma(2-\gamma) \frac{\varepsilon_{k+1}^2}{m_k}. \]

Summing from \(0\) to \(N\) gives

\[ \gamma(2-\gamma)\sum_{k=0}^{N} \frac{\varepsilon_{k+1}^2}{m_k} \leq V_0-V_{N+1}\leq V_0. \]

Because \(V_{N+1}\geq0\) and the partial sums are monotone, their infinite limit exists and is no larger than \(V_0/[\gamma(2-\gamma)]\).

Problem 4 (One-sample delay): Suppose the actuator applies \(u_{k-1}\) during interval \(k\). What estimator regressor should be used, and why is \([y_k,u_k,1]^T\) incorrect?

Solution: The sampled plant is

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

Therefore the correct regressor is \(\boldsymbol{\psi}_k=[y_k,u_{k-1},1]^T\). Using \(u_k\) pairs the output with an input that has not yet affected it. The resulting regression error is correlated with the omitted delayed input and generally biases the estimate.

Problem 5 (Tracking without identification): A simulation shows small tracking error, but \(\hat b_k\) remains far from the true \(b\). Is the adaptive controller failing? Give a rigorous explanation and propose a diagnostic experiment.

Solution:

Not necessarily. Small tracking error requires the combined prediction mismatch entering the stable error filter to be small; it does not require every component of the parameter error to vanish. If the regressor is not persistently exciting, parameter-error components in unexcited directions are unobservable from the collected trajectory. To diagnose this, compute the minimum eigenvalue of a moving-window information matrix

\[ \mathbf{G}_k=\sum_{i=k}^{k+N_e-1} \boldsymbol{\psi}_i\boldsymbol{\psi}_i^T. \]

If its smallest eigenvalue remains near zero, add a safe multilevel or multisine reference segment, repeat the experiment, and compare parameter convergence. Any probing signal must remain within actuator, velocity, and safety constraints.

A sampled velocity loop was modeled as a first-order affine regression, identified online with a normalized-gradient law, and controlled by certainty-equivalent one-step model matching. The parameter-error calculation proved bounded estimates and summability of normalized prediction error for constant nominal parameters and \(0<\gamma<2\). Projection preserves this property while enforcing the known positive input direction. The tracking error is a stable filtered version of prediction mismatch, so tracking can converge even when individual physical parameters do not. Persistent excitation is needed for parameter convergence. Practical implementation additionally requires correct sample indexing, use of the actual saturated input, dead-zone logic, quantization awareness, delay modeling, and a safe fallback strategy.