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

• derive an exact zero-order-hold discrete model from a continuous-time linear plant;
• write causal SISO difference equations using backward-shift polynomials;
• convert a plant equation into the linear regression form used by online estimators;
• construct and verify an equivalent discrete state-space realization;
• check delay, stability, model order, and data-rank assumptions before adaptation.

In adaptive control, the plant model has two simultaneous roles. It is a prediction model that explains measured data and a control-oriented model from which controller coefficients will later be computed or directly adjusted. A model that predicts well but hides input delay, unstable zeros, or the sign of the leading input gain may still be unsuitable for adaptive controller design.

flowchart TD
  CT["Continuous plant and sampler"] --> DT["Discrete-time dynamics"]
  DT --> IO["Polynomial input-output model"]
  DT --> SS["State-space model"]
  IO --> REG["Linear regression: \noutput = regressor times parameters"]
  SS --> PROP["Poles, controllability, observability"]
  REG --> EST["Online estimation in later lessons"]
  PROP --> CTRL["Adaptive controller construction in later lessons"]
        

Consider the continuous-time LTI plant \( \dot{\mathbf{x} }(t)=\mathbf{A}_c\mathbf{x}(t)+\mathbf{B}_c\mathbf{u}(t) \) with sampled output \( \mathbf{y}(k)=\mathbf{C}_c\mathbf{x}(kT_s)+\mathbf{D}_c\mathbf{u}(kT_s) \). A zero-order hold keeps the input constant on every sampling interval:

\[ \mathbf{u}(t)=\mathbf{u}(k),\qquad kT_s\leq t < (k+1)T_s. \]

Exact discretization theorem. The sampled state satisfies

\[ \mathbf{x}(k+1)=\mathbf{A}_d\mathbf{x}(k)+\mathbf{B}_d\mathbf{u}(k), \qquad \mathbf{y}(k)=\mathbf{C}_d\mathbf{x}(k)+\mathbf{D}_d\mathbf{u}(k), \]

\[ \mathbf{A}_d=e^{\mathbf{A}_cT_s},\qquad \mathbf{B}_d=\int_0^{T_s}e^{\mathbf{A}_c\tau}\mathbf{B}_c\,d\tau, \qquad \mathbf{C}_d=\mathbf{C}_c,\qquad \mathbf{D}_d=\mathbf{D}_c. \]

Proof. Variation of constants over one interval gives

\[ \mathbf{x}((k+1)T_s)=e^{\mathbf{A}_cT_s}\mathbf{x}(kT_s) +\int_{kT_s}^{(k+1)T_s}e^{\mathbf{A}_c((k+1)T_s-\sigma)} \mathbf{B}_c\mathbf{u}(k)\,d\sigma. \]

Set \( \tau=(k+1)T_s-\sigma \). Because the held input is constant, it leaves the integral, producing exactly \( \mathbf{A}_d \) and \( \mathbf{B}_d \) above. If \( \mathbf{A}_c \) is nonsingular, integration of the matrix exponential also yields

\[ \mathbf{B}_d=\mathbf{A}_c^{-1} \left(\mathbf{A}_d-\mathbf{I}\right)\mathbf{B}_c. \]

For the scalar plant \( \dot y(t)=-\alpha y(t)+\beta u(t) \), with \( \alpha > 0 \), the exact sampled model is

\[ y(k+1)=e^{-\alpha T_s}y(k) +\frac{\beta}{\alpha}\left(1-e^{-\alpha T_s}\right)u(k). \]

Thus sampling changes the numerical parameterization. The unknown continuous parameters \( \alpha,\beta \) do not enter this exact discrete model linearly, whereas the discrete coefficients can be treated as independent unknown parameters and estimated linearly.

Let \( q^{-1} \) denote the backward-shift operator: \( q^{-1}y(k)=y(k-1) \). A causal SISO plant model is

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

\[ A(q^{-1})=1+a_1q^{-1}+\cdots+a_{n_a}q^{-n_a}, \qquad B(q^{-1})=b_0+b_1q^{-1}+\cdots+b_{n_b}q^{-n_b}. \]

Here \(d\geq 1\) is the known integer input delay and \(v(k)\) collects disturbances, measurement noise, and unmodeled effects. Expanding the equation gives

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

The normalization of the leading coefficient of \(A\) to one removes a scale ambiguity. The quantity \(b_0\) is the leading or high-frequency input gain in this convention. Its nonzero value and, in many adaptive designs, its known sign are structural assumptions.

The corresponding rational transfer description for the deterministic part is

\[ G(q^{-1})=\frac{q^{-d}B(q^{-1})}{A(q^{-1})}. \]

No Z-transform derivation is required here: the expression is shorthand for the same causal recurrence. The roots of \(z^{n_a}+a_1z^{n_a-1}+\cdots+a_{n_a}\) are the plant poles. The deterministic model is asymptotically stable when every pole obeys \( |z_i|<1 \).

The polynomial model is linear in its unknown coefficients. Define

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

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

Then the plant equation becomes

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

This equation is the bridge from plant dynamics to adaptation. The regressor contains measured past inputs and outputs; the parameter vector contains constant but unknown plant coefficients. An online estimator can therefore update \(\hat{\boldsymbol{\theta} }(k)\) without nonlinear optimization.

Stacked-data form. For samples \(k=k_0,\ldots,N\),

\[ \mathbf{Y}_N=\mathbf{\Phi}_N\boldsymbol{\theta}+\mathbf{V}_N, \qquad \mathbf{\Phi}_N= \begin{bmatrix}\boldsymbol{\phi}^{T}(k_0)\\ \boldsymbol{\phi}^{T}(k_0+1)\\ \vdots\\ \boldsymbol{\phi}^{T}(N)\end{bmatrix}. \]

Uniqueness proposition. In the noise-free case, the parameter vector is unique if

\[ \mathbf{\Phi}_N^{T}\mathbf{\Phi}_N > 0. \]

Proof. Suppose two vectors generate the same output data.

\[ \mathbf{\Phi}_N\boldsymbol{\theta}_1= \mathbf{\Phi}_N\boldsymbol{\theta}_2. \]

With \(\boldsymbol{\delta}=\boldsymbol{\theta}_1- \boldsymbol{\theta}_2\), one has \(\mathbf{\Phi}_N\boldsymbol{\delta}=\mathbf{0}\), hence

\[ 0=\|\mathbf{\Phi}_N\boldsymbol{\delta}\|_2^2 =\boldsymbol{\delta}^{T}\mathbf{\Phi}_N^{T}\mathbf{\Phi}_N \boldsymbol{\delta}. \]

Positive definiteness implies \(\boldsymbol{\delta}=\mathbf{0}\). Therefore \(\boldsymbol{\theta}_1=\boldsymbol{\theta}_2\). This is the finite-data counterpart of the persistent-excitation ideas developed in Chapter 10.

The general discrete state-space plant is

\[ \mathbf{x}(k+1)=\mathbf{A}_d\mathbf{x}(k)+\mathbf{B}_d u(k),\qquad y(k)=\mathbf{C}_d\mathbf{x}(k)+D_d u(k). \]

State-space form is preferable when states are measured, the plant is multivariable, or structural properties such as controllability and observability must be explicit. Input-output form is preferable for output-feedback self-tuning regulators because it directly produces a measurable regressor.

For the strictly proper one-step-delay transfer model

\[ G(z)=\frac{b_0z^{n-1}+b_1z^{n-2}+\cdots+b_{n-1} } {z^n+a_1z^{n-1}+\cdots+a_n}, \]

a controllable companion realization is

\[ \begin{aligned} \mathbf{A}_d &= \begin{bmatrix} -a_1&-a_2&\cdots&-a_{n-1}&-a_n\\ 1&0&\cdots&0&0\\ 0&1&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&1&0 \end{bmatrix}, \qquad \mathbf{B}_d = \begin{bmatrix} 1\\0\\\vdots\\0 \end{bmatrix}, \\[8pt] \mathbf{C}_d &= \begin{bmatrix} b_0&b_1&\cdots&b_{n-1} \end{bmatrix}, \qquad D_d=0. \end{aligned} \]

In the displayed pattern, the last row contains a one in the next-to-last column and zero elsewhere. For additional pure delay \(d>1\), append input-delay states.

Equivalence proof. The first column of

\[ (z\mathbf{I}-\mathbf{A}_d)^{-1}\mathbf{B}_d \]

has denominator equal to the characteristic polynomial

\[ \det(z\mathbf{I}-\mathbf{A}_d)=z^n+a_1z^{n-1}+\cdots+a_n. \]

Its numerator vector is \(\begin{bmatrix}z^{n-1}&z^{n-2}&\cdots&1\end{bmatrix}^{T}\). Left multiplication by \(\mathbf{C}_d\) gives the stated numerator polynomial, so \(\mathbf{C}_d(z\mathbf{I}-\mathbf{A}_d)^{-1}\mathbf{B}_d=G(z)\).

Adaptive control needs an explicit distinction between plant dynamics and stochastic uncertainty. With white innovation \(e(k)\), common model structures include:

ARX:

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

ARMAX:

\[ A(q^{-1})y(k)=q^{-d}B(q^{-1})u(k)+C(q^{-1})e(k). \]

Output-error:

\[ y(k)=\frac{q^{-d}B(q^{-1})}{F(q^{-1})}u(k)+e(k). \]

Box-Jenkins:

\[ y(k)=\frac{q^{-d}B(q^{-1})}{F(q^{-1})}u(k) +\frac{C(q^{-1})}{D(q^{-1})}e(k). \]

ARX is especially convenient for introductory adaptive control because it is immediately linear in the plant parameters. Its limitation is that the same polynomial \(A\) shapes both deterministic dynamics and noise. Output-error and Box-Jenkins structures can describe noise more faithfully, but their prediction-error parameterization is generally more involved. The correct choice is the simplest model that preserves the dynamics needed for controller design and leaves a sufficiently uncorrelated prediction residual.

Causality. A controller can only use data available at sample \(k\). Therefore a strictly proper plant model normally has \(d\geq1\).

Known delay. An incorrect delay shifts input columns of the regressor. Even with abundant data, the estimated polynomial then represents the wrong control channel.

Stability. For

\[ A(q^{-1})y(k)=w(k), \]

if all roots of the associated polynomial in \(z\) satisfy \(|z_i|<1\), the impulse response is absolutely summable. Consequently, every bounded input \(w(k)\) produces a bounded zero-state output. This is the discrete BIBO-stability condition.

Model order. Under-modeling moves omitted dynamics into \(v(k)\) and can bias estimates; over-modeling creates weakly excited directions and poor numerical conditioning.

Minimum phase. Some later certainty-equivalent control laws require the zeros of \(B\) to lie inside the unit disk. This is not required merely to simulate or estimate the plant, but it must be recorded before selecting the controller architecture.

Sampling period. A very large \(T_s\) hides fast dynamics and can create aliasing; a very small \(T_s\) makes adjacent regressors nearly collinear and can worsen finite-precision estimation. Sampling is therefore part of the model design.

flowchart TD
  A["Choose sampling period and signal units"] --> B["Determine input-output delay"]
  B --> C["Select state-space or polynomial structure"]
  C --> D["Choose orders na and nb"]
  D --> E["Form regressor from available past data"]
  E --> F["Check boundedness, rank, and conditioning"]
  F --> G["Estimate or verify coefficients"]
  G --> H["Inspect poles, zeros, residuals, and delay"]
  H --> I{"Adequate for control?"}
  I -->|"no"| C
  I -->|"yes"| J["Freeze assumptions for adaptive-law design"]
        

A plant model should not be accepted only because its simulation curve appears plausible. The delay, regressor rank, pole locations, residual correlation, and leading input gain must all be checked. These checks prevent the adaptive law from compensating for a modeling error that it cannot structurally remove.

All implementations use the stable plant

\[ y(k)=1.50y(k-1)-0.56y(k-2)+0.12u(k-1)+0.05u(k-2). \]

Its normalized polynomial parameters are

\[ A(q^{-1})=1-1.50q^{-1}+0.56q^{-2},\qquad B(q^{-1})=0.12+0.05q^{-1},\qquad d=1. \]

The parameter vector and regressor are

\[ \boldsymbol{\theta}= \begin{bmatrix}-1.50&0.56&0.12&0.05\end{bmatrix}^{T},\qquad \boldsymbol{\phi}(k)= \begin{bmatrix}-y(k-1)&-y(k-2)&u(k-1)&u(k-2)\end{bmatrix}^{T} \]

An equivalent realization is

\[ \mathbf{A}_d=\begin{bmatrix}1.50&-0.56\\1&0\end{bmatrix},\quad \mathbf{B}_d=\begin{bmatrix}1\\0\end{bmatrix},\quad \mathbf{C}_d=\begin{bmatrix}0.12&0.05\end{bmatrix},\quad D_d=0. \]

The poles are \(0.8\) and \(0.7\), so the deterministic plant is stable. Every program below performs three tests: exact agreement of input-output and state-space simulations, full-rank regression, and recovery of the four coefficients from noise-free data.

NumPy supplies arrays, least squares, eigenvalues, and polynomial roots. SciPy's signal module and the third-party python-control package are useful for larger discrete state-space and transfer-function studies, but this implementation keeps the plant recurrence explicit.

Chapter14_Lesson1.py

"""Chapter14_Lesson1.py

Discrete-time plant-model laboratory for Adaptive Control, Chapter 14, Lesson 1.
The script verifies equivalence between an ARX difference equation and a
controllable-canonical state-space realization, then recovers the linearly
parameterized model from noise-free data.

Required library: NumPy
Optional visualization: Matplotlib
"""

from __future__ import annotations

import numpy as np


def input_signal(n_samples: int) -> np.ndarray:
    """Create a deterministic, sufficiently rich bounded input."""
    k = np.arange(n_samples, dtype=float)
    square = np.where((k.astype(int) % 40) < 20, 0.3, -0.3)
    return 0.7 * np.sin(0.11 * k) + 0.4 * np.sin(0.037 * k) + square


def simulate_input_output(u: np.ndarray) -> np.ndarray:
    """Simulate A(q^-1)y(k) = q^-1 B(q^-1)u(k)."""
    y = np.zeros_like(u, dtype=float)
    for k in range(1, len(u)):
        y_k_minus_2 = y[k - 2] if k >= 2 else 0.0
        u_k_minus_2 = u[k - 2] if k >= 2 else 0.0
        y[k] = (
            1.50 * y[k - 1]
            - 0.56 * y_k_minus_2
            + 0.12 * u[k - 1]
            + 0.05 * u_k_minus_2
        )
    return y


def simulate_state_space(u: np.ndarray) -> np.ndarray:
    """Simulate the equivalent controllable-canonical realization."""
    a_d = np.array([[1.50, -0.56], [1.0, 0.0]], dtype=float)
    b_d = np.array([1.0, 0.0], dtype=float)
    c_d = np.array([0.12, 0.05], dtype=float)

    x = np.zeros(2, dtype=float)
    y = np.zeros_like(u, dtype=float)
    for k, uk in enumerate(u):
        y[k] = c_d @ x
        x = a_d @ x + b_d * uk
    return y


def build_regression(u: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Build Y = Phi theta for theta = [a1, a2, b0, b1]^T."""
    phi = []
    target = []
    for k in range(2, len(y)):
        phi.append([-y[k - 1], -y[k - 2], u[k - 1], u[k - 2]])
        target.append(y[k])
    return np.asarray(phi, dtype=float), np.asarray(target, dtype=float)


def main() -> None:
    n_samples = 240
    u = input_signal(n_samples)
    y_io = simulate_input_output(u)
    y_ss = simulate_state_space(u)

    phi, target = build_regression(u, y_io)
    theta_hat, residuals, rank, singular_values = np.linalg.lstsq(
        phi, target, rcond=None
    )

    true_theta = np.array([-1.50, 0.56, 0.12, 0.05])
    gram = phi.T @ phi
    gram_eigenvalues = np.linalg.eigvalsh(gram)
    poles = np.roots([1.0, -1.50, 0.56])

    print("Maximum input-output/state-space mismatch:", np.max(np.abs(y_io - y_ss)))
    print("True theta:     ", true_theta)
    print("Estimated theta:", theta_hat)
    print("Regression rank:", rank)
    print("Singular values:", singular_values)
    print("Minimum eigenvalue of Phi^T Phi:", gram_eigenvalues.min())
    print("Plant poles:", poles)
    print("Stable:", bool(np.all(np.abs(poles) < 1.0)))
    print("Residual sum of squares:", float(residuals[0]) if residuals.size else 0.0)

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(u, label="u(k)")
        plt.plot(y_io, label="y(k)")
        plt.xlabel("Sample k")
        plt.ylabel("Amplitude")
        plt.title("Chapter 14 Lesson 1: discrete-time plant data")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.savefig("Chapter14_Lesson1_response.png", dpi=160)
        print("Saved optional plot: Chapter14_Lesson1_response.png")
    except ImportError:
        print("Matplotlib is not installed; numerical verification is complete.")


if __name__ == "__main__":
    main()

The example uses only the C++17 standard library and implements the four-parameter normal-equation solve from scratch. For larger systems, Eigen is the customary header-only linear algebra library; QR or SVD should replace normal equations when conditioning is poor.

Chapter14_Lesson1.cpp

// Chapter14_Lesson1.cpp
// Discrete-time plant model, state-space equivalence, and batch regression.
// Build: g++ -std=c++17 -O2 Chapter14_Lesson1.cpp -o Chapter14_Lesson1

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

using Vector4 = std::array<double, 4>;
using Matrix4 = std::array<std::array<double, 4>, 4>;

std::vector<double> inputSignal(std::size_t n) {
    std::vector<double> u(n, 0.0);
    for (std::size_t k = 0; k < n; ++k) {
        const double square = (k % 40 < 20) ? 0.3 : -0.3;
        u[k] = 0.7 * std::sin(0.11 * static_cast<double>(k))
             + 0.4 * std::sin(0.037 * static_cast<double>(k)) + square;
    }
    return u;
}

std::vector<double> simulateInputOutput(const std::vector<double>& u) {
    std::vector<double> y(u.size(), 0.0);
    for (std::size_t k = 1; k < u.size(); ++k) {
        const double yKMinus2 = (k >= 2) ? y[k - 2] : 0.0;
        const double uKMinus2 = (k >= 2) ? u[k - 2] : 0.0;
        y[k] = 1.50 * y[k - 1] - 0.56 * yKMinus2
             + 0.12 * u[k - 1] + 0.05 * uKMinus2;
    }
    return y;
}

std::vector<double> simulateStateSpace(const std::vector<double>& u) {
    std::vector<double> y(u.size(), 0.0);
    std::array<double, 2> x{0.0, 0.0};
    for (std::size_t k = 0; k < u.size(); ++k) {
        y[k] = 0.12 * x[0] + 0.05 * x[1];
        const std::array<double, 2> next{
            1.50 * x[0] - 0.56 * x[1] + u[k],
            x[0]
        };
        x = next;
    }
    return y;
}

Vector4 solveLinearSystem(Matrix4 a, Vector4 b) {
    for (std::size_t pivot = 0; pivot < 4; ++pivot) {
        std::size_t best = pivot;
        for (std::size_t row = pivot + 1; row < 4; ++row) {
            if (std::abs(a[row][pivot]) > std::abs(a[best][pivot])) {
                best = row;
            }
        }
        if (std::abs(a[best][pivot]) < 1e-12) {
            throw std::runtime_error("Normal matrix is singular or ill-conditioned.");
        }
        std::swap(a[pivot], a[best]);
        std::swap(b[pivot], b[best]);

        const double diagonal = a[pivot][pivot];
        for (std::size_t col = pivot; col < 4; ++col) {
            a[pivot][col] /= diagonal;
        }
        b[pivot] /= diagonal;

        for (std::size_t row = 0; row < 4; ++row) {
            if (row == pivot) {
                continue;
            }
            const double factor = a[row][pivot];
            for (std::size_t col = pivot; col < 4; ++col) {
                a[row][col] -= factor * a[pivot][col];
            }
            b[row] -= factor * b[pivot];
        }
    }
    return b;
}

int main() {
    constexpr std::size_t nSamples = 240;
    const auto u = inputSignal(nSamples);
    const auto yIO = simulateInputOutput(u);
    const auto ySS = simulateStateSpace(u);

    Matrix4 gram{};
    Vector4 rhs{};
    double maxMismatch = 0.0;

    for (std::size_t k = 0; k < nSamples; ++k) {
        maxMismatch = std::max(maxMismatch, std::abs(yIO[k] - ySS[k]));
    }

    for (std::size_t k = 2; k < nSamples; ++k) {
        const Vector4 phi{-yIO[k - 1], -yIO[k - 2], u[k - 1], u[k - 2]};
        for (std::size_t i = 0; i < 4; ++i) {
            rhs[i] += phi[i] * yIO[k];
            for (std::size_t j = 0; j < 4; ++j) {
                gram[i][j] += phi[i] * phi[j];
            }
        }
    }

    const Vector4 thetaHat = solveLinearSystem(gram, rhs);
    const double discriminant = 1.50 * 1.50 - 4.0 * 0.56;
    const double pole1 = (1.50 + std::sqrt(discriminant)) / 2.0;
    const double pole2 = (1.50 - std::sqrt(discriminant)) / 2.0;

    std::cout << std::setprecision(12);
    std::cout << "Maximum input-output/state-space mismatch: " << maxMismatch << '\n';
    std::cout << "Estimated theta: [";
    for (std::size_t i = 0; i < thetaHat.size(); ++i) {
        std::cout << thetaHat[i] << (i + 1 == thetaHat.size() ? "]\n" : ", ");
    }
    std::cout << "Plant poles: " << pole1 << ", " << pole2 << '\n';
    std::cout << "Stable: " << std::boolalpha
              << (std::abs(pole1) < 1.0 && std::abs(pole2) < 1.0) << '\n';
    return 0;
}

The Java version also uses no external dependency. Apache Commons Math or EJML is recommended for production-quality QR, SVD, eigenvalue, and state-space computations.

Chapter14_Lesson1.java

// Chapter14_Lesson1.java
// Discrete-time plant model, state-space equivalence, and batch regression.
// Build and run: javac Chapter14_Lesson1.java && java Chapter14_Lesson1

import java.util.Arrays;

public final class Chapter14_Lesson1 {
    private static final int PARAMETER_COUNT = 4;

    private Chapter14_Lesson1() {
    }

    private static double[] inputSignal(int n) {
        double[] u = new double[n];
        for (int k = 0; k < n; k++) {
            double square = (k % 40 < 20) ? 0.3 : -0.3;
            u[k] = 0.7 * Math.sin(0.11 * k) + 0.4 * Math.sin(0.037 * k) + square;
        }
        return u;
    }

    private static double[] simulateInputOutput(double[] u) {
        double[] y = new double[u.length];
        for (int k = 1; k < u.length; k++) {
            double yKMinus2 = (k >= 2) ? y[k - 2] : 0.0;
            double uKMinus2 = (k >= 2) ? u[k - 2] : 0.0;
            y[k] = 1.50 * y[k - 1] - 0.56 * yKMinus2
                    + 0.12 * u[k - 1] + 0.05 * uKMinus2;
        }
        return y;
    }

    private static double[] simulateStateSpace(double[] u) {
        double[] y = new double[u.length];
        double[] x = {0.0, 0.0};
        for (int k = 0; k < u.length; k++) {
            y[k] = 0.12 * x[0] + 0.05 * x[1];
            double[] next = {
                    1.50 * x[0] - 0.56 * x[1] + u[k],
                    x[0]
            };
            x = next;
        }
        return y;
    }

    private static double[] solveLinearSystem(double[][] aInput, double[] bInput) {
        double[][] a = new double[PARAMETER_COUNT][PARAMETER_COUNT];
        double[] b = Arrays.copyOf(bInput, PARAMETER_COUNT);
        for (int i = 0; i < PARAMETER_COUNT; i++) {
            a[i] = Arrays.copyOf(aInput[i], PARAMETER_COUNT);
        }

        for (int pivot = 0; pivot < PARAMETER_COUNT; pivot++) {
            int best = pivot;
            for (int row = pivot + 1; row < PARAMETER_COUNT; row++) {
                if (Math.abs(a[row][pivot]) > Math.abs(a[best][pivot])) {
                    best = row;
                }
            }
            if (Math.abs(a[best][pivot]) < 1e-12) {
                throw new IllegalStateException("Normal matrix is singular or ill-conditioned.");
            }

            double[] rowTemp = a[pivot];
            a[pivot] = a[best];
            a[best] = rowTemp;
            double scalarTemp = b[pivot];
            b[pivot] = b[best];
            b[best] = scalarTemp;

            double diagonal = a[pivot][pivot];
            for (int col = pivot; col < PARAMETER_COUNT; col++) {
                a[pivot][col] /= diagonal;
            }
            b[pivot] /= diagonal;

            for (int row = 0; row < PARAMETER_COUNT; row++) {
                if (row == pivot) {
                    continue;
                }
                double factor = a[row][pivot];
                for (int col = pivot; col < PARAMETER_COUNT; col++) {
                    a[row][col] -= factor * a[pivot][col];
                }
                b[row] -= factor * b[pivot];
            }
        }
        return b;
    }

    public static void main(String[] args) {
        int nSamples = 240;
        double[] u = inputSignal(nSamples);
        double[] yIO = simulateInputOutput(u);
        double[] ySS = simulateStateSpace(u);

        double[][] gram = new double[PARAMETER_COUNT][PARAMETER_COUNT];
        double[] rhs = new double[PARAMETER_COUNT];
        double maxMismatch = 0.0;

        for (int k = 0; k < nSamples; k++) {
            maxMismatch = Math.max(maxMismatch, Math.abs(yIO[k] - ySS[k]));
        }

        for (int k = 2; k < nSamples; k++) {
            double[] phi = {-yIO[k - 1], -yIO[k - 2], u[k - 1], u[k - 2]};
            for (int i = 0; i < PARAMETER_COUNT; i++) {
                rhs[i] += phi[i] * yIO[k];
                for (int j = 0; j < PARAMETER_COUNT; j++) {
                    gram[i][j] += phi[i] * phi[j];
                }
            }
        }

        double[] thetaHat = solveLinearSystem(gram, rhs);
        double discriminant = 1.50 * 1.50 - 4.0 * 0.56;
        double pole1 = (1.50 + Math.sqrt(discriminant)) / 2.0;
        double pole2 = (1.50 - Math.sqrt(discriminant)) / 2.0;

        System.out.printf("Maximum input-output/state-space mismatch: %.12e%n", maxMismatch);
        System.out.println("Estimated theta: " + Arrays.toString(thetaHat));
        System.out.printf("Plant poles: %.12f, %.12f%n", pole1, pole2);
        System.out.println("Stable: " + (Math.abs(pole1) < 1.0 && Math.abs(pole2) < 1.0));
    }
}

MATLAB's matrix operators make regression and pole checks direct. Control System Toolbox provides ss, tf, c2d, and dlsim. The script below performs the core calculations from first principles and, when Simulink is licensed, creates a discrete transfer-function verification model programmatically.

Chapter14_Lesson1.m

% Chapter14_Lesson1.m
% Discrete-time plant model, state-space equivalence, batch regression,
% and optional programmatic construction of a Simulink verification model.

clear; clc; close all;

N = 240;
k = (0:N-1)';
squarePart = 0.3 * ones(N,1);
squarePart(mod(k,40) >= 20) = -0.3;
u = 0.7*sin(0.11*k) + 0.4*sin(0.037*k) + squarePart;

% A(q^-1)y(k) = q^-1 B(q^-1)u(k)
% A(q^-1) = 1 - 1.50 q^-1 + 0.56 q^-2
% B(q^-1) = 0.12 + 0.05 q^-1
yIO = zeros(N,1);
yIO(2) = 0.12*u(1); % zero prehistory: y(-1)=0 and u(-1)=0
for index = 3:N
    yIO(index) = 1.50*yIO(index-1) - 0.56*yIO(index-2) ...
               + 0.12*u(index-1) + 0.05*u(index-2);
end

% Equivalent controllable-canonical state-space realization.
Ad = [1.50 -0.56; 1.0 0.0];
Bd = [1.0; 0.0];
Cd = [0.12 0.05];
Dd = 0.0;

x = zeros(2,1);
ySS = zeros(N,1);
for index = 1:N
    ySS(index) = Cd*x + Dd*u(index);
    x = Ad*x + Bd*u(index);
end

% Batch regression Y = Phi*theta.
Phi = [-yIO(2:N-1), -yIO(1:N-2), u(2:N-1), u(1:N-2)];
Y = yIO(3:N);
thetaHat = Phi \ Y;
trueTheta = [-1.50; 0.56; 0.12; 0.05];
plantPoles = roots([1 -1.50 0.56]);

fprintf('Maximum input-output/state-space mismatch: %.12e\n', max(abs(yIO-ySS)));
fprintf('True theta:      [% .8f % .8f % .8f % .8f]^T\n', trueTheta);
fprintf('Estimated theta: [% .8f % .8f % .8f % .8f]^T\n', thetaHat);
fprintf('rank(Phi) = %d\n', rank(Phi));
fprintf('minimum eigenvalue of Phi''*Phi = %.8e\n', min(eig(Phi'*Phi)));
fprintf('plant poles = %.8f, %.8f\n', plantPoles(1), plantPoles(2));
fprintf('stable = %d\n', all(abs(plantPoles) < 1));

figure('Name','Chapter14 Lesson1');
plot(k,u,'DisplayName','u(k)'); hold on;
plot(k,yIO,'DisplayName','y(k)');
grid on; xlabel('Sample k'); ylabel('Amplitude');
title('Discrete-time plant input and output'); legend('Location','best');

% Optional Simulink construction. The Discrete Transfer Fcn block uses the
% explicit one-sample delay through the leading zero in the numerator.
if license('test','Simulink')
    modelName = 'Chapter14_Lesson1_Model';
    if bdIsLoaded(modelName)
        close_system(modelName,0);
    end
    if exist([modelName '.slx'],'file')
        delete([modelName '.slx']);
    end

    new_system(modelName);
    open_system(modelName);
    set_param(modelName,'StopTime',num2str(N-1),'Solver','FixedStepDiscrete', ...
        'FixedStep','1');

    inputData = timeseries(u,k);
    assignin('base','inputData',inputData);

    add_block('simulink/Sources/From Workspace',[modelName '/Input'], ...
        'VariableName','inputData','Position',[40 70 150 100]);
    add_block('simulink/Discrete/Discrete Transfer Fcn',[modelName '/Plant'], ...
        'Numerator','[0 0.12 0.05]', ...
        'Denominator','[1 -1.50 0.56]', ...
        'SampleTime','1','Position',[220 55 390 115]);
    add_block('simulink/Sinks/To Workspace',[modelName '/Output'], ...
        'VariableName','simulinkOutput','SaveFormat','Timeseries', ...
        'Position',[470 70 580 100]);
    add_line(modelName,'Input/1','Plant/1');
    add_line(modelName,'Plant/1','Output/1');
    save_system(modelName);

    simulationResult = sim(modelName);
    ySimulink = simulationResult.simulinkOutput.Data;
    fprintf('Maximum MATLAB/Simulink mismatch: %.12e\n', ...
        max(abs(yIO-ySimulink(1:N))));
else
    fprintf('Simulink license unavailable; MATLAB verification completed.\n');
end

Mathematica provides exact symbolic polynomial manipulation together with numerical LeastSquares, MatrixRank, and Eigenvalues. The downloadable file is a notebook expression containing executable input cells.

Chapter14_Lesson1.nb

Notebook[{
  Cell[
    "Chapter 14, Lesson 1: Discrete-Time Plant Models for Adaptive Control",
    "Title"
  ],

  Cell[
    "This notebook verifies an ARX model, an equivalent state-space realization, and linear regression parameterization.",
    "Text"
  ],

  Cell[BoxData[
    "ClearAll[\"Global`*\"]"
  ], "Input"],

  Cell[BoxData[
    "n = 240;\nk = Range[0, n - 1];"
  ], "Input"],

  Cell[BoxData[
    "u = Table[\n  0.7 Sin[0.11 i] + 0.4 Sin[0.037 i] +\n    If[Mod[i, 40] < 20, 0.3, -0.3],\n  {i, 0, n - 1}\n];"
  ], "Input"],

  Cell[BoxData[
    "yIO = ConstantArray[0., n];\nyIO[[2]] = 0.12 u[[1]];\nDo[\n  yIO[[i]] = 1.50 yIO[[i - 1]] - 0.56 yIO[[i - 2]] +\n    0.12 u[[i - 1]] + 0.05 u[[i - 2]],\n  {i, 3, n}\n];"
  ], "Input"],

  Cell[BoxData[
    "Ad = { {1.50, -0.56}, {1.0, 0.0} };\nBd = {1.0, 0.0};\nCd = {0.12, 0.05};"
  ], "Input"],

  Cell[BoxData[
    "x = {0., 0.};\nySS = ConstantArray[0., n];\nDo[\n  ySS[[i]] = Cd.x;\n  x = Ad.x + Bd u[[i]],\n  {i, 1, n}\n];"
  ], "Input"],

  Cell[BoxData[
    "phi = Table[\n  {-yIO[[i - 1]], -yIO[[i - 2]], u[[i - 1]], u[[i - 2]]},\n  {i, 3, n}\n];\ntarget = yIO[[3 ;; n]];"
  ], "Input"],

  Cell[BoxData[
    "thetaHat = LeastSquares[phi, target];\ntrueTheta = {-1.50, 0.56, 0.12, 0.05};"
  ], "Input"],

  Cell[BoxData[
    "poles = z /. Solve[z^2 - 1.50 z + 0.56 == 0, z];"
  ], "Input"],

  Cell[BoxData[
    "Print[\"Maximum input-output/state-space mismatch: \", Max[Abs[yIO - ySS]]];\nPrint[\"True theta: \", trueTheta];\nPrint[\"Estimated theta: \", thetaHat];\nPrint[\"Rank of Phi: \", MatrixRank[phi]];\nPrint[\"Minimum eigenvalue of Transpose[Phi].Phi: \", Min[Eigenvalues[Transpose[phi].phi]]];\nPrint[\"Plant poles: \", poles];\nPrint[\"Stable: \", And @@ Thread[Abs[poles] < 1]];"
  ], "Input"],

  Cell[BoxData[
    "ListLinePlot[\n  {u, yIO},\n  PlotLegends -> {\"u(k)\", \"y(k)\"},\n  GridLines -> Automatic,\n  AxesLabel -> {\"k\", \"amplitude\"},\n  PlotLabel -> \"Discrete-time plant data\"\n]"
  ], "Input"]
},
  WindowTitle -> "Chapter14_Lesson1.nb",
  StyleDefinitions -> "Default.nb"
]        

In exact arithmetic, the two simulations coincide because they are realizations of the same rational operator. The batch estimate also equals the true coefficient vector when the data matrix has full column rank. Small nonzero discrepancies in software are caused by floating-point roundoff.

Before using a discrete model in an adaptive controller, verify:

• the sample index and initial-condition conventions are consistent;
• the input delay is represented explicitly;
• the regressor contains only data available at the current sample;
• the normal matrix is sufficiently well conditioned;
• plant poles and zeros satisfy the assumptions of the intended controller;
• input and output units, scaling, saturation limits, and sign conventions are documented.

Problem 1 — Exact scalar discretization. For

\[ \dot y(t)=-2y(t)+3u(t),\qquad T_s=0.1, \]

derive the exact zero-order-hold recurrence.

Solution. Apply the scalar formula with \(\alpha=2\) and \(\beta=3\):

\[ y(k+1)=e^{-0.2}y(k)+\frac{3}{2}\left(1-e^{-0.2}\right)u(k). \]

Numerically,

\[ y(k+1)\approx0.81873075y(k)+0.27190387u(k). \]

The discrete pole is \(e^{-0.2}\), which lies inside the unit disk, preserving stability.

Problem 2 — Difference equation and regressor. Given

\[ \left(1-1.2q^{-1}+0.32q^{-2}\right)y(k) =q^{-2}\left(0.4-0.1q^{-1}\right)u(k)+v(k), \]

write the scalar recurrence, parameter vector, and regressor.

Solution.

\[ y(k)=1.2y(k-1)-0.32y(k-2)+0.4u(k-2)-0.1u(k-3)+v(k). \]

\[ \boldsymbol{\theta}= \begin{bmatrix}-1.2&0.32&0.4&-0.1\end{bmatrix}^{T},\qquad \boldsymbol{\phi}(k)= \begin{bmatrix}-y(k-1)&-y(k-2)&u(k-2)&u(k-3)\end{bmatrix}^{T}. \]

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

Problem 3 — Companion realization. Construct a state-space realization for

\[ G(z)=\frac{0.3z+0.1}{z^2-1.1z+0.28}. \]

Solution. Here \(a_1=-1.1\), \(a_2=0.28\), \(b_0=0.3\), and \(b_1=0.1\). Therefore

\[ \mathbf{A}_d=\begin{bmatrix}1.1&-0.28\\1&0\end{bmatrix},\quad \mathbf{B}_d=\begin{bmatrix}1\\0\end{bmatrix},\quad \mathbf{C}_d=\begin{bmatrix}0.3&0.1\end{bmatrix},\quad D_d=0. \]

The poles solve \(z^2-1.1z+0.28=0\), giving \(0.7\) and \(0.4\); hence the realization is stable.

Problem 4 — Parameter nonuniqueness under poor excitation. Suppose \(u(k)=0\) and the plant starts at rest, so \(y(k)=0\) for every measured sample. Can the four parameters of a second-order, one-delay ARX model be recovered?

Solution. Every regressor is the zero vector:

\[ \boldsymbol{\phi}(k)=\mathbf{0},\qquad \mathbf{\Phi}_N^{T}\mathbf{\Phi}_N=\mathbf{0}. \]

The matrix is not positive definite and every parameter vector predicts the same zero data. The plant is therefore unidentifiable from this experiment. Nonzero, sufficiently rich excitation or informative initial conditions are necessary.

Problem 5 — Effect of an incorrect delay. The true plant has \(d=2\), but a model is fitted with \(d=1\). Explain why increasing the number of samples does not necessarily remove the error.

Solution. The true input terms are \(u(k-2),u(k-3),\ldots\), whereas the fitted regressor uses \(u(k-1),u(k-2),\ldots\). Thus the true data-generating equation is generally outside the selected model class. More samples reduce variance around the best approximation within the wrong class, but they do not correct the structural shift. Residual-input correlation and a biased leading input coefficient typically remain.

Problem 6 — Stability proof for the lesson example. Show that the poles of

\[ z^2-1.50z+0.56 \]

are inside the unit disk.

Solution. Factorization gives

\[ z^2-1.50z+0.56=(z-0.8)(z-0.7). \]

Both magnitudes are less than one, so the homogeneous response is a linear combination of \(0.8^k\) and \(0.7^k\), each converging to zero. The plant is asymptotically and BIBO stable.

A discrete-time adaptive controller must begin from a causal, control-oriented plant model. Exact zero-order-hold discretization maps continuous state dynamics to sampled state dynamics. Polynomial input-output equations expose delay and produce the linear regression form \(y(k)=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}+v(k)\). State-space and polynomial realizations are equivalent when their transfer operators match, but they emphasize different information. Before any adaptation law is introduced, model order, delay, stability, leading-gain assumptions, signal scaling, and regressor rank must be verified.