After completing this lesson, students should be able to:

• test a candidate regressor using its finite-window excitation Gramian;
• explain why a scalar constant can be PE while a constant vector is usually non-PE;
• prove that sine/cosine quadrature provides two independent excitation directions;
• construct discrete-time PE sequences and rank-deficient counterexamples;
• identify the difference between finite excitation and persistent excitation; and
• estimate the smallest Gramian eigenvalue numerically without mistaking numerical noise for excitation.

Lessons 1 and 2 introduced persistent excitation and its role in parameter convergence. Here the object being tested is the regressor \( \boldsymbol{\phi}(t) \), not merely the externally commanded input. This distinction is essential: plant dynamics and filters can transform an apparently rich input into a lower-rank internal signal.

A bounded regressor \( \boldsymbol{\phi}:\mathbb{R}_{\ge 0}\to\mathbb{R}^{p} \) is persistently exciting when there exist constants \( T>0 \) and \( \alpha>0 \) such that, for every starting time \( t\ge 0 \),

\[ \mathbf{G}(t,T) :=\int_{t}^{t+T}\boldsymbol{\phi}(\tau)\boldsymbol{\phi}^{\mathsf T}(\tau)\,d\tau \succeq \alpha\mathbf{I}_{p}. \]

The matrix \( \mathbf{G}(t,T) \) is the excitation Gramian. The most direct numerical certificate is its smallest eigenvalue:

\[ \lambda_{\min}\!\left(\mathbf{G}(t,T)\right)\ge\alpha \quad\text{uniformly for all }t\ge 0. \]

The equivalent directional statement is especially useful. For every nonzero \( \mathbf{v}\in\mathbb{R}^{p} \),

\[ \mathbf{v}^{\mathsf T}\mathbf{G}(t,T)\mathbf{v} =\int_t^{t+T}\left(\mathbf{v}^{\mathsf T}\boldsymbol{\phi}(\tau)\right)^2d\tau \ge\alpha\lVert\mathbf{v}\rVert_2^2. \]

Thus, a regressor is PE only if every parameter-space direction receives a uniformly positive amount of energy in every admissible window.

flowchart TD
  A["Candidate input or regressor"] --> B["Choose window length T"]
  B --> C["Compute sliding Gramian G"]
  C --> D["Find lambda_min of G"]
  D --> E{"Positive lower bound \nfor every window?"}
  E -->|"yes"| F["PE for the tested regressor"]
  E -->|"no"| G["Non-PE or insufficient observation horizon"]
  G --> H["Check rank, decay, window choice, and numerical tolerance"]
        

3.1 Nonzero Constant Scalar — PE

Let \( p=1 \) and \( \phi(t)=c \), where \( c\ne 0 \). Then, for any selected \( T>0 \),

\[ G(t,T)=\int_t^{t+T}c^2d\tau=c^2T. \]

Hence \( \alpha=c^2T \) is independent of \( t \), so a nonzero constant is PE for a one-dimensional parameter. No change of direction is needed because only one direction exists.

3.2 Scalar Sinusoid — PE

Let \( \phi(t)=A\sin(\omega t+\psi) \), with \( A\ne 0 \) and \( \omega\ne 0 \). Choose one full period \( T=2\pi/|\omega| \). Then

\[ G(t,T)=A^2\int_t^{t+T}\sin^2(\omega\tau+\psi)d\tau =\frac{A^2T}{2}. \]

The result does not depend on phase or starting time. Therefore the scalar sinusoid is PE with \( \alpha=A^2T/2 \).

3.3 Scalar Pulse of Finite Duration — Non-PE

Suppose \( \phi(t)=1 \) for \( 0\le t\le t_f \) and zero afterward. Early windows contain information, but any window starting after \( t_f \) has zero Gramian. Thus no positive uniform lower bound exists. The signal is informative on a finite interval, but it is not persistently exciting over an infinite operating horizon.

Consider the two-dimensional regressor

\[ \boldsymbol{\phi}(t)= \begin{bmatrix}\sin(\omega t)\\\cos(\omega t)\end{bmatrix}, \qquad \omega\ne 0. \]

Over one period \( T=2\pi/|\omega| \), the Gramian is

\[ \mathbf{G}(t,T)= \int_t^{t+T} \begin{bmatrix} \sin^2(\omega\tau) & \sin(\omega\tau)\cos(\omega\tau)\\ \sin(\omega\tau)\cos(\omega\tau) & \cos^2(\omega\tau) \end{bmatrix}d\tau. \]

Orthogonality over a full period gives

\[ \int_t^{t+T}\sin^2(\omega\tau)d\tau=\frac{T}{2},\qquad \int_t^{t+T}\cos^2(\omega\tau)d\tau=\frac{T}{2}, \]

\[ \int_t^{t+T}\sin(\omega\tau)\cos(\omega\tau)d\tau=0. \]

Therefore

\[ \mathbf{G}(t,T)=\frac{T}{2}\mathbf{I}_2, \qquad \lambda_{\min}\!\left(\mathbf{G}(t,T)\right)=\frac{T}{2}. \]

The regressor is PE with \( \alpha=T/2 \). Geometrically, its direction rotates continuously and spans the complete plane during each period.

4.1 A Three-Dimensional Periodic PE Regressor

Adding a constant component produces

\[ \boldsymbol{\phi}(t)= \begin{bmatrix}1 & \sin(\omega t) & \cos(\omega t)\end{bmatrix}^{\mathsf T}. \]

Over one period, all cross terms integrate to zero and

\[ \mathbf{G}(t,T)= \operatorname{diag}\!\left(T,\frac{T}{2},\frac{T}{2}\right) \succeq\frac{T}{2}\mathbf{I}_3. \]

This example shows how an offset and a quadrature pair can independently excite three parameter directions.

5.1 Constant Vector — Rank-One Gramian

Let \( p>1 \) and \( \boldsymbol{\phi}(t)=\mathbf{c}\ne\mathbf{0} \). Then

\[ \mathbf{G}(t,T)=T\mathbf{c}\mathbf{c}^{\mathsf T}. \]

This outer product has rank one. Any nonzero vector \( \mathbf{v} \) satisfying \( \mathbf{v}^{\mathsf T}\mathbf{c}=0 \) yields

\[ \mathbf{v}^{\mathsf T}\mathbf{G}(t,T)\mathbf{v} =T\left(\mathbf{v}^{\mathsf T}\mathbf{c}\right)^2=0. \]

Hence the Gramian cannot dominate a positive multiple of \( \mathbf{I}_p \). A constant vector excites only its own fixed direction.

5.2 Time-Varying but Collinear Components — Non-PE

Consider

\[ \boldsymbol{\phi}(t)= \begin{bmatrix}\sin t\\2\sin t\end{bmatrix} =\sin t\begin{bmatrix}1\\2\end{bmatrix}. \]

Although both entries vary, the regressor always lies on one line. Over a full period,

\[ \mathbf{G}(t,2\pi)=\pi \begin{bmatrix}1&2\\2&4\end{bmatrix}, \]

whose eigenvalues are \( 0 \) and \( 5\pi \). Time variation alone is not enough; independent directions are required.

5.3 Exponentially Decaying Regressor — Nonuniform Information

Let

\[ \boldsymbol{\phi}(t)= \begin{bmatrix}e^{-t}\\e^{-2t}\end{bmatrix}. \]

For every fixed window length, all entries of the Gramian approach zero as the window moves to late time. Using \( \lVert\boldsymbol{\phi}(\tau)\rVert_2^2 \le e^{-2\tau}+e^{-4\tau} \),

\[ 0\le\lambda_{\min}\!\left(\mathbf{G}(t,T)\right) \le\operatorname{tr}\!\left(\mathbf{G}(t,T)\right) \le\int_t^{t+T}\left(e^{-2\tau}+e^{-4\tau}\right)d\tau. \]

The upper bound tends to zero as \( t\to\infty \). Therefore no constant \( \alpha>0 \) can hold uniformly.

A discrete regressor \( \boldsymbol{\phi}[k]\in\mathbb{R}^{p} \) is PE when there are an integer \( N\ge 1 \) and \( \alpha>0 \) such that

\[ \sum_{i=k}^{k+N-1}\boldsymbol{\phi}[i]\boldsymbol{\phi}^{\mathsf T}[i] \succeq\alpha\mathbf{I}_p \qquad\text{for every }k\ge 0. \]

6.1 Alternating Basis Vectors — PE

Let the period-two sequence be

\[ \boldsymbol{\phi}[2m]=\begin{bmatrix}1\\0\end{bmatrix},\qquad \boldsymbol{\phi}[2m+1]=\begin{bmatrix}0\\1\end{bmatrix}. \]

Every two-sample window contains both directions, so

\[ \sum_{i=k}^{k+1}\boldsymbol{\phi}[i]\boldsymbol{\phi}^{\mathsf T}[i] =\mathbf{I}_2. \]

The sequence is PE with \( N=2 \) and \( \alpha=1 \).

6.2 Repeated Direction — Non-PE

If \( \boldsymbol{\phi}[k]=[1\;1]^{\mathsf T} \) for all \( k \), then over an \( N \)-sample window,

\[ \sum_{i=k}^{k+N-1}\boldsymbol{\phi}[i]\boldsymbol{\phi}^{\mathsf T}[i] =N\begin{bmatrix}1&1\\1&1\end{bmatrix}. \]

Its eigenvalues are \( 0 \) and \( 2N \); therefore it is non-PE for two unknown parameters.

Adaptive laws typically see filtered states, outputs, and commands rather than the raw command alone. In symbolic form,

\[ u(t)\quad\mathop{\longmapsto}^{\text{plant and filters} }\quad \boldsymbol{\phi}(t)\quad\mathop{\longmapsto}^{\text{adaptive law} }\quad \dot{\hat{\boldsymbol{\theta} } }(t). \]

A multisine command can fail to produce a PE regressor if a relevant frequency is cancelled, strongly attenuated, or repeated through dependent channels. Conversely, a simple command may produce a richer regressor when coupled with measured states that have independent transient or periodic content. The correct test is therefore applied to the actual \( \boldsymbol{\phi}(t) \) entering the update law.

For safe experiment design, excitation amplitude must also respect actuator limits, bandwidth, noise sensitivity, and unmodeled dynamics. PE is an information condition, not permission to inject arbitrarily large or high-frequency probing signals.

flowchart TD
  U["Command u"] --> P["Plant dynamics"]
  P --> M["Measured states and outputs"]
  U --> F["Reference or command filters"]
  M --> R["Build regressor phi"]
  F --> R
  R --> G["Sliding Gramian test"]
  G --> Q{"All directions excited?"}
  Q -->|"yes"| C["Parameter convergence may be possible"]
  Q -->|"no"| S["Only a parameter subspace is identifiable"]
        

With samples \( t_j=t+j\Delta t \), a practical approximation is

\[ \widehat{\mathbf{G} }(t,T) =\sum_{j=0}^{M}w_j \boldsymbol{\phi}(t_j)\boldsymbol{\phi}^{\mathsf T}(t_j)\Delta t, \]

where trapezoidal weights use \( w_0=w_M=1/2 \) and \( w_j=1 \) otherwise. Define the empirical score

\[ \widehat{\alpha}_{\mathrm{obs} } =\min_{t\in\mathcal{T}_{\mathrm{obs} } } \lambda_{\min}\!\left(\widehat{\mathbf{G} }(t,T)\right). \]

A positive score over finite data is evidence of PE on the observed horizon, not a proof for all future time. Use a tolerance scaled to signal energy, inspect the Gramian condition number, and repeat the test for several window lengths. A tiny positive eigenvalue caused by floating-point roundoff should be treated as zero.

Chapter10_Lesson3.py

import math
from typing import Callable

import numpy as np


def regressor(case_name: str, t: float) -> np.ndarray:
    """Return the regressor used in each PE/non-PE example."""
    if case_name == "scalar_constant":
        return np.array([1.0])
    if case_name == "sin_cos":
        return np.array([math.sin(t), math.cos(t)])
    if case_name == "constant_vector":
        return np.array([1.0, 1.0])
    if case_name == "collinear_sine":
        return np.array([math.sin(t), 2.0 * math.sin(t)])
    if case_name == "decaying":
        return np.array([math.exp(-t), math.exp(-2.0 * t)])
    raise ValueError(f"Unknown case: {case_name}")


def window_gramian(
    phi: Callable[[float], np.ndarray],
    start: float,
    window: float,
    dt: float,
) -> np.ndarray:
    """Approximate integral(phi phi^T) over one sliding window."""
    sample_count = max(2, int(round(window / dt)) + 1)
    times = np.linspace(start, start + window, sample_count)
    values = np.vstack([phi(t) for t in times])
    # Trapezoidal integration, applied entry by entry.
    products = np.einsum("ni,nj->nij", values, values)
    return np.trapezoid(products, times, axis=0)


def scan_case(
    case_name: str,
    horizon: float = 20.0,
    window: float = 2.0 * math.pi,
    dt: float = 0.002,
    start_step: float = 0.25,
) -> tuple[float, float]:
    starts = np.arange(0.0, horizon - window + 1e-12, start_step)
    smallest = []
    for start in starts:
        gramian = window_gramian(
            lambda tau: regressor(case_name, tau), start, window, dt
        )
        smallest.append(float(np.linalg.eigvalsh(gramian)[0]))
    return min(smallest), max(smallest)


def discrete_window_gramian(sequence: list[np.ndarray], start: int, length: int) -> np.ndarray:
    dimension = sequence[0].size
    gramian = np.zeros((dimension, dimension))
    for k in range(start, start + length):
        vector = sequence[k % len(sequence)]
        gramian += np.outer(vector, vector)
    return gramian


def main() -> None:
    cases = [
        "scalar_constant",
        "sin_cos",
        "constant_vector",
        "collinear_sine",
        "decaying",
    ]

    print("Continuous-time sliding-Gramian test")
    for case_name in cases:
        worst, best = scan_case(case_name)
        classification = "numerically PE" if worst > 1e-3 else "non-PE"
        print(
            f"{case_name:18s} min(lambda_min)={worst: .6e}, "
            f"max(lambda_min)={best: .6e} -> {classification}"
        )

    alternating_basis = [np.array([1.0, 0.0]), np.array([0.0, 1.0])]
    repeated_direction = [np.array([1.0, 1.0])]

    print("\nDiscrete-time two-sample tests")
    for name, sequence in [
        ("alternating_basis", alternating_basis),
        ("repeated_direction", repeated_direction),
    ]:
        gramian = discrete_window_gramian(sequence, start=0, length=2)
        eigenvalues = np.linalg.eigvalsh(gramian)
        print(f"{name:18s} G=\n{gramian}\neigenvalues={eigenvalues}\n")


if __name__ == "__main__":
    main()

Chapter10_Lesson3.cpp

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

struct Matrix2 {
    double a11{0.0};
    double a12{0.0};
    double a22{0.0};
};

std::vector<double> regressor(const std::string& caseName, double t) {
    if (caseName == "scalar_constant") {
        return {1.0};
    }
    if (caseName == "sin_cos") {
        return {std::sin(t), std::cos(t)};
    }
    if (caseName == "constant_vector") {
        return {1.0, 1.0};
    }
    if (caseName == "collinear_sine") {
        return {std::sin(t), 2.0 * std::sin(t)};
    }
    if (caseName == "decaying") {
        return {std::exp(-t), std::exp(-2.0 * t)};
    }
    throw std::invalid_argument("Unknown case: " + caseName);
}

double smallestEigenvalue(const std::vector<double>& gramian) {
    if (gramian.size() == 1U) {
        return gramian[0];
    }
    const double a = gramian[0];
    const double b = gramian[1];
    const double d = gramian[3];
    const double trace = a + d;
    const double discriminant = std::sqrt((a - d) * (a - d) + 4.0 * b * b);
    return 0.5 * (trace - discriminant);
}

std::vector<double> windowGramian(
    const std::string& caseName,
    double start,
    double window,
    double dt
) {
    const std::vector<double> first = regressor(caseName, start);
    const std::size_t dimension = first.size();
    std::vector<double> gramian(dimension * dimension, 0.0);
    const int intervals = std::max(1, static_cast<int>(std::round(window / dt)));
    const double actualDt = window / static_cast<double>(intervals);

    for (int k = 0; k <= intervals; ++k) {
        const double t = start + actualDt * static_cast<double>(k);
        const double weight = (k == 0 || k == intervals) ? 0.5 : 1.0;
        const std::vector<double> phi = regressor(caseName, t);
        for (std::size_t i = 0; i < dimension; ++i) {
            for (std::size_t j = 0; j < dimension; ++j) {
                gramian[i * dimension + j] += weight * phi[i] * phi[j] * actualDt;
            }
        }
    }
    return gramian;
}

std::pair<double, double> scanCase(
    const std::string& caseName,
    double horizon = 20.0,
    double window = 2.0 * std::acos(-1.0),
    double dt = 0.002,
    double startStep = 0.25
) {
    double worst = 1.0e300;
    double best = -1.0e300;
    for (double start = 0.0; start + window <= horizon + 1.0e-12; start += startStep) {
        const double value = smallestEigenvalue(windowGramian(caseName, start, window, dt));
        worst = std::min(worst, value);
        best = std::max(best, value);
    }
    return {worst, best};
}

Matrix2 discreteGramian(
    const std::vector<std::vector<double>>& sequence,
    int start,
    int length
) {
    Matrix2 gramian;
    for (int k = start; k < start + length; ++k) {
        const auto& v = sequence[static_cast<std::size_t>(k) % sequence.size()];
        gramian.a11 += v[0] * v[0];
        gramian.a12 += v[0] * v[1];
        gramian.a22 += v[1] * v[1];
    }
    return gramian;
}

void printDiscrete(const std::string& name, const std::vector<std::vector<double>>& sequence) {
    const Matrix2 g = discreteGramian(sequence, 0, 2);
    const std::vector<double> flat = {g.a11, g.a12, g.a12, g.a22};
    const double lambdaMin = smallestEigenvalue(flat);
    const double lambdaMax = g.a11 + g.a22 - lambdaMin;
    std::cout << name << "\nG = [[" << g.a11 << ", " << g.a12 << "], ["
              << g.a12 << ", " << g.a22 << "]]\n"
              << "eigenvalues = [" << lambdaMin << ", " << lambdaMax << "]\n\n";
}

int main() {
    const std::vector<std::string> cases = {
        "scalar_constant",
        "sin_cos",
        "constant_vector",
        "collinear_sine",
        "decaying"
    };

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "Continuous-time sliding-Gramian test\n";
    for (const auto& caseName : cases) {
        const auto [worst, best] = scanCase(caseName);
        const std::string classification = worst > 1.0e-3 ? "numerically PE" : "non-PE";
        std::cout << std::setw(18) << std::left << caseName
                  << " min(lambda_min)=" << worst
                  << ", max(lambda_min)=" << best
                  << " -> " << classification << '\n';
    }

    std::cout << "\nDiscrete-time two-sample tests\n";
    printDiscrete("alternating_basis", { {1.0, 0.0}, {0.0, 1.0} });
    printDiscrete("repeated_direction", { {1.0, 1.0} });
    return 0;
}

Chapter10_Lesson3.java

import java.util.List;

public final class Chapter10_Lesson3 {
    private Chapter10_Lesson3() {
    }

    private static double[] regressor(String caseName, double t) {
        return switch (caseName) {
            case "scalar_constant" -> new double[] {1.0};
            case "sin_cos" -> new double[] {Math.sin(t), Math.cos(t)};
            case "constant_vector" -> new double[] {1.0, 1.0};
            case "collinear_sine" -> new double[] {Math.sin(t), 2.0 * Math.sin(t)};
            case "decaying" -> new double[] {Math.exp(-t), Math.exp(-2.0 * t)};
            default -> throw new IllegalArgumentException("Unknown case: " + caseName);
        };
    }

    private static double smallestEigenvalue(double[] gramian) {
        if (gramian.length == 1) {
            return gramian[0];
        }
        double a = gramian[0];
        double b = gramian[1];
        double d = gramian[3];
        double trace = a + d;
        double discriminant = Math.sqrt((a - d) * (a - d) + 4.0 * b * b);
        return 0.5 * (trace - discriminant);
    }

    private static double[] windowGramian(
            String caseName,
            double start,
            double window,
            double dt) {
        int dimension = regressor(caseName, start).length;
        double[] gramian = new double[dimension * dimension];
        int intervals = Math.max(1, (int) Math.round(window / dt));
        double actualDt = window / intervals;

        for (int k = 0; k <= intervals; k++) {
            double t = start + actualDt * k;
            double weight = (k == 0 || k == intervals) ? 0.5 : 1.0;
            double[] phi = regressor(caseName, t);
            for (int i = 0; i < dimension; i++) {
                for (int j = 0; j < dimension; j++) {
                    gramian[i * dimension + j] += weight * phi[i] * phi[j] * actualDt;
                }
            }
        }
        return gramian;
    }

    private static double[] scanCase(
            String caseName,
            double horizon,
            double window,
            double dt,
            double startStep) {
        double worst = Double.POSITIVE_INFINITY;
        double best = Double.NEGATIVE_INFINITY;
        for (double start = 0.0; start + window <= horizon + 1.0e-12; start += startStep) {
            double value = smallestEigenvalue(windowGramian(caseName, start, window, dt));
            worst = Math.min(worst, value);
            best = Math.max(best, value);
        }
        return new double[] {worst, best};
    }

    private static double[] discreteGramian(List<double[]> sequence, int start, int length) {
        double[] gramian = new double[4];
        for (int k = start; k < start + length; k++) {
            double[] v = sequence.get(k % sequence.size());
            gramian[0] += v[0] * v[0];
            gramian[1] += v[0] * v[1];
            gramian[2] += v[1] * v[0];
            gramian[3] += v[1] * v[1];
        }
        return gramian;
    }

    private static void printDiscrete(String name, List<double[]> sequence) {
        double[] g = discreteGramian(sequence, 0, 2);
        double lambdaMin = smallestEigenvalue(g);
        double lambdaMax = g[0] + g[3] - lambdaMin;
        System.out.printf(
                "%s%nG = [[%.6f, %.6f], [%.6f, %.6f]]%n"
                        + "eigenvalues = [%.6f, %.6f]%n%n",
                name, g[0], g[1], g[2], g[3], lambdaMin, lambdaMax);
    }

    public static void main(String[] args) {
        List<String> cases = List.of(
                "scalar_constant",
                "sin_cos",
                "constant_vector",
                "collinear_sine",
                "decaying");

        System.out.println("Continuous-time sliding-Gramian test");
        for (String caseName : cases) {
            double[] range = scanCase(caseName, 20.0, 2.0 * Math.PI, 0.002, 0.25);
            String classification = range[0] > 1.0e-3 ? "numerically PE" : "non-PE";
            System.out.printf(
                    "%-18s min(lambda_min)=% .6e, max(lambda_min)=% .6e -> %s%n",
                    caseName, range[0], range[1], classification);
        }

        System.out.println("\nDiscrete-time two-sample tests");
        printDiscrete(
                "alternating_basis",
                List.of(new double[] {1.0, 0.0}, new double[] {0.0, 1.0}));
        printDiscrete("repeated_direction", List.of(new double[] {1.0, 1.0}));
    }
}

The script runs the MATLAB tests directly and, when Simulink is licensed, programmatically creates Chapter10_Lesson3_PE_Simulink.slx.

Chapter10_Lesson3.m

%% Chapter10_Lesson3.m
% Examples of persistently exciting and non-persistently exciting regressors.
clear; clc;

caseNames = {
    'scalar_constant', ...
    'sin_cos', ...
    'constant_vector', ...
    'collinear_sine', ...
    'decaying'};

horizon = 20.0;
window = 2*pi;
dt = 0.002;
startStep = 0.25;

fprintf('Continuous-time sliding-Gramian test\n');
for i = 1:numel(caseNames)
    [worstValue, bestValue] = scanCase(...
        caseNames{i}, horizon, window, dt, startStep);
    if worstValue > 1e-3
        classification = 'numerically PE';
    else
        classification = 'non-PE';
    end
    fprintf('%-18s min(lambda_min)=% .6e, max(lambda_min)=% .6e -> %s\n', ...
        caseNames{i}, worstValue, bestValue, classification);
end

fprintf('\nDiscrete-time two-sample tests\n');
alternatingBasis = [1 0; 0 1];
repeatedDirection = [1 1; 1 1];
printDiscrete('alternating_basis', alternatingBasis);
printDiscrete('repeated_direction', repeatedDirection);

%% Optional Simulink demonstration
% This section creates a model that generates phi(t) = [sin(t); cos(t)],
% exports the regressor, and checks the Gramian over one period.
if license('test', 'Simulink')
    modelName = 'Chapter10_Lesson3_PE_Simulink';
    buildSimulinkDemo(modelName);
    simulationOutput = sim(modelName, ...
        'StopTime', num2str(4*pi), ...
        'ReturnWorkspaceOutputs', 'on');
    phiStructure = simulationOutput.get('phi_sim');
    timeSim = phiStructure.time;
    phiSim = phiStructure.signals.values;

    onePeriod = timeSim <= 2*pi + 10*eps;
    selectedTime = timeSim(onePeriod);
    selectedPhi = phiSim(onePeriod, :);
    productsSim = zeros(2, 2, numel(selectedTime));
    for k = 1:numel(selectedTime)
        vector = selectedPhi(k, :).';
        productsSim(:, :, k) = vector*vector.';
    end
    gramianSim = trapz(selectedTime, productsSim, 3);
    fprintf('\nSimulink sin/cos Gramian eigenvalues: [%g, %g]\n', ...
        eig(gramianSim));
    open_system(modelName);
else
    fprintf('\nSimulink is not licensed; the MATLAB tests still ran.\n');
end

function vector = regressor(caseName, t)
    switch caseName
        case 'scalar_constant'
            vector = 1.0;
        case 'sin_cos'
            vector = [sin(t); cos(t)];
        case 'constant_vector'
            vector = [1.0; 1.0];
        case 'collinear_sine'
            vector = [sin(t); 2.0*sin(t)];
        case 'decaying'
            vector = [exp(-t); exp(-2.0*t)];
        otherwise
            error('Unknown case: %s', caseName);
    end
end

function gramian = windowGramian(caseName, startTime, window, dt)
    time = startTime:dt:(startTime + window);
    if time(end) < startTime + window
        time(end + 1) = startTime + window;
    end
    firstVector = regressor(caseName, time(1));
    dimension = numel(firstVector);
    products = zeros(dimension, dimension, numel(time));
    for k = 1:numel(time)
        vector = regressor(caseName, time(k));
        products(:, :, k) = vector*vector.';
    end
    gramian = trapz(time, products, 3);
end

function [worstValue, bestValue] = scanCase(...
        caseName, horizon, window, dt, startStep)
    starts = 0:startStep:(horizon - window);
    values = zeros(size(starts));
    for k = 1:numel(starts)
        gramian = windowGramian(caseName, starts(k), window, dt);
        eigenvalues = eig((gramian + gramian.')/2);
        values(k) = min(eigenvalues);
    end
    worstValue = min(values);
    bestValue = max(values);
end

function printDiscrete(name, samples)
    gramian = samples.'*samples;
    eigenvalues = sort(eig(gramian));
    fprintf('%s\n', name);
    disp(gramian);
    fprintf('eigenvalues = [%g, %g]\n\n', eigenvalues(1), eigenvalues(2));
end

function buildSimulinkDemo(modelName)
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if exist([modelName '.slx'], 'file')
        delete([modelName '.slx']);
    end

    new_system(modelName);
    set_param(modelName, 'Solver', 'ode45');

    add_block('simulink/Sources/Sine Wave', [modelName '/sin_t'], ...
        'Amplitude', '1', 'Frequency', '1', 'Phase', '0', ...
        'Position', [40 40 110 70]);
    add_block('simulink/Sources/Sine Wave', [modelName '/cos_t'], ...
        'Amplitude', '1', 'Frequency', '1', 'Phase', 'pi/2', ...
        'Position', [40 110 110 140]);
    add_block('simulink/Signal Routing/Mux', [modelName '/Regressor'], ...
        'Inputs', '2', 'Position', [170 55 175 125]);
    add_block('simulink/Sinks/To Workspace', [modelName '/Export_phi'], ...
        'VariableName', 'phi_sim', 'SaveFormat', 'Structure With Time', ...
        'Position', [240 75 340 105]);

    add_line(modelName, 'sin_t/1', 'Regressor/1');
    add_line(modelName, 'cos_t/1', 'Regressor/2');
    add_line(modelName, 'Regressor/1', 'Export_phi/1');
    save_system(modelName);
end

Chapter10_Lesson3.nb

Notebook[{
 Cell["Chapter 10, Lesson 3: Examples of PE and Non-PE Inputs", "Title"],
 Cell["Sliding excitation-Gramian calculations for continuous- and discrete-time regressors.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];
regressor[\"scalar_constant\", t_] := {1.0};
regressor[\"sin_cos\", t_] := {Sin[t], Cos[t]};
regressor[\"constant_vector\", t_] := {1.0, 1.0};
regressor[\"collinear_sine\", t_] := {Sin[t], 2.0 Sin[t]};
regressor[\"decaying\", t_] := {Exp[-t], Exp[-2.0 t]};"], "Input"],
 Cell[BoxData[
"windowGramian[case_, start_, window_] := NIntegrate[
  Outer[Times, regressor[case, tau], regressor[case, tau]],
  {tau, start, start + window},
  Method -> {\"GlobalAdaptive\", \"SymbolicProcessing\" -> 0}
];
minimumWindowEigenvalue[case_, horizon_: 20.0, window_: 2.0 Pi, step_: 0.25] :=
 Module[{starts, values},
  starts = Range[0.0, horizon - window, step];
  values = Min[Eigenvalues[(# + Transpose[#])/2]] & /@
    (windowGramian[case, #, window] & /@ starts);
  {Min[values], Max[values]}
];"], "Input"],
 Cell[BoxData[
"cases = {\"scalar_constant\", \"sin_cos\", \"constant_vector\",
   \"collinear_sine\", \"decaying\"};
continuousResults = AssociationMap[minimumWindowEigenvalue, cases];
Dataset[KeyValueMap[
  <|\"Case\" -> #1, \"WorstLambdaMin\" -> #2[[1]],
    \"BestLambdaMin\" -> #2[[2]],
    \"Classification\" -> If[#2[[1]] > 10^-3, \"numerically PE\", \"non-PE\"]|> &,
  continuousResults
]]"], "Input"],
 Cell[BoxData[
"discreteGramian[sequence_, start_, length_] := Total[
  Outer[Times, #, #] & /@ Table[
    sequence[[Mod[k, Length[sequence]] + 1]],
    {k, start, start + length - 1}
  ]
];
alternatingBasis = { {1.0, 0.0}, {0.0, 1.0} };
repeatedDirection = { {1.0, 1.0} };
{
  discreteGramian[alternatingBasis, 0, 2],
  Eigenvalues[discreteGramian[alternatingBasis, 0, 2]],
  discreteGramian[repeatedDirection, 0, 2],
  Eigenvalues[discreteGramian[repeatedDirection, 0, 2]]
}"], "Input"],
 Cell[BoxData[
"Plot[
 Evaluate[{
   Min[Eigenvalues[windowGramian[\"sin_cos\", s, 2 Pi]]],
   Min[Eigenvalues[windowGramian[\"decaying\", s, 2 Pi]]]
 }],
 {s, 0, 12},
 PlotLegends -> {\"sin/cos\", \"decaying\"},
 AxesLabel -> {\"window start\", \"minimum eigenvalue\"},
 PlotRange -> All
]"], "Input"]
}, WindowTitle -> "Chapter10_Lesson3"]

With \( T=2\pi \), the scalar constant has Gramian \( 2\pi \), while the sine/cosine regressor has the exact Gramian \( \pi\mathbf{I}_2 \). Their smallest eigenvalues remain uniformly positive as the window slides.

The constant vector and collinear-sine cases produce a zero eigenvalue because their Gramians are rank one. The decaying regressor can have a positive smallest eigenvalue in early windows, but its worst-window value approaches zero as the observation horizon expands. This is the numerical signature of finite early information without persistent excitation.

The discrete alternating-basis sequence gives \( \mathbf{G}=\mathbf{I}_2 \) in every two-sample window. Repeating \( [1\;1]^{\mathsf T} \) instead gives a rank-one matrix, regardless of window length.

Problem 1: Constant Scalar Versus Constant Vector

Determine whether \( \phi(t)=3 \) is PE for one unknown parameter and whether \( \boldsymbol{\phi}(t)=[3\;4]^{\mathsf T} \) is PE for two unknown parameters.

Solution: For the scalar signal,

\[ G(t,T)=\int_t^{t+T}9d\tau=9T, \]

so it is PE for every \( T>0 \). For the vector signal,

\[ \mathbf{G}(t,T)=T \begin{bmatrix}9&12\\12&16\end{bmatrix}. \]

Its determinant is \( 9\cdot16-12^2=0 \), so one eigenvalue is zero. The vector is non-PE because it excites only the direction \( [3\;4]^{\mathsf T} \).

Problem 2: Frequency Scaling of the Quadrature Regressor

For \( \boldsymbol{\phi}(t)=[\sin(5t)\;\cos(5t)]^{\mathsf T} \), find a valid window and lower bound.

Solution: One period is \( T=2\pi/5 \). Orthogonality gives

\[ \mathbf{G}(t,T)=\frac{T}{2}\mathbf{I}_2 =\frac{\pi}{5}\mathbf{I}_2. \]

Therefore a valid PE pair is \( (T,\alpha)=(2\pi/5,\pi/5) \).

Problem 3: Phase-Shifted Components

Analyze \( \boldsymbol{\phi}(t)=[\sin t\;\sin(t+\delta)]^{\mathsf T} \) over \( T=2\pi \). For which values of \( \delta \) is it PE?

Solution: The Gramian is

\[ \mathbf{G}(t,2\pi)=\pi \begin{bmatrix}1&\cos\delta\\\cos\delta&1\end{bmatrix}. \]

Its eigenvalues are

\[ \lambda_{1,2}=\pi\left(1\pm|\cos\delta|\right). \]

The smallest eigenvalue is positive exactly when \( |\cos\delta|<1 \). Thus PE fails for \( \delta=m\pi \), where the two components are equal or opposite and therefore collinear. The strongest conditioning occurs at \( \delta=\pi/2+m\pi \).

Problem 4: A Discrete Period-Three Sequence

Let \( \boldsymbol{\phi}[3m]=[1\;0]^{\mathsf T} \), \( \boldsymbol{\phi}[3m+1]=[1\;1]^{\mathsf T} \), and \( \boldsymbol{\phi}[3m+2]=[0\;1]^{\mathsf T} \). Show that it is PE.

Solution: Every length-three window contains all three vectors, possibly in cyclic order. Therefore

\[ \mathbf{G}= \begin{bmatrix}1\\0\end{bmatrix}\begin{bmatrix}1&0\end{bmatrix} +\begin{bmatrix}1\\1\end{bmatrix}\begin{bmatrix}1&1\end{bmatrix} +\begin{bmatrix}0\\1\end{bmatrix}\begin{bmatrix}0&1\end{bmatrix} =\begin{bmatrix}2&1\\1&2\end{bmatrix}. \]

Its eigenvalues are \( 1 \) and \( 3 \). Hence the sequence is PE with \( N=3 \) and \( \alpha=1 \).

Problem 5: Finite Excitation Is Not PE

Let \( \boldsymbol{\phi}(t)=[\sin t\;\cos t]^{\mathsf T} \) for \( 0\le t\le 4\pi \) and \( \boldsymbol{\phi}(t)=\mathbf{0} \) afterward. Explain why accurate parameter estimates may be obtained during the initial interval although the signal is not PE on \( [0,\infty) \).

Solution: The first two periods generate the positive information matrix

\[ \int_0^{4\pi}\boldsymbol{\phi}(\tau) \boldsymbol{\phi}^{\mathsf T}(\tau)d\tau=2\pi\mathbf{I}_2. \]

This can substantially reduce parameter error during that interval. However, any window beginning after \( 4\pi \) has zero Gramian, so the uniform all-time PE inequality fails. The example separates information accumulated over one experiment from persistent information supplied indefinitely.

Problem 6: Choosing a Numerical Tolerance

A computed Gramian has eigenvalues \( 3.14 \) and \( 2.0\times10^{-13} \). Should the regressor be classified as PE in double precision?

Solution: Not from this computation. The eigenvalue ratio is

\[ \frac{2.0\times10^{-13} }{3.14}\approx6.4\times10^{-14}, \]

which is consistent with rank deficiency contaminated by floating-point and quadrature errors. The tolerance should be relative to the largest eigenvalue and expected measurement noise. The regressor should be treated as numerically non-PE unless independent analysis establishes a positive lower bound.

Persistent excitation is a uniform full-rank condition on a sliding information matrix. A nonzero scalar constant is PE for one parameter, but a constant vector is rank deficient when several parameters must be separated. Quadrature sinusoids are PE because their Gramian is a positive multiple of the identity. Collinear, finite-duration, and decaying regressors are non-PE for distinct reasons: missing directions, missing late-time excitation, or a vanishing lower bound. In discrete time, periodically cycling through basis directions gives a transparent PE construction. Numerical eigenvalue tests are useful diagnostics, but finite data cannot by itself prove an all-time property.