After completing this lesson, students should be able to:

• derive the augmented state-estimation-parameter error system;
• explain why the classical linear separation principle does not automatically extend to adaptive systems;
• identify when a Lyapunov cross-term can be generated from the measurable output innovation;
• distinguish signal boundedness, state convergence, tracking convergence, and parameter convergence;
• derive ultimate bounds produced by leakage or \(\sigma\)-modification;
• recognize peaking, bursting, saturation mismatch, and unmodeled-dynamics failure mechanisms.

In a fixed linear system, a stabilizing state-feedback gain and a convergent observer can often be designed separately. In an adaptive output-feedback system, the estimated state drives the controller, the controller changes the plant trajectory, the plant trajectory supplies the regressor used for adaptation, and the parameter estimate enters the observer. The resulting loop is nonlinear and time varying even when the physical plant is linear in its state.

\[ \text{observer error} \longrightarrow \text{control error} \longrightarrow \text{regressor change} \longrightarrow \\ \text{parameter update} \longrightarrow \text{observer error}. \]

Therefore, proving that \(A-bK\) and \(A-Lc^{T}\) are Hurwitz is necessary in many designs, but it is generally not sufficient for stability of the adaptive interconnection.

Consider the uncertain SISO plant

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

where \(x\in\mathbb{R}^{n}\) is not measured, \(y\in\mathbb{R}\) is measured, \(\theta\in\mathbb{R}^{p}\) is an unknown constant parameter vector, and \(\phi(y)\in\mathbb{R}^{p}\) is a known regressor. A certainty-equivalent observer-controller is

\[ \begin{aligned} u &= -K\hat{x}+k_r r-\hat{\theta}^{T}\phi(y),\\ \dot{\hat{x} } &=A\hat{x}+b\left[u+\hat{\theta}^{T}\phi(y)\right] +L\left(y-c^{T}\hat{x}\right),\\ \hat{y}&=c^{T}\hat{x}. \end{aligned} \]

Define the observer and parameter errors by

\[ e_o=x-\hat{x}, \qquad \widetilde{\theta}=\theta-\hat{\theta}, \qquad e_y=y-\hat{y}=c^{T}e_o. \]

Subtraction of the observer from the plant gives

\[ \dot{e}_o=(A-Lc^{T})e_o+b\widetilde{\theta}^{T}\phi(y). \]

Although this equation resembles a stable observer driven by an input, that input is endogenous: it depends on the adaptive parameter error and on the plant output generated by the same closed loop.

flowchart TD
  R["Reference r"] --> C1["Adaptive controller"]
  C1 --> U["Control u"]
  U --> P["Uncertain plant"]
  P --> Y["Measured output y"]
  Y --> O["State observer"]
  O --> XH["Estimated state xhat"]
  XH --> C1
  Y --> A1["Parameter update"]
  O --> E["Output innovation y - yhat"]
  E --> A1
  A1 --> TH["Parameter estimate thetahat"]
  TH --> C1
  TH --> O
        

For a known linear plant, substituting \(u=-K\hat{x}\) yields a triangular system in suitable coordinates. The controller poles and observer poles can then be assigned independently. In the adaptive case, the augmented state includes \(\widetilde{\theta}\), and its update law depends on measured closed-loop signals.

\[ \frac{d}{dt} \begin{bmatrix} x\\ e_o\\ \widetilde{\theta} \end{bmatrix} = \begin{bmatrix} (A-bK)x+bKe_o+b\widetilde{\theta}^{T}\phi(y)+bk_r r\\ (A-Lc^{T})e_o+b\widetilde{\theta}^{T}\phi(y)\\ -\dot{\hat{\theta} }(e_y,y,\hat{x},r) \end{bmatrix}. \]

The Jacobian of this system is not block triangular in general because the regressor and adaptive law couple all three subsystems. Three distinct issues follow.

1. Observer-controller coupling: estimation error enters the control through \(bKe_o\).
2. Observer-parameter coupling: \(b\widetilde{\theta}^{T}\phi(y)\) drives the observer error.
3. Trajectory-dependent adaptation: the controller and observer change \(y\), hence change the excitation available to the parameter estimator.

Consequently, local pole placement does not by itself establish global boundedness, and convergence of the observer for a frozen parameter does not imply convergence when the parameter estimate is moving.

Let \(A_o=A-Lc^{T}\) be Hurwitz. For any \(Q=Q^{T}>0\), there is a unique \(P=P^{T}>0\) satisfying

\[ A_o^{T}P+PA_o=-Q. \]

Use the candidate

\[ V=e_o^{T}Pe_o+\widetilde{\theta}^{T}\Gamma^{-1}\widetilde{\theta}, \qquad \Gamma=\Gamma^{T}>0. \]

Its derivative is

\[ \dot{V} =-e_o^{T}Qe_o +2\widetilde{\theta}^{T}\phi(y)e_o^{T}Pb -2\widetilde{\theta}^{T}\Gamma^{-1}\dot{\hat{\theta} }. \]

Full-state adaptive control can often use \(e_o^{T}Pb\) directly. Output feedback cannot, because \(e_o\) is unknown. The cancellation becomes implementable only if the cross-term is reconstructible from the measured innovation. A sufficient scalar-output condition is

\[ Pb=c\rho, \qquad \rho\neq 0. \]

Then \(e_o^{T}Pb=\rho e_y\), so the measurable law

\[ \dot{\hat{\theta} }=\Gamma\phi(y)\rho e_y \]

produces

\[ \dot{V}=-e_o^{T}Qe_o\le 0. \]

This matching condition is restrictive. It is not a generic property of an arbitrary observer, input direction, and output map. Its importance is conceptual: output-feedback adaptation requires not merely a stable observer, but an error structure whose Lyapunov cross-term is measurable or can be generated by additional filters.

Theorem 1 (matched adaptive observer-controller, regulation).

Assume:

1. \(A_o=A-Lc^{T}\) and \(A_c=A-bK\) are Hurwitz;
2. \(P=P^{T}>0\) solves \(A_o^{T}P+PA_o=-Q\) for \(Q=Q^{T}>0\);
3. \(Pb=c\rho\) for a known nonzero scalar \(\rho\);
4. \(\phi(y)\) is locally Lipschitz and bounded whenever \(y\) is bounded;
5. the command is \(r=0\), and there is no saturation, measurement noise, or unmodeled dynamics.

Under

\[ \begin{aligned} u&=-K\hat{x}-\hat{\theta}^{T}\phi(y),\\ \dot{\hat{x} }&=A\hat{x}+b\left[u+\hat{\theta}^{T}\phi(y)\right]+Le_y,\\ \dot{\hat{\theta} }&=\Gamma\phi(y)\rho e_y, \end{aligned} \]

all closed-loop signals are bounded, \(e_o(t)\to 0\), \(\hat{x}(t)\to 0\), and \(x(t)\to 0\). Parameter convergence is not guaranteed.

Proof. From Section 4,

\[ \dot{V}=-e_o^{T}Qe_o. \]

Hence \(V(t)\le V(0)\), so \(e_o\in\mathcal{L}_{\infty}\) and \(\widetilde{\theta}\in\mathcal{L}_{\infty}\). Also,

\[ \int_{0}^{\infty}\|e_o(t)\|^{2}dt \le \frac{V(0)}{\lambda_{\min}(Q)}, \]

so \(e_o\in\mathcal{L}_{2}\). Substitution of the control law into the observer equation gives the triangular relation

\[ \dot{\hat{x} }=(A-bK)\hat{x}+L c^{T}e_o. \]

Since \(A-bK\) is Hurwitz and \(e_o\in\mathcal{L}_{\infty}\), the observer state is bounded. Therefore \(x=\hat{x}+e_o\) and \(y=c^{T}x\) are bounded; the assumed property of \(\phi\) implies a bounded regressor. The error equation then gives bounded \(\dot e_o\). Barbalat's lemma yields \(e_o(t)\to 0\). The stable linear system for \(\hat{x}\) is driven by an input tending to zero, hence \(\hat{x}(t)\to 0\). Finally, \(x=\hat{x}+e_o\to 0\). The derivative contains no negative-definite term in \(\widetilde{\theta}\), so the proof does not establish \(\widetilde{\theta}\to 0\). \(\square\)

Suppose the exact condition is replaced by

\[ Pb=c\rho+\delta, \qquad \delta\neq 0. \]

Using the same output-error adaptation law leaves

\[ \dot{V}=-e_o^{T}Qe_o +2\widetilde{\theta}^{T}\phi(y)e_o^{T}\delta. \]

Young's inequality gives

\[ 2\left|\widetilde{\theta}^{T}\phi(y)e_o^{T}\delta\right| \le \frac{\lambda_{\min}(Q)}{2}\|e_o\|^{2} +\frac{2\|\delta\|^{2} }{\lambda_{\min}(Q)} \|\widetilde{\theta}\|^{2}\|\phi(y)\|^{2}. \]

Thus

\[ \dot{V} \le -\frac{\lambda_{\min}(Q)}{2}\|e_o\|^{2} +\frac{2\|\delta\|^{2} }{\lambda_{\min}(Q)} \|\widetilde{\theta}\|^{2}\|\phi(y)\|^{2}. \]

There is no negative parameter-error term to absorb the residual. Consequently, an arbitrarily small structural mismatch cannot simply be dismissed as a small perturbation unless an additional robustness mechanism supplies damping, a known parameter bound, or an input-to-state stability estimate.

Output-feedback adaptive designs often use a high-gain observer to make estimation dynamics faster than plant dynamics. In observer canonical coordinates, gains may contain inverse powers of a small parameter \(\varepsilon\):

\[ L_{\varepsilon} =\begin{bmatrix} \ell_1/\varepsilon & \ell_2/\varepsilon^{2} & \cdots & \ell_n/\varepsilon^{n} \end{bmatrix}^{T}. \]

Fast asymptotic convergence does not imply a small transient. With inconsistent initial conditions, internal state estimates may exhibit peaking whose magnitude increases as \(\varepsilon\) decreases. Since \(u=-K\hat{x}+\cdots\), observer peaking can produce a large control pulse. If the adaptive law uses the same transient signals, the parameter estimate may move significantly before the observer enters its asymptotic regime.

\[ \|D_{\varepsilon}^{-1}e_o(t)\| \le M e^{-\alpha t/\varepsilon} \|D_{\varepsilon}^{-1}e_o(0)\| +\text{forced-response terms}, \]

where the scaling matrix \(D_{\varepsilon}\) contains powers of \(\varepsilon\). A practical design therefore combines high-gain observation with control saturation, parameter projection, staged adaptation, or a gain schedule that starts with a moderate observer bandwidth.

Theorem 1 proves convergence of the estimation error but not necessarily of the parameter error. A standard sufficient condition for parameter convergence is persistent excitation of the effective regressor:

\[ \exists\,T>0,\;\alpha>0 \quad\text{such that}\quad \int_{t}^{t+T}\phi(\tau)\phi^{T}(\tau)d\tau \ge \alpha I, \qquad \forall t\ge 0. \]

Regulation often destroys excitation because the output and regressor approach zero. A nonzero parameter error can then remain hidden while the system is near equilibrium. Later disturbances or reference changes can reactivate the poorly identified direction and produce a burst of oscillation. This is why a quiet parameter trace is not evidence that the estimate is correct.

Measurement bias creates a second problem. If the innovation contains a nonzero mean component, a pure integral adaptation law can accumulate the bias and cause parameter drift even when the plant state is small.

9.1 Sigma-modification

Replace the ideal law by

\[ \dot{\hat{\theta} } =\Gamma\left[\phi(y)\rho e_y -\sigma\left(\hat{\theta}-\theta_0\right)\right], \qquad \sigma>0. \]

Under exact matching and no other disturbances,

\[ \begin{aligned} \dot{V} &=-e_o^{T}Qe_o +2\sigma\widetilde{\theta}^{T} \left(\hat{\theta}-\theta_0\right)\\ &=-e_o^{T}Qe_o -2\sigma\|\widetilde{\theta}\|^{2} +2\sigma\widetilde{\theta}^{T}(\theta-\theta_0)\\ &\le -\lambda_{\min}(Q)\|e_o\|^{2} -\sigma\|\widetilde{\theta}\|^{2} +\sigma\|\theta-\theta_0\|^{2}. \end{aligned} \]

The estimate is pulled toward \(\theta_0\). Exact parameter convergence is sacrificed, but drift is prevented and the combined error becomes uniformly ultimately bounded.

9.2 Projection

If prior knowledge gives a compact convex parameter set \(\Omega\), use

\[ \dot{\hat{\theta} } =\operatorname{Proj}_{\Omega} \left(\hat{\theta},\Gamma\phi(y)\rho e_y\right). \]

The projection operator is selected so that, for the true parameter in \(\Omega\), its correction does not increase the standard Lyapunov derivative. Projection provides a hard parameter bound but does not by itself remove measurement-noise sensitivity inside the admissible set.

9.3 Dead zone and normalization

\[ \dot{\hat{\theta} }= \begin{cases} \Gamma\dfrac{\phi(y)\rho e_y}{1+\nu\|\phi(y)\|^{2} }, & |e_y|>\varepsilon_d,\\ 0, & |e_y|\le\varepsilon_d. \end{cases} \]

The dead zone prevents adaptation to innovations indistinguishable from noise, while normalization limits updates under large regressors. Their thresholds must be chosen using noise and modeling-error bounds rather than solely by trial and error.

Suppose the commanded input is \(u_c\), but the plant receives

\[ u=\operatorname{sat}(u_c), \qquad \Delta_u=u-u_c. \]

If the observer uses \(u_c\) while the plant receives \(u\), the error equation acquires the unmatched term \(b\Delta_u\):

\[ \dot{e}_o=A_o e_o+b\widetilde{\theta}^{T}\phi(y)+b\Delta_u. \]

Then

\[ \dot{V}=-e_o^{T}Qe_o+2e_o^{T}Pb\Delta_u \]

even when the parameter cross-term is cancelled. The observer should use the actual applied actuator signal whenever it is available. If it is not available, an anti-windup or saturation-error model is required, and the claimed stability result should be weakened to an ultimate-bound or regional result.

flowchart TD
  S["Start with plant, measured output, and uncertainty model"] --> O["Choose observer structure"]
  O --> H["Verify nominal observer matrix is Hurwitz"]
  H --> C["Choose nominal controller and verify closed-loop matrix is Hurwitz"]
  C --> M["Derive exact augmented error dynamics"]
  M --> X["Can the Lyapunov cross-term be computed from measured signals?"]
  X -->|"No"| F["Add filters, dynamic extension, or redesign parameterization"]
  X -->|"Yes"| R["Add projection, leakage, normalization, and dead zone as needed"]
  F --> R
  R --> P["Check observer peaking and actuator saturation"]
  P --> U["Include noise and unmodeled dynamics in the proof or bound"]
  U --> E["State the actual result: asymptotic, regional, or ultimate boundedness"]
  E --> V["Validate with worst-case initial conditions and bandwidth sweeps"]
        

A rigorous design report should explicitly identify the domain of attraction, assumed parameter set, noise class, saturation model, required excitation, and whether the result is asymptotic stability, uniform asymptotic stability, input-to-state stability, or uniform ultimate boundedness.

The downloadable programs use

\[ A=\begin{bmatrix}0&1\\-2&-0.8\end{bmatrix}, \quad c=\begin{bmatrix}1\\0\end{bmatrix}, \quad L=\begin{bmatrix}4\\4\end{bmatrix}. \]

For \(Q=I\), the Lyapunov equation gives approximately

\[ P=\begin{bmatrix} 0.51902174&-0.26268116\\ -0.26268116&0.29664855 \end{bmatrix}. \]

Choose the input direction as

\[ b=P^{-1}c \approx \begin{bmatrix}3.49139367\\3.09161577\end{bmatrix}, \]

so that \(Pb=c\) by construction. The feedback gain

\[ K=\begin{bmatrix}0.65194504&0.29880919\end{bmatrix} \]

places the eigenvalues of \(A-bK\) at approximately \(-1.5\) and \(-2.5\). The plant uses \(\theta=1.2\) and \(\phi(y)=\tanh(y)\). Two cases are compared:

1. Ideal matched case: no noise, no leakage, and effectively no saturation.
2. Noisy robust case: high-frequency measurement noise, control limiting, and \(\sigma\)-modification.

The ideal simulation drives the plant and observer states to zero, but the estimate need not converge to \(1.2\). This is the expected non-PE behavior, not a numerical error.

This implementation uses NumPy, Matplotlib, and a from-scratch fourth-order Runge-Kutta integrator. It verifies the Lyapunov and matching residuals, writes CSV files, and plots state estimates, parameter estimates, and control.

Chapter15_Lesson5.py

"""Chapter15_Lesson5.py

Output-feedback adaptive-control stability laboratory.

The example implements a structured adaptive observer/controller for

    x_dot = A x + b [u + theta * tanh(y)],    y = c^T x,

with an observer satisfying the measurable matching condition P b = c.
It compares an ideal case with a noisy, saturated, sigma-modified case.
Only NumPy and Matplotlib are required.
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


A = np.array([[0.0, 1.0], [-2.0, -0.8]], dtype=float)
C = np.array([1.0, 0.0], dtype=float)
L = np.array([4.0, 4.0], dtype=float)
P = np.array(
    [[0.5190217391304348, -0.2626811594202899],
     [-0.2626811594202899, 0.2966485507246377]],
    dtype=float,
)
B = np.array([3.49139367, 3.09161577], dtype=float)  # approximately P^{-1} C
K = np.array([0.65194504, 0.29880919], dtype=float)
TRUE_THETA = 1.2


@dataclass(frozen=True)
class Case:
    name: str
    gamma: float
    sigma: float
    noise_amplitude: float
    control_limit: float


def phi(value: float) -> float:
    """Bounded regressor used by plant, observer, and adaptive controller."""
    return math.tanh(value)


def measurement_noise(t: float, amplitude: float) -> float:
    """Deterministic high-frequency measurement disturbance."""
    return amplitude * (math.sin(37.0 * t) + 0.5 * math.sin(91.0 * t))


def saturate(value: float, limit: float) -> float:
    return max(-limit, min(limit, value))


def dynamics(t: float, z: np.ndarray, case: Case) -> np.ndarray:
    """Combined plant-observer-parameter dynamics.

    z = [x1, x2, xhat1, xhat2, theta_hat].
    """
    x = z[0:2]
    x_hat = z[2:4]
    theta_hat = float(z[4])

    y = float(C @ x)
    y_measured = y + measurement_noise(t, case.noise_amplitude)
    y_hat = float(C @ x_hat)
    innovation = y_measured - y_hat

    regressor_measured = phi(y_measured)
    control_unsaturated = -float(K @ x_hat) - theta_hat * regressor_measured
    control = saturate(control_unsaturated, case.control_limit)

    x_dot = A @ x + B * (control + TRUE_THETA * phi(y))
    x_hat_dot = (
        A @ x_hat
        + B * (control + theta_hat * regressor_measured)
        + L * innovation
    )

    # The first term is the ideal Lyapunov cancellation because P B = C.
    # The second is sigma-modification, producing ultimate boundedness under noise.
    theta_hat_dot = (
        case.gamma * regressor_measured * innovation
        - case.sigma * case.gamma * theta_hat
    )

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


def rk4_step(
    f: Callable[[float, np.ndarray, Case], np.ndarray],
    t: float,
    z: np.ndarray,
    step: float,
    case: Case,
) -> np.ndarray:
    k1 = f(t, z, case)
    k2 = f(t + 0.5 * step, z + 0.5 * step * k1, case)
    k3 = f(t + 0.5 * step, z + 0.5 * step * k2, case)
    k4 = f(t + step, z + step * k3, case)
    return z + (step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(case: Case, final_time: float = 20.0, step: float = 0.001) -> np.ndarray:
    count = int(round(final_time / step)) + 1
    history = np.zeros((count, 9), dtype=float)
    z = np.array([1.0, -0.5, 0.0, 0.0, 0.0], dtype=float)

    for index in range(count):
        t = index * step
        x = z[0:2]
        x_hat = z[2:4]
        y = float(C @ x)
        y_measured = y + measurement_noise(t, case.noise_amplitude)
        innovation = y_measured - float(C @ x_hat)
        control_raw = -float(K @ x_hat) - float(z[4]) * phi(y_measured)
        control = saturate(control_raw, case.control_limit)

        history[index, :] = [
            t,
            x[0],
            x[1],
            x_hat[0],
            x_hat[1],
            z[4],
            innovation,
            control,
            control_raw,
        ]
        if index + 1 < count:
            z = rk4_step(dynamics, t, z, step, case)

    return history


def report(case: Case, history: np.ndarray) -> None:
    final = history[-1]
    peak_control = float(np.max(np.abs(history[:, 7])))
    peak_unsaturated = float(np.max(np.abs(history[:, 8])))
    rms_innovation = float(np.sqrt(np.mean(history[:, 6] ** 2)))
    print(f"\n{case.name}")
    print(f"  final plant state       = [{final[1]: .6f}, {final[2]: .6f}]")
    print(f"  final observer state    = [{final[3]: .6f}, {final[4]: .6f}]")
    print(f"  final theta_hat         = {final[5]: .6f}")
    print(f"  RMS output innovation   = {rms_innovation: .6f}")
    print(f"  peak applied control    = {peak_control: .6f}")
    print(f"  peak requested control  = {peak_unsaturated: .6f}")


def save_results(case: Case, history: np.ndarray, output_directory: Path) -> None:
    output_directory.mkdir(parents=True, exist_ok=True)
    safe_name = case.name.lower().replace(" ", "_").replace("-", "_")
    csv_path = output_directory / f"Chapter15_Lesson5_{safe_name}.csv"
    np.savetxt(
        csv_path,
        history,
        delimiter=",",
        header="t,x1,x2,xhat1,xhat2,theta_hat,innovation,u,u_raw",
        comments="",
    )


def plot_cases(results: list[tuple[Case, np.ndarray]], output_directory: Path) -> None:
    output_directory.mkdir(parents=True, exist_ok=True)

    plt.figure(figsize=(10, 6))
    for case, history in results:
        plt.plot(history[:, 0], history[:, 1], label=f"x1: {case.name}")
        plt.plot(history[:, 0], history[:, 3], "--", label=f"xhat1: {case.name}")
    plt.xlabel("time [s]")
    plt.ylabel("measured state and estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_directory / "Chapter15_Lesson5_output_estimation.png", dpi=180)

    plt.figure(figsize=(10, 6))
    for case, history in results:
        plt.plot(history[:, 0], history[:, 5], label=case.name)
    plt.axhline(TRUE_THETA, linestyle=":", label="true theta")
    plt.xlabel("time [s]")
    plt.ylabel("parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_directory / "Chapter15_Lesson5_parameter_estimate.png", dpi=180)

    plt.figure(figsize=(10, 6))
    for case, history in results:
        plt.plot(history[:, 0], history[:, 7], label=case.name)
    plt.xlabel("time [s]")
    plt.ylabel("applied control")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_directory / "Chapter15_Lesson5_control.png", dpi=180)
    plt.show()


def verify_design_identities() -> None:
    observer_matrix = A - np.outer(L, C)
    controller_matrix = A - np.outer(B, K)
    lyapunov_residual = observer_matrix.T @ P + P @ observer_matrix + np.eye(2)
    matching_residual = P @ B - C

    print("Observer eigenvalues:", np.linalg.eigvals(observer_matrix))
    print("Controller eigenvalues:", np.linalg.eigvals(controller_matrix))
    print("||A_o^T P + P A_o + I||_2 =", np.linalg.norm(lyapunov_residual, 2))
    print("||P B - C||_2             =", np.linalg.norm(matching_residual, 2))


def main() -> None:
    verify_design_identities()

    cases = [
        Case(
            name="Ideal matched case",
            gamma=8.0,
            sigma=0.0,
            noise_amplitude=0.0,
            control_limit=1.0e9,
        ),
        Case(
            name="Noisy robust case",
            gamma=8.0,
            sigma=0.03,
            noise_amplitude=0.01,
            control_limit=8.0,
        ),
    ]

    output_directory = Path("Chapter15_Lesson5_results")
    results: list[tuple[Case, np.ndarray]] = []
    for case in cases:
        history = simulate(case)
        report(case, history)
        save_results(case, history, output_directory)
        results.append((case, history))

    plot_cases(results, output_directory)


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. It implements fixed-size vectors, RK4 integration, saturation, deterministic sensor noise, CSV output, and finite-value checks.

Chapter15_Lesson5.cpp

// Chapter15_Lesson5.cpp
// Output-feedback adaptive-control stability laboratory.
// Compile: g++ -std=c++17 -O2 Chapter15_Lesson5.cpp -o Chapter15_Lesson5

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

namespace {

using State = std::array<double, 5>;
using Vector2 = std::array<double, 2>;

constexpr std::array<std::array<double, 2>, 2> A{ {
    { {0.0, 1.0} },
    { {-2.0, -0.8} },
} };
constexpr Vector2 C{ {1.0, 0.0} };
constexpr Vector2 L{ {4.0, 4.0} };
constexpr Vector2 B{ {3.49139367, 3.09161577} };
constexpr Vector2 K{ {0.65194504, 0.29880919} };
constexpr double TRUE_THETA = 1.2;

struct Case {
    std::string name;
    double gamma;
    double sigma;
    double noiseAmplitude;
    double controlLimit;
};

double dot(const Vector2& first, const Vector2& second) {
    return first[0] * second[0] + first[1] * second[1];
}

Vector2 matVec(const std::array<std::array<double, 2>, 2>& matrix,
               const Vector2& vector) {
    return { {
        matrix[0][0] * vector[0] + matrix[0][1] * vector[1],
        matrix[1][0] * vector[0] + matrix[1][1] * vector[1],
    } };
}

double phi(double value) {
    return std::tanh(value);
}

double measurementNoise(double time, double amplitude) {
    return amplitude * (std::sin(37.0 * time) + 0.5 * std::sin(91.0 * time));
}

double saturate(double value, double limit) {
    return std::max(-limit, std::min(limit, value));
}

State addScaled(const State& state, const State& increment, double scale) {
    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = state[i] + scale * increment[i];
    }
    return result;
}

State dynamics(double time, const State& state, const Case& simulationCase) {
    const Vector2 x{ {state[0], state[1]} };
    const Vector2 xHat{ {state[2], state[3]} };
    const double thetaHat = state[4];

    const double y = dot(C, x);
    const double yMeasured = y + measurementNoise(time, simulationCase.noiseAmplitude);
    const double yHat = dot(C, xHat);
    const double innovation = yMeasured - yHat;
    const double regressorMeasured = phi(yMeasured);

    const double rawControl = -dot(K, xHat) - thetaHat * regressorMeasured;
    const double control = saturate(rawControl, simulationCase.controlLimit);

    Vector2 xDot = matVec(A, x);
    Vector2 xHatDot = matVec(A, xHat);
    for (std::size_t i = 0; i < 2; ++i) {
        xDot[i] += B[i] * (control + TRUE_THETA * phi(y));
        xHatDot[i] += B[i] * (control + thetaHat * regressorMeasured)
                      + L[i] * innovation;
    }

    const double thetaHatDot =
        simulationCase.gamma * regressorMeasured * innovation
        - simulationCase.sigma * simulationCase.gamma * thetaHat;

    return { {xDot[0], xDot[1], xHatDot[0], xHatDot[1], thetaHatDot} };
}

State rk4Step(double time, const State& state, double step, const Case& simulationCase) {
    const State k1 = dynamics(time, state, simulationCase);
    const State k2 = dynamics(time + 0.5 * step,
                              addScaled(state, k1, 0.5 * step), simulationCase);
    const State k3 = dynamics(time + 0.5 * step,
                              addScaled(state, k2, 0.5 * step), simulationCase);
    const State k4 = dynamics(time + step, addScaled(state, k3, step), simulationCase);

    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = state[i] + (step / 6.0)
            * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(result[i])) {
            throw std::runtime_error("Non-finite state encountered.");
        }
    }
    return result;
}

void simulate(const Case& simulationCase,
              const std::string& fileName,
              double finalTime = 20.0,
              double step = 0.001) {
    std::ofstream output(fileName);
    if (!output) {
        throw std::runtime_error("Cannot open output file: " + fileName);
    }

    output << "t,x1,x2,xhat1,xhat2,theta_hat,innovation,u,u_raw\n";
    output << std::setprecision(12);

    State state{ {1.0, -0.5, 0.0, 0.0, 0.0} };
    const int sampleCount = static_cast<int>(std::llround(finalTime / step)) + 1;
    double peakControl = 0.0;
    double innovationSquaredSum = 0.0;

    for (int index = 0; index < sampleCount; ++index) {
        const double time = index * step;
        const Vector2 x{ {state[0], state[1]} };
        const Vector2 xHat{ {state[2], state[3]} };
        const double y = dot(C, x);
        const double yMeasured = y + measurementNoise(time, simulationCase.noiseAmplitude);
        const double innovation = yMeasured - dot(C, xHat);
        const double rawControl = -dot(K, xHat) - state[4] * phi(yMeasured);
        const double control = saturate(rawControl, simulationCase.controlLimit);

        peakControl = std::max(peakControl, std::abs(control));
        innovationSquaredSum += innovation * innovation;

        output << time << ',' << state[0] << ',' << state[1] << ','
               << state[2] << ',' << state[3] << ',' << state[4] << ','
               << innovation << ',' << control << ',' << rawControl << '\n';

        if (index + 1 < sampleCount) {
            state = rk4Step(time, state, step, simulationCase);
        }
    }

    const double rmsInnovation = std::sqrt(innovationSquaredSum / sampleCount);
    std::cout << "\n" << simulationCase.name << "\n"
              << "  final plant state     = [" << state[0] << ", " << state[1] << "]\n"
              << "  final observer state  = [" << state[2] << ", " << state[3] << "]\n"
              << "  final theta_hat       = " << state[4] << "\n"
              << "  RMS innovation        = " << rmsInnovation << "\n"
              << "  peak applied control  = " << peakControl << "\n";
}

}  // namespace

int main() {
    try {
        const Case ideal{"Ideal matched case", 8.0, 0.0, 0.0, 1.0e9};
        const Case robust{"Noisy robust case", 8.0, 0.03, 0.01, 8.0};
        simulate(ideal, "Chapter15_Lesson5_ideal_cpp.csv");
        simulate(robust, "Chapter15_Lesson5_robust_cpp.csv");
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses a record for simulation settings and standard Java file I/O. It produces trajectories numerically consistent with the Python and C++ programs.

Chapter15_Lesson5.java

// Chapter15_Lesson5.java
// Output-feedback adaptive-control stability laboratory.
// Compile: javac Chapter15_Lesson5.java
// Run:     java Chapter15_Lesson5

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 Chapter15_Lesson5 {
    private static final double[][] A = {
        {0.0, 1.0},
        {-2.0, -0.8}
    };
    private static final double[] C = {1.0, 0.0};
    private static final double[] L = {4.0, 4.0};
    private static final double[] B = {3.49139367, 3.09161577};
    private static final double[] K = {0.65194504, 0.29880919};
    private static final double TRUE_THETA = 1.2;

    private record SimulationCase(
        String name,
        double gamma,
        double sigma,
        double noiseAmplitude,
        double controlLimit
    ) {}

    private Chapter15_Lesson5() {}

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

    private static double[] matrixVector(double[][] matrix, double[] vector) {
        return new double[] {
            matrix[0][0] * vector[0] + matrix[0][1] * vector[1],
            matrix[1][0] * vector[0] + matrix[1][1] * vector[1]
        };
    }

    private static double phi(double value) {
        return Math.tanh(value);
    }

    private static double measurementNoise(double time, double amplitude) {
        return amplitude * (Math.sin(37.0 * time) + 0.5 * Math.sin(91.0 * time));
    }

    private static double saturate(double value, double limit) {
        return Math.max(-limit, Math.min(limit, value));
    }

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

    private static double[] dynamics(double time, double[] state, SimulationCase simulationCase) {
        double[] x = {state[0], state[1]};
        double[] xHat = {state[2], state[3]};
        double thetaHat = state[4];

        double y = dot(C, x);
        double yMeasured = y + measurementNoise(time, simulationCase.noiseAmplitude());
        double yHat = dot(C, xHat);
        double innovation = yMeasured - yHat;
        double regressorMeasured = phi(yMeasured);

        double rawControl = -dot(K, xHat) - thetaHat * regressorMeasured;
        double control = saturate(rawControl, simulationCase.controlLimit());

        double[] xDot = matrixVector(A, x);
        double[] xHatDot = matrixVector(A, xHat);
        for (int i = 0; i < 2; i++) {
            xDot[i] += B[i] * (control + TRUE_THETA * phi(y));
            xHatDot[i] += B[i] * (control + thetaHat * regressorMeasured)
                + L[i] * innovation;
        }

        double thetaHatDot = simulationCase.gamma() * regressorMeasured * innovation
            - simulationCase.sigma() * simulationCase.gamma() * thetaHat;

        return new double[] {xDot[0], xDot[1], xHatDot[0], xHatDot[1], thetaHatDot};
    }

    private static double[] rk4Step(
        double time,
        double[] state,
        double step,
        SimulationCase simulationCase
    ) {
        double[] k1 = dynamics(time, state, simulationCase);
        double[] k2 = dynamics(
            time + 0.5 * step,
            addScaled(state, k1, 0.5 * step),
            simulationCase
        );
        double[] k3 = dynamics(
            time + 0.5 * step,
            addScaled(state, k2, 0.5 * step),
            simulationCase
        );
        double[] k4 = dynamics(
            time + step,
            addScaled(state, k3, step),
            simulationCase
        );

        double[] result = new double[state.length];
        for (int i = 0; i < result.length; i++) {
            result[i] = state[i] + (step / 6.0)
                * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(result[i])) {
                throw new IllegalStateException("Non-finite state encountered.");
            }
        }
        return result;
    }

    private static void simulate(
        SimulationCase simulationCase,
        String fileName,
        double finalTime,
        double step
    ) throws IOException {
        int sampleCount = (int) Math.round(finalTime / step) + 1;
        double[] state = {1.0, -0.5, 0.0, 0.0, 0.0};
        double innovationSquaredSum = 0.0;
        double peakControl = 0.0;

        try (BufferedWriter output = Files.newBufferedWriter(Path.of(fileName))) {
            output.write("t,x1,x2,xhat1,xhat2,theta_hat,innovation,u,u_raw\n");

            for (int index = 0; index < sampleCount; index++) {
                double time = index * step;
                double[] x = {state[0], state[1]};
                double[] xHat = {state[2], state[3]};
                double y = dot(C, x);
                double yMeasured = y + measurementNoise(time, simulationCase.noiseAmplitude());
                double innovation = yMeasured - dot(C, xHat);
                double rawControl = -dot(K, xHat) - state[4] * phi(yMeasured);
                double control = saturate(rawControl, simulationCase.controlLimit());

                innovationSquaredSum += innovation * innovation;
                peakControl = Math.max(peakControl, Math.abs(control));

                output.write(String.format(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    time, state[0], state[1], state[2], state[3], state[4],
                    innovation, control, rawControl
                ));

                if (index + 1 < sampleCount) {
                    state = rk4Step(time, state, step, simulationCase);
                }
            }
        }

        double rmsInnovation = Math.sqrt(innovationSquaredSum / sampleCount);
        System.out.printf(Locale.US, "%n%s%n", simulationCase.name());
        System.out.printf(Locale.US, "  final plant state     = [%.6f, %.6f]%n", state[0], state[1]);
        System.out.printf(Locale.US, "  final observer state  = [%.6f, %.6f]%n", state[2], state[3]);
        System.out.printf(Locale.US, "  final theta_hat       = %.6f%n", state[4]);
        System.out.printf(Locale.US, "  RMS innovation        = %.6f%n", rmsInnovation);
        System.out.printf(Locale.US, "  peak applied control  = %.6f%n", peakControl);
    }

    public static void main(String[] args) {
        SimulationCase ideal = new SimulationCase(
            "Ideal matched case", 8.0, 0.0, 0.0, 1.0e9
        );
        SimulationCase robust = new SimulationCase(
            "Noisy robust case", 8.0, 0.03, 0.01, 8.0
        );

        try {
            simulate(ideal, "Chapter15_Lesson5_ideal_java.csv", 20.0, 0.001);
            simulate(robust, "Chapter15_Lesson5_robust_java.csv", 20.0, 0.001);
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script uses base MATLAB functions and local functions in one file. A corresponding Simulink realization can be assembled from two vector Integrator blocks for \(x\) and \(\hat{x}\), one scalar Integrator for \(\hat{\theta}\), Matrix Gain blocks for \(A,b,K,L,c^{T}\), a Trigonometric Function block set to tanh, a Saturation block, and Sum/Product blocks implementing the innovation and adaptation law. The observer must receive the actual saturated control signal, not merely the unsaturated command.

Chapter15_Lesson5.m

% Chapter15_Lesson5.m
% Output-feedback adaptive-control stability laboratory.
% The script compares an ideal matched design with a noisy, saturated,
% sigma-modified implementation. It requires only base MATLAB.

clear; clc; close all;

A = [0, 1; -2, -0.8];
C = [1, 0];
L = [4; 4];
P = [0.5190217391304348, -0.2626811594202899; ...
    -0.2626811594202899, 0.2966485507246377];
B = [3.49139367; 3.09161577];
K = [0.65194504, 0.29880919];
thetaTrue = 1.2;

Ao = A - L*C;
Ac = A - B*K;
observerPoles = eig(Ao);
controllerPoles = eig(Ac);
fprintf('Observer eigenvalues: %g%+gi, %g%+gi\n', ...
    real(observerPoles(1)), imag(observerPoles(1)), ...
    real(observerPoles(2)), imag(observerPoles(2)));
fprintf('Matching residual norm: %.3e\n', norm(P*B - C'));
fprintf('Lyapunov residual norm: %.3e\n', norm(Ao'*P + P*Ao + eye(2)));
fprintf('Controller poles: %g, %g\n', controllerPoles(1), controllerPoles(2));

cases(1) = struct('name', 'Ideal matched case', ...
    'gamma', 8.0, 'sigma', 0.0, 'noiseAmplitude', 0.0, 'controlLimit', 1.0e9);
cases(2) = struct('name', 'Noisy robust case', ...
    'gamma', 8.0, 'sigma', 0.03, 'noiseAmplitude', 0.01, 'controlLimit', 8.0);

finalTime = 20.0;
step = 0.001;
time = (0:step:finalTime)';
results = cell(numel(cases), 1);

for caseIndex = 1:numel(cases)
    currentCase = cases(caseIndex);
    state = [1.0; -0.5; 0.0; 0.0; 0.0];
    history = zeros(numel(time), 9);

    for index = 1:numel(time)
        t = time(index);
        x = state(1:2);
        xHat = state(3:4);
        thetaHat = state(5);

        y = C*x;
        yMeasured = y + measurementNoise(t, currentCase.noiseAmplitude);
        innovation = yMeasured - C*xHat;
        regressorMeasured = tanh(yMeasured);
        rawControl = -K*xHat - thetaHat*regressorMeasured;
        control = saturate(rawControl, currentCase.controlLimit);

        history(index, :) = [t, x', xHat', thetaHat, innovation, control, rawControl];

        if index < numel(time)
            state = rk4Step(@combinedDynamics, t, state, step, currentCase, ...
                A, C, L, B, K, thetaTrue);
        end
    end

    results{caseIndex} = history;
    rmsInnovation = sqrt(mean(history(:, 7).^2));
    fprintf('\n%s\n', currentCase.name);
    fprintf('  final plant state     = [% .6f, % .6f]\n', history(end, 2), history(end, 3));
    fprintf('  final observer state  = [% .6f, % .6f]\n', history(end, 4), history(end, 5));
    fprintf('  final theta_hat       = % .6f\n', history(end, 6));
    fprintf('  RMS innovation        = % .6f\n', rmsInnovation);
    fprintf('  peak applied control  = % .6f\n', max(abs(history(:, 8))));

    safeName = lower(strrep(strrep(currentCase.name, ' ', '_'), '-', '_'));
    writematrix(history, ['Chapter15_Lesson5_', safeName, '.csv']);
end

figure;
for caseIndex = 1:numel(cases)
    history = results{caseIndex};
    plot(history(:, 1), history(:, 2), 'DisplayName', ['x_1: ', cases(caseIndex).name]);
    hold on;
    plot(history(:, 1), history(:, 4), '--', ...
        'DisplayName', ['xhat_1: ', cases(caseIndex).name]);
end
xlabel('time [s]'); ylabel('output state and estimate'); grid on; legend('Location', 'best');

figure;
for caseIndex = 1:numel(cases)
    history = results{caseIndex};
    plot(history(:, 1), history(:, 6), 'DisplayName', cases(caseIndex).name);
    hold on;
end
yline(thetaTrue, ':', 'true theta');
xlabel('time [s]'); ylabel('parameter estimate'); grid on; legend('Location', 'best');

figure;
for caseIndex = 1:numel(cases)
    history = results{caseIndex};
    plot(history(:, 1), history(:, 8), 'DisplayName', cases(caseIndex).name);
    hold on;
end
xlabel('time [s]'); ylabel('applied control'); grid on; legend('Location', 'best');

function derivative = combinedDynamics(t, state, currentCase, A, C, L, B, K, thetaTrue)
    x = state(1:2);
    xHat = state(3:4);
    thetaHat = state(5);

    y = C*x;
    yMeasured = y + measurementNoise(t, currentCase.noiseAmplitude);
    yHat = C*xHat;
    innovation = yMeasured - yHat;
    regressorMeasured = tanh(yMeasured);

    rawControl = -K*xHat - thetaHat*regressorMeasured;
    control = saturate(rawControl, currentCase.controlLimit);

    xDot = A*x + B*(control + thetaTrue*tanh(y));
    xHatDot = A*xHat + B*(control + thetaHat*regressorMeasured) + L*innovation;
    thetaHatDot = currentCase.gamma*regressorMeasured*innovation ...
        - currentCase.sigma*currentCase.gamma*thetaHat;

    derivative = [xDot; xHatDot; thetaHatDot];
end

function nextState = rk4Step(functionHandle, t, state, step, currentCase, A, C, L, B, K, thetaTrue)
    k1 = functionHandle(t, state, currentCase, A, C, L, B, K, thetaTrue);
    k2 = functionHandle(t + step/2, state + step*k1/2, currentCase, A, C, L, B, K, thetaTrue);
    k3 = functionHandle(t + step/2, state + step*k2/2, currentCase, A, C, L, B, K, thetaTrue);
    k4 = functionHandle(t + step, state + step*k3, currentCase, A, C, L, B, K, thetaTrue);
    nextState = state + step*(k1 + 2*k2 + 2*k3 + k4)/6;
end

function value = measurementNoise(t, amplitude)
    value = amplitude*(sin(37*t) + 0.5*sin(91*t));
end

function value = saturate(inputValue, limit)
    value = max(-limit, min(limit, inputValue));
end

The notebook uses NDSolveValue, verifies the matrix identities, simulates both cases, and constructs output, parameter, and control plots.

Chapter15_Lesson5.nb

Notebook[{
 Cell["Chapter 15, Lesson 5: Stability Challenges in Output-Feedback Adaptive Control", "Title"],
 Cell["This notebook simulates the structured output-feedback adaptive controller used in the lesson. The ideal case satisfies P.b = c, while the robust case adds measurement noise, control saturation, and sigma-modification.", "Text"],
 Cell[BoxData @ ToBoxes[
   Defer[CompoundExpression[
    ClearAll["Global`*"],
    a = { {0., 1.}, {-2., -0.8} },
    c = {1., 0.},
    l = {4., 4.},
    p = { {0.5190217391304348, -0.2626811594202899},
         {-0.2626811594202899, 0.2966485507246377} },
    b = {3.49139367, 3.09161577},
    k = {0.65194504, 0.29880919},
    thetaTrue = 1.2,
    phi[value_] := Tanh[value],
    measurementNoise[t_, amplitude_] :=
      amplitude (Sin[37. t] + 0.5 Sin[91. t]),
    saturate[value_, limit_] := Clip[value, {-limit, limit}],
    observerMatrix = a - Outer[Times, l, c],
    controllerMatrix = a - Outer[Times, b, k],
    Print["Observer eigenvalues: ", Eigenvalues[observerMatrix]],
    Print["Controller eigenvalues: ", Eigenvalues[controllerMatrix]],
    Print["Matching residual norm: ", Norm[p.b - c]],
    Print["Lyapunov residual norm: ",
      Norm[Transpose[observerMatrix].p + p.observerMatrix + IdentityMatrix[2]]],
    simulateCase[name_, gamma_, sigma_, noiseAmplitude_, controlLimit_] :=
      Module[{equations, initialConditions, solution, timeGrid, data},
       equations = {
         x1'[t] == x2[t] + b[[1]] (u[t] + thetaTrue phi[x1[t]]),
         x2'[t] == -2. x1[t] - 0.8 x2[t]
           + b[[2]] (u[t] + thetaTrue phi[x1[t]]),
         xh1'[t] == xh2[t] + b[[1]] (u[t] + th[t] phi[ym[t]])
           + l[[1]] (ym[t] - xh1[t]),
         xh2'[t] == -2. xh1[t] - 0.8 xh2[t]
           + b[[2]] (u[t] + th[t] phi[ym[t]])
           + l[[2]] (ym[t] - xh1[t]),
         th'[t] == gamma phi[ym[t]] (ym[t] - xh1[t]) - sigma gamma th[t],
         ym[t] == x1[t] + measurementNoise[t, noiseAmplitude],
         u[t] == saturate[-k.{xh1[t], xh2[t]} - th[t] phi[ym[t]], controlLimit]
       };
       initialConditions = {
         x1[0] == 1., x2[0] == -0.5,
         xh1[0] == 0., xh2[0] == 0., th[0] == 0.
       };
       solution = NDSolveValue[
         Join[equations, initialConditions],
         {x1, x2, xh1, xh2, th, u}, {t, 0., 20.},
         Method -> {"EquationSimplification" -> "Residual"},
         MaxStepFraction -> 1/5000
       ];
       timeGrid = Subdivide[0., 20., 4000];
       data = Table[
         {time, solution[[1]][time], solution[[2]][time],
          solution[[3]][time], solution[[4]][time],
          solution[[5]][time], solution[[6]][time]},
         {time, timeGrid}
       ];
       <|"Name" -> name, "Solution" -> solution, "Data" -> data|>
      ],
    ideal = simulateCase["Ideal matched case", 8., 0., 0., 10.^9],
    robust = simulateCase["Noisy robust case", 8., 0.03, 0.01, 8.],
    Print["Ideal final state and parameter: ", ideal["Data"][[-1, {2, 3, 6}]]],
    Print["Robust final state and parameter: ", robust["Data"][[-1, {2, 3, 6}]]],
    outputPlot = ListLinePlot[
      {ideal["Data"][[All, {1, 2}]], ideal["Data"][[All, {1, 4}]],
       robust["Data"][[All, {1, 2}]], robust["Data"][[All, {1, 4}]]},
      PlotLegends -> {"ideal x1", "ideal xhat1", "robust x1", "robust xhat1"},
      Frame -> True, FrameLabel -> {"time [s]", "output and estimate"},
      PlotRange -> All, ImageSize -> Large
    ],
    parameterPlot = ListLinePlot[
      {ideal["Data"][[All, {1, 6}]], robust["Data"][[All, {1, 6}]],
       Table[{time, thetaTrue}, {time, Subdivide[0., 20., 200]}]},
      PlotLegends -> {"ideal thetaHat", "robust thetaHat", "true theta"},
      Frame -> True, FrameLabel -> {"time [s]", "parameter estimate"},
      PlotRange -> All, ImageSize -> Large
    ],
    controlPlot = ListLinePlot[
      {ideal["Data"][[All, {1, 7}]], robust["Data"][[All, {1, 7}]]},
      PlotLegends -> {"ideal control", "robust control"},
      Frame -> True, FrameLabel -> {"time [s]", "control"},
      PlotRange -> All, ImageSize -> Large
    ],
    Column[{outputPlot, parameterPlot, controlPlot}]
   ]], StandardForm], "Input"]
}, WindowTitle -> "Chapter15_Lesson5", Saveable -> True]

18.1 Ideal matched case

The output innovation decays, the observer state approaches the plant state, and both states converge to the origin. The parameter estimate typically settles at a constant different from the true parameter. Once the state and regressor become small, the data no longer distinguish the remaining parameter error.

18.2 Noisy robust case

Leakage prevents unbounded drift, but it biases the estimate toward its nominal center. The state approaches a residual neighborhood rather than the exact origin because the noisy innovation continuously drives the observer and adaptive law. This is consistent with uniform ultimate boundedness, not with exact asymptotic convergence.

18.3 Experiments students should perform

1. increase \(\Gamma\) and record noise amplification;
2. set \(\sigma=0\) in the noisy case and observe parameter wandering;
3. reduce the control limit until saturation destroys the nominal transient;
4. replace \(b=P^{-1}c\) with a mismatched vector and compare the Lyapunov residual;
5. apply a persistently exciting command and compare parameter convergence with regulation.

Problem 1 (derive the observer error): For the plant and observer in Section 2, derive \(\dot e_o\) and identify every cancellation.

Solution: Since \(e_o=x-\hat{x}\),

\[ \begin{aligned} \dot e_o &=Ax+b[u+\theta^{T}\phi(y)]\\ &\quad-A\hat{x}-b[u+\hat{\theta}^{T}\phi(y)]-L(y-c^{T}\hat{x})\\ &=(A-Lc^{T})(x-\hat{x}) +b(\theta-\hat{\theta})^{T}\phi(y)\\ &=(A-Lc^{T})e_o+b\widetilde{\theta}^{T}\phi(y). \end{aligned} \]

The known control input cancels because the same applied input is used in the plant and observer. The estimated uncertainty cancels only partially, leaving the parameter-error term.

Problem 2 (measurable matching): Show that \(Pb=c\rho\) makes the Lyapunov cross-term measurable.

Solution:

\[ e_o^{T}Pb=e_o^{T}c\rho=\rho c^{T}e_o=\rho(y-\hat y)=\rho e_y. \]

Therefore the adaptive law can use only measured output and observer output while reproducing the cross-term required by the proof.

Problem 3 (sigma-modification bound): Starting from the leakage law in Section 9, derive an ultimate-bound inequality.

Solution: Exact matching cancels the nominal cross-term. Since \(\hat\theta-\theta_0=(\theta-\theta_0)-\widetilde\theta\),

\[ \begin{aligned} \dot V &=-e_o^{T}Qe_o +2\sigma\widetilde\theta^{T}(\hat\theta-\theta_0)\\ &=-e_o^{T}Qe_o-2\sigma\|\widetilde\theta\|^2 +2\sigma\widetilde\theta^{T}(\theta-\theta_0)\\ &\le-\lambda_{\min}(Q)\|e_o\|^2 -\sigma\|\widetilde\theta\|^2 +\sigma\|\theta-\theta_0\|^2. \end{aligned} \]

Outside a ball determined by the distance between the true parameter and the leakage center, \(\dot V<0\). Hence the errors are uniformly ultimately bounded.

Problem 4 (saturation mismatch): The observer uses \(u_c\), but the plant receives \(u=\operatorname{sat}(u_c)\). Derive the new error equation and explain the risk.

Solution: Let \(\Delta_u=u-u_c\). Then

\[ \dot e_o=A_o e_o+b\widetilde\theta^{T}\phi(y)+b\Delta_u. \]

The Lyapunov derivative contains \(2e_o^{T}Pb\Delta_u\), which has no fixed sign. Severe or persistent saturation can invalidate the asymptotic proof. Feeding the actual actuator signal to the observer removes this particular mismatch.

Problem 5 (state convergence without parameter convergence): Explain mathematically how \(x\to0\) can occur while \(\widetilde\theta\not\to0\).

Solution: For \(\phi(0)=0\), regulation drives the regressor to zero. The uncertain term then satisfies

\[ b\widetilde\theta^{T}\phi(y)\to0 \]

for any bounded constant parameter error. The closed loop no longer contains enough information to distinguish different parameter values. Persistent excitation is needed to infer parameter convergence from the vanishing prediction or tracking error.

Problem 6 (numerical design verification): For the matrices in Section 12, verify that \(Pb\approx c\) and that both nominal matrices are Hurwitz.

Solution: Direct multiplication gives

\[ Pb-c\approx \begin{bmatrix} 4.9\times10^{-10}\\0 \end{bmatrix}, \]

up to rounded constants. The observer eigenvalues are approximately \(-2.4\pm1.8547j\), and the controller eigenvalues are approximately \(-2.5\) and \(-1.5\). All have negative real parts.

Output-feedback adaptive stability is not obtained by independently combining a stable observer, a stable state-feedback controller, and an adaptive law. The augmented nonlinear interconnection must be analyzed. A measurable Lyapunov matching condition can recover a clean proof for a special structured class, but mismatch, high-gain peaking, insufficient excitation, measurement noise, unmodeled dynamics, and saturation create additional terms that may invalidate asymptotic stability. Projection, leakage, dead zones, normalization, saturation-aware observers, and explicit ultimate-bound analysis are therefore central design tools.