After completing this lesson, students should be able to:

• form a control-oriented linear regression for a simple dynamic process;
• derive RLS with exponential forgetting from a weighted least-squares cost;
• derive a certainty-equivalent pole-placement control law from current parameter estimates;
• obtain the closed-loop tracking-error equation in terms of parameter error;
• explain why tracking can be satisfactory even when parameters have not fully converged;
• implement projection, denominator protection, covariance symmetrization, probing, and input saturation;
• interpret online estimates, prediction errors, tracking errors, and control effort together.

The lesson deliberately uses a first-order process so that every estimation and control equation can be derived explicitly. Higher-order polynomial STR design is not introduced here; it belongs after the discrete-time adaptive control foundations developed in later chapters.

An indirect self-tuning regulator separates adaptation into two connected computations. The identifier estimates a plant model, while the controller design map converts that model into feedback and feedforward gains. This is the certainty-equivalence principle: at each sample, the estimated model is treated as though it were the true model for the purpose of control design.

flowchart TD
  R["Reference r(k)"] --> C["Controller redesign from aHat(k), bHat(k)"]
  C --> U["Bounded input u(k)"]
  U --> P["Unknown process"]
  P --> Y["Measured output y(k+1)"]
  Y --> I["RLS identifier"]
  U --> I
  I --> E["Updated estimates aHat(k+1), bHat(k+1)"]
  E --> C
  Y --> C
        

The loop is nonlinear and time varying even though the process model is linear, because the controller coefficients depend on recursively updated estimates. Consequently, fixed-parameter pole-placement intuition alone is not a complete stability proof for the adaptive closed loop.

Consider the unknown first-order discrete-time process

\[ y(k+1)=a\,y(k)+b\,u(k)+v(k+1), \qquad b\neq 0, \]

where \( y(k) \) is the measured output, \( u(k) \) is the control input, and \( v(k+1) \) collects measurement noise, disturbances, and small unmodeled effects. Define

\[ \boldsymbol{\theta}=\begin{bmatrix}a\\b\end{bmatrix}, \qquad \boldsymbol{\phi}(k)=\begin{bmatrix}y(k)\\u(k)\end{bmatrix}. \]

The one-step model is linear in the unknown parameters:

\[ y(k+1)=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}+v(k+1). \]

At sample \( k \), the available estimate is \( \hat{\boldsymbol{\theta}}(k) \). The one-step prediction and innovation are

\[ \hat{y}(k+1\mid k)=\boldsymbol{\phi}^{T}(k) \hat{\boldsymbol{\theta}}(k), \qquad \varepsilon(k+1)=y(k+1)-\hat{y}(k+1\mid k). \]

The sign of \( b \) is assumed known in this introductory STR example. Its magnitude is unknown. Unknown control direction requires a different design and is treated in Chapter 16.

Let \( 0<\lambda\leq 1 \) be the forgetting factor. After observing data through sample \( k \), define the regularized weighted cost

\[ \begin{aligned} J_{k+1}(\boldsymbol{\vartheta}) &=\lambda^{k+1} (\boldsymbol{\vartheta}-\hat{\boldsymbol{\theta}}(0))^{T} \mathbf{P}^{-1}(0) (\boldsymbol{\vartheta}-\hat{\boldsymbol{\theta}}(0)) \\ &\quad+\sum_{i=0}^{k}\lambda^{k-i} \left[y(i+1)-\boldsymbol{\phi}^{T}(i)\boldsymbol{\vartheta}\right]^{2}. \end{aligned} \]

Differentiating with respect to \( \boldsymbol{\vartheta} \) and setting the gradient to zero gives the normal equation

\[ \mathbf{R}(k+1)\hat{\boldsymbol{\theta}}(k+1)=\mathbf{q}(k+1), \]

with recursions

\[ \begin{aligned} \mathbf{R}(k+1)&=\lambda\mathbf{R}(k)+ \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k),\\ \mathbf{q}(k+1)&=\lambda\mathbf{q}(k)+ \boldsymbol{\phi}(k)y(k+1). \end{aligned} \]

Directly inverting \( \mathbf{R}(k+1) \) at every sample is unnecessary. Define \( \mathbf{P}(k)=\mathbf{R}^{-1}(k) \). Applying the matrix inversion lemma to \( \lambda\mathbf{R}(k)+\boldsymbol{\phi}(k) \boldsymbol{\phi}^{T}(k) \) yields the covariance recursion used by RLS.

\[ \left(\mathbf{A}+\mathbf{u}\mathbf{v}^{T}\right)^{-1} =\mathbf{A}^{-1}- \frac{\mathbf{A}^{-1}\mathbf{u}\mathbf{v}^{T}\mathbf{A}^{-1}} {1+\mathbf{v}^{T}\mathbf{A}^{-1}\mathbf{u}}. \]

Use \( \mathbf{A}=\lambda\mathbf{R}(k) \) and \( \mathbf{u}=\mathbf{v}=\boldsymbol{\phi}(k) \). Because \( \mathbf{A}^{-1}=\lambda^{-1}\mathbf{P}(k) \), the result is

\[ \mathbf{P}(k+1)=\frac{1}{\lambda}\left[ \mathbf{P}(k)- \frac{\mathbf{P}(k)\boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k) \mathbf{P}(k)} {\lambda+\boldsymbol{\phi}^{T}(k)\mathbf{P}(k) \boldsymbol{\phi}(k)}\right]. \]

If \( \mathbf{P}(0) \) is symmetric positive definite and \( \lambda>0 \), the exact-arithmetic information matrix remains positive definite. Finite-precision arithmetic can slowly destroy symmetry, so practical implementations replace \( \mathbf{P} \) by \( (\mathbf{P}+\mathbf{P}^{T})/2 \) after every update.

The implementable RLS recursion is

\[ \mathbf{K}(k)= \frac{\mathbf{P}(k)\boldsymbol{\phi}(k)} {\lambda+\boldsymbol{\phi}^{T}(k)\mathbf{P}(k) \boldsymbol{\phi}(k)}, \]

\[ \hat{\boldsymbol{\theta}}(k+1)= \hat{\boldsymbol{\theta}}(k)+\mathbf{K}(k)\varepsilon(k+1), \]

\[ \mathbf{P}(k+1)=\frac{1}{\lambda} \left[\mathbf{P}(k)-\mathbf{K}(k) \boldsymbol{\phi}^{T}(k)\mathbf{P}(k)\right]. \]

The gain \( \mathbf{K}(k) \) is large in parameter-space directions that remain uncertain and that are excited by the current regressor. The denominator is strictly positive in exact arithmetic when \( \mathbf{P}(k) \) is positive semidefinite and \( \lambda>0 \).

For \( \lambda=1 \), all historical samples retain equal weight. For \( \lambda<1 \), a sample that is \( j \) steps old has relative weight \( \lambda^{j} \). Useful memory approximations are

\[ N_{\mathrm{eff}}\approx\frac{1}{1-\lambda}, \qquad N_{1/2}=\frac{\ln(1/2)}{\ln(\lambda)}. \]

Smaller forgetting factors adapt more rapidly to parameter changes but increase estimate variance and sensitivity to noise. This is the same speed-versus-noise trade-off emphasized when RLS was compared with gradient adaptation in Lesson 4.

Select a desired scalar closed-loop pole \( p_m \) satisfying \( |p_m|<1 \). The desired first-order command response is

\[ y_d(k+1)=p_m y_d(k)+(1-p_m)r(k). \]

If the true parameters were known, equating the plant update to the desired update would give

\[ a y(k)+b u(k)=p_m y(k)+(1-p_m)r(k). \]

Solving for the ideal input:

\[ u^{\star}(k)=\frac{p_m-a}{b}y(k)+ \frac{1-p_m}{b}r(k). \]

The indirect STR replaces the unknown parameters with their current RLS estimates:

\[ u_{\mathrm{ce}}(k)= \frac{p_m-\hat{a}(k)}{\hat{b}(k)}y(k)+ \frac{1-p_m}{\hat{b}(k)}r(k). \]

The denominator makes the controller highly sensitive to an estimate near zero. With known positive control direction, the practical implementation uses a projected estimate \( \hat{b}(k)\in[b_{\min},b_{\max}] \), where \( b_{\min}>0 \).

6.1 Closed-Loop Error Equation

First ignore probing and saturation. Define the parameter error \( \tilde{\boldsymbol{\theta}}(k)= \boldsymbol{\theta}-\hat{\boldsymbol{\theta}}(k) \). Because the controller enforces

\[ \hat{a}(k)y(k)+\hat{b}(k)u_{\mathrm{ce}}(k) =p_m y(k)+(1-p_m)r(k), \]

the actual process satisfies

\[ y(k+1)=p_m y(k)+(1-p_m)r(k)+ \tilde{\boldsymbol{\theta}}^{T}(k) \boldsymbol{\phi}(k)+v(k+1). \]

Let the reference-model state obey the same desired recursion and define \( e_c(k)=y(k)-y_d(k) \). Then

\[ e_c(k+1)=p_m e_c(k)+ \tilde{\boldsymbol{\theta}}^{T}(k) \boldsymbol{\phi}(k)+v(k+1). \]

This equation shows the central STR mechanism. The nominal tracking dynamics are stable because \( |p_m|<1 \); parameter error enters as an additive forcing term. Exact parameter convergence drives that term toward zero, but good tracking may occur earlier if the particular combination \( \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} \) is small.

In the implementation, the applied input is \( u(k)=\operatorname{sat}(u_{\mathrm{ce}}(k)+d_p(k)) \), where \( d_p(k) \) is a small probing signal. These modifications add a bounded input-mismatch term to the ideal error equation. They improve identifiability and safety but prevent exact equality with the nominal pole-placement equation during probing or saturation.

flowchart TD
  A["Start sample k with y(k), thetaHat(k), P(k)"] --> B["Read reference r(k)"]
  B --> C["Protect bHat and compute certainty-equivalent input"]
  C --> D["Add small probe and apply input saturation"]
  D --> E["Apply u(k) to process and measure y(k+1)"]
  E --> F["Build regressor phi(k) = [y(k), u(k)]"]
  F --> G["Compute prediction error and RLS gain"]
  G --> H["Update thetaHat(k+1) and P(k+1)"]
  H --> I["Project parameters and symmetrize covariance"]
  I --> J["Store diagnostics and advance to k+1"]
  J --> B
        

7.1 Initialization

\[ \hat{\boldsymbol{\theta}}(0)= \begin{bmatrix}\hat{a}(0)\\\hat{b}(0)\end{bmatrix}, \qquad \mathbf{P}(0)=\alpha\mathbf{I}, \quad \alpha>0. \]

A large \( \alpha \) represents substantial initial uncertainty and permits rapid initial adaptation. It can also create large early gains, so controller saturation and parameter projection are essential.

7.2 Practical Safeguards

• Known-sign projection: constrain \( \hat{b} \) away from zero.
• Input saturation: enforce actuator limits before applying the command.
• Covariance symmetrization: remove roundoff-induced asymmetry.
• Small probing signal: preserve information when regulation would otherwise remove excitation.
• Signal monitoring: log prediction error, tracking error, estimates, covariance, and saturation events.
• Finite-value checks: stop or reset safely if the RLS denominator or estimates become non-finite.

These safeguards do not constitute a general stability theorem. They are implementation measures consistent with the projection, normalization, excitation, and robustness concepts already developed in Chapters 8–10.

All implementations use the same deterministic benchmark:

The reproducible disturbance is

\[ v(k+1)=0.01\sin(0.17k)+0.005\cos(0.043k). \]

In the verified Python, C++, and Java runs, the final estimates were approximately \( \hat{a}=0.820780 \) and \( \hat{b}=0.179812 \). The tracking RMSE after sample 50 was approximately \( 0.135492 \), and the one-step prediction RMSE was approximately \( 0.008011 \). Exact floating-point output may vary slightly by platform.

This version uses NumPy for vector and matrix operations and Matplotlib for diagnostics. The RLS and controller equations are implemented from scratch rather than hidden inside an identification package.

Chapter12_Lesson5.py

"""Chapter 12, Lesson 5: RLS-based self-tuning regulator example.

Requires:
    numpy
    matplotlib

The benchmark is deliberately deterministic so its results can be compared with
C++, Java, MATLAB, Simulink, and Wolfram Mathematica implementations.
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    samples: int = 600
    a_true: float = 0.82
    b_true: float = 0.18
    desired_pole: float = 0.35
    forgetting_factor: float = 0.995
    initial_covariance: float = 1000.0
    u_min: float = -4.0
    u_max: float = 4.0
    a_min: float = -0.98
    a_max: float = 0.98
    b_min: float = 0.05
    b_max: float = 0.80
    probe_amplitude: float = 0.04


def reference(k: int) -> float:
    """Piecewise-constant command used by all implementations."""
    if k < 40:
        return 0.0
    if k < 200:
        return 1.0
    if k < 340:
        return -0.5
    if k < 470:
        return 1.5
    return 0.25


def deterministic_probe(k: int, amplitude: float) -> float:
    """Small binary probing signal that prevents immediate loss of excitation."""
    return amplitude if (37 * k) % 23 < 11 else -amplitude


def deterministic_disturbance(k: int) -> float:
    """Small reproducible output disturbance."""
    return 0.01 * math.sin(0.17 * k) + 0.005 * math.cos(0.043 * k)


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


def simulate(cfg: Config = Config()) -> dict[str, np.ndarray | float]:
    n = cfg.samples

    y = np.zeros(n + 1)
    u = np.zeros(n)
    r = np.zeros(n)
    prediction_error = np.zeros(n)
    a_hat_history = np.zeros(n + 1)
    b_hat_history = np.zeros(n + 1)

    theta_hat = np.array([0.20, 0.40], dtype=float)
    covariance = cfg.initial_covariance * np.eye(2)
    a_hat_history[0], b_hat_history[0] = theta_hat

    for k in range(n):
        r[k] = reference(k)

        # Certainty-equivalent pole-placement controller:
        # a_hat*y + b_hat*u = p_m*y + (1-p_m)*r.
        safe_b_hat = max(theta_hat[1], cfg.b_min)
        u_ce = (
            (cfg.desired_pole - theta_hat[0]) * y[k]
            + (1.0 - cfg.desired_pole) * r[k]
        ) / safe_b_hat

        u[k] = clip(
            u_ce + deterministic_probe(k, cfg.probe_amplitude),
            cfg.u_min,
            cfg.u_max,
        )

        # True process: y(k+1) = a*y(k) + b*u(k) + v(k+1).
        y[k + 1] = (
            cfg.a_true * y[k]
            + cfg.b_true * u[k]
            + deterministic_disturbance(k)
        )

        # RLS uses phi(k) = [y(k), u(k)]^T to predict y(k+1).
        phi = np.array([y[k], u[k]])
        denominator = cfg.forgetting_factor + phi @ covariance @ phi
        gain = covariance @ phi / denominator
        prediction_error[k] = y[k + 1] - phi @ theta_hat
        theta_hat = theta_hat + gain * prediction_error[k]
        covariance = (
            covariance - np.outer(gain, phi @ covariance)
        ) / cfg.forgetting_factor

        # Numerical safeguards already introduced in earlier lessons.
        covariance = 0.5 * (covariance + covariance.T)
        theta_hat[0] = clip(theta_hat[0], cfg.a_min, cfg.a_max)
        theta_hat[1] = clip(theta_hat[1], cfg.b_min, cfg.b_max)

        a_hat_history[k + 1], b_hat_history[k + 1] = theta_hat

    tracking_rmse = float(np.sqrt(np.mean((y[50:n] - r[50:n]) ** 2)))
    prediction_rmse = float(np.sqrt(np.mean(prediction_error[50:] ** 2)))

    return {
        "y": y,
        "u": u,
        "r": r,
        "prediction_error": prediction_error,
        "a_hat": a_hat_history,
        "b_hat": b_hat_history,
        "tracking_rmse": tracking_rmse,
        "prediction_rmse": prediction_rmse,
    }


def main() -> None:
    cfg = Config()
    data = simulate(cfg)

    print(f"Final a_hat: {data['a_hat'][-1]:.6f} (true {cfg.a_true:.6f})")
    print(f"Final b_hat: {data['b_hat'][-1]:.6f} (true {cfg.b_true:.6f})")
    print(f"Tracking RMSE after sample 50: {data['tracking_rmse']:.6f}")
    print(f"Prediction RMSE after sample 50: {data['prediction_rmse']:.6f}")

    k = np.arange(cfg.samples)

    plt.figure()
    plt.step(k, data["r"], where="post", label="reference r(k)")
    plt.plot(np.arange(cfg.samples + 1), data["y"], label="output y(k)")
    plt.xlabel("sample k")
    plt.ylabel("output")
    plt.title("RLS-based STR tracking")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(data["a_hat"], label="a_hat")
    plt.plot(data["b_hat"], label="b_hat")
    plt.axhline(cfg.a_true, linestyle="--", label="true a")
    plt.axhline(cfg.b_true, linestyle="--", label="true b")
    plt.xlabel("sample k")
    plt.ylabel("parameter value")
    plt.title("Online parameter estimates")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(k, data["u"])
    plt.xlabel("sample k")
    plt.ylabel("u(k)")
    plt.title("Control input")
    plt.grid(True)

    plt.show()


if __name__ == "__main__":
    main()

This C++17 implementation uses only the standard library, explicit two-dimensional matrix operations, finite-value checks, and CSV export. It is suitable as a starting point for embedded or real-time restructuring.

Chapter12_Lesson5.cpp

// Chapter 12, Lesson 5: RLS-based self-tuning regulator example.
// Build: g++ -std=c++17 -O2 Chapter12_Lesson5.cpp -o Chapter12_Lesson5

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

namespace {

struct Config {
    int samples = 600;
    double aTrue = 0.82;
    double bTrue = 0.18;
    double desiredPole = 0.35;
    double forgettingFactor = 0.995;
    double initialCovariance = 1000.0;
    double uMin = -4.0;
    double uMax = 4.0;
    double aMin = -0.98;
    double aMax = 0.98;
    double bMin = 0.05;
    double bMax = 0.80;
    double probeAmplitude = 0.04;
};

double reference(const int k) {
    if (k < 40) return 0.0;
    if (k < 200) return 1.0;
    if (k < 340) return -0.5;
    if (k < 470) return 1.5;
    return 0.25;
}

double deterministicProbe(const int k, const double amplitude) {
    return ((37 * k) % 23 < 11) ? amplitude : -amplitude;
}

double deterministicDisturbance(const int k) {
    return 0.01 * std::sin(0.17 * k) + 0.005 * std::cos(0.043 * k);
}

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

}  // namespace

int main() {
    try {
        const Config cfg;
        const int n = cfg.samples;

        std::vector<double> y(n + 1, 0.0);
        std::vector<double> u(n, 0.0);
        std::vector<double> r(n, 0.0);
        std::vector<double> predictionError(n, 0.0);
        std::vector<double> aHatHistory(n + 1, 0.0);
        std::vector<double> bHatHistory(n + 1, 0.0);

        std::array<double, 2> thetaHat{0.20, 0.40};
        std::array<std::array<double, 2>, 2> P{{
            {cfg.initialCovariance, 0.0},
            {0.0, cfg.initialCovariance}
        }};

        aHatHistory[0] = thetaHat[0];
        bHatHistory[0] = thetaHat[1];

        for (int k = 0; k < n; ++k) {
            r[k] = reference(k);

            const double safeBHat = std::max(thetaHat[1], cfg.bMin);
            const double uCE = ((cfg.desiredPole - thetaHat[0]) * y[k]
                              + (1.0 - cfg.desiredPole) * r[k]) / safeBHat;
            u[k] = clip(uCE + deterministicProbe(k, cfg.probeAmplitude),
                        cfg.uMin, cfg.uMax);

            y[k + 1] = cfg.aTrue * y[k] + cfg.bTrue * u[k]
                     + deterministicDisturbance(k);

            const std::array<double, 2> phi{y[k], u[k]};
            const std::array<double, 2> Pphi{
                P[0][0] * phi[0] + P[0][1] * phi[1],
                P[1][0] * phi[0] + P[1][1] * phi[1]
            };
            const double phiTPphi = phi[0] * Pphi[0] + phi[1] * Pphi[1];
            const double denominator = cfg.forgettingFactor + phiTPphi;
            if (!(denominator > 0.0) || !std::isfinite(denominator)) {
                throw std::runtime_error("RLS denominator became invalid.");
            }

            const std::array<double, 2> gain{
                Pphi[0] / denominator,
                Pphi[1] / denominator
            };
            const double prediction = phi[0] * thetaHat[0] + phi[1] * thetaHat[1];
            predictionError[k] = y[k + 1] - prediction;

            thetaHat[0] += gain[0] * predictionError[k];
            thetaHat[1] += gain[1] * predictionError[k];

            const std::array<double, 2> phiTP{
                phi[0] * P[0][0] + phi[1] * P[1][0],
                phi[0] * P[0][1] + phi[1] * P[1][1]
            };

            std::array<std::array<double, 2>, 2> Pnew{};
            for (int i = 0; i < 2; ++i) {
                for (int j = 0; j < 2; ++j) {
                    Pnew[i][j] = (P[i][j] - gain[i] * phiTP[j])
                               / cfg.forgettingFactor;
                }
            }

            // Enforce symmetry to suppress roundoff-induced covariance drift.
            const double offDiagonal = 0.5 * (Pnew[0][1] + Pnew[1][0]);
            Pnew[0][1] = offDiagonal;
            Pnew[1][0] = offDiagonal;
            P = Pnew;

            thetaHat[0] = clip(thetaHat[0], cfg.aMin, cfg.aMax);
            thetaHat[1] = clip(thetaHat[1], cfg.bMin, cfg.bMax);

            aHatHistory[k + 1] = thetaHat[0];
            bHatHistory[k + 1] = thetaHat[1];
        }

        double trackingSse = 0.0;
        double predictionSse = 0.0;
        int trackingCount = 0;
        for (int k = 50; k < n; ++k) {
            const double trackingError = y[k] - r[k];
            trackingSse += trackingError * trackingError;
            predictionSse += predictionError[k] * predictionError[k];
            ++trackingCount;
        }

        const double trackingRmse = std::sqrt(trackingSse / trackingCount);
        const double predictionRmse = std::sqrt(predictionSse / trackingCount);

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Final a_hat: " << thetaHat[0]
                  << " (true " << cfg.aTrue << ")\n";
        std::cout << "Final b_hat: " << thetaHat[1]
                  << " (true " << cfg.bTrue << ")\n";
        std::cout << "Tracking RMSE after sample 50: " << trackingRmse << '\n';
        std::cout << "Prediction RMSE after sample 50: " << predictionRmse << '\n';

        std::ofstream csv("Chapter12_Lesson5_cpp_results.csv");
        if (!csv) {
            throw std::runtime_error("Could not create CSV output file.");
        }
        csv << "k,r,y,u,a_hat,b_hat,prediction_error\n";
        for (int k = 0; k < n; ++k) {
            csv << k << ',' << r[k] << ',' << y[k] << ',' << u[k] << ','
                << aHatHistory[k] << ',' << bHatHistory[k] << ','
                << predictionError[k] << '\n';
        }
        csv << n << ",," << y[n] << ",," << aHatHistory[n] << ','
            << bHatHistory[n] << ",\n";

        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

This implementation uses only the Java standard library. Arrays represent the two-parameter vector and covariance matrix, while a CSV file records the complete adaptation history.

Chapter12_Lesson5.java

// Chapter 12, Lesson 5: RLS-based self-tuning regulator example.
// Build and run:
//   javac Chapter12_Lesson5.java
//   java Chapter12_Lesson5

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

public final class Chapter12_Lesson5 {
    private static final int SAMPLES = 600;
    private static final double A_TRUE = 0.82;
    private static final double B_TRUE = 0.18;
    private static final double DESIRED_POLE = 0.35;
    private static final double FORGETTING_FACTOR = 0.995;
    private static final double INITIAL_COVARIANCE = 1000.0;
    private static final double U_MIN = -4.0;
    private static final double U_MAX = 4.0;
    private static final double A_MIN = -0.98;
    private static final double A_MAX = 0.98;
    private static final double B_MIN = 0.05;
    private static final double B_MAX = 0.80;
    private static final double PROBE_AMPLITUDE = 0.04;

    private Chapter12_Lesson5() {
        // Utility class.
    }

    private static double reference(int k) {
        if (k < 40) return 0.0;
        if (k < 200) return 1.0;
        if (k < 340) return -0.5;
        if (k < 470) return 1.5;
        return 0.25;
    }

    private static double deterministicProbe(int k) {
        return ((37 * k) % 23 < 11) ? PROBE_AMPLITUDE : -PROBE_AMPLITUDE;
    }

    private static double deterministicDisturbance(int k) {
        return 0.01 * Math.sin(0.17 * k) + 0.005 * Math.cos(0.043 * k);
    }

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

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

        double[] y = new double[SAMPLES + 1];
        double[] u = new double[SAMPLES];
        double[] r = new double[SAMPLES];
        double[] predictionError = new double[SAMPLES];
        double[] aHatHistory = new double[SAMPLES + 1];
        double[] bHatHistory = new double[SAMPLES + 1];

        double[] thetaHat = {0.20, 0.40};
        double[][] p = {
            {INITIAL_COVARIANCE, 0.0},
            {0.0, INITIAL_COVARIANCE}
        };

        aHatHistory[0] = thetaHat[0];
        bHatHistory[0] = thetaHat[1];

        for (int k = 0; k < SAMPLES; k++) {
            r[k] = reference(k);

            double safeBHat = Math.max(thetaHat[1], B_MIN);
            double uCE = ((DESIRED_POLE - thetaHat[0]) * y[k]
                    + (1.0 - DESIRED_POLE) * r[k]) / safeBHat;
            u[k] = clip(uCE + deterministicProbe(k), U_MIN, U_MAX);

            y[k + 1] = A_TRUE * y[k] + B_TRUE * u[k]
                    + deterministicDisturbance(k);

            double[] phi = {y[k], u[k]};
            double[] pPhi = {
                p[0][0] * phi[0] + p[0][1] * phi[1],
                p[1][0] * phi[0] + p[1][1] * phi[1]
            };
            double phiTPPhi = phi[0] * pPhi[0] + phi[1] * pPhi[1];
            double denominator = FORGETTING_FACTOR + phiTPPhi;
            if (!(denominator > 0.0) || !Double.isFinite(denominator)) {
                throw new IllegalStateException("RLS denominator became invalid.");
            }

            double[] gain = {
                pPhi[0] / denominator,
                pPhi[1] / denominator
            };
            double prediction = phi[0] * thetaHat[0] + phi[1] * thetaHat[1];
            predictionError[k] = y[k + 1] - prediction;

            thetaHat[0] += gain[0] * predictionError[k];
            thetaHat[1] += gain[1] * predictionError[k];

            double[] phiTP = {
                phi[0] * p[0][0] + phi[1] * p[1][0],
                phi[0] * p[0][1] + phi[1] * p[1][1]
            };

            double[][] pNew = new double[2][2];
            for (int i = 0; i < 2; i++) {
                for (int j = 0; j < 2; j++) {
                    pNew[i][j] = (p[i][j] - gain[i] * phiTP[j])
                            / FORGETTING_FACTOR;
                }
            }

            double offDiagonal = 0.5 * (pNew[0][1] + pNew[1][0]);
            pNew[0][1] = offDiagonal;
            pNew[1][0] = offDiagonal;
            p = pNew;

            thetaHat[0] = clip(thetaHat[0], A_MIN, A_MAX);
            thetaHat[1] = clip(thetaHat[1], B_MIN, B_MAX);

            aHatHistory[k + 1] = thetaHat[0];
            bHatHistory[k + 1] = thetaHat[1];
        }

        double trackingSse = 0.0;
        double predictionSse = 0.0;
        int count = 0;
        for (int k = 50; k < SAMPLES; k++) {
            double trackingError = y[k] - r[k];
            trackingSse += trackingError * trackingError;
            predictionSse += predictionError[k] * predictionError[k];
            count++;
        }

        double trackingRmse = Math.sqrt(trackingSse / count);
        double predictionRmse = Math.sqrt(predictionSse / count);

        System.out.printf("Final a_hat: %.6f (true %.6f)%n", thetaHat[0], A_TRUE);
        System.out.printf("Final b_hat: %.6f (true %.6f)%n", thetaHat[1], B_TRUE);
        System.out.printf("Tracking RMSE after sample 50: %.6f%n", trackingRmse);
        System.out.printf("Prediction RMSE after sample 50: %.6f%n", predictionRmse);

        writeCsv(y, u, r, predictionError, aHatHistory, bHatHistory);
    }

    private static void writeCsv(
            double[] y,
            double[] u,
            double[] r,
            double[] predictionError,
            double[] aHatHistory,
            double[] bHatHistory) {
        Path output = Path.of("Chapter12_Lesson5_java_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output)) {
            writer.write("k,r,y,u,a_hat,b_hat,prediction_error\n");
            for (int k = 0; k < SAMPLES; k++) {
                writer.write(String.format(
                        Locale.US,
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        k, r[k], y[k], u[k], aHatHistory[k],
                        bHatHistory[k], predictionError[k]));
            }
            writer.write(String.format(
                    Locale.US,
                    "%d,,%.12f,,%.12f,%.12f,%n",
                    SAMPLES, y[SAMPLES], aHatHistory[SAMPLES],
                    bHatHistory[SAMPLES]));
        } catch (IOException ex) {
            throw new RuntimeException("Could not write CSV output.", ex);
        }
    }
}

The MATLAB script uses core matrix operations and local functions. System Identification Toolbox is not required because the complete forgetting-factor RLS recursion is written explicitly.

Chapter12_Lesson5.m

%% Chapter 12, Lesson 5: RLS-Based STR Example for a Simple Process
% This script uses only core MATLAB functionality.
% It implements the same deterministic benchmark as the Python, C++, Java,
% Simulink, and Wolfram Mathematica versions.

clear; clc; close all;

%% Configuration
N = 600;
aTrue = 0.82;
bTrue = 0.18;
pm = 0.35;
lambda = 0.995;
P = 1000 * eye(2);
thetaHat = [0.20; 0.40];

uMin = -4.0;
uMax = 4.0;
aMin = -0.98;
aMax = 0.98;
bMin = 0.05;
bMax = 0.80;
probeAmplitude = 0.04;

%% Storage
y = zeros(N + 1, 1);
u = zeros(N, 1);
r = zeros(N, 1);
predictionError = zeros(N, 1);
aHatHistory = zeros(N + 1, 1);
bHatHistory = zeros(N + 1, 1);
aHatHistory(1) = thetaHat(1);
bHatHistory(1) = thetaHat(2);

%% Adaptive control loop
for k = 0:(N - 1)
    idx = k + 1;
    r(idx) = referenceSignal(k);

    % Certainty-equivalent pole placement:
    % aHat*y(k) + bHat*u(k) = pm*y(k) + (1-pm)*r(k).
    safeBHat = max(thetaHat(2), bMin);
    uCE = ((pm - thetaHat(1)) * y(idx) + (1 - pm) * r(idx)) / safeBHat;
    probe = deterministicProbe(k, probeAmplitude);
    u(idx) = clipScalar(uCE + probe, uMin, uMax);

    % True process.
    disturbance = 0.01 * sin(0.17 * k) + 0.005 * cos(0.043 * k);
    y(idx + 1) = aTrue * y(idx) + bTrue * u(idx) + disturbance;

    % RLS update using phi(k) = [y(k); u(k)] to predict y(k+1).
    phi = [y(idx); u(idx)];
    denominator = lambda + phi' * P * phi;
    if ~(isfinite(denominator) && denominator > 0)
        error('RLS denominator became invalid.');
    end

    K = (P * phi) / denominator;
    predictionError(idx) = y(idx + 1) - phi' * thetaHat;
    thetaHat = thetaHat + K * predictionError(idx);
    P = (P - K * (phi' * P)) / lambda;

    % Numerical safeguards.
    P = 0.5 * (P + P');
    thetaHat(1) = clipScalar(thetaHat(1), aMin, aMax);
    thetaHat(2) = clipScalar(thetaHat(2), bMin, bMax);

    aHatHistory(idx + 1) = thetaHat(1);
    bHatHistory(idx + 1) = thetaHat(2);
end

%% Metrics
metricIndices = 51:N;
trackingRmse = sqrt(mean((y(metricIndices) - r(metricIndices)).^2));
predictionRmse = sqrt(mean(predictionError(metricIndices).^2));

fprintf('Final a_hat: %.6f (true %.6f)\n', thetaHat(1), aTrue);
fprintf('Final b_hat: %.6f (true %.6f)\n', thetaHat(2), bTrue);
fprintf('Tracking RMSE after sample 50: %.6f\n', trackingRmse);
fprintf('Prediction RMSE after sample 50: %.6f\n', predictionRmse);

%% Plots
k = (0:(N - 1))';

figure;
stairs(k, r, 'LineWidth', 1.2); hold on;
plot((0:N)', y, 'LineWidth', 1.2);
grid on;
xlabel('sample k'); ylabel('output');
title('RLS-based STR tracking');
legend('reference r(k)', 'output y(k)', 'Location', 'best');

figure;
plot((0:N)', aHatHistory, 'LineWidth', 1.2); hold on;
plot((0:N)', bHatHistory, 'LineWidth', 1.2);
yline(aTrue, '--');
yline(bTrue, '--');
grid on;
xlabel('sample k'); ylabel('parameter value');
title('Online parameter estimates');
legend('a hat', 'b hat', 'true a', 'true b', 'Location', 'best');

figure;
plot(k, u, 'LineWidth', 1.2);
grid on;
xlabel('sample k'); ylabel('u(k)');
title('Control input');

%% Save results
results = table(k, r, y(1:N), u, aHatHistory(1:N), bHatHistory(1:N), ...
    predictionError, 'VariableNames', ...
    {'k', 'r', 'y', 'u', 'a_hat', 'b_hat', 'prediction_error'});
writetable(results, 'Chapter12_Lesson5_matlab_results.csv');

%% Local functions
function r = referenceSignal(k)
    if k < 40
        r = 0.0;
    elseif k < 200
        r = 1.0;
    elseif k < 340
        r = -0.5;
    elseif k < 470
        r = 1.5;
    else
        r = 0.25;
    end
end

function value = deterministicProbe(k, amplitude)
    if mod(37 * k, 23) < 11
        value = amplitude;
    else
        value = -amplitude;
    end
end

function value = clipScalar(value, lowerBound, upperBound)
    value = min(max(value, lowerBound), upperBound);
end

This script requires Simulink and creates a discrete closed-loop model programmatically. A MATLAB Function block maintains persistent RLS states, computes the controller, and logs the estimated parameters.

Chapter12_Lesson5_Simulink.m

%% Chapter 12, Lesson 5: Programmatic Simulink Model
% Requires Simulink. The script creates Chapter12_Lesson5_Simulink.slx.
% The adaptive controller is a discrete MATLAB Function block with persistent
% RLS states. The plant implements y(k+1) = 0.82*y(k) + 0.18*u(k).

clear; clc;

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

N = 600;
t = (0:(N - 1))';
r = zeros(N, 1);
r(t >= 40 & t < 200) = 1.0;
r(t >= 200 & t < 340) = -0.5;
r(t >= 340 & t < 470) = 1.5;
r(t >= 470) = 0.25;
r_ts = timeseries(r, t);
assignin('base', 'r_ts', r_ts);

new_system(model);
open_system(model);
set_param(model, ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '1', ...
    'StopTime', num2str(N - 1), ...
    'SaveTime', 'on', ...
    'TimeSaveName', 'tout');

add_block('simulink/Sources/From Workspace', [model '/Reference'], ...
    'VariableName', 'r_ts', ...
    'Position', [40 120 150 150]);

add_block('simulink/Discrete/Discrete Transfer Fcn', [model '/Plant'], ...
    'Numerator', '[0 0.18]', ...
    'Denominator', '[1 -0.82]', ...
    'InitialStates', '0', ...
    'SampleTime', '1', ...
    'Position', [500 80 650 130]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/RLS_STR_Controller'], ...
    'Position', [235 70 410 180]);

add_block('simulink/Sinks/To Workspace', [model '/y_out'], ...
    'VariableName', 'y_out', 'SaveFormat', 'Structure With Time', ...
    'Position', [720 70 830 100]);
add_block('simulink/Sinks/To Workspace', [model '/u_out'], ...
    'VariableName', 'u_out', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 195 610 225]);
add_block('simulink/Sinks/To Workspace', [model '/a_hat_out'], ...
    'VariableName', 'a_hat_out', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 245 610 275]);
add_block('simulink/Sinks/To Workspace', [model '/b_hat_out'], ...
    'VariableName', 'b_hat_out', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 295 610 325]);

% Configure the MATLAB Function block through the Stateflow API.
root = sfroot;
chart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [model '/RLS_STR_Controller']);
if isempty(chart)
    error('Could not locate the MATLAB Function block chart.');
end

chart.Script = sprintf([ ...
'function [u, aHat, bHat] = controller(y, r)\n' ...
'%%#codegen\n' ...
'persistent thetaHat P previousY previousU sampleIndex\n' ...
'if isempty(thetaHat)\n' ...
'    thetaHat = [0.20; 0.40];\n' ...
'    P = 1000.0 * eye(2);\n' ...
'    previousY = 0.0;\n' ...
'    previousU = 0.0;\n' ...
'    sampleIndex = 0.0;\n' ...
'end\n' ...
'lambda = 0.995;\n' ...
'pm = 0.35;\n' ...
'bMin = 0.05;\n' ...
'phi = [previousY; previousU];\n' ...
'denominator = lambda + phi'' * P * phi;\n' ...
'K = (P * phi) / denominator;\n' ...
'predictionError = y - phi'' * thetaHat;\n' ...
'thetaHat = thetaHat + K * predictionError;\n' ...
'P = (P - K * (phi'' * P)) / lambda;\n' ...
'P = 0.5 * (P + P'');\n' ...
'if thetaHat(1) < -0.98, thetaHat(1) = -0.98; end\n' ...
'if thetaHat(1) >  0.98, thetaHat(1) =  0.98; end\n' ...
'if thetaHat(2) <  0.05, thetaHat(2) =  0.05; end\n' ...
'if thetaHat(2) >  0.80, thetaHat(2) =  0.80; end\n' ...
'safeBHat = max(thetaHat(2), bMin);\n' ...
'uCE = ((pm - thetaHat(1)) * y + (1.0 - pm) * r) / safeBHat;\n' ...
'if mod(37.0 * sampleIndex, 23.0) < 11.0\n' ...
'    probe = 0.04;\n' ...
'else\n' ...
'    probe = -0.04;\n' ...
'end\n' ...
'u = min(max(uCE + probe, -4.0), 4.0);\n' ...
'previousY = y;\n' ...
'previousU = u;\n' ...
'sampleIndex = sampleIndex + 1.0;\n' ...
'aHat = thetaHat(1);\n' ...
'bHat = thetaHat(2);\n' ...
'end\n']);

add_line(model, 'Plant/1', 'RLS_STR_Controller/1', 'autorouting', 'on');
add_line(model, 'Reference/1', 'RLS_STR_Controller/2', 'autorouting', 'on');
add_line(model, 'RLS_STR_Controller/1', 'Plant/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'y_out/1', 'autorouting', 'on');
add_line(model, 'RLS_STR_Controller/1', 'u_out/1', 'autorouting', 'on');
add_line(model, 'RLS_STR_Controller/2', 'a_hat_out/1', 'autorouting', 'on');
add_line(model, 'RLS_STR_Controller/3', 'b_hat_out/1', 'autorouting', 'on');

save_system(model);
set_param(model, 'SimulationCommand', 'update');

simOut = sim(model, 'ReturnWorkspaceOutputs', 'on');

% Retrieve logged signals from the SimulationOutput object when available.
yLog = simOut.get('y_out');
uLog = simOut.get('u_out');
aLog = simOut.get('a_hat_out');
bLog = simOut.get('b_hat_out');

figure;
stairs(t, r, 'LineWidth', 1.2); hold on;
plot(yLog.time, yLog.signals.values, 'LineWidth', 1.2);
grid on;
xlabel('sample k'); ylabel('output');
title('Simulink RLS-based STR tracking');
legend('reference r(k)', 'output y(k)', 'Location', 'best');

figure;
plot(aLog.time, aLog.signals.values, 'LineWidth', 1.2); hold on;
plot(bLog.time, bLog.signals.values, 'LineWidth', 1.2);
yline(0.82, '--'); yline(0.18, '--');
grid on;
xlabel('sample k'); ylabel('parameter value');
title('Simulink online parameter estimates');
legend('a hat', 'b hat', 'true a', 'true b', 'Location', 'best');

fprintf('Created and simulated %s.slx\n', model);
fprintf('Final a_hat: %.6f\n', aLog.signals.values(end));
fprintf('Final b_hat: %.6f\n', bLog.signals.values(end));

The notebook contains the following Wolfram Language implementation. It performs the same deterministic loop, creates diagnostic plots, prints the final metrics, and exports a CSV file.

Chapter12_Lesson5.nb

(* Chapter 12, Lesson 5: RLS-based self-tuning regulator example. *)

ClearAll["Global`*"];

n = 600;
aTrue = 0.82;
bTrue = 0.18;
pm = 0.35;
lambda = 0.995;
p = 1000.0 IdentityMatrix[2];
thetaHat = {0.20, 0.40};

uMin = -4.0;
uMax = 4.0;
aMin = -0.98;
aMax = 0.98;
bMin = 0.05;
bMax = 0.80;
probeAmplitude = 0.04;

reference[k_Integer] := Piecewise[{
    {0.0, k < 40},
    {1.0, k < 200},
    {-0.5, k < 340},
    {1.5, k < 470}
}, 0.25];

deterministicProbe[k_Integer] :=
    If[Mod[37 k, 23] < 11, probeAmplitude, -probeAmplitude];

deterministicDisturbance[k_Integer] :=
    0.01 Sin[0.17 k] + 0.005 Cos[0.043 k];

y = ConstantArray[0.0, n + 1];
u = ConstantArray[0.0, n];
r = ConstantArray[0.0, n];
predictionError = ConstantArray[0.0, n];
aHatHistory = ConstantArray[0.0, n + 1];
bHatHistory = ConstantArray[0.0, n + 1];
aHatHistory[[1]] = thetaHat[[1]];
bHatHistory[[1]] = thetaHat[[2]];

Do[
    idx = k + 1;
    r[[idx]] = reference[k];

    safeBHat = Max[thetaHat[[2]], bMin];
    uCE = ((pm - thetaHat[[1]]) y[[idx]]
        + (1.0 - pm) r[[idx]])/safeBHat;
    u[[idx]] = Clip[uCE + deterministicProbe[k], {uMin, uMax}];

    y[[idx + 1]] = aTrue y[[idx]] + bTrue u[[idx]]
        + deterministicDisturbance[k];

    phi = {y[[idx]], u[[idx]]};
    denominator = lambda + phi . p . phi;
    If[!NumericQ[denominator] || denominator <= 0,
        Print["Invalid RLS denominator at sample ", k]; Abort[]
    ];

    gain = p . phi/denominator;
    predictionError[[idx]] = y[[idx + 1]] - phi . thetaHat;
    thetaHat = thetaHat + gain predictionError[[idx]];
    p = (p - Outer[Times, gain, phi . p])/lambda;

    p = 0.5 (p + Transpose[p]);
    thetaHat[[1]] = Clip[thetaHat[[1]], {aMin, aMax}];
    thetaHat[[2]] = Clip[thetaHat[[2]], {bMin, bMax}];

    aHatHistory[[idx + 1]] = thetaHat[[1]];
    bHatHistory[[idx + 1]] = thetaHat[[2]],
    {k, 0, n - 1}
];

trackingRmse = Sqrt[Mean[(y[[51 ;; n]] - r[[51 ;; n]])^2]];
predictionRmse = Sqrt[Mean[predictionError[[51 ;; n]]^2]];

Print["Final a_hat: ", NumberForm[thetaHat[[1]], {8, 6}],
    " (true ", NumberForm[aTrue, {8, 6}], ")"];
Print["Final b_hat: ", NumberForm[thetaHat[[2]], {8, 6}],
    " (true ", NumberForm[bTrue, {8, 6}], ")"];
Print["Tracking RMSE after sample 50: ",
    NumberForm[trackingRmse, {8, 6}]];
Print["Prediction RMSE after sample 50: ",
    NumberForm[predictionRmse, {8, 6}]];

trackingPlot = ListLinePlot[
    {
        Transpose[{Range[0, n - 1], r}],
        Transpose[{Range[0, n], y}]
    },
    PlotLegends -> {"reference r(k)", "output y(k)"},
    PlotLabel -> "RLS-based STR tracking",
    AxesLabel -> {"sample k", "output"},
    GridLines -> Automatic,
    ImageSize -> Large
];

parameterPlot = ListLinePlot[
    {
        Transpose[{Range[0, n], aHatHistory}],
        Transpose[{Range[0, n], bHatHistory}],
        Transpose[{Range[0, n], ConstantArray[aTrue, n + 1]}],
        Transpose[{Range[0, n], ConstantArray[bTrue, n + 1]}]
    },
    PlotLegends -> {"a hat", "b hat", "true a", "true b"},
    PlotLabel -> "Online parameter estimates",
    AxesLabel -> {"sample k", "parameter value"},
    GridLines -> Automatic,
    ImageSize -> Large
];

controlPlot = ListLinePlot[
    Transpose[{Range[0, n - 1], u}],
    PlotLabel -> "Control input",
    AxesLabel -> {"sample k", "u(k)"},
    GridLines -> Automatic,
    ImageSize -> Large
];

Print[trackingPlot];
Print[parameterPlot];
Print[controlPlot];

csvRows = Prepend[
    Table[
        {k, r[[k + 1]], y[[k + 1]], u[[k + 1]],
         aHatHistory[[k + 1]], bHatHistory[[k + 1]],
         predictionError[[k + 1]]},
        {k, 0, n - 1}
    ],
    {"k", "r", "y", "u", "a_hat", "b_hat", "prediction_error"}
];
Export["Chapter12_Lesson5_mathematica_results.csv", csvRows];

• Python: NumPy supplies reliable dense linear algebra, while Matplotlib visualizes estimates, tracking, and input. The python-control package can represent and simulate discrete LTI systems, but the adaptive recursion still needs to be constructed explicitly for this lesson.
• C++: The standard library is sufficient for a two-parameter example. Larger models should use a tested linear-algebra library such as Eigen and a factorized square-root RLS implementation when conditioning is critical.
• Java: The standard library is sufficient here. Larger implementations benefit from a numerical package that provides matrix factorizations and condition-number diagnostics.
• MATLAB: The from-scratch script exposes every equation. In licensed environments, the System Identification Toolbox provides the recursiveLS System object for online least-squares models.
• Simulink: The supplied model uses a MATLAB Function block for transparency. The System Identification Toolbox also provides a Recursive Least Squares Estimator block for models that are linear in their parameters.
• Wolfram Mathematica: built-in list, matrix, plotting, and export functions are sufficient for symbolic inspection and numerical experimentation.

Library implementations reduce coding effort, but students should first be able to map every library input and output to \( \boldsymbol{\phi}(k) \), \( \varepsilon(k+1) \), \( \mathbf{K}(k) \), \( \hat{\boldsymbol{\theta}}(k) \), and \( \mathbf{P}(k) \).

16.1 Parameter Convergence Is Not the Same as Tracking

A regulator can track a constant command while the regressor loses rank. After transients, both \( y(k) \) and \( u(k) \) may approach constants, so repeated regressors provide information mainly in one direction of the two-dimensional parameter space. The output can remain well regulated even though \( \hat{a} \) and \( \hat{b} \) are not individually accurate.

16.2 Prediction Error Must Be Interpreted with Excitation

A small innovation does not by itself prove that every parameter is correct. It may mean that the current input-output trajectory does not distinguish among several parameter vectors. Covariance eigenvalues, regressor richness, and response to reference changes should be inspected together with the innovation.

16.3 Aggressive Desired Poles Increase Control Sensitivity

Choosing \( p_m \) close to zero requests a rapid response. From the controller equation, this may require large gains when \( |\hat{b}| \) is small. Saturation then breaks the exact pole-placement relation and can also bias identification because the applied input differs from an unconstrained design assumption.

16.4 Forgetting Prevents Covariance Collapse but Raises Variance

When \( \lambda<1 \), old information decays and \( \mathbf{P}(k) \) can remain responsive to parameter changes. The cost is larger steady-state variation in estimates under noise. Therefore, the forgetting factor should be selected together with signal scaling, noise level, and expected parameter variation rate.

16.5 Closed-Loop Identification Is Endogenous

The input is generated from measured output and current estimates, so the regression data are not externally prescribed. Noise, estimation, and control are statistically coupled. Classical open-loop least-squares intuition must therefore be applied carefully in adaptive regulation.

Problem 1 (Controller Derivation): For \( y(k+1)=a y(k)+b u(k) \), derive a control law that gives \( y(k+1)=p_m y(k)+(1-p_m)r(k) \) when \( a \) and \( b \) are known. Then write the certainty-equivalent law.

Solution: Equate the plant and desired updates:

\[ a y(k)+b u(k)=p_m y(k)+(1-p_m)r(k). \]

Solving for the input gives

\[ u^{\star}(k)=\frac{p_m-a}{b}y(k)+ \frac{1-p_m}{b}r(k). \]

Replacing the unknown parameters by current estimates yields

\[ u_{\mathrm{ce}}(k)=\frac{p_m-\hat{a}(k)}{\hat{b}(k)}y(k)+ \frac{1-p_m}{\hat{b}(k)}r(k). \]

Problem 2 (One Numerical RLS Update): Let \( \hat{\boldsymbol{\theta}}(k)=[0.5,\;0.25]^T \), \( \mathbf{P}(k)=\operatorname{diag}(10,4) \), \( \boldsymbol{\phi}(k)=[1.2,\;-0.5]^T \), \( y(k+1)=0.9 \), and \( \lambda=0.98 \). Compute the gain, innovation, updated estimate, and covariance.

Solution: First,

\[ \mathbf{P}\boldsymbol{\phi}= \begin{bmatrix}12\\-2\end{bmatrix}, \qquad \lambda+\boldsymbol{\phi}^{T}\mathbf{P}\boldsymbol{\phi} =0.98+15.4=16.38. \]

\[ \mathbf{K}(k)=\frac{1}{16.38} \begin{bmatrix}12\\-2\end{bmatrix} =\begin{bmatrix}0.732601\\-0.122100\end{bmatrix}. \]

The predicted output and innovation are

\[ \hat{y}(k+1\mid k)=1.2(0.5)-0.5(0.25)=0.475, \qquad \varepsilon(k+1)=0.9-0.475=0.425. \]

Therefore,

\[ \hat{\boldsymbol{\theta}}(k+1)= \begin{bmatrix}0.5\\0.25\end{bmatrix} +\begin{bmatrix}0.732601\\-0.122100\end{bmatrix}(0.425) =\begin{bmatrix}0.811355\\0.198107\end{bmatrix}. \]

The covariance becomes approximately

\[ \mathbf{P}(k+1)= \begin{bmatrix} 1.233460 & 1.495104\\ 1.495104 & 3.832449 \end{bmatrix}. \]

Problem 3 (Tracking-Error Dynamics): Assume no saturation or probing. Show that the certainty-equivalent controller gives \( e_c(k+1)=p_m e_c(k)+ \tilde{\boldsymbol{\theta}}^{T}(k)\boldsymbol{\phi}(k)+v(k+1) \).

Solution: Add and subtract the estimated model:

\[ \begin{aligned} y(k+1) &=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}+v(k+1)\\ &=\boldsymbol{\phi}^{T}(k)\hat{\boldsymbol{\theta}}(k) +\boldsymbol{\phi}^{T}(k)\tilde{\boldsymbol{\theta}}(k)+v(k+1). \end{aligned} \]

The controller makes the first term equal to the desired update:

\[ y(k+1)=p_m y(k)+(1-p_m)r(k)+ \tilde{\boldsymbol{\theta}}^{T}(k)\boldsymbol{\phi}(k)+v(k+1). \]

Subtracting \( y_d(k+1)=p_m y_d(k)+(1-p_m)r(k) \) proves the result.

Problem 4 (Forgetting-Factor Memory): For \( \lambda=0.995 \), estimate the effective memory and half-weight age.

Solution:

\[ N_{\mathrm{eff}}\approx\frac{1}{1-0.995}=200. \]

\[ N_{1/2}=\frac{\ln(0.5)}{\ln(0.995)}\approx138.3. \]

Thus, the estimator behaves roughly as though it retained about 200 equally weighted samples, and a sample about 138 steps old has one-half the weight of the newest sample.

Problem 5 (Why Constant Regulation Is Not Persistently Exciting): Explain why a constant command can permit good tracking without uniquely identifying \( a \) and \( b \).

Solution: At steady state,

\[ y(k)\approx\bar{y}, \qquad u(k)\approx\bar{u}, \qquad \boldsymbol{\phi}(k)\approx \begin{bmatrix}\bar{y}\\\bar{u}\end{bmatrix}. \]

Repeating essentially the same regressor contributes an information matrix proportional to \( \boldsymbol{\phi}\boldsymbol{\phi}^{T} \), which has rank one. A two-parameter vector cannot be uniquely recovered from one regressor direction. Nevertheless, the parameter combination \( a\bar{y}+b\bar{u} \) can be predicted accurately, so regulation may remain good. Reference changes or probing create additional regressor directions.

Problem 6 (Input-Gain Protection): Suppose the raw RLS update produces \( \hat{b}(k)=0.002 \) while the known input direction is positive. Explain the control risk and give a safe projected denominator when \( b_{\min}=0.05 \).

Solution: The controller contains terms proportional to \( 1/\hat{b}(k) \). Using 0.002 would amplify the computed input by a factor of 500, causing severe peaking and immediate saturation. With known positive sign, use

\[ \hat{b}_{\mathrm{safe}}(k)= \max\left(\hat{b}(k),b_{\min}\right)=0.05. \]

Parameter projection should also return the stored estimate to the admissible interval. Saturation remains necessary because even a protected denominator does not guarantee that the requested input lies within actuator limits.

This lesson constructed a complete indirect self-tuning regulator for a first-order discrete process. The process was written as a linear regression, exponentially weighted RLS was derived through the matrix inversion lemma, and the current estimates were mapped to a pole-placement controller by certainty equivalence. The resulting tracking-error equation exposed the role of parameter error, while the implementation showed why projection, covariance symmetrization, excitation, saturation, and diagnostics are indispensable. The cross-language benchmark provides a controlled basis for later study of higher-order discrete-time adaptive controllers.