A conventional observer assumes that the plant matrices and nonlinear terms used by the observer are known. For a linear time-invariant plant,

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}, \qquad \mathbf{y}=\mathbf{C}\mathbf{x}, \]

a full-order Luenberger observer has the familiar form

\[ \dot{\hat{\mathbf{x} } } =\mathbf{A}\hat{\mathbf{x} }+\mathbf{B}\mathbf{u} +\mathbf{L}\left(\mathbf{y}-\hat{\mathbf{y} }\right), \qquad \hat{\mathbf{y} }=\mathbf{C}\hat{\mathbf{x} }. \]

With the state-estimation error \( \mathbf{e}=\mathbf{x}-\hat{\mathbf{x} } \), the error dynamics are autonomous:

\[ \dot{\mathbf{e} }=(\mathbf{A}-\mathbf{L}\mathbf{C})\mathbf{e}. \]

The standard design problem is therefore to choose \( \mathbf{L} \) so that \( \mathbf{A}-\mathbf{L}\mathbf{C} \) is Hurwitz. When the dynamics contain an unknown parameter vector, however, copying the plant model into the observer introduces a forcing term proportional to the parameter error. The observer and identifier become coupled, and pole placement alone no longer proves convergence.

An adaptive observer combines three mechanisms:

1. a dynamic state estimator,
2. an output innovation or prediction error, and
3. an online parameter-update law driven by measurable signals.

Two goals must be kept logically separate:

• State observation: \( \mathbf{e}(t)\longrightarrow\mathbf{0} \).
• Parameter identification: \( \tilde{\boldsymbol{\theta} }(t) =\boldsymbol{\theta}-\hat{\boldsymbol{\theta} }(t) \longrightarrow\mathbf{0} \).

The first goal may be achieved without the second. A plant can be observed accurately even when the data do not uniquely reveal every unknown parameter.

Consider the single-output system

\[ \dot{\mathbf{x} } =\mathbf{A}\mathbf{x} +\mathbf{b}\left[u+\boldsymbol{\theta}^{T}\boldsymbol{\phi}(t)\right], \qquad y=\mathbf{c}^{T}\mathbf{x}, \]

where \( \mathbf{x}\in\mathbb{R}^{n} \) is not fully measured, \( y\in\mathbb{R} \) is measured, \( \boldsymbol{\theta}\in\mathbb{R}^{p} \) is an unknown constant parameter vector, and \( \boldsymbol{\phi}(t)\in\mathbb{R}^{p} \) is a known bounded regressor constructed from measured or filtered signals.

A natural adaptive observer is

\[ \dot{\hat{\mathbf{x} } } =\mathbf{A}\hat{\mathbf{x} } +\mathbf{b}\left[u+\hat{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}(t)\right] +\mathbf{L}(y-\hat{y}), \qquad \hat{y}=\mathbf{c}^{T}\hat{\mathbf{x} }. \]

Define \( \mathbf{A}_{o}=\mathbf{A}-\mathbf{L}\mathbf{c}^{T} \), \( e_y=y-\hat{y}=\mathbf{c}^{T}\mathbf{e} \), and \( \tilde{\boldsymbol{\theta} } =\boldsymbol{\theta}-\hat{\boldsymbol{\theta} } \). Since the true parameter is constant,

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

flowchart TD
  U["Known input u"] --> P["Plant with unknown parameter theta"]
  P --> Y["Measured output y"]
  U --> O["Observer model"]
  Y --> E["Output error: ey = y - yhat"]
  O --> E
  E --> I["Innovation term L times ey"]
  I --> O
  E --> A["Parameter update law"]
  R["Known regressor phi"] --> A
  A --> T["Parameter estimate thetaHat"]
  T --> O
  O --> XH["Estimated state xHat"]
        

The difficulty is visible in the error equation: even if \( \mathbf{A}_{o} \) is Hurwitz, the unknown forcing \( \mathbf{b}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi} \) prevents a conventional observer proof. The adaptation law must cancel this coupling in the derivative of a composite Lyapunov function.

Let the transfer function from the scalar mismatch signal \( v=\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \) to the measurable output error \( e_y \) be

\[ H(s)=\mathbf{c}^{T}(s\mathbf{I}-\mathbf{A}_{o})^{-1}\mathbf{b}. \]

A central classical construction requires \( H(s) \) to be strictly positive real (SPR). Informally, SPR means that the stable linear error channel has a strictly passive input-output map. For a minimal realization, the Kalman-Yakubovich-Popov lemma gives symmetric matrices \( \mathbf{P}\succ 0 \) and \( \mathbf{Q}\succ 0 \) satisfying

\[ \mathbf{A}_{o}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{o} =-\mathbf{Q}, \qquad \mathbf{P}\mathbf{b}=\mathbf{c}. \]

The second identity is the essential matching relation. It converts the unmeasured state cross term \( \mathbf{e}^{T}\mathbf{P}\mathbf{b} \) into the measurable output error:

\[ \mathbf{e}^{T}\mathbf{P}\mathbf{b} =\mathbf{e}^{T}\mathbf{c}=e_y. \]

Consequently, the adaptation law can use only \( e_y \) and the known regressor; it does not require the unavailable full state-estimation error.

The SPR requirement is structural, not merely a gain-tuning detail. A stable transfer function of relative degree two or more cannot be SPR in the standard scalar sense. Classical adaptive-observer designs often introduce stable filters, nonminimal realizations, output injections, or canonical transformations to produce an SPR error model.

Choose a symmetric adaptation-gain matrix \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ 0 \) and the gradient update law

\[ \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}(t)e_y. \]

Consider the composite Lyapunov function

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

Differentiating along the estimation-error dynamics gives

\[ \begin{aligned} \dot{V} &=\mathbf{e}^{T} (\mathbf{A}_{o}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{o})\mathbf{e} +2\mathbf{e}^{T}\mathbf{P}\mathbf{b} \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \\ &\quad +2\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }. \end{aligned} \]

Use the KYP identities and \( \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_y \):

\[ \begin{aligned} \dot{V} &=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} +2e_y\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} -2\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}e_y \\ &=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}. \end{aligned} \]

Therefore \( V(t) \) is nonincreasing, and both \( \mathbf{e}(t) \) and \( \tilde{\boldsymbol{\theta} }(t) \) are bounded. Also,

\[ \int_{0}^{\infty} \mathbf{e}^{T}(t)\mathbf{Q}\mathbf{e}(t)\,dt \le V(0), \]

so \( \mathbf{e}\in\mathcal{L}_{2} \). If \( u \) and \( \boldsymbol{\phi} \) are bounded, the error dynamics imply bounded \( \dot{\mathbf{e} } \). Hence \( \mathbf{e}(t) \) is uniformly continuous, and Barbalat's lemma yields

\[ \lim_{t\to\infty}\mathbf{e}(t)=\mathbf{0}, \qquad \lim_{t\to\infty}e_y(t)=0. \]

This proves asymptotic state observation. The proof does not yet prove \( \tilde{\boldsymbol{\theta} }(t)\longrightarrow\mathbf{0} \).

Parameter convergence requires that the measured experiment contain enough independent information. A bounded regressor \( \boldsymbol{\phi}(t) \) is persistently exciting if there exist constants \( T>0 \) and \( \alpha>0 \) such that, for every \( t\ge 0 \),

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

In a filtered adaptive-observer derivation, the relevant excitation may be the original regressor, a filtered regressor, or an observability sensitivity signal. Under the required PE condition and regularity assumptions, the zero output-error limit forces \( \tilde{\boldsymbol{\theta} } \) to zero.

To see why PE matters, suppose \( \boldsymbol{\phi}(t)=\mathbf{0} \) after some time. Then the update law stops. The stable observer error can still decay to zero, but every parameter estimate compatible with the final output behavior remains possible. More generally, if the regressor explores only a lower-dimensional subspace, only the projection of the parameter vector onto that subspace can be identified.

Three notions should not be conflated:

• Observability: whether the state is reconstructible when the model is known.
• Structural identifiability: whether distinct parameters can generate identical ideal input-output behavior.
• Excitation/informativity: whether the particular experiment activates the directions needed to distinguish parameters.

No adaptive law can recover a parameter that is structurally absent from the output map or never excited by the data. Increasing \( \boldsymbol{\Gamma} \) changes adaptation speed; it does not manufacture information.

The ideal proof assumes exact measurements, an exact linearly parameterized model, and constant parameters. Let the measured output be

\[ y_m=\mathbf{c}^{T}\mathbf{x}+n(t), \]

and let the plant contain an unmatched residual \( \mathbf{d}(t) \). The adaptive law then reacts to noise and approximation error. A large adaptation gain may amplify high-frequency measurement noise and cause parameter drift.

Standard robustifications from earlier chapters can be inserted:

\[ \dot{\hat{\boldsymbol{\theta} } } =\operatorname{Proj}_{\Omega} \left( \hat{\boldsymbol{\theta} }, \boldsymbol{\Gamma}\boldsymbol{\phi}e_y -\sigma\boldsymbol{\Gamma}\hat{\boldsymbol{\theta} } \right), \qquad \sigma>0. \]

Here projection confines the estimate to a known admissible set \( \Omega \), while the leakage term produces an ultimate bound instead of exact convergence. A dead zone is another option:

\[ \dot{\hat{\boldsymbol{\theta} } } =\begin{cases} \boldsymbol{\Gamma}\boldsymbol{\phi}e_y, & |e_y|>\delta,\\ \mathbf{0}, & |e_y|\le\delta. \end{cases} \]

Normalization limits updates when the regressor becomes large:

\[ \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma} \frac{\boldsymbol{\phi}e_y} {1+\boldsymbol{\phi}^{T}\boldsymbol{\phi} }. \]

In practice, the observer bandwidth, adaptation bandwidth, sensor bandwidth, sampling rate, and unmodeled dynamics must be designed as one coupled system. A fast state observer does not justify an arbitrarily fast parameter update.

Deterministic adaptive observers usually seek stability and asymptotic convergence for unknown but fixed quantities. Stochastic filters model uncertainty probabilistically. Consider

\[ \mathbf{x}_{k+1} =\mathbf{f}(\mathbf{x}_{k},\mathbf{u}_{k},\boldsymbol{\theta}_{k}) +\mathbf{w}_{k}, \qquad \mathbf{y}_{k} =\mathbf{h}(\mathbf{x}_{k},\boldsymbol{\theta}_{k})+\mathbf{v}_{k}, \]

with process noise \( \mathbf{w}_{k}\sim\mathcal{N}(\mathbf{0},\mathbf{Q}_{k}) \) and measurement noise \( \mathbf{v}_{k}\sim\mathcal{N}(\mathbf{0},\mathbf{R}_{k}) \). A common joint-estimation device augments the state with the unknown parameter:

\[ \mathbf{z}_{k} =\begin{bmatrix}\mathbf{x}_{k}\\ \boldsymbol{\theta}_{k}\end{bmatrix}, \qquad \boldsymbol{\theta}_{k+1} =\boldsymbol{\theta}_{k}+\boldsymbol{\eta}_{k}. \]

The artificial random walk \( \boldsymbol{\eta}_{k} \) gives the filter a tunable parameter process covariance. For an augmented extended Kalman filter, define the nonlinear augmented transition \( \mathbf{z}_{k+1}=\mathbf{F}_{a}(\mathbf{z}_{k},\mathbf{u}_{k}) \) and measurement \( \mathbf{y}_{k}=\mathbf{H}_{a}(\mathbf{z}_{k})+\mathbf{v}_{k} \). Linearizing about the current estimate yields Jacobians

\[ \mathbf{F}_{k} =\left.\frac{\partial\mathbf{F}_{a} }{\partial\mathbf{z} } \right|_{\hat{\mathbf{z} }_{k|k} }, \qquad \mathbf{H}_{k} =\left.\frac{\partial\mathbf{H}_{a} }{\partial\mathbf{z} } \right|_{\hat{\mathbf{z} }_{k|k-1} }. \]

The prediction equations are

\[ \hat{\mathbf{z} }_{k|k-1} =\mathbf{F}_{a}(\hat{\mathbf{z} }_{k-1|k-1},\mathbf{u}_{k-1}), \]

\[ \mathbf{P}_{k|k-1} =\mathbf{F}_{k-1}\mathbf{P}_{k-1|k-1}\mathbf{F}_{k-1}^{T} +\mathbf{Q}_{a,k-1}. \]

The correction equations are

\[ \mathbf{S}_{k} =\mathbf{H}_{k}\mathbf{P}_{k|k-1}\mathbf{H}_{k}^{T} +\mathbf{R}_{k}, \qquad \mathbf{K}_{k} =\mathbf{P}_{k|k-1}\mathbf{H}_{k}^{T}\mathbf{S}_{k}^{-1}, \]

\[ \hat{\mathbf{z} }_{k|k} =\hat{\mathbf{z} }_{k|k-1} +\mathbf{K}_{k} \left[\mathbf{y}_{k}-\mathbf{H}_{a}(\hat{\mathbf{z} }_{k|k-1})\right], \]

\[ \mathbf{P}_{k|k} =(\mathbf{I}-\mathbf{K}_{k}\mathbf{H}_{k}) \mathbf{P}_{k|k-1} (\mathbf{I}-\mathbf{K}_{k}\mathbf{H}_{k})^{T} +\mathbf{K}_{k}\mathbf{R}_{k}\mathbf{K}_{k}^{T}. \]

The last expression is the Joseph covariance update, preferred in finite precision because it better preserves symmetry and positive semidefiniteness.

State augmentation does not eliminate identifiability requirements. If the output is insensitive to a parameter, the state-parameter cross covariance cannot reliably update it. Poor scaling, weak excitation, incorrect noise covariances, and linearization error can make the joint EKF appear confident while being wrong.

This lesson does not replace a Stochastic Control course. A complete treatment there should derive conditional expectations, innovations, Riccati equations, Bayesian filtering, smoothing, covariance consistency, and nonlinear filters. Similarly, a complete Observer Design course should cover reduced-order observers, unknown-input observers, high-gain observers, sliding-mode observers, nonlinear observability, and geometric observer forms.

flowchart TD
  S["Specify plant, measured outputs, and unknown parameters"] --> O["Check observability for fixed parameters"]
  O --> I["Check structural identifiability and signal availability"]
  I --> R["Construct measurable or filtered regressor"]
  R --> C["Choose observer injection and stable error matrix"]
  C --> P["Verify SPR or solve an equivalent Lyapunov/LMI condition"]
  P --> A["Derive adaptation law from a composite Lyapunov function"]
  A --> B["Add projection, leakage, normalization, or dead zone"]
  B --> E["Design excitation and sampling experiment"]
  E --> V["Validate state error, parameter error, residuals, and numerical conditioning"]
        

A credible adaptive-observer design should explicitly answer the following questions:

1. Which states and parameters are being estimated?
2. Which signals in the observer and update law are actually measured?
3. What is the exact estimation-error system?
4. Which matrix or passivity condition makes the Lyapunov proof work?
5. What assumptions guarantee state convergence?
6. What additional excitation assumptions guarantee parameter convergence?
7. How are noise, saturation, sampling, and parameter bounds handled?
8. Which residuals reveal model mismatch or estimator inconsistency?

Consider the plant

\[ \dot{\mathbf{x} } =\underbrace{\begin{bmatrix}-1&1\\-2&-3\end{bmatrix} }_{\mathbf{A} } \mathbf{x} +\underbrace{\begin{bmatrix}1\\0\end{bmatrix} }_{\mathbf{b} } \left[u+\theta\phi(t)\right], \qquad y=\underbrace{\begin{bmatrix}1&0\end{bmatrix} }_{\mathbf{c}^{T} } \mathbf{x}. \]

The scalar parameter \( \theta \) is unknown. Choose

\[ \mathbf{L}=\begin{bmatrix}1\\0\end{bmatrix}, \qquad \mathbf{A}_{o}=\mathbf{A}-\mathbf{L}\mathbf{c}^{T} =\begin{bmatrix}-2&1\\-2&-3\end{bmatrix}. \]

The mismatch-to-output-error transfer function is

\[ H(s)=\mathbf{c}^{T}(s\mathbf{I}-\mathbf{A}_{o})^{-1}\mathbf{b} =\frac{s+3}{s^{2}+5s+8}. \]

For \( s=j\omega \), direct calculation gives

\[ \operatorname{Re}\{H(j\omega)\} =\frac{24+2\omega^{2} } {(8-\omega^{2})^{2}+25\omega^{2} }, \]

which is strictly positive for every real \( \omega \). Thus \( H(s) \) is SPR. The matrices

\[ \mathbf{P}=\begin{bmatrix}1&0\\0&\tfrac{1}{2}\end{bmatrix}, \qquad \mathbf{Q}=\begin{bmatrix}4&0\\0&3\end{bmatrix} \]

satisfy exactly

\[ \mathbf{A}_{o}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{o}=-\mathbf{Q}, \qquad \mathbf{P}\mathbf{b}=\mathbf{c}. \]

The adaptive observer is therefore

\[ \dot{\hat{\mathbf{x} } } =\mathbf{A}\hat{\mathbf{x} } +\mathbf{b}\left[u+\hat{\theta}\phi(t)\right] +\mathbf{L}(y-\hat{y}), \qquad \dot{\hat{\theta} }=\gamma\phi(t)(y-\hat{y}). \]

With \( V=\mathbf{e}^{T}\mathbf{P}\mathbf{e} +\tilde{\theta}^{2}/\gamma \), the derivative is the exact quadratic expression

\[ \dot{V}=-4e_{1}^{2}-3e_{2}^{2}. \]

The implementations below use \( \theta=1.4 \), \( \gamma=8 \),

\[ u(t)=\sin(0.55t)+0.35\sin(1.70t), \qquad \phi(t)=0.70\sin(0.90t)+0.45u(t), \]

and projection to the interval \( [-3,3] \). The multisine regressor provides rich excitation for this scalar parameter.

The Python implementation uses NumPy for matrix operations, SciPy for adaptive-step integration, and optional Matplotlib visualization. For larger filtering projects, relevant libraries include filterpy, python-control, CasADi, and probabilistic programming tools; the present Lyapunov observer is implemented directly so that every equation remains visible.

Chapter28_Lesson3.py

"""Chapter28_Lesson3.py

Deterministic adaptive observer for a stable second-order SISO system.
The example implements the SPR/Lyapunov construction developed in the lesson.

Dependencies:
    pip install numpy scipy matplotlib
"""

from __future__ import annotations

import numpy as np
from scipy.integrate import solve_ivp


A = np.array([[-1.0, 1.0], [-2.0, -3.0]])
b = np.array([1.0, 0.0])
C = np.array([1.0, 0.0])
L = np.array([1.0, 0.0])
P = np.diag([1.0, 0.5])
Q = np.diag([4.0, 3.0])

THETA_TRUE = 1.40
GAMMA = 8.0
THETA_MIN = -3.0
THETA_MAX = 3.0


def command(t: float) -> float:
    """Known plant command."""
    return np.sin(0.55 * t) + 0.35 * np.sin(1.70 * t)


def regressor(t: float, u: float) -> float:
    """Known persistently rich scalar regressor."""
    return 0.70 * np.sin(0.90 * t) + 0.45 * u


def project_scalar(theta_hat: float, raw_update: float) -> float:
    """Continuous-time projection that prevents outward motion at the bounds."""
    if theta_hat <= THETA_MIN and raw_update < 0.0:
        return 0.0
    if theta_hat >= THETA_MAX and raw_update > 0.0:
        return 0.0
    return raw_update


def dynamics(t: float, z: np.ndarray) -> np.ndarray:
    """Combined plant, observer, and adaptation dynamics."""
    x = z[0:2]
    x_hat = z[2:4]
    theta_hat = float(z[4])

    u = command(t)
    phi = regressor(t, u)
    y = float(C @ x)
    y_hat = float(C @ x_hat)
    output_error = y - y_hat

    x_dot = A @ x + b * (u + THETA_TRUE * phi)
    x_hat_dot = A @ x_hat + b * (u + theta_hat * phi) + L * output_error

    raw_theta_dot = GAMMA * phi * output_error
    theta_hat_dot = project_scalar(theta_hat, raw_theta_dot)

    return np.concatenate((x_dot, x_hat_dot, [theta_hat_dot]))


def main() -> None:
    t_span = (0.0, 30.0)
    t_eval = np.linspace(t_span[0], t_span[1], 6001)
    z0 = np.array([1.2, -0.8, 0.0, 0.0, -0.5])

    solution = solve_ivp(
        dynamics,
        t_span,
        z0,
        t_eval=t_eval,
        rtol=1e-9,
        atol=1e-11,
        max_step=0.02,
    )
    if not solution.success:
        raise RuntimeError(solution.message)

    x = solution.y[0:2, :]
    x_hat = solution.y[2:4, :]
    theta_hat = solution.y[4, :]
    error = x - x_hat

    rmse = np.sqrt(np.mean(error**2, axis=1))
    final_error = error[:, -1]
    final_theta_error = THETA_TRUE - theta_hat[-1]

    # Verify the exact matrices used in the proof.
    A_o = A - np.outer(L, C)
    lyapunov_residual = A_o.T @ P + P @ A_o + Q
    matching_residual = P @ b - C

    print("Adaptive observer simulation")
    print(f"RMSE state 1: {rmse[0]:.6e}")
    print(f"RMSE state 2: {rmse[1]:.6e}")
    print(f"Final state error: {final_error}")
    print(f"True theta: {THETA_TRUE:.6f}")
    print(f"Final theta_hat: {theta_hat[-1]:.6f}")
    print(f"Final parameter error: {final_theta_error:.6e}")
    print(f"Lyapunov identity residual norm: {np.linalg.norm(lyapunov_residual):.3e}")
    print(f"SPR matching residual norm: {np.linalg.norm(matching_residual):.3e}")

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(solution.t, x[0], label="x1")
        plt.plot(solution.t, x_hat[0], "--", label="x1_hat")
        plt.plot(solution.t, x[1], label="x2")
        plt.plot(solution.t, x_hat[1], "--", label="x2_hat")
        plt.xlabel("Time [s]")
        plt.ylabel("State")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        plt.figure()
        plt.plot(solution.t, theta_hat, label="theta_hat")
        plt.axhline(THETA_TRUE, linestyle="--", label="theta_true")
        plt.xlabel("Time [s]")
        plt.ylabel("Parameter")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()
    except ImportError:
        print("matplotlib is not installed; numerical simulation completed without plots.")


if __name__ == "__main__":
    main()

The C++ version implements fourth-order Runge-Kutta integration using the standard library. In production systems, Eigen is the usual choice for fixed- and dynamic-size matrix algebra; Boost.Odeint provides reusable numerical integrators, and automatic-differentiation libraries can support EKF Jacobians. The from-scratch form below avoids hiding the observer equations.

Chapter28_Lesson3.cpp

// Chapter28_Lesson3.cpp
// From-scratch RK4 simulation of the deterministic adaptive observer.
// Build: g++ -std=c++17 -O2 Chapter28_Lesson3.cpp -o Chapter28_Lesson3

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

using State = std::array<double, 5>; // x1, x2, xhat1, xhat2, theta_hat

constexpr double theta_true = 1.40;
constexpr double gamma_gain = 8.0;
constexpr double theta_min = -3.0;
constexpr double theta_max = 3.0;

static double command(double t) {
    return std::sin(0.55 * t) + 0.35 * std::sin(1.70 * t);
}

static double regressor(double t, double u) {
    return 0.70 * std::sin(0.90 * t) + 0.45 * u;
}

static double projected_update(double theta_hat, double raw_update) {
    if (theta_hat <= theta_min && raw_update < 0.0) return 0.0;
    if (theta_hat >= theta_max && raw_update > 0.0) return 0.0;
    return raw_update;
}

static State derivative(double t, const State& z) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xh1 = z[2];
    const double xh2 = z[3];
    const double theta_hat = z[4];

    const double u = command(t);
    const double phi = regressor(t, u);
    const double output_error = x1 - xh1;

    State dz{};
    dz[0] = -x1 + x2 + u + theta_true * phi;
    dz[1] = -2.0 * x1 - 3.0 * x2;
    dz[2] = -xh1 + xh2 + u + theta_hat * phi + output_error;
    dz[3] = -2.0 * xh1 - 3.0 * xh2;
    dz[4] = projected_update(theta_hat, gamma_gain * phi * output_error);
    return dz;
}

static State add_scaled(const State& a, const State& b, double scale) {
    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = a[i] + scale * b[i];
    }
    return result;
}

static State rk4_step(double t, double h, const State& z) {
    const State k1 = derivative(t, z);
    const State k2 = derivative(t + 0.5 * h, add_scaled(z, k1, 0.5 * h));
    const State k3 = derivative(t + 0.5 * h, add_scaled(z, k2, 0.5 * h));
    const State k4 = derivative(t + h, add_scaled(z, k3, h));

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + h * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
    }
    next[4] = std::clamp(next[4], theta_min, theta_max);
    return next;
}

int main() {
    constexpr double h = 0.002;
    constexpr double final_time = 30.0;
    const int steps = static_cast<int>(final_time / h);

    State z{1.2, -0.8, 0.0, 0.0, -0.5};
    double sum_sq_e1 = 0.0;
    double sum_sq_e2 = 0.0;

    for (int k = 0; k < steps; ++k) {
        const double t = k * h;
        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        sum_sq_e1 += e1 * e1;
        sum_sq_e2 += e2 * e2;
        z = rk4_step(t, h, z);
    }

    const double e1 = z[0] - z[2];
    const double e2 = z[1] - z[3];
    const double rmse1 = std::sqrt(sum_sq_e1 / steps);
    const double rmse2 = std::sqrt(sum_sq_e2 / steps);

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "Adaptive observer simulation\n";
    std::cout << "RMSE state 1: " << rmse1 << '\n';
    std::cout << "RMSE state 2: " << rmse2 << '\n';
    std::cout << "Final state error: [" << e1 << ", " << e2 << "]\n";
    std::cout << "True theta: " << theta_true << '\n';
    std::cout << "Final theta_hat: " << z[4] << '\n';
    std::cout << "Final parameter error: " << theta_true - z[4] << '\n';

    if (!std::isfinite(z[4])) {
        throw std::runtime_error("Non-finite parameter estimate detected.");
    }
    return 0;
}

The Java implementation also uses explicit RK4 and has no external dependency. For larger observer/filter projects, EJML or Apache Commons Math can provide matrices and decompositions. Real-time Java deployments should preallocate arrays, avoid per-step object creation, and separate estimator state from visualization and logging threads.

Chapter28_Lesson3.java

// Chapter28_Lesson3.java
// From-scratch RK4 simulation of the deterministic adaptive observer.
// Build and run: javac Chapter28_Lesson3.java && java Chapter28_Lesson3

import java.util.Locale;

public final class Chapter28_Lesson3 {
    private static final double THETA_TRUE = 1.40;
    private static final double GAMMA = 8.0;
    private static final double THETA_MIN = -3.0;
    private static final double THETA_MAX = 3.0;

    private Chapter28_Lesson3() {}

    private static double command(double t) {
        return Math.sin(0.55 * t) + 0.35 * Math.sin(1.70 * t);
    }

    private static double regressor(double t, double u) {
        return 0.70 * Math.sin(0.90 * t) + 0.45 * u;
    }

    private static double projectedUpdate(double thetaHat, double rawUpdate) {
        if (thetaHat <= THETA_MIN && rawUpdate < 0.0) return 0.0;
        if (thetaHat >= THETA_MAX && rawUpdate > 0.0) return 0.0;
        return rawUpdate;
    }

    private static double[] derivative(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double xh1 = z[2];
        double xh2 = z[3];
        double thetaHat = z[4];

        double u = command(t);
        double phi = regressor(t, u);
        double outputError = x1 - xh1;

        return new double[] {
            -x1 + x2 + u + THETA_TRUE * phi,
            -2.0 * x1 - 3.0 * x2,
            -xh1 + xh2 + u + thetaHat * phi + outputError,
            -2.0 * xh1 - 3.0 * xh2,
            projectedUpdate(thetaHat, GAMMA * phi * outputError)
        };
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] + scale * b[i];
        }
        return result;
    }

    private static double[] rk4Step(double t, double h, double[] z) {
        double[] k1 = derivative(t, z);
        double[] k2 = derivative(t + 0.5 * h, addScaled(z, k1, 0.5 * h));
        double[] k3 = derivative(t + 0.5 * h, addScaled(z, k2, 0.5 * h));
        double[] k4 = derivative(t + h, addScaled(z, k3, h));

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + h * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
        }
        next[4] = Math.max(THETA_MIN, Math.min(THETA_MAX, next[4]));
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        final double h = 0.002;
        final double finalTime = 30.0;
        final int steps = (int) Math.round(finalTime / h);

        double[] z = {1.2, -0.8, 0.0, 0.0, -0.5};
        double sumSqE1 = 0.0;
        double sumSqE2 = 0.0;

        for (int k = 0; k < steps; k++) {
            double t = k * h;
            double e1 = z[0] - z[2];
            double e2 = z[1] - z[3];
            sumSqE1 += e1 * e1;
            sumSqE2 += e2 * e2;
            z = rk4Step(t, h, z);
        }

        double e1 = z[0] - z[2];
        double e2 = z[1] - z[3];
        double rmse1 = Math.sqrt(sumSqE1 / steps);
        double rmse2 = Math.sqrt(sumSqE2 / steps);

        System.out.println("Adaptive observer simulation");
        System.out.printf("RMSE state 1: %.6e%n", rmse1);
        System.out.printf("RMSE state 2: %.6e%n", rmse2);
        System.out.printf("Final state error: [%.6e, %.6e]%n", e1, e2);
        System.out.printf("True theta: %.6f%n", THETA_TRUE);
        System.out.printf("Final theta_hat: %.6f%n", z[4]);
        System.out.printf("Final parameter error: %.6e%n", THETA_TRUE - z[4]);

        if (!Double.isFinite(z[4])) {
            throw new IllegalStateException("Non-finite parameter estimate detected.");
        }
    }
}

MATLAB provides ode45, matrix operations, and the Control System Toolbox functions ss, lyap, and isstable. In Simulink, realize the plant and observer as parallel state-space subsystems, form the output innovation, integrate the parameter update, and use Saturation or MATLAB Function blocks for projection. Solver step size must be selected relative to both observer and adaptation bandwidths.

Chapter28_Lesson3.m

% Chapter28_Lesson3.m
% Deterministic adaptive observer using ODE45.
% Related tools: Control System Toolbox (ss, lsim, lyap) and Simulink.

clear; clc; close all;

A = [-1 1; -2 -3];
b = [1; 0];
C = [1 0];
L = [1; 0];
P = diag([1 0.5]);
Q = diag([4 3]);
Ao = A - L*C;

thetaTrue = 1.40;
gamma = 8.0;
thetaMin = -3.0;
thetaMax = 3.0;

z0 = [1.2; -0.8; 0.0; 0.0; -0.5];
tspan = linspace(0, 30, 6001);
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11, 'MaxStep', 0.02);

[t, z] = ode45(@combinedDynamics, tspan, z0, options);

x = z(:, 1:2);
xHat = z(:, 3:4);
thetaHat = z(:, 5);
e = x - xHat;
rmse = sqrt(mean(e.^2, 1));

fprintf('Adaptive observer simulation\n');
fprintf('RMSE state 1: %.6e\n', rmse(1));
fprintf('RMSE state 2: %.6e\n', rmse(2));
fprintf('Final state error: [%.6e, %.6e]\n', e(end,1), e(end,2));
fprintf('True theta: %.6f\n', thetaTrue);
fprintf('Final thetaHat: %.6f\n', thetaHat(end));
fprintf('Final parameter error: %.6e\n', thetaTrue - thetaHat(end));
fprintf('Lyapunov identity residual: %.3e\n', norm(Ao'*P + P*Ao + Q));
fprintf('SPR matching residual: %.3e\n', norm(P*b - C'));

figure;
plot(t, x(:,1), t, xHat(:,1), '--', t, x(:,2), t, xHat(:,2), '--');
grid on;
xlabel('Time [s]'); ylabel('State');
legend('x_1', 'xHat_1', 'x_2', 'xHat_2', 'Location', 'best');
title('Adaptive state observation');

figure;
plot(t, thetaHat, 'LineWidth', 1.2); hold on;
yline(thetaTrue, '--');
grid on;
xlabel('Time [s]'); ylabel('Parameter');
legend('thetaHat', 'theta true', 'Location', 'best');
title('Online parameter estimate');

% Simulink realization:
% 1. Use one State-Space block for the plant (A, b, C, 0).
% 2. Use a second State-Space block for the observer model.
% 3. Form e_y = y - yHat with a Sum block.
% 4. Add L*e_y to the observer state derivative.
% 5. Form gamma*phi*e_y, pass it through saturation/projection logic,
%    and integrate it to obtain thetaHat.
% 6. Feed thetaHat*phi into the observer input channel.

function dz = combinedDynamics(t, z)
    A = [-1 1; -2 -3];
    b = [1; 0];
    C = [1 0];
    L = [1; 0];

    thetaTrue = 1.40;
    gamma = 8.0;
    thetaMin = -3.0;
    thetaMax = 3.0;

    x = z(1:2);
    xHat = z(3:4);
    thetaHat = z(5);

    u = sin(0.55*t) + 0.35*sin(1.70*t);
    phi = 0.70*sin(0.90*t) + 0.45*u;
    y = C*x;
    yHat = C*xHat;
    outputError = y - yHat;

    xDot = A*x + b*(u + thetaTrue*phi);
    xHatDot = A*xHat + b*(u + thetaHat*phi) + L*outputError;

    rawThetaDot = gamma*phi*outputError;
    if (thetaHat <= thetaMin && rawThetaDot < 0) || ...
       (thetaHat >= thetaMax && rawThetaDot > 0)
        thetaHatDot = 0;
    else
        thetaHatDot = rawThetaDot;
    end

    dz = [xDot; xHatDot; thetaHatDot];
end

Mathematica is useful for symbolic verification of Lyapunov identities, transfer-function algebra, and high-accuracy numerical integration. The downloadable notebook is represented as a valid textual Notebook[...] expression containing one executable input cell.

Chapter28_Lesson3.nb

Notebook[{
  Cell["Chapter 28, Lesson 3: Adaptive Observers and Estimation", "Title"],
  Cell["Deterministic SPR adaptive-observer simulation", "Section"],
  Cell["ClearAll[\"Global`*\"];
A = { {-1., 1.}, {-2., -3.} }; b = {1., 0.}; c = {1., 0.}; l = {1., 0.};
p = DiagonalMatrix[{1., 0.5}]; q = DiagonalMatrix[{4., 3.}];
ao = A - Outer[Times, l, c];
thetaTrue = 1.40; gamma = 8.0; thetaMin = -3.0; thetaMax = 3.0;
u[t_] := Sin[0.55 t] + 0.35 Sin[1.70 t];
phi[t_] := 0.70 Sin[0.90 t] + 0.45 u[t];
projectedUpdate[theta_, raw_] := Which[
  theta <= thetaMin && raw < 0., 0.,
  theta >= thetaMax && raw > 0., 0.,
  True, raw
];
solution = NDSolveValue[{
  x1'[t] == -x1[t] + x2[t] + u[t] + thetaTrue phi[t],
  x2'[t] == -2 x1[t] - 3 x2[t],
  xh1'[t] == -xh1[t] + xh2[t] + u[t] + th[t] phi[t] +
    (x1[t] - xh1[t]),
  xh2'[t] == -2 xh1[t] - 3 xh2[t],
  th'[t] == projectedUpdate[
    th[t], gamma phi[t] (x1[t] - xh1[t])
  ],
  x1[0] == 1.2, x2[0] == -0.8,
  xh1[0] == 0., xh2[0] == 0., th[0] == -0.5
}, {x1, x2, xh1, xh2, th}, {t, 0, 30},
  Method -> {\"TimeIntegration\" -> {
    \"ExplicitRungeKutta\", \"DifferenceOrder\" -> 4
  } },
  MaxStepFraction -> 1/1500
];
{x1f, x2f, xh1f, xh2f, thf} = solution;
Print[\"Final state error = \",
  {x1f[30] - xh1f[30], x2f[30] - xh2f[30]}];
Print[\"Final theta estimate = \", thf[30]];
Print[\"Final parameter error = \", thetaTrue - thf[30]];
Print[\"Lyapunov identity residual norm = \",
  Norm[Transpose[ao].p + p.ao + q]];
Print[\"SPR matching residual norm = \", Norm[p.b - c]];
statePlot = Plot[
  Evaluate[{x1f[t], xh1f[t], x2f[t], xh2f[t]}], {t, 0, 30},
  PlotLegends -> {\"x1\", \"x1Hat\", \"x2\", \"x2Hat\"},
  Frame -> True, FrameLabel -> {\"Time [s]\", \"State\"},
  PlotRange -> All
];
parameterPlot = Plot[
  Evaluate[{thf[t], thetaTrue}], {t, 0, 30},
  PlotLegends -> {\"thetaHat\", \"theta true\"},
  Frame -> True, FrameLabel -> {\"Time [s]\", \"Parameter\"},
  PlotRange -> All
];
Column[{statePlot, parameterPlot}]", "Input"]
},
WindowTitle -> "Chapter28_Lesson3",
StyleDefinitions -> "Default.nb"
]        

A simulation should report more than a visually convincing plot. For the worked example, verify:

1. the algebraic residual \( \|\mathbf{A}_{o}^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_{o}+\mathbf{Q}\| \);
2. the matching residual \( \|\mathbf{P}\mathbf{b}-\mathbf{c}\| \);
3. state-error root-mean-square values;
4. the final parameter error;
5. whether projection is active for a significant time;
6. the response to reduced excitation, noise, and model mismatch.

With the nominal multisine experiment, the supplied Python, C++, and Java programs produce final state errors near machine precision and recover the true scalar parameter. That numerical result illustrates the theorem; it is not itself a proof. The proof rests on the SPR/KYP identities and the excitation assumptions.

A useful stress test is to replace \( \phi(t) \) by a signal that becomes zero. The state error still decays because the observer error matrix is stable, while the parameter estimate generally freezes at a nontrue value. This experiment demonstrates the distinction between observation and identification.

Problem 1 (Derive the Error Dynamics): For \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} +\mathbf{b}(u+\boldsymbol{\theta}^{T}\boldsymbol{\phi}) \) and the adaptive observer in Section 2, derive the state- and parameter-error equations.

Solution:

Subtract the observer equation from the plant equation. The known input cancels, and \( \boldsymbol{\theta}^{T}\boldsymbol{\phi} -\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} =\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \). Since \( y-\hat{y}=\mathbf{c}^{T}\mathbf{e} \),

\[ \dot{\mathbf{e} } =(\mathbf{A}-\mathbf{L}\mathbf{c}^{T})\mathbf{e} +\mathbf{b}\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]

Because the true parameter is constant,

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\dot{\hat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_y. \]

Problem 2 (Verify the KYP Identities): For the worked example, verify that \( \mathbf{P}=\operatorname{diag}(1,1/2) \) and \( \mathbf{Q}=\operatorname{diag}(4,3) \) satisfy the Lyapunov and matching equations.

Solution:

\[ \mathbf{A}_{o} =\begin{bmatrix}-2&1\\-2&-3\end{bmatrix}, \qquad \mathbf{P} =\begin{bmatrix}1&0\\0&1/2\end{bmatrix}. \]

Direct multiplication gives

\[ \mathbf{A}_{o}^{T}\mathbf{P} =\begin{bmatrix}-2&-1\\1&-3/2\end{bmatrix}, \qquad \mathbf{P}\mathbf{A}_{o} =\begin{bmatrix}-2&1\\-1&-3/2\end{bmatrix}. \]

Therefore

\[ \mathbf{A}_{o}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{o} =\begin{bmatrix}-4&0\\0&-3\end{bmatrix}=-\mathbf{Q}. \]

Also,

\[ \mathbf{P}\mathbf{b} =\begin{bmatrix}1&0\\0&1/2\end{bmatrix} \begin{bmatrix}1\\0\end{bmatrix} =\begin{bmatrix}1\\0\end{bmatrix}=\mathbf{c}. \]

Problem 3 (Why State Convergence Does Not Imply Parameter Convergence): Give a mathematical argument for why \( \mathbf{e}(t)\longrightarrow\mathbf{0} \) does not alone imply \( \tilde{\boldsymbol{\theta} }(t) \longrightarrow\mathbf{0} \).

Solution:

The Lyapunov derivative \( \dot{V}=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \) contains no negative-definite term in \( \tilde{\boldsymbol{\theta} } \). Thus the invariance set associated with \( \dot{V}=0 \) initially includes \( \mathbf{e}=0 \) and arbitrary bounded parameter errors. The error dynamics restricted to \( \mathbf{e}=0 \) impose

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

Without excitation, a nonzero parameter error can satisfy this equality. PE eliminates such nonzero indistinguishable directions and shrinks the largest invariant set to the origin in both state and parameter errors.

Problem 4 (Leakage and Ultimate Boundedness): Replace the ideal update law by \( \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}e_y -\sigma\boldsymbol{\Gamma}\hat{\boldsymbol{\theta} } \), where \( \sigma>0 \). Show the additional term in \( \dot{V} \) and derive an ultimate-bound inequality.

Solution:

Since \( \hat{\boldsymbol{\theta} } =\boldsymbol{\theta}-\tilde{\boldsymbol{\theta} } \),

\[ \dot{\tilde{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e_y +\sigma\boldsymbol{\Gamma} (\boldsymbol{\theta}-\tilde{\boldsymbol{\theta} }). \]

After cancellation of the state-parameter cross term,

\[ \dot{V} =-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} -2\sigma\|\tilde{\boldsymbol{\theta} }\|^{2} +2\sigma\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\theta}. \]

By Young's inequality, \( 2\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\theta} \le \|\tilde{\boldsymbol{\theta} }\|^{2} +\|\boldsymbol{\theta}\|^{2} \), hence

\[ \dot{V} \le -\lambda_{\min}(\mathbf{Q})\|\mathbf{e}\|^{2} -\sigma\|\tilde{\boldsymbol{\theta} }\|^{2} +\sigma\|\boldsymbol{\theta}\|^{2}. \]

The positive constant term prevents a proof of exact convergence but establishes uniform ultimate boundedness under standard comparison arguments.

Problem 5 (Augmented EKF Jacobian): Consider the scalar nonlinear model

\[ x_{k+1}=x_k+T_s[-a x_k+\theta u_k], \qquad \theta_{k+1}=\theta_k, \qquad y_k=x_k+v_k. \]

Form the augmented state and derive the transition and measurement Jacobians.

Solution:

Let \( \mathbf{z}_k=[x_k,\theta_k]^{T} \). Then

\[ \mathbf{F}_{a}(\mathbf{z}_k,u_k) =\begin{bmatrix} (1-aT_s)x_k+T_s\theta_k u_k\\ \theta_k \end{bmatrix}. \]

The Jacobian is

\[ \mathbf{F}_k =\frac{\partial\mathbf{F}_{a} }{\partial\mathbf{z} } =\begin{bmatrix} 1-aT_s&T_su_k\\ 0&1 \end{bmatrix}. \]

Since \( h_a(\mathbf{z}_k)=x_k \),

\[ \mathbf{H}_k =\frac{\partial h_a}{\partial\mathbf{z} } =\begin{bmatrix}1&0\end{bmatrix}. \]

The parameter is updated only through the state-parameter covariance. When \( u_k \) lacks excitation, the sensitivity term \( T_su_k \) is weak and parameter learning becomes ill-conditioned.

An adaptive observer estimates unmeasured states and unknown model parameters in one coupled dynamical system. For the linearly parameterized output-error model, an SPR transfer function supplies the KYP identity \( \mathbf{P}\mathbf{b}=\mathbf{c} \), which turns the unmeasured state cross term into a measurable output error. The gradient law then cancels the state-parameter coupling in a composite Lyapunov derivative and proves boundedness and asymptotic state convergence. Parameter convergence is stronger and requires identifiability plus persistent excitation of the relevant regressor.

Robust modifications are needed when noise, disturbances, unmodeled dynamics, sampling, and numerical limitations invalidate the ideal assumptions. Stochastic joint estimators approach the same state-parameter problem through conditional means and covariances; augmented EKF, UKF, dual filtering, and moving-horizon estimation belong to the broader Stochastic Control and Estimation curriculum. The unifying principle is that neither adaptation nor covariance tuning can replace missing information in the measured experiment.