After this lesson, students should be able to:

• separate the identifier state from the physical plant state and controller state;
• write an STR as an estimator followed by a controller-design map;
• derive the certainty-equivalent control law for a simple uncertain plant;
• show explicitly how parameter errors perturb the desired closed-loop model;
• implement a safe sample-by-sample update order with projection and saturation.

An indirect STR does not adapt controller gains as independent free variables. It first estimates a plant parameter vector \( \hat{\boldsymbol{\theta}}(k) \), then computes the controller parameters through a design map \( \mathcal{C} \):

\[ \hat{\boldsymbol{\theta}}(k) \xrightarrow{\text{controller synthesis}} \boldsymbol{\kappa}(k)=\mathcal{C}\!\left(\hat{\boldsymbol{\theta}}(k)\right), \qquad u(k)=\mathcal{U}\!\left(\boldsymbol{\kappa}(k),y(k),r(k)\right). \]

The word self-tuning refers to this repeated model estimation and controller redesign. The architecture is adaptive because the control law changes as measured input-output data change the estimated model.

The STR contains three conceptually distinct dynamical objects:

1. Plant: produces measured data and is governed by unknown parameters.
2. Identifier: updates a model from the measured transition data.
3. Controller and redesign mechanism: maps the current model estimate to implementable gains.

flowchart TD
  R["Reference r(k)"] --> CL["Control law"]
  CL --> U["Actuator command u(k)"]
  U --> P["Unknown plant"]
  P --> Y["Measured output y(k)"]
  Y --> CL
  U --> ID["Online identifier"]
  Y --> ID
  ID --> TH["Estimated model parameters"]
  TH --> DS["Controller-design map"]
  DS --> KG["Updated controller gains"]
  KG --> CL
        

This is a closed-loop learning system. The controller determines the input that generates identification data, while the identifier determines the model used to redesign the controller. Consequently, estimator convergence cannot generally be analyzed as though the input were externally fixed.

Consider the first-order discrete-time SISO plant

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

where \( a \) and \( b \) are constant or slowly varying unknown parameters and \( v(k) \) is a bounded disturbance or measurement-equivalent innovation. Define

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

The one-step predictor used by the identifier is

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

The prediction error, available when \( y(k+1) \) has been measured, is

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

with parameter error \( \tilde{\boldsymbol{\theta}}(k)= \boldsymbol{\theta}-\hat{\boldsymbol{\theta}}(k) \). Thus the identifier receives information about parameter error only through the regressor directions excited by the closed-loop signals.

Detailed recursive least-squares design is postponed to Chapter 12. For the present architectural lesson, a normalized gradient law is sufficient:

\[ \hat{\boldsymbol{\theta}}(k+1)= \operatorname{Proj}_{\Omega}\!\left[ \hat{\boldsymbol{\theta}}(k)+ \gamma\frac{\boldsymbol{\phi}(k)\varepsilon(k+1)} {\delta+\boldsymbol{\phi}^T(k)\boldsymbol{\phi}(k)}\right], \]

where \( \delta > 0 \) prevents division by a small regressor norm, \( 0 < \gamma \le 1 \) is an adaptation gain, and \( \Omega \) is a physically meaningful compact parameter set. For example, if the sign of the high-frequency gain is known to be positive, choose

\[ \Omega = \left\{ (\hat a,\hat b): a_{\min} \le \hat a \le a_{\max}, \; 0 < b_{\min} \le \hat b \le b_{\max} \right\}. \]

4.1 A descent calculation without disturbance

Ignore projection temporarily and set \( v(k)=0 \). Let \( D(k)=\delta+\|\boldsymbol{\phi}(k)\|^2 \) and \( V_\theta(k)=\|\tilde{\boldsymbol{\theta}}(k)\|^2 \). Because

\[ \tilde{\boldsymbol{\theta}}(k+1)= \tilde{\boldsymbol{\theta}}(k)- \gamma\frac{\boldsymbol{\phi}(k)\varepsilon(k+1)}{D(k)}, \]

direct expansion gives

\[ \begin{aligned} \Delta V_\theta(k) &= -2\gamma\frac{\varepsilon^2(k+1)}{D(k)} +\gamma^2\frac{\|\boldsymbol{\phi}(k)\|^2\varepsilon^2(k+1)}{D^2(k)} \\ &\le -\gamma(2-\gamma) \frac{\varepsilon^2(k+1)}{D(k)}\le 0. \end{aligned} \]

The inequality uses \( \|\boldsymbol{\phi}(k)\|^2/D(k)\le 1 \). Projection onto a convex set is nonexpansive relative to any true parameter inside that set, so it does not invalidate the basic boundedness argument. This establishes prediction-error descent, not automatic convergence of every parameter component; persistent excitation remains relevant.

Specify the desired first-order closed-loop relation

\[ y_{d}(k+1)=a_m y(k)+b_m r(k), \qquad |a_m|<1. \]

If \( a \) and \( b \) were known, the algebraic model-matching controller would be

\[ u^*(k)=\frac{a_m-a}{b}y(k)+\frac{b_m}{b}r(k). \]

The indirect STR applies the same design formula to the current estimates:

\[ u(k)=K_y(k)y(k)+K_r(k)r(k), \]

\[ K_y(k)=\frac{a_m-\hat a(k)}{\hat b(k)}, \qquad K_r(k)=\frac{b_m}{\hat b(k)}. \]

This substitution is the certainty-equivalence principle: the estimated model is treated as though it were the true model when the controller is synthesized. Certainty equivalence is an architecture and design principle; by itself it is not a complete proof of closed-loop stability.

Define \( \tilde a(k)=a-\hat a(k) \) and \( \tilde b(k)=b-\hat b(k) \). Substitution of the certainty-equivalent controller into the true plant gives

\[ \begin{aligned} y(k+1) &=a_m y(k)+b_m r(k)+d_\theta(k)+v(k),\\ d_\theta(k) &=\left[\tilde a(k)+ \frac{\tilde b(k)}{\hat b(k)} \left(a_m-\hat a(k)\right)\right]y(k) +\frac{\tilde b(k)}{\hat b(k)}b_m r(k). \end{aligned} \]

Therefore parameter mismatch appears as an endogenous disturbance \( d_\theta(k) \). It is amplified when \( |\hat b(k)| \) becomes small, which explains why a lower bound on the estimated input gain is an architectural necessity, not merely a numerical convenience.

6.1 Conditional tracking result

Let the reference-model state satisfy \( y_m(k+1)=a_m y_m(k)+b_m r(k) \) and define \( e(k)=y(k)-y_m(k) \). Then

\[ e(k+1)=a_m e(k)+d_\theta(k)+v(k). \]

Repeated substitution yields

\[ e(k)=a_m^k e(0)+ \sum_{i=0}^{k-1}a_m^{k-1-i} \left[d_\theta(i)+v(i)\right]. \]

If \( |d_\theta(k)+v(k)|\le\bar d \), then

\[ |e(k)|\le |a_m|^k|e(0)|+ \bar d\frac{1-|a_m|^k}{1-|a_m|}, \qquad \limsup_{k\to\infty}|e(k)|\le \frac{\bar d}{1-|a_m|}. \]

In the disturbance-free case, if the estimates converge to the true parameters and \( \hat b(k) \) stays away from zero, then \( d_\theta(k)\to 0 \) and the stable error recursion gives \( e(k)\to 0 \). This is a conditional result: proving that the coupled adaptive loop actually supplies those conditions is the deeper STR convergence problem.

At sample \( k \), the transition generated by \( y(k-1),u(k-1) \) has just produced \( y(k) \). A causally correct implementation performs the following sequence:

flowchart TD
  A["Measure y(k)"] --> B["Form previous regressor: y(k-1), u(k-1)"]
  B --> C["Compute one-step prediction error"]
  C --> D["Update and project plant estimates"]
  D --> E["Validate estimated model"]
  E --> F["Recompute controller gains"]
  F --> G["Apply gain limits and actuator saturation"]
  G --> H["Compute and apply u(k)"]
  H --> I["Store data for the next sample"]
        

Using \( y(k+1) \) to compute \( u(k) \) would be noncausal. Updating the controller before processing the newest completed transition introduces an avoidable one-sample delay in adaptation. In multirate implementations the identifier may execute every sample while redesign occurs less frequently to reduce computational load and gain chatter.

8.1 Parameter projection

Choose bounds from physical knowledge. For the example, \( \hat b(k)\ge b_{\min}>0 \) guarantees that the controller map is defined and preserves the known control direction.

8.2 Controller validity test

A real redesign module should reject an estimated model if it violates assumptions used by the synthesis method. A generic acceptance rule is

\[ \text{accept}(k)=1 \quad\Longleftrightarrow\quad \hat{\boldsymbol{\theta}}(k)\in\Omega, \;|\hat b(k)|\ge b_{\min}, \;\boldsymbol{\kappa}(k)\in\mathcal{K}_{\text{safe}}. \]

If validation fails, retain the last accepted controller rather than applying a mathematically undefined or excessively aggressive redesign.

8.3 Rate limiting and saturation

Controller gains may be filtered or rate limited:

\[ \boldsymbol{\kappa}_{\text{applied}}(k)= \boldsymbol{\kappa}_{\text{applied}}(k-1)+ \operatorname{sat}_{\Delta\kappa_{\max}} \!\left(\boldsymbol{\kappa}_{\text{new}}(k)- \boldsymbol{\kappa}_{\text{applied}}(k-1)\right). \]

Actuator saturation should also be explicit. Later chapters examine why saturation must be communicated back to the adaptive law; otherwise the identifier may interpret saturation-induced behavior as plant uncertainty.

Good regulation tends to reduce output and input variation. Identification, however, needs informative variation in the regressor. This creates a structural tension:

\[ \text{strong regulation}\quad\Longrightarrow\quad \boldsymbol{\phi}(k)\text{ may lose richness} \quad\Longrightarrow\quad \hat{\boldsymbol{\theta}}(k)\text{ may not converge uniquely}. \]

Closed-loop tracking can therefore be satisfactory while parameter estimates remain biased or nonunique. The controller may need only a combination of parameters sufficient for the desired input-output behavior. This distinction between state/tracking convergence and parameter convergence was established in Chapter 10 and is central to STR analysis.

Deliberate probing can improve identification but increases control cost. The resulting dual-control issue is not solved by the basic certainty-equivalent architecture: the input is chosen for regulation, not explicitly for future information value.

A practical high-level algorithm is:

Initialize parameter estimate, admissible set, safe controller, and data buffers.
For each sample k:
    1. Measure y(k).
    2. Use y(k-1), u(k-1), and y(k) to compute the prediction error.
    3. Update and project the plant-parameter estimate.
    4. Test whether the estimated model is valid for controller synthesis.
    5. If valid, redesign gains; otherwise retain the last accepted gains.
    6. Apply gain-rate limits and calculate u(k) from y(k) and r(k).
    7. Saturate u(k), send it to the plant, and store the transition data.
      

In software, the identifier and controller-design map should be separate modules with explicit interfaces. This prevents accidental algebraic loops, makes fallback behavior testable, and allows the estimator to be replaced later by RLS without changing the control-law interface.

Required libraries: NumPy for numerical arrays and Matplotlib for plots. The estimator is implemented from scratch.

Chapter11_Lesson2.py

# Chapter11_Lesson2.py
# Self-Tuning Regulator architecture: normalized-gradient identifier
# plus certainty-equivalent controller redesign for a first-order plant.

from __future__ import annotations

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


@dataclass
class STRConfig:
    samples: int = 300
    gamma: float = 0.35
    delta: float = 0.10
    a_min: float = -1.50
    a_max: float = 1.50
    b_min: float = 0.10
    b_max: float = 2.00
    u_limit: float = 5.00
    a_model: float = 0.25
    b_model: float = 0.75
    noise_std: float = 0.01
    seed: int = 7


def reference_signal(k: int) -> float:
    if k < 40:
        return 0.0
    if k < 160:
        return 1.0
    return -0.5


def true_a(k: int) -> float:
    # A parameter change demonstrates why retuning is useful.
    return 0.75 if k < 150 else 0.45


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


def run_str(config: STRConfig) -> dict[str, np.ndarray]:
    rng = np.random.default_rng(config.seed)
    n = config.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)
    b_hat_history = np.zeros(n)
    a_true_history = np.zeros(n)

    # Initial model: deliberately inaccurate but with the known positive gain sign.
    a_hat = 0.0
    b_hat = 0.8

    for k in range(n):
        # Identifier update uses the completed transition from k-1 to k.
        if k > 0:
            phi = np.array([y[k - 1], u[k - 1]])
            y_prediction = a_hat * phi[0] + b_hat * phi[1]
            epsilon = y[k] - y_prediction
            denominator = config.delta + float(phi @ phi)

            theta_update = config.gamma * phi * epsilon / denominator
            a_hat = project(a_hat + theta_update[0], config.a_min, config.a_max)
            b_hat = project(b_hat + theta_update[1], config.b_min, config.b_max)
            prediction_error[k] = epsilon

        # Certainty-equivalent redesign from the current estimates.
        r[k] = reference_signal(k)
        b_design = max(b_hat, config.b_min)
        k_y = (config.a_model - a_hat) / b_design
        k_r = config.b_model / b_design

        u_unsaturated = k_y * y[k] + k_r * r[k]
        u[k] = float(np.clip(u_unsaturated, -config.u_limit, config.u_limit))

        a_plant = true_a(k)
        process_noise = rng.normal(0.0, config.noise_std)
        y[k + 1] = a_plant * y[k] + 1.2 * u[k] + process_noise

        a_hat_history[k] = a_hat
        b_hat_history[k] = b_hat
        a_true_history[k] = a_plant

    return {
        "y": y,
        "u": u,
        "r": r,
        "prediction_error": prediction_error,
        "a_hat": a_hat_history,
        "b_hat": b_hat_history,
        "a_true": a_true_history,
    }


def main() -> None:
    config = STRConfig()
    data = run_str(config)
    k = np.arange(config.samples)

    print(f"Final a_hat = {data['a_hat'][-1]:.6f}")
    print(f"Final b_hat = {data['b_hat'][-1]:.6f}")
    print(f"RMS prediction error = {np.sqrt(np.mean(data['prediction_error'][1:] ** 2)):.6f}")

    plt.figure(figsize=(10, 5))
    plt.step(k, data["r"], where="post", label="reference r(k)")
    plt.plot(np.arange(config.samples + 1), data["y"], label="output y(k)")
    plt.xlabel("sample k")
    plt.ylabel("signal")
    plt.title("Indirect STR: reference tracking")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(10, 5))
    plt.plot(k, data["a_true"], label="true a(k)")
    plt.plot(k, data["a_hat"], label="estimated a_hat(k)")
    plt.plot(k, data["b_hat"], label="estimated b_hat(k)")
    plt.xlabel("sample k")
    plt.ylabel("parameter value")
    plt.title("Identifier state used by the controller redesign")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()
      

This version uses only the C++17 standard library. The controller and normalized-gradient identifier are implemented directly.

Chapter11_Lesson2.cpp

// Chapter11_Lesson2.cpp
// Build: g++ -std=c++17 -O2 Chapter11_Lesson2.cpp -o str_demo

#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>

struct STRConfig {
    int samples = 300;
    double gamma = 0.35;
    double delta = 0.10;
    double aMin = -1.50;
    double aMax = 1.50;
    double bMin = 0.10;
    double bMax = 2.00;
    double uLimit = 5.00;
    double aModel = 0.25;
    double bModel = 0.75;
    double noiseStd = 0.01;
};

double referenceSignal(int k) {
    if (k < 40) return 0.0;
    if (k < 160) return 1.0;
    return -0.5;
}

double trueA(int k) {
    return (k < 150) ? 0.75 : 0.45;
}

double project(double value, double lower, double upper) {
    return std::clamp(value, lower, upper);
}

int main() {
    const STRConfig cfg;
    std::mt19937 generator(7);
    std::normal_distribution<double> noise(0.0, cfg.noiseStd);

    std::vector<double> y(cfg.samples + 1, 0.0);
    std::vector<double> u(cfg.samples, 0.0);
    std::vector<double> predictionError(cfg.samples, 0.0);
    std::vector<double> aHatHistory(cfg.samples, 0.0);
    std::vector<double> bHatHistory(cfg.samples, 0.0);

    double aHat = 0.0;
    double bHat = 0.8;

    for (int k = 0; k < cfg.samples; ++k) {
        if (k > 0) {
            const std::array<double, 2> phi{y[k - 1], u[k - 1]};
            const double yPrediction = aHat * phi[0] + bHat * phi[1];
            const double epsilon = y[k] - yPrediction;
            const double denominator = cfg.delta + phi[0] * phi[0] + phi[1] * phi[1];

            aHat = project(aHat + cfg.gamma * phi[0] * epsilon / denominator,
                           cfg.aMin, cfg.aMax);
            bHat = project(bHat + cfg.gamma * phi[1] * epsilon / denominator,
                           cfg.bMin, cfg.bMax);
            predictionError[k] = epsilon;
        }

        const double reference = referenceSignal(k);
        const double bDesign = std::max(bHat, cfg.bMin);
        const double kY = (cfg.aModel - aHat) / bDesign;
        const double kR = cfg.bModel / bDesign;
        const double uUnbounded = kY * y[k] + kR * reference;
        u[k] = std::clamp(uUnbounded, -cfg.uLimit, cfg.uLimit);

        y[k + 1] = trueA(k) * y[k] + 1.2 * u[k] + noise(generator);
        aHatHistory[k] = aHat;
        bHatHistory[k] = bHat;
    }

    double squaredErrorSum = 0.0;
    for (int k = 1; k < cfg.samples; ++k) {
        squaredErrorSum += predictionError[k] * predictionError[k];
    }
    const double rmsError = std::sqrt(squaredErrorSum / (cfg.samples - 1));

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Final a_hat = " << aHatHistory.back() << '\n';
    std::cout << "Final b_hat = " << bHatHistory.back() << '\n';
    std::cout << "RMS prediction error = " << rmsError << '\n';
    std::cout << "Final output y(N) = " << y.back() << '\n';
    return 0;
}
      

This implementation uses the Java standard library and preserves the same update ordering and safeguards as the Python and C++ versions.

Chapter11_Lesson2.java

// Chapter11_Lesson2.java
// Build and run: javac Chapter11_Lesson2.java && java Chapter11_Lesson2

import java.util.Random;

public final class Chapter11_Lesson2 {
    private static final int SAMPLES = 300;
    private static final double GAMMA = 0.35;
    private static final double DELTA = 0.10;
    private static final double A_MIN = -1.50;
    private static final double A_MAX = 1.50;
    private static final double B_MIN = 0.10;
    private static final double B_MAX = 2.00;
    private static final double U_LIMIT = 5.00;
    private static final double A_MODEL = 0.25;
    private static final double B_MODEL = 0.75;
    private static final double NOISE_STD = 0.01;

    private Chapter11_Lesson2() {
    }

    private static double referenceSignal(int k) {
        if (k < 40) return 0.0;
        if (k < 160) return 1.0;
        return -0.5;
    }

    private static double trueA(int k) {
        return (k < 150) ? 0.75 : 0.45;
    }

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

    public static void main(String[] args) {
        Random random = new Random(7L);
        double[] y = new double[SAMPLES + 1];
        double[] u = new double[SAMPLES];
        double[] predictionError = new double[SAMPLES];
        double[] aHatHistory = new double[SAMPLES];
        double[] bHatHistory = new double[SAMPLES];

        double aHat = 0.0;
        double bHat = 0.8;

        for (int k = 0; k < SAMPLES; k++) {
            if (k > 0) {
                double phiY = y[k - 1];
                double phiU = u[k - 1];
                double yPrediction = aHat * phiY + bHat * phiU;
                double epsilon = y[k] - yPrediction;
                double denominator = DELTA + phiY * phiY + phiU * phiU;

                aHat = clamp(aHat + GAMMA * phiY * epsilon / denominator, A_MIN, A_MAX);
                bHat = clamp(bHat + GAMMA * phiU * epsilon / denominator, B_MIN, B_MAX);
                predictionError[k] = epsilon;
            }

            double reference = referenceSignal(k);
            double bDesign = Math.max(bHat, B_MIN);
            double kY = (A_MODEL - aHat) / bDesign;
            double kR = B_MODEL / bDesign;
            double uUnbounded = kY * y[k] + kR * reference;
            u[k] = clamp(uUnbounded, -U_LIMIT, U_LIMIT);

            double processNoise = NOISE_STD * random.nextGaussian();
            y[k + 1] = trueA(k) * y[k] + 1.2 * u[k] + processNoise;
            aHatHistory[k] = aHat;
            bHatHistory[k] = bHat;
        }

        double squaredErrorSum = 0.0;
        for (int k = 1; k < SAMPLES; k++) {
            squaredErrorSum += predictionError[k] * predictionError[k];
        }
        double rmsError = Math.sqrt(squaredErrorSum / (SAMPLES - 1));

        System.out.printf("Final a_hat = %.6f%n", aHatHistory[SAMPLES - 1]);
        System.out.printf("Final b_hat = %.6f%n", bHatHistory[SAMPLES - 1]);
        System.out.printf("RMS prediction error = %.6f%n", rmsError);
        System.out.printf("Final output y(N) = %.6f%n", y[SAMPLES]);
    }
}
      

The MATLAB script simulates the STR directly. Its final local function can create a high-level Simulink architecture model when Simulink is installed. The estimator equations remain visible rather than hidden in a toolbox block.

Chapter11_Lesson2.m

% Chapter11_Lesson2.m
% Indirect STR with a normalized-gradient identifier and certainty-equivalent
% pole-placement controller. The last function optionally creates a Simulink
% architecture diagram when Simulink is installed.

clear; clc; close all;
rng(7);

N = 300;
gamma = 0.35;
delta = 0.10;
aBounds = [-1.50, 1.50];
bBounds = [0.10, 2.00];
uLimit = 5.00;
aModel = 0.25;
bModel = 0.75;
noiseStd = 0.01;

y = zeros(N + 1, 1);
u = zeros(N, 1);
r = zeros(N, 1);
predictionError = zeros(N, 1);
aHatHistory = zeros(N, 1);
bHatHistory = zeros(N, 1);
aTrueHistory = zeros(N, 1);

aHat = 0.0;
bHat = 0.8;

for k = 1:N
    if k > 1
        phi = [y(k - 1); u(k - 1)];
        yPrediction = [aHat, bHat] * phi;
        epsilon = y(k) - yPrediction;
        denominator = delta + phi.' * phi;
        thetaUpdate = gamma * phi * epsilon / denominator;

        aHat = min(max(aHat + thetaUpdate(1), aBounds(1)), aBounds(2));
        bHat = min(max(bHat + thetaUpdate(2), bBounds(1)), bBounds(2));
        predictionError(k) = epsilon;
    end

    if k <= 40
        r(k) = 0.0;
    elseif k <= 160
        r(k) = 1.0;
    else
        r(k) = -0.5;
    end

    bDesign = max(bHat, bBounds(1));
    kY = (aModel - aHat) / bDesign;
    kR = bModel / bDesign;
    uUnbounded = kY * y(k) + kR * r(k);
    u(k) = min(max(uUnbounded, -uLimit), uLimit);

    if k <= 150
        aPlant = 0.75;
    else
        aPlant = 0.45;
    end
    y(k + 1) = aPlant * y(k) + 1.2 * u(k) + noiseStd * randn;

    aHatHistory(k) = aHat;
    bHatHistory(k) = bHat;
    aTrueHistory(k) = aPlant;
end

fprintf('Final a_hat = %.6f\n', aHatHistory(end));
fprintf('Final b_hat = %.6f\n', bHatHistory(end));
fprintf('RMS prediction error = %.6f\n', ...
    sqrt(mean(predictionError(2:end) .^ 2)));

figure;
stairs(0:N-1, r, 'DisplayName', 'reference r(k)'); hold on;
plot(0:N, y, 'DisplayName', 'output y(k)');
grid on; xlabel('sample k'); ylabel('signal');
title('Indirect STR: reference tracking'); legend('Location', 'best');

figure;
plot(0:N-1, aTrueHistory, 'DisplayName', 'true a(k)'); hold on;
plot(0:N-1, aHatHistory, 'DisplayName', 'estimated aHat(k)');
plot(0:N-1, bHatHistory, 'DisplayName', 'estimated bHat(k)');
grid on; xlabel('sample k'); ylabel('parameter value');
title('Identifier state used by controller redesign'); legend('Location', 'best');

% Uncomment to create a high-level Simulink architecture model.
% buildSTRSimulinkArchitecture('Chapter11_Lesson2_STR_Architecture');

function buildSTRSimulinkArchitecture(modelName)
    if ~license('test', 'Simulink')
        error('Simulink is not available in this MATLAB installation.');
    end

    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/In1', [modelName '/Reference'], ...
        'Position', [30 70 60 90]);
    add_block('simulink/Ports & Subsystems/Subsystem', [modelName '/Identifier'], ...
        'Position', [250 175 380 245]);
    add_block('simulink/Ports & Subsystems/Subsystem', [modelName '/Controller Redesign'], ...
        'Position', [460 175 620 245]);
    add_block('simulink/Ports & Subsystems/Subsystem', [modelName '/Control Law'], ...
        'Position', [460 45 600 125]);
    add_block('simulink/Ports & Subsystems/Subsystem', [modelName '/Plant'], ...
        'Position', [720 50 830 120]);
    add_block('simulink/Sinks/Out1', [modelName '/Output'], ...
        'Position', [930 70 960 90]);

    % Add the extra subsystem input ports required by the architecture.
    add_block('simulink/Ports & Subsystems/In1', ...
        [modelName '/Identifier/Input u'], 'Port', '2', ...
        'Position', [40 75 70 95]);
    add_block('simulink/Ports & Subsystems/In1', ...
        [modelName '/Control Law/Input y'], 'Port', '2', ...
        'Position', [40 75 70 95]);
    add_block('simulink/Ports & Subsystems/In1', ...
        [modelName '/Control Law/Input gains'], 'Port', '3', ...
        'Position', [40 125 70 145]);

    add_line(modelName, 'Reference/1', 'Control Law/1', 'autorouting', 'on');
    add_line(modelName, 'Control Law/1', 'Plant/1', 'autorouting', 'on');
    add_line(modelName, 'Control Law/1', 'Identifier/2', 'autorouting', 'on');
    add_line(modelName, 'Plant/1', 'Output/1', 'autorouting', 'on');
    add_line(modelName, 'Plant/1', 'Identifier/1', 'autorouting', 'on');
    add_line(modelName, 'Plant/1', 'Control Law/2', 'autorouting', 'on');
    add_line(modelName, 'Identifier/1', 'Controller Redesign/1', 'autorouting', 'on');
    add_line(modelName, 'Controller Redesign/1', 'Control Law/3', 'autorouting', 'on');

    set_param(modelName, 'StopTime', '300');
    save_system(modelName);
end
      

The following is a textual Wolfram Notebook expression. Saving it as the indicated .nb file produces an openable notebook containing the simulation and plots.

Chapter11_Lesson2.nb

Notebook[{
 Cell["Chapter 11, Lesson 2: Self-Tuning Regulator Architecture", "Title"],
 Cell["Normalized-gradient identifier plus certainty-equivalent controller redesign for a first-order discrete-time plant.", "Text"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"n", "=", "300"}], ";",
   RowBox[{"gamma", "=", "0.35"}], ";",
   RowBox[{"delta", "=", "0.10"}], ";",
   RowBox[{"aModel", "=", "0.25"}], ";",
   RowBox[{"bModel", "=", "0.75"}], ";",
   RowBox[{"uLimit", "=", "5.0"}], ";",
   RowBox[{"SeedRandom", "[", "7", "]"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"reference", "[", "k_", "]", ":=", "Which", "[",
    RowBox[{"k", "<", "40", ",", "0.0", ",", "k", "<", "160", ",", "1.0", ",", "True", ",", RowBox[{"-", "0.5"}]}], "]"}], ";",
   RowBox[{"trueA", "[", "k_", "]", ":=", "If", "[", RowBox[{"k", "<", "150", ",", "0.75", ",", "0.45"}], "]"}], ";",
   RowBox[{"clip", "[", RowBox[{"x_", ",", "lo_", ",", "hi_"}], "]", ":=", "Min", "[", RowBox[{"hi", ",", RowBox[{"Max", "[", RowBox[{"lo", ",", "x"}], "]"}]}], "]"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"y", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", RowBox[{"n", "+", "1"}]}], "]"}]}], ";",
   RowBox[{"u", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"eps", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"aHistory", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"bHistory", "=", RowBox[{"ConstantArray", "[", RowBox[{"0.0", ",", "n"}], "]"}]}], ";",
   RowBox[{"aHat", "=", "0.0"}], ";",
   RowBox[{"bHat", "=", "0.8"}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Do", "[",
   RowBox[{
    RowBox[{"If", "[", RowBox[{"k", ">", "1", ",",
     RowBox[{
      RowBox[{"phi", "=", RowBox[{"{", RowBox[{RowBox[{"y", "[[", RowBox[{"k", "-", "1"}], "]]"}], ",", RowBox[{"u", "[[", RowBox[{"k", "-", "1"}], "]]"}]}], "}"}]}], ";",
      RowBox[{"prediction", "=", RowBox[{RowBox[{"{", RowBox[{"aHat", ",", "bHat"}], "}"}], ".", "phi"}]}], ";",
      RowBox[{RowBox[{"eps", "[[", "k", "]]"}], "=", RowBox[{RowBox[{"y", "[[", "k", "]]"}], "-", "prediction"}]}], ";",
      RowBox[{"update", "=", RowBox[{"gamma", "*", "phi", "*", RowBox[{"eps", "[[", "k", "]]"}], "/", RowBox[{"(", RowBox[{"delta", "+", RowBox[{"phi", ".", "phi"}]}], ")"}]}]}], ";",
      RowBox[{"aHat", "=", RowBox[{"clip", "[", RowBox[{RowBox[{"aHat", "+", RowBox[{"update", "[[", "1", "]]"}]}], ",", RowBox[{"-", "1.5"}], ",", "1.5"}], "]"}]}], ";",
      RowBox[{"bHat", "=", RowBox[{"clip", "[", RowBox[{RowBox[{"bHat", "+", RowBox[{"update", "[[", "2", "]]"}]}], ",", "0.1", ",", "2.0"}], "]"}]}], ";"
     }]}], "]"}], ";",
    RowBox[{"bDesign", "=", RowBox[{"Max", "[", RowBox[{"bHat", ",", "0.1"}], "]"}]}], ";",
    RowBox[{"kY", "=", RowBox[{RowBox[{"(", RowBox[{"aModel", "-", "aHat"}], ")"}], "/", "bDesign"}]}], ";",
    RowBox[{"kR", "=", RowBox[{"bModel", "/", "bDesign"}]}], ";",
    RowBox[{RowBox[{"u", "[[", "k", "]]"}], "=", RowBox[{"clip", "[", RowBox[{RowBox[{RowBox[{"kY", "*", RowBox[{"y", "[[", "k", "]]"}]}], "+", RowBox[{"kR", "*", RowBox[{"reference", "[", RowBox[{"k", "-", "1"}], "]"}]}]}], ",", RowBox[{"-", "uLimit"}], ",", "uLimit"}], "]"}]}], ";",
    RowBox[{RowBox[{"y", "[[", RowBox[{"k", "+", "1"}], "]]"}], "=", RowBox[{RowBox[{RowBox[{"trueA", "[", RowBox[{"k", "-", "1"}], "]"}], "*", RowBox[{"y", "[[", "k", "]]"}]}], "+", RowBox[{"1.2", "*", RowBox[{"u", "[[", "k", "]]"}]}], "+", RowBox[{"RandomVariate", "[", RowBox[{"NormalDistribution", "[", RowBox[{"0", ",", "0.01"}], "]"}], "]"}]}]}], ";",
    RowBox[{RowBox[{"aHistory", "[[", "k", "]]"}], "=", "aHat"}], ";",
    RowBox[{RowBox[{"bHistory", "[[", "k", "]]"}], "=", "bHat"}], ";"
   }], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "n"}], "}"}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Print", "[", RowBox[{"\"Final aHat = \"", ",", RowBox[{"Last", "[", "aHistory", "]"}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Print", "[", RowBox[{"\"Final bHat = \"", ",", RowBox[{"Last", "[", "bHistory", "]"}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"ListLinePlot", "[",
   RowBox[{RowBox[{"{", RowBox[{"y", ",", RowBox[{"PadRight", "[", RowBox[{RowBox[{"Table", "[", RowBox[{RowBox[{"reference", "[", "k", "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "0", ",", RowBox[{"n", "-", "1"}]}], "}"}]}], "]"}], ",", RowBox[{"n", "+", "1"}], ",", RowBox[{"reference", "[", RowBox[{"n", "-", "1"}], "]"}]}], "]"}]}], "}"}], ",", RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"y(k)\"", ",", "\"r(k)\""}], "}"}]}], ",", RowBox[{"PlotRange", "->", "All"}]}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"ListLinePlot", "[",
   RowBox[{RowBox[{"{", RowBox[{"aHistory", ",", "bHistory", ",", RowBox[{"Table", "[", RowBox[{RowBox[{"trueA", "[", "k", "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "0", ",", RowBox[{"n", "-", "1"}]}], "}"}]}], "]"}]}], "}"}], ",", RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"aHat(k)\"", ",", "\"bHat(k)\"", ",", "\"true a(k)\""}], "}"}]}], ",", RowBox[{"PlotRange", "->", "All"}]}], "]"}]], "Input"]
}, WindowTitle -> "Chapter11_Lesson2", StyleDefinitions -> "Default.nb"]
      

Problem 1 (Certainty-Equivalent Gain Derivation): For \( y(k+1)=a y(k)+b u(k) \), derive a feedback law \( u(k)=K_y y(k)+K_r r(k) \) that produces \( y(k+1)=a_m y(k)+b_m r(k) \) when the parameters are known.

Solution: Substitution gives

\[ y(k+1)=\left(a+bK_y\right)y(k)+bK_r r(k). \]

Coefficient matching requires

\[ a+bK_y=a_m, \qquad bK_r=b_m. \]

Because \( b\neq0 \),

\[ K_y=\frac{a_m-a}{b}, \qquad K_r=\frac{b_m}{b}. \]

The STR replaces \( a,b \) with \( \hat a(k),\hat b(k) \).

Problem 2 (Why the Input-Gain Bound Is Necessary): Suppose \( \hat b(k) \) approaches zero while \( a_m-\hat a(k) \) and \( b_m \) remain nonzero. Explain the consequence for the control law.

Solution: The gains satisfy

\[ |K_y(k)|=\frac{|a_m-\hat a(k)|}{|\hat b(k)|}, \qquad |K_r(k)|=\frac{|b_m|}{|\hat b(k)|}. \]

Both can become arbitrarily large. The commanded input can therefore saturate or destabilize unmodeled dynamics. Projection enforcing \( |\hat b(k)|\ge b_{\min}>0 \), together with gain and actuator limits, keeps the redesign map well posed.

Problem 3 (Prediction-Error Descent): In the disturbance-free case, prove that the unprojected normalized gradient update has nonincreasing \( V_\theta(k)=\|\tilde{\boldsymbol{\theta}}(k)\|^2 \) for \( 0<\gamma\le1 \).

Solution: With \( \varepsilon=\tilde{\boldsymbol{\theta}}^T\boldsymbol{\phi} \) and \( D=\delta+\|\boldsymbol{\phi}\|^2 \), expansion yields

\[ \Delta V_\theta=-2\gamma\frac{\varepsilon^2}{D}+ \gamma^2\frac{\|\boldsymbol{\phi}\|^2\varepsilon^2}{D^2}. \]

Since \( \|\boldsymbol{\phi}\|^2/D\le1 \),

\[ \Delta V_\theta\le- \gamma(2-\gamma)\frac{\varepsilon^2}{D}\le0. \]

Hence the parameter-error norm is bounded and the normalized prediction error is square summable. Unique parameter convergence still requires sufficient excitation.

Problem 4 (Tracking Bound): Given \( e(k+1)=a_m e(k)+w(k) \), with \( |a_m|<1 \) and \( |w(k)|\le\bar w \), derive a bound on \( e(k) \).

Solution: Iterating the recursion gives

\[ e(k)=a_m^k e(0)+\sum_{i=0}^{k-1}a_m^{k-1-i}w(i). \]

Taking absolute values and summing the geometric series gives

\[ |e(k)|\le|a_m|^k|e(0)|+ \bar w\frac{1-|a_m|^k}{1-|a_m|}. \]

Thus the error is ultimately bounded by \( \bar w/(1-|a_m|) \). In the STR, \( w(k)=d_\theta(k)+v(k) \).

Problem 5 (Causal Update Ordering): At sample \( k \), identify the newest transition that can be used by the identifier and state why the controller cannot use \( y(k+1) \) to compute \( u(k) \).

Solution: The newest completed transition is \( (y(k-1),u(k-1))\mapsto y(k) \). The value \( y(k+1) \) is produced only after \( u(k) \) has been applied, so using it in \( u(k) \) would require future information and would be noncausal.

An indirect self-tuning regulator is a coupled identifier–redesign–control loop. The identifier converts completed input-output transitions into a plant estimate; the controller-design map converts that estimate into gains; and the certainty-equivalent controller closes the physical loop. For the first-order example, parameter errors enter the desired closed-loop dynamics as a calculable disturbance. Projection, model validation, gain-rate limits, and actuator saturation keep the redesign mathematically and physically well posed. The next lesson examines how online estimates drive controller updates in greater detail, while Chapter 12 develops RLS as the principal estimator.