After completing this lesson, students should be able to:

• distinguish plant-parameter adaptation from controller-parameter adaptation;
• derive the certainty-equivalent gain map used by an indirect STR;
• derive a direct controller-gain update from a filtered tracking error;
• prove the ideal asymptotic equivalence of compatible direct and indirect designs;
• explain why their transient behavior is generally not equivalent;
• quantify how plant-estimation covariance is propagated into controller gains;
• select an architecture using stability, computation, interpretability, and constraint criteria.

The lesson builds on the RLS equations introduced in Chapter 12 and on controller-gain parameterization and Lyapunov reasoning from Lessons 2 and 3 of this chapter. Students are assumed to know discrete-time linear systems, pole placement, feedback stability, and state-space or transfer- function representations from linear control.

An indirect STR identifies a plant model and then computes controller coefficients from that model. A direct STR updates the controller coefficients themselves without requiring an explicit online estimate of every physical plant parameter.

flowchart TD
  R1["Reference r(k)"] --> C1["Indirect controller with gains from F(theta_p_hat)"]
  C1 --> P1["Unknown plant"]
  P1 --> Y1["Output y(k)"]
  Y1 --> C1
  Y1 --> ID["Online plant estimator"]
  C1 --> ID
  ID --> SYN["Controller synthesis map F"]
  SYN --> C1

  R2["Reference r(k)"] --> C2["Direct controller with theta_c_hat"]
  C2 --> P2["Unknown plant"]
  P2 --> Y2["Output y(k)"]
  Y2 --> C2
  Y2 --> ERR["Tracking-error model"]
  R2 --> ERR
  ERR --> ADP["Direct gain update"]
  ADP --> C2
        

Let the unknown plant-parameter vector be \( \boldsymbol{\theta}_p \in \mathbb{R}^{n_p} \) and the controller-parameter vector be \( \boldsymbol{\theta}_c \in \mathbb{R}^{n_c} \). The indirect architecture contains an explicit synthesis map \( \mathcal{F} \):

\[ \widehat{\boldsymbol{\theta} }_p(k) \xrightarrow{\text{controller synthesis }\mathcal{F} } \widehat{\boldsymbol{\theta} }_c(k) = \mathcal{F}\!\left(\widehat{\boldsymbol{\theta} }_p(k)\right). \]

The direct architecture instead applies an adaptive law directly in controller coordinates:

\[ \widehat{\boldsymbol{\theta} }_c(k+1) = \mathcal{A}_c\!\left( \widehat{\boldsymbol{\theta} }_c(k),\,\boldsymbol{\phi}_c(k),\,e(k) \right). \]

Therefore, “direct” does not mean that no model assumptions are used. It means that the online adaptive state is expressed in controller coordinates rather than as an explicit plant model.

Consider the uncertain discrete-time plant

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

and the stable reference model

\[ y_m(k+1)=a_m y_m(k)+b_m r(k), \qquad |a_m| < 1. \]

The controller is parameterized directly by two gains:

\[ u(k)=\theta_y(k)y(k)+\theta_r(k)r(k) =\boldsymbol{\theta}_c^{\mathsf T}(k)\boldsymbol{\phi}_c(k), \] \[ \boldsymbol{\theta}_c(k)= \begin{bmatrix}\theta_y(k)\\\theta_r(k)\end{bmatrix},\qquad \boldsymbol{\phi}_c(k)= \begin{bmatrix}y(k)\\r(k)\end{bmatrix}. \]

Substitution into the plant gives the closed-loop equation

\[ y(k+1)=\left(a+b\theta_y(k)\right)y(k) +b\theta_r(k)r(k)+d(k). \]

Exact reference-model matching is achieved by the ideal controller gains

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

These equations are central to the comparison: the indirect STR computes them from plant estimates, whereas the direct STR attempts to approach them without explicitly estimating \(a\) and \(b\) separately.

4.1 Plant regression and RLS

The plant is written as a linear regression:

\[ y(k+1)=\boldsymbol{\phi}_p^{\mathsf T}(k) \boldsymbol{\theta}_p+d(k), \] \[ \boldsymbol{\phi}_p(k)= \begin{bmatrix}y(k)\\u(k)\end{bmatrix},\qquad \boldsymbol{\theta}_p= \begin{bmatrix}a\\b\end{bmatrix}. \]

With forgetting factor \(0 < \lambda \le 1\), RLS is

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

4.2 Certainty-equivalent synthesis map

The controller gains are recomputed as

\[ \widehat{\theta}_y(k)= \frac{a_m-\widehat{a}(k)}{\widehat{b}(k)},\qquad \widehat{\theta}_r(k)= \frac{b_m}{\widehat{b}(k)}. \]

This is called certainty equivalence: the estimated plant is treated as though it were the true plant when the controller is synthesized. In practice, a protected denominator is required:

\[ \widehat{b}_{\mathrm{safe} }(k)= \operatorname{sgn}\!\left(\widehat b(k)\right) \max\!\left(|\widehat b(k)|,b_{\min}\right), \qquad b_{\min} > 0. \]

4.3 General polynomial RST interpretation

For a SISO plant

\[ A(q^{-1})y(k)=q^{-d}B(q^{-1})u(k)+v(k), \]

an RST controller is

\[ R(q^{-1})u(k)=T(q^{-1})r(k)-S(q^{-1})y(k). \]

An indirect STR estimates the coefficients of \(A(q^{-1})\) and \(B(q^{-1})\), then solves a Diophantine equation such as

\[ A(q^{-1})R(q^{-1})+q^{-d}B(q^{-1})S(q^{-1}) =A_c(q^{-1}), \]

where \(A_c(q^{-1})\) is the desired closed-loop characteristic polynomial. Thus an indirect STR separates the online tasks into identification and controller redesign.

Define the tracking error and controller-parameter error as

\[ e(k)=y(k)-y_m(k),\qquad \widetilde{\boldsymbol{\theta} }_c(k)= \boldsymbol{\theta}_c(k)-\boldsymbol{\theta}_c^{\star}. \]

Using the matching identities \(a+b\theta_y^{\star}=a_m\) and \(b\theta_r^{\star}=b_m\), the error dynamics become

\[ e(k+1)=a_m e(k)+b\, \widetilde{\boldsymbol{\theta} }_c^{\mathsf T}(k) \boldsymbol{\phi}_c(k)+d(k). \]

Introduce the filtered adaptation error

\[ \varepsilon(k+1)=e(k+1)-a_m e(k). \]

In the noise-free case,

\[ \varepsilon(k+1)=b\, \widetilde{\boldsymbol{\theta} }_c^{\mathsf T}(k) \boldsymbol{\phi}_c(k). \]

If the control direction \(\operatorname{sgn}(b)\) is known, a normalized direct update is

\[ \boxed{ \boldsymbol{\theta}_c(k+1)= \operatorname{Proj}_{\Omega_c}\!\left[ \boldsymbol{\theta}_c(k)- \gamma\,\operatorname{sgn}(b) \frac{\boldsymbol{\phi}_c(k)\varepsilon(k+1)} {\delta+\boldsymbol{\phi}_c^{\mathsf T}(k)\boldsymbol{\phi}_c(k)} \right]}, \]

where \(\delta > 0\) prevents division by a small regressor norm and \(\Omega_c\) is an admissible set of controller gains. The magnitude of \(b\) is not explicitly estimated; it is absorbed into the effective adaptation rate.

5.1 A one-step parameter-energy argument

Ignore projection and disturbance temporarily. Let \(s(k)=\boldsymbol{\phi}_c^{\mathsf T}(k) \boldsymbol{\phi}_c(k)\). Then

\[ \widetilde{\boldsymbol{\theta} }_c(k+1)= \left[\mathbf{I}- \frac{\gamma|b|}{\delta+s(k)} \boldsymbol{\phi}_c(k)\boldsymbol{\phi}_c^{\mathsf T}(k) \right] \widetilde{\boldsymbol{\theta} }_c(k). \]

Therefore, with \(z(k)=\boldsymbol{\phi}_c^{\mathsf T}(k) \widetilde{\boldsymbol{\theta} }_c(k)\),

\[ \begin{aligned} \|\widetilde{\boldsymbol{\theta} }_c(k+1)\|^2 -\|\widetilde{\boldsymbol{\theta} }_c(k)\|^2 &=-\alpha(k)\left[2-\alpha(k)s(k)\right]z^2(k),\\ \alpha(k)&=\frac{\gamma|b|}{\delta+s(k)}. \end{aligned} \]

Since \(0 \le \alpha(k)s(k) < \gamma|b|\), the parameter-error norm is nonincreasing whenever \(0 < \gamma|b| < 2\). Persistent excitation is still required to force every component of the parameter error to zero; without it, only the component observed through the regressor is reduced.

6.1 Proposition: common ideal controller

Suppose the plant belongs exactly to the assumed model class, \(b \neq 0\), and both architectures use the same controller structure and reference model. Then the ideal direct gains are exactly the certainty-equivalent gains evaluated at the true plant:

\[ \boldsymbol{\theta}_c^{\star} =\mathcal{F}(\boldsymbol{\theta}_p) =\begin{bmatrix} (a_m-a)/b\\ b_m/b \end{bmatrix}. \]

6.2 Theorem: asymptotic closed-loop equivalence

Assume:

1. \(\widehat{\boldsymbol{\theta} }_p(k)\to \boldsymbol{\theta}_p\) for the indirect estimator;
2. \(\mathcal{F}\) is continuous at the true plant parameters;
3. \(\widehat{\boldsymbol{\theta} }_c(k)\to \boldsymbol{\theta}_c^{\star}\) for the direct law;
4. the limiting reference-model dynamics are stable.

Then

\[ \lim_{k\to\infty} \widehat{\boldsymbol{\theta} }_{c,\mathrm{ind} }(k) =\lim_{k\to\infty} \widehat{\boldsymbol{\theta} }_{c,\mathrm{dir} }(k) =\boldsymbol{\theta}_c^{\star}, \]

and both limiting closed loops reproduce the same reference-model coefficients.

Proof. By continuity of the synthesis map,

\[ \widehat{\boldsymbol{\theta} }_{c,\mathrm{ind} }(k) =\mathcal{F}\!\left(\widehat{\boldsymbol{\theta} }_p(k)\right) \to\mathcal{F}(\boldsymbol{\theta}_p) =\boldsymbol{\theta}_c^{\star}. \]

The direct convergence assumption gives the same limit. Substitution of either limiting gain vector into the plant yields \(a+b\theta_y^{\star}=a_m\) and \(b\theta_r^{\star}=b_m\). Hence the limiting nominal tracking-error equation is \(e(k+1)=a_m e(k)\), whose solution converges to zero because \(|a_m| < 1\). This proves the result. \(\square\)

6.3 Why the trajectories are generally different

At finite time, the indirect gain increment is induced by a nonlinear coordinate transformation:

\[ \Delta\widehat{\boldsymbol{\theta} }_{c,\mathrm{ind} }(k) =\mathcal{F}\!\left( \widehat{\boldsymbol{\theta} }_p(k)+ \Delta\widehat{\boldsymbol{\theta} }_p(k)\right) -\mathcal{F}\!\left(\widehat{\boldsymbol{\theta} }_p(k)\right). \]

The direct increment is generated by a controller-space adaptation law. Except in special cases where that law is deliberately chosen as the transformed plant-estimator update, the two increments are unequal. Consequently, equal asymptotic targets do not imply equal transients, equal robustness, or equal control effort.

The first-order certainty-equivalent map is

\[ \mathcal{F}(a,b)= \begin{bmatrix}(a_m-a)/b\\b_m/b\end{bmatrix}. \]

Its Jacobian is

\[ \mathbf{J}_{\mathcal{F} }(a,b)= \frac{\partial\boldsymbol{\theta}_c}{\partial\boldsymbol{\theta}_p} =\begin{bmatrix} -1/b & -(a_m-a)/b^2\\ 0 & -b_m/b^2 \end{bmatrix}. \]

For a small plant-estimation error \(\Delta\boldsymbol{\theta}_p\), first-order propagation gives

\[ \Delta\boldsymbol{\theta}_c \approx\mathbf{J}_{\mathcal{F} } \Delta\boldsymbol{\theta}_p. \]

In scalar form,

\[ \Delta\theta_y\approx- rac{\Delta a}{b} -\frac{a_m-a}{b^2}\Delta b,\qquad \Delta\theta_r\approx- rac{b_m}{b^2}\Delta b. \]

Thus uncertainty in the input gain is amplified proportionally to \(1/b^2\). This explains why an indirect STR can exhibit large gain excursions when \(\widehat b(k)\) passes near zero, even if the plant estimate itself changes only modestly.

7.1 Covariance propagation

If the plant estimate has local covariance \(\mathbf{P}_p(k)\), the delta method gives the approximate controller-gain covariance

\[ \boxed{\mathbf{P}_c(k)\approx \mathbf{J}_{\mathcal{F} } \mathbf{P}_p(k) \mathbf{J}_{\mathcal{F} }^{\mathsf T} }. \]

A direct STR avoids this particular nonlinear covariance transformation, but it is not noise-free. Its controller regressor is generated inside the feedback loop, so regressor-noise correlation and unmodeled dynamics can bias the direct adaptation error. Direct adaptation removes one sensitivity mechanism; it does not remove all statistical difficulties.

8.1 Indirect STR

The indirect loop combines an estimator and a time-varying controller synthesis map. Convergence of the plant estimate is useful but does not, by itself, guarantee safe transient control. During adaptation, the estimated model can generate unstable or excessively aggressive controllers. Typical safeguards include:

• projection of estimated plant parameters onto a physically admissible set;
• lower bounds on estimated high-frequency gain magnitude;
• closed-loop pole checks before accepting a newly synthesized controller;
• rate limits or filtering on controller-coefficient changes;
• a fixed fallback controller during estimator loss of confidence.

8.2 Direct STR

A direct law can often be designed together with an error model so that a Lyapunov or parameter-energy argument constrains the update. This can make the stability proof structurally shorter because no separate plant-to-controller map appears. However, the proof normally depends on matching conditions, known control direction, suitable normalization, and a controller parameterization rich enough to contain the ideal controller.

8.3 Tracking convergence versus parameter convergence

For either architecture, small tracking error does not imply unique parameter convergence. If the reference does not excite all regressor directions, many parameter vectors can produce similar closed-loop behavior. Persistent excitation is needed when exact convergence of all estimated coefficients is required.

\[ \exists\,N\in\mathbb{N},\;\alpha > 0:\qquad \sum_{i=k}^{k+N-1}\boldsymbol{\phi}(i) \boldsymbol{\phi}^{\mathsf T}(i)\succeq\alpha\mathbf{I} \quad\forall k. \]

The relevant regressor is different: an indirect STR needs excitation of the plant-identification regressor, whereas a direct STR needs excitation of the controller-error regressor.

flowchart TD
  S["Start with adaptive-control objective"] --> M["Is an explicit online plant model needed for monitoring or redesign?"]
  M -->|"yes"| I["Prefer indirect STR"]
  M -->|"no"| C["Can the ideal controller be \nlinearly parameterized with \nknown control direction?"]
  C -->|"yes"| D["Direct STR is a strong candidate"]
  C -->|"no"| I2["Use indirect STR or \nredesign the parameterization"]
  I --> G["Protect synthesis map, \nvalidate poles, and \nkeep fallback control"]
  D --> P["Use normalization, projection, \nsaturation handling, \nand stability proof"]
  I2 --> G
  G --> V["Validate excitation, noise sensitivity, and plant-change transients"]
  P --> V
        

A useful engineering rule is: choose indirect adaptation when the model is valuable beyond the immediate controller, and choose direct adaptation when a well-justified controller parameterization and error model are available and implementation simplicity is dominant.

The code examples use the piecewise plant

\[ (a,b)= \begin{cases} (0.65,0.55), & 0\le k < 400,\\ (0.78,0.42), & 400\le k < 800, \end{cases} \]

with the reference model

\[ y_m(k+1)=0.35y_m(k)+0.65r(k). \]

The indirect STR uses two-parameter RLS with \(\lambda=0.995\), denominator protection, controller-gain projection, and input saturation. The direct STR uses the filtered error \(\varepsilon(k)=e(k)-0.35e(k-1)\), normalized gradient adaptation, projection, and the same input saturation.

With the deterministic excitation and disturbance in the supplied code, the implementations produce approximately the following metrics:

These numbers demonstrate one tuning, not a universal ranking. Different adaptation gains, forgetting factors, excitation spectra, disturbances, constraints, or plant structures can reverse the result. The experiment is intended to make the architectural pathways observable and comparable.

Libraries: NumPy performs vector and matrix operations; Matplotlib is optional for visualization. The estimator and adaptive laws are implemented from scratch rather than through a specialized adaptive- control package.

Chapter13_Lesson4.py

"""Chapter 13, Lesson 4: Direct versus indirect self-tuning regulators.

The program compares two adaptive controllers on the same first-order,
piecewise-constant discrete-time plant:

    y(k+1) = a(k) y(k) + b(k) u(k) + d(k)

Indirect STR:
    1. Estimate [a, b] with recursive least squares.
    2. Compute certainty-equivalent gains from the reference-model matching map.

Direct STR:
    1. Parameterize u(k) = theta_y(k) y(k) + theta_r(k) r(k).
    2. Update the controller gains directly from a filtered tracking error.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Dict

import numpy as np


@dataclass
class SimulationResult:
    reference: np.ndarray
    model_output: np.ndarray
    indirect_output: np.ndarray
    direct_output: np.ndarray
    indirect_control: np.ndarray
    direct_control: np.ndarray
    plant_estimates: np.ndarray
    indirect_gains: np.ndarray
    direct_gains: np.ndarray


def clipped(value: float, lower: float, upper: float) -> float:
    """Return value projected onto the closed interval [lower, upper]."""
    return float(np.minimum(np.maximum(value, lower), upper))


def rmse(signal: np.ndarray) -> float:
    """Root-mean-square value."""
    return float(np.sqrt(np.mean(np.square(signal))))


def build_reference(samples: int) -> np.ndarray:
    """Piecewise command plus deterministic excitation for identifiability."""
    k = np.arange(samples)
    reference = np.select(
        [k < 120, k < 300, k < 460, k < 620],
        [0.0, 1.0, -0.6, 0.8],
        default=-0.2,
    ).astype(float)
    reference += 0.08 * np.sin(0.071 * k) + 0.04 * np.sin(0.19 * k)
    return reference


def simulate(samples: int = 800) -> SimulationResult:
    """Simulate direct and indirect STR architectures."""
    if samples < 10:
        raise ValueError("samples must be at least 10")

    k_axis = np.arange(samples)
    reference = build_reference(samples)

    # The plant changes at k = samples/2 to test tracking of new dynamics.
    change_index = samples // 2
    plant_a = np.where(k_axis < change_index, 0.65, 0.78)
    plant_b = np.where(k_axis < change_index, 0.55, 0.42)
    disturbance = 0.012 * np.sin(0.37 * k_axis) + 0.006 * np.cos(0.11 * k_axis)

    # Stable reference model with unit steady-state gain.
    a_m = 0.35
    b_m = 1.0 - a_m

    model_output = np.zeros(samples + 1)
    indirect_output = np.zeros(samples + 1)
    direct_output = np.zeros(samples + 1)
    indirect_control = np.zeros(samples)
    direct_control = np.zeros(samples)

    plant_estimates = np.zeros((samples, 2))
    indirect_gains = np.zeros((samples, 2))
    direct_gains = np.zeros((samples, 2))

    # Indirect STR state: RLS estimate theta_p = [a_hat, b_hat].
    theta_p = np.array([0.0, 0.5], dtype=float)
    covariance = 100.0 * np.eye(2)
    forgetting_factor = 0.995
    minimum_input_gain = 0.08

    # Direct STR state: theta_c = [theta_y, theta_r].
    theta_c = np.array([0.0, 1.0], dtype=float)
    direct_step = 1.2
    normalization = 0.2

    for k in range(samples):
        model_output[k + 1] = a_m * model_output[k] + b_m * reference[k]

        # ----------------------- Indirect STR -----------------------
        # At time k, update the plant model using the completed transition
        # y(k-1), u(k-1) -> y(k).
        if k > 0:
            phi_p = np.array([indirect_output[k - 1], indirect_control[k - 1]])
            denominator = forgetting_factor + phi_p @ covariance @ phi_p
            kalman_gain = covariance @ phi_p / denominator
            prediction_error = indirect_output[k] - phi_p @ theta_p
            theta_p = theta_p + kalman_gain * prediction_error
            covariance = (
                covariance - np.outer(kalman_gain, phi_p) @ covariance
            ) / forgetting_factor
            covariance = 0.5 * (covariance + covariance.T)

        a_hat, b_hat = theta_p
        # The known control direction is positive. The lower bound prevents the
        # certainty-equivalent gain map from dividing by a value near zero.
        b_safe = max(abs(b_hat), minimum_input_gain)
        gain_y = clipped((a_m - a_hat) / b_safe, -4.0, 4.0)
        gain_r = clipped(b_m / b_safe, 0.0, 4.0)

        indirect_control[k] = clipped(
            gain_y * indirect_output[k] + gain_r * reference[k], -5.0, 5.0
        )
        indirect_output[k + 1] = (
            plant_a[k] * indirect_output[k]
            + plant_b[k] * indirect_control[k]
            + disturbance[k]
        )

        plant_estimates[k, :] = theta_p
        indirect_gains[k, :] = [gain_y, gain_r]

        # ------------------------- Direct STR -----------------------
        # From e(k) = a_m e(k-1) + b * theta_tilde^T phi(k-1) + d(k-1),
        # use epsilon(k) = e(k) - a_m e(k-1) as the direct adaptation error.
        # Only the sign of b is used; b itself is not estimated.
        if k > 0:
            error_k = direct_output[k] - model_output[k]
            error_previous = direct_output[k - 1] - model_output[k - 1]
            filtered_error = error_k - a_m * error_previous
            phi_c = np.array([direct_output[k - 1], reference[k - 1]])
            denominator = normalization + phi_c @ phi_c
            theta_c = theta_c - direct_step * phi_c * filtered_error / denominator

            # Projection is a practical safeguard against gain drift.
            theta_c[0] = clipped(theta_c[0], -4.0, 4.0)
            theta_c[1] = clipped(theta_c[1], 0.0, 4.0)

        direct_control[k] = clipped(
            theta_c @ np.array([direct_output[k], reference[k]]), -5.0, 5.0
        )
        direct_output[k + 1] = (
            plant_a[k] * direct_output[k]
            + plant_b[k] * direct_control[k]
            + disturbance[k]
        )
        direct_gains[k, :] = theta_c

    return SimulationResult(
        reference=reference,
        model_output=model_output,
        indirect_output=indirect_output,
        direct_output=direct_output,
        indirect_control=indirect_control,
        direct_control=direct_control,
        plant_estimates=plant_estimates,
        indirect_gains=indirect_gains,
        direct_gains=direct_gains,
    )


def summarize(result: SimulationResult) -> Dict[str, float]:
    """Compute comparable performance indices."""
    samples = result.reference.size
    midpoint = samples // 2
    indirect_error = result.indirect_output[1:] - result.model_output[1:]
    direct_error = result.direct_output[1:] - result.model_output[1:]

    metrics = {
        "indirect_rmse_all": rmse(indirect_error),
        "direct_rmse_all": rmse(direct_error),
        "indirect_rmse_before_change": rmse(indirect_error[:midpoint]),
        "direct_rmse_before_change": rmse(direct_error[:midpoint]),
        "indirect_rmse_after_change": rmse(indirect_error[midpoint:]),
        "direct_rmse_after_change": rmse(direct_error[midpoint:]),
        "indirect_control_energy": float(np.mean(result.indirect_control**2)),
        "direct_control_energy": float(np.mean(result.direct_control**2)),
    }
    return metrics


def main() -> None:
    result = simulate()
    metrics = summarize(result)

    print("Direct versus indirect STR comparison")
    print("-" * 45)
    for name, value in metrics.items():
        print(f"{name:36s}: {value:.6f}")

    print("\nFinal indirect plant estimate [a_hat, b_hat]:")
    print(result.plant_estimates[-1])
    print("Final indirect gains [gain_y, gain_r]:")
    print(result.indirect_gains[-1])
    print("Final direct gains [theta_y, theta_r]:")
    print(result.direct_gains[-1])

    # Plotting is optional so that the simulation still runs in a minimal
    # NumPy-only environment.
    try:
        import matplotlib.pyplot as plt

        k = np.arange(result.reference.size)

        plt.figure(figsize=(10, 5))
        plt.plot(k, result.reference, "--", label="reference")
        plt.plot(k, result.model_output[1:], label="reference model")
        plt.plot(k, result.indirect_output[1:], label="indirect STR")
        plt.plot(k, result.direct_output[1:], label="direct STR")
        plt.xlabel("sample k")
        plt.ylabel("output")
        plt.title("Tracking comparison")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        plt.figure(figsize=(10, 5))
        plt.plot(k, result.indirect_gains[:, 0], label="indirect gain_y")
        plt.plot(k, result.indirect_gains[:, 1], label="indirect gain_r")
        plt.plot(k, result.direct_gains[:, 0], "--", label="direct theta_y")
        plt.plot(k, result.direct_gains[:, 1], "--", label="direct theta_r")
        plt.xlabel("sample k")
        plt.ylabel("controller gain")
        plt.title("Controller-parameter trajectories")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        plt.figure(figsize=(10, 5))
        plt.plot(k, result.plant_estimates[:, 0], label="a_hat")
        plt.plot(k, result.plant_estimates[:, 1], label="b_hat")
        plt.axvline(result.reference.size // 2, linestyle="--", label="plant change")
        plt.xlabel("sample k")
        plt.ylabel("estimated plant parameter")
        plt.title("Plant estimates used only by the indirect STR")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()
    except ImportError:
        print("\nMatplotlib is not installed; numerical simulation completed without plots.")


if __name__ == "__main__":
    main()

Libraries: only the ISO C++ standard library is used. The two-dimensional RLS covariance update is written explicitly, making the implementation suitable for later migration to Eigen, Armadillo, or an embedded target.

Chapter13_Lesson4.cpp

// Chapter 13, Lesson 4: Direct versus indirect self-tuning regulators.
// Standard-library-only C++17 implementation.

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

namespace {

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

double referenceAt(int k) {
    double base = -0.2;
    if (k < 120) {
        base = 0.0;
    } else if (k < 300) {
        base = 1.0;
    } else if (k < 460) {
        base = -0.6;
    } else if (k < 620) {
        base = 0.8;
    }
    return base + 0.08 * std::sin(0.071 * k) + 0.04 * std::sin(0.19 * k);
}

double rmse(const std::vector<double>& error, int begin, int end) {
    if (begin < 0 || end > static_cast<int>(error.size()) || begin >= end) {
        throw std::invalid_argument("invalid RMSE interval");
    }
    double sum = 0.0;
    for (int i = begin; i < end; ++i) {
        sum += error[i] * error[i];
    }
    return std::sqrt(sum / static_cast<double>(end - begin));
}

}  // namespace

int main() {
    constexpr int samples = 800;
    constexpr int changeIndex = samples / 2;
    constexpr double aModel = 0.35;
    constexpr double bModel = 1.0 - aModel;

    std::vector<double> reference(samples, 0.0);
    std::vector<double> modelOutput(samples + 1, 0.0);
    std::vector<double> indirectOutput(samples + 1, 0.0);
    std::vector<double> directOutput(samples + 1, 0.0);
    std::vector<double> indirectControl(samples, 0.0);
    std::vector<double> directControl(samples, 0.0);

    std::vector<std::array<double, 2>> plantEstimates(samples);
    std::vector<std::array<double, 2>> indirectGains(samples);
    std::vector<std::array<double, 2>> directGains(samples);

    // Indirect STR: thetaPlant = [aHat, bHat], P is the RLS covariance matrix.
    std::array<double, 2> thetaPlant{0.0, 0.5};
    std::array<std::array<double, 2>, 2> P{ { {100.0, 0.0}, {0.0, 100.0} } };
    constexpr double forgettingFactor = 0.995;
    constexpr double minimumInputGain = 0.08;

    // Direct STR: thetaController = [thetaY, thetaR].
    std::array<double, 2> thetaController{0.0, 1.0};
    constexpr double directStep = 1.2;
    constexpr double normalization = 0.2;

    for (int k = 0; k < samples; ++k) {
        reference[k] = referenceAt(k);
        const double plantA = (k < changeIndex) ? 0.65 : 0.78;
        const double plantB = (k < changeIndex) ? 0.55 : 0.42;
        const double disturbance =
            0.012 * std::sin(0.37 * k) + 0.006 * std::cos(0.11 * k);

        modelOutput[k + 1] = aModel * modelOutput[k] + bModel * reference[k];

        // ----------------------- Indirect STR -----------------------
        if (k > 0) {
            const std::array<double, 2> phi{
                indirectOutput[k - 1], indirectControl[k - 1]};

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

            const double prediction =
                phi[0] * thetaPlant[0] + phi[1] * thetaPlant[1];
            const double predictionError = indirectOutput[k] - prediction;
            thetaPlant[0] += gain[0] * predictionError;
            thetaPlant[1] += gain[1] * predictionError;

            const double phiTP0 = phi[0] * P[0][0] + phi[1] * P[1][0];
            const double phiTP1 = phi[0] * P[0][1] + phi[1] * P[1][1];
            std::array<std::array<double, 2>, 2> newP{};
            newP[0][0] = (P[0][0] - gain[0] * phiTP0) / forgettingFactor;
            newP[0][1] = (P[0][1] - gain[0] * phiTP1) / forgettingFactor;
            newP[1][0] = (P[1][0] - gain[1] * phiTP0) / forgettingFactor;
            newP[1][1] = (P[1][1] - gain[1] * phiTP1) / forgettingFactor;

            // Symmetrize to suppress roundoff asymmetry.
            const double offDiagonal = 0.5 * (newP[0][1] + newP[1][0]);
            newP[0][1] = offDiagonal;
            newP[1][0] = offDiagonal;
            P = newP;
        }

        const double aHat = thetaPlant[0];
        const double bHatSafe = std::max(std::abs(thetaPlant[1]), minimumInputGain);
        const double gainY = clipped((aModel - aHat) / bHatSafe, -4.0, 4.0);
        const double gainR = clipped(bModel / bHatSafe, 0.0, 4.0);

        indirectControl[k] = clipped(
            gainY * indirectOutput[k] + gainR * reference[k], -5.0, 5.0);
        indirectOutput[k + 1] = plantA * indirectOutput[k]
                              + plantB * indirectControl[k]
                              + disturbance;

        plantEstimates[k] = thetaPlant;
        indirectGains[k] = {gainY, gainR};

        // ------------------------- Direct STR -----------------------
        if (k > 0) {
            const double errorK = directOutput[k] - modelOutput[k];
            const double errorPrevious = directOutput[k - 1] - modelOutput[k - 1];
            const double filteredError = errorK - aModel * errorPrevious;
            const std::array<double, 2> phi{
                directOutput[k - 1], reference[k - 1]};
            const double denominator =
                normalization + phi[0] * phi[0] + phi[1] * phi[1];

            thetaController[0] -= directStep * phi[0] * filteredError / denominator;
            thetaController[1] -= directStep * phi[1] * filteredError / denominator;
            thetaController[0] = clipped(thetaController[0], -4.0, 4.0);
            thetaController[1] = clipped(thetaController[1], 0.0, 4.0);
        }

        directControl[k] = clipped(
            thetaController[0] * directOutput[k]
                + thetaController[1] * reference[k],
            -5.0,
            5.0);
        directOutput[k + 1] = plantA * directOutput[k]
                            + plantB * directControl[k]
                            + disturbance;
        directGains[k] = thetaController;
    }

    std::vector<double> indirectError(samples);
    std::vector<double> directError(samples);
    for (int k = 0; k < samples; ++k) {
        indirectError[k] = indirectOutput[k + 1] - modelOutput[k + 1];
        directError[k] = directOutput[k + 1] - modelOutput[k + 1];
    }

    const double indirectEnergy = std::inner_product(
        indirectControl.begin(), indirectControl.end(), indirectControl.begin(), 0.0)
        / samples;
    const double directEnergy = std::inner_product(
        directControl.begin(), directControl.end(), directControl.begin(), 0.0)
        / samples;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Direct versus indirect STR comparison\n";
    std::cout << "---------------------------------------------\n";
    std::cout << "Indirect RMSE, all samples      : "
              << rmse(indirectError, 0, samples) << '\n';
    std::cout << "Direct RMSE, all samples        : "
              << rmse(directError, 0, samples) << '\n';
    std::cout << "Indirect RMSE, before change    : "
              << rmse(indirectError, 0, changeIndex) << '\n';
    std::cout << "Direct RMSE, before change      : "
              << rmse(directError, 0, changeIndex) << '\n';
    std::cout << "Indirect RMSE, after change     : "
              << rmse(indirectError, changeIndex, samples) << '\n';
    std::cout << "Direct RMSE, after change       : "
              << rmse(directError, changeIndex, samples) << '\n';
    std::cout << "Indirect mean control energy    : " << indirectEnergy << '\n';
    std::cout << "Direct mean control energy      : " << directEnergy << '\n';
    std::cout << "Final plant estimate [a,b]      : ["
              << plantEstimates.back()[0] << ", " << plantEstimates.back()[1] << "]\n";
    std::cout << "Final indirect gains [y,r]      : ["
              << indirectGains.back()[0] << ", " << indirectGains.back()[1] << "]\n";
    std::cout << "Final direct gains [y,r]        : ["
              << directGains.back()[0] << ", " << directGains.back()[1] << "]\n";

    std::ofstream csv("Chapter13_Lesson4_results.csv");
    if (!csv) {
        std::cerr << "Warning: could not create CSV output.\n";
        return 0;
    }
    csv << "k,reference,model_output,indirect_output,direct_output,"
           "indirect_control,direct_control,a_hat,b_hat,"
           "indirect_gain_y,indirect_gain_r,direct_gain_y,direct_gain_r\n";
    for (int k = 0; k < samples; ++k) {
        csv << k << ',' << reference[k] << ',' << modelOutput[k + 1] << ','
            << indirectOutput[k + 1] << ',' << directOutput[k + 1] << ','
            << indirectControl[k] << ',' << directControl[k] << ','
            << plantEstimates[k][0] << ',' << plantEstimates[k][1] << ','
            << indirectGains[k][0] << ',' << indirectGains[k][1] << ','
            << directGains[k][0] << ',' << directGains[k][1] << '\n';
    }

    return 0;
}

Libraries: the Java standard library is sufficient. For larger adaptive-control projects, EJML or Apache Commons Math can replace the explicit small-matrix arithmetic.

Chapter13_Lesson4.java

// Chapter 13, Lesson 4: Direct versus indirect self-tuning regulators.
// Standard-library-only Java implementation.

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 Chapter13_Lesson4 {
    private static final int SAMPLES = 800;
    private static final int CHANGE_INDEX = SAMPLES / 2;
    private static final double A_MODEL = 0.35;
    private static final double B_MODEL = 1.0 - A_MODEL;

    private Chapter13_Lesson4() {
        // Utility class.
    }

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

    private static double referenceAt(int k) {
        double base = -0.2;
        if (k < 120) {
            base = 0.0;
        } else if (k < 300) {
            base = 1.0;
        } else if (k < 460) {
            base = -0.6;
        } else if (k < 620) {
            base = 0.8;
        }
        return base + 0.08 * Math.sin(0.071 * k) + 0.04 * Math.sin(0.19 * k);
    }

    private static double rmse(double[] error, int begin, int end) {
        if (begin < 0 || end > error.length || begin >= end) {
            throw new IllegalArgumentException("invalid RMSE interval");
        }
        double sum = 0.0;
        for (int i = begin; i < end; i++) {
            sum += error[i] * error[i];
        }
        return Math.sqrt(sum / (end - begin));
    }

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

        double[] reference = new double[SAMPLES];
        double[] modelOutput = new double[SAMPLES + 1];
        double[] indirectOutput = new double[SAMPLES + 1];
        double[] directOutput = new double[SAMPLES + 1];
        double[] indirectControl = new double[SAMPLES];
        double[] directControl = new double[SAMPLES];

        double[][] plantEstimates = new double[SAMPLES][2];
        double[][] indirectGains = new double[SAMPLES][2];
        double[][] directGains = new double[SAMPLES][2];

        // Indirect STR state.
        double[] thetaPlant = {0.0, 0.5};
        double[][] covariance = { {100.0, 0.0}, {0.0, 100.0} };
        final double forgettingFactor = 0.995;
        final double minimumInputGain = 0.08;

        // Direct STR state.
        double[] thetaController = {0.0, 1.0};
        final double directStep = 1.2;
        final double normalization = 0.2;

        for (int k = 0; k < SAMPLES; k++) {
            reference[k] = referenceAt(k);
            double plantA = (k < CHANGE_INDEX) ? 0.65 : 0.78;
            double plantB = (k < CHANGE_INDEX) ? 0.55 : 0.42;
            double disturbance =
                    0.012 * Math.sin(0.37 * k) + 0.006 * Math.cos(0.11 * k);

            modelOutput[k + 1] = A_MODEL * modelOutput[k] + B_MODEL * reference[k];

            // ----------------------- Indirect STR -----------------------
            if (k > 0) {
                double phi0 = indirectOutput[k - 1];
                double phi1 = indirectControl[k - 1];
                double pPhi0 = covariance[0][0] * phi0 + covariance[0][1] * phi1;
                double pPhi1 = covariance[1][0] * phi0 + covariance[1][1] * phi1;
                double denominator = forgettingFactor + phi0 * pPhi0 + phi1 * pPhi1;
                double gain0 = pPhi0 / denominator;
                double gain1 = pPhi1 / denominator;

                double prediction = phi0 * thetaPlant[0] + phi1 * thetaPlant[1];
                double predictionError = indirectOutput[k] - prediction;
                thetaPlant[0] += gain0 * predictionError;
                thetaPlant[1] += gain1 * predictionError;

                double phiTP0 = phi0 * covariance[0][0] + phi1 * covariance[1][0];
                double phiTP1 = phi0 * covariance[0][1] + phi1 * covariance[1][1];
                double[][] nextCovariance = new double[2][2];
                nextCovariance[0][0] =
                        (covariance[0][0] - gain0 * phiTP0) / forgettingFactor;
                nextCovariance[0][1] =
                        (covariance[0][1] - gain0 * phiTP1) / forgettingFactor;
                nextCovariance[1][0] =
                        (covariance[1][0] - gain1 * phiTP0) / forgettingFactor;
                nextCovariance[1][1] =
                        (covariance[1][1] - gain1 * phiTP1) / forgettingFactor;

                double offDiagonal =
                        0.5 * (nextCovariance[0][1] + nextCovariance[1][0]);
                nextCovariance[0][1] = offDiagonal;
                nextCovariance[1][0] = offDiagonal;
                covariance = nextCovariance;
            }

            double aHat = thetaPlant[0];
            double bHatSafe = Math.max(Math.abs(thetaPlant[1]), minimumInputGain);
            double gainY = clipped((A_MODEL - aHat) / bHatSafe, -4.0, 4.0);
            double gainR = clipped(B_MODEL / bHatSafe, 0.0, 4.0);

            indirectControl[k] = clipped(
                    gainY * indirectOutput[k] + gainR * reference[k], -5.0, 5.0);
            indirectOutput[k + 1] = plantA * indirectOutput[k]
                    + plantB * indirectControl[k]
                    + disturbance;

            plantEstimates[k][0] = thetaPlant[0];
            plantEstimates[k][1] = thetaPlant[1];
            indirectGains[k][0] = gainY;
            indirectGains[k][1] = gainR;

            // ------------------------- Direct STR -----------------------
            if (k > 0) {
                double errorK = directOutput[k] - modelOutput[k];
                double errorPrevious = directOutput[k - 1] - modelOutput[k - 1];
                double filteredError = errorK - A_MODEL * errorPrevious;
                double phi0 = directOutput[k - 1];
                double phi1 = reference[k - 1];
                double denominator = normalization + phi0 * phi0 + phi1 * phi1;

                thetaController[0] -= directStep * phi0 * filteredError / denominator;
                thetaController[1] -= directStep * phi1 * filteredError / denominator;
                thetaController[0] = clipped(thetaController[0], -4.0, 4.0);
                thetaController[1] = clipped(thetaController[1], 0.0, 4.0);
            }

            directControl[k] = clipped(
                    thetaController[0] * directOutput[k]
                            + thetaController[1] * reference[k],
                    -5.0,
                    5.0);
            directOutput[k + 1] = plantA * directOutput[k]
                    + plantB * directControl[k]
                    + disturbance;
            directGains[k][0] = thetaController[0];
            directGains[k][1] = thetaController[1];
        }

        double[] indirectError = new double[SAMPLES];
        double[] directError = new double[SAMPLES];
        double indirectEnergy = 0.0;
        double directEnergy = 0.0;
        for (int k = 0; k < SAMPLES; k++) {
            indirectError[k] = indirectOutput[k + 1] - modelOutput[k + 1];
            directError[k] = directOutput[k + 1] - modelOutput[k + 1];
            indirectEnergy += indirectControl[k] * indirectControl[k];
            directEnergy += directControl[k] * directControl[k];
        }
        indirectEnergy /= SAMPLES;
        directEnergy /= SAMPLES;

        System.out.printf("Direct versus indirect STR comparison%n");
        System.out.printf("---------------------------------------------%n");
        System.out.printf("Indirect RMSE, all samples      : %.6f%n",
                rmse(indirectError, 0, SAMPLES));
        System.out.printf("Direct RMSE, all samples        : %.6f%n",
                rmse(directError, 0, SAMPLES));
        System.out.printf("Indirect RMSE, before change    : %.6f%n",
                rmse(indirectError, 0, CHANGE_INDEX));
        System.out.printf("Direct RMSE, before change      : %.6f%n",
                rmse(directError, 0, CHANGE_INDEX));
        System.out.printf("Indirect RMSE, after change     : %.6f%n",
                rmse(indirectError, CHANGE_INDEX, SAMPLES));
        System.out.printf("Direct RMSE, after change       : %.6f%n",
                rmse(directError, CHANGE_INDEX, SAMPLES));
        System.out.printf("Indirect mean control energy    : %.6f%n", indirectEnergy);
        System.out.printf("Direct mean control energy      : %.6f%n", directEnergy);
        System.out.printf("Final plant estimate [a,b]      : [%.6f, %.6f]%n",
                plantEstimates[SAMPLES - 1][0], plantEstimates[SAMPLES - 1][1]);
        System.out.printf("Final indirect gains [y,r]      : [%.6f, %.6f]%n",
                indirectGains[SAMPLES - 1][0], indirectGains[SAMPLES - 1][1]);
        System.out.printf("Final direct gains [y,r]        : [%.6f, %.6f]%n",
                directGains[SAMPLES - 1][0], directGains[SAMPLES - 1][1]);

        Path outputPath = Path.of("Chapter13_Lesson4_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(outputPath)) {
            writer.write("k,reference,model_output,indirect_output,direct_output,"
                    + "indirect_control,direct_control,a_hat,b_hat,"
                    + "indirect_gain_y,indirect_gain_r,direct_gain_y,direct_gain_r\n");
            for (int k = 0; k < SAMPLES; k++) {
                writer.write(String.format(Locale.US,
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,"
                                + "%.12f,%.12f,%.12f,%.12f%n",
                        k,
                        reference[k],
                        modelOutput[k + 1],
                        indirectOutput[k + 1],
                        directOutput[k + 1],
                        indirectControl[k],
                        directControl[k],
                        plantEstimates[k][0],
                        plantEstimates[k][1],
                        indirectGains[k][0],
                        indirectGains[k][1],
                        directGains[k][0],
                        directGains[k][1]));
            }
        }
    }
}

Libraries: the script uses base MATLAB. Simulink can realize the same equations using Unit Delay, Sum, Gain, Saturation, and MATLAB Function blocks.

Chapter13_Lesson4.m

% Chapter 13, Lesson 4: Direct versus indirect self-tuning regulators.
% This script uses base MATLAB. Simulink implementation guidance is included
% in the lesson HTML after this code listing.

clear; clc; close all;

samples = 800;
changeIndex = samples / 2;
kAxis = 0:(samples - 1);

reference = -0.2 * ones(1, samples);
reference(kAxis < 120) = 0.0;
reference(kAxis >= 120 & kAxis < 300) = 1.0;
reference(kAxis >= 300 & kAxis < 460) = -0.6;
reference(kAxis >= 460 & kAxis < 620) = 0.8;
reference = reference + 0.08 * sin(0.071 * kAxis) ...
                      + 0.04 * sin(0.19 * kAxis);

plantA = 0.65 * ones(1, samples);
plantB = 0.55 * ones(1, samples);
plantA((changeIndex + 1):end) = 0.78;
plantB((changeIndex + 1):end) = 0.42;
disturbance = 0.012 * sin(0.37 * kAxis) + 0.006 * cos(0.11 * kAxis);

% Stable reference model with unit steady-state gain.
aModel = 0.35;
bModel = 1.0 - aModel;

modelOutput = zeros(1, samples + 1);
indirectOutput = zeros(1, samples + 1);
directOutput = zeros(1, samples + 1);
indirectControl = zeros(1, samples);
directControl = zeros(1, samples);

plantEstimates = zeros(2, samples);
indirectGains = zeros(2, samples);
directGains = zeros(2, samples);

% Indirect STR state: RLS estimate thetaPlant = [aHat; bHat].
thetaPlant = [0.0; 0.5];
covariance = 100.0 * eye(2);
forgettingFactor = 0.995;
minimumInputGain = 0.08;

% Direct STR state: thetaController = [thetaY; thetaR].
thetaController = [0.0; 1.0];
directStep = 1.2;
normalization = 0.2;

for index = 1:samples
    modelOutput(index + 1) = aModel * modelOutput(index) ...
                           + bModel * reference(index);

    % ----------------------- Indirect STR -----------------------
    if index > 1
        phiPlant = [indirectOutput(index - 1); indirectControl(index - 1)];
        denominator = forgettingFactor ...
                    + phiPlant' * covariance * phiPlant;
        rlsGain = covariance * phiPlant / denominator;
        predictionError = indirectOutput(index) - phiPlant' * thetaPlant;
        thetaPlant = thetaPlant + rlsGain * predictionError;
        covariance = (covariance - rlsGain * phiPlant' * covariance) ...
                   / forgettingFactor;
        covariance = 0.5 * (covariance + covariance');
    end

    aHat = thetaPlant(1);
    bHatSafe = max(abs(thetaPlant(2)), minimumInputGain);
    gainY = min(max((aModel - aHat) / bHatSafe, -4.0), 4.0);
    gainR = min(max(bModel / bHatSafe, 0.0), 4.0);

    indirectControl(index) = min(max( ...
        gainY * indirectOutput(index) + gainR * reference(index), ...
        -5.0), 5.0);
    indirectOutput(index + 1) = plantA(index) * indirectOutput(index) ...
                              + plantB(index) * indirectControl(index) ...
                              + disturbance(index);

    plantEstimates(:, index) = thetaPlant;
    indirectGains(:, index) = [gainY; gainR];

    % ------------------------- Direct STR -----------------------
    if index > 1
        errorCurrent = directOutput(index) - modelOutput(index);
        errorPrevious = directOutput(index - 1) - modelOutput(index - 1);
        filteredError = errorCurrent - aModel * errorPrevious;
        phiController = [directOutput(index - 1); reference(index - 1)];
        denominator = normalization + phiController' * phiController;
        thetaController = thetaController ...
            - directStep * phiController * filteredError / denominator;

        % Projection safeguards.
        thetaController(1) = min(max(thetaController(1), -4.0), 4.0);
        thetaController(2) = min(max(thetaController(2), 0.0), 4.0);
    end

    directControl(index) = min(max( ...
        thetaController' * [directOutput(index); reference(index)], ...
        -5.0), 5.0);
    directOutput(index + 1) = plantA(index) * directOutput(index) ...
                            + plantB(index) * directControl(index) ...
                            + disturbance(index);
    directGains(:, index) = thetaController;
end

indirectError = indirectOutput(2:end) - modelOutput(2:end);
directError = directOutput(2:end) - modelOutput(2:end);
rmseValue = @(signal) sqrt(mean(signal .^ 2));

fprintf('Direct versus indirect STR comparison\n');
fprintf('---------------------------------------------\n');
fprintf('Indirect RMSE, all samples      : %.6f\n', rmseValue(indirectError));
fprintf('Direct RMSE, all samples        : %.6f\n', rmseValue(directError));
fprintf('Indirect RMSE, before change    : %.6f\n', ...
    rmseValue(indirectError(1:changeIndex)));
fprintf('Direct RMSE, before change      : %.6f\n', ...
    rmseValue(directError(1:changeIndex)));
fprintf('Indirect RMSE, after change     : %.6f\n', ...
    rmseValue(indirectError((changeIndex + 1):end)));
fprintf('Direct RMSE, after change       : %.6f\n', ...
    rmseValue(directError((changeIndex + 1):end)));
fprintf('Indirect mean control energy    : %.6f\n', mean(indirectControl .^ 2));
fprintf('Direct mean control energy      : %.6f\n', mean(directControl .^ 2));
fprintf('Final plant estimate [a,b]      : [%.6f, %.6f]\n', ...
    plantEstimates(1, end), plantEstimates(2, end));
fprintf('Final indirect gains [y,r]      : [%.6f, %.6f]\n', ...
    indirectGains(1, end), indirectGains(2, end));
fprintf('Final direct gains [y,r]        : [%.6f, %.6f]\n', ...
    directGains(1, end), directGains(2, end));

figure('Name', 'Tracking comparison');
plot(kAxis, reference, '--', 'DisplayName', 'reference'); hold on;
plot(kAxis, modelOutput(2:end), 'DisplayName', 'reference model');
plot(kAxis, indirectOutput(2:end), 'DisplayName', 'indirect STR');
plot(kAxis, directOutput(2:end), 'DisplayName', 'direct STR');
grid on; xlabel('sample k'); ylabel('output');
title('Tracking comparison'); legend('Location', 'best');

figure('Name', 'Controller gains');
plot(kAxis, indirectGains(1, :), 'DisplayName', 'indirect gainY'); hold on;
plot(kAxis, indirectGains(2, :), 'DisplayName', 'indirect gainR');
plot(kAxis, directGains(1, :), '--', 'DisplayName', 'direct thetaY');
plot(kAxis, directGains(2, :), '--', 'DisplayName', 'direct thetaR');
grid on; xlabel('sample k'); ylabel('controller gain');
title('Controller-parameter trajectories'); legend('Location', 'best');

figure('Name', 'Indirect plant estimates');
plot(kAxis, plantEstimates(1, :), 'DisplayName', 'aHat'); hold on;
plot(kAxis, plantEstimates(2, :), 'DisplayName', 'bHat');
xline(changeIndex, '--', 'DisplayName', 'plant change');
grid on; xlabel('sample k'); ylabel('estimated plant parameter');
title('Plant estimates used only by the indirect STR');
legend('Location', 'best');

resultTable = table(kAxis', reference', modelOutput(2:end)', ...
    indirectOutput(2:end)', directOutput(2:end)', indirectControl', ...
    directControl', plantEstimates(1, :)', plantEstimates(2, :)', ...
    indirectGains(1, :)', indirectGains(2, :)', ...
    directGains(1, :)', directGains(2, :)', ...
    'VariableNames', {'k', 'reference', 'modelOutput', 'indirectOutput', ...
    'directOutput', 'indirectControl', 'directControl', 'aHat', 'bHat', ...
    'indirectGainY', 'indirectGainR', 'directGainY', 'directGainR'});
writetable(resultTable, 'Chapter13_Lesson4_results.csv');

15.1 Simulink realization

Construct two parallel closed-loop branches driven by one reference:

1. Implement the reference model with a Unit Delay and the equation \(y_m(k+1)=a_m y_m(k)+b_m r(k)\).
2. Implement each plant with a Unit Delay, two Gain blocks for \(a(k)\) and \(b(k)\), a Sum block, and a disturbance input.
3. In the indirect branch, place the RLS recursion in a MATLAB Function block. Feed its \(\widehat a\) and \(\widehat b\) outputs to a second MATLAB Function block implementing the protected gain map.
4. In the direct branch, use a MATLAB Function block with persistent variables for \(\theta_y\), \(\theta_r\), and the previous tracking error.
5. Apply identical Saturation blocks to both controller outputs and log outputs, controls, estimates, and gains with To Workspace blocks.
6. Use a fixed-step discrete solver and a sample time of one simulation unit so that the block diagram matches the difference equations exactly.

Libraries: the notebook uses built-in Wolfram Language list, matrix, plotting, association, and export functionality.

Chapter13_Lesson4.nb

Notebook[{
 Cell["Chapter 13, Lesson 4: Direct versus Indirect Self-Tuning Regulators", "Title"],
 Cell["The notebook compares certainty-equivalent indirect adaptation with direct controller-gain adaptation on the same time-varying first-order plant.", "Text"],
 Cell["ClearAll[\"Global`*\"];

samples = 800; changeIndex = Floor[samples/2]; kAxis = Range[0, samples - 1];
clip[x_, lower_, upper_] := Min[Max[x, lower], upper]; rmse[x_List] := Sqrt[Mean[x^2]];
referenceAt[k_] := Module[{base}, base = Which[k < 120, 0.0, k < 300, 1.0, k < 460, -0.6, k < 620, 0.8, True, -0.2];
  base + 0.08 Sin[0.071 k] + 0.04 Sin[0.19 k]];
reference = referenceAt /@ kAxis;
plantA = Table[If[k < changeIndex, 0.65, 0.78], {k, kAxis}]; plantB = Table[If[k < changeIndex, 0.55, 0.42], {k, kAxis}];
disturbance = 0.012 Sin[0.37 kAxis] + 0.006 Cos[0.11 kAxis];
aModel = 0.35; bModel = 1.0 - aModel;
modelOutput = ConstantArray[0.0, samples + 1]; indirectOutput = ConstantArray[0.0, samples + 1]; directOutput = ConstantArray[0.0, samples + 1];
indirectControl = ConstantArray[0.0, samples]; directControl = ConstantArray[0.0, samples];
plantEstimates = ConstantArray[{0.0, 0.0}, samples]; indirectGains = ConstantArray[{0.0, 0.0}, samples]; directGains = ConstantArray[{0.0, 0.0}, samples];
thetaPlant = {0.0, 0.5}; covariance = 100.0 IdentityMatrix[2]; forgettingFactor = 0.995; minimumInputGain = 0.08;
thetaController = {0.0, 1.0}; directStep = 1.2; normalization = 0.2;
Do[
  modelOutput[[index + 1]] = aModel modelOutput[[index]] + bModel reference[[index]];
  (* Indirect STR: estimate plant parameters, then synthesize gains. *)
  If[index > 1, phiPlant = {indirectOutput[[index - 1]], indirectControl[[index - 1]]};
    denominator = forgettingFactor + phiPlant.covariance.phiPlant; rlsGain = covariance.phiPlant/denominator;
    predictionError = indirectOutput[[index]] - phiPlant.thetaPlant; thetaPlant = thetaPlant + rlsGain predictionError;
    covariance = (covariance - Outer[Times, rlsGain, phiPlant].covariance)/forgettingFactor;
    covariance = 0.5 (covariance + Transpose[covariance]);];
  aHat = thetaPlant[[1]]; bHatSafe = Max[Abs[thetaPlant[[2]]], minimumInputGain];
  gainY = clip[(aModel - aHat)/bHatSafe, -4.0, 4.0]; gainR = clip[bModel/bHatSafe, 0.0, 4.0];
  indirectControl[[index]] = clip[gainY indirectOutput[[index]] + gainR reference[[index]], -5.0, 5.0];
  indirectOutput[[index + 1]] = plantA[[index]] indirectOutput[[index]] + plantB[[index]] indirectControl[[index]] + disturbance[[index]];
  plantEstimates[[index]] = thetaPlant; indirectGains[[index]] = {gainY, gainR};
  (* Direct STR: update controller gains from filtered tracking error. *)
  If[index > 1, errorCurrent = directOutput[[index]] - modelOutput[[index]];
    errorPrevious = directOutput[[index - 1]] - modelOutput[[index - 1]]; filteredError = errorCurrent - aModel errorPrevious;
    phiController = {directOutput[[index - 1]], reference[[index - 1]]}; denominator = normalization + phiController.phiController;
    thetaController = thetaController - directStep phiController filteredError/denominator;
    thetaController[[1]] = clip[thetaController[[1]], -4.0, 4.0]; thetaController[[2]] = clip[thetaController[[2]], 0.0, 4.0];];
  directControl[[index]] = clip[thetaController.{directOutput[[index]], reference[[index]]}, -5.0, 5.0];
  directOutput[[index + 1]] = plantA[[index]] directOutput[[index]] + plantB[[index]] directControl[[index]] + disturbance[[index]];
  directGains[[index]] = thetaController;, {index, 1, samples}];
indirectError = Rest[indirectOutput] - Rest[modelOutput]; directError = Rest[directOutput] - Rest[modelOutput];
metrics = <|\"IndirectRMSEAll\" -> rmse[indirectError], \"DirectRMSEAll\" -> rmse[directError],
  \"IndirectRMSEBeforeChange\" -> rmse[Take[indirectError, changeIndex]], \"DirectRMSEBeforeChange\" -> rmse[Take[directError, changeIndex]],
  \"IndirectRMSEAfterChange\" -> rmse[Drop[indirectError, changeIndex]], \"DirectRMSEAfterChange\" -> rmse[Drop[directError, changeIndex]],
  \"IndirectControlEnergy\" -> Mean[indirectControl^2], \"DirectControlEnergy\" -> Mean[directControl^2]|>;
Print[Dataset[metrics]]; Print[\"Final plant estimate: \", Last[plantEstimates]]; Print[\"Final indirect gains: \", Last[indirectGains]]; Print[\"Final direct gains: \", Last[directGains]];
trackingPlot = ListLinePlot[{reference, Rest[modelOutput], Rest[indirectOutput], Rest[directOutput]},
  PlotLegends -> {\"reference\", \"reference model\", \"indirect STR\", \"direct STR\"}, Frame -> True,
  FrameLabel -> {\"sample k\", \"output\"}, PlotLabel -> \"Tracking comparison\", ImageSize -> Large];
gainPlot = ListLinePlot[{indirectGains[[All, 1]], indirectGains[[All, 2]], directGains[[All, 1]], directGains[[All, 2]]},
  PlotLegends -> {\"indirect gainY\", \"indirect gainR\", \"direct thetaY\", \"direct thetaR\"}, Frame -> True,
  FrameLabel -> {\"sample k\", \"controller gain\"}, PlotLabel -> \"Controller-parameter trajectories\", ImageSize -> Large];
Print[trackingPlot]; Print[gainPlot];
Export[\"Chapter13_Lesson4_results.csv\", Prepend[MapThread[Join, {List /@ kAxis, List /@ reference,
    List /@ Rest[modelOutput], List /@ Rest[indirectOutput], List /@ Rest[directOutput], List /@ indirectControl,
    List /@ directControl, plantEstimates, indirectGains, directGains}],
  {\"k\", \"reference\", \"modelOutput\", \"indirectOutput\", \"directOutput\", \"indirectControl\", \"directControl\",
   \"aHat\", \"bHat\", \"indirectGainY\", \"indirectGainR\", \"directGainY\", \"directGainR\"}]];", "Input"]
}, WindowSize -> {1100, 800}, StyleDefinitions -> "Default.nb"]        

Problem 1: Derive the ideal gains and tracking-error equation

For \(y(k+1)=ay(k)+bu(k)\), reference model \(y_m(k+1)=a_my_m(k)+b_mr(k)\), and controller \(u(k)=\theta_yy(k)+\theta_rr(k)\), derive the ideal gains and the error equation for arbitrary gain estimates.

Solution. Substituting the controller gives

\[ y(k+1)=(a+b\theta_y)y(k)+b\theta_r r(k). \]

Coefficient matching with the reference model requires

\[ a+b\theta_y^{\star}=a_m,\qquad b\theta_r^{\star}=b_m. \]

Since \(b\neq0\),

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

Writing \(\theta_y=\theta_y^{\star}+\widetilde\theta_y\) and \(\theta_r=\theta_r^{\star}+\widetilde\theta_r\),

\[ \begin{aligned} e(k+1)&=y(k+1)-y_m(k+1)\\ &=a_m e(k)+b\widetilde\theta_y(k)y(k) +b\widetilde\theta_r(k)r(k)\\ &=a_m e(k)+b\widetilde{\boldsymbol\theta}_c^{\mathsf T}(k) \boldsymbol\phi_c(k). \end{aligned} \]

Problem 2: Quantify covariance amplification

Let \(a=0.7\), \(b=0.5\), \(a_m=0.3\), \(b_m=0.7\), and suppose the plant-estimation covariance is

\[ \mathbf{P}_p=\begin{bmatrix}0.0025&0\\0&0.0004\end{bmatrix}. \]

Compute the approximate controller-gain covariance.

Solution. The Jacobian is

\[ \mathbf{J}_{\mathcal F}= \begin{bmatrix} -1/0.5 & -(0.3-0.7)/0.5^2\\ 0 & -0.7/0.5^2 \end{bmatrix} =\begin{bmatrix}-2&1.6\\0&-2.8\end{bmatrix}. \]

Therefore,

\[ \mathbf{P}_c\approx \mathbf{J}_{\mathcal F}\mathbf{P}_p \mathbf{J}_{\mathcal F}^{\mathsf T} =\begin{bmatrix} 0.011024 & -0.001792\\ -0.001792 & 0.003136 \end{bmatrix}. \]

The standard deviations are approximately \(0.1050\) and \(0.0560\). Even a small variance in the estimated input gain creates noticeable controller-gain uncertainty because the map contains \(1/b^2\) terms.

Problem 3: Prove nonincrease of the direct parameter-error norm

In the noise-free case, prove that the normalized direct update is nonexpansive when \(0 < \gamma|b| < 2\).

Solution. The error update is

\[ \widetilde{\boldsymbol\theta}^{+} =\widetilde{\boldsymbol\theta} -\alpha\boldsymbol\phi\boldsymbol\phi^{\mathsf T} \widetilde{\boldsymbol\theta},\qquad \alpha=\frac{\gamma|b|}{\delta+\|\boldsymbol\phi\|^2}. \]

Let \(z=\boldsymbol\phi^{\mathsf T} \widetilde{\boldsymbol\theta}\) and \(s=\|\boldsymbol\phi\|^2\). Then

\[ \begin{aligned} \|\widetilde{\boldsymbol\theta}^{+}\|^2 &=\|\widetilde{\boldsymbol\theta}\|^2 -2\alpha z^2+\alpha^2 s z^2,\\ \|\widetilde{\boldsymbol\theta}^{+}\|^2 -\|\widetilde{\boldsymbol\theta}\|^2 &=-\alpha(2-\alpha s)z^2. \end{aligned} \]

Since \(\delta > 0\), \(\alpha s < \gamma|b| < 2\). Hence \(\alpha(2-\alpha s) > 0\), and the difference is nonpositive. Equality occurs when the current regressor is orthogonal to the parameter error or when the regressor is zero.

Problem 4: Compare gain errors produced by a biased plant estimate

The true plant has \(a=0.72\) and \(b=0.40\). The desired model has \(a_m=0.30\) and \(b_m=0.70\). An indirect estimator temporarily reports \(\widehat a=0.70\) and \(\widehat b=0.32\). Find the true ideal gains, the indirect gains, and the gain errors.

Solution. The ideal gains are

\[ \theta_y^{\star}=\frac{0.30-0.72}{0.40}=-1.05,\qquad \theta_r^{\star}=\frac{0.70}{0.40}=1.75. \]

The certainty-equivalent gains are

\[ \widehat\theta_y=\frac{0.30-0.70}{0.32}=-1.25,\qquad \widehat\theta_r=\frac{0.70}{0.32}=2.1875. \]

Thus

\[ \Delta\theta_y=-0.20,\qquad \Delta\theta_r=0.4375. \]

Although the input-gain estimate is wrong by only \(-0.08\), the reference gain changes substantially because it is inversely proportional to \(\widehat b\).

Problem 5: Select an architecture

A process-control system must provide online estimates of process gain and time constant for operator diagnostics, periodically switch between regulation and set-point tracking designs, and maintain a certified fixed fallback controller. Which STR architecture is more appropriate? State the required safeguards.

Solution. An indirect STR is preferable because the explicit plant estimate serves diagnostics and can be reused by multiple controller-design modes. Required safeguards include:

• projection of plant estimates onto known physical bounds;
• minimum magnitude and fixed sign constraints for the input gain;
• validation of every synthesized controller before activation;
• rate limiting or bumpless transfer between controller coefficients;
• estimator confidence monitoring and automatic fallback control;
• excitation monitoring so that poorly identified models are not trusted.

Direct and indirect STRs can converge to the same ideal controller when their model class, controller structure, and reference model are compatible. Their adaptation pathways are nevertheless different. Indirect STRs estimate plant parameters and transform their uncertainty through a controller-synthesis map; this improves interpretability and redesign flexibility but introduces synthesis singularities and additional transient-stability checks. Direct STRs adapt controller gains through a tracking-error model; this can reduce computation and integrate naturally with a Lyapunov argument, but it requires a valid controller parameterization, known control direction, and careful robustness modifications. Architecture selection should therefore be based on the required information flow and guarantees, not on the labels “direct” or “indirect” alone.