After completing this lesson, students should be able to:

• form a control-oriented regression model from sampled plant data;
• derive a normalized online parameter update from prediction error;
• map estimated plant parameters into updated feedback and feedforward gains;
• distinguish tracking convergence from full parameter convergence;
• explain why identifier and controller dynamics cannot be analyzed independently;
• implement projection, denominator protection, saturation, and excitation safeguards.

The lesson assumes the concepts developed earlier in the course: linear parameterization, gradient adaptation, projection, normalization, persistent excitation, and the identifier-controller structure of a self-tuning regulator. Detailed recursive least-squares derivations are intentionally postponed to Chapter 12.

Consider a sampled first-order SISO plant whose parameters are constant or slowly varying over the estimator time scale:

\[ y(k+1)=\alpha(k)y(k)+\beta(k)u(k)+d(k), \]

where \( y(k) \) is the measured output, \( u(k) \) is the applied input, and \( d(k) \) collects measurement noise, process disturbance, and unmodeled dynamics. Introduce the regression form

\[ y(k+1)=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}(k)+d(k), \quad \boldsymbol{\phi}(k)=\begin{bmatrix}y(k)\\u(k)\end{bmatrix}, \quad \boldsymbol{\theta}(k)=\begin{bmatrix}\alpha(k)\\\beta(k)\end{bmatrix}. \]

The control direction is assumed known, so the admissible parameter set is

\[ \Omega=\left\{(\alpha,\beta):\alpha_{\min}\le\alpha\le\alpha_{\max},\; 0<\beta_{\min}\le\beta\le\beta_{\max}\right\}. \]

The strict lower bound on \( \beta \) is not merely a numerical convenience. The controller divides by its estimate, so allowing \( \hat{\beta}(k) \) to approach zero can produce arbitrarily large gains.

flowchart TD
  R["Reference r(k)"] --> C["Controller using current estimates"]
  C --> U["Applied input u(k)"]
  U --> P["Unknown plant"]
  P --> Y["Measured output y(k+1)"]
  Y --> E["Online estimator"]
  U --> E
  E --> T["Estimated alpha and beta"]
  T --> D["Controller gain calculation"]
  D --> C
        

Before observing \( y(k+1) \), the current model predicts

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

The one-step prediction error is

\[ \varepsilon(k+1)=y(k+1)-\hat{y}(k+1\mid k). \]

Using the instantaneous quadratic loss \( J(k+1)=\tfrac{1}{2}\varepsilon^{2}(k+1) \), its gradient with respect to the parameter estimate is

\[ \frac{\partial J(k+1)}{\partial\hat{\boldsymbol{\theta}}(k)} =-\boldsymbol{\phi}(k)\varepsilon(k+1). \]

A normalized gradient step followed by projection onto the admissible set is

\[ \hat{\boldsymbol{\theta}}(k+1)=\Pi_{\Omega}\!\left[ \hat{\boldsymbol{\theta}}(k)+ \frac{\gamma\boldsymbol{\phi}(k)\varepsilon(k+1)} {\delta+\boldsymbol{\phi}^{T}(k)\boldsymbol{\phi}(k)}\right], \quad \delta>0. \]

Normalization prevents a large regressor from creating a disproportionately large parameter jump. Projection enforces prior physical information and keeps the controller calculation well posed. The scalar component updates before projection are

\[ \begin{aligned} \hat{\alpha}^{+} &= \hat{\alpha}(k)+ \frac{\gamma y(k)\varepsilon(k+1)}{\delta+y^{2}(k)+u^{2}(k)},\\ \hat{\beta}^{+} &= \hat{\beta}(k)+ \frac{\gamma u(k)\varepsilon(k+1)}{\delta+y^{2}(k)+u^{2}(k)}. \end{aligned} \]

Here, superscript \( + \) denotes the unconstrained candidate. The projected values become the estimates used at the next sample.

Choose a desired stable one-step response with pole \( p_d \):

\[ y_d(k+1)=p_d y(k)+(1-p_d)r(k), \quad |p_d|<1. \]

Certainty equivalence treats the current estimates as though they were the true parameters when solving the nominal control problem. Setting the estimated next output equal to the desired next output gives

\[ \hat{\alpha}(k)y(k)+\hat{\beta}(k)u_c(k)=y_d(k+1), \]

and therefore

\[ u_c(k)=\frac{y_d(k+1)-\hat{\alpha}(k)y(k)}{\hat{\beta}(k)}. \]

The same controller can be written in gain form:

\[ u_c(k)=-\hat{k}_y(k)y(k)+\hat{k}_r(k)r(k), \]

\[ \hat{k}_y(k)=\frac{\hat{\alpha}(k)-p_d}{\hat{\beta}(k)}, \qquad \hat{k}_r(k)=\frac{1-p_d}{\hat{\beta}(k)}. \]

Thus, online parameter estimation drives controller update through the nonlinear mapping \( (\hat{\alpha},\hat{\beta})\mapsto (\hat{k}_y,\hat{k}_r) \). The applied input may include bounded excitation and saturation:

\[ u(k)=\operatorname{sat}_{[-u_{\max},u_{\max}]} \left(u_c(k)+q(k)\right), \]

where \( q(k) \) is a small probing signal. The probe is not required for the algebraic controller design, but it may be required to keep the regression informative.

First ignore saturation, probing, and disturbance. If the estimates equal the true parameters, substitution of the certainty-equivalence law into the plant gives

\[ \begin{aligned} y(k+1) &=\alpha y(k)+\beta \frac{p_d y(k)+(1-p_d)r(k)-\alpha y(k)}{\beta}\\ &=p_d y(k)+(1-p_d)r(k). \end{aligned} \]

For a constant reference, define \( e(k)=y(k)-r \). Then

\[ e(k+1)=p_d e(k), \qquad e(k)=p_d^{k}e(0), \]

so the nominal tracking error converges geometrically when \( |p_d|<1 \).

Now define the parameter error as \( \tilde{\boldsymbol{\theta}}(k)= \boldsymbol{\theta}(k)-\hat{\boldsymbol{\theta}}(k) \). The actual deviation from the desired one-step response is

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

This equation is fundamental: controller performance depends on parameter error only through the regressor-weighted combination \( \boldsymbol{\phi}^{T}\tilde{\boldsymbol{\theta}} \). Consequently, excellent tracking can coexist with incorrect individual parameter estimates when the closed-loop data do not excite every parameter direction.

Consider the noise-free, constant-parameter case and temporarily omit projection. Then

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

Let

\[ V(k)=\frac{1}{2}\tilde{\boldsymbol{\theta}}^{T}(k) \tilde{\boldsymbol{\theta}}(k). \]

Because the true parameter is fixed, the normalized update implies

\[ \tilde{\boldsymbol{\theta}}(k+1)= \tilde{\boldsymbol{\theta}}(k)- \frac{\gamma\boldsymbol{\phi}(k)\varepsilon(k+1)} {\delta+\|\boldsymbol{\phi}(k)\|^{2}}. \]

Expanding the squared norm yields

\[ \begin{aligned} V(k+1)-V(k) &=-\frac{\gamma\varepsilon^{2}(k+1)} {\delta+\|\boldsymbol{\phi}(k)\|^{2}}\\ &\quad+\frac{\gamma^{2}\|\boldsymbol{\phi}(k)\|^{2} \varepsilon^{2}(k+1)} {2\left(\delta+\|\boldsymbol{\phi}(k)\|^{2}\right)^{2}}\\ &=-\frac{\gamma\varepsilon^{2}(k+1)} {\delta+\|\boldsymbol{\phi}(k)\|^{2}} \left(1-\frac{\gamma\|\boldsymbol{\phi}(k)\|^{2}} {2\left(\delta+\|\boldsymbol{\phi}(k)\|^{2}\right)}\right). \end{aligned} \]

Since \( 0\le\|\boldsymbol{\phi}\|^{2}/ (\delta+\|\boldsymbol{\phi}\|^{2})<1 \), choosing \( 0<\gamma<2 \) gives

\[ V(k+1)-V(k)\le 0. \]

Therefore the parameter-error norm is nonincreasing in this idealized case, and the normalized prediction error is square summable. This proof alone does not establish \( \tilde{\boldsymbol{\theta}}(k)→\mathbf{0} \). Persistent excitation is still required for convergence of every parameter component. Projection preserves boundedness and, for a convex admissible set containing the true parameter, does not increase the distance to that set.

In offline identification, the experimenter may choose an input primarily to reveal the plant dynamics. In an STR, the controller chooses the input primarily to regulate the plant. Hence

\[ \hat{\boldsymbol{\theta}}(k) →u(k)→\boldsymbol{\phi}(k) →\varepsilon(k+1)→ \hat{\boldsymbol{\theta}}(k+1). \]

This loop invalidates a naive separation argument. Better tracking often reduces signal variation, which can weaken excitation and slow identification. Conversely, aggressive probing improves parameter information but worsens regulation and may violate actuator constraints. This is the adaptive-control form of the regulation-versus-learning trade-off.

From Chapter 10, full convergence requires a positive-definite information sum. For this two-parameter regression, a finite-window form is

\[ \sum_{i=k}^{k+N-1}\boldsymbol{\phi}(i) \boldsymbol{\phi}^{T}(i)\succeq \mu\mathbf{I}, \qquad \mu>0. \]

A nearly constant output and a nearly constant input can make this matrix poorly conditioned. A bounded multi-sine probe, reference changes, or naturally varying disturbances may supply the missing excitation. The probe amplitude must be selected according to performance and safety constraints rather than solely for estimator speed.

The controller map is particularly sensitive to the estimated input gain. For

\[ u_c=\frac{y_d-\hat{\alpha}y}{\hat{\beta}}, \]

the local sensitivities are

\[ \frac{\partial u_c}{\partial\hat{\alpha}}=-\frac{y}{\hat{\beta}}, \qquad \frac{\partial u_c}{\partial\hat{\beta}}=-\frac{u_c}{\hat{\beta}}. \]

As \( \hat{\beta} \) approaches zero, both sensitivities can become unbounded. A practical implementation therefore uses several layers of protection:

• Projection: enforce physically meaningful parameter bounds.
• Denominator floor: use \( \delta>0 \) in the normalized update.
• Input saturation: bound the command sent to the actuator.
• Rate limiting: restrict changes in controller gains or input.
• Update gating: freeze adaptation when sensors are invalid or saturation dominates.
• Excitation monitoring: inspect the information matrix condition number.
• Fallback controller: revert to validated conservative gains after a fault.

These mechanisms modify the ideal equations, so their effects must be included in implementation-level verification. Stability of the nominal certainty- equivalence formula does not automatically guarantee stability under arbitrary saturation, estimator transients, or unmodeled dynamics.

flowchart TD
  A["Read y(k), reference, and stored estimates"] --> B["Compute desired next output"]
  B --> C["Compute controller gains from alphaHat and betaHat"]
  C --> D["Add bounded probe and apply input saturation"]
  D --> E["Apply u(k) to plant"]
  E --> F["Measure y(k+1)"]
  F --> G["Compute one-step prediction error"]
  G --> H["Normalized parameter update"]
  H --> I["Project estimates into safe bounds"]
  I --> J["Store estimates for next sample"]
  J --> A
        

The corresponding algorithm is:

1. Initialize \( \hat{\alpha}(0) \) and \( \hat{\beta}(0) \) inside \( \Omega \).
2. Measure \( y(k) \) and obtain \( r(k) \).
3. Compute \( y_d(k+1) \) from the selected closed-loop pole.
4. Compute \( u_c(k) \) from the current parameter estimates.
5. Add a bounded probe if required and saturate the actuator command.
6. After sampling \( y(k+1) \), calculate the prediction error.
7. Perform the normalized gradient update and project the result.
8. Repeat using the newly stored estimate.

The implementations below use the plant

\[ y(k+1)=\alpha(k)y(k)+\beta(k)u(k)+v(k), \]

with an abrupt parameter change:

\[ (\alpha,\beta)=\begin{cases} (0.82,0.18), & k<350,\\ (0.68,0.26), & k\ge350. \end{cases} \]

The desired pole is \( p_d=0.55 \), the adaptation gain is \( \gamma=0.45 \), and the estimates are projected onto

\[ 0.40\le\hat{\alpha}\le0.98, \qquad 0.05\le\hat{\beta}\le0.50. \]

A low-amplitude two-frequency probe is added to the certainty-equivalence input. This experiment is intentionally more demanding than convergence to a fixed model: the controller must recover after a change in both the open-loop pole and input gain. The random-noise sequences differ slightly across programming languages, so numerical values are expected to be close rather than identical.

This implementation uses NumPy for arrays, Matplotlib for plots, and the standard csv module for exporting the sample history.

Chapter11_Lesson3.py

"""Chapter11_Lesson3.py

Online parameter estimation driving a certainty-equivalence controller update
for a sampled first-order plant:
    y[k+1] = alpha[k] y[k] + beta[k] u[k] + noise[k]

The identifier uses a projected normalized-gradient update. The controller
recomputes its gains from the latest estimates and adds a small probing signal.
"""

from __future__ import annotations

import csv
import math
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def clamp(value: float, lower: float, upper: float) -> float:
    """Project a scalar onto a closed interval."""
    return max(lower, min(upper, value))


def reference_at(k: int) -> float:
    """Piecewise-constant command used in the numerical experiment."""
    if k < 80:
        return 0.0
    if k < 240:
        return 1.0
    if k < 400:
        return -0.6
    if k < 560:
        return 0.8
    return 0.25


def true_parameters(k: int) -> tuple[float, float]:
    """Unknown plant parameters, including an abrupt change at sample 350."""
    if k < 350:
        return 0.82, 0.18
    return 0.68, 0.26


def rmse(values: np.ndarray) -> float:
    return float(np.sqrt(np.mean(np.square(values))))


def main() -> None:
    samples = 700
    desired_pole = 0.55
    adaptation_gain = 0.45
    normalization = 0.05
    input_limit = 4.0

    alpha_bounds = (0.40, 0.98)
    beta_bounds = (0.05, 0.50)

    rng = np.random.default_rng(7)

    y = np.zeros(samples + 1)
    u = np.zeros(samples)
    r = np.zeros(samples)
    y_desired_next = np.zeros(samples)
    prediction_error = np.zeros(samples)

    alpha_true = np.zeros(samples)
    beta_true = np.zeros(samples)
    alpha_hat = np.zeros(samples + 1)
    beta_hat = np.zeros(samples + 1)

    alpha_hat[0] = 0.60
    beta_hat[0] = 0.10

    for k in range(samples):
        r[k] = reference_at(k)
        alpha_true[k], beta_true[k] = true_parameters(k)

        # Desired one-step response: exact parameters would produce
        # y[k+1] = desired_pole*y[k] + (1-desired_pole)*r[k].
        y_desired_next[k] = desired_pole * y[k] + (1.0 - desired_pole) * r[k]

        # Small deterministic excitation improves identifiability while retaining
        # closed-loop regulation. In practice it must respect performance limits.
        probe = 0.03 * math.sin(2.0 * math.pi * k / 47.0)
        probe += 0.02 * math.sin(2.0 * math.pi * k / 19.0)

        control_nominal = (
            y_desired_next[k] - alpha_hat[k] * y[k]
        ) / beta_hat[k]
        u[k] = clamp(control_nominal + probe, -input_limit, input_limit)

        process_noise = 0.004 * rng.standard_normal()
        y[k + 1] = (
            alpha_true[k] * y[k]
            + beta_true[k] * u[k]
            + process_noise
        )

        y_predicted = alpha_hat[k] * y[k] + beta_hat[k] * u[k]
        prediction_error[k] = y[k + 1] - y_predicted

        denominator = normalization + y[k] ** 2 + u[k] ** 2
        alpha_candidate = (
            alpha_hat[k]
            + adaptation_gain * y[k] * prediction_error[k] / denominator
        )
        beta_candidate = (
            beta_hat[k]
            + adaptation_gain * u[k] * prediction_error[k] / denominator
        )

        alpha_hat[k + 1] = clamp(alpha_candidate, *alpha_bounds)
        beta_hat[k + 1] = clamp(beta_candidate, *beta_bounds)

    tracking_error = y[:-1] - r
    pre_change_rmse = rmse(tracking_error[250:340])
    post_change_rmse = rmse(tracking_error[430:650])

    print(f"Pre-change tracking RMSE:  {pre_change_rmse:.6f}")
    print(f"Post-change tracking RMSE: {post_change_rmse:.6f}")
    print(f"Final alpha estimate:      {alpha_hat[-1]:.6f}")
    print(f"Final beta estimate:       {beta_hat[-1]:.6f}")

    csv_path = Path(__file__).with_name("Chapter11_Lesson3_results.csv")
    with csv_path.open("w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(
            [
                "k",
                "reference",
                "output",
                "input",
                "alpha_true",
                "alpha_hat",
                "beta_true",
                "beta_hat",
                "prediction_error",
            ]
        )
        for k in range(samples):
            writer.writerow(
                [
                    k,
                    r[k],
                    y[k],
                    u[k],
                    alpha_true[k],
                    alpha_hat[k],
                    beta_true[k],
                    beta_hat[k],
                    prediction_error[k],
                ]
            )

    plt.figure(figsize=(10, 5))
    plt.plot(r, label="reference")
    plt.plot(y[:-1], label="output")
    plt.xlabel("sample k")
    plt.ylabel("signal")
    plt.title("Adaptive certainty-equivalence tracking")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(10, 5))
    plt.plot(alpha_true, label="true alpha")
    plt.plot(alpha_hat[:-1], label="estimated alpha")
    plt.plot(beta_true, label="true beta")
    plt.plot(beta_hat[:-1], label="estimated beta")
    plt.xlabel("sample k")
    plt.ylabel("parameter value")
    plt.title("Online parameter estimates")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(10, 4))
    plt.plot(u)
    plt.xlabel("sample k")
    plt.ylabel("u[k]")
    plt.title("Control input")
    plt.grid(True)
    plt.tight_layout()

    plt.show()


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. It writes a CSV file suitable for plotting in Python, MATLAB, GNU Octave, or a spreadsheet.

Chapter11_Lesson3.cpp

// Chapter11_Lesson3.cpp
// Online parameter estimation driving a certainty-equivalence controller update.

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

namespace {

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

double referenceAt(int k) {
    if (k < 80) {
        return 0.0;
    }
    if (k < 240) {
        return 1.0;
    }
    if (k < 400) {
        return -0.6;
    }
    if (k < 560) {
        return 0.8;
    }
    return 0.25;
}

std::pair<double, double> trueParameters(int k) {
    if (k < 350) {
        return {0.82, 0.18};
    }
    return {0.68, 0.26};
}

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

}  // namespace

int main() {
    constexpr int samples = 700;
    constexpr double pi = 3.14159265358979323846;
    constexpr double desiredPole = 0.55;
    constexpr double adaptationGain = 0.45;
    constexpr double normalization = 0.05;
    constexpr double inputLimit = 4.0;

    constexpr double alphaMin = 0.40;
    constexpr double alphaMax = 0.98;
    constexpr double betaMin = 0.05;
    constexpr double betaMax = 0.50;

    std::vector<double> y(samples + 1, 0.0);
    std::vector<double> u(samples, 0.0);
    std::vector<double> reference(samples, 0.0);
    std::vector<double> alphaTrue(samples, 0.0);
    std::vector<double> betaTrue(samples, 0.0);
    std::vector<double> alphaHat(samples + 1, 0.0);
    std::vector<double> betaHat(samples + 1, 0.0);
    std::vector<double> predictionError(samples, 0.0);
    std::vector<double> trackingError(samples, 0.0);

    alphaHat[0] = 0.60;
    betaHat[0] = 0.10;

    std::mt19937 generator(7);
    std::normal_distribution<double> standardNormal(0.0, 1.0);

    for (int k = 0; k < samples; ++k) {
        reference[k] = referenceAt(k);
        const auto [alpha, beta] = trueParameters(k);
        alphaTrue[k] = alpha;
        betaTrue[k] = beta;

        const double desiredNext = desiredPole * y[k]
                                   + (1.0 - desiredPole) * reference[k];
        const double probe = 0.03 * std::sin(2.0 * pi * k / 47.0)
                             + 0.02 * std::sin(2.0 * pi * k / 19.0);
        const double nominalControl = (desiredNext - alphaHat[k] * y[k])
                                      / betaHat[k];
        u[k] = clampValue(nominalControl + probe, -inputLimit, inputLimit);

        const double processNoise = 0.004 * standardNormal(generator);
        y[k + 1] = alphaTrue[k] * y[k] + betaTrue[k] * u[k] + processNoise;

        const double predictedOutput = alphaHat[k] * y[k] + betaHat[k] * u[k];
        predictionError[k] = y[k + 1] - predictedOutput;

        const double denominator = normalization + y[k] * y[k] + u[k] * u[k];
        const double alphaCandidate = alphaHat[k]
                                      + adaptationGain * y[k]
                                            * predictionError[k] / denominator;
        const double betaCandidate = betaHat[k]
                                     + adaptationGain * u[k]
                                           * predictionError[k] / denominator;

        alphaHat[k + 1] = clampValue(alphaCandidate, alphaMin, alphaMax);
        betaHat[k + 1] = clampValue(betaCandidate, betaMin, betaMax);
        trackingError[k] = y[k] - reference[k];
    }

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Pre-change tracking RMSE:  "
              << rmse(trackingError, 250, 340) << '\n';
    std::cout << "Post-change tracking RMSE: "
              << rmse(trackingError, 430, 650) << '\n';
    std::cout << "Final alpha estimate:      " << alphaHat.back() << '\n';
    std::cout << "Final beta estimate:       " << betaHat.back() << '\n';

    std::ofstream csv("Chapter11_Lesson3_results_cpp.csv");
    if (!csv) {
        std::cerr << "Could not create output CSV file.\n";
        return 1;
    }
    csv << "k,reference,output,input,alpha_true,alpha_hat,beta_true,beta_hat,prediction_error\n";
    csv << std::setprecision(12);
    for (int k = 0; k < samples; ++k) {
        csv << k << ',' << reference[k] << ',' << y[k] << ',' << u[k] << ','
            << alphaTrue[k] << ',' << alphaHat[k] << ',' << betaTrue[k] << ','
            << betaHat[k] << ',' << predictionError[k] << '\n';
    }

    return 0;
}

The Java implementation uses standard arrays, java.util.Random, and buffered CSV output. No external dependency is required.

Chapter11_Lesson3.java

// Chapter11_Lesson3.java
// Online parameter estimation driving a certainty-equivalence controller update.

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

public final class Chapter11_Lesson3 {
    private static final int SAMPLES = 700;
    private static final double DESIRED_POLE = 0.55;
    private static final double ADAPTATION_GAIN = 0.45;
    private static final double NORMALIZATION = 0.05;
    private static final double INPUT_LIMIT = 4.0;

    private static final double ALPHA_MIN = 0.40;
    private static final double ALPHA_MAX = 0.98;
    private static final double BETA_MIN = 0.05;
    private static final double BETA_MAX = 0.50;

    private Chapter11_Lesson3() {
        // Utility class.
    }

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

    private static double referenceAt(int k) {
        if (k < 80) {
            return 0.0;
        }
        if (k < 240) {
            return 1.0;
        }
        if (k < 400) {
            return -0.6;
        }
        if (k < 560) {
            return 0.8;
        }
        return 0.25;
    }

    private static double[] trueParameters(int k) {
        if (k < 350) {
            return new double[] {0.82, 0.18};
        }
        return new double[] {0.68, 0.26};
    }

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

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

        double[] y = new double[SAMPLES + 1];
        double[] u = new double[SAMPLES];
        double[] reference = new double[SAMPLES];
        double[] alphaTrue = new double[SAMPLES];
        double[] betaTrue = new double[SAMPLES];
        double[] alphaHat = new double[SAMPLES + 1];
        double[] betaHat = new double[SAMPLES + 1];
        double[] predictionError = new double[SAMPLES];
        double[] trackingError = new double[SAMPLES];

        alphaHat[0] = 0.60;
        betaHat[0] = 0.10;

        Random random = new Random(7L);

        for (int k = 0; k < SAMPLES; k++) {
            reference[k] = referenceAt(k);
            double[] parameters = trueParameters(k);
            alphaTrue[k] = parameters[0];
            betaTrue[k] = parameters[1];

            double desiredNext = DESIRED_POLE * y[k]
                    + (1.0 - DESIRED_POLE) * reference[k];
            double probe = 0.03 * Math.sin(2.0 * Math.PI * k / 47.0)
                    + 0.02 * Math.sin(2.0 * Math.PI * k / 19.0);
            double nominalControl = (desiredNext - alphaHat[k] * y[k])
                    / betaHat[k];
            u[k] = clamp(nominalControl + probe, -INPUT_LIMIT, INPUT_LIMIT);

            double processNoise = 0.004 * random.nextGaussian();
            y[k + 1] = alphaTrue[k] * y[k]
                    + betaTrue[k] * u[k]
                    + processNoise;

            double predictedOutput = alphaHat[k] * y[k] + betaHat[k] * u[k];
            predictionError[k] = y[k + 1] - predictedOutput;

            double denominator = NORMALIZATION + y[k] * y[k] + u[k] * u[k];
            double alphaCandidate = alphaHat[k]
                    + ADAPTATION_GAIN * y[k] * predictionError[k] / denominator;
            double betaCandidate = betaHat[k]
                    + ADAPTATION_GAIN * u[k] * predictionError[k] / denominator;

            alphaHat[k + 1] = clamp(alphaCandidate, ALPHA_MIN, ALPHA_MAX);
            betaHat[k + 1] = clamp(betaCandidate, BETA_MIN, BETA_MAX);
            trackingError[k] = y[k] - reference[k];
        }

        System.out.printf("Pre-change tracking RMSE:  %.6f%n",
                rmse(trackingError, 250, 340));
        System.out.printf("Post-change tracking RMSE: %.6f%n",
                rmse(trackingError, 430, 650));
        System.out.printf("Final alpha estimate:      %.6f%n", alphaHat[SAMPLES]);
        System.out.printf("Final beta estimate:       %.6f%n", betaHat[SAMPLES]);

        Path csvPath = Path.of("Chapter11_Lesson3_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                csvPath, StandardCharsets.UTF_8)) {
            writer.write("k,reference,output,input,alpha_true,alpha_hat,"
                    + "beta_true,beta_hat,prediction_error");
            writer.newLine();
            for (int k = 0; k < SAMPLES; k++) {
                writer.write(String.format(
                        Locale.US,
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                        k,
                        reference[k],
                        y[k],
                        u[k],
                        alphaTrue[k],
                        alphaHat[k],
                        betaTrue[k],
                        betaHat[k],
                        predictionError[k]));
                writer.newLine();
            }
        }
    }
}

The MATLAB script simulates, plots, and exports the complete history. For a Simulink realization, use Unit Delay blocks for the plant output and parameter estimates; Sum, Product, and Divide blocks for the prediction and update laws; and Saturation blocks for the parameters and actuator command. A MATLAB Function block may alternatively contain the estimator and certainty-equivalence equations.

Chapter11_Lesson3.m

% Chapter11_Lesson3.m
% Online parameter estimation driving a certainty-equivalence controller update.

clear;
clc;
close all;

samples = 700;
desiredPole = 0.55;
adaptationGain = 0.45;
normalization = 0.05;
inputLimit = 4.0;

alphaBounds = [0.40, 0.98];
betaBounds = [0.05, 0.50];

rng(7, 'twister');

y = zeros(samples + 1, 1);
u = zeros(samples, 1);
r = zeros(samples, 1);
alphaTrue = zeros(samples, 1);
betaTrue = zeros(samples, 1);
alphaHat = zeros(samples + 1, 1);
betaHat = zeros(samples + 1, 1);
predictionError = zeros(samples, 1);

alphaHat(1) = 0.60;
betaHat(1) = 0.10;

for k = 1:samples
    index = k - 1;

    if index < 80
        r(k) = 0.0;
    elseif index < 240
        r(k) = 1.0;
    elseif index < 400
        r(k) = -0.6;
    elseif index < 560
        r(k) = 0.8;
    else
        r(k) = 0.25;
    end

    if index < 350
        alphaTrue(k) = 0.82;
        betaTrue(k) = 0.18;
    else
        alphaTrue(k) = 0.68;
        betaTrue(k) = 0.26;
    end

    desiredNext = desiredPole * y(k) + (1 - desiredPole) * r(k);
    probe = 0.03 * sin(2 * pi * index / 47) ...
          + 0.02 * sin(2 * pi * index / 19);

    nominalControl = (desiredNext - alphaHat(k) * y(k)) / betaHat(k);
    u(k) = min(inputLimit, max(-inputLimit, nominalControl + probe));

    processNoise = 0.004 * randn();
    y(k + 1) = alphaTrue(k) * y(k) ...
             + betaTrue(k) * u(k) ...
             + processNoise;

    predictedOutput = alphaHat(k) * y(k) + betaHat(k) * u(k);
    predictionError(k) = y(k + 1) - predictedOutput;

    denominator = normalization + y(k)^2 + u(k)^2;
    alphaCandidate = alphaHat(k) ...
                   + adaptationGain * y(k) * predictionError(k) / denominator;
    betaCandidate = betaHat(k) ...
                  + adaptationGain * u(k) * predictionError(k) / denominator;

    alphaHat(k + 1) = min(alphaBounds(2), max(alphaBounds(1), alphaCandidate));
    betaHat(k + 1) = min(betaBounds(2), max(betaBounds(1), betaCandidate));
end

trackingError = y(1:end-1) - r;
preChangeRmse = sqrt(mean(trackingError(251:340).^2));
postChangeRmse = sqrt(mean(trackingError(431:650).^2));

fprintf('Pre-change tracking RMSE:  %.6f\n', preChangeRmse);
fprintf('Post-change tracking RMSE: %.6f\n', postChangeRmse);
fprintf('Final alpha estimate:      %.6f\n', alphaHat(end));
fprintf('Final beta estimate:       %.6f\n', betaHat(end));

results = table((0:samples-1)', r, y(1:end-1), u, alphaTrue, ...
    alphaHat(1:end-1), betaTrue, betaHat(1:end-1), predictionError, ...
    'VariableNames', {'k', 'reference', 'output', 'input', 'alpha_true', ...
    'alpha_hat', 'beta_true', 'beta_hat', 'prediction_error'});
writetable(results, 'Chapter11_Lesson3_results_matlab.csv');

figure;
plot(0:samples-1, r, '--', 'LineWidth', 1.2);
hold on;
plot(0:samples-1, y(1:end-1), 'LineWidth', 1.2);
grid on;
xlabel('sample k');
ylabel('signal');
title('Adaptive certainty-equivalence tracking');
legend('reference', 'output', 'Location', 'best');

figure;
plot(0:samples-1, alphaTrue, '--', 'LineWidth', 1.2);
hold on;
plot(0:samples-1, alphaHat(1:end-1), 'LineWidth', 1.2);
plot(0:samples-1, betaTrue, '--', 'LineWidth', 1.2);
plot(0:samples-1, betaHat(1:end-1), 'LineWidth', 1.2);
grid on;
xlabel('sample k');
ylabel('parameter value');
title('Online parameter estimates');
legend('true alpha', 'estimated alpha', 'true beta', 'estimated beta', ...
    'Location', 'best');

figure;
plot(0:samples-1, u, 'LineWidth', 1.2);
grid on;
xlabel('sample k');
ylabel('u[k]');
title('Control input');

% Simulink implementation guidance:
% 1. Use a Unit Delay block for y[k] and for each parameter estimate.
% 2. Compute predictionError = y[k+1] - alphaHat*y[k] - betaHat*u[k].
% 3. Implement the normalized-gradient law with Product, Sum, and Divide blocks.
% 4. Apply Saturation blocks to alphaHat, betaHat, and u.
% 5. Recompute nominalControl from the projected estimates at every sample.

The notebook contains the same adaptive loop and produces tracking, parameter, and input plots using native Wolfram Language functions.

Chapter11_Lesson3.nb

Notebook[{
  Cell["Chapter 11 - Lesson 3", "Title"],
  Cell["Online Parameter Estimation Driving Controller Update", "Subtitle"],
  Cell[BoxData["(* Chapter11_Lesson3.wl
   Online parameter estimation driving a certainty-equivalence controller update. *)
ClearAll[\"Global`*\"];
samples = 700; desiredPole = 0.55; adaptationGain = 0.45;
normalization = 0.05; inputLimit = 4.0; SeedRandom[7];
alphaBounds = {0.40, 0.98}; betaBounds = {0.05, 0.50};
referenceAt[k_Integer] := Piecewise[{
  {0.0, k < 80}, {1.0, k < 240}, {-0.6, k < 400},
  {0.8, k < 560}}, 0.25];
trueParameters[k_Integer] := If[k < 350, {0.82, 0.18}, {0.68, 0.26}];
clip[x_, {lower_, upper_}] := Clip[x, {lower, upper}];
y = ConstantArray[0.0, samples + 1]; u = ConstantArray[0.0, samples];
r = ConstantArray[0.0, samples]; alphaTrue = ConstantArray[0.0, samples];
betaTrue = ConstantArray[0.0, samples]; predictionError = ConstantArray[0.0, samples];
alphaHat = ConstantArray[0.0, samples + 1]; betaHat = ConstantArray[0.0, samples + 1];
alphaHat[[1]] = 0.60; betaHat[[1]] = 0.10;
Do[
  Module[{index, desiredNext, probe, nominalControl, processNoise,
    predictedOutput, denominator, alphaCandidate, betaCandidate, parameters},
   index = k - 1; r[[k]] = referenceAt[index];
   parameters = trueParameters[index];
   alphaTrue[[k]] = parameters[[1]]; betaTrue[[k]] = parameters[[2]];
   desiredNext = desiredPole*y[[k]] + (1 - desiredPole)*r[[k]];
   probe = 0.03*Sin[2*Pi*index/47] + 0.02*Sin[2*Pi*index/19];
   nominalControl = (desiredNext - alphaHat[[k]]*y[[k]])/betaHat[[k]];
   u[[k]] = clip[nominalControl + probe, {-inputLimit, inputLimit}];
   processNoise = RandomVariate[NormalDistribution[0, 0.004]];
   y[[k + 1]] = alphaTrue[[k]]*y[[k]] + betaTrue[[k]]*u[[k]] + processNoise;
   predictedOutput = alphaHat[[k]]*y[[k]] + betaHat[[k]]*u[[k]];
   predictionError[[k]] = y[[k + 1]] - predictedOutput;
   denominator = normalization + y[[k]]^2 + u[[k]]^2;
   alphaCandidate = alphaHat[[k]] +
     adaptationGain*y[[k]]*predictionError[[k]]/denominator;
   betaCandidate = betaHat[[k]] +
     adaptationGain*u[[k]]*predictionError[[k]]/denominator;
   alphaHat[[k + 1]] = clip[alphaCandidate, alphaBounds];
   betaHat[[k + 1]] = clip[betaCandidate, betaBounds];],
  {k, 1, samples}];
trackingError = Most[y] - r; rmse[data_] := Sqrt[Mean[data^2]];
preChangeRmse = rmse[trackingError[[251 ;; 340]]];
postChangeRmse = rmse[trackingError[[431 ;; 650]]];
Print[\"Pre-change tracking RMSE:  \", NumberForm[preChangeRmse, {8, 6}]];
Print[\"Post-change tracking RMSE: \", NumberForm[postChangeRmse, {8, 6}]];
Print[\"Final alpha estimate:      \", NumberForm[Last[alphaHat], {8, 6}]];
Print[\"Final beta estimate:       \", NumberForm[Last[betaHat], {8, 6}]];
trackingPlot = ListLinePlot[{r, Most[y]},
  PlotLegends -> {\"reference\", \"output\"}, Frame -> True,
  FrameLabel -> {\"sample k\", \"signal\"},
  PlotLabel -> \"Adaptive certainty-equivalence tracking\", ImageSize -> Large];
parameterPlot = ListLinePlot[{alphaTrue, Most[alphaHat], betaTrue, Most[betaHat]},
  PlotLegends -> {\"true alpha\", \"estimated alpha\", \"true beta\", \"estimated beta\"},
  Frame -> True, FrameLabel -> {\"sample k\", \"parameter value\"},
  PlotLabel -> \"Online parameter estimates\", ImageSize -> Large];
inputPlot = ListLinePlot[u, Frame -> True, FrameLabel -> {\"sample k\", \"u[k]\"},
  PlotLabel -> \"Control input\", ImageSize -> Large];
Column[{trackingPlot, parameterPlot, inputPlot}]
"], "Input"]}, WindowTitle -> "Chapter11_Lesson3"]        

During the first operating regime, the prediction error corrects the initially inaccurate model and the controller gains move toward values consistent with the desired pole. At sample 350, the old model becomes incorrect. The prediction error increases, producing renewed parameter motion and automatic gain retuning.

The estimates need not settle exactly at the true values. The reference becomes mostly piecewise constant, the feedback action suppresses output variation, and the probing signal is deliberately small. Therefore the regressor may not satisfy a strong persistent-excitation condition over every interval. Nevertheless, the controller can achieve small tracking error because the relevant combination \( \boldsymbol{\phi}^{T}\tilde{\boldsymbol{\theta}} \) becomes small along the realized closed-loop trajectory.

This distinction is central to adaptive control: identification quality should be evaluated in relation to the control objective, but inaccurate estimates must not be allowed to generate unsafe gains. Projection, saturation, monitoring, and a fallback policy therefore remain necessary even when nominal tracking is good.

Problem 1 (Controller-Gain Derivation): For the estimated plant \( \hat{y}(k+1)=\hat{\alpha}y(k)+\hat{\beta}u(k) \), derive feedback and feedforward gains that impose \( y_d(k+1)=p_d y(k)+(1-p_d)r(k) \).

Solution: Set the estimated response equal to the desired response:

\[ \hat{\alpha}y+\hat{\beta}u=p_d y+(1-p_d)r. \]

Solving for the input gives

\[ u=-\frac{\hat{\alpha}-p_d}{\hat{\beta}}y +\frac{1-p_d}{\hat{\beta}}r. \]

Hence

\[ \hat{k}_y=\frac{\hat{\alpha}-p_d}{\hat{\beta}}, \qquad \hat{k}_r=\frac{1-p_d}{\hat{\beta}}. \]

Problem 2 (Prediction-Error Relation): Assume the true plant is noise free and define \( \tilde{\boldsymbol{\theta}}= \boldsymbol{\theta}-\hat{\boldsymbol{\theta}} \). Show that the prediction error is a projection of parameter error onto the regressor.

Solution:

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

Therefore zero prediction error constrains only the component of parameter error visible in the current regressor direction. Repeatedly varying independent regressor directions are needed to identify the full vector.

Problem 3 (Allowed Adaptation-Gain Range): Using the expression for \( V(k+1)-V(k) \) in Section 6, find a simple constant range of \( \gamma \) that guarantees nonincrease of \( V \) for every regressor.

Solution: The factor that must remain nonnegative is

\[ 1-\frac{\gamma\|\boldsymbol{\phi}\|^{2}} {2(\delta+\|\boldsymbol{\phi}\|^{2})}. \]

Because

\[ 0\le\frac{\|\boldsymbol{\phi}\|^{2}} {\delta+\|\boldsymbol{\phi}\|^{2}}<1, \]

any constant choice \( 0<\gamma<2 \) is sufficient. This is a descent condition for the ideal estimator; it is not by itself a complete closed-loop robustness guarantee.

Problem 4 (Why Projection of the Input Gain Is Essential): Suppose \( \hat{\beta}(k)=0.005 \), \( y_d(k+1)-\hat{\alpha}(k)y(k)=0.5 \), and the actuator limit is \( 4 \). Compute the nominal input and explain the conflict.

Solution:

\[ u_c(k)=\frac{0.5}{0.005}=100. \]

The actuator can apply only the saturated value \( 4 \), so the identifier-controller equations used in the nominal derivation no longer describe the actual loop. Enforcing \( \hat{\beta}\ge\beta_{\min} \), selecting physically meaningful bounds, and gating adaptation during severe saturation reduce this risk.

Problem 5 (Excitation Matrix): Suppose over a window the data satisfy \( u(i)=c y(i) \) for every sample. Show why the two parameters cannot be independently identified from that window.

Solution: Every regressor has the form

\[ \boldsymbol{\phi}(i)=y(i) \begin{bmatrix}1\\c\end{bmatrix}. \]

Hence the information matrix is

\[ \sum_i\boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i) =\left(\sum_i y^{2}(i)\right) \begin{bmatrix}1\\c\end{bmatrix} \begin{bmatrix}1&c\end{bmatrix}. \]

This matrix has rank one, whereas two independent parameter directions require rank two. Only the combination \( \alpha+c\beta \) is visible.

Problem 6 (Slow Parameter Variation): Let \( \boldsymbol{\theta}(k+1)=\boldsymbol{\theta}(k)+ \boldsymbol{\nu}(k) \). Derive the parameter-error recursion before projection.

Solution: Since

\[ \tilde{\boldsymbol{\theta}}(k+1)= \boldsymbol{\theta}(k+1)-\hat{\boldsymbol{\theta}}(k+1), \]

substitution of the estimator update gives

\[ \tilde{\boldsymbol{\theta}}(k+1)= \tilde{\boldsymbol{\theta}}(k)+\boldsymbol{\nu}(k)- \frac{\gamma\boldsymbol{\phi}(k)\varepsilon(k+1)} {\delta+\|\boldsymbol{\phi}(k)\|^{2}}. \]

The variation term continually injects parameter error. Consequently, a tracking estimator generally reaches a bounded error neighborhood rather than exact convergence; increasing adaptation speed can reduce lag but increases sensitivity to noise.

An indirect adaptive controller closes two interacting loops. The prediction-error loop updates a plant model, while the certainty-equivalence loop maps that model into controller gains. For the sampled first-order plant, a normalized projected gradient estimator provides bounded parameter motion, and the estimated parameters determine feedback and reference gains that impose a desired nominal pole.

The descent proof establishes nonincrease of an ideal parameter-error energy for \( 0<\gamma<2 \), but full parameter convergence still requires persistent excitation. Closed-loop regulation may remove that excitation, so good tracking does not imply exact identification. Safe implementation therefore combines model-based gain updates with projection, saturation, excitation management, update gating, monitoring, and a fallback controller.