After completing this lesson, students should be able to:

• distinguish convergence of the tracking error \( e(t)=x(t)-x_m(t) \) from convergence of the parameter error \( \tilde{\theta}(t)=\hat{\theta}(t)-\theta^* \);
• prove tracking-error convergence with a state-parameter Lyapunov function;
• identify the logical gap between \( \dot V\le0 \) and \( \tilde{\theta}\to0 \);
• use the excitation Gramian to determine identifiable directions; and
• explain why PE can yield simultaneous state and parameter convergence.

\[ e(t)\to0 \quad \not\Rightarrow \quad \tilde{\theta}(t)\to0. \]

A controller only needs the current estimate to generate the correct input-output behavior along the observed trajectory. Several parameter vectors may generate the same behavior when the regressor does not explore all independent directions.

flowchart TD
  A["Adaptive closed loop"] --> B["Lyapunov derivative is negative semidefinite"]
  B --> C["All error signals remain bounded"]
  C --> D["Tracking error tends to zero"]
  D --> E{"Is the regressor \npersistently exciting?"}
  E -->|No| F["Only excited parameter \ncombinations are learned"]
  F --> G["Parameter error may approach \na nonzero invariant set"]
  E -->|Yes| H["Every parameter direction \nis repeatedly excited"]
  H --> I["Parameter error tends to zero"]
        

Consider the linearly parameterized scalar error model developed in the preceding MRAC lessons:

\[ \dot e(t)=-a e(t)+\tilde{\theta}^T(t)\phi(t), \qquad a>0, \]

where \( e\in\mathbb{R} \) is a tracking error, \( \phi(t)\in\mathbb{R}^p \) is the measurable regressor, \( \theta^* \) is the constant ideal parameter, and

\[ \tilde{\theta}(t)=\hat{\theta}(t)-\theta^*. \]

Use the gradient adaptive law

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

Since \( \theta^* \) is constant, \( \dot{\tilde{\theta} }=\dot{\hat{\theta} } \). The combined linear time-varying system is

\[ \begin{bmatrix} \dot e \\ \dot{\tilde{\theta} } \end{bmatrix} = \begin{bmatrix} -a & \phi^T(t) \\ -\Gamma\phi(t) & 0 \end{bmatrix} \begin{bmatrix} e \\ \tilde{\theta} \end{bmatrix}. \]

Choose the composite Lyapunov function

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

Differentiating along the trajectories gives

\[ \begin{aligned} \dot V &= e\dot e+\tilde{\theta}^T\Gamma^{-1}\dot{\tilde{\theta} } \\ &= e\left(-ae+\tilde{\theta}^T\phi\right) + \tilde{\theta}^T\Gamma^{-1}\left(-\Gamma\phi e\right) \\ &= -ae^2 + e\tilde{\theta}^T\phi - \tilde{\theta}^T\phi e \\ &= -ae^2\le0. \end{aligned} \]

Hence \(V(t)\le V(0)\), so \(e(t)\) and \(\tilde{\theta}(t)\) are bounded. Furthermore,

\[ a\int_0^t e^2(\tau)\,d\tau = V(0)-V(t) \le V(0), \]

which yields \(e\in L_2\cap L_\infty\). If \(\phi(t)\) is bounded, the error equation makes \(\dot e\) bounded, and therefore \(e(t)\) is uniformly continuous. Barbalat's lemma gives

\[ \boxed{\lim_{t\to\infty}e(t)=0}. \]

In MRAC this means \(x(t)-x_m(t)\to0\). It does not necessarily mean \(x(t)\to0\); the plant follows the reference-model trajectory.

The derivative \(\dot V=-ae^2\) contains no strictly negative term in \(\tilde{\theta}\). It proves parameter boundedness, not decay. The zero-derivative set is

\[ \mathcal{S}=\left\{(e,\tilde{\theta}):e=0\right\}. \]

To remain in this set, the error derivative must also be zero:

\[ 0=\dot e=\tilde{\theta}^T\phi(t). \]

Therefore, the largest invariant subset of \(\mathcal{S}\) is characterized by

\[ \mathcal{M} = \left\{ (0,\tilde{\theta}): \tilde{\theta}^T\phi(t)=0 \text{ for all future }t \right\}. \]

If the regressor remains in a proper subspace \(\mathcal{R}\subsetneq\mathbb{R}^p\), every nonzero parameter error in \(\mathcal{R}^{\perp}\) is invisible:

\[ \tilde{\theta}\in\mathcal{R}^{\perp} \quad\Rightarrow\quad \tilde{\theta}^T\phi(t)=0. \]

The controller can consequently achieve zero tracking error with an incorrect parameter vector. This is lack of identifiability, not instability.

For a window of length \(T\), define

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

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

\[ v^TG(t,T)v = \int_t^{t+T}\left(v^T\phi(\tau)\right)^2d\tau. \]

Thus, a zero value means no information in direction \(v\); a small eigenvalue means weak information; and a uniformly positive smallest eigenvalue means every parameter direction is repeatedly visible.

Persistent excitation requires constants \(T>0\) and \(\alpha>0\) such that

\[ \int_t^{t+T}\phi(\tau)\phi^T(\tau)\,d\tau \ge\alpha I \qquad \text{for every }t\ge0. \]

flowchart TD
  A["Regressor samples over a moving window"] --> B["Compute G = integral(phi phi^T)"]
  B --> C{"Is the smallest eigenvalue \nuniformly positive?"}
  C -->|No| D["Rank-deficient information"]
  D --> E["Invisible parameter directions remain"]
  C -->|Yes| F["Full directional information"]
  F --> G["Unique parameter vector is identifiable"]
        

From the error dynamics,

\[ \tilde{\theta}^T(t)\phi(t)=\dot e(t)+ae(t). \]

Under the boundedness and regularity assumptions used above, \(e(t)\to0\) and \(\dot e(t)\to0\), so

\[ \tilde{\theta}^T(t)\phi(t)\to0. \]

This is only a projection result. PE turns it into a vector result because

\[ \int_t^{t+T} \left(\tilde{\theta}^T(t)\phi(\tau)\right)^2d\tau = \tilde{\theta}^T(t)G(t,T)\tilde{\theta}(t) \ge \alpha\left\|\tilde{\theta}(t)\right\|^2. \]

Moreover, \(\dot{\tilde{\theta} }=-\Gamma\phi e\to0\), so the parameter error varies increasingly slowly over a fixed window. Its projected energy can approach zero only when its norm approaches zero.

Theorem (simultaneous convergence under PE).

Suppose \(\phi\) and \(\dot\phi\) are bounded, the regressor is persistently exciting, \(a>0\), and \(\Gamma>0\). Then the origin of the combined error system is uniformly asymptotically stable; under the standard uniform-PE hypotheses it is exponentially stable. Consequently,

\[ e(t)\to0, \qquad \tilde{\theta}(t)\to0. \]

Let

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

7.1 Non-PE regressor

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

Every sample lies on one line, so the Gramian has rank one. The observed behavior depends only on \(\theta_1+2\theta_2\). Since \(2+2(-0.5)=1\), every estimate satisfying

\[ \hat{\theta}_1+2\hat{\theta}_2=1 \]

is indistinguishable in this experiment. With a zero initial estimate, the gradient update approaches the minimum-norm point on this line:

\[ \hat{\theta}_{\infty} = \frac{1}{1^2+2^2} \begin{bmatrix} 1\\ 2 \end{bmatrix} = \begin{bmatrix} 0.2\\ 0.4 \end{bmatrix}. \]

The tracking error converges to zero, but

\[ \tilde{\theta}_{\infty} = \begin{bmatrix} -1.8\\ 0.9 \end{bmatrix}, \qquad \left\|\tilde{\theta}_{\infty}\right\| \approx2.0125. \]

7.2 PE regressor

\[ \phi_{PE}(t) = \begin{bmatrix} \sin t\\ \cos(1.7t) \end{bmatrix}. \]

Independent frequency content repeatedly spans two directions. The same adaptive law now drives both \(e(t)\to0\) and \(\hat{\theta}(t)\to\theta^*\).

1. Tracking convergence: \(e=x-x_m\to0\).
2. Plant-state convergence: \(x\to0\), which also requires \(x_m\to0\).
3. Parameter convergence: \(\hat{\theta}\to\theta^*\), which requires identifiability and adequate excitation.

Exact identification matters when parameters will be interpreted physically, transferred to another operating condition, used for fault detection, or inserted into an indirect controller. For pure tracking, exact physical identification may be unnecessary.

A simulation should not infer PE merely because an estimate settles. Compute the moving-window Gramian and monitor

\[ \lambda_{\min}\left(G(t,T)\right). \]

• tracking metric: \(|e(t)|\);
• parameter metric: \(\|\tilde{\theta}(t)\|\) when the true value is known;
• excitation metric: the smallest moving-window Gramian eigenvalue;
• prediction metric: \(\tilde{\theta}^T\phi\), which can vanish even when the full parameter error does not.

NumPy supplies vector operations and Gramian eigenvalues, Matplotlib creates the convergence plots, and an explicit RK4 method keeps the adaptive dynamics transparent.

Chapter10_Lesson4.py

"""
Chapter10_Lesson4.py
State convergence versus parameter convergence in an adaptive error model.

The simulated system is
    e_dot         = -a e + tilde_theta^T phi(t)
    hat_theta_dot = -gamma phi(t) e
where tilde_theta = hat_theta - theta_star.

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

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    a: float = 2.0
    gamma: float = 5.0
    theta_star_1: float = 2.0
    theta_star_2: float = -0.5
    t_final: float = 40.0
    dt: float = 0.002

    @property
    def theta_star(self) -> np.ndarray:
        return np.array([self.theta_star_1, self.theta_star_2], dtype=float)


@dataclass
class SimulationResult:
    name: str
    t: np.ndarray
    z: np.ndarray
    min_gramian_eigenvalue: float

    @property
    def error(self) -> np.ndarray:
        return self.z[:, 0]

    @property
    def theta_hat(self) -> np.ndarray:
        return self.z[:, 1:3]


def regressor(t: float, persistent_excitation: bool) -> np.ndarray:
    """Return the two-dimensional regressor."""
    if persistent_excitation:
        return np.array([np.sin(t), np.cos(1.7 * t)], dtype=float)
    return np.array([np.sin(t), 2.0 * np.sin(t)], dtype=float)


def dynamics(
    t: float,
    z: np.ndarray,
    p: Parameters,
    persistent_excitation: bool,
) -> np.ndarray:
    """Evaluate [e_dot, theta_hat_1_dot, theta_hat_2_dot]."""
    e = z[0]
    theta_hat = z[1:3]
    phi = regressor(t, persistent_excitation)
    theta_tilde = theta_hat - p.theta_star

    e_dot = -p.a * e + float(theta_tilde @ phi)
    theta_hat_dot = -p.gamma * phi * e
    return np.concatenate(([e_dot], theta_hat_dot))


def rk4_step(
    rhs: Callable[[float, np.ndarray], np.ndarray],
    t: float,
    z: np.ndarray,
    dt: float,
) -> np.ndarray:
    """One classical fourth-order Runge-Kutta step."""
    k1 = rhs(t, z)
    k2 = rhs(t + 0.5 * dt, z + 0.5 * dt * k1)
    k3 = rhs(t + 0.5 * dt, z + 0.5 * dt * k2)
    k4 = rhs(t + dt, z + dt * k3)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def minimum_window_gramian_eigenvalue(
    t: np.ndarray,
    persistent_excitation: bool,
    window_length: float = 8.0,
    stride: int = 250,
) -> float:
    """Numerically estimate min_t lambda_min(int_t^(t+T) phi phi^T d tau)."""
    dt = float(t[1] - t[0])
    window_samples = int(round(window_length / dt))
    minimum = np.inf

    for start in range(0, len(t) - window_samples, stride):
        gramian = np.zeros((2, 2), dtype=float)
        for k in range(start, start + window_samples):
            phi = regressor(float(t[k]), persistent_excitation)
            gramian += np.outer(phi, phi) * dt
        minimum = min(minimum, float(np.linalg.eigvalsh(gramian)[0]))

    return float(minimum)


def simulate(
    name: str,
    persistent_excitation: bool,
    p: Parameters,
) -> SimulationResult:
    """Run one experiment."""
    sample_count = int(round(p.t_final / p.dt)) + 1
    t = np.linspace(0.0, p.t_final, sample_count)
    z = np.zeros((sample_count, 3), dtype=float)
    z[0] = np.array([1.0, 0.0, 0.0], dtype=float)

    rhs = lambda time, state: dynamics(time, state, p, persistent_excitation)
    for k in range(sample_count - 1):
        z[k + 1] = rk4_step(rhs, float(t[k]), z[k], p.dt)

    min_eigenvalue = minimum_window_gramian_eigenvalue(
        t, persistent_excitation
    )
    return SimulationResult(name, t, z, min_eigenvalue)


def save_csv(result: SimulationResult, p: Parameters) -> Path:
    """Write simulation data to a CSV file."""
    theta_error = result.theta_hat - p.theta_star
    data = np.column_stack(
        (
            result.t,
            result.error,
            result.theta_hat,
            theta_error,
            np.linalg.norm(theta_error, axis=1),
        )
    )
    filename = Path(f"Chapter10_Lesson4_{result.name}.csv")
    header = (
        "time,e,theta_hat_1,theta_hat_2,"
        "theta_tilde_1,theta_tilde_2,theta_tilde_norm"
    )
    np.savetxt(filename, data, delimiter=",", header=header, comments="")
    return filename


def report(result: SimulationResult, p: Parameters) -> None:
    """Print final convergence diagnostics."""
    final_theta = result.theta_hat[-1]
    final_theta_error = final_theta - p.theta_star
    print(f"\nCase: {result.name}")
    print(f"  final tracking error       = {result.error[-1]: .6e}")
    print(
        "  final parameter estimate  = "
        f"[{final_theta[0]: .6f}, {final_theta[1]: .6f}]"
    )
    print(f"  final parameter-error norm = {np.linalg.norm(final_theta_error): .6e}")
    print(
        "  minimum window Gramian eigenvalue "
        f"= {result.min_gramian_eigenvalue: .6e}"
    )


def plot_results(
    non_pe: SimulationResult,
    pe: SimulationResult,
    p: Parameters,
) -> None:
    """Create separate figures for tracking and parameter errors."""
    plt.figure()
    plt.plot(non_pe.t, non_pe.error, label="non-PE")
    plt.plot(pe.t, pe.error, label="PE")
    plt.xlabel("Time (s)")
    plt.ylabel("Tracking error e(t)")
    plt.title("State/tracking-error convergence")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter10_Lesson4_tracking_error.png", dpi=180)

    plt.figure()
    non_pe_norm = np.linalg.norm(non_pe.theta_hat - p.theta_star, axis=1)
    pe_norm = np.linalg.norm(pe.theta_hat - p.theta_star, axis=1)
    plt.semilogy(non_pe.t, np.maximum(non_pe_norm, 1e-14), label="non-PE")
    plt.semilogy(pe.t, np.maximum(pe_norm, 1e-14), label="PE")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter-error norm")
    plt.title("Parameter convergence requires sufficient excitation")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter10_Lesson4_parameter_error.png", dpi=180)

    plt.show()


def main() -> None:
    p = Parameters()
    non_pe = simulate("non_pe", False, p)
    pe = simulate("pe", True, p)

    report(non_pe, p)
    report(pe, p)
    save_csv(non_pe, p)
    save_csv(pe, p)
    plot_results(non_pe, pe, p)


if __name__ == "__main__":
    main()
      

The C++17 version uses no third-party package. It writes CSV files and computes the eigenvalues of each symmetric two-by-two Gramian analytically.

Chapter10_Lesson4.cpp

/*
Chapter10_Lesson4.cpp
State convergence versus parameter convergence in an adaptive error model.

Build:
    g++ -std=c++17 -O2 Chapter10_Lesson4.cpp -o Chapter10_Lesson4
Run:
    ./Chapter10_Lesson4
*/

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

struct Parameters {
    double a = 2.0;
    double gamma = 5.0;
    std::array<double, 2> thetaStar{2.0, -0.5};
    double tFinal = 40.0;
    double dt = 0.002;
};

struct State {
    double e = 0.0;
    double theta1 = 0.0;
    double theta2 = 0.0;
};

State operator+(const State& x, const State& y) {
    return {x.e + y.e, x.theta1 + y.theta1, x.theta2 + y.theta2};
}

State operator*(double c, const State& x) {
    return {c * x.e, c * x.theta1, c * x.theta2};
}

std::array<double, 2> regressor(double t, bool persistentExcitation) {
    if (persistentExcitation) {
        return {std::sin(t), std::cos(1.7 * t)};
    }
    return {std::sin(t), 2.0 * std::sin(t)};
}

State derivative(
    double t,
    const State& z,
    const Parameters& p,
    bool persistentExcitation
) {
    const auto phi = regressor(t, persistentExcitation);
    const double thetaTilde1 = z.theta1 - p.thetaStar[0];
    const double thetaTilde2 = z.theta2 - p.thetaStar[1];

    State dz;
    dz.e = -p.a * z.e
         + thetaTilde1 * phi[0]
         + thetaTilde2 * phi[1];
    dz.theta1 = -p.gamma * phi[0] * z.e;
    dz.theta2 = -p.gamma * phi[1] * z.e;
    return dz;
}

State rk4Step(
    double t,
    const State& z,
    const Parameters& p,
    bool persistentExcitation
) {
    const double h = p.dt;
    const State k1 = derivative(t, z, p, persistentExcitation);
    const State k2 = derivative(
        t + 0.5 * h, z + (0.5 * h) * k1, p, persistentExcitation
    );
    const State k3 = derivative(
        t + 0.5 * h, z + (0.5 * h) * k2, p, persistentExcitation
    );
    const State k4 = derivative(
        t + h, z + h * k3, p, persistentExcitation
    );

    return z + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

double smallestEigenvalue2x2(
    double a11,
    double a12,
    double a22
) {
    const double trace = a11 + a22;
    const double discriminant =
        std::sqrt((a11 - a22) * (a11 - a22) + 4.0 * a12 * a12);
    return 0.5 * (trace - discriminant);
}

double minimumWindowGramianEigenvalue(
    bool persistentExcitation,
    const Parameters& p,
    double windowLength = 8.0,
    int stride = 250
) {
    const int totalSteps = static_cast<int>(std::lround(p.tFinal / p.dt));
    const int windowSteps = static_cast<int>(std::lround(windowLength / p.dt));
    double minimum = std::numeric_limits<double>::infinity();

    for (int start = 0; start + windowSteps <= totalSteps; start += stride) {
        double g11 = 0.0;
        double g12 = 0.0;
        double g22 = 0.0;

        for (int k = start; k < start + windowSteps; ++k) {
            const double t = k * p.dt;
            const auto phi = regressor(t, persistentExcitation);
            g11 += phi[0] * phi[0] * p.dt;
            g12 += phi[0] * phi[1] * p.dt;
            g22 += phi[1] * phi[1] * p.dt;
        }
        minimum = std::min(minimum, smallestEigenvalue2x2(g11, g12, g22));
    }
    return minimum;
}

void runCase(
    const std::string& name,
    bool persistentExcitation,
    const Parameters& p
) {
    const int steps = static_cast<int>(std::lround(p.tFinal / p.dt));
    State z{1.0, 0.0, 0.0};

    std::ofstream csv("Chapter10_Lesson4_" + name + ".csv");
    if (!csv) {
        throw std::runtime_error("Could not create CSV output.");
    }
    csv << "time,e,theta_hat_1,theta_hat_2,"
           "theta_tilde_1,theta_tilde_2,theta_tilde_norm\n";
    csv << std::setprecision(12);

    for (int k = 0; k <= steps; ++k) {
        const double t = k * p.dt;
        const double d1 = z.theta1 - p.thetaStar[0];
        const double d2 = z.theta2 - p.thetaStar[1];
        csv << t << ',' << z.e << ',' << z.theta1 << ',' << z.theta2 << ','
            << d1 << ',' << d2 << ',' << std::hypot(d1, d2) << '\n';

        if (k < steps) {
            z = rk4Step(t, z, p, persistentExcitation);
        }
    }

    const double d1 = z.theta1 - p.thetaStar[0];
    const double d2 = z.theta2 - p.thetaStar[1];
    const double gramianMinimum =
        minimumWindowGramianEigenvalue(persistentExcitation, p);

    std::cout << "\nCase: " << name << '\n'
              << "  final tracking error       = " << std::scientific
              << z.e << '\n'
              << "  final parameter estimate  = [" << std::fixed
              << z.theta1 << ", " << z.theta2 << "]\n"
              << "  final parameter-error norm = " << std::scientific
              << std::hypot(d1, d2) << '\n'
              << "  minimum Gramian eigenvalue = "
              << gramianMinimum << '\n';
}

int main() {
    try {
        const Parameters p;
        runCase("non_pe", false, p);
        runCase("pe", true, p);
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}
      

The Java implementation mirrors the C++ simulation and writes portable CSV output without a third-party numerical dependency.

Chapter10_Lesson4.java

/*
Chapter10_Lesson4.java
State convergence versus parameter convergence in an adaptive error model.

Build:
    javac Chapter10_Lesson4.java
Run:
    java Chapter10_Lesson4
*/

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_Lesson4 {
    private Chapter10_Lesson4() {
    }

    private static final class Parameters {
        final double a = 2.0;
        final double gamma = 5.0;
        final double[] thetaStar = {2.0, -0.5};
        final double tFinal = 40.0;
        final double dt = 0.002;
    }

    private static final class State {
        final double e;
        final double theta1;
        final double theta2;

        State(double e, double theta1, double theta2) {
            this.e = e;
            this.theta1 = theta1;
            this.theta2 = theta2;
        }

        State add(State other) {
            return new State(
                e + other.e,
                theta1 + other.theta1,
                theta2 + other.theta2
            );
        }

        State scale(double value) {
            return new State(
                value * e,
                value * theta1,
                value * theta2
            );
        }
    }

    private static double[] regressor(
        double t,
        boolean persistentExcitation
    ) {
        if (persistentExcitation) {
            return new double[]{Math.sin(t), Math.cos(1.7 * t)};
        }
        return new double[]{Math.sin(t), 2.0 * Math.sin(t)};
    }

    private static State derivative(
        double t,
        State z,
        Parameters p,
        boolean persistentExcitation
    ) {
        double[] phi = regressor(t, persistentExcitation);
        double thetaTilde1 = z.theta1 - p.thetaStar[0];
        double thetaTilde2 = z.theta2 - p.thetaStar[1];

        double eDot = -p.a * z.e
            + thetaTilde1 * phi[0]
            + thetaTilde2 * phi[1];
        double theta1Dot = -p.gamma * phi[0] * z.e;
        double theta2Dot = -p.gamma * phi[1] * z.e;
        return new State(eDot, theta1Dot, theta2Dot);
    }

    private static State rk4Step(
        double t,
        State z,
        Parameters p,
        boolean persistentExcitation
    ) {
        double h = p.dt;
        State k1 = derivative(t, z, p, persistentExcitation);
        State k2 = derivative(
            t + 0.5 * h,
            z.add(k1.scale(0.5 * h)),
            p,
            persistentExcitation
        );
        State k3 = derivative(
            t + 0.5 * h,
            z.add(k2.scale(0.5 * h)),
            p,
            persistentExcitation
        );
        State k4 = derivative(
            t + h,
            z.add(k3.scale(h)),
            p,
            persistentExcitation
        );

        return z.add(
            k1.add(k2.scale(2.0))
              .add(k3.scale(2.0))
              .add(k4)
              .scale(h / 6.0)
        );
    }

    private static double smallestEigenvalue2x2(
        double a11,
        double a12,
        double a22
    ) {
        double trace = a11 + a22;
        double discriminant = Math.sqrt(
            (a11 - a22) * (a11 - a22) + 4.0 * a12 * a12
        );
        return 0.5 * (trace - discriminant);
    }

    private static double minimumWindowGramianEigenvalue(
        boolean persistentExcitation,
        Parameters p,
        double windowLength,
        int stride
    ) {
        int totalSteps = (int) Math.round(p.tFinal / p.dt);
        int windowSteps = (int) Math.round(windowLength / p.dt);
        double minimum = Double.POSITIVE_INFINITY;

        for (
            int start = 0;
            start + windowSteps <= totalSteps;
            start += stride
        ) {
            double g11 = 0.0;
            double g12 = 0.0;
            double g22 = 0.0;

            for (int k = start; k < start + windowSteps; k++) {
                double t = k * p.dt;
                double[] phi = regressor(t, persistentExcitation);
                g11 += phi[0] * phi[0] * p.dt;
                g12 += phi[0] * phi[1] * p.dt;
                g22 += phi[1] * phi[1] * p.dt;
            }
            minimum = Math.min(
                minimum,
                smallestEigenvalue2x2(g11, g12, g22)
            );
        }
        return minimum;
    }

    private static void runCase(
        String name,
        boolean persistentExcitation,
        Parameters p
    ) throws IOException {
        int steps = (int) Math.round(p.tFinal / p.dt);
        State z = new State(1.0, 0.0, 0.0);
        Path csvPath = Path.of("Chapter10_Lesson4_" + name + ".csv");

        try (BufferedWriter writer = Files.newBufferedWriter(csvPath)) {
            writer.write(
                "time,e,theta_hat_1,theta_hat_2,"
                + "theta_tilde_1,theta_tilde_2,theta_tilde_norm\n"
            );

            for (int k = 0; k <= steps; k++) {
                double t = k * p.dt;
                double d1 = z.theta1 - p.thetaStar[0];
                double d2 = z.theta2 - p.thetaStar[1];
                writer.write(String.format(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t,
                    z.e,
                    z.theta1,
                    z.theta2,
                    d1,
                    d2,
                    Math.hypot(d1, d2)
                ));

                if (k < steps) {
                    z = rk4Step(t, z, p, persistentExcitation);
                }
            }
        }

        double d1 = z.theta1 - p.thetaStar[0];
        double d2 = z.theta2 - p.thetaStar[1];
        double gramianMinimum = minimumWindowGramianEigenvalue(
            persistentExcitation, p, 8.0, 250
        );

        System.out.printf(Locale.US, "%nCase: %s%n", name);
        System.out.printf(
            Locale.US,
            "  final tracking error       = %.6e%n",
            z.e
        );
        System.out.printf(
            Locale.US,
            "  final parameter estimate  = [%.6f, %.6f]%n",
            z.theta1,
            z.theta2
        );
        System.out.printf(
            Locale.US,
            "  final parameter-error norm = %.6e%n",
            Math.hypot(d1, d2)
        );
        System.out.printf(
            Locale.US,
            "  minimum Gramian eigenvalue = %.6e%n",
            gramianMinimum
        );
    }

    public static void main(String[] args) {
        Parameters p = new Parameters();
        try {
            runCase("non_pe", false, p);
            runCase("pe", true, p);
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }
    }
}
      

The MATLAB script performs the RK4 comparison and includes an optional programmatic Simulink builder. Set buildSimulink = true to generate a model with state \(z=[e,\hat{\theta}_1,\hat{\theta}_2]^T\). The CaseFlag constant selects non-PE or PE excitation.

Chapter10_Lesson4.m

% Chapter10_Lesson4.m
% State convergence versus parameter convergence.
%
% The script simulates
%   e_dot         = -a*e + theta_tilde'*phi(t)
%   theta_hat_dot = -gamma*phi(t)*e
% for non-PE and PE regressors. It can also build an equivalent
% continuous-time Simulink model when Simulink is installed.

clear;
clc;
close all;

p.a = 2.0;
p.gamma = 5.0;
p.thetaStar = [2.0; -0.5];
p.tFinal = 40.0;
p.dt = 0.002;

nonPE = runCase(false, p);
PE = runCase(true, p);

reportCase("non_pe", nonPE, p);
reportCase("pe", PE, p);

writematrix( ...
    [nonPE.t, nonPE.z, nonPE.z(:, 2:3) - p.thetaStar.', ...
     vecnorm(nonPE.z(:, 2:3) - p.thetaStar.', 2, 2)], ...
    "Chapter10_Lesson4_non_pe.csv" ...
);
writematrix( ...
    [PE.t, PE.z, PE.z(:, 2:3) - p.thetaStar.', ...
     vecnorm(PE.z(:, 2:3) - p.thetaStar.', 2, 2)], ...
    "Chapter10_Lesson4_pe.csv" ...
);

figure;
plot(nonPE.t, nonPE.z(:, 1), "LineWidth", 1.2);
hold on;
plot(PE.t, PE.z(:, 1), "LineWidth", 1.2);
grid on;
xlabel("Time (s)");
ylabel("Tracking error e(t)");
title("State/tracking-error convergence");
legend("non-PE", "PE", "Location", "best");

figure;
semilogy( ...
    nonPE.t, ...
    max(vecnorm(nonPE.z(:, 2:3) - p.thetaStar.', 2, 2), 1e-14), ...
    "LineWidth", 1.2 ...
);
hold on;
semilogy( ...
    PE.t, ...
    max(vecnorm(PE.z(:, 2:3) - p.thetaStar.', 2, 2), 1e-14), ...
    "LineWidth", 1.2 ...
);
grid on;
xlabel("Time (s)");
ylabel("Parameter-error norm");
title("Parameter convergence requires sufficient excitation");
legend("non-PE", "PE", "Location", "best");

% Set this flag to true to generate Chapter10_Lesson4_Simulink.slx.
buildSimulink = false;
if buildSimulink
    buildSimulinkModel(p);
end


function result = runCase(persistentExcitation, p)
    numberOfSamples = round(p.tFinal / p.dt) + 1;
    t = linspace(0.0, p.tFinal, numberOfSamples).';
    z = zeros(numberOfSamples, 3);
    z(1, :) = [1.0, 0.0, 0.0];

    for k = 1:numberOfSamples - 1
        z(k + 1, :) = rk4Step( ...
            t(k), z(k, :).', p, persistentExcitation ...
        ).';
    end

    result.t = t;
    result.z = z;
    result.minGramianEigenvalue = minimumGramianEigenvalue( ...
        t, persistentExcitation, 8.0, 250 ...
    );
end


function dz = dynamics(t, z, p, persistentExcitation)
    phi = regressor(t, persistentExcitation);
    thetaTilde = z(2:3) - p.thetaStar;

    dz = zeros(3, 1);
    dz(1) = -p.a * z(1) + thetaTilde.' * phi;
    dz(2:3) = -p.gamma * phi * z(1);
end


function zNext = rk4Step(t, z, p, persistentExcitation)
    h = p.dt;
    k1 = dynamics(t, z, p, persistentExcitation);
    k2 = dynamics(t + 0.5 * h, z + 0.5 * h * k1, ...
                  p, persistentExcitation);
    k3 = dynamics(t + 0.5 * h, z + 0.5 * h * k2, ...
                  p, persistentExcitation);
    k4 = dynamics(t + h, z + h * k3, p, persistentExcitation);

    zNext = z + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
end


function phi = regressor(t, persistentExcitation)
    if persistentExcitation
        phi = [sin(t); cos(1.7 * t)];
    else
        phi = [sin(t); 2.0 * sin(t)];
    end
end


function minimum = minimumGramianEigenvalue( ...
    t, persistentExcitation, windowLength, stride ...
)
    dt = t(2) - t(1);
    windowSamples = round(windowLength / dt);
    minimum = inf;

    for startIndex = 1:stride:(numel(t) - windowSamples)
        gramian = zeros(2, 2);
        for k = startIndex:(startIndex + windowSamples - 1)
            phi = regressor(t(k), persistentExcitation);
            gramian = gramian + (phi * phi.') * dt;
        end
        eigenvalues = eig(gramian);
        minimum = min(minimum, min(eigenvalues));
    end
end


function reportCase(name, result, p)
    finalTheta = result.z(end, 2:3).';
    finalParameterError = finalTheta - p.thetaStar;

    fprintf("\nCase: %s\n", name);
    fprintf("  final tracking error       = %.6e\n", result.z(end, 1));
    fprintf( ...
        "  final parameter estimate  = [%.6f, %.6f]\n", ...
        finalTheta(1), finalTheta(2) ...
    );
    fprintf( ...
        "  final parameter-error norm = %.6e\n", ...
        norm(finalParameterError) ...
    );
    fprintf( ...
        "  minimum Gramian eigenvalue = %.6e\n", ...
        result.minGramianEigenvalue ...
    );
end


function buildSimulinkModel(p)
    % Build a three-state model z = [e; thetaHat1; thetaHat2].
    % CaseFlag = 0 selects non-PE; CaseFlag = 1 selects PE.
    modelName = "Chapter10_Lesson4_Simulink";

    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if isfile(modelName + ".slx")
        delete(modelName + ".slx");
    end

    new_system(modelName);
    open_system(modelName);

    add_block( ...
        "simulink/Sources/Clock", ...
        modelName + "/Clock", ...
        "Position", [40, 55, 70, 75] ...
    );
    add_block( ...
        "simulink/Sources/Constant", ...
        modelName + "/CaseFlag", ...
        "Value", "1", ...
        "Position", [35, 125, 70, 145] ...
    );
    add_block( ...
        "simulink/User-Defined Functions/MATLAB Function", ...
        modelName + "/AdaptiveDynamics", ...
        "Position", [145, 50, 315, 155] ...
    );
    add_block( ...
        "simulink/Continuous/Integrator", ...
        modelName + "/StateIntegrator", ...
        "InitialCondition", "[1; 0; 0]", ...
        "Position", [375, 75, 405, 125] ...
    );
    add_block( ...
        "simulink/Sinks/To Workspace", ...
        modelName + "/StateLog", ...
        "VariableName", "zSim", ...
        "SaveFormat", "Structure With Time", ...
        "Position", [475, 75, 565, 105] ...
    );

    root = sfroot;
    chart = root.find( ...
        "-isa", "Stateflow.EMChart", ...
        "Path", modelName + "/AdaptiveDynamics" ...
    );

    chart.Script = sprintf([ ...
        "function dz = fcn(t,z,caseFlag)\n" ...
        "a = %.17g;\n" ...
        "gamma = %.17g;\n" ...
        "thetaStar = [%.17g; %.17g];\n" ...
        "if caseFlag >= 0.5\n" ...
        "    phi = [sin(t); cos(1.7*t)];\n" ...
        "else\n" ...
        "    phi = [sin(t); 2*sin(t)];\n" ...
        "end\n" ...
        "thetaTilde = z(2:3) - thetaStar;\n" ...
        "dz = zeros(3,1);\n" ...
        "dz(1) = -a*z(1) + thetaTilde'*phi;\n" ...
        "dz(2:3) = -gamma*phi*z(1);\n" ...
        "end\n" ...
    ], p.a, p.gamma, p.thetaStar(1), p.thetaStar(2));

    add_line(modelName, "Clock/1", "AdaptiveDynamics/1");
    add_line(modelName, "StateIntegrator/1", "AdaptiveDynamics/2");
    add_line(modelName, "CaseFlag/1", "AdaptiveDynamics/3");
    add_line(modelName, "AdaptiveDynamics/1", "StateIntegrator/1");
    add_line(modelName, "StateIntegrator/1", "StateLog/1");

    set_param( ...
        modelName, ...
        "StopTime", num2str(p.tFinal), ...
        "Solver", "ode45" ...
    );
    save_system(modelName);
    fprintf("Created %s.slx\n", modelName);
end
      

The notebook uses NDSolveValue for the differential equations, NIntegrate for the excitation Gramian, and Eigenvalues for the PE diagnostic.

Chapter10_Lesson4.nb

Notebook[{
  Cell["Chapter 10, Lesson 4: State Convergence vs Parameter Convergence", "Title"],
  Cell["Adaptive error-model simulation with non-PE and PE regressors.", "Text"], Cell[BoxData[
"ClearAll[\"Global`*\"]; a = 2.0; gamma = 5.0;
thetaStar = {2.0, -0.5}; tFinal = 40.0;

regressor[t_, persistentExcitation_] := If[persistentExcitation,
  {Sin[t], Cos[1.7 t]}, {Sin[t], 2 Sin[t]}];
simulate[persistentExcitation_] := Module[{solution, e, theta1, theta2},
  solution = NDSolveValue[{
    e'[t] == -a e[t] + ({theta1[t], theta2[t]} - thetaStar).
      regressor[t, persistentExcitation],
    theta1'[t] == -gamma regressor[t, persistentExcitation][[1]] e[t],
    theta2'[t] == -gamma regressor[t, persistentExcitation][[2]] e[t],
    e[0] == 1.0, theta1[0] == 0.0, theta2[0] == 0.0},
    {e, theta1, theta2}, {t, 0, tFinal},
    Method -> {\"TimeIntegration\" -> {\"ExplicitRungeKutta\"} }];
  solution
];
minimumGramianEigenvalue[persistentExcitation_, windowLength_: 8.0] :=
 Module[{starts, gramian, values},
  starts = Range[0, tFinal - windowLength, 0.5];
  values = Table[
    gramian = NIntegrate[
      With[{phi = regressor[tau, persistentExcitation]},
        Outer[Times, phi, phi]], {tau, start, start + windowLength}];
    Min[Eigenvalues[gramian]], {start, starts}];
  Min[values]
];

nonPE = simulate[False]; pe = simulate[True];
finalReport[name_, solution_, persistentExcitation_] := Module[
  {finalError, finalTheta, parameterError},
  finalError = solution[[1]][tFinal];
  finalTheta = {solution[[2]][tFinal], solution[[3]][tFinal]};
  parameterError = finalTheta - thetaStar;
  Print[\"Case: \", name];
  Print[\"  final tracking error       = \", N[finalError, 8]];
  Print[\"  final parameter estimate  = \", N[finalTheta, 8]];
  Print[\"  final parameter-error norm = \", N[Norm[parameterError], 8]];
  Print[\"  minimum Gramian eigenvalue = \",
    N[minimumGramianEigenvalue[persistentExcitation], 8]];
];

finalReport[\"non_pe\", nonPE, False]; finalReport[\"pe\", pe, True];
trackingPlot = Plot[Evaluate[{nonPE[[1]][t], pe[[1]][t]}], {t, 0, tFinal},
  PlotLegends -> {\"non-PE\", \"PE\"}, AxesLabel -> {\"Time (s)\", \"e(t)\"},
  PlotLabel -> \"State/tracking-error convergence\",
  GridLines -> Automatic, PlotRange -> All];
parameterPlot = LogPlot[Evaluate[{
    Norm[{nonPE[[2]][t], nonPE[[3]][t]} - thetaStar],
    Norm[{pe[[2]][t], pe[[3]][t]} - thetaStar]}], {t, 0, tFinal},
  PlotLegends -> {\"non-PE\", \"PE\"},
  AxesLabel -> {\"Time (s)\", \"Parameter-error norm\"},
  PlotLabel -> \"Parameter convergence requires excitation\",
  GridLines -> Automatic, PlotRange -> All];
Print[trackingPlot]; Print[parameterPlot];
Export[\"Chapter10_Lesson4_non_pe.csv\", Table[{
  time, nonPE[[1]][time], nonPE[[2]][time], nonPE[[3]][time],
  nonPE[[2]][time] - thetaStar[[1]], nonPE[[3]][time] - thetaStar[[2]]},
  {time, 0, tFinal, 0.01}]];
Export[\"Chapter10_Lesson4_pe.csv\", Table[{
  time, pe[[1]][time], pe[[2]][time], pe[[3]][time],
  pe[[2]][time] - thetaStar[[1]], pe[[3]][time] - thetaStar[[2]]},
  {time, 0, tFinal, 0.01}]];"], "Input"]
}, WindowSize -> {1200, 800},
WindowMargins -> { {Automatic, 40}, {Automatic, 40} }, StyleDefinitions -> "Default.nb"]
      

For \(a=2\), \(\Gamma=5I\), \(e(0)=1\), and \(\hat{\theta}(0)=0\), the implementations give approximately:

The non-PE case is the decisive counterexample: perfect asymptotic tracking does not certify correct parameter identification.

Misinterpretation 1: Since \(V\) converges, every error converges to zero.

Correction: A negative-semidefinite derivative may leave nontrivial invariant directions.

Misinterpretation 2: A settled estimate must be the true parameter.

Correction: It may be one member of an input-output-equivalent parameter set.

Misinterpretation 3: A sinusoid is always PE.

Correction: PE depends on regressor dimension and rank.

Misinterpretation 4: Increasing the adaptation gain creates missing information.

Correction: Gain changes transient speed and noise sensitivity, not Gramian rank.

Misinterpretation 5: Zero tracking error proves the model is physically correct.

Correction: It proves compatibility with the observed trajectory.

Problem 1 (Lyapunov conclusion): Prove \(e(t)\to0\) for the model in Section 2 under bounded \(\phi(t)\). What can be concluded about \(\tilde{\theta}(t)\) without PE?

Solution:

The Lyapunov function gives \(\dot V=-ae^2\), so \(e,\tilde{\theta}\in L_\infty\) and \(e\in L_2\). Bounded \(\phi\) makes \(\dot e\) bounded. Barbalat's lemma yields \(e\to0\). Without PE, only boundedness of \(\tilde{\theta}\) follows generally; its limit may be nonzero.

Problem 2 (Unidentifiable subspace): For \(\phi(t)=[\sin t,\;2\sin t]^T\), determine the Gramian null direction and all estimates indistinguishable from \(\theta^*=[2,-0.5]^T\).

Solution:

Write \(\phi=q\sin t\) with \(q=[1,2]^T\). Then

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

so the rank is one and a null vector is \([-2,1]^T\). The observable combination is \(\theta_1+2\theta_2\). Therefore,

\[ \left\{ \hat{\theta}: \hat{\theta}_1+2\hat{\theta}_2=1 \right\} = \theta^*+\operatorname{span}\left\{[-2,1]^T\right\}. \]

Problem 3 (Exact PE Gramian): For \(\phi(t)=[\sin(\omega t),\cos(\omega t)]^T\), evaluate the Gramian over \(T=2\pi/\omega\).

Solution:

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

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

Hence

\[ G\left(t,\frac{2\pi}{\omega}\right) = \frac{\pi}{\omega}I, \]

uniformly in \(t\).

Problem 4 (Tracking versus state convergence): If \(e=x-x_m\to0\) while a stable reference model approaches a nonzero periodic response, does \(x\to0\)?

Solution:

No. Since \(x=x_m+e\), the plant approaches the same periodic reference-model trajectory; only their difference approaches zero.

Problem 5 (Adaptation gain): Can replacing \(\Gamma=5I\) by \(\Gamma=50I\) make the non-PE regressor identify the true two-parameter vector?

Solution:

No. The Gramian depends on \(\phi\), not on \(\Gamma\). A larger gain changes transients but does not remove the Gramian nullspace.

The standard adaptive Lyapunov function proves boundedness of state and parameter errors and, with Barbalat's lemma, convergence of the tracking error. Because its derivative is only negative semidefinite in the combined state-parameter vector, exact parameter convergence does not follow. Non-PE regressors identify only certain parameter combinations and leave an invariant unidentifiable subspace. PE makes the moving-window excitation Gramian uniformly positive definite and permits simultaneous tracking and parameter convergence.