A controller is adaptive when some internal quantity changes online in response to measured data. A controller is learning-based when that change is produced by a statistical, optimization-based, function-approximation, or trajectory-learning mechanism. A controller is data-driven when controller synthesis or prediction is performed directly from measured trajectories, possibly without first producing a conventional parametric plant model.

These labels overlap. Recursive least squares in a self-tuning regulator is already an online learning algorithm. A neural adaptive controller is an adaptive controller with a nonlinear function approximator. Virtual Reference Feedback Tuning is a direct batch data-based design method. Behavioral methods such as data-enabled predictive control use measured trajectory matrices as an implicit model. The essential control question remains:

\[ \begin{aligned} &\text{How can data-dependent updates improve performance} \\ &\text{while preserving boundedness, constraint satisfaction,} \\ &\text{and closed-loop stability?} \end{aligned} \]

flowchart TD
  A["Measured input, state, \nand output data"] --> B["Learning mechanism"]
  B --> C1["Online parameter \nestimator"]
  B --> C2["Direct controller \ntuning"]
  B --> C3["Function \napproximator"]
  B --> C4["Trajectory or \nbehavioral representation"]
  C1 --> D["Adaptive controller"]
  C2 --> D
  C3 --> D
  C4 --> D
  D --> E["Plant and constraints"]
  E --> A
  D --> F["Safety layer: projection, \nsaturation, monitoring"]
        

Let the measured information vector at sample \( k \) be \( z_k \), let the controller be \( u_k=\pi_{\eta_k}(z_k) \), and let \( \eta_k \) denote learned controller parameters, model parameters, neural weights, or coefficients of a data representation. A broad class of online updates can be written as

\[ \eta_{k+1}=\operatorname{Proj}_{\Omega}\!\left( \eta_k-\alpha_k\widehat{g}_k\right), \qquad \widehat{g}_k \approx \nabla_{\eta}\ell_k(\eta_k). \]

Here \( \ell_k \) may be prediction error, tracking error, a finite-horizon control cost, or a Bellman residual. The projection set \( \Omega \) incorporates prior engineering knowledge. Classical gradient MRAC, RLS-based STR, neural adaptive control, and many policy-gradient updates all fit this structure, but their stability analyses differ because \( \widehat{g}_k \) is generated inside a feedback loop rather than from independent data.

For the lesson laboratory, consider the unknown discrete-time nonlinear plant

\[ y_{k+1}=a_*y_k+b_*u_k+c_*y_k^3+d_*+w_k =\theta_*^{\mathsf T}\phi_k+w_k, \]

\[ \theta_*=[a_*,b_*,c_*,d_*]^{\mathsf T}, \qquad \phi_k=[y_k,u_k,y_k^3,1]^{\mathsf T}. \]

The dynamics are nonlinear in the measured output but linear in the unknown parameters. This structure allows an online learner with a transparent Lyapunov and least-squares interpretation while still illustrating a learning-based nonlinear adaptive controller.

Batch learning uses a fixed dataset

\[ \mathcal D_N=\{(z_i,u_i,y_{i+1})\}_{i=0}^{N-1} \]

and minimizes an empirical objective such as

\[ \widehat J_N(\eta)=\frac{1}{N}\sum_{i=0}^{N-1} \ell\!\left(y_{i+1},\widehat y_{i+1}(\eta)\right) +\lambda_{\eta}\mathcal R(\eta). \]

The result is deployed after training. Online learning instead updates the parameter after each new transition. In control, online operation creates a dual effect: the current input regulates the plant and determines the quality of future data. Consequently, a purely greedy controller may suppress the excitation needed for learning.

Direct data-based controller design bypasses an explicit parametric plant model. Examples include Iterative Feedback Tuning, Virtual Reference Feedback Tuning, data-driven LMIs, and behavioral trajectory optimization. “No explicit model” does not mean “no assumptions.” Typical assumptions concern controllability, data richness, bounded noise, consistent input-output dimensions, or existence of a stabilizing controller class.

Given transitions up to sample \( k \), define the exponentially weighted regularized loss

\[ J_k(\theta)=\sum_{i=0}^{k-1}\lambda^{k-1-i} \left(y_{i+1}-\theta^{\mathsf T}\phi_i\right)^2 +(\theta-\theta_0)^{\mathsf T}P_0^{-1}(\theta-\theta_0), \]

with \( 0<\lambda\leq 1 \). The minimizer can be updated recursively. Define the one-step prediction error

\[ \varepsilon_{k+1}=y_{k+1}-\widehat\theta_k^{\mathsf T}\phi_k. \]

The RLS recursion is

\[ K_k=\frac{P_k\phi_k}{\lambda+\phi_k^{\mathsf T}P_k\phi_k}, \]

\[ \widehat\theta_{k+1}=\operatorname{Proj}_{\Omega} \left(\widehat\theta_k+K_k\varepsilon_{k+1}\right), \]

\[ P_{k+1}=\frac{1}{\lambda}\left( P_k-K_k\phi_k^{\mathsf T}P_k\right). \]

The matrix \( P_k \) acts as an inverse information matrix. Large eigenvalues imply that some parameter directions are weakly informed by the data. Forgetting \( \lambda<1 \) allows slow parameter variation to be tracked but increases sensitivity to noise. Projection keeps the estimate inside a physically plausible compact set.

Introduce the information matrix and information vector

\[ S_k=P_k^{-1}, \qquad q_k=S_k\widehat\theta_k. \]

From the weighted least-squares objective,

\[ S_{k+1}=\lambda S_k+\phi_k\phi_k^{\mathsf T}, \qquad q_{k+1}=\lambda q_k+\phi_k y_{k+1}. \]

Applying the matrix inversion lemma to \( S_{k+1} \) yields

\[ (\lambda S_k+\phi_k\phi_k^{\mathsf T})^{-1} =\frac{1}{\lambda}\left[P_k- \frac{P_k\phi_k\phi_k^{\mathsf T}P_k} {\lambda+\phi_k^{\mathsf T}P_k\phi_k}\right]. \]

Substitution into \( \widehat\theta_{k+1}=P_{k+1}q_{k+1} \) gives the gain and parameter recursions of Section 4. Thus RLS is not an arbitrary adaptation rule; it is an exact recursive optimizer for a quadratic online learning objective.

In the noise-free, constant-parameter case with \( \lambda=1 \) and no projection, the parameter error \( \widetilde\theta_k=\theta_*-\widehat\theta_k \) satisfies a non-increasing information-weighted energy. A standard identity is

\[ V_{k+1}-V_k=-\frac{ (\phi_k^{\mathsf T}\widetilde\theta_k)^2} {1+\phi_k^{\mathsf T}P_k\phi_k}\leq 0, \qquad V_k=\widetilde\theta_k^{\mathsf T}P_k^{-1}\widetilde\theta_k. \]

This proves boundedness of the estimation error in the corresponding metric; convergence to zero additionally requires sufficient excitation.

Let \( \widehat\theta_k=[\widehat a_k,\widehat b_k, \widehat c_k,\widehat d_k]^{\mathsf T} \). The nominal inverse control that would force the learned one-step predictor to equal the next reference is

\[ u_k^{\mathrm{nom} }=\frac{1}{\widehat b_k} \left(r_{k+1}-\widehat a_k y_k-\widehat c_k y_k^3- \widehat d_k\right). \]

Direct division is unsafe when \( \widehat b_k \) is small. We therefore use a regularized inverse, probing excitation, and saturation:

\[ v_k=r_{k+1}-\widehat a_k y_k- \widehat c_k y_k^3-\widehat d_k, \]

\[ u_k=\operatorname{sat}_{u_{\max} }\!\left( \frac{\widehat b_k}{\widehat b_k^2+\rho}v_k+d_k^{\mathrm{exc} } \right), \qquad \rho>0. \]

Projection enforces a known control direction by requiring \( \widehat b_k\geq b_{\min}>0 \). The excitation \( d_k^{\mathrm{exc} } \) is applied only during early learning in the laboratory. Saturation prevents an imperfect predictor from commanding excessive input.

Assume temporarily that saturation is inactive and define \( \widetilde\theta_k=\theta_*-\widehat\theta_k \). The learned predictor is

\[ \widehat y_{k+1}=\widehat\theta_k^{\mathsf T}\phi_k. \]

Since the regularized inverse gives

\[ \widehat y_{k+1}-r_{k+1} =-\frac{\rho}{\widehat b_k^2+\rho}v_k +\widehat b_k d_k^{\mathrm{exc} }, \]

the actual tracking error satisfies

\[ e_{k+1}=y_{k+1}-r_{k+1} =\widetilde\theta_k^{\mathsf T}\phi_k -\frac{\rho}{\widehat b_k^2+\rho}v_k +\widehat b_kd_k^{\mathrm{exc} }+w_k. \]

If \( \widehat b_k\geq b_{\min}>0 \), then

\[ |e_{k+1}|\leq \|\widetilde\theta_k\|\,\|\phi_k\| +\frac{\rho}{b_{\min}^2+\rho}|v_k| +|\widehat b_k|\,|d_k^{\mathrm{exc} }|+|w_k|. \]

This inequality exposes four performance terms: learning error, regularization bias, deliberate exploration, and disturbance/noise. It also explains why reducing \( \rho \) is not automatically beneficial: a very small regularizer amplifies errors in the estimated input gain and may activate saturation.

Parameter convergence requires the regressor to explore all parameter directions. A discrete-time persistent-excitation condition is: there exist \( N\in\mathbb N \) and \( \alpha>0 \) such that

\[ \sum_{i=k}^{k+N-1}\phi_i\phi_i^{\mathsf T}\succeq \alpha I \quad \text{for every }k. \]

If the reference is constant and regulation becomes nearly perfect, then \( y_k \) and \( u_k \) may become almost constant. The information matrix becomes poorly conditioned, and different combinations of \( a_* \) and \( c_* \) can explain the same narrow operating region. Good tracking can therefore coexist with inaccurate individual parameter estimates.

Practical mechanisms for improving data richness include finite probing signals, reference variation, command shaping, safe exploration, concurrent learning from stored informative data, and composite adaptation using both tracking and prediction errors. Excitation must remain compatible with input, state, and safety constraints.

For a nonlinear uncertainty \( \Delta(z) \), a neural or basis-function approximation can be written as

\[ \Delta(z)=W_*^{\mathsf T}\sigma(z)+\epsilon(z), \qquad \|\epsilon(z)\|\leq\bar\epsilon \text{ on a compact domain}. \]

In a reference-model error system \( \dot e=A_m e+B(\widetilde W^{\mathsf T}\sigma+ \epsilon) \), choose

\[ V=e^{\mathsf T}Pe+ \operatorname{tr}\!\left(\widetilde W^{\mathsf T} \Gamma^{-1}\widetilde W\right), \qquad A_m^{\mathsf T}P+PA_m=-Q, \]

and the weight law

\[ \dot{\widehat W}=\Gamma\sigma(z)e^{\mathsf T}PB. \]

The cross term cancels, giving the bound

\[ \dot V\leq-\lambda_{\min}(Q)\|e\|^2 +2\|e\|\,\|PB\|\bar\epsilon. \]

Therefore the tracking error is ultimately bounded by a radius proportional to the approximation error. Deep networks complicate this proof because the model is nonlinear in all weights. Common safety-oriented architectures keep a proven baseline controller, learn only a bounded residual, use projection or leakage, and reject learned actions that violate a control-barrier, Lyapunov, or constraint filter.

Let \( u^d=\{u_0^d,\ldots,u_{T-1}^d\} \) be a recorded input sequence. Its depth-\( L \) block Hankel matrix is

\[ H_L(u^d)=\begin{bmatrix} u_0^d & u_1^d & \cdots & u_{T-L}^d\\ u_1^d & u_2^d & \cdots & u_{T-L+1}^d\\ \vdots & \vdots & \ddots & \vdots\\ u_{L-1}^d & u_L^d & \cdots & u_{T-1}^d \end{bmatrix}. \]

For a controllable LTI system, sufficiently persistently exciting data allow every length-\( L \) input-output trajectory to be expressed as

\[ \begin{bmatrix}u\\y\end{bmatrix} =\begin{bmatrix}H_L(u^d)\\H_L(y^d)\end{bmatrix}g \]

for some coefficient vector \( g \). This is an implicit model obtained from trajectories rather than identified state matrices. Data-enabled predictive control optimizes future inputs and outputs subject to such trajectory constraints. Data-driven LMI methods similarly replace unknown system matrices by algebraic relations constructed from measured state and input data.

These methods are linked to, but not identical with, adaptive control. A batch trajectory controller becomes adaptive when the dataset, regularization, uncertainty description, or controller is updated online in response to new measurements.

Three distinct claims are often confused:

Prediction generalization: the learned predictor performs well on trajectories not used for fitting.

Closed-loop stability: the feedback interconnection keeps states, signals, and learned parameters bounded and produces an appropriate convergence property.

Constraint safety: inputs and states remain inside prescribed admissible sets during learning and deployment.

Small one-step prediction error does not by itself prove stability. Likewise, a Lyapunov proof for an unconstrained adaptive law does not guarantee that actuator saturation or state limits are respected. A deployable design should state assumptions, define the certified operating region, bound the learning update, monitor excitation and covariance conditioning, and provide a fallback controller.

\[ u_k=\mathcal S\!\left(u_k^{\mathrm{learn} },x_k,\mathcal X, \mathcal U,V\right), \]

where \( \mathcal S \) is a safety filter that modifies or rejects the learned command when constraints or a Lyapunov decrease condition would be violated.

flowchart TD
  A["Define tracking, constraints, and uncertainty"] --> B["Choose a proven baseline controller"]
  B --> C["Select learned quantity: model, residual, gains, or trajectories"]
  C --> D["Define loss, update law, projection, and regularization"]
  D --> E["Plan informative but safe excitation"]
  E --> F["Prove boundedness or state explicit guarantees"]
  F --> G["Test noise, delay, saturation, and distribution shift"]
  G --> H["Deploy with monitoring and fallback"]
  H --> I["Log new data and validate before updates"]
  I --> D
        

The workflow deliberately begins with control requirements rather than a choice of machine-learning architecture. The learned component should be the smallest component that materially improves performance while preserving a tractable safety argument.

Python: NumPy and SciPy support linear algebra and online estimation; python-control supports classical control models; PyTorch and JAX support differentiable neural controllers. The laboratory uses NumPy and Matplotlib so the adaptive recursion remains explicit.

C++: Eigen is commonly used for fixed- and dynamic-size matrices; LibTorch can execute PyTorch networks; dlib provides optimization and learning utilities. The laboratory uses only the C++17 standard library and explicit four-dimensional matrix operations.

Java: EJML supports numerical linear algebra, while ND4J and Deeplearning4j support tensor and neural-network computations. The laboratory uses standard arrays to expose every RLS operation.

MATLAB/Simulink: Control System Toolbox, System Identification Toolbox, Deep Learning Toolbox, Reinforcement Learning Toolbox, and Simulink are relevant. The script includes a local adaptive-step function that can be moved into a Simulink MATLAB Function block with Unit Delay, Saturation, plant, and Scope blocks.

Wolfram Mathematica: built-in matrix operations, LinearModelFit, NDSolve, optimization functions, and neural-network functions such as NetTrain support symbolic and numerical analysis. The notebook implements the recursion directly.

The Python experiment learns the four plant coefficients online. Two sinusoidal probing components are applied only during the first 250 samples. The controller then relies on the learned model, projection, regularized inversion, and saturation.

Chapter28_Lesson2.py

"""Chapter28_Lesson2.py
Learning-based/data-driven adaptive control using online RLS and a
certainty-equivalent one-step controller with projection and saturation.
"""

from __future__ import annotations

import math
from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    steps: int = 700
    forgetting: float = 0.998
    regularization: float = 0.05
    input_limit: float = 3.0
    noise_std: float = 0.005
    seed: int = 7


def project(theta: np.ndarray) -> np.ndarray:
    """Componentwise projection onto conservative parameter bounds."""
    lower = np.array([-0.95, 0.10, -0.40, -0.50], dtype=float)
    upper = np.array([0.95, 1.50, 0.40, 0.50], dtype=float)
    return np.clip(theta, lower, upper)


def reference(k: int) -> float:
    if k < 100:
        return 0.0
    if k < 220:
        return 1.2
    if k < 340:
        return -1.0
    if k < 500:
        return 0.7 * math.sin(0.06 * (k - 340))
    return 0.8 * math.sin(0.035 * (k - 500)) + 0.3 * math.sin(0.11 * (k - 500))


def run_experiment(cfg: Config = Config()) -> dict[str, np.ndarray]:
    rng = np.random.default_rng(cfg.seed)

    # Unknown plant coefficients: y(k+1) = theta_true^T phi(k) + disturbance.
    theta_true = np.array([0.72, 0.55, 0.12, 0.0], dtype=float)
    theta_hat = np.array([0.20, 0.80, 0.00, 0.0], dtype=float)
    covariance = 30.0 * np.eye(4)

    y = np.zeros(cfg.steps + 1)
    u = np.zeros(cfg.steps)
    r = np.array([reference(k) for k in range(cfg.steps + 1)])
    theta_history = np.zeros((cfg.steps + 1, 4))
    theta_history[0] = theta_hat
    prediction_error = np.zeros(cfg.steps)

    previous_phi: np.ndarray | None = None

    for k in range(cfg.steps):
        # Update the online predictor using the transition observed at time k.
        if previous_phi is not None:
            innovation = y[k] - float(previous_phi @ theta_hat)
            denominator = cfg.forgetting + float(previous_phi @ covariance @ previous_phi)
            gain = (covariance @ previous_phi) / denominator
            theta_hat = project(theta_hat + gain * innovation)
            covariance = (
                covariance - np.outer(gain, previous_phi) @ covariance
            ) / cfg.forgetting
            covariance = 0.5 * (covariance + covariance.T)
            eigvals, eigvecs = np.linalg.eigh(covariance)
            covariance = eigvecs @ np.diag(np.clip(eigvals, 1.0e-6, 1.0e4)) @ eigvecs.T
            prediction_error[k] = innovation

        # One-step certainty-equivalent tracking control.
        a_hat, b_hat, c_hat, d_hat = theta_hat
        target = r[k + 1] - a_hat * y[k] - c_hat * y[k] ** 3 - d_hat
        inverse_b = b_hat / (b_hat * b_hat + cfg.regularization)
        probing = (0.12 * math.sin(0.31 * k) + 0.08 * math.sin(0.073 * k)) if k < 250 else 0.0
        u[k] = float(np.clip(inverse_b * target + probing, -cfg.input_limit, cfg.input_limit))

        phi = np.array([y[k], u[k], y[k] ** 3, 1.0], dtype=float)
        disturbance = 0.02 * math.sin(0.025 * k)
        y[k + 1] = float(
            phi @ theta_true + disturbance + rng.normal(0.0, cfg.noise_std)
        )
        previous_phi = phi
        theta_history[k + 1] = theta_hat

    return {
        "y": y,
        "u": u,
        "r": r,
        "theta_history": theta_history,
        "prediction_error": prediction_error,
        "theta_true": theta_true,
    }


def main() -> None:
    data = run_experiment()
    y, u, r = data["y"], data["u"], data["r"]
    theta_history = data["theta_history"]
    theta_true = data["theta_true"]

    tracking_rmse = float(np.sqrt(np.mean((r[150:] - y[150:]) ** 2)))
    print(f"Tracking RMSE after initial learning: {tracking_rmse:.6f}")
    print("Final parameter estimate:", np.array2string(theta_history[-1], precision=5))
    print("True parameter vector:  ", np.array2string(theta_true, precision=5))

    output_dir = Path(__file__).resolve().parent

    plt.figure(figsize=(10, 5))
    plt.plot(r, label="reference")
    plt.plot(y, label="output")
    plt.xlabel("sample k")
    plt.ylabel("output")
    plt.title("Learning-based adaptive tracking")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter28_Lesson2_tracking.png", dpi=160)
    plt.close()

    plt.figure(figsize=(10, 5))
    for index, label in enumerate(("a_hat", "b_hat", "c_hat", "d_hat")):
        plt.plot(theta_history[:, index], label=label)
        plt.axhline(theta_true[index], linestyle="--", linewidth=1)
    plt.xlabel("sample k")
    plt.ylabel("parameter estimate")
    plt.title("Online RLS parameter evolution")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter28_Lesson2_parameters.png", dpi=160)
    plt.close()


if __name__ == "__main__":
    main()

The C++17 implementation uses fixed-size arrays, writes a CSV file, and avoids external dependencies. In larger projects, the explicit matrix routines can be replaced by Eigen while preserving the same recursion.

Chapter28_Lesson2.cpp

// Chapter28_Lesson2.cpp
// Learning-based/data-driven adaptive control using online RLS and a
// certainty-equivalent one-step controller with projection and saturation.

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

namespace {
constexpr int N = 4;
using Vec = std::array<double, N>;
using Mat = std::array<std::array<double, N>, N>;

Vec matVec(const Mat& a, const Vec& x) {
    Vec out{};
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) out[i] += a[i][j] * x[j];
    }
    return out;
}

double dot(const Vec& a, const Vec& b) {
    double result = 0.0;
    for (int i = 0; i < N; ++i) result += a[i] * b[i];
    return result;
}

Vec project(Vec theta) {
    const Vec lo{-0.95, 0.10, -0.40, -0.50};
    const Vec hi{ 0.95, 1.50,  0.40,  0.50};
    for (int i = 0; i < N; ++i) theta[i] = std::clamp(theta[i], lo[i], hi[i]);
    return theta;
}

double reference(int k) {
    if (k < 100) return 0.0;
    if (k < 220) return 1.2;
    if (k < 340) return -1.0;
    if (k < 500) return 0.7 * std::sin(0.06 * (k - 340));
    return 0.8 * std::sin(0.035 * (k - 500)) + 0.3 * std::sin(0.11 * (k - 500));
}

void symmetrizeAndBound(Mat& p) {
    for (int i = 0; i < N; ++i) {
        for (int j = i + 1; j < N; ++j) {
            const double value = 0.5 * (p[i][j] + p[j][i]);
            p[i][j] = value;
            p[j][i] = value;
        }
        p[i][i] = std::clamp(p[i][i], 1.0e-6, 1.0e4);
    }
}
}  // namespace

int main() {
    try {
        constexpr int steps = 700;
        constexpr double forgetting = 0.998;
        constexpr double regularization = 0.05;
        constexpr double inputLimit = 3.0;
        constexpr double noiseStd = 0.005;

        const Vec thetaTrue{0.72, 0.55, 0.12, 0.0};
        Vec thetaHat{0.20, 0.80, 0.00, 0.0};
        Mat covariance{};
        for (int i = 0; i < N; ++i) covariance[i][i] = 30.0;

        std::vector<double> y(steps + 1, 0.0);
        std::vector<double> u(steps, 0.0);
        std::vector<double> r(steps + 1, 0.0);
        std::vector<Vec> thetaHistory(steps + 1);
        thetaHistory[0] = thetaHat;
        for (int k = 0; k <= steps; ++k) r[k] = reference(k);

        std::mt19937 generator(7);
        std::normal_distribution<double> noise(0.0, noiseStd);
        Vec previousPhi{};
        bool havePrevious = false;

        for (int k = 0; k < steps; ++k) {
            if (havePrevious) {
                const double innovation = y[k] - dot(previousPhi, thetaHat);
                const Vec pPhi = matVec(covariance, previousPhi);
                const double denominator = forgetting + dot(previousPhi, pPhi);
                Vec gain{};
                for (int i = 0; i < N; ++i) gain[i] = pPhi[i] / denominator;
                for (int i = 0; i < N; ++i) thetaHat[i] += gain[i] * innovation;
                thetaHat = project(thetaHat);

                Mat updated{};
                for (int i = 0; i < N; ++i) {
                    for (int j = 0; j < N; ++j) {
                        double correction = 0.0;
                        for (int m = 0; m < N; ++m) {
                            correction += gain[i] * previousPhi[m] * covariance[m][j];
                        }
                        updated[i][j] = (covariance[i][j] - correction) / forgetting;
                    }
                }
                covariance = updated;
                symmetrizeAndBound(covariance);
            }

            const double aHat = thetaHat[0];
            const double bHat = thetaHat[1];
            const double cHat = thetaHat[2];
            const double dHat = thetaHat[3];
            const double target = r[k + 1] - aHat * y[k] - cHat * std::pow(y[k], 3) - dHat;
            const double inverseB = bHat / (bHat * bHat + regularization);
            const double probing = k < 250
                ? 0.12 * std::sin(0.31 * k) + 0.08 * std::sin(0.073 * k)
                : 0.0;
            u[k] = std::clamp(inverseB * target + probing, -inputLimit, inputLimit);

            const Vec phi{y[k], u[k], std::pow(y[k], 3), 1.0};
            const double disturbance = 0.02 * std::sin(0.025 * k);
            y[k + 1] = dot(phi, thetaTrue) + disturbance + noise(generator);
            previousPhi = phi;
            havePrevious = true;
            thetaHistory[k + 1] = thetaHat;
        }

        double squaredError = 0.0;
        int count = 0;
        for (int k = 150; k <= steps; ++k) {
            const double e = r[k] - y[k];
            squaredError += e * e;
            ++count;
        }
        const double rmse = std::sqrt(squaredError / count);

        std::cout << std::fixed << std::setprecision(6)
                  << "Tracking RMSE after initial learning: " << rmse << '\n';
        std::cout << "Final parameter estimate:";
        for (double value : thetaHat) std::cout << ' ' << value;
        std::cout << "\nTrue parameter vector:  ";
        for (double value : thetaTrue) std::cout << ' ' << value;
        std::cout << '\n';

        std::ofstream csv("Chapter28_Lesson2_cpp_results.csv");
        if (!csv) throw std::runtime_error("Could not create output CSV file.");
        csv << "k,reference,output,input,a_hat,b_hat,c_hat,d_hat\n";
        for (int k = 0; k <= steps; ++k) {
            const double input = k < steps ? u[k] : 0.0;
            csv << k << ',' << r[k] << ',' << y[k] << ',' << input;
            for (double value : thetaHistory[k]) csv << ',' << value;
            csv << '\n';
        }
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

The Java implementation follows the same sample ordering as the Python and C++ programs. It uses standard arrays and exports parameter histories for plotting or further analysis.

Chapter28_Lesson2.java

// Chapter28_Lesson2.java
// Learning-based/data-driven adaptive control using online RLS and a
// certainty-equivalent one-step controller with projection and saturation.

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;
import java.util.Random;

public final class Chapter28_Lesson2 {
    private static final int N = 4;

    private Chapter28_Lesson2() {
    }

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

    private static double[] matVec(double[][] matrix, double[] vector) {
        double[] out = new double[N];
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                out[i] += matrix[i][j] * vector[j];
            }
        }
        return out;
    }

    private static double dot(double[] left, double[] right) {
        double result = 0.0;
        for (int i = 0; i < N; i++) {
            result += left[i] * right[i];
        }
        return result;
    }

    private static void project(double[] theta) {
        double[] lower = {-0.95, 0.10, -0.40, -0.50};
        double[] upper = {0.95, 1.50, 0.40, 0.50};
        for (int i = 0; i < N; i++) {
            theta[i] = clamp(theta[i], lower[i], upper[i]);
        }
    }

    private static double reference(int k) {
        if (k < 100) {
            return 0.0;
        }
        if (k < 220) {
            return 1.2;
        }
        if (k < 340) {
            return -1.0;
        }
        if (k < 500) {
            return 0.7 * Math.sin(0.06 * (k - 340));
        }
        return 0.8 * Math.sin(0.035 * (k - 500))
                + 0.3 * Math.sin(0.11 * (k - 500));
    }

    private static void symmetrizeAndBound(double[][] covariance) {
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) {
                double value = 0.5 * (covariance[i][j] + covariance[j][i]);
                covariance[i][j] = value;
                covariance[j][i] = value;
            }
            covariance[i][i] = clamp(covariance[i][i], 1.0e-6, 1.0e4);
        }
    }

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

        final int steps = 700;
        final double forgetting = 0.998;
        final double regularization = 0.05;
        final double inputLimit = 3.0;
        final double noiseStd = 0.005;

        double[] thetaTrue = {0.72, 0.55, 0.12, 0.0};
        double[] thetaHat = {0.20, 0.80, 0.00, 0.0};
        double[][] covariance = new double[N][N];
        for (int i = 0; i < N; i++) {
            covariance[i][i] = 30.0;
        }

        double[] y = new double[steps + 1];
        double[] u = new double[steps];
        double[] r = new double[steps + 1];
        double[][] thetaHistory = new double[steps + 1][N];
        System.arraycopy(thetaHat, 0, thetaHistory[0], 0, N);
        for (int k = 0; k <= steps; k++) {
            r[k] = reference(k);
        }

        Random random = new Random(7);
        double[] previousPhi = null;

        for (int k = 0; k < steps; k++) {
            if (previousPhi != null) {
                double innovation = y[k] - dot(previousPhi, thetaHat);
                double[] pPhi = matVec(covariance, previousPhi);
                double denominator = forgetting + dot(previousPhi, pPhi);
                double[] gain = new double[N];
                for (int i = 0; i < N; i++) {
                    gain[i] = pPhi[i] / denominator;
                    thetaHat[i] += gain[i] * innovation;
                }
                project(thetaHat);

                double[][] updated = new double[N][N];
                for (int i = 0; i < N; i++) {
                    for (int j = 0; j < N; j++) {
                        double correction = 0.0;
                        for (int m = 0; m < N; m++) {
                            correction += gain[i] * previousPhi[m] * covariance[m][j];
                        }
                        updated[i][j] = (covariance[i][j] - correction) / forgetting;
                    }
                }
                covariance = updated;
                symmetrizeAndBound(covariance);
            }

            double aHat = thetaHat[0];
            double bHat = thetaHat[1];
            double cHat = thetaHat[2];
            double dHat = thetaHat[3];
            double target = r[k + 1] - aHat * y[k]
                    - cHat * Math.pow(y[k], 3) - dHat;
            double inverseB = bHat / (bHat * bHat + regularization);
            double probing = k < 250
                    ? 0.12 * Math.sin(0.31 * k) + 0.08 * Math.sin(0.073 * k)
                    : 0.0;
            u[k] = clamp(inverseB * target + probing, -inputLimit, inputLimit);

            double[] phi = {y[k], u[k], Math.pow(y[k], 3), 1.0};
            double disturbance = 0.02 * Math.sin(0.025 * k);
            y[k + 1] = dot(phi, thetaTrue) + disturbance
                    + noiseStd * random.nextGaussian();
            previousPhi = phi;
            System.arraycopy(thetaHat, 0, thetaHistory[k + 1], 0, N);
        }

        double squaredError = 0.0;
        int count = 0;
        for (int k = 150; k <= steps; k++) {
            double error = r[k] - y[k];
            squaredError += error * error;
            count++;
        }
        double rmse = Math.sqrt(squaredError / count);

        System.out.printf("Tracking RMSE after initial learning: %.6f%n", rmse);
        System.out.printf("Final parameter estimate: %.6f %.6f %.6f %.6f%n",
                thetaHat[0], thetaHat[1], thetaHat[2], thetaHat[3]);
        System.out.printf("True parameter vector:   %.6f %.6f %.6f %.6f%n",
                thetaTrue[0], thetaTrue[1], thetaTrue[2], thetaTrue[3]);

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter28_Lesson2_java_results.csv"))) {
            writer.write("k,reference,output,input,a_hat,b_hat,c_hat,d_hat\n");
            for (int k = 0; k <= steps; k++) {
                double input = k < steps ? u[k] : 0.0;
                writer.write(String.format(Locale.US,
                        "%d,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                        k, r[k], y[k], input,
                        thetaHistory[k][0], thetaHistory[k][1],
                        thetaHistory[k][2], thetaHistory[k][3]));
            }
        } catch (IOException exception) {
            System.err.println("Could not write CSV: " + exception.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script runs the complete experiment and documents a direct mapping to Simulink. The RLS state consists of \( \widehat\theta_k \) and \( P_k \); these can be persistent variables in a MATLAB Function block. Unit Delay blocks preserve the previous regressor and measured transition.

Chapter28_Lesson2.m

%% Chapter28_Lesson2.m
% Learning-based/data-driven adaptive control using online RLS and a
% certainty-equivalent one-step controller with projection and saturation.
% The local function learningAdaptiveStep can also be placed in a Simulink
% MATLAB Function block, with thetaHat and covariance stored as persistent
% states and y(k), r(k+1) supplied through input ports.

clear; clc; close all;
rng(7, 'twister');

steps = 700;
forgetting = 0.998;
regularization = 0.05;
inputLimit = 3.0;
noiseStd = 0.005;

thetaTrue = [0.72; 0.55; 0.12; 0.0];
thetaHat = [0.20; 0.80; 0.00; 0.0];
covariance = 30.0 * eye(4);

y = zeros(steps + 1, 1);
u = zeros(steps, 1);
r = arrayfun(@referenceSignal, (0:steps)');
thetaHistory = zeros(steps + 1, 4);
thetaHistory(1, :) = thetaHat';
previousPhi = [];

for kIndex = 1:steps
    k = kIndex - 1;

    if ~isempty(previousPhi)
        innovation = y(kIndex) - previousPhi' * thetaHat;
        denominator = forgetting + previousPhi' * covariance * previousPhi;
        gain = covariance * previousPhi / denominator;
        thetaHat = projectParameters(thetaHat + gain * innovation);
        covariance = (covariance - gain * previousPhi' * covariance) / forgetting;
        covariance = 0.5 * (covariance + covariance');
        [vectors, values] = eig(covariance);
        boundedValues = min(max(diag(values), 1.0e-6), 1.0e4);
        covariance = vectors * diag(boundedValues) * vectors';
    end

    [u(kIndex), phi] = learningAdaptiveStep( ...
        y(kIndex), r(kIndex + 1), thetaHat, regularization, ...
        inputLimit, k);

    disturbance = 0.02 * sin(0.025 * k);
    y(kIndex + 1) = phi' * thetaTrue + disturbance + noiseStd * randn;
    previousPhi = phi;
    thetaHistory(kIndex + 1, :) = thetaHat';
end

trackingRmse = sqrt(mean((r(151:end) - y(151:end)).^2));
fprintf('Tracking RMSE after initial learning: %.6f\n', trackingRmse);
fprintf('Final parameter estimate: %.6f %.6f %.6f %.6f\n', thetaHat);
fprintf('True parameter vector:   %.6f %.6f %.6f %.6f\n', thetaTrue);

figure('Name', 'Learning-based adaptive tracking');
plot(0:steps, r, 'LineWidth', 1.2); hold on;
plot(0:steps, y, 'LineWidth', 1.1);
xlabel('sample k'); ylabel('output'); grid on;
legend('reference', 'output', 'Location', 'best');
title('Learning-based adaptive tracking');

figure('Name', 'Online parameter estimates');
plot(0:steps, thetaHistory, 'LineWidth', 1.1); hold on;
yline(thetaTrue(1), '--'); yline(thetaTrue(2), '--');
yline(thetaTrue(3), '--'); yline(thetaTrue(4), '--');
xlabel('sample k'); ylabel('parameter estimate'); grid on;
legend('a hat', 'b hat', 'c hat', 'd hat', 'Location', 'best');
title('Online RLS parameter evolution');

results = table((0:steps)', r, y, [u; 0], thetaHistory(:, 1), ...
    thetaHistory(:, 2), thetaHistory(:, 3), thetaHistory(:, 4), ...
    'VariableNames', {'k', 'reference', 'output', 'input', ...
    'a_hat', 'b_hat', 'c_hat', 'd_hat'});
writetable(results, 'Chapter28_Lesson2_matlab_results.csv');

%% Simulink mapping
% A discrete Simulink realization uses:
% 1. Unit Delay blocks for y(k-1), u(k-1), and the previous regressor.
% 2. A MATLAB Function block implementing learningAdaptiveStep and the RLS
%    update above, with persistent thetaHat and covariance variables.
% 3. A Saturation block with limits [-3, 3].
% 4. A nonlinear plant subsystem implementing
%       y(k+1) = 0.72*y(k) + 0.55*u(k) + 0.12*y(k)^3 + disturbance.
% 5. Scope blocks for reference, output, input, and parameter estimates.

function [control, phi] = learningAdaptiveStep( ...
    yCurrent, referenceNext, thetaHat, regularization, inputLimit, k)
    aHat = thetaHat(1);
    bHat = thetaHat(2);
    cHat = thetaHat(3);
    dHat = thetaHat(4);

    target = referenceNext - aHat * yCurrent - cHat * yCurrent^3 - dHat;
    inverseB = bHat / (bHat^2 + regularization);
    if k < 250
        probing = 0.12 * sin(0.31 * k) + 0.08 * sin(0.073 * k);
    else
        probing = 0.0;
    end
    control = min(max(inverseB * target + probing, -inputLimit), inputLimit);
    phi = [yCurrent; control; yCurrent^3; 1.0];
end

function theta = projectParameters(theta)
    lower = [-0.95; 0.10; -0.40; -0.50];
    upper = [ 0.95; 1.50;  0.40;  0.50];
    theta = min(max(theta, lower), upper);
end

function value = referenceSignal(k)
    if k < 100
        value = 0.0;
    elseif k < 220
        value = 1.2;
    elseif k < 340
        value = -1.0;
    elseif k < 500
        value = 0.7 * sin(0.06 * (k - 340));
    else
        value = 0.8 * sin(0.035 * (k - 500)) ...
            + 0.3 * sin(0.11 * (k - 500));
    end
end

The notebook is represented as a textual Mathematica Notebook expression. Opening the downloaded Chapter28_Lesson2.nb file creates executable input cells for the learner, controller, simulation, and plots.

Chapter28_Lesson2.nb

Notebook[{
  Cell["Chapter 28, Lesson 2: Learning-Based / Data-Driven Adaptive Control", "Title"],
  Cell["Online recursive least squares with certainty-equivalent control, projection, probing excitation, and saturation.", "Text"],
  Cell[BoxData[StringJoin[
    "ClearAll[\"Global`*\"];\nSeedRandom[7];\n\n",
    "steps = 700;\nforgetting = 0.998;\nregularization = 0.05;\n",
    "inputLimit = 3.0;\nnoiseStd = 0.005;\n\n",
    "thetaTrue = {0.72, 0.55, 0.12, 0.0};\n", "thetaHat = {0.20, 0.80, 0.00, 0.0};\n",
    "covariance = 30.0 IdentityMatrix[4];\n",
    "lower = {-0.95, 0.10, -0.40, -0.50};\n",
    "upper = {0.95, 1.50, 0.40, 0.50};\n\n",
    "project[theta_] := MapThread[Clip[#1, { #2, #3 }] &, {theta, lower, upper}];\n\n",
    "reference[k_] := Piecewise[{\n",
    "   {0.0, k < 100},\n   {1.2, k < 220},\n   {-1.0, k < 340},\n",
    "   {0.7 Sin[0.06 (k - 340)], k < 500} },\n",
    "  0.8 Sin[0.035 (k - 500)] + 0.3 Sin[0.11 (k - 500)]];\n\n",
    "y = ConstantArray[0.0, steps + 1];\nu = ConstantArray[0.0, steps];\n",
    "r = Table[reference[k], {k, 0, steps}];\n",
    "thetaHistory = ConstantArray[0.0, {steps + 1, 4}];\n",
    "thetaHistory[[1]] = thetaHat;\npreviousPhi = None;\n\n",
    "Do[\n  If[previousPhi =!= None,\n",
    "    innovation = y[[k + 1]] - previousPhi.thetaHat;\n",
    "    gain = (covariance.previousPhi)/(\n",
    "      forgetting + previousPhi.covariance.previousPhi);\n",
    "    thetaHat = project[thetaHat + gain innovation];\n",
    "    covariance = (covariance -\n",
    "       Outer[Times, gain, previousPhi].covariance)/forgetting;\n",
    "    covariance = (covariance + Transpose[covariance])/2.0;\n  ];\n\n",
    "  target = r[[k + 2]] - thetaHat[[1]] y[[k + 1]] -\n",
    "    thetaHat[[3]] y[[k + 1]]^3 - thetaHat[[4]];\n",
    "  inverseB = thetaHat[[2]]/(thetaHat[[2]]^2 + regularization);\n",
    "  probing = If[k < 250,\n",
    "    0.12 Sin[0.31 k] + 0.08 Sin[0.073 k], 0.0];\n",
    "  u[[k + 1]] = Clip[inverseB target + probing,\n",
    "    {-inputLimit, inputLimit}];\n\n",
    "  phi = {y[[k + 1]], u[[k + 1]], y[[k + 1]]^3, 1.0};\n",
    "  disturbance = 0.02 Sin[0.025 k];\n",
    "  y[[k + 2]] = phi.thetaTrue + disturbance +\n",
    "    RandomVariate[NormalDistribution[0, noiseStd]];\n",
    "  previousPhi = phi;\n  thetaHistory[[k + 2]] = thetaHat;\n",
    ", {k, 0, steps - 1}];\n\n",
    "trackingRMSE = Sqrt[Mean[(r[[151 ;;]] - y[[151 ;;]])^2]];\n",
    "Print[\"Tracking RMSE after initial learning: \", N[trackingRMSE, 8]];\n",
    "Print[\"Final parameter estimate: \", N[thetaHat, 8]];\n",
    "Print[\"True parameter vector: \", thetaTrue];\n\n",
    "trackingPlot = ListLinePlot[{r, y},\n",
    "  PlotLegends -> {\"reference\", \"output\"},\n  AxesLabel -> {\"sample k\", \"output\"},\n",
    "  PlotLabel -> \"Learning-based adaptive tracking\",\n  ImageSize -> Large];\n\n",
    "parameterPlot = ListLinePlot[Transpose[thetaHistory],\n",
    "  PlotLegends -> {\"a hat\", \"b hat\", \"c hat\", \"d hat\"},\n",
    "  AxesLabel -> {\"sample k\", \"parameter estimate\"},\n",
    "  PlotLabel -> \"Online RLS parameter evolution\",\n  ImageSize -> Large];\n\n",
    "trackingPlot\nparameterPlot\n\n",
    "results = Prepend[\n  Table[{k, r[[k + 1]], y[[k + 1]],\n",
    "    If[k < steps, u[[k + 1]], 0.0],\n",
    "    Sequence @@ thetaHistory[[k + 1]]}, {k, 0, steps}],\n",
    "  {\"k\", \"reference\", \"output\", \"input\", \"a_hat\", \"b_hat\", \"c_hat\", \"d_hat\"}];\n",
    "Export[\"Chapter28_Lesson2_mathematica_results.csv\", results];"
  ]], "Input"]
}, WindowTitle -> "Chapter28_Lesson2", StyleDefinitions -> "Default.nb"]        

The true coefficients in the supplied programs are \( [0.72,0.55,0.12,0]^{\mathsf T} \), whereas the estimator begins at \( [0.20,0.80,0,0]^{\mathsf T} \). With the fixed random seed, the tested implementations achieve a post-transient tracking RMSE of approximately \( 0.04 \). Exact estimates differ slightly across languages because their Gaussian random-number generators differ.

The parameter history is as important as the tracking plot. During step-only operation, the regressors may be correlated. The later multi-sinusoidal reference broadens the operating range and improves separation of the linear and cubic terms. Students should repeat the experiment with no probing, a constant reference, increased noise, and tighter saturation to observe how learning and tracking degrade for different reasons.

Problem 1 (RLS from Information Form): Starting from \( S_{k+1}=\lambda S_k+\phi_k\phi_k^{\mathsf T} \), derive the covariance recursion and the RLS gain.

Solution: Apply the matrix inversion lemma with \( A=\lambda S_k \), \( U=\phi_k \), and \( C=1 \):

\[ P_{k+1}=S_{k+1}^{-1}=\frac{1}{\lambda}P_k- \frac{1}{\lambda}P_k\phi_k \frac{1}{\lambda+\phi_k^{\mathsf T}P_k\phi_k} \phi_k^{\mathsf T}P_k. \]

Defining

\[ K_k=\frac{P_k\phi_k}{\lambda+ \phi_k^{\mathsf T}P_k\phi_k} \]

gives

\[ P_{k+1}=\lambda^{-1} (P_k-K_k\phi_k^{\mathsf T}P_k). \]

Combining the information-vector recursion with \( \widehat\theta_{k+1}=P_{k+1}q_{k+1} \) yields \( \widehat\theta_{k+1}=\widehat\theta_k+K_k (y_{k+1}-\phi_k^{\mathsf T}\widehat\theta_k) \) before projection.

Problem 2 (Regularization Bias): Assume no probing, no disturbance, and inactive saturation. Derive the tracking error caused by model error and regularized inversion.

Solution: The learned prediction under the applied input is

\[ \widehat y_{k+1}=\widehat a_k y_k+ \widehat b_k\frac{\widehat b_k}{\widehat b_k^2+\rho}v_k+ \widehat c_k y_k^3+\widehat d_k. \]

Since the non-input terms plus \( v_k \) equal the reference,

\[ \widehat y_{k+1}-r_{k+1} =-\frac{\rho}{\widehat b_k^2+\rho}v_k. \]

Adding the prediction mismatch gives

\[ e_{k+1}=\widetilde\theta_k^{\mathsf T}\phi_k- \frac{\rho}{\widehat b_k^2+\rho}v_k. \]

Thus regularization prevents singular inversion at the price of a known bias that vanishes as \( \rho \) decreases, provided the learned gain remains safely separated from zero.

Problem 3 (Failure of Excitation): Suppose regulation reaches the constant equilibrium \( y_k=\bar y \) and \( u_k=\bar u \). Determine the rank of the accumulated regressor matrix.

Solution: Every regressor is identical:

\[ \phi_k=\bar\phi=[\bar y,\bar u,\bar y^3,1]^{\mathsf T}. \]

Therefore

\[ \sum_{i=k}^{k+N-1}\phi_i\phi_i^{\mathsf T} =N\bar\phi\bar\phi^{\mathsf T}, \]

which has rank one when \( \bar\phi\neq0 \). Four independent parameter directions cannot be identified. Tracking at that operating point may still be excellent because only the scalar combination \( \theta_*^{\mathsf T}\bar\phi \) matters locally.

Problem 4 (Projection and Control Direction): Explain why the projection interval for \( \widehat b_k \) in the programs is \( [0.10,1.50] \) rather than an interval containing zero.

Solution: The inverse controller assumes a known positive control direction. If the estimate crosses zero, the control command may become arbitrarily large or reverse sign, creating positive feedback. The projection encodes prior knowledge \( b_*>0 \) and supplies a lower bound used in the tracking-error inequality. If the sign is genuinely unknown, this controller is inappropriate; Chapter 16 methods based on Nussbaum-type gains or another sign-robust architecture are required.

Problem 5 (Neural Approximation Bound): From \( \dot V\leq-q\|e\|^2+2c\bar\epsilon\|e\| \), where \( q>0 \), derive an ultimate bound.

Solution: Factor the right-hand side:

\[ \dot V\leq-\|e\|\left(q\|e\|-2c\bar\epsilon\right). \]

Hence \( \dot V<0 \) whenever \( \|e\|>2c\bar\epsilon/q \). Standard ultimate-boundedness arguments imply convergence to a compact set whose error radius is no larger than a constant proportional to \( 2c\bar\epsilon/q \). A richer approximator may reduce \( \bar\epsilon \), but it may also increase training complexity and uncertainty outside the approximation domain.

Problem 6 (Hankel-Matrix Dimensions): A scalar input record has length \( T \). What are the dimensions of \( H_L(u^d) \), and what condition is necessary for it to have full row rank?

Solution: The matrix contains \( L \) shifted rows and \( T-L+1 \) columns, so

\[ H_L(u^d)\in\mathbb R^{L\times(T-L+1)}. \]

Full row rank requires at least as many columns as rows, \( T-L+1\geq L \), and the samples must be sufficiently diverse. The dimension inequality is necessary but not sufficient; repeated or highly correlated inputs can still produce rank deficiency.

Learning-based adaptive control extends classical adaptation by using online optimization, statistical learning, nonlinear function approximation, or trajectory representations. RLS is an exact online optimizer for a weighted least-squares loss and can drive a certainty-equivalent nonlinear controller. Its performance depends on excitation, estimator conditioning, regularization, projection, and saturation. Neural approximators can represent broader uncertainties but normally provide ultimate rather than exact convergence when approximation error remains. Behavioral methods show that sufficiently rich measured trajectories can serve as implicit models. Across all approaches, prediction quality, closed-loop stability, and constraint safety must be analyzed separately.