After completing this lesson, students should be able to:

• distinguish tracking-error or prediction-error convergence from convergence of the parameter error;
• express PE through a sliding-window regressor Gramian and interpret its minimum eigenvalue;
• prove parameter convergence for a continuous-time gradient estimator under bounded PE regressors;
• identify the invariant parameter directions that remain under a rank-deficient regressor;
• explain how PE enters the parameter-convergence argument for MRAC;
• numerically compare PE and non-PE experiments and evaluate an empirical excitation metric.

The lesson uses the Lyapunov and Barbalat arguments introduced in Chapter 3, the MRAC error structures developed in Chapters 4–7, and the PE definition from Lesson 1 of this chapter. No new adaptive-controller architecture is required.

Consider an exact, noise-free, linearly parameterized relation

\[ y(t)=\boldsymbol{\phi}^{T}(t)\boldsymbol{\theta}^{*}, \qquad \boldsymbol{\theta}^{*}\in\mathbb{R}^{p}, \qquad \boldsymbol{\phi}(t)\in\mathbb{R}^{p}, \]

where \( \boldsymbol{\theta}^{*} \) is an unknown constant parameter vector. An estimate \( \hat{\boldsymbol{\theta} }(t) \) produces

\[ \hat{y}(t)=\boldsymbol{\phi}^{T}(t)\hat{\boldsymbol{\theta} }(t), \qquad \varepsilon(t)=\hat{y}(t)-y(t). \]

Define the parameter error as \( \tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*} \). Then

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

The crucial logical question is: does \( \varepsilon(t)\to 0 \) force \( \tilde{\boldsymbol{\theta} }(t)\to\mathbf{0} \)? In general, the answer is no. A scalar inner product can vanish while the vector multiplying it remains nonzero.

flowchart TD
  A["Adaptive law is stable"] --> B["Prediction or tracking error tends to zero"]
  B --> C{"Does the regressor excite \nevery parameter direction?"}
  C -->|"Yes: window Gramian is uniformly positive definite"| D["No indistinguishable direction remains"]
  D --> E["Parameter error tends to zero"]
  C -->|"No: window Gramian is rank deficient"| F["An unobservable parameter subspace remains"]
  F --> G["Correct output may coexist with incorrect parameters"]
        

Use the continuous-time gradient law

\[ \dot{\hat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}(t)\varepsilon(t), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ 0. \]

Because the true parameter is constant,

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}(t) \boldsymbol{\phi}^{T}(t)\tilde{\boldsymbol{\theta} }. \]

Choose the quadratic Lyapunov function

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

Its derivative is

\[ \begin{aligned} \dot{V} &= \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } } \\ &= -\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}\boldsymbol{\phi}^{T} \tilde{\boldsymbol{\theta} } \\ &= -\left( \boldsymbol{\phi}^{T}\tilde{\boldsymbol{\theta} } \right)^{2} =-\varepsilon^{2} \le 0. \end{aligned} \]

Therefore, the parameter estimate is bounded and \( V(t) \) is nonincreasing. Under the usual boundedness and uniform-continuity conditions, Barbalat's lemma also yields \( \varepsilon(t)\to 0 \). However, \( \dot V=-\varepsilon^{2} \) is only negative semidefinite with respect to the full parameter error. It does not yet prove \( \tilde{\boldsymbol{\theta} }\to\mathbf{0} \).

Over a time window of length \( T \), define the regressor information matrix, or window Gramian,

\[ \mathbf{G}_{\phi}(t,T) = \int_{t}^{t+T} \boldsymbol{\phi}(\tau) \boldsymbol{\phi}^{T}(\tau)\,d\tau. \]

For any fixed direction \( \mathbf{v}\in\mathbb{R}^{p} \),

\[ \mathbf{v}^{T}\mathbf{G}_{\phi}(t,T)\mathbf{v} = \int_{t}^{t+T} \left( \boldsymbol{\phi}^{T}(\tau)\mathbf{v} \right)^{2}d\tau. \]

Hence, \( \mathbf{G}_{\phi}(t,T) \) is singular exactly when at least one nonzero direction is invisible over that window. If \( \mathbf{v}\neq\mathbf{0} \) satisfies

\[ \boldsymbol{\phi}^{T}(\tau)\mathbf{v}=0 \quad\text{for almost every }\tau\in[t,t+T], \]

then the parameter vectors \( \boldsymbol{\theta}^{*} \) and \( \boldsymbol{\theta}^{*}+\mathbf{v} \) generate the same output throughout that window. No estimator using only those data can distinguish them.

4.1 Windowed least-squares interpretation

Freeze a candidate parameter vector \( \boldsymbol{\vartheta} \) over one window and define

\[ J_{t,T}(\boldsymbol{\vartheta}) = \frac{1}{2} \int_{t}^{t+T} \left[ y(\tau) -\boldsymbol{\phi}^{T}(\tau)\boldsymbol{\vartheta} \right]^{2}d\tau. \]

In the exact regression,

\[ J_{t,T}(\boldsymbol{\vartheta}) = \frac{1}{2} \left( \boldsymbol{\vartheta}-\boldsymbol{\theta}^{*} \right)^{T} \mathbf{G}_{\phi}(t,T) \left( \boldsymbol{\vartheta}-\boldsymbol{\theta}^{*} \right). \]

A positive-definite Gramian therefore gives a unique minimizer. A singular Gramian gives a flat direction and a continuum of parameters that fit the data equally well.

A bounded regressor is persistently exciting of order \( p \) when there exist constants \( T>0 \), \( \alpha>0 \), and \( \beta\ge\alpha \) such that, for every \( t\ge t_{0} \),

\[ \alpha\mathbf{I}_{p} \le \int_{t}^{t+T} \boldsymbol{\phi}(\tau) \boldsymbol{\phi}^{T}(\tau)\,d\tau \le \beta\mathbf{I}_{p}. \]

Equivalently,

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

The lower bound is the essential part for identifiability. It means that every unit parameter direction contributes at least \( \alpha \) units of squared regressor projection during every window. In terms of the windowed cost,

\[ J_{t,T}(\boldsymbol{\vartheta}) \ge \frac{\alpha}{2} \left\| \boldsymbol{\vartheta}-\boldsymbol{\theta}^{*} \right\|^{2}. \]

Thus, PE is a uniform strong-identifiability condition, not merely a statement that the signal is nonzero or oscillatory.

Theorem. Consider

\[ \dot{\tilde{\boldsymbol{\theta} } } = -\boldsymbol{\Gamma} \boldsymbol{\phi}(t) \boldsymbol{\phi}^{T}(t) \tilde{\boldsymbol{\theta} }, \qquad \boldsymbol{\Gamma}\succ 0. \]

Suppose that \( \boldsymbol{\phi}(t) \) is bounded, piecewise continuous, and PE. Then the origin \( \tilde{\boldsymbol{\theta} }=\mathbf{0} \) is uniformly exponentially stable. Consequently, constants \( C\ge 1 \) and \( \lambda>0 \) exist such that

\[ \left\| \tilde{\boldsymbol{\theta} }(t) \right\| \le C e^{-\lambda(t-t_{0})} \left\| \tilde{\boldsymbol{\theta} }(t_{0}) \right\|, \qquad t\ge t_{0}. \]

Therefore, \( \hat{\boldsymbol{\theta} }(t) \to\boldsymbol{\theta}^{*} \).

The statement is stronger than asymptotic convergence: PE supplies a uniform contraction over repeated windows. The exact rate depends on the excitation level, the window length, the regressor bound, and the adaptation-gain matrix.

7.1 Symmetric normalization

Define

\[ \mathbf{z} = \boldsymbol{\Gamma}^{-1/2} \tilde{\boldsymbol{\theta} }, \qquad \boldsymbol{\psi}(t) = \boldsymbol{\Gamma}^{1/2} \boldsymbol{\phi}(t). \]

The parameter-error system becomes

\[ \dot{\mathbf{z} } = -\boldsymbol{\psi}(t) \boldsymbol{\psi}^{T}(t)\mathbf{z}. \]

Let \( W=\frac{1}{2}\|\mathbf{z}\|^{2} \). Then

\[ \dot W = -\left| \boldsymbol{\psi}^{T}(t)\mathbf{z}(t) \right|^{2}. \]

7.2 PE after the coordinate transformation

If the original PE lower bound is \( \alpha\mathbf{I} \), then

\[ \int_{t}^{t+T} \boldsymbol{\psi}(\tau) \boldsymbol{\psi}^{T}(\tau)\,d\tau \ge \alpha_{\psi}\mathbf{I}, \qquad \alpha_{\psi} = \alpha\lambda_{\min}(\boldsymbol{\Gamma}). \]

Boundedness of the regressor gives a constant \( L>0 \) such that \( \|\boldsymbol{\psi}(t)\|\le L \).

7.3 Energy collected over one PE window

For a fixed starting time, define

\[ E_{t} = \int_{t}^{t+T} \left| \boldsymbol{\psi}^{T}(\tau)\mathbf{z}(\tau) \right|^{2}d\tau. \]

Since \( \dot{\mathbf{z} } =-\boldsymbol{\psi}\boldsymbol{\psi}^{T}\mathbf{z} \),

\[ \begin{aligned} \|\mathbf{z}(\tau)-\mathbf{z}(t)\| &\le \int_{t}^{\tau} \|\boldsymbol{\psi}(s)\| \left| \boldsymbol{\psi}^{T}(s)\mathbf{z}(s) \right|ds \\ &\le L\sqrt{T}\sqrt{E_{t} }. \end{aligned} \]

The triangle inequality in the \( L_{2}[t,t+T] \) norm yields

\[ \left( \int_{t}^{t+T} \left| \boldsymbol{\psi}^{T}(\tau)\mathbf{z}(t) \right|^{2}d\tau \right)^{1/2} \le \left(1+L^{2}T\right)\sqrt{E_{t} }. \]

On the other hand, PE gives

\[ \int_{t}^{t+T} \left| \boldsymbol{\psi}^{T}(\tau)\mathbf{z}(t) \right|^{2}d\tau \ge \alpha_{\psi}\|\mathbf{z}(t)\|^{2}. \]

Combining the two inequalities gives

\[ E_{t} \ge \frac{\alpha_{\psi} } {\left(1+L^{2}T\right)^{2} } \|\mathbf{z}(t)\|^{2}. \]

7.4 Contraction and exponential decay

Integration of \( \dot W \) over one window gives

\[ \begin{aligned} W(t+T) &= W(t)-E_{t} \\ &\le \left[ 1- \frac{2\alpha_{\psi} } {\left(1+L^{2}T\right)^{2} } \right]W(t) \\ &= qW(t), \qquad 0\le q<1. \end{aligned} \]

Iteration at \( t=t_{0}+kT \) gives

\[ W(t_{0}+kT) \le q^{k}W(t_{0}). \]

Because \( W(t) \) is nonincreasing between window endpoints and the transformation through \( \boldsymbol{\Gamma}^{1/2} \) is nonsingular, this geometric window-by-window decrease implies uniform exponential convergence of the original parameter error. A conservative window-based norm rate is proportional to

\[ \lambda = -\frac{\log q}{2T}. \]

Suppose the regressor always lies in one fixed direction:

\[ \boldsymbol{\phi}(t)=s(t)\mathbf{a}, \qquad \mathbf{a}\in\mathbb{R}^{p}. \]

Then

\[ \mathbf{G}_{\phi}(t,T) = \left( \int_{t}^{t+T}s^{2}(\tau)d\tau \right) \mathbf{a}\mathbf{a}^{T}, \]

whose rank is at most one. The output contains only the scalar parameter combination \( \mathbf{a}^{T}\boldsymbol{\theta}^{*} \). Every vector \( \mathbf{v} \) in the nullspace \( \mathbf{a}^{T}\mathbf{v}=0 \) is unidentifiable.

With \( \boldsymbol{\Gamma}=\gamma\mathbf{I} \), the gradient law moves the estimate only along \( \mathbf{a} \). Therefore, the orthogonal component of the initial parameter error remains unchanged. When enough scalar excitation is present along \( \mathbf{a} \), the limiting estimate is

\[ \hat{\boldsymbol{\theta} }_{\infty} = \hat{\boldsymbol{\theta} }(0) - \mathbf{a} \frac{ \mathbf{a}^{T} \left[ \hat{\boldsymbol{\theta} }(0)-\boldsymbol{\theta}^{*} \right] }{ \|\mathbf{a}\|^{2} }. \]

This estimate reproduces the measured output, but it need not equal the physical parameter vector.

Let

\[ \boldsymbol{\theta}^{*} = \begin{bmatrix} 1.5\\-0.5 \end{bmatrix}, \qquad \hat{\boldsymbol{\theta} }(0) = \begin{bmatrix} -1\\1 \end{bmatrix}. \]

9.1 PE regressor

\[ \boldsymbol{\phi}_{\mathrm{PE} }(t) = \begin{bmatrix} \sin t\\ \cos t \end{bmatrix}. \]

Over any interval of length \( 2\pi \),

\[ \int_{t}^{t+2\pi} \boldsymbol{\phi}_{\mathrm{PE} }(\tau) \boldsymbol{\phi}_{\mathrm{PE} }^{T}(\tau)d\tau = \pi\mathbf{I}_{2}. \]

Both parameter directions receive equal information, and the parameter error converges to zero.

9.2 Non-PE regressor

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

Its window Gramian is

\[ \int_{t}^{t+2\pi} \boldsymbol{\phi}_{\mathrm{NPE} }(\tau) \boldsymbol{\phi}_{\mathrm{NPE} }^{T}(\tau)d\tau = \pi \begin{bmatrix} 1&2\\ 2&4 \end{bmatrix}, \]

with eigenvalues \( 0 \) and \( 5\pi \). Only \( \theta_{1}+2\theta_{2} \) can be learned. The analytical limiting estimate is

\[ \hat{\boldsymbol{\theta} }_{\infty} = \begin{bmatrix} -1.1\\0.8 \end{bmatrix} \neq \begin{bmatrix} 1.5\\-0.5 \end{bmatrix}, \]

but both vectors satisfy

\[ \theta_{1}+2\theta_{2}=0.5, \]

so the final prediction error is zero despite a nonzero parameter error.

flowchart TD
  S["Choose a sliding window of length T"] --> G["Accumulate G = integral of phi phi-transpose"]
  G --> E["Compute the smallest eigenvalue"]
  E --> Q{"Is the smallest eigenvalue uniformly \nabove a positive threshold?"}
  Q -->|"Yes"| PE["All parameter directions are repeatedly visible"]
  PE --> PC["Gradient parameter error contracts \nover every window"]
  Q -->|"No"| NP["At least one direction is weak or invisible"]
  NP --> Z["Prediction error may vanish without \ntrue parameter recovery"]
        

A common MRAC error model has the form

\[ \dot{\mathbf{e} } = \mathbf{A}_{m}\mathbf{e} + \mathbf{b} \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(t), \qquad \mathbf{A}_{m}\text{ Hurwitz}. \]

Let \( \mathbf{P}\succ0 \) solve

\[ \mathbf{A}_{m}^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_{m} = -\mathbf{Q}, \qquad \mathbf{Q}\succ0. \]

With the standard cancellation-based adaptive law, a composite Lyapunov function can satisfy

\[ V = \mathbf{e}^{T}\mathbf{P}\mathbf{e} + \tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }, \qquad \dot V = -\mathbf{e}^{T}\mathbf{Q}\mathbf{e}. \]

This proves boundedness and, under standard regularity assumptions, \( \mathbf{e}(t)\to\mathbf{0} \). It does not by itself prove parameter convergence because the derivative is not negative definite in \( \tilde{\boldsymbol{\theta} } \).

Parameter convergence requires an additional identifiability argument. In many standard MRAC structures, PE of the relevant regressor, together with the stable error-model dynamics, rules out a nonzero constant parameter error that could keep producing zero tracking error. The precise regressor must be checked: PE of the command signal alone does not automatically guarantee PE of every internally generated regressor.

11.1 PE is stronger than nonzero activity

A regressor can be nonzero for all time and still be rank deficient. PE concerns directional richness over windows, not signal magnitude alone.

11.2 Adaptation gain cannot create missing information

Increasing \( \boldsymbol{\Gamma} \) can accelerate motion in excited directions, but it cannot remove the nullspace of the information matrix. Structural lack of excitation is not a tuning problem.

11.3 Redundant parameterizations prevent full recovery

If two columns of a regressor are always linearly dependent, the corresponding parameters are not separately identifiable. Rewriting the model with a minimal parameter vector may be necessary before asking for full convergence.

11.4 Noise changes the conclusion

The exact theorem above assumes a correct regression with constant parameters. With measurement noise, disturbances, or unmodeled dynamics, persistent excitation preserves information but exact convergence may be replaced by convergence to a neighborhood, statistical consistency under additional assumptions, or a biased limit.

11.5 Regulation often destroys excitation

A successful regulator can drive states and inputs toward constants, causing the regressor to lose richness. This explains why good tracking can coexist with poor identification. Practical excitation design must balance parameter learning against control performance and actuator limitations.

For sampled regressor data \( \boldsymbol{\phi}_{k} =\boldsymbol{\phi}(t_{k}) \) with sampling interval \( \Delta t \), approximate the window Gramian by

\[ \widehat{\mathbf{G} }_{\phi}(k,N) = \Delta t \sum_{i=k}^{k+N} w_{i} \boldsymbol{\phi}_{i} \boldsymbol{\phi}_{i}^{T}, \]

where trapezoidal weights use \( w_{k}=w_{k+N}=\frac{1}{2} \) and unit interior weights. A practical diagnostic is

\[ \mu_{k} = \lambda_{\min} \left( \widehat{\mathbf{G} }_{\phi}(k,N) \right). \]

A consistently positive \( \mu_{k} \) supports a PE claim for the observed experiment. A value near zero indicates a weak or invisible parameter direction. Finite data can only diagnose the tested windows; it cannot, by itself, prove an infinite-horizon condition.

A scale-independent supplementary metric is the Gramian condition number, but the minimum eigenvalue remains the direct measure of the weakest excited direction.

This implementation uses NumPy, an explicit fourth-order Runge–Kutta solver, and a sliding-window Gramian calculation. Matplotlib is optional.

Chapter10_Lesson2.py

"""
Chapter10_Lesson2.py
Persistent excitation (PE) and parameter convergence for a continuous-time
gradient estimator.

Model:
    y(t) = phi(t)^T theta_star
    d(theta_hat)/dt = -Gamma phi(t) (phi(t)^T theta_hat - y(t))

Two regressors are compared:
1. PE:     phi(t) = [sin(t), cos(t)]^T
2. Non-PE: phi(t) = [sin(t), 2 sin(t)]^T

Dependencies:
    numpy
    matplotlib (optional, only for plots)
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import numpy as np


Vector = np.ndarray
Regressor = Callable[[float], Vector]


@dataclass
class SimulationResult:
    time: Vector
    theta_hat: np.ndarray
    prediction_error: Vector
    min_window_eigenvalue: Vector


THETA_STAR = np.array([1.5, -0.5], dtype=float)
THETA_HAT_0 = np.array([-1.0, 1.0], dtype=float)
GAMMA = np.diag([3.0, 3.0])
WINDOW = 2.0 * np.pi


def phi_pe(t: float) -> Vector:
    """A two-dimensional PE regressor."""
    return np.array([np.sin(t), np.cos(t)], dtype=float)


def phi_non_pe(t: float) -> Vector:
    """A rank-one regressor: all samples are parallel to [1, 2]^T."""
    return np.array([np.sin(t), 2.0 * np.sin(t)], dtype=float)


def derivative(t: float, theta_hat: Vector, regressor: Regressor) -> Vector:
    phi = regressor(t)
    y = float(phi @ THETA_STAR)
    prediction_error = float(phi @ theta_hat - y)
    return -GAMMA @ phi * prediction_error


def rk4_step(
    t: float,
    state: Vector,
    dt: float,
    regressor: Regressor,
) -> Vector:
    k1 = derivative(t, state, regressor)
    k2 = derivative(t + 0.5 * dt, state + 0.5 * dt * k1, regressor)
    k3 = derivative(t + 0.5 * dt, state + 0.5 * dt * k2, regressor)
    k4 = derivative(t + dt, state + dt * k3, regressor)
    return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def window_gramian(
    regressor: Regressor,
    start_time: float,
    window: float,
    samples: int = 2001,
) -> np.ndarray:
    """Numerically approximate integral phi(t) phi(t)^T dt."""
    tau = np.linspace(start_time, start_time + window, samples)
    phis = np.vstack([regressor(float(t)) for t in tau])
    outer_products = np.einsum("ni,nj->nij", phis, phis)
    return np.trapezoid(outer_products, tau, axis=0)


def simulate(
    regressor: Regressor,
    final_time: float = 30.0,
    dt: float = 0.002,
    diagnostic_stride: int = 250,
) -> SimulationResult:
    steps = int(round(final_time / dt))
    time = np.linspace(0.0, final_time, steps + 1)
    theta_hat = np.zeros((steps + 1, 2), dtype=float)
    prediction_error = np.zeros(steps + 1, dtype=float)
    min_window_eigenvalue = np.full(steps + 1, np.nan, dtype=float)
    theta_hat[0] = THETA_HAT_0

    for k, t in enumerate(time[:-1]):
        phi = regressor(float(t))
        prediction_error[k] = float(phi @ (theta_hat[k] - THETA_STAR))
        theta_hat[k + 1] = rk4_step(float(t), theta_hat[k], dt, regressor)

        if k % diagnostic_stride == 0:
            gramian = window_gramian(regressor, float(t), WINDOW)
            min_window_eigenvalue[k] = float(np.linalg.eigvalsh(gramian)[0])

    final_phi = regressor(float(time[-1]))
    prediction_error[-1] = float(final_phi @ (theta_hat[-1] - THETA_STAR))
    gramian = window_gramian(regressor, float(time[-1]), WINDOW)
    min_window_eigenvalue[-1] = float(np.linalg.eigvalsh(gramian)[0])

    return SimulationResult(
        time=time,
        theta_hat=theta_hat,
        prediction_error=prediction_error,
        min_window_eigenvalue=min_window_eigenvalue,
    )


def save_csv(path: Path, result: SimulationResult) -> None:
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(
            ["time", "theta_hat_1", "theta_hat_2", "prediction_error"]
        )
        for t, theta, error in zip(
            result.time, result.theta_hat, result.prediction_error
        ):
            writer.writerow([t, theta[0], theta[1], error])


def print_report(name: str, regressor: Regressor, result: SimulationResult) -> None:
    gramian = window_gramian(regressor, 0.0, WINDOW)
    eigenvalues = np.linalg.eigvalsh(gramian)
    parameter_error = result.theta_hat[-1] - THETA_STAR

    print(f"\n{name}")
    print("-" * len(name))
    print("Window Gramian:")
    print(gramian)
    print("Window-Gramian eigenvalues:", eigenvalues)
    print("Final theta_hat:", result.theta_hat[-1])
    print("True theta:", THETA_STAR)
    print("Final parameter error norm:", np.linalg.norm(parameter_error))
    print("Final prediction error:", result.prediction_error[-1])


def create_plots(pe: SimulationResult, non_pe: SimulationResult) -> None:
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print("\nmatplotlib is not installed; CSV files were still generated.")
        return

    plt.figure()
    plt.plot(pe.time, pe.theta_hat[:, 0], label="theta_hat_1")
    plt.plot(pe.time, pe.theta_hat[:, 1], label="theta_hat_2")
    plt.axhline(THETA_STAR[0], linestyle="--", label="theta_1 true")
    plt.axhline(THETA_STAR[1], linestyle="--", label="theta_2 true")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.title("PE regressor: both parameters converge")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter10_Lesson2_PE.png", dpi=180)

    plt.figure()
    plt.plot(non_pe.time, non_pe.theta_hat[:, 0], label="theta_hat_1")
    plt.plot(non_pe.time, non_pe.theta_hat[:, 1], label="theta_hat_2")
    plt.axhline(THETA_STAR[0], linestyle="--", label="theta_1 true")
    plt.axhline(THETA_STAR[1], linestyle="--", label="theta_2 true")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.title("Non-PE regressor: only one parameter combination is learned")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter10_Lesson2_NonPE.png", dpi=180)


def main() -> None:
    pe = simulate(phi_pe)
    non_pe = simulate(phi_non_pe)

    print_report("Persistently exciting case", phi_pe, pe)
    print_report("Non-persistently exciting case", phi_non_pe, non_pe)

    save_csv(Path("Chapter10_Lesson2_PE.csv"), pe)
    save_csv(Path("Chapter10_Lesson2_NonPE.csv"), non_pe)
    create_plots(pe, non_pe)


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library, performs RK4 integration, computes the two Gramian eigenvalues analytically, and writes CSV output.

Chapter10_Lesson2.cpp

/*
Chapter10_Lesson2.cpp
Continuous-time gradient parameter estimator under PE and non-PE regressors.

Compile:
    g++ -std=c++17 -O2 Chapter10_Lesson2.cpp -o Chapter10_Lesson2
Run:
    ./Chapter10_Lesson2
*/

#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>

namespace {

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;
using Regressor = Vec2 (*)(double);

constexpr double PI = 3.14159265358979323846;
constexpr Vec2 THETA_STAR{1.5, -0.5};
constexpr Vec2 THETA_HAT_0{-1.0, 1.0};
constexpr double GAMMA = 3.0;
constexpr double WINDOW = 2.0 * PI;

Vec2 add(const Vec2& a, const Vec2& b) {
    return {a[0] + b[0], a[1] + b[1]};
}

Vec2 scale(double s, const Vec2& a) {
    return {s * a[0], s * a[1]};
}

double dot(const Vec2& a, const Vec2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

double norm(const Vec2& a) {
    return std::sqrt(dot(a, a));
}

Vec2 phiPE(double t) {
    return {std::sin(t), std::cos(t)};
}

Vec2 phiNonPE(double t) {
    return {std::sin(t), 2.0 * std::sin(t)};
}

Vec2 derivative(double t, const Vec2& thetaHat, Regressor phiFunction) {
    const Vec2 phi = phiFunction(t);
    const double y = dot(phi, THETA_STAR);
    const double predictionError = dot(phi, thetaHat) - y;
    return {
        -GAMMA * phi[0] * predictionError,
        -GAMMA * phi[1] * predictionError
    };
}

Vec2 rk4Step(double t, const Vec2& state, double dt, Regressor phiFunction) {
    const Vec2 k1 = derivative(t, state, phiFunction);
    const Vec2 k2 = derivative(
        t + 0.5 * dt,
        add(state, scale(0.5 * dt, k1)),
        phiFunction
    );
    const Vec2 k3 = derivative(
        t + 0.5 * dt,
        add(state, scale(0.5 * dt, k2)),
        phiFunction
    );
    const Vec2 k4 = derivative(
        t + dt,
        add(state, scale(dt, k3)),
        phiFunction
    );

    return add(
        state,
        scale(
            dt / 6.0,
            add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))
        )
    );
}

Mat2 windowGramian(
    Regressor phiFunction,
    double startTime,
    double window,
    int intervals = 20000
) {
    if (intervals <= 0) {
        throw std::invalid_argument("intervals must be positive");
    }

    Mat2 gramian{ { {0.0, 0.0}, {0.0, 0.0} } };
    const double h = window / static_cast<double>(intervals);

    for (int i = 0; i <= intervals; ++i) {
        const double t = startTime + i * h;
        const double weight = (i == 0 || i == intervals) ? 0.5 : 1.0;
        const Vec2 phi = phiFunction(t);

        gramian[0][0] += weight * phi[0] * phi[0] * h;
        gramian[0][1] += weight * phi[0] * phi[1] * h;
        gramian[1][0] += weight * phi[1] * phi[0] * h;
        gramian[1][1] += weight * phi[1] * phi[1] * h;
    }

    return gramian;
}

Vec2 symmetricEigenvalues(const Mat2& matrix) {
    const double trace = matrix[0][0] + matrix[1][1];
    const double difference = matrix[0][0] - matrix[1][1];
    const double discriminant = std::sqrt(
        difference * difference + 4.0 * matrix[0][1] * matrix[1][0]
    );
    return {0.5 * (trace - discriminant), 0.5 * (trace + discriminant)};
}

void simulate(
    const std::string& name,
    const std::string& csvFile,
    Regressor phiFunction,
    double finalTime = 30.0,
    double dt = 0.002
) {
    std::ofstream csv(csvFile);
    if (!csv) {
        throw std::runtime_error("Cannot open output file: " + csvFile);
    }
    csv << "time,theta_hat_1,theta_hat_2,prediction_error\n";
    csv << std::setprecision(12);

    Vec2 thetaHat = THETA_HAT_0;
    const int steps = static_cast<int>(std::llround(finalTime / dt));

    for (int k = 0; k <= steps; ++k) {
        const double t = k * dt;
        const Vec2 phi = phiFunction(t);
        const Vec2 parameterError{
            thetaHat[0] - THETA_STAR[0],
            thetaHat[1] - THETA_STAR[1]
        };
        const double predictionError = dot(phi, parameterError);

        csv << t << "," << thetaHat[0] << "," << thetaHat[1] << ","
            << predictionError << "\n";

        if (k < steps) {
            thetaHat = rk4Step(t, thetaHat, dt, phiFunction);
        }
    }

    const Mat2 gramian = windowGramian(phiFunction, 0.0, WINDOW);
    const Vec2 eigenvalues = symmetricEigenvalues(gramian);
    const Vec2 finalError{
        thetaHat[0] - THETA_STAR[0],
        thetaHat[1] - THETA_STAR[1]
    };
    const Vec2 finalPhi = phiFunction(finalTime);

    std::cout << "\n" << name << "\n";
    std::cout << std::string(name.size(), '-') << "\n";
    std::cout << "Window Gramian = [[" << gramian[0][0] << ", "
              << gramian[0][1] << "], [" << gramian[1][0] << ", "
              << gramian[1][1] << "]]\n";
    std::cout << "Gramian eigenvalues = [" << eigenvalues[0] << ", "
              << eigenvalues[1] << "]\n";
    std::cout << "Final theta_hat = [" << thetaHat[0] << ", "
              << thetaHat[1] << "]\n";
    std::cout << "Final parameter-error norm = " << norm(finalError) << "\n";
    std::cout << "Final prediction error = " << dot(finalPhi, finalError)
              << "\n";
}

}  // namespace

int main() {
    try {
        simulate(
            "Persistently exciting case",
            "Chapter10_Lesson2_PE_cpp.csv",
            phiPE
        );
        simulate(
            "Non-persistently exciting case",
            "Chapter10_Lesson2_NonPE_cpp.csv",
            phiNonPE
        );
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << "\n";
        return 1;
    }
    return 0;
}

The Java implementation mirrors the same estimator, regressor cases, RK4 integration, Gramian calculation, and CSV export.

Chapter10_Lesson2.java

/*
Chapter10_Lesson2.java
Continuous-time gradient parameter estimator under PE and non-PE regressors.

Compile:
    javac Chapter10_Lesson2.java
Run:
    java Chapter10_Lesson2
*/

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

public final class Chapter10_Lesson2 {
    private static final double PI = Math.PI;
    private static final double[] THETA_STAR = {1.5, -0.5};
    private static final double[] THETA_HAT_0 = {-1.0, 1.0};
    private static final double GAMMA = 3.0;
    private static final double WINDOW = 2.0 * PI;

    @FunctionalInterface
    private interface Regressor {
        double[] value(double time);
    }

    private Chapter10_Lesson2() {
    }

    private static double[] phiPE(double time) {
        return new double[] {Math.sin(time), Math.cos(time)};
    }

    private static double[] phiNonPE(double time) {
        return new double[] {Math.sin(time), 2.0 * Math.sin(time)};
    }

    private static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1];
    }

    private static double norm(double[] a) {
        return Math.sqrt(dot(a, a));
    }

    private static double[] add(double[] a, double[] b) {
        return new double[] {a[0] + b[0], a[1] + b[1]};
    }

    private static double[] scale(double scalar, double[] vector) {
        return new double[] {scalar * vector[0], scalar * vector[1]};
    }

    private static double[] derivative(
        double time,
        double[] thetaHat,
        Regressor regressor
    ) {
        double[] phi = regressor.value(time);
        double y = dot(phi, THETA_STAR);
        double predictionError = dot(phi, thetaHat) - y;

        return new double[] {
            -GAMMA * phi[0] * predictionError,
            -GAMMA * phi[1] * predictionError
        };
    }

    private static double[] rk4Step(
        double time,
        double[] state,
        double dt,
        Regressor regressor
    ) {
        double[] k1 = derivative(time, state, regressor);
        double[] k2 = derivative(
            time + 0.5 * dt,
            add(state, scale(0.5 * dt, k1)),
            regressor
        );
        double[] k3 = derivative(
            time + 0.5 * dt,
            add(state, scale(0.5 * dt, k2)),
            regressor
        );
        double[] k4 = derivative(
            time + dt,
            add(state, scale(dt, k3)),
            regressor
        );

        return add(
            state,
            scale(
                dt / 6.0,
                add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))
            )
        );
    }

    private static double[][] windowGramian(
        Regressor regressor,
        double startTime,
        double window,
        int intervals
    ) {
        if (intervals <= 0) {
            throw new IllegalArgumentException("intervals must be positive");
        }

        double[][] gramian = { {0.0, 0.0}, {0.0, 0.0} };
        double h = window / intervals;

        for (int i = 0; i <= intervals; i++) {
            double time = startTime + i * h;
            double weight = (i == 0 || i == intervals) ? 0.5 : 1.0;
            double[] phi = regressor.value(time);

            gramian[0][0] += weight * phi[0] * phi[0] * h;
            gramian[0][1] += weight * phi[0] * phi[1] * h;
            gramian[1][0] += weight * phi[1] * phi[0] * h;
            gramian[1][1] += weight * phi[1] * phi[1] * h;
        }

        return gramian;
    }

    private static double[] symmetricEigenvalues(double[][] matrix) {
        double trace = matrix[0][0] + matrix[1][1];
        double difference = matrix[0][0] - matrix[1][1];
        double discriminant = Math.sqrt(
            difference * difference
                + 4.0 * matrix[0][1] * matrix[1][0]
        );

        return new double[] {
            0.5 * (trace - discriminant),
            0.5 * (trace + discriminant)
        };
    }

    private static void simulate(
        String name,
        String csvFile,
        Regressor regressor,
        double finalTime,
        double dt
    ) throws IOException {
        int steps = (int) Math.round(finalTime / dt);
        double[] thetaHat = THETA_HAT_0.clone();

        try (BufferedWriter writer = Files.newBufferedWriter(Path.of(csvFile))) {
            writer.write("time,theta_hat_1,theta_hat_2,prediction_error");
            writer.newLine();

            for (int k = 0; k <= steps; k++) {
                double time = k * dt;
                double[] phi = regressor.value(time);
                double[] parameterError = {
                    thetaHat[0] - THETA_STAR[0],
                    thetaHat[1] - THETA_STAR[1]
                };
                double predictionError = dot(phi, parameterError);

                writer.write(String.format(
                    Locale.ROOT,
                    "%.12f,%.12f,%.12f,%.12f",
                    time,
                    thetaHat[0],
                    thetaHat[1],
                    predictionError
                ));
                writer.newLine();

                if (k < steps) {
                    thetaHat = rk4Step(time, thetaHat, dt, regressor);
                }
            }
        }

        double[][] gramian = windowGramian(regressor, 0.0, WINDOW, 20000);
        double[] eigenvalues = symmetricEigenvalues(gramian);
        double[] finalError = {
            thetaHat[0] - THETA_STAR[0],
            thetaHat[1] - THETA_STAR[1]
        };
        double[] finalPhi = regressor.value(finalTime);

        System.out.println();
        System.out.println(name);
        System.out.println("-".repeat(name.length()));
        System.out.printf(
            Locale.ROOT,
            "Window Gramian = [[%.8f, %.8f], [%.8f, %.8f]]%n",
            gramian[0][0],
            gramian[0][1],
            gramian[1][0],
            gramian[1][1]
        );
        System.out.printf(
            Locale.ROOT,
            "Gramian eigenvalues = [%.8f, %.8f]%n",
            eigenvalues[0],
            eigenvalues[1]
        );
        System.out.printf(
            Locale.ROOT,
            "Final theta_hat = [%.8f, %.8f]%n",
            thetaHat[0],
            thetaHat[1]
        );
        System.out.printf(
            Locale.ROOT,
            "Final parameter-error norm = %.8e%n",
            norm(finalError)
        );
        System.out.printf(
            Locale.ROOT,
            "Final prediction error = %.8e%n",
            dot(finalPhi, finalError)
        );
    }

    public static void main(String[] args) {
        try {
            simulate(
                "Persistently exciting case",
                "Chapter10_Lesson2_PE_java.csv",
                Chapter10_Lesson2::phiPE,
                30.0,
                0.002
            );
            simulate(
                "Non-persistently exciting case",
                "Chapter10_Lesson2_NonPE_java.csv",
                Chapter10_Lesson2::phiNonPE,
                30.0,
                0.002
            );
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script uses ode45, computes the information matrices, produces comparison plots, and includes a function that programmatically creates a Simulink model when Simulink is installed.

Chapter10_Lesson2.m

% Chapter10_Lesson2.m
% Persistent excitation and parameter convergence for a continuous-time
% gradient estimator. The script compares PE and non-PE regressors and can
% also construct a Simulink model programmatically.
%
% Required:
%   MATLAB
% Optional:
%   Simulink (for buildPEParameterConvergenceModel)

clear; close all; clc;

thetaStar = [1.5; -0.5];
thetaHat0 = [-1.0; 1.0];
Gamma = 3.0 * eye(2);
timeSpan = [0, 30];
window = 2*pi;

phiPE = @(t) [sin(t); cos(t)];
phiNonPE = @(t) [sin(t); 2*sin(t)];

rhsPE = @(t, thetaHat) adaptiveLaw(t, thetaHat, thetaStar, Gamma, phiPE);
rhsNonPE = @(t, thetaHat) adaptiveLaw( ...
    t, thetaHat, thetaStar, Gamma, phiNonPE);

options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[tPE, thetaPE] = ode45(rhsPE, timeSpan, thetaHat0, options);
[tNonPE, thetaNonPE] = ode45(rhsNonPE, timeSpan, thetaHat0, options);

gramianPE = windowGramian(phiPE, 0, window);
gramianNonPE = windowGramian(phiNonPE, 0, window);

fprintf('\nPersistently exciting case\n');
fprintf('Window-Gramian eigenvalues: %.10f, %.10f\n', eig(gramianPE));
fprintf('Final thetaHat: [%.10f, %.10f]^T\n', ...
    thetaPE(end, 1), thetaPE(end, 2));
fprintf('Final parameter-error norm: %.6e\n', ...
    norm(thetaPE(end, :)' - thetaStar));

fprintf('\nNon-persistently exciting case\n');
fprintf('Window-Gramian eigenvalues: %.10f, %.10f\n', ...
    eig(gramianNonPE));
fprintf('Final thetaHat: [%.10f, %.10f]^T\n', ...
    thetaNonPE(end, 1), thetaNonPE(end, 2));
fprintf('Final parameter-error norm: %.6e\n', ...
    norm(thetaNonPE(end, :)' - thetaStar));

figure('Name', 'PE parameter convergence');
plot(tPE, thetaPE(:, 1), 'LineWidth', 1.2); hold on;
plot(tPE, thetaPE(:, 2), 'LineWidth', 1.2);
yline(thetaStar(1), '--', 'theta_1 true');
yline(thetaStar(2), '--', 'theta_2 true');
xlabel('Time (s)');
ylabel('Parameter estimate');
title('PE regressor: both parameters converge');
grid on;
legend('thetaHat_1', 'thetaHat_2', 'Location', 'best');

figure('Name', 'Non-PE parameter behavior');
plot(tNonPE, thetaNonPE(:, 1), 'LineWidth', 1.2); hold on;
plot(tNonPE, thetaNonPE(:, 2), 'LineWidth', 1.2);
yline(thetaStar(1), '--', 'theta_1 true');
yline(thetaStar(2), '--', 'theta_2 true');
xlabel('Time (s)');
ylabel('Parameter estimate');
title('Non-PE regressor: only one parameter combination is learned');
grid on;
legend('thetaHat_1', 'thetaHat_2', 'Location', 'best');

peTable = table(tPE, thetaPE(:, 1), thetaPE(:, 2), ...
    'VariableNames', {'time', 'thetaHat1', 'thetaHat2'});
nonPeTable = table(tNonPE, thetaNonPE(:, 1), thetaNonPE(:, 2), ...
    'VariableNames', {'time', 'thetaHat1', 'thetaHat2'});

writetable(peTable, 'Chapter10_Lesson2_PE_matlab.csv');
writetable(nonPeTable, 'Chapter10_Lesson2_NonPE_matlab.csv');

% Set this flag to true when Simulink is installed.
buildSimulinkModel = false;
if buildSimulinkModel
    buildPEParameterConvergenceModel;
end


function dThetaHat = adaptiveLaw( ...
    time, thetaHat, thetaStar, Gamma, phiFunction)
    phi = phiFunction(time);
    measuredOutput = phi' * thetaStar;
    predictedOutput = phi' * thetaHat;
    predictionError = predictedOutput - measuredOutput;
    dThetaHat = -Gamma * phi * predictionError;
end


function gramian = windowGramian(phiFunction, startTime, window)
    time = linspace(startTime, startTime + window, 4001);
    phiSamples = zeros(2, numel(time));

    for index = 1:numel(time)
        phiSamples(:, index) = phiFunction(time(index));
    end

    gramian = zeros(2);
    gramian(1, 1) = trapz(time, phiSamples(1, :).^2);
    gramian(1, 2) = trapz( ...
        time, phiSamples(1, :) .* phiSamples(2, :));
    gramian(2, 1) = gramian(1, 2);
    gramian(2, 2) = trapz(time, phiSamples(2, :).^2);
end


function buildPEParameterConvergenceModel
    % Creates a Simulink model with a vector Integrator and a MATLAB
    % Function block. Mode = 1 selects PE; Mode = 0 selects non-PE.

    model = 'Chapter10_Lesson2_PE_Model';

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

    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Clock', [model '/Clock'], ...
        'Position', [40, 55, 70, 85]);
    add_block('simulink/Sources/Constant', [model '/Mode'], ...
        'Value', '1', ...
        'Position', [40, 140, 80, 170]);
    add_block('simulink/Continuous/Integrator', [model '/ThetaHat'], ...
        'InitialCondition', '[-1; 1]', ...
        'Position', [390, 55, 425, 105]);
    add_block( ...
        'simulink/User-Defined Functions/MATLAB Function', ...
        [model '/AdaptiveLaw'], ...
        'Position', [155, 40, 315, 185]);
    add_block('simulink/Signal Routing/Mux', [model '/Mux'], ...
        'Inputs', '3', ...
        'Position', [495, 50, 500, 155]);
    add_block('simulink/Sinks/Scope', [model '/Scope'], ...
        'Position', [555, 75, 590, 115]);
    add_block('simulink/Sinks/To Workspace', [model '/ThetaHatLog'], ...
        'VariableName', 'thetaHatSimulink', ...
        'SaveFormat', 'Structure With Time', ...
        'Position', [540, 145, 650, 175]);

    root = sfroot;
    chart = root.find( ...
        '-isa', 'Stateflow.EMChart', ...
        'Path', [model '/AdaptiveLaw']);

    chart.Script = sprintf([ ...
        'function [dThetaHat, predictionError, minEigExpected] = ' ...
        'fcn(time, thetaHat, mode)\n' ...
        'thetaStar = [1.5; -0.5];\n' ...
        'Gamma = 3.0 * eye(2);\n' ...
        'if mode >= 0.5\n' ...
        '    phi = [sin(time); cos(time)];\n' ...
        '    minEigExpected = pi;\n' ...
        'else\n' ...
        '    phi = [sin(time); 2*sin(time)];\n' ...
        '    minEigExpected = 0.0;\n' ...
        'end\n' ...
        'measuredOutput = phi'' * thetaStar;\n' ...
        'predictedOutput = phi'' * thetaHat;\n' ...
        'predictionError = predictedOutput - measuredOutput;\n' ...
        'dThetaHat = -Gamma * phi * predictionError;\n' ...
        'end\n']);

    add_line(model, 'Clock/1', 'AdaptiveLaw/1', 'autorouting', 'on');
    add_line(model, 'ThetaHat/1', 'AdaptiveLaw/2', 'autorouting', 'on');
    add_line(model, 'Mode/1', 'AdaptiveLaw/3', 'autorouting', 'on');
    add_line(model, 'AdaptiveLaw/1', 'ThetaHat/1', 'autorouting', 'on');
    add_line(model, 'ThetaHat/1', 'Mux/1', 'autorouting', 'on');
    add_line(model, 'AdaptiveLaw/2', 'Mux/2', 'autorouting', 'on');
    add_line(model, 'AdaptiveLaw/3', 'Mux/3', 'autorouting', 'on');
    add_line(model, 'Mux/1', 'Scope/1', 'autorouting', 'on');
    add_line(model, 'ThetaHat/1', 'ThetaHatLog/1', 'autorouting', 'on');

    set_param(model, ...
        'StopTime', '30', ...
        'Solver', 'ode45', ...
        'SaveTime', 'on', ...
        'SaveOutput', 'on');

    save_system(model);
    fprintf('Created %s.slx. Set Mode to 1 for PE or 0 for non-PE.\n', model);
end

The notebook uses NDSolveValue for the adaptive law and NIntegrate for the two window Gramians.

Chapter10_Lesson2.nb

Notebook[{
 Cell["Chapter 10, Lesson 2: PE as a Condition for Parameter Convergence", "Title"],
 Cell["This notebook compares a persistently exciting regressor with a rank-one non-PE regressor for a continuous-time gradient estimator.", "Text"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";",
   RowBox[{"thetaStar", "=", RowBox[{"{", RowBox[{"1.5", ",", RowBox[{"-", "0.5"}]}], "}"}]}], ";",
   RowBox[{"thetaHat0", "=", RowBox[{"{", RowBox[{RowBox[{"-", "1.0"}], ",", "1.0"}], "}"}]}], ";",
   RowBox[{"gamma", "=", "3.0"}], ";",
   RowBox[{"window", "=", RowBox[{"2", " ", "Pi"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"phiPE", "[", "t_", "]"}], ":=",
   RowBox[{"{", RowBox[{RowBox[{"Sin", "[", "t", "]"}], ",", RowBox[{"Cos", "[", "t", "]"}]}], "}"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"phiNonPE", "[", "t_", "]"}], ":=",
   RowBox[{"{", RowBox[{RowBox[{"Sin", "[", "t", "]"}], ",", RowBox[{"2", " ", RowBox[{"Sin", "[", "t", "]"}]}]}], "}"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"simulate", "[", "phi_", "]"}], ":=",
   RowBox[{"NDSolveValue", "[",
    RowBox[{
     RowBox[{"{",
      RowBox[{
       RowBox[{
        RowBox[{
         RowBox[{"thetaHat", "'"}], "[", "t", "]"}], "==",
        RowBox[{
         RowBox[{"-", "gamma"}], " ",
         RowBox[{"phi", "[", "t", "]"}], " ",
         RowBox[{"(",
          RowBox[{
           RowBox[{"phi", "[", "t", "]"}], ".",
           RowBox[{"(",
            RowBox[{
             RowBox[{"thetaHat", "[", "t", "]"}], "-", "thetaStar"}], ")"}]
          }], ")"}]
        }]}], ",",
       RowBox[{
        RowBox[{"thetaHat", "[", "0", "]"}], "==", "thetaHat0"}]
      }], "}"}], ",",
     "thetaHat", ",",
     RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}],
     ",",
     RowBox[{"Method", "->", "\"StiffnessSwitching\""}]
    }], "]"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"windowGramian", "[", "phi_", "]"}], ":=",
   RowBox[{"NIntegrate", "[",
    RowBox[{
     RowBox[{"Outer", "[",
      RowBox[{"Times", ",", RowBox[{"phi", "[", "tau", "]"}], ",", RowBox[{"phi", "[", "tau", "]"}]}],
     "]"}], ",",
     RowBox[{"{", RowBox[{"tau", ",", "0", ",", "window"}], "}"}]
    }], "]"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"solutionPE", "=", RowBox[{"simulate", "[", "phiPE", "]"}]}], ";",
   RowBox[{"solutionNonPE", "=", RowBox[{"simulate", "[", "phiNonPE", "]"}]}], ";",
   RowBox[{"gramianPE", "=", RowBox[{"windowGramian", "[", "phiPE", "]"}]}], ";",
   RowBox[{"gramianNonPE", "=", RowBox[{"windowGramian", "[", "phiNonPE", "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Column", "[",
   RowBox[{"{",
    RowBox[{
     RowBox[{"Row", "[",
      RowBox[{"{", RowBox[{"\"PE Gramian eigenvalues: \"", ",", RowBox[{"Eigenvalues", "[", "gramianPE", "]"}]}], "}"}],
     "]"}], ",",
     RowBox[{"Row", "[",
      RowBox[{"{", RowBox[{"\"Non-PE Gramian eigenvalues: \"", ",", RowBox[{"Eigenvalues", "[", "gramianNonPE", "]"}]}], "}"}],
     "]"}], ",",
     RowBox[{"Row", "[",
      RowBox[{"{", RowBox[{"\"Final PE estimate: \"", ",", RowBox[{"solutionPE", "[", "30", "]"}]}], "}"}],
     "]"}], ",",
     RowBox[{"Row", "[",
      RowBox[{"{", RowBox[{"\"Final non-PE estimate: \"", ",", RowBox[{"solutionNonPE", "[", "30", "]"}]}], "}"}],
     "]"}]
    }], "}"}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[",
   RowBox[{
    RowBox[{"Evaluate", "[", RowBox[{"solutionPE", "[", "t", "]"}], "]"}], ",",
    RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}], ",",
    RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"thetaHat1\"", ",", "\"thetaHat2\""}], "}"}]}], ",",
    RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"time\"", ",", "\"estimate\""}], "}"}]}], ",",
    RowBox[{"PlotLabel", "->", "\"PE: both estimates converge\""}]
   }], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[",
   RowBox[{
    RowBox[{"Evaluate", "[", RowBox[{"solutionNonPE", "[", "t", "]"}], "]"}], ",",
    RowBox[{"{", RowBox[{"t", ",", "0", ",", "30"}], "}"}], ",",
    RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"thetaHat1\"", ",", "\"thetaHat2\""}], "}"}]}], ",",
    RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"time\"", ",", "\"estimate\""}], "}"}]}], ",",
    RowBox[{"PlotLabel", "->", "\"Non-PE: an unidentifiable direction remains\""}]
   }], "]"}]], "Input"],
 Cell["For the PE case, the window Gramian is positive definite. For the non-PE case, it has one zero eigenvalue, so only one linear combination of the two parameters can be identified.", "Text"]
},
WindowTitle -> "Chapter10_Lesson2",
StyleDefinitions -> "Default.nb"
]

For the PE regressor, the numerical Gramian eigenvalues should both be approximately \( \pi \). The two parameter estimates approach \( 1.5 \) and \( -0.5 \), respectively.

For the non-PE regressor, the eigenvalues should be approximately \( 0 \) and \( 5\pi \). The estimates approach approximately \( -1.1 \) and \( 0.8 \). The parameter-error norm remains nonzero, while the prediction error approaches zero.

This side-by-side result is the central lesson: error convergence demonstrates agreement with the observed trajectory; PE is what makes that trajectory informative enough to select the unique true parameter vector.

Problem 1 (Exact PE Gramian): For \( \boldsymbol{\phi}(t) =[\sin(\omega t),\cos(\omega t)]^{T} \), where \( \omega\neq0 \), prove that the regressor is PE and find a valid window and lower bound.

Solution: Choose

\[ T=\frac{2\pi}{|\omega|}. \]

Over one complete period, the squared sine and cosine each integrate to half the period, and the cross term integrates to zero:

\[ \begin{aligned} \int_{t}^{t+T}\sin^{2}(\omega\tau)d\tau &=\frac{\pi}{|\omega|},\\ \int_{t}^{t+T}\cos^{2}(\omega\tau)d\tau &=\frac{\pi}{|\omega|},\\ \int_{t}^{t+T}\sin(\omega\tau)\cos(\omega\tau)d\tau &=0. \end{aligned} \]

Therefore,

\[ \mathbf{G}_{\phi}(t,T) = \frac{\pi}{|\omega|}\mathbf{I}_{2}. \]

A valid PE lower and upper bound is \( \alpha=\beta=\pi/|\omega| \).

Problem 2 (Rank-One Non-PE Regressor): Let \( \boldsymbol{\phi}(t)=s(t)\mathbf{a} \) with \( \mathbf{a}\neq\mathbf{0} \) and \( p>1 \). Show that this regressor cannot be PE of order \( p \).

Solution: The Gramian is

\[ \mathbf{G}_{\phi}(t,T) = \left( \int_{t}^{t+T}s^{2}(\tau)d\tau \right) \mathbf{a}\mathbf{a}^{T}. \]

Since \( \mathbf{a}\mathbf{a}^{T} \) has rank one, the Gramian has at least \( p-1 \) zero eigenvalues. Choose any nonzero \( \mathbf{v} \) satisfying \( \mathbf{a}^{T}\mathbf{v}=0 \). Then

\[ \mathbf{v}^{T}\mathbf{G}_{\phi}(t,T)\mathbf{v}=0, \]

which contradicts the PE lower bound \( \alpha\|\mathbf{v}\|^{2}>0 \).

Problem 3 (Why Barbalat Is Insufficient): The gradient estimator satisfies \( \dot V=-\varepsilon^{2} \). Assume the standard conditions needed for Barbalat's lemma. State exactly what follows and explain why parameter convergence does not follow without PE.

Solution: Integration gives

\[ \int_{0}^{\infty}\varepsilon^{2}(t)dt \le V(0)<\infty. \]

Uniform continuity then gives \( \varepsilon(t)\to0 \). But \( \varepsilon=\boldsymbol{\phi}^{T} \tilde{\boldsymbol{\theta} } \) is one scalar projection. If the regressor does not span the entire parameter space over time, a nonzero parameter error can remain in the common nullspace of those projections. PE is the condition that excludes such a persistent nullspace uniformly.

Problem 4 (PE Under a Constant Coordinate Change): Suppose \( \boldsymbol{\phi}(t) \) is PE and \( \mathbf{H} \) is a constant nonsingular matrix. Show that \( \boldsymbol{\psi}(t)=\mathbf{H}\boldsymbol{\phi}(t) \) is also PE.

Solution:

\[ \begin{aligned} \int_{t}^{t+T} \boldsymbol{\psi}\boldsymbol{\psi}^{T}d\tau &= \mathbf{H} \left( \int_{t}^{t+T} \boldsymbol{\phi}\boldsymbol{\phi}^{T}d\tau \right) \mathbf{H}^{T}\\ &\ge \alpha\mathbf{H}\mathbf{H}^{T}\\ &\ge \alpha\sigma_{\min}^{2}(\mathbf{H})\mathbf{I}. \end{aligned} \]

Thus, \( \boldsymbol{\psi} \) is PE with the same window and lower bound \( \alpha_{\psi} =\alpha\sigma_{\min}^{2}(\mathbf{H}) \).

Problem 5 (Analytical Limit in the Non-PE Example): For \( \mathbf{a}=[1,2]^{T} \), \( \boldsymbol{\theta}^{*}=[1.5,-0.5]^{T} \), and \( \hat{\boldsymbol{\theta} }(0)=[-1,1]^{T} \), compute the limiting estimate under a scalar-gain gradient law.

Solution: First compute

\[ \mathbf{a}^{T} \left[ \hat{\boldsymbol{\theta} }(0)-\boldsymbol{\theta}^{*} \right] = \begin{bmatrix}1&2\end{bmatrix} \begin{bmatrix}-2.5\\1.5\end{bmatrix} =0.5, \qquad \|\mathbf{a}\|^{2}=5. \]

Hence,

\[ \hat{\boldsymbol{\theta} }_{\infty} = \begin{bmatrix}-1\\1\end{bmatrix} - \begin{bmatrix}1\\2\end{bmatrix} \frac{0.5}{5} = \begin{bmatrix}-1.1\\0.8\end{bmatrix}. \]

The final estimate is incorrect, but \( [-1.1]+2[0.8]=0.5 =1.5+2[-0.5] \), so it produces the same output for this rank-one regressor.

Problem 6 (Finite-Data Excitation Test): A sampled two-parameter experiment gives the approximate window Gramian

\[ \widehat{\mathbf{G} } = \begin{bmatrix} 4.0&1.2\\ 1.2&0.5 \end{bmatrix}. \]

Compute its eigenvalues and state what the smaller eigenvalue says about the weakest parameter direction.

Solution: For a symmetric two-by-two matrix,

\[ \lambda_{\pm} = \frac{ \operatorname{tr}(\widehat{\mathbf{G} }) \pm \sqrt{ \left( \widehat{G}_{11}-\widehat{G}_{22} \right)^{2} +4\widehat{G}_{12}^{2} } }{2}. \]

Substitution gives

\[ \lambda_{\min}\approx0.128, \qquad \lambda_{\max}\approx4.372. \]

The matrix is positive definite, so both directions are represented in this window. However, the small minimum eigenvalue and the ratio \( \lambda_{\max}/\lambda_{\min}\approx34.1 \) indicate strongly unequal excitation and potentially slow or noise-sensitive convergence in the weak direction.

The standard adaptive Lyapunov argument guarantees bounded estimates and drives a tracking or prediction error toward zero, but it does not automatically identify the true parameters. Persistent excitation adds the missing information condition: every parameter direction must be observed with a uniform positive energy over every sufficiently long window.

For the gradient parameter-error dynamics, PE converts negative semidefinite Lyapunov decay into a strict window-by-window contraction, yielding uniform exponential parameter convergence. Without PE, the estimator can converge to an equivalence set of parameters that reproduce the observed output while retaining a nonzero error in unexcited directions.