After completing this lesson, the student should be able to:

• formulate a basic direct STR in which controller gains, rather than plant parameters, are the adaptive variables;
• derive the ideal controller gains from a reference-model matching condition;
• prove boundedness and asymptotic tracking for the ideal, noise-free, unsaturated closed loop;
• translate the continuous-time adaptive law into a sampled numerical implementation;
• evaluate nominal operation, abrupt plant variation, measurement noise, projection, and actuator saturation; and
• distinguish tracking convergence from controller-parameter convergence.

The implementation is intentionally based on a transparent scalar plant. This isolates the essential direct-STR mechanism before the course proceeds to the more general discrete-time and output-feedback developments of later chapters.

A direct self-tuning regulator updates the controller parameters themselves. It does not first construct estimates \( \hat a \) and \( \hat b \) of physical plant coefficients and then redesign a controller from those estimates. The adjustable variables in this lesson are the feedback gain \( k_y(t) \) and the command gain \( k_r(t) \).

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> YM["Desired output ym"]
  R --> C["Controller u = ky*y + kr*r"]
  Y["Measured output y"] --> C
  C --> P["Unknown plant"]
  P --> Y
  Y --> E["Tracking error e = y - ym"]
  YM --> E
  E --> A["Direct adaptation law"]
  Y --> A
  R --> A
  A --> G["Updated gains ky and kr"]
  G --> C
        

The key signal path is therefore \( e,y,r \mapsto \dot k_y,\dot k_r \). The adaptation mechanism is driven by the tracking objective, not by a separate prediction-error identification objective.

Consider the uncertain first-order plant

\[ \dot y(t)=-a(t)y(t)+b(t)u(t), \qquad b(t)>0. \]

The controller knows the sign of the high-frequency gain \( b(t) \), but it does not know its magnitude or the value of \( a(t) \). For the theoretical derivation, first suppose that \( a \) and \( b \) are unknown constants. The desired closed-loop response is specified by the stable reference model

\[ \dot y_m(t)=-a_m y_m(t)+b_m r(t), \qquad a_m>0. \]

The direct controller is parameterized as

\[ u(t)=k_y(t)y(t)+k_r(t)r(t). \]

This two-gain parameterization is minimal for matching the state coefficient and the command coefficient of the first-order reference model.

If the adjustable gains were equal to fixed ideal values \( k_y^\ast \) and \( k_r^\ast \), the plant would satisfy

\[ \dot y = \left(-a+b k_y^\ast\right)y + b k_r^\ast r. \]

Exact model matching requires coefficient equality:

\[ -a+b k_y^\ast=-a_m, \qquad b k_r^\ast=b_m. \]

Consequently,

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

These equations are used only to prove that suitable ideal controller gains exist. The online algorithm never evaluates them because \( a \) and \( b \) are not estimated. This is the defining direct-adaptation feature.

Define tracking and gain errors by

\[ e=y-y_m, \qquad \tilde k_y=k_y-k_y^\ast, \qquad \tilde k_r=k_r-k_r^\ast. \]

Substituting the controller into the plant, subtracting the reference model, and using the matching equations gives

\[ \begin{aligned} \dot e &= -a y+b k_y y+b k_r r+a_m y_m-b_m r \\ &= -a_m e+b\tilde k_y y+b\tilde k_r r. \end{aligned} \]

With the regressor and gain-error vector

\[ \boldsymbol{\phi} = \begin{bmatrix} y \\ r \end{bmatrix}, \qquad \tilde{\boldsymbol{k} } = \begin{bmatrix} \tilde k_y \\ \tilde k_r \end{bmatrix}, \]

the error system becomes

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

Choose positive adaptation gains \( \gamma_y \) and \( \gamma_r \). The baseline direct adaptive law is

\[ \dot k_y=-\gamma_y e y, \qquad \dot k_r=-\gamma_r e r. \]

In vector form, with \( \boldsymbol{\Gamma}= \operatorname{diag}(\gamma_y,\gamma_r) \),

\[ \dot{\boldsymbol{k} } = -\boldsymbol{\Gamma}\boldsymbol{\phi}e. \]

Assume that \( a \) and \( b \) are constant, \( b>0 \), the reference input is bounded, and neither actuator saturation nor measurement noise is present. Consider

\[ V = \frac{1}{2}e^2 + \frac{b}{2\gamma_y}\tilde k_y^2 + \frac{b}{2\gamma_r}\tilde k_r^2. \]

Since \( b>0 \) and the adaptation gains are positive, \( V \) is positive definite in \( e,\tilde k_y,\tilde k_r \). Because the ideal gains are constant, \( \dot{\tilde k}_y=\dot k_y \) and \( \dot{\tilde k}_r=\dot k_r \). Therefore,

\[ \begin{aligned} \dot V &= e\dot e + \frac{b}{\gamma_y}\tilde k_y\dot k_y + \frac{b}{\gamma_r}\tilde k_r\dot k_r \\ &= e\left(-a_m e+b\tilde k_y y+b\tilde k_r r\right) - b\tilde k_y e y - b\tilde k_r e r \\ &= -a_m e^2 \leq 0. \end{aligned} \]

Hence \( V(t)\leq V(0) \). It follows that \( e,\tilde k_y,\tilde k_r \) are bounded. A bounded command and a stable reference model imply bounded \( y_m \); because \( y=e+y_m \), the plant output is also bounded. Consequently, the regressor, control input, and \( \dot e \) are bounded.

Integrating the Lyapunov derivative yields

\[ a_m\int_{0}^{T}e^2(\tau)\,d\tau = V(0)-V(T) \leq V(0). \]

Thus \( e\in L_2\cap L_\infty \). Since \( \dot e\in L_\infty \), Barbalat's lemma gives

\[ \lim_{t→\infty}e(t)=0. \]

This proof establishes asymptotic tracking, but not necessarily \( k_y→k_y^\ast \) and \( k_r→k_r^\ast \). Gain convergence requires sufficient excitation of the regressor, as developed in Chapter 10.

The programs use a small integration interval \( h \) to realize the continuous-time equations. This is a numerical implementation of the continuous direct law, not yet a general discrete-time adaptive-control derivation.

\[ \begin{aligned} y_{k+1} &= y_k+h\left(-a_k y_k+b_k u_k\right), \\ y_{m,k+1} &= y_{m,k}+h\left(-a_m y_{m,k}+b_m r_k\right), \\ k_{y,k+1} &= k_{y,k}-h\gamma_y e_k y_k, \\ k_{r,k+1} &= k_{r,k}-h\gamma_r e_k r_k. \end{aligned} \]

The controller is evaluated from the current sampled values:

\[ u_k=k_{y,k}y_k+k_{r,k}r_k. \]

flowchart TD
  S["Initialize y, ym, ky, kr, step h"] --> M["Read command r and measured output y"]
  M --> E["Compute e = y - ym"]
  E --> U["Compute u = ky*y + kr*r"]
  U --> SAT["Apply actuator limit"]
  SAT --> A["Compute ky_dot = -gamma_y*e*y and kr_dot = -gamma_r*e*r"]
  A --> PR["Integrate gains and apply projection"]
  PR --> P["Integrate plant and reference model"]
  P --> L["Log y, ym, e, u, ky, kr"]
  L --> Q{"Final time reached?"}
  Q -->|No| M
  Q -->|Yes| R["Compute performance metrics and plots"]
        

8.1 Projection of the Direct Gains

The implementations constrain each gain to the interval \( [-k_{\max},k_{\max}] \). In scalar form,

\[ k_{i,k+1} = \operatorname{clip} \left( k_{i,k}+h\dot k_{i,k}, -k_{\max}, k_{\max} \right). \]

If the ideal gain lies strictly inside the projection set, a standard projection operator can preserve the Lyapunov non-increase property. The simple clipping implementation is a numerical approximation to that operator.

8.2 Actuator Saturation

The practical controller uses

\[ u_{\mathrm{sat} } = \operatorname{clip} \left( k_y y+k_r r, -u_{\max}, u_{\max} \right). \]

Saturation destroys the exact error-model equality used in the Lyapunov cancellation whenever the unsaturated command exceeds the actuator limit. Therefore, the proof of Section 6 applies only while the saturation is inactive. A complete treatment of saturation-aware adaptive laws appears in Chapter 19.

8.3 Measurement Noise

With a measured output \( y_{\mathrm{m} }=y+n \), both the controller and the gain update receive noise. The implemented update becomes

\[ \dot k_y = -\gamma_y \left(y_{\mathrm{m} }-y_m\right)y_{\mathrm{m} }, \qquad \dot k_r = -\gamma_r \left(y_{\mathrm{m} }-y_m\right)r. \]

Large adaptation gains may then convert high-frequency measurement noise into gain motion and control chatter. The robust modifications, dead zones, and normalization techniques studied in Chapters 8 and 9 are appropriate extensions when noise is substantial.

The common design values are

\[ a_m=2.5, \quad b_m=2.5, \quad \gamma_y=12, \quad \gamma_r=8, \quad h=10^{-3}\;\mathrm{s}. \]

The command is piecewise constant:

\[ r(t)= \begin{cases} 1, & 0\leq t<5, \\ -0.5, & 5\leq t<10, \\ 0.8, & 10\leq t<18, \\ 0.2, & 18\leq t<24, \\ -0.8, & 24\leq t\leq 30. \end{cases} \]

9.1 Case A: Constant Unknown Plant

\[ a(t)=1.2, \qquad b(t)=1.0. \]

The ideal matching gains for interpretation are

\[ k_y^\ast=-1.3, \qquad k_r^\ast=2.5. \]

9.2 Case B: Abrupt Plant Change

At \( t=15\;\mathrm{s} \), the physical plant changes:

\[ (a,b) = \begin{cases} (1.2,1.0), & t<15, \\ (2.0,0.65), & t\geq 15. \end{cases} \]

The corresponding ideal gains change from \( (-1.3,2.5) \) to

\[ k_y^\ast = \frac{2.0-2.5}{0.65} \approx -0.7692, \qquad k_r^\ast = \frac{2.5}{0.65} \approx 3.8462. \]

Because the baseline theorem assumes constant ideal gains, the abrupt change is a robustness and tracking experiment rather than a direct consequence of the constant-parameter Lyapunov proof.

Tracking and control behavior are evaluated with several complementary indices. The root-mean-square tracking error is

\[ \operatorname{RMSE} = \sqrt{ \frac{1}{N} \sum_{k=0}^{N-1}e_k^2 }. \]

The integral of absolute error is approximated by

\[ \operatorname{IAE} \approx h\sum_{k=0}^{N-1}|e_k|. \]

The RMS control effort is

\[ u_{\mathrm{RMS} } = \sqrt{ \frac{1}{N} \sum_{k=0}^{N-1}u_k^2 }. \]

A simple measure of control activity is the discrete total variation

\[ \operatorname{TV}(u) = \sum_{k=1}^{N-1}|u_k-u_{k-1}|. \]

Low tracking error alone is not sufficient. Excessive \( \operatorname{TV}(u) \), repeated saturation, or gains resting on their projection limits indicate an implementation problem even when the output visually follows the command.

This implementation uses NumPy for numerical arrays and Matplotlib for plotting. The adaptive law itself is coded from first principles; no specialized adaptive-control package is required.

Chapter13_Lesson5.py

"""
Chapter13_Lesson5.py
Basic direct self-tuning regulator (direct STR) case study.

The controller gains are adapted directly; no online estimate of the plant
parameters a(t) and b(t) is formed.

Dependencies:
    numpy
    matplotlib
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Dict
import numpy as np
import matplotlib.pyplot as plt


@dataclass
class Config:
    dt: float = 0.001
    final_time: float = 30.0
    a_model: float = 2.5
    b_model: float = 2.5
    gamma_y: float = 12.0
    gamma_r: float = 8.0
    u_limit: float = 6.0
    gain_limit: float = 12.0
    noise_std: float = 0.01
    seed: int = 7


def reference(t: float) -> float:
    """Bounded, piecewise-constant command with several excitation levels."""
    if t < 5.0:
        return 1.0
    if t < 10.0:
        return -0.5
    if t < 18.0:
        return 0.8
    if t < 24.0:
        return 0.2
    return -0.8


def plant_parameters(t: float, parameter_jump: bool) -> tuple[float, float]:
    """Unknown-to-controller plant coefficients."""
    if parameter_jump and t >= 15.0:
        return 2.0, 0.65
    return 1.2, 1.0


def project(value: float, bound: float) -> float:
    return float(np.clip(value, -bound, bound))


def simulate(config: Config, parameter_jump: bool) -> Dict[str, np.ndarray]:
    rng = np.random.default_rng(config.seed + int(parameter_jump))
    steps = int(round(config.final_time / config.dt)) + 1
    t = np.linspace(0.0, config.final_time, steps)

    y = np.zeros(steps)
    y_model = np.zeros(steps)
    error = np.zeros(steps)
    control = np.zeros(steps)
    k_y = np.zeros(steps)
    k_r = np.zeros(steps)
    command = np.zeros(steps)
    a_hist = np.zeros(steps)
    b_hist = np.zeros(steps)

    k_y[0] = 0.0
    k_r[0] = 0.5

    for i in range(steps - 1):
        ti = t[i]
        command[i] = reference(ti)
        a_i, b_i = plant_parameters(ti, parameter_jump)
        a_hist[i], b_hist[i] = a_i, b_i

        measured_y = y[i] + rng.normal(0.0, config.noise_std)
        error[i] = measured_y - y_model[i]

        # Direct controller parameterization: u = k_y*y + k_r*r.
        unsaturated_u = k_y[i] * measured_y + k_r[i] * command[i]
        control[i] = float(np.clip(unsaturated_u, -config.u_limit, config.u_limit))

        # Baseline direct update from the Lyapunov derivation.
        # Projection prevents numerical gain escape if a bound is reached.
        k_y_dot = -config.gamma_y * error[i] * measured_y
        k_r_dot = -config.gamma_r * error[i] * command[i]

        k_y[i + 1] = project(k_y[i] + config.dt * k_y_dot, config.gain_limit)
        k_r[i + 1] = project(k_r[i] + config.dt * k_r_dot, config.gain_limit)

        # Forward-Euler realization of the continuous-time plant/reference model.
        y_dot = -a_i * y[i] + b_i * control[i]
        y_model_dot = -config.a_model * y_model[i] + config.b_model * command[i]
        y[i + 1] = y[i] + config.dt * y_dot
        y_model[i + 1] = y_model[i] + config.dt * y_model_dot

    command[-1] = reference(t[-1])
    a_hist[-1], b_hist[-1] = plant_parameters(t[-1], parameter_jump)
    error[-1] = y[-1] - y_model[-1]
    control[-1] = np.clip(k_y[-1] * y[-1] + k_r[-1] * command[-1],
                          -config.u_limit, config.u_limit)

    return {
        "t": t,
        "r": command,
        "y": y,
        "ym": y_model,
        "e": error,
        "u": control,
        "ky": k_y,
        "kr": k_r,
        "a": a_hist,
        "b": b_hist,
    }


def metrics(result: Dict[str, np.ndarray]) -> Dict[str, float]:
    t = result["t"]
    dt = t[1] - t[0]
    e = result["e"]
    u = result["u"]
    return {
        "rmse": float(np.sqrt(np.mean(e**2))),
        "iae": float(np.sum(np.abs(e)) * dt),
        "max_abs_error": float(np.max(np.abs(e))),
        "control_rms": float(np.sqrt(np.mean(u**2))),
        "control_total_variation": float(np.sum(np.abs(np.diff(u)))),
    }


def plot_case(result: Dict[str, np.ndarray], title: str, file_name: str) -> None:
    t = result["t"]

    plt.figure(figsize=(10, 7))
    plt.subplot(3, 1, 1)
    plt.plot(t, result["r"], "--", label="reference r")
    plt.plot(t, result["ym"], label="reference model ym")
    plt.plot(t, result["y"], label="plant output y")
    plt.ylabel("output")
    plt.legend()
    plt.grid(True)

    plt.subplot(3, 1, 2)
    plt.plot(t, result["u"], label="control u")
    plt.plot(t, result["e"], label="tracking error e")
    plt.ylabel("u, e")
    plt.legend()
    plt.grid(True)

    plt.subplot(3, 1, 3)
    plt.plot(t, result["ky"], label="k_y")
    plt.plot(t, result["kr"], label="k_r")
    plt.xlabel("time [s]")
    plt.ylabel("direct gains")
    plt.legend()
    plt.grid(True)

    plt.suptitle(title)
    plt.tight_layout()
    plt.savefig(file_name, dpi=160)
    plt.close()


def main() -> None:
    config = Config()
    nominal = simulate(config, parameter_jump=False)
    changed = simulate(config, parameter_jump=True)

    print("Nominal case metrics:")
    for key, value in metrics(nominal).items():
        print(f"  {key:26s}: {value:.6f}")

    print("\nAbrupt plant-change case metrics:")
    for key, value in metrics(changed).items():
        print(f"  {key:26s}: {value:.6f}")

    plot_case(
        nominal,
        "Direct STR: nominal plant",
        "Chapter13_Lesson5_nominal.png",
    )
    plot_case(
        changed,
        "Direct STR: plant parameters change at t = 15 s",
        "Chapter13_Lesson5_parameter_jump.png",
    )


if __name__ == "__main__":
    main()

The C++17 program uses only the standard library. It executes both simulation cases, prints the metrics, and writes the plant-change case to a CSV file for plotting in another environment.

Chapter13_Lesson5.cpp

/*
Chapter13_Lesson5.cpp
Basic direct self-tuning regulator (direct STR) case study.

Build:
    g++ -std=c++17 -O2 Chapter13_Lesson5.cpp -o Chapter13_Lesson5
Run:
    ./Chapter13_Lesson5

The program writes Chapter13_Lesson5_cpp_results.csv.
*/

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

struct Config {
    double dt = 0.001;
    double finalTime = 30.0;
    double aModel = 2.5;
    double bModel = 2.5;
    double gammaY = 12.0;
    double gammaR = 8.0;
    double uLimit = 6.0;
    double gainLimit = 12.0;
    double noiseStd = 0.01;
    unsigned int seed = 7;
};

double reference(double t) {
    if (t < 5.0) return 1.0;
    if (t < 10.0) return -0.5;
    if (t < 18.0) return 0.8;
    if (t < 24.0) return 0.2;
    return -0.8;
}

void plantParameters(double t, bool parameterJump, double& a, double& b) {
    if (parameterJump && t >= 15.0) {
        a = 2.0;
        b = 0.65;
    } else {
        a = 1.2;
        b = 1.0;
    }
}

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

struct Sample {
    double t{};
    double r{};
    double y{};
    double ym{};
    double e{};
    double u{};
    double ky{};
    double kr{};
    double a{};
    double b{};
};

std::vector<Sample> simulate(const Config& cfg, bool parameterJump) {
    const std::size_t steps =
        static_cast<std::size_t>(std::llround(cfg.finalTime / cfg.dt)) + 1U;
    std::vector<Sample> data(steps);

    double y = 0.0;
    double ym = 0.0;
    double ky = 0.0;
    double kr = 0.5;

    std::mt19937 generator(cfg.seed + static_cast<unsigned int>(parameterJump));
    std::normal_distribution<double> noise(0.0, cfg.noiseStd);

    for (std::size_t i = 0; i < steps; ++i) {
        const double t = static_cast<double>(i) * cfg.dt;
        const double r = reference(t);
        double a = 0.0;
        double b = 0.0;
        plantParameters(t, parameterJump, a, b);

        const double measuredY = y + noise(generator);
        const double e = measuredY - ym;
        const double rawU = ky * measuredY + kr * r;
        const double u = clip(rawU, -cfg.uLimit, cfg.uLimit);

        data[i] = {t, r, y, ym, e, u, ky, kr, a, b};

        if (i + 1U == steps) break;

        const double kyDot = -cfg.gammaY * e * measuredY;
        const double krDot = -cfg.gammaR * e * r;

        ky = clip(ky + cfg.dt * kyDot, -cfg.gainLimit, cfg.gainLimit);
        kr = clip(kr + cfg.dt * krDot, -cfg.gainLimit, cfg.gainLimit);

        const double yDot = -a * y + b * u;
        const double ymDot = -cfg.aModel * ym + cfg.bModel * r;
        y += cfg.dt * yDot;
        ym += cfg.dt * ymDot;
    }

    return data;
}

void reportMetrics(const std::vector<Sample>& data, const Config& cfg,
                   const std::string& label) {
    double squaredErrorSum = 0.0;
    double absoluteErrorSum = 0.0;
    double maxAbsError = 0.0;
    double squaredControlSum = 0.0;
    double totalVariation = 0.0;

    for (std::size_t i = 0; i < data.size(); ++i) {
        squaredErrorSum += data[i].e * data[i].e;
        absoluteErrorSum += std::abs(data[i].e) * cfg.dt;
        maxAbsError = std::max(maxAbsError, std::abs(data[i].e));
        squaredControlSum += data[i].u * data[i].u;
        if (i > 0U) totalVariation += std::abs(data[i].u - data[i - 1U].u);
    }

    std::cout << label << "\n"
              << "  RMSE: " << std::sqrt(squaredErrorSum / data.size()) << "\n"
              << "  IAE: " << absoluteErrorSum << "\n"
              << "  Maximum absolute error: " << maxAbsError << "\n"
              << "  Control RMS: " << std::sqrt(squaredControlSum / data.size()) << "\n"
              << "  Control total variation: " << totalVariation << "\n\n";
}

void writeCsv(const std::vector<Sample>& data, const std::string& fileName) {
    std::ofstream file(fileName);
    if (!file) throw std::runtime_error("Cannot open output CSV file.");

    file << "t,r,y,ym,e,u,ky,kr,a,b\n";
    file << std::setprecision(12);
    for (const auto& s : data) {
        file << s.t << ',' << s.r << ',' << s.y << ',' << s.ym << ','
             << s.e << ',' << s.u << ',' << s.ky << ',' << s.kr << ','
             << s.a << ',' << s.b << '\n';
    }
}

int main() {
    try {
        const Config cfg;
        const auto nominal = simulate(cfg, false);
        const auto changed = simulate(cfg, true);

        reportMetrics(nominal, cfg, "Nominal case metrics:");
        reportMetrics(changed, cfg, "Abrupt plant-change case metrics:");
        writeCsv(changed, "Chapter13_Lesson5_cpp_results.csv");
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses standard records, collections, random numbers, and file I/O. It follows the same equations and parameter values as the Python and C++ versions.

Chapter13_Lesson5.java

/*
Chapter13_Lesson5.java
Basic direct self-tuning regulator (direct STR) case study.

Build:
    javac Chapter13_Lesson5.java
Run:
    java Chapter13_Lesson5

The program writes Chapter13_Lesson5_java_results.csv.
*/

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.Random;

public final class Chapter13_Lesson5 {
    private Chapter13_Lesson5() {}

    private static final class Config {
        final double dt = 0.001;
        final double finalTime = 30.0;
        final double aModel = 2.5;
        final double bModel = 2.5;
        final double gammaY = 12.0;
        final double gammaR = 8.0;
        final double uLimit = 6.0;
        final double gainLimit = 12.0;
        final double noiseStd = 0.01;
        final long seed = 7L;
    }

    private record Sample(
        double t, double r, double y, double ym, double e,
        double u, double ky, double kr, double a, double b
    ) {}

    private static double reference(double t) {
        if (t < 5.0) return 1.0;
        if (t < 10.0) return -0.5;
        if (t < 18.0) return 0.8;
        if (t < 24.0) return 0.2;
        return -0.8;
    }

    private static double[] plantParameters(double t, boolean parameterJump) {
        if (parameterJump && t >= 15.0) {
            return new double[] {2.0, 0.65};
        }
        return new double[] {1.2, 1.0};
    }

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

    private static List<Sample> simulate(Config cfg, boolean parameterJump) {
        int steps = (int) Math.round(cfg.finalTime / cfg.dt) + 1;
        List<Sample> data = new ArrayList<>(steps);

        double y = 0.0;
        double ym = 0.0;
        double ky = 0.0;
        double kr = 0.5;
        Random random = new Random(cfg.seed + (parameterJump ? 1L : 0L));

        for (int i = 0; i < steps; i++) {
            double t = i * cfg.dt;
            double r = reference(t);
            double[] parameters = plantParameters(t, parameterJump);
            double a = parameters[0];
            double b = parameters[1];

            double measuredY = y + cfg.noiseStd * random.nextGaussian();
            double e = measuredY - ym;
            double rawU = ky * measuredY + kr * r;
            double u = clip(rawU, -cfg.uLimit, cfg.uLimit);

            data.add(new Sample(t, r, y, ym, e, u, ky, kr, a, b));
            if (i + 1 == steps) break;

            double kyDot = -cfg.gammaY * e * measuredY;
            double krDot = -cfg.gammaR * e * r;

            ky = clip(ky + cfg.dt * kyDot, -cfg.gainLimit, cfg.gainLimit);
            kr = clip(kr + cfg.dt * krDot, -cfg.gainLimit, cfg.gainLimit);

            double yDot = -a * y + b * u;
            double ymDot = -cfg.aModel * ym + cfg.bModel * r;
            y += cfg.dt * yDot;
            ym += cfg.dt * ymDot;
        }
        return data;
    }

    private static void reportMetrics(
        List<Sample> data, Config cfg, String label
    ) {
        double squaredErrorSum = 0.0;
        double absoluteErrorSum = 0.0;
        double maxAbsError = 0.0;
        double squaredControlSum = 0.0;
        double totalVariation = 0.0;

        for (int i = 0; i < data.size(); i++) {
            Sample sample = data.get(i);
            squaredErrorSum += sample.e() * sample.e();
            absoluteErrorSum += Math.abs(sample.e()) * cfg.dt;
            maxAbsError = Math.max(maxAbsError, Math.abs(sample.e()));
            squaredControlSum += sample.u() * sample.u();
            if (i > 0) {
                totalVariation += Math.abs(sample.u() - data.get(i - 1).u());
            }
        }

        System.out.printf(Locale.US, "%s%n", label);
        System.out.printf(Locale.US, "  RMSE: %.6f%n",
            Math.sqrt(squaredErrorSum / data.size()));
        System.out.printf(Locale.US, "  IAE: %.6f%n", absoluteErrorSum);
        System.out.printf(Locale.US, "  Maximum absolute error: %.6f%n", maxAbsError);
        System.out.printf(Locale.US, "  Control RMS: %.6f%n",
            Math.sqrt(squaredControlSum / data.size()));
        System.out.printf(Locale.US, "  Control total variation: %.6f%n%n",
            totalVariation);
    }

    private static void writeCsv(List<Sample> data, Path path)
        throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(path)) {
            writer.write("t,r,y,ym,e,u,ky,kr,a,b");
            writer.newLine();
            for (Sample s : data) {
                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    s.t(), s.r(), s.y(), s.ym(), s.e(),
                    s.u(), s.ky(), s.kr(), s.a(), s.b()
                ));
            }
        }
    }

    public static void main(String[] args) {
        Config cfg = new Config();
        List<Sample> nominal = simulate(cfg, false);
        List<Sample> changed = simulate(cfg, true);

        reportMetrics(nominal, cfg, "Nominal case metrics:");
        reportMetrics(changed, cfg, "Abrupt plant-change case metrics:");

        try {
            writeCsv(changed, Path.of("Chapter13_Lesson5_java_results.csv"));
        } catch (IOException error) {
            System.err.println("Could not write CSV: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script requires only base MATLAB. It also contains a Simulink realization map. A Simulink model can be assembled with four Integrator blocks for \( y,y_m,k_y,k_r \), a MATLAB Function block for the controller and adaptive law, a Saturation block for the actuator, and Scope or To Workspace blocks for logging.

Chapter13_Lesson5.m

% Chapter13_Lesson5.m
% Basic direct self-tuning regulator (direct STR) case study.
%
% The controller gains are adapted directly:
%       u = k_y*y + k_r*r
% No estimate of the unknown plant coefficients a(t), b(t) is formed.
%
% This script requires only base MATLAB. The final section explains the
% equivalent Simulink block realization.

clear; clc; close all;

cfg.dt = 1e-3;
cfg.tf = 30;
cfg.am = 2.5;
cfg.bm = 2.5;
cfg.gammaY = 12;
cfg.gammaR = 8;
cfg.uLimit = 6;
cfg.gainLimit = 12;
cfg.noiseStd = 0.01;
cfg.seed = 7;

nominal = simulateDirectSTR(cfg, false);
changed = simulateDirectSTR(cfg, true);

fprintf('Nominal case:\n');
printMetrics(nominal, cfg.dt);

fprintf('\nAbrupt parameter-change case:\n');
printMetrics(changed, cfg.dt);

figure('Name', 'Chapter 13 Lesson 5 - Direct STR');
tiledlayout(3,1);

nexttile;
plot(changed.t, changed.r, '--', 'LineWidth', 1.0); hold on;
plot(changed.t, changed.ym, 'LineWidth', 1.1);
plot(changed.t, changed.y, 'LineWidth', 1.1);
grid on;
ylabel('Output');
legend('r', 'y_m', 'y', 'Location', 'best');
title('Plant parameters change at t = 15 s');

nexttile;
plot(changed.t, changed.u, 'LineWidth', 1.0); hold on;
plot(changed.t, changed.e, 'LineWidth', 1.0);
grid on;
ylabel('u, e');
legend('u', 'e', 'Location', 'best');

nexttile;
plot(changed.t, changed.ky, 'LineWidth', 1.0); hold on;
plot(changed.t, changed.kr, 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Direct gains');
legend('k_y', 'k_r', 'Location', 'best');

% Simulink realization:
% 1. Use two Integrator blocks for y and y_m.
% 2. Plant derivative: -a(t)*y + b(t)*u.
% 3. Reference-model derivative: -am*y_m + bm*r.
% 4. MATLAB Function block inputs: y, y_m, r, k_y, k_r.
% 5. Two additional Integrator blocks integrate k_y_dot and k_r_dot.
% 6. Apply Saturation blocks to u, k_y, and k_r.
% 7. The MATLAB Function block equations are:
%
% function [u, kyDot, krDot] = directSTR(y, ym, r, ky, kr)
% %#codegen
% gammaY = 12; gammaR = 8; uLimit = 6;
% e = y - ym;
% u = min(max(ky*y + kr*r, -uLimit), uLimit);
% kyDot = -gammaY*e*y;
% krDot = -gammaR*e*r;
% end

function result = simulateDirectSTR(cfg, parameterJump)
    rng(cfg.seed + double(parameterJump));

    t = (0:cfg.dt:cfg.tf)';
    n = numel(t);
    y = zeros(n,1);
    ym = zeros(n,1);
    e = zeros(n,1);
    u = zeros(n,1);
    ky = zeros(n,1);
    kr = zeros(n,1);
    r = zeros(n,1);
    aHist = zeros(n,1);
    bHist = zeros(n,1);

    ky(1) = 0;
    kr(1) = 0.5;

    for k = 1:n-1
        r(k) = command(t(k));
        [a, b] = plantParameters(t(k), parameterJump);
        aHist(k) = a;
        bHist(k) = b;

        measuredY = y(k) + cfg.noiseStd*randn;
        e(k) = measuredY - ym(k);
        rawU = ky(k)*measuredY + kr(k)*r(k);
        u(k) = min(max(rawU, -cfg.uLimit), cfg.uLimit);

        kyDot = -cfg.gammaY*e(k)*measuredY;
        krDot = -cfg.gammaR*e(k)*r(k);

        ky(k+1) = min(max(ky(k) + cfg.dt*kyDot, ...
                          -cfg.gainLimit), cfg.gainLimit);
        kr(k+1) = min(max(kr(k) + cfg.dt*krDot, ...
                          -cfg.gainLimit), cfg.gainLimit);

        yDot = -a*y(k) + b*u(k);
        ymDot = -cfg.am*ym(k) + cfg.bm*r(k);
        y(k+1) = y(k) + cfg.dt*yDot;
        ym(k+1) = ym(k) + cfg.dt*ymDot;
    end

    r(end) = command(t(end));
    [aHist(end), bHist(end)] = plantParameters(t(end), parameterJump);
    e(end) = y(end) - ym(end);
    u(end) = min(max(ky(end)*y(end) + kr(end)*r(end), ...
                     -cfg.uLimit), cfg.uLimit);

    result = table(t, r, y, ym, e, u, ky, kr, aHist, bHist, ...
        'VariableNames', {'t','r','y','ym','e','u','ky','kr','a','b'});
end

function r = command(t)
    if t < 5
        r = 1;
    elseif t < 10
        r = -0.5;
    elseif t < 18
        r = 0.8;
    elseif t < 24
        r = 0.2;
    else
        r = -0.8;
    end
end

function [a, b] = plantParameters(t, parameterJump)
    if parameterJump && t >= 15
        a = 2.0;
        b = 0.65;
    else
        a = 1.2;
        b = 1.0;
    end
end

function printMetrics(data, dt)
    rmse = sqrt(mean(data.e.^2));
    iae = sum(abs(data.e))*dt;
    maxError = max(abs(data.e));
    controlRms = sqrt(mean(data.u.^2));
    totalVariation = sum(abs(diff(data.u)));

    fprintf('  RMSE: %.6f\n', rmse);
    fprintf('  IAE: %.6f\n', iae);
    fprintf('  Maximum absolute error: %.6f\n', maxError);
    fprintf('  Control RMS: %.6f\n', controlRms);
    fprintf('  Control total variation: %.6f\n', totalVariation);
end

Mathematica's NDSolveValue integrates the plant, reference-model, and gain dynamics as a coupled nonlinear differential system. The notebook computes RMSE and IAE and produces tracking, control-error, and gain plots.

Chapter13_Lesson5.nb

(* Chapter13_Lesson5.nb
   Basic direct self-tuning regulator case study.
   The gains ky(t) and kr(t) are adapted directly. *)

ClearAll["Global`*"];

tf = 30.;
am = 2.5;
bm = 2.5;
gammaY = 12.;
gammaR = 8.;
uLimit = 6.;

r[t_] := Piecewise[{
    {1., t < 5.},
    {-0.5, t < 10.},
    {0.8, t < 18.},
    {0.2, t < 24.}
    }, -0.8];

a[t_] := Piecewise[{ {1.2, t < 15.} }, 2.0];
b[t_] := Piecewise[{ {1.0, t < 15.} }, 0.65];

solution = NDSolveValue[
   {
    y'[t] == -a[t] y[t] +
      b[t] Clip[ky[t] y[t] + kr[t] r[t], {-uLimit, uLimit}],
    ym'[t] == -am ym[t] + bm r[t],
    ky'[t] == -gammaY (y[t] - ym[t]) y[t],
    kr'[t] == -gammaR (y[t] - ym[t]) r[t],
    y[0] == 0,
    ym[0] == 0,
    ky[0] == 0,
    kr[0] == 0.5
    },
   {y, ym, ky, kr},
   {t, 0, tf},
   MaxStepFraction -> 1/5000
   ];

{ySol, ymSol, kySol, krSol} = solution;

uSol[t_] := Clip[
   kySol[t] ySol[t] + krSol[t] r[t],
   {-uLimit, uLimit}
   ];

error[t_] := ySol[t] - ymSol[t];

rmse = Sqrt[NIntegrate[error[t]^2, {t, 0, tf}]/tf];
iae = NIntegrate[Abs[error[t]], {t, 0, tf}];

Print["RMSE = ", N[rmse, 6]];
Print["IAE = ", N[iae, 6]];

trackingPlot = Plot[
   {r[t], ymSol[t], ySol[t]},
   {t, 0, tf},
   PlotLegends -> {"r", "ym", "y"},
   PlotRange -> All,
   AxesLabel -> {"t", "output"},
   GridLines -> Automatic
   ];

gainPlot = Plot[
   {kySol[t], krSol[t]},
   {t, 0, tf},
   PlotLegends -> {"ky", "kr"},
   PlotRange -> All,
   AxesLabel -> {"t", "gain"},
   GridLines -> Automatic
   ];

controlPlot = Plot[
   {uSol[t], error[t]},
   {t, 0, tf},
   PlotLegends -> {"u", "e"},
   PlotRange -> All,
   AxesLabel -> {"t", "signal"},
   GridLines -> Automatic
   ];

Column[{trackingPlot, controlPlot, gainPlot}]

16.1 Initial Adaptation Transient

The gains begin at values that do not satisfy the matching equations. Therefore, the plant initially differs from the reference model and the nonzero tracking error drives gain adaptation. Larger adaptation gains generally shorten this transient, but they also increase sensitivity to noise and numerical integration errors.

16.2 Response to Command Changes

Each command step temporarily changes the tracking error and therefore injects information into the gain update. A constant command held for a long interval is usually insufficient to identify both ideal gains uniquely. Nevertheless, the controller can still produce small tracking error.

16.3 Response to the Abrupt Plant Change

When \( a \) and \( b \) change at 15 seconds, the previously learned gains no longer match the reference model. The tracking error reappears and redirects the gains toward values appropriate for the new plant. The output should recover, provided the plant retains the assumed control direction and the actuator has enough authority.

16.4 Why the Gains Need Not Equal the Ideal Values Exactly

The Lyapunov proof forces the tracking error toward zero but does not impose unique gain identification. If the regressor is not persistently exciting, several gain combinations may generate nearly identical trajectories over the observed command sequence. Small measurement noise, projection, saturation, and finite integration steps further prevent exact equality with the theoretical ideal gains.

17.1 Incorrect Control-Direction Assumption

The adaptive-law sign in Section 5 relies on \( b>0 \). If the true input gain is negative, the cancellation used in the Lyapunov proof has the wrong sign and the adaptive loop can destabilize. Unknown control direction is treated separately in Chapter 16.

17.2 Excessive Integration Step

The continuous-time proof does not automatically guarantee stability of a coarse forward-Euler implementation. The gain increments \( h\gamma_y e_k y_k \) and \( h\gamma_r e_k r_k \) should remain small relative to the gain scale. Reducing \( h \) is the first numerical diagnostic when oscillations appear.

17.3 Persistent Saturation

If the actuator remains saturated, the controller parameterization cannot realize the reference model. Continuing to adapt from the resulting error may push gains to their limits. A practical monitor should record the fraction of samples at saturation and pause or modify adaptation when that fraction becomes large.

17.4 Noise-Driven Gain Motion

A small output error does not imply that the gains are stationary. Plotting \( \dot k_y \), \( \dot k_r \), and \( \operatorname{TV}(u) \) reveals noise-driven motion that may be hidden in the output plot.

Problem 1 — Matching Conditions: For \( \dot y=-a y+b u \), \( u=k_y y+k_r r \), and \( \dot y_m=-a_m y_m+b_m r \), derive the ideal direct controller gains.

Solution:

With constant ideal gains, the plant becomes \( \dot y=(-a+b k_y^\ast)y+b k_r^\ast r \). Matching the coefficient of \( y \) to \( -a_m \) and the coefficient of \( r \) to \( b_m \) gives

\[ -a+b k_y^\ast=-a_m, \qquad b k_r^\ast=b_m. \]

Therefore,

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

Problem 2 — Lyapunov Cancellation: Starting from \( \dot e=-a_m e+b\tilde k_y y+b\tilde k_r r \), find adaptive laws that make \( \dot V=-a_m e^2 \) for

\[ V= \frac{1}{2}e^2+ \frac{b}{2\gamma_y}\tilde k_y^2+ \frac{b}{2\gamma_r}\tilde k_r^2. \]

Solution:

Differentiation gives

\[ \dot V = -a_m e^2 + b e\tilde k_y y + b e\tilde k_r r + \frac{b}{\gamma_y}\tilde k_y\dot k_y + \frac{b}{\gamma_r}\tilde k_r\dot k_r. \]

Select each gain derivative to cancel the associated cross term:

\[ \dot k_y=-\gamma_y e y, \qquad \dot k_r=-\gamma_r e r. \]

Substitution cancels both products and leaves \( \dot V=-a_m e^2\leq0 \).

Problem 3 — Effect of a Plant-Gain Change: Let \( a_m=b_m=2.5 \). Compute the ideal gains before and after the change \( (a,b):(1.2,1.0)→(2.0,0.65) \).

Solution:

Before the change,

\[ k_y^\ast=\frac{1.2-2.5}{1.0}=-1.3, \qquad k_r^\ast=\frac{2.5}{1.0}=2.5. \]

After the change,

\[ k_y^\ast=\frac{2.0-2.5}{0.65}\approx-0.7692, \qquad k_r^\ast=\frac{2.5}{0.65}\approx3.8462. \]

The command gain must increase because the plant input gain decreases. The feedback gain becomes less negative because the new open-loop decay coefficient is closer to the desired model decay coefficient.

Problem 4 — Sampled Gain Update: Suppose \( h=0.001 \), \( \gamma_y=12 \), \( \gamma_r=8 \), \( e_k=0.2 \), \( y_k=0.6 \), and \( r_k=1 \). Compute one gain-update step.

Solution:

\[ \Delta k_y = -h\gamma_y e_k y_k = -(0.001)(12)(0.2)(0.6) = -0.00144, \]

\[ \Delta k_r = -h\gamma_r e_k r_k = -(0.001)(8)(0.2)(1) = -0.0016. \]

Thus \( k_{y,k+1}=k_{y,k}-0.00144 \) and \( k_{r,k+1}=k_{r,k}-0.0016 \), before projection.

Problem 5 — Tracking Versus Gain Convergence: Explain why \( e(t)→0 \) does not by itself prove \( \tilde k_y(t)→0 \) and \( \tilde k_r(t)→0 \).

Solution:

The error equation depends on the scalar combination \( \tilde{\boldsymbol{k} }^T\boldsymbol{\phi} \). If the regressor explores only a restricted set of directions, a nonzero gain-error vector can remain orthogonal to the observed regressor and produce zero tracking error. Unique gain convergence requires enough independent excitation over time. Persistent excitation supplies this condition; boundedness and asymptotic tracking alone do not.

Problem 6 — Saturation and the Proof: Let \( u_{\mathrm{sat} }=u+\Delta_u \), where \( \Delta_u \) is the saturation discrepancy. Derive the additional term in the tracking-error and Lyapunov derivatives.

Solution:

The plant receives \( u+\Delta_u \), so the error equation becomes

\[ \dot e = -a_m e + b\tilde k_y y + b\tilde k_r r + b\Delta_u. \]

With the same adaptive laws, the gain-error cross terms still cancel, but

\[ \dot V = -a_m e^2+b e\Delta_u. \]

The final term has no fixed sign. This demonstrates precisely why actuator saturation invalidates the simple global proof whenever \( \Delta_u\neq0 \).

This lesson implemented a basic direct self-tuning regulator for an uncertain first-order plant. The controller gains were parameterized directly, ideal matching gains were derived only for analysis, and a Lyapunov argument established boundedness and asymptotic tracking for the ideal constant-parameter, noise-free, unsaturated system. The implementation then added numerical integration, gain projection, actuator limits, measurement noise, performance metrics, and an abrupt plant change. The case study also emphasized a central adaptive-control distinction: output tracking can converge even when the adaptive gains do not converge to unique ideal values.

1. Åström, K.J., & Wittenmark, B. (1973). On self-tuning regulators. Automatica, 9(2), 185–199.
2. Clarke, D.W., & Gawthrop, P.J. (1975). Self-tuning controller. Proceedings of the Institution of Electrical Engineers, 122(9), 929–934.
3. Landau, I.D. (1974). A survey of model reference adaptive techniques—Theory and applications. Automatica, 10(4), 353–379.
4. Narendra, K.S., & Valavani, L.S. (1978). Stable adaptive controller design—Direct control. IEEE Transactions on Automatic Control, 23(4), 570–583.
5. Morse, A.S. (1980). Global stability of parameter-adaptive control systems. IEEE Transactions on Automatic Control, 25(3), 433–439.
6. Goodwin, G.C., Ramadge, P.J., & Caines, P.E. (1980). Discrete-time multivariable adaptive control. IEEE Transactions on Automatic Control, 25(3), 449–456.
7. Ioannou, P.A., & Tsakalis, K.S. (1986). A robust direct adaptive controller. IEEE Transactions on Automatic Control, 31(11), 1033–1043.
8. Rohrs, C.E., Valavani, L., Athans, M., & Stein, G. (1985). Robustness of continuous-time adaptive control algorithms in the presence of unmodeled dynamics. IEEE Transactions on Automatic Control, 30(9), 881–889.