After completing this lesson, a student should be able to:

• distinguish linearity in unknown parameters from linearity in measured signals;
• construct a causal regressor vector from sampled input-output data;
• stack scalar regression equations into a matrix regression;
• derive the normal equations and state the rank condition for a unique estimate;
• explain why excitation, scaling, noise correlation, and delay alignment matter;
• implement and validate the regression model in five computational environments.

In an indirect adaptive controller, the identifier receives measured plant data, estimates physical or input-output parameters, and passes those estimates to a controller-design map. The identifier cannot operate until the plant equation has been written in a suitable parametric form.

flowchart TD
  U["Measured input u(k)"] --> R["Construct causal regressor phi(k)"]
  Y["Measured output y(k)"] --> R
  R --> P["Predict output using parameter estimate"]
  P --> E["Prediction error"]
  E --> ID["Parameter estimator in later lessons"]
  ID --> C["Controller update in an indirect adaptive loop"]
        

Let \( y(t) \in \mathbb{R} \) be a measured scalar signal, \( \boldsymbol{\phi}(t) \in \mathbb{R}^{p} \) a known regressor vector, \( \boldsymbol{\theta} \in \mathbb{R}^{p} \) an unknown constant parameter vector, and \( v(t) \) a disturbance, measurement error, or model mismatch. The regression equation is

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

Expanding the inner product gives

\[ y(t)=\sum_{i=1}^{p}\phi_i(t)\theta_i+v(t). \]

The model is linear in the unknown parameters because every \( \theta_i \) appears to the first power and parameters do not multiply one another. The regressor components may themselves be nonlinear functions of measured variables. For example,

\[ y(t)=\theta_1 x(t)+\theta_2 x^2(t)+\theta_3\sin(u(t))+v(t) \]

is still linear in the parameter vector

\[ \boldsymbol{\theta}=\begin{bmatrix}\theta_1&\theta_2&\theta_3\end{bmatrix}^{\mathsf{T}},\qquad \boldsymbol{\phi}(t)=\begin{bmatrix}x(t)&x^2(t)&\sin(u(t))\end{bmatrix}^{\mathsf{T}}. \]

By contrast, \( y=\theta_1 e^{\theta_2 x} \) is not linear in the unknown parameter vector and cannot be used directly by ordinary RLS.

Given a current parameter estimate \( \widehat{\boldsymbol{\theta}}(t) \), define the predicted output as

\[ \widehat{y}(t)=\boldsymbol{\phi}^{\mathsf{T}}(t)\widehat{\boldsymbol{\theta}}(t). \]

The prediction error is

\[ \varepsilon(t)=y(t)-\widehat{y}(t). \]

Define the parameter error using the convention

\[ \widetilde{\boldsymbol{\theta}}(t)=\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}}(t). \]

Substitution gives the fundamental error decomposition

\[ \varepsilon(t)=\boldsymbol{\phi}^{\mathsf{T}}(t)\widetilde{\boldsymbol{\theta}}(t)+v(t). \]

This identity separates error caused by incorrect parameters from error caused by noise or structural mismatch. Even when the parameters are exact, \( \varepsilon(t)=v(t) \); therefore, zero instantaneous prediction error is generally not a sensible requirement in noisy systems.

Consider the sampled first-order input-output model

\[ y(k)=a\,y(k-1)+b\,u(k-1)+v(k). \]

Choose

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

Then the plant is exactly in regression form:

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

The one-sample delay is essential. At sample \( k \), the values \( y(k-1) \) and \( u(k-1) \) are already available, so the regressor is causal. Using a future signal such as \( u(k+1) \) would create an unrealizable online estimator.

A more general finite-order model may be written as

\[ y(k)=\sum_{i=1}^{n_a}a_i y(k-i)+\sum_{j=1}^{n_b}b_j u(k-j)+v(k). \]

with

\[ \boldsymbol{\theta}=\begin{bmatrix}a_1&\cdots&a_{n_a}&b_1&\cdots&b_{n_b}\end{bmatrix}^{\mathsf{T}}, \]

\[ \boldsymbol{\phi}(k)=\begin{bmatrix}y(k-1)&\cdots&y(k-n_a)&u(k-1)&\cdots&u(k-n_b)\end{bmatrix}^{\mathsf{T}}. \]

Some identification texts place minus signs before past-output terms. Either convention is valid, but the parameter definition and regressor definition must use the same sign convention throughout the estimator and controller.

Suppose samples are collected at instants \( t_1,\ldots,t_N \). Stack the scalar equations:

\[ \mathbf{Y}_N= \begin{bmatrix}y(t_1)\\y(t_2)\\\vdots\\y(t_N)\end{bmatrix},\qquad \mathbf{V}_N= \begin{bmatrix}v(t_1)\\v(t_2)\\\vdots\\v(t_N)\end{bmatrix}. \]

The regression matrix is formed by placing each transposed regressor in one row:

\[ \mathbf{\Phi}_N= \begin{bmatrix} \boldsymbol{\phi}^{\mathsf{T}}(t_1)\\ \boldsymbol{\phi}^{\mathsf{T}}(t_2)\\ \vdots\\ \boldsymbol{\phi}^{\mathsf{T}}(t_N) \end{bmatrix}\in\mathbb{R}^{N\times p}. \]

The complete data equation is

\[ \mathbf{Y}_N=\mathbf{\Phi}_N\boldsymbol{\theta}+\mathbf{V}_N. \]

Dimension checks are a reliable way to detect implementation mistakes: \( \mathbf{\Phi}_N \) has one row per observation and one column per unknown parameter, while \( \boldsymbol{\theta} \) is a column vector.

Define the finite-data information or Gram matrix

\[ \mathbf{G}_N=\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N =\sum_{i=1}^{N}\boldsymbol{\phi}(t_i)\boldsymbol{\phi}^{\mathsf{T}}(t_i). \]

For any vector \( \mathbf{z}\in\mathbb{R}^{p} \),

\[ \mathbf{z}^{\mathsf{T}}\mathbf{G}_N\mathbf{z} =\sum_{i=1}^{N}\left(\boldsymbol{\phi}^{\mathsf{T}}(t_i)\mathbf{z}\right)^2. \]

Hence \( \mathbf{G}_N \) is symmetric positive semidefinite. It is positive definite exactly when no nonzero direction \( \mathbf{z} \) is orthogonal to every collected regressor. The equivalent uniqueness conditions are

\[ \operatorname{rank}(\mathbf{\Phi}_N)=p \quad\Longleftrightarrow\quad \mathbf{G}_N\succ 0 \quad\Longleftrightarrow\quad \ker(\mathbf{\Phi}_N)=\{\mathbf{0}\}. \]

If the rank is deficient, two or more parameter vectors generate the same predictions on the collected data. This is a data-based identifiability failure, not merely a numerical inconvenience.

flowchart TD
  D["Collected regressors"] --> M["Build data matrix Phi"]
  M --> Q["Check rank and singular values"]
  Q -->|"full column rank"| U["Unique finite-data parameter fit"]
  Q -->|"rank deficient"| N["Unobservable parameter direction"]
  U --> K["Inspect conditioning before online estimation"]
  N --> X["Change input, model order, or measured signals"]
        

Persistent excitation was introduced in Chapter 10. In this chapter it supplies the long-window version of the same idea: the regressor must keep exciting all parameter directions rather than becoming confined to a lower- dimensional subspace.

Before implementing a recursive estimator, it is useful to validate the regression offline. Consider the quadratic criterion

\[ J_N(\boldsymbol{\vartheta})= \frac{1}{2}\left\|\mathbf{Y}_N-\mathbf{\Phi}_N\boldsymbol{\vartheta}\right\|_2^2. \]

Expanding the criterion yields

\[ J_N(\boldsymbol{\vartheta})= \frac{1}{2}\mathbf{Y}_N^{\mathsf{T}}\mathbf{Y}_N -\boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{Y}_N +\frac{1}{2}\boldsymbol{\vartheta}^{\mathsf{T}} \mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N\boldsymbol{\vartheta}. \]

The gradient and Hessian are

\[ \nabla_{\boldsymbol{\vartheta}}J_N =-\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{Y}_N +\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N\boldsymbol{\vartheta}, \]

\[ \nabla_{\boldsymbol{\vartheta}}^2 J_N =\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N=\mathbf{G}_N. \]

Setting the gradient to zero gives the normal equations

\[ \mathbf{G}_N\widehat{\boldsymbol{\theta}}_N =\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{Y}_N. \]

If \( \mathbf{G}_N\succ 0 \), the unique minimizer is

\[ \widehat{\boldsymbol{\theta}}_N =\left(\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N\right)^{-1} \mathbf{\Phi}_N^{\mathsf{T}}\mathbf{Y}_N. \]

Numerically, one should solve the least-squares system using QR or singular- value decomposition rather than explicitly forming the inverse. The formula is used here for analysis. Lesson 2 will show how the same information is updated recursively without resolving the entire batch problem after every sample.

Let the residual vector at the least-squares solution be

\[ \mathbf{e}_N=\mathbf{Y}_N-\mathbf{\Phi}_N\widehat{\boldsymbol{\theta}}_N. \]

The normal equations imply

\[ \mathbf{\Phi}_N^{\mathsf{T}}\mathbf{e}_N=\mathbf{0}. \]

Thus the residual is orthogonal to every column of the regression matrix. Geometrically, \( \mathbf{\Phi}_N\widehat{\boldsymbol{\theta}}_N \) is the orthogonal projection of \( \mathbf{Y}_N \) onto the column space of \( \mathbf{\Phi}_N \).

If the data satisfy

\[ \mathbf{Y}_N=\mathbf{\Phi}_N\boldsymbol{\theta}+\mathbf{V}_N, \]

then the estimation error is

\[ \widehat{\boldsymbol{\theta}}_N-\boldsymbol{\theta} =\left(\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N\right)^{-1} \mathbf{\Phi}_N^{\mathsf{T}}\mathbf{V}_N. \]

Conditional unbiasedness follows when the noise has zero conditional mean relative to the regressor data:

\[ \mathbb{E}\!\left[\mathbf{V}_N\mid\mathbf{\Phi}_N\right]=\mathbf{0}. \]

Under conditionally white noise with covariance

\[ \operatorname{Cov}\!\left(\mathbf{V}_N\mid\mathbf{\Phi}_N\right) =\sigma_v^2\mathbf{I}, \]

the conditional parameter covariance is

\[ \operatorname{Cov}\!\left(\widehat{\boldsymbol{\theta}}_N\mid\mathbf{\Phi}_N\right) =\sigma_v^2\left(\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{\Phi}_N\right)^{-1}. \]

Large information eigenvalues reduce parameter uncertainty. A tiny minimum eigenvalue indicates that at least one parameter combination is weakly informed and will be sensitive to noise.

9.1 Causality and timing

Every component of \( \boldsymbol{\phi}(k) \) must be known when the estimate at sample \( k \) is updated. Input-output records must be aligned with physical delays and acquisition timestamps.

9.2 Parameter constancy over the estimation window

The baseline model assumes a fixed \( \boldsymbol{\theta} \). Slowly varying parameters motivate forgetting mechanisms, but those mechanisms do not repair an incorrect regressor structure.

9.3 Offset and bias terms

A constant unknown offset is represented by appending a component equal to one:

\[ y(k)=\begin{bmatrix}\boldsymbol{\phi}^{\mathsf{T}}(k)&1\end{bmatrix} \begin{bmatrix}\boldsymbol{\theta}\\d\end{bmatrix}+v(k). \]

9.4 Scaling and conditioning

If one regressor component is many orders of magnitude larger than another, the Gram matrix can be ill-conditioned even when it has full rank. Physical normalization or consistent engineering units improve numerical reliability.

9.5 Closed-loop correlation and bias

In feedback operation, past outputs may contain noise that also influences future inputs. Therefore the regressor can become statistically correlated with the disturbance. The simple least-squares unbiasedness condition may fail, even though the algebraic regression equation appears correct.

9.6 Overparameterization

Adding unnecessary delays or basis functions increases dimension and can create nearly dependent columns. A more complicated model is not necessarily more useful for adaptive control; it may reduce estimator speed and amplify noise sensitivity.

Suppose the true sampled plant is

\[ y(k)=0.78\,y(k-1)+0.42\,u(k-1)+v(k). \]

The unknown parameter vector and measured regressor are

\[ \boldsymbol{\theta}=\begin{bmatrix}0.78\\0.42\end{bmatrix},\qquad \boldsymbol{\phi}(k)=\begin{bmatrix}y(k-1)\\u(k-1)\end{bmatrix}. \]

For four usable samples, assume

\[ \mathbf{\Phi}_4= \begin{bmatrix} 0.00&1.00\\ 0.42&-0.50\\ 0.12&0.80\\ 0.43&-0.20 \end{bmatrix},\qquad \mathbf{Y}_4= \begin{bmatrix}0.42\\0.118\\0.4296\\0.2514\end{bmatrix}. \]

Then

\[ \mathbf{\Phi}_4^{\mathsf{T}}\mathbf{\Phi}_4= \begin{bmatrix}0.3757&-0.2000\\-0.2000&1.9300\end{bmatrix}. \]

Its determinant is

\[ \det\!\left(\mathbf{\Phi}_4^{\mathsf{T}}\mathbf{\Phi}_4\right) =0.685101. \]

Since the determinant is nonzero, the two columns are linearly independent and the finite-data fit is unique. With noise-free values generated by the stated plant, solving the normal equations recovers the true coefficients.

This NumPy implementation generates a dynamic data record, constructs \( \mathbf{\Phi} \), solves a batch least-squares validation problem using a stable library routine, and checks the orthogonality condition. Required library: numpy.

Chapter12_Lesson1.py

"""Chapter 12, Lesson 1: linear-in-the-parameters regression model.

Generates data from y(k) = a*y(k-1) + b*u(k-1) + v(k), constructs
phi(k) = [y(k-1), u(k-1)]^T, and solves the batch least-squares
problem only to validate that the chosen regression is identifiable.
RLS recursion is intentionally deferred to Lesson 2.
"""

from __future__ import annotations

import numpy as np


def build_regression(y: np.ndarray, u: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Return Phi and Y for k = 1, ..., N-1."""
    if y.ndim != 1 or u.ndim != 1 or y.size != u.size:
        raise ValueError("y and u must be one-dimensional arrays of equal length")
    if y.size < 3:
        raise ValueError("at least three samples are required")

    phi = np.column_stack((y[:-1], u[:-1]))
    target = y[1:]
    return phi, target


def main() -> None:
    rng = np.random.default_rng(12)
    n_samples = 300
    theta_true = np.array([0.78, 0.42])
    noise_std = 0.03

    # A persistently varying input produces a well-conditioned data matrix.
    u = rng.normal(0.0, 1.0, n_samples)
    y = np.zeros(n_samples)
    for k in range(1, n_samples):
        regressor = np.array([y[k - 1], u[k - 1]])
        y[k] = regressor @ theta_true + rng.normal(0.0, noise_std)

    phi, target = build_regression(y, u)
    gram = phi.T @ phi
    theta_hat, residuals, rank, singular_values = np.linalg.lstsq(
        phi, target, rcond=None
    )

    prediction = phi @ theta_hat
    residual = target - prediction

    print("True theta:     ", theta_true)
    print("Estimated theta:", theta_hat)
    print("rank(Phi):      ", rank)
    print("singular values:", singular_values)
    print("cond(Phi^T Phi):", np.linalg.cond(gram))
    print("residual RMSE:  ", np.sqrt(np.mean(residual**2)))
    print("normal-equation residual:", np.linalg.norm(phi.T @ residual))
    if residuals.size:
        print("least-squares SSE:", residuals[0])


if __name__ == "__main__":
    main()

The C++ version uses only the standard library. Because this lesson's example has two parameters, the normal equations are solved analytically as a two-by-two system. For larger models, use a tested linear-algebra library such as Eigen and solve by QR or SVD.

Chapter12_Lesson1.cpp

// Chapter 12, Lesson 1: linear-in-the-parameters regression model.
// Standard-library-only demonstration for a two-parameter dynamic regression.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <vector>

struct Estimate {
    double theta1{};
    double theta2{};
    double determinant{};
};

Estimate solveTwoParameterLeastSquares(const std::vector<double>& y,
                                       const std::vector<double>& u) {
    if (y.size() != u.size() || y.size() < 3) {
        throw std::invalid_argument("y and u must have equal length >= 3");
    }

    double s11 = 0.0;
    double s12 = 0.0;
    double s22 = 0.0;
    double q1 = 0.0;
    double q2 = 0.0;

    for (std::size_t k = 1; k < y.size(); ++k) {
        const double phi1 = y[k - 1];
        const double phi2 = u[k - 1];
        const double target = y[k];

        s11 += phi1 * phi1;
        s12 += phi1 * phi2;
        s22 += phi2 * phi2;
        q1 += phi1 * target;
        q2 += phi2 * target;
    }

    const double det = s11 * s22 - s12 * s12;
    if (std::abs(det) < 1.0e-12) {
        throw std::runtime_error("Phi^T Phi is singular: parameters are not identifiable");
    }

    return {
        (s22 * q1 - s12 * q2) / det,
        (-s12 * q1 + s11 * q2) / det,
        det
    };
}

int main() {
    constexpr std::size_t nSamples = 300;
    constexpr double aTrue = 0.78;
    constexpr double bTrue = 0.42;
    constexpr double noiseStd = 0.03;

    std::mt19937 generator(12);
    std::normal_distribution<double> inputDistribution(0.0, 1.0);
    std::normal_distribution<double> noiseDistribution(0.0, noiseStd);

    std::vector<double> u(nSamples, 0.0);
    std::vector<double> y(nSamples, 0.0);
    for (double& sample : u) {
        sample = inputDistribution(generator);
    }

    for (std::size_t k = 1; k < nSamples; ++k) {
        y[k] = aTrue * y[k - 1] + bTrue * u[k - 1]
             + noiseDistribution(generator);
    }

    try {
        const Estimate estimate = solveTwoParameterLeastSquares(y, u);

        double squaredError = 0.0;
        double normal1 = 0.0;
        double normal2 = 0.0;
        for (std::size_t k = 1; k < nSamples; ++k) {
            const double prediction = estimate.theta1 * y[k - 1]
                                    + estimate.theta2 * u[k - 1];
            const double residual = y[k] - prediction;
            squaredError += residual * residual;
            normal1 += y[k - 1] * residual;
            normal2 += u[k - 1] * residual;
        }

        const double rmse = std::sqrt(squaredError / static_cast<double>(nSamples - 1));
        const double normalResidual = std::sqrt(normal1 * normal1 + normal2 * normal2);

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "True theta:      [" << aTrue << ", " << bTrue << "]\n";
        std::cout << "Estimated theta: [" << estimate.theta1 << ", "
                  << estimate.theta2 << "]\n";
        std::cout << "det(Phi^T Phi):  " << estimate.determinant << '\n';
        std::cout << "residual RMSE:   " << rmse << '\n';
        std::cout << "normal residual: " << normalResidual << '\n';
    } catch (const std::exception& exception) {
        std::cerr << "Estimation failed: " << exception.what() << '\n';
        return 1;
    }

    return 0;
}

The Java implementation also uses only the standard library and performs the same dimensional, singularity, residual, and normal-equation checks. For larger regression matrices, suitable libraries include EJML and Apache Commons Math, preferably with QR or singular-value decomposition rather than an explicitly formed matrix inverse.

Chapter12_Lesson1.java

// Chapter 12, Lesson 1: linear-in-the-parameters regression model.
// Standard-library-only demonstration for a two-parameter dynamic regression.

import java.util.Locale;
import java.util.Random;

public final class Chapter12_Lesson1 {
    private Chapter12_Lesson1() {
    }

    private static final class Estimate {
        final double theta1;
        final double theta2;
        final double determinant;

        Estimate(double theta1, double theta2, double determinant) {
            this.theta1 = theta1;
            this.theta2 = theta2;
            this.determinant = determinant;
        }
    }

    private static Estimate solveTwoParameterLeastSquares(double[] y, double[] u) {
        if (y.length != u.length || y.length < 3) {
            throw new IllegalArgumentException("y and u must have equal length >= 3");
        }

        double s11 = 0.0;
        double s12 = 0.0;
        double s22 = 0.0;
        double q1 = 0.0;
        double q2 = 0.0;

        for (int k = 1; k < y.length; k++) {
            double phi1 = y[k - 1];
            double phi2 = u[k - 1];
            double target = y[k];

            s11 += phi1 * phi1;
            s12 += phi1 * phi2;
            s22 += phi2 * phi2;
            q1 += phi1 * target;
            q2 += phi2 * target;
        }

        double determinant = s11 * s22 - s12 * s12;
        if (Math.abs(determinant) < 1.0e-12) {
            throw new IllegalStateException(
                "Phi^T Phi is singular: parameters are not identifiable"
            );
        }

        double theta1 = (s22 * q1 - s12 * q2) / determinant;
        double theta2 = (-s12 * q1 + s11 * q2) / determinant;
        return new Estimate(theta1, theta2, determinant);
    }

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

        final int nSamples = 300;
        final double aTrue = 0.78;
        final double bTrue = 0.42;
        final double noiseStd = 0.03;

        Random random = new Random(12L);
        double[] u = new double[nSamples];
        double[] y = new double[nSamples];

        for (int k = 0; k < nSamples; k++) {
            u[k] = random.nextGaussian();
        }
        for (int k = 1; k < nSamples; k++) {
            y[k] = aTrue * y[k - 1] + bTrue * u[k - 1]
                 + noiseStd * random.nextGaussian();
        }

        Estimate estimate = solveTwoParameterLeastSquares(y, u);

        double squaredError = 0.0;
        double normal1 = 0.0;
        double normal2 = 0.0;
        for (int k = 1; k < nSamples; k++) {
            double prediction = estimate.theta1 * y[k - 1]
                              + estimate.theta2 * u[k - 1];
            double residual = y[k] - prediction;
            squaredError += residual * residual;
            normal1 += y[k - 1] * residual;
            normal2 += u[k - 1] * residual;
        }

        double rmse = Math.sqrt(squaredError / (nSamples - 1));
        double normalResidual = Math.hypot(normal1, normal2);

        System.out.printf("True theta:      [%.6f, %.6f]%n", aTrue, bTrue);
        System.out.printf(
            "Estimated theta: [%.6f, %.6f]%n",
            estimate.theta1, estimate.theta2
        );
        System.out.printf("det(Phi^T Phi):  %.6f%n", estimate.determinant);
        System.out.printf("residual RMSE:   %.6f%n", rmse);
        System.out.printf("normal residual: %.6e%n", normalResidual);
    }
}

MATLAB constructs the stacked regression and uses the backslash operator, which selects an appropriate numerical factorization. When Simulink is licensed, the script also creates a model implementing \( \widehat{y}=\widehat{\theta}_1\phi_1+\widehat{\theta}_2\phi_2 \) with two Gain blocks and a Sum block.

Chapter12_Lesson1.m

%% Chapter12_Lesson1.m
% Linear-in-the-parameters model for online parameter estimation.
% The batch least-squares result is used only to validate the regression
% structure. Recursive least squares is introduced in Lesson 2.

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

N = 300;
thetaTrue = [0.78; 0.42];
noiseStd = 0.03;

u = randn(N, 1);
y = zeros(N, 1);
for k = 2:N
    phiK = [y(k - 1); u(k - 1)];
    y(k) = phiK.' * thetaTrue + noiseStd * randn();
end

Phi = [y(1:N - 1), u(1:N - 1)];
Y = y(2:N);
G = Phi.' * Phi;

if rcond(G) < 1.0e-12
    error('Phi''*Phi is numerically singular; parameters are not identifiable.');
end

thetaHat = Phi \ Y;
residual = Y - Phi * thetaHat;

fprintf('True theta:      [%.6f, %.6f]\n', thetaTrue(1), thetaTrue(2));
fprintf('Estimated theta: [%.6f, %.6f]\n', thetaHat(1), thetaHat(2));
fprintf('rank(Phi):       %d\n', rank(Phi));
fprintf('cond(Phi''*Phi): %.6f\n', cond(G));
fprintf('residual RMSE:   %.6f\n', sqrt(mean(residual.^2)));
fprintf('normal residual: %.6e\n', norm(Phi.' * residual));

figure('Name', 'Chapter 12 Lesson 1');
plot(2:N, Y, 'DisplayName', 'Measured output');
hold on;
plot(2:N, Phi * thetaHat, '--', 'DisplayName', 'Linear-regression prediction');
grid on;
xlabel('Sample k');
ylabel('Output');
legend('Location', 'best');
title('Measured output and fitted parametric model');

%% Optional Simulink construction
% This creates a static regression block y_hat = theta1*phi1 + theta2*phi2.
% It requires Simulink. The model is saved but not opened automatically.
if license('test', 'Simulink')
    modelName = 'Chapter12_Lesson1_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if exist([modelName '.slx'], 'file')
        delete([modelName '.slx']);
    end

    new_system(modelName);
    add_block('simulink/Sources/In1', [modelName '/phi1'], ...
        'Position', [40 55 70 75]);
    add_block('simulink/Sources/In1', [modelName '/phi2'], ...
        'Position', [40 125 70 145], 'Port', '2');
    add_block('simulink/Math Operations/Gain', [modelName '/theta1'], ...
        'Gain', num2str(thetaHat(1), 17), 'Position', [120 45 190 85]);
    add_block('simulink/Math Operations/Gain', [modelName '/theta2'], ...
        'Gain', num2str(thetaHat(2), 17), 'Position', [120 115 190 155]);
    add_block('simulink/Math Operations/Sum', [modelName '/sum'], ...
        'Inputs', '++', 'Position', [240 75 270 125]);
    add_block('simulink/Sinks/Out1', [modelName '/y_hat'], ...
        'Position', [330 90 360 110]);

    add_line(modelName, 'phi1/1', 'theta1/1');
    add_line(modelName, 'phi2/1', 'theta2/1');
    add_line(modelName, 'theta1/1', 'sum/1');
    add_line(modelName, 'theta2/1', 'sum/2');
    add_line(modelName, 'sum/1', 'y_hat/1');

    set_param(modelName, 'Solver', 'FixedStepDiscrete', 'StopTime', '10');
    save_system(modelName);
    close_system(modelName);
    fprintf('Created Simulink model: %s.slx\n', modelName);
else
    fprintf('Simulink is unavailable; skipped model generation.\n');
end

The notebook uses LeastSquares, MatrixRank, and Eigenvalues to validate the same parametric model and inspect its information matrix.

Chapter12_Lesson1.nb

Notebook[{
 Cell["Chapter 12, Lesson 1: Parametric Model with Linear Regression", "Title"],
 Cell["This notebook constructs phi(k) and theta for y(k)=phi(k)^T theta+v(k). Batch least squares is used only to validate the model; recursive updates are deferred to Lesson 2.", "Text"],
 Cell[BoxData["ClearAll[\"Global`*\"]\nSeedRandom[12];\nnSamples = 300;\nthetaTrue = {0.78, 0.42};\nnoiseStd = 0.03;\nu = RandomVariate[NormalDistribution[0, 1], nSamples];\ny = ConstantArray[0., nSamples];\nDo[\n  phiK = {y[[k - 1]], u[[k - 1]]};\n  y[[k]] = phiK.thetaTrue + RandomVariate[NormalDistribution[0, noiseStd]],\n  {k, 2, nSamples}\n];"], "Input"],
 Cell[BoxData["phi = Table[{y[[k - 1]], u[[k - 1]]}, {k, 2, nSamples}];\ntarget = y[[2 ;;]];\ngram = Transpose[phi].phi;\nIf[MatrixRank[phi] < Length[thetaTrue],\n  Print[\"Regression matrix is rank deficient.\"]; Abort[]\n];\nthetaHat = LeastSquares[phi, target];\nprediction = phi.thetaHat;\nresidual = target - prediction;"], "Input"],
 Cell[BoxData["Print[\"True theta:      \", thetaTrue];\nPrint[\"Estimated theta: \", thetaHat];\nPrint[\"rank(Phi):       \", MatrixRank[phi]];\nPrint[\"eigenvalues(Phi^T Phi): \", Eigenvalues[gram]];\nPrint[\"residual RMSE:   \", Sqrt[Mean[residual^2]]];\nPrint[\"normal residual: \", Norm[Transpose[phi].residual]];"], "Input"],
 Cell[BoxData["ListLinePlot[\n  {target, prediction},\n  PlotLegends -> {\"Measured output\", \"Regression prediction\"},\n  AxesLabel -> {\"sample\", \"output\"},\n  PlotLabel -> \"Linear-in-the-parameters model validation\",\n  ImageSize -> Large\n]"], "Input"]
}, WindowTitle -> "Chapter12_Lesson1"]

Problem 1 (Regression Construction): Consider \( y(k)+a_1y(k-1)+a_2y(k-2)=b_1u(k-1)+v(k) \). Construct a regression of the form \( y(k)=\boldsymbol{\phi}^{\mathsf{T}}(k)\boldsymbol{\theta}+v(k) \).

Solution: Isolate the measured output:

\[ y(k)=-a_1y(k-1)-a_2y(k-2)+b_1u(k-1)+v(k). \]

A consistent choice is

\[ \boldsymbol{\theta}=\begin{bmatrix}a_1\\a_2\\b_1\end{bmatrix},\qquad \boldsymbol{\phi}(k)=\begin{bmatrix}-y(k-1)\\-y(k-2)\\u(k-1)\end{bmatrix}. \]

The regressor uses only past measurements and is therefore causal. An equally valid convention absorbs the minus signs into the parameter definitions.

Problem 2 (Rank Deficiency): Let every collected regressor be \( \boldsymbol{\phi}(k)=\begin{bmatrix}1&2\end{bmatrix}^{\mathsf{T}} \). Can two parameters be uniquely identified?

Solution: The data matrix has identical rows:

\[ \mathbf{\Phi}_N= \begin{bmatrix}1&2\\1&2\\\vdots&\vdots\\1&2\end{bmatrix}. \]

Its rank is one, not two. Only the combination \( \theta_1+2\theta_2 \) affects the predicted output. For any scalar \( c \), the parameter pairs

\[ \begin{bmatrix}\theta_1\\\theta_2\end{bmatrix} \quad\text{and}\quad \begin{bmatrix}\theta_1-2c\\\theta_2+c\end{bmatrix} \]

produce identical predictions. The experiment must generate non-collinear regressors before a unique estimate is possible.

Problem 3 (Proof of the Orthogonality Condition): Show that a minimizer of the least-squares criterion has a residual orthogonal to every column of \( \mathbf{\Phi}_N \).

Solution: The criterion is

\[ J_N(\boldsymbol{\vartheta})= \frac{1}{2}\left(\mathbf{Y}_N-\mathbf{\Phi}_N\boldsymbol{\vartheta}\right)^{\mathsf{T}} \left(\mathbf{Y}_N-\mathbf{\Phi}_N\boldsymbol{\vartheta}\right). \]

At a differentiable unconstrained minimum,

\[ \mathbf{0}=\nabla J_N =-\mathbf{\Phi}_N^{\mathsf{T}} \left(\mathbf{Y}_N-\mathbf{\Phi}_N\widehat{\boldsymbol{\theta}}_N\right). \]

Defining the residual vector gives

\[ \mathbf{\Phi}_N^{\mathsf{T}}\mathbf{e}_N=\mathbf{0}. \]

The \( j \)th entry is the inner product between the residual and the \( j \)th regression column, proving orthogonality to the complete column space.

Problem 4 (Bias from Regressor-Noise Correlation): Starting from the least-squares error expression, explain why \( \mathbb{E}[\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{V}_N]\neq\mathbf{0} \) can create biased estimates.

Solution: The parameter estimation error is

\[ \widehat{\boldsymbol{\theta}}_N-\boldsymbol{\theta} =\mathbf{G}_N^{-1}\mathbf{\Phi}_N^{\mathsf{T}}\mathbf{V}_N. \]

If the regressor-noise cross term has nonzero mean, the expected right-hand side generally does not vanish. In a closed loop, this can occur because measurement noise influences the controller and therefore future input values, while delayed noisy outputs are also included in the regressor. Increasing the sample count alone does not automatically remove structural correlation.

Problem 5 (Conditioning and Rescaling): Suppose \( \phi_1 \) is measured in millivolts and has magnitude near \( 10^3 \), while \( \phi_2 \) has magnitude near \( 10^{-2} \). Explain how to rescale the regression without changing its predictions.

Solution: Introduce a nonsingular diagonal scaling matrix

\[ \mathbf{S}=\operatorname{diag}(s_1,s_2),\qquad \overline{\boldsymbol{\phi}}=\mathbf{S}^{-1}\boldsymbol{\phi},\qquad \overline{\boldsymbol{\theta}}=\mathbf{S}\boldsymbol{\theta}. \]

Then

\[ \overline{\boldsymbol{\phi}}^{\mathsf{T}}\overline{\boldsymbol{\theta}} =\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\theta}. \]

Choose \( s_1 \) and \( s_2 \) so that the scaled regressor components have comparable typical magnitudes. After estimation, recover the physical parameters through \( \boldsymbol{\theta}=\mathbf{S}^{-1}\overline{\boldsymbol{\theta}} \).

A control-oriented parameter estimator begins with the model \( y(t)=\boldsymbol{\phi}^{\mathsf{T}}(t)\boldsymbol{\theta}+v(t) \). The regressor must contain known, causal, correctly aligned signals, while the unknown quantities are collected in a constant parameter vector. Stacking data produces \( \mathbf{Y}=\mathbf{\Phi}\boldsymbol{\theta}+\mathbf{V} \). Full column rank of \( \mathbf{\Phi} \) is necessary for a unique finite-data fit; singular values and the information matrix reveal weak parameter directions. Batch least squares provides the normal equations, residual orthogonality, and a useful structural test before recursive operation. The next lesson converts this batch information into RLS update equations and introduces forgetting factors for online tracking.