Consider the known linear time-invariant plant \( \dot{x}=Ax+Bu,\; y=Cx \). Let the controller and observer be

\[ u=-K\hat{x}, \qquad \dot{\hat{x} }=A\hat{x}+Bu+L(y-C\hat{x}). \]

Define the observation error \( e_o=x-\hat{x} \). Direct subtraction gives

\[ \dot e_o=(A-LC)e_o. \]

Since \( \hat{x}=x-e_o \), the plant state obeys

\[ \dot x=(A-BK)x+BK e_o. \]

The combined dynamics are therefore block upper triangular:

\[ \begin{bmatrix}\dot x\\ \dot e_o\end{bmatrix} = \begin{bmatrix} A-BK & BK\\ 0 & A-LC \end{bmatrix} \begin{bmatrix}x\\e_o\end{bmatrix}. \]

The characteristic polynomial factors exactly:

\[ \det\!\left(sI-\mathcal A_{cl}\right) =\det\!\left(sI-(A-BK)\right) \det\!\left(sI-(A-LC)\right). \]

Thus the controller poles and observer poles can be selected independently. This is the deterministic linear separation principle. Notice that the proof depends on a structural fact—not merely on both subsystems being stable: the lower-left block is exactly zero.

flowchart TD
  R["Reference"] --> C1["Fixed controller K"]
  Y["Measured output y"] --> O1["Fixed observer L"]
  O1 --> XH["State estimate xhat"]
  XH --> C1
  C1 --> U["Control u"]
  U --> P["Known linear plant"]
  P --> Y
  P --> X["Plant state x"]
  X --> EO["Observer error dynamics"]
  XH --> EO
  EO --> NOTE["Triangular combined matrix: independent pole design"]
        

An adaptive loop contains not only the plant and observer, but also a parameter estimator. For a matched, linearly parameterized uncertainty, consider

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

where \( \theta^*\in\mathbb R^p \) is an unknown constant vector and \( \phi(x)\in\mathbb R^p \) is a known regressor map. A certainty-equivalent observer and controller may be written as

\[ \begin{aligned} \dot{\hat x}&=A\hat x+B\left(u+\hat\theta^T\phi(\hat x)\right) +L(y-C\hat x),\\ u&=-K\hat x-\hat\theta^T\phi(\hat x),\\ \dot{\hat\theta}&=\Gamma\Psi(e_y,\hat x,u), \qquad e_y=y-C\hat x. \end{aligned} \]

Define \( e_o=x-\hat x \) and \( \tilde\theta=\theta^*-\hat\theta \). The complete adaptive state is \( z=[x^T,e_o^T,\tilde\theta^T]^T \). Even when the original plant is linear in its state, the augmented dynamics are nonlinear because products such as \( \tilde\theta^T\phi(x) \) and \( \hat\theta^T[\phi(x)-\phi(\hat x)] \) occur.

Substituting the adaptive control law into the plant gives

\[ \begin{aligned} \dot x &=(A-BK)x+BK e_o\\ &\quad+B\left[\tilde\theta^T\phi(x) +\hat\theta^T\left(\phi(x)-\phi(\hat x)\right)\right]. \end{aligned} \]

Subtracting the observer equation from the plant equation gives

\[ \begin{aligned} \dot e_o &=(A-LC)e_o\\ &\quad+B\left[\tilde\theta^T\phi(x) +\hat\theta^T\left(\phi(x)-\phi(\hat x)\right)\right]. \end{aligned} \]

The parameter error satisfies

\[ \dot{\tilde\theta}=-\Gamma\Psi(Ce_o,\hat x,u). \]

These equations reveal three bidirectional paths. First, the observer error enters the controller through \( \hat x=x-e_o \). Second, the parameter error drives both plant and observer dynamics. Third, the adaptation law is driven by signals generated by the plant, observer, and controller. Consequently, the combined Jacobian generally has the pattern

\[ J(z)= \begin{bmatrix} J_{xx} & J_{xe} & J_{x\theta}\\ J_{ex} & J_{ee} & J_{e\theta}\\ J_{\theta x} & J_{\theta e} & J_{\theta\theta} \end{bmatrix}, \]

with several nonzero off-diagonal blocks. It is not block triangular, so separate pole placement for \( K \) and \( L \) does not establish stability of the adaptive interconnection.

flowchart TD
  P["Plant state x"] --> Y["Measured output y"]
  Y --> O["Adaptive observer"]
  O --> XH["State estimate xhat"]
  XH --> C["Certainty-equivalent controller"]
  C --> U["Control u"]
  U --> P
  Y --> A["Parameter update law"]
  XH --> A
  U --> A
  A --> TH["Parameter estimate thetaHat"]
  TH --> O
  TH --> C
  P --> REG["Regressor mismatch \nphi(x) - phi(xhat)"]
  XH --> REG
  REG --> O
  REG --> P
  NOTE["Closed loop must be analyzed \nas one nonlinear system"] --> A
        

Suppose \( A-BK \) and \( A-LC \) are Hurwitz. This only proves stability of two frozen linear subsystems. It does not bound the nonlinear mismatch

\[ d(x,e_o,\tilde\theta) =B\left[\tilde\theta^T\phi(x) +\hat\theta^T\left(\phi(x)-\phi(x-e_o)\right)\right]. \]

If \( \phi \) is locally Lipschitz on a region \( \mathcal D \), then

\[ \|\phi(x)-\phi(\hat x)\| \leq \ell_\phi\|e_o\|, \qquad x,\hat x\in\mathcal D. \]

Projection can ensure \( \|\hat\theta\|\leq\bar\theta \), but the resulting term is still proportional to the observer error:

\[ \|B\hat\theta^T(\phi(x)-\phi(\hat x))\| \leq \|B\|\bar\theta\ell_\phi\|e_o\|. \]

Therefore, increasing observer speed, controller bandwidth, or adaptation gain changes the gain of the whole nonlinear feedback cycle. The familiar engineering instruction “make the observer several times faster than the controller” is only a heuristic here; it is not a proof.

Let positive-definite matrices \( P_c \) and \( P_o \) satisfy

\[ (A-BK)^TP_c+P_c(A-BK)=-Q_c, \qquad (A-LC)^TP_o+P_o(A-LC)=-Q_o, \]

with \( Q_c>0 \) and \( Q_o>0 \). A natural composite candidate is

\[ V=x^TP_cx+\rho e_o^TP_oe_o +\tilde\theta^T\Gamma^{-1}\tilde\theta, \qquad \rho>0. \]

After applying the adaptation law and bounding terms that do not cancel, one often obtains an inequality of the form

\[ \dot V \leq-\alpha_x\|x\|^2-\alpha_e\|e_o\|^2 +2\beta\|x\|\|e_o\|+\Delta(z), \]

where \( \Delta(z) \) contains residual parameter, disturbance, normalization, or leakage terms. Ignoring \( \Delta \) momentarily, define

\[ M=\begin{bmatrix}\alpha_x & -\beta\\-\beta & \alpha_e\end{bmatrix}. \]

The quadratic part is negative definite precisely when

\[ \alpha_x>0,\qquad \alpha_e>0, \qquad \alpha_x\alpha_e>\beta^2. \]

This is a small-gain-type condition. The quantities \( \alpha_x \), \( \alpha_e \), and \( \beta \) depend jointly on controller gains, observer gains, adaptation gains, regressor bounds, and parameter bounds. Hence the design is coupled even when separate Lyapunov equations are used.

If the residual satisfies \( \Delta(z)\leq d_0\|[x^T,e_o^T]^T\| \), then

\[ \dot V\leq-\lambda_M\|[x^T,e_o^T]^T\|^2 +d_0\|[x^T,e_o^T]^T\|, \]

where \( \lambda_M=\lambda_{\min}(M) \). This proves uniform ultimate boundedness outside the ball whose radius is proportional to \( d_0/\lambda_M \), not exact asymptotic separation.

Consider the uncertain oscillator-like plant

\[ \dot x_1=x_2,\qquad \dot x_2=-\theta^*x_1+u,\qquad y=x_1. \]

Use the observer, controller, and output-error update

\[ \begin{aligned} \dot{\hat x}_1&=\hat x_2+l_1(y-\hat x_1),\\ \dot{\hat x}_2&=-\hat\theta\hat x_1+u+l_2(y-\hat x_1),\\ u&=-k_1\hat x_1-k_2\hat x_2+\hat\theta\hat x_1,\\ \dot{\hat\theta}&=\gamma\hat x_1(y-\hat x_1). \end{aligned} \]

With \( e_1=x_1-\hat x_1 \), \( e_2=x_2-\hat x_2 \), and \( \tilde\theta=\theta^*-\hat\theta \), the observer errors satisfy

\[ \dot e_1=e_2-l_1e_1, \qquad \dot e_2=-(\theta^*+l_2)e_1-\tilde\theta\hat x_1. \]

Choose the candidate

\[ V_o=\frac{1}{2}q_1e_1^2+\frac{1}{2}q_2e_2^2 +\frac{1}{2\gamma}\tilde\theta^2. \]

Its derivative is

\[ \begin{aligned} \dot V_o &=-q_1l_1e_1^2 +\left(q_1-q_2(\theta^*+l_2)\right)e_1e_2\\ &\quad-q_2e_2\tilde\theta\hat x_1 -e_1\tilde\theta\hat x_1. \end{aligned} \]

Selecting \( q_1=q_2(\theta^*+l_2) \) removes the \( e_1e_2 \) term, but two indefinite parameter-error products remain. The update based only on the measurable output error \( e_1 \) cannot cancel the term containing the unmeasured error \( e_2 \). This is the essential output-feedback obstruction: a state-feedback Lyapunov cancellation may require a signal that is unavailable from the measured output.

Adaptive output-feedback designs therefore introduce additional filtered errors, augmented errors, dynamic regressor extensions, high-gain mechanisms, or specially structured adaptive observers. Those techniques will be developed conceptually in the next lesson; here the important conclusion is that naïve substitution of estimated states is not a separation theorem.

A defensible workflow treats the observer, controller, and estimator as one nonlinear system. The following checks are especially important:

1. Derive the exact augmented error model before tuning gains. Do not infer the adaptive dynamics from the frozen linear poles.
2. Identify which Lyapunov cross terms can be cancelled by measurable signals and which remain as perturbations.
3. Establish bounded regressors and parameter estimates. Projection and normalization from Chapter 9 are often structural requirements, not cosmetic numerical options.
4. Increase \( \gamma \) only after checking the coupled gain inequality. Faster adaptation can amplify observer transients and noise instead of improving tracking.
5. Use leakage, dead zones, or projection to obtain boundedness under imperfect cancellation. These modifications do not recreate exact linear separation; they reshape the residual term \( \Delta(z) \).
6. Test large initial observer and parameter errors. Local simulations near the equilibrium can hide peaking and finite-escape behavior.

In singular-perturbation terminology, a fast observer may sometimes justify an approximate two-time-scale argument. Such an argument requires explicit bounds, a stable boundary-layer system, and a reduced adaptive system that remains stable. Merely placing observer poles far left is not sufficient.

The implementations below simulate the second-order example using a from-scratch fourth-order Runge–Kutta method. Three scenarios use the same stabilizing controller and observer gains:

\[ k_1=4,\quad k_2=3,\quad l_1=4,\quad l_2=4, \quad \theta^*=2. \]

The moderate case uses \( \gamma=10 \). The aggressive case uses \( \gamma=100 \) and demonstrates coupled instability. The robustified case retains \( \gamma=100 \) but adds \( \sigma=1 \) leakage and projection to \( 0\leq\hat\theta\leq4 \). For the supplied initial condition, the aggressive run reaches the numerical blow-up threshold at approximately 2.11 seconds, whereas the other two runs remain bounded.

The robustified estimate does not converge to the true parameter because regulation removes excitation and leakage biases the estimate toward its nominal value. This is consistent with Chapter 10: bounded tracking does not imply parameter convergence.

This version uses NumPy for arrays and Matplotlib for plots. The RK4 integrator, projection logic, stopping test, and metrics are implemented explicitly.

Chapter15_Lesson3.py

"""Chapter 15, Lesson 3: nonlinear coupling in adaptive output feedback.

The script uses a from-scratch RK4 solver to compare three cases:
1) moderate adaptation,
2) aggressive adaptation that may trigger coupled instability,
3) the same aggressive gain with sigma modification and projection.

Dependencies: numpy, matplotlib.
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    name: str
    gamma: float
    sigma: float
    projection: bool
    theta_min: float = 0.0
    theta_max: float = 4.0


@dataclass
class Result:
    time: np.ndarray
    state: np.ndarray
    stopped_early: bool
    stop_reason: str


THETA_TRUE = 2.0
THETA_NOMINAL = 0.5
K1, K2 = 4.0, 3.0
L1, L2 = 4.0, 4.0
DT = 1.0e-4
T_FINAL = 10.0
BLOWUP_LIMIT = 1.0e6
INITIAL_STATE = np.array([1.0, 0.0, 0.0, 0.0, THETA_NOMINAL], dtype=float)


def derivative(z: np.ndarray, cfg: Config) -> np.ndarray:
    """Return [x1_dot, x2_dot, xhat1_dot, xhat2_dot, theta_hat_dot]."""
    x1, x2, xhat1, xhat2, theta_hat = z
    output_error = x1 - xhat1

    # Certainty-equivalent controller.
    control = -K1 * xhat1 - K2 * xhat2 + theta_hat * xhat1

    theta_hat_dot = (
        cfg.gamma * xhat1 * output_error
        - cfg.sigma * (theta_hat - THETA_NOMINAL)
    )

    # Tangent-cone projection: block only outward motion at a boundary.
    if cfg.projection:
        if theta_hat <= cfg.theta_min and theta_hat_dot < 0.0:
            theta_hat_dot = 0.0
        if theta_hat >= cfg.theta_max and theta_hat_dot > 0.0:
            theta_hat_dot = 0.0

    return np.array(
        [
            x2,
            -THETA_TRUE * x1 + control,
            xhat2 + L1 * output_error,
            -theta_hat * xhat1 + control + L2 * output_error,
            theta_hat_dot,
        ],
        dtype=float,
    )


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


def simulate(cfg: Config) -> Result:
    steps = int(round(T_FINAL / DT))
    time = np.empty(steps + 1, dtype=float)
    state = np.empty((steps + 1, 5), dtype=float)
    time[0] = 0.0
    state[0] = INITIAL_STATE

    stopped_early = False
    stop_reason = "completed"
    last = steps

    for i in range(steps):
        z_next = rk4_step(state[i], DT, lambda z: derivative(z, cfg))

        if cfg.projection:
            z_next[4] = np.clip(z_next[4], cfg.theta_min, cfg.theta_max)

        time[i + 1] = time[i] + DT
        state[i + 1] = z_next

        if not np.all(np.isfinite(z_next)):
            stopped_early = True
            stop_reason = "non-finite state"
            last = i + 1
            break
        if np.max(np.abs(z_next)) > BLOWUP_LIMIT:
            stopped_early = True
            stop_reason = f"state magnitude exceeded {BLOWUP_LIMIT:.1e}"
            last = i + 1
            break

    return Result(time[: last + 1], state[: last + 1], stopped_early, stop_reason)


def summarize(cfg: Config, result: Result) -> None:
    y = result.state[:, 0]
    output_error = result.state[:, 0] - result.state[:, 2]
    rms_y = float(np.sqrt(np.mean(y * y)))
    rms_error = float(np.sqrt(np.mean(output_error * output_error)))
    max_norm = float(np.max(np.abs(result.state)))
    final_theta = float(result.state[-1, 4])

    print(f"\nScenario: {cfg.name}")
    print(f"  stopped early : {result.stopped_early} ({result.stop_reason})")
    print(f"  final time    : {result.time[-1]:.6f} s")
    print(f"  RMS output    : {rms_y:.6e}")
    print(f"  RMS obs. error: {rms_error:.6e}")
    print(f"  max |state|   : {max_norm:.6e}")
    print(f"  final thetaHat: {final_theta:.6f}")


def save_csv(cfg: Config, result: Result, output_dir: Path) -> None:
    output = np.column_stack((result.time, result.state))
    header = "time,x1,x2,xhat1,xhat2,theta_hat"
    np.savetxt(output_dir / f"Chapter15_Lesson3_{cfg.name}.csv", output, delimiter=",", header=header, comments="")


def make_plots(results: list[tuple[Config, Result]], output_dir: Path) -> None:
    plt.figure(figsize=(10, 6))
    for cfg, result in results:
        plt.plot(result.time, result.state[:, 0], label=cfg.name)
    plt.xlabel("Time (s)")
    plt.ylabel("Measured output y = x1")
    plt.title("Adaptive output-feedback coupling benchmark")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter15_Lesson3_output.png", dpi=180)

    plt.figure(figsize=(10, 6))
    for cfg, result in results:
        plt.plot(result.time, result.state[:, 4], label=cfg.name)
    plt.axhline(THETA_TRUE, linestyle="--", label="true theta")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.title("Parameter evolution under observer-controller coupling")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter15_Lesson3_parameter.png", dpi=180)
    plt.show()


def main() -> None:
    output_dir = Path(__file__).resolve().parent
    scenarios = [
        Config("moderate", gamma=10.0, sigma=0.0, projection=False),
        Config("aggressive", gamma=100.0, sigma=0.0, projection=False),
        Config("robustified", gamma=100.0, sigma=1.0, projection=True),
    ]

    results: list[tuple[Config, Result]] = []
    for cfg in scenarios:
        result = simulate(cfg)
        summarize(cfg, result)
        save_csv(cfg, result, output_dir)
        results.append((cfg, result))

    make_plots(results, output_dir)


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library. It writes one CSV file per scenario and stops safely when a state exceeds the configured blow-up limit.

Chapter15_Lesson3.cpp

// Chapter 15, Lesson 3: nonlinear coupling in adaptive output feedback.
// Build: g++ -std=c++17 -O2 Chapter15_Lesson3.cpp -o Chapter15_Lesson3

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

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

struct Config {
    std::string name;
    double gamma;
    double sigma;
    bool projection;
    double thetaMin{0.0};
    double thetaMax{4.0};
};

struct Sample {
    double time;
    State state;
};

struct Result {
    std::vector<Sample> samples;
    bool stoppedEarly{false};
    std::string stopReason{"completed"};
};

constexpr double THETA_TRUE = 2.0;
constexpr double THETA_NOMINAL = 0.5;
constexpr double K1 = 4.0;
constexpr double K2 = 3.0;
constexpr double L1 = 4.0;
constexpr double L2 = 4.0;
constexpr double DT = 1.0e-4;
constexpr double T_FINAL = 10.0;
constexpr double BLOWUP_LIMIT = 1.0e6;

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

State derivative(const State& z, const Config& cfg) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xhat1 = z[2];
    const double xhat2 = z[3];
    const double thetaHat = z[4];
    const double outputError = x1 - xhat1;
    const double control = -K1 * xhat1 - K2 * xhat2 + thetaHat * xhat1;

    double thetaHatDot = cfg.gamma * xhat1 * outputError
                       - cfg.sigma * (thetaHat - THETA_NOMINAL);

    if (cfg.projection) {
        if (thetaHat <= cfg.thetaMin && thetaHatDot < 0.0) thetaHatDot = 0.0;
        if (thetaHat >= cfg.thetaMax && thetaHatDot > 0.0) thetaHatDot = 0.0;
    }

    return State{
        x2,
        -THETA_TRUE * x1 + control,
        xhat2 + L1 * outputError,
        -thetaHat * xhat1 + control + L2 * outputError,
        thetaHatDot
    };
}

State rk4Step(const State& z, const Config& cfg) {
    const State k1 = derivative(z, cfg);
    const State k2 = derivative(addScaled(z, k1, 0.5 * DT), cfg);
    const State k3 = derivative(addScaled(z, k2, 0.5 * DT), cfg);
    const State k4 = derivative(addScaled(z, k3, DT), cfg);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    if (cfg.projection) {
        next[4] = std::clamp(next[4], cfg.thetaMin, cfg.thetaMax);
    }
    return next;
}

bool finiteState(const State& z) {
    for (double value : z) {
        if (!std::isfinite(value)) return false;
    }
    return true;
}

double maxAbs(const State& z) {
    double result = 0.0;
    for (double value : z) result = std::max(result, std::abs(value));
    return result;
}

Result simulate(const Config& cfg) {
    const std::size_t steps = static_cast<std::size_t>(std::llround(T_FINAL / DT));
    Result result;
    result.samples.reserve(steps + 1);

    State z{1.0, 0.0, 0.0, 0.0, THETA_NOMINAL};
    result.samples.push_back({0.0, z});

    for (std::size_t i = 0; i < steps; ++i) {
        z = rk4Step(z, cfg);
        const double t = (i + 1) * DT;
        result.samples.push_back({t, z});

        if (!finiteState(z)) {
            result.stoppedEarly = true;
            result.stopReason = "non-finite state";
            break;
        }
        if (maxAbs(z) > BLOWUP_LIMIT) {
            result.stoppedEarly = true;
            result.stopReason = "state magnitude exceeded blow-up limit";
            break;
        }
    }
    return result;
}

void saveCsv(const Config& cfg, const Result& result) {
    std::ofstream file("Chapter15_Lesson3_" + cfg.name + ".csv");
    if (!file) throw std::runtime_error("Cannot create CSV output.");

    file << "time,x1,x2,xhat1,xhat2,theta_hat\n";
    file << std::setprecision(12);
    for (const auto& sample : result.samples) {
        file << sample.time;
        for (double value : sample.state) file << ',' << value;
        file << '\n';
    }
}

void summarize(const Config& cfg, const Result& result) {
    double sumY2 = 0.0;
    double sumE2 = 0.0;
    double maximum = 0.0;

    for (const auto& sample : result.samples) {
        const double y = sample.state[0];
        const double error = sample.state[0] - sample.state[2];
        sumY2 += y * y;
        sumE2 += error * error;
        maximum = std::max(maximum, maxAbs(sample.state));
    }

    const double count = static_cast<double>(result.samples.size());
    const State& finalState = result.samples.back().state;
    std::cout << "\nScenario: " << cfg.name << '\n'
              << "  stopped early : " << std::boolalpha << result.stoppedEarly
              << " (" << result.stopReason << ")\n"
              << "  final time    : " << result.samples.back().time << " s\n"
              << "  RMS output    : " << std::sqrt(sumY2 / count) << '\n'
              << "  RMS obs. error: " << std::sqrt(sumE2 / count) << '\n'
              << "  max |state|   : " << maximum << '\n'
              << "  final thetaHat: " << finalState[4] << '\n';
}

int main() {
    try {
        const std::vector<Config> scenarios{
            {"moderate", 10.0, 0.0, false},
            {"aggressive", 100.0, 0.0, false},
            {"robustified", 100.0, 1.0, true}
        };

        for (const Config& cfg : scenarios) {
            const Result result = simulate(cfg);
            summarize(cfg, result);
            saveCsv(cfg, result);
        }
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses records for configuration and samples, standard arrays for the state, and java.nio.file for CSV output. It requires a modern JDK supporting records.

Chapter15_Lesson3.java

// Chapter 15, Lesson 3: nonlinear coupling in adaptive output feedback.
// Build and run: javac Chapter15_Lesson3.java && java Chapter15_Lesson3

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

public final class Chapter15_Lesson3 {
    private static final double THETA_TRUE = 2.0;
    private static final double THETA_NOMINAL = 0.5;
    private static final double K1 = 4.0;
    private static final double K2 = 3.0;
    private static final double L1 = 4.0;
    private static final double L2 = 4.0;
    private static final double DT = 1.0e-4;
    private static final double T_FINAL = 10.0;
    private static final double BLOWUP_LIMIT = 1.0e6;

    private record Config(
            String name,
            double gamma,
            double sigma,
            boolean projection,
            double thetaMin,
            double thetaMax) {
        Config(String name, double gamma, double sigma, boolean projection) {
            this(name, gamma, sigma, projection, 0.0, 4.0);
        }
    }

    private record Sample(double time, double[] state) {}

    private static final class Result {
        final List<Sample> samples = new ArrayList<>();
        boolean stoppedEarly = false;
        String stopReason = "completed";
    }

    private static double[] derivative(double[] z, Config cfg) {
        double x1 = z[0];
        double x2 = z[1];
        double xhat1 = z[2];
        double xhat2 = z[3];
        double thetaHat = z[4];
        double outputError = x1 - xhat1;
        double control = -K1 * xhat1 - K2 * xhat2 + thetaHat * xhat1;

        double thetaHatDot = cfg.gamma() * xhat1 * outputError
                - cfg.sigma() * (thetaHat - THETA_NOMINAL);

        if (cfg.projection()) {
            if (thetaHat <= cfg.thetaMin() && thetaHatDot < 0.0) thetaHatDot = 0.0;
            if (thetaHat >= cfg.thetaMax() && thetaHatDot > 0.0) thetaHatDot = 0.0;
        }

        return new double[] {
                x2,
                -THETA_TRUE * x1 + control,
                xhat2 + L1 * outputError,
                -thetaHat * xhat1 + control + L2 * outputError,
                thetaHatDot
        };
    }

    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[] z, Config cfg) {
        double[] k1 = derivative(z, cfg);
        double[] k2 = derivative(addScaled(z, k1, 0.5 * DT), cfg);
        double[] k3 = derivative(addScaled(z, k2, 0.5 * DT), cfg);
        double[] k4 = derivative(addScaled(z, k3, DT), cfg);

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        if (cfg.projection()) {
            next[4] = Math.max(cfg.thetaMin(), Math.min(cfg.thetaMax(), next[4]));
        }
        return next;
    }

    private static boolean finiteState(double[] z) {
        for (double value : z) if (!Double.isFinite(value)) return false;
        return true;
    }

    private static double maxAbs(double[] z) {
        double maximum = 0.0;
        for (double value : z) maximum = Math.max(maximum, Math.abs(value));
        return maximum;
    }

    private static Result simulate(Config cfg) {
        int steps = (int) Math.round(T_FINAL / DT);
        Result result = new Result();
        double[] z = {1.0, 0.0, 0.0, 0.0, THETA_NOMINAL};
        result.samples.add(new Sample(0.0, z.clone()));

        for (int i = 0; i < steps; i++) {
            z = rk4Step(z, cfg);
            result.samples.add(new Sample((i + 1) * DT, z.clone()));

            if (!finiteState(z)) {
                result.stoppedEarly = true;
                result.stopReason = "non-finite state";
                break;
            }
            if (maxAbs(z) > BLOWUP_LIMIT) {
                result.stoppedEarly = true;
                result.stopReason = "state magnitude exceeded blow-up limit";
                break;
            }
        }
        return result;
    }

    private static void summarize(Config cfg, Result result) {
        double sumY2 = 0.0;
        double sumE2 = 0.0;
        double maximum = 0.0;

        for (Sample sample : result.samples) {
            double[] z = sample.state();
            double y = z[0];
            double error = z[0] - z[2];
            sumY2 += y * y;
            sumE2 += error * error;
            maximum = Math.max(maximum, maxAbs(z));
        }

        double count = result.samples.size();
        Sample last = result.samples.get(result.samples.size() - 1);
        System.out.printf(Locale.US,
                "%nScenario: %s%n  stopped early : %s (%s)%n  final time    : %.6f s%n"
                        + "  RMS output    : %.6e%n  RMS obs. error: %.6e%n"
                        + "  max |state|   : %.6e%n  final thetaHat: %.6f%n",
                cfg.name(), result.stoppedEarly, result.stopReason, last.time(),
                Math.sqrt(sumY2 / count), Math.sqrt(sumE2 / count), maximum, last.state()[4]);
    }

    private static void saveCsv(Config cfg, Result result) throws IOException {
        Path path = Path.of("Chapter15_Lesson3_" + cfg.name() + ".csv");
        try (BufferedWriter writer = Files.newBufferedWriter(path)) {
            writer.write("time,x1,x2,xhat1,xhat2,theta_hat\n");
            for (Sample sample : result.samples) {
                double[] z = sample.state();
                writer.write(String.format(Locale.US,
                        "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                        sample.time(), z[0], z[1], z[2], z[3], z[4]));
            }
        }
    }

    public static void main(String[] args) {
        List<Config> scenarios = List.of(
                new Config("moderate", 10.0, 0.0, false),
                new Config("aggressive", 100.0, 0.0, false),
                new Config("robustified", 100.0, 1.0, true));

        try {
            for (Config cfg : scenarios) {
                Result result = simulate(cfg);
                summarize(cfg, result);
                saveCsv(cfg, result);
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script requires only base MATLAB. Comments at the end map the five-state differential equations to Integrator, Sum, Product, Saturation, and Gain blocks for an equivalent Simulink model. The key modeling point is that all five states must remain in one nonlinear loop.

Chapter15_Lesson3.m

% Chapter 15, Lesson 3: nonlinear coupling in adaptive output feedback.
% This script compares moderate, aggressive, and robustified adaptation.
% It uses a local, from-scratch RK4 implementation and requires base MATLAB.

clear; clc; close all;

scenarios = struct( ...
    'name',       {'moderate', 'aggressive', 'robustified'}, ...
    'gamma',      {10.0,       100.0,        100.0}, ...
    'sigma',      {0.0,        0.0,          1.0}, ...
    'projection', {false,      false,        true}, ...
    'thetaMin',   {0.0,        0.0,          0.0}, ...
    'thetaMax',   {4.0,        4.0,          4.0});

results = cell(size(scenarios));
for i = 1:numel(scenarios)
    results{i} = simulateScenario(scenarios(i));
    printSummary(scenarios(i), results{i});
    writeResult(scenarios(i), results{i});
end

figure;
hold on;
for i = 1:numel(scenarios)
    plot(results{i}.time, results{i}.state(:, 1), 'DisplayName', scenarios(i).name);
end
xlabel('Time (s)'); ylabel('Measured output y = x_1');
title('Adaptive output-feedback coupling benchmark');
grid on; legend('Location', 'best');

figure;
hold on;
for i = 1:numel(scenarios)
    plot(results{i}.time, results{i}.state(:, 5), 'DisplayName', scenarios(i).name);
end
yline(2.0, '--', 'true theta');
xlabel('Time (s)'); ylabel('Parameter estimate');
title('Parameter evolution under observer-controller coupling');
grid on; legend('Location', 'best');

% Simulink mapping:
% 1. Use two Integrator blocks for x_1 and x_2.
% 2. Use two Integrator blocks for xhat_1 and xhat_2.
% 3. Form e_y = x_1 - xhat_1 with a Sum block.
% 4. Implement u = -k_1*xhat_1 - k_2*xhat_2 + thetaHat*xhat_1.
% 5. Integrate thetaHatDot; add Saturation [0,4] for the projected case.
% 6. Add a leakage branch -sigma*(thetaHat-thetaNominal).
% The complete five-state loop must be simulated as one nonlinear model.

function result = simulateScenario(cfg)
    dt = 1.0e-4;
    tFinal = 10.0;
    blowupLimit = 1.0e6;
    steps = round(tFinal / dt);

    time = zeros(steps + 1, 1);
    state = zeros(steps + 1, 5);
    state(1, :) = [1.0, 0.0, 0.0, 0.0, 0.5];
    stoppedEarly = false;
    stopReason = 'completed';
    last = steps + 1;

    for k = 1:steps
        z = state(k, :).';
        k1 = rhs(z, cfg);
        k2 = rhs(z + 0.5 * dt * k1, cfg);
        k3 = rhs(z + 0.5 * dt * k2, cfg);
        k4 = rhs(z + dt * k3, cfg);
        zNext = z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);

        if cfg.projection
            zNext(5) = min(max(zNext(5), cfg.thetaMin), cfg.thetaMax);
        end

        time(k + 1) = time(k) + dt;
        state(k + 1, :) = zNext.';

        if any(~isfinite(zNext))
            stoppedEarly = true;
            stopReason = 'non-finite state';
            last = k + 1;
            break;
        end
        if max(abs(zNext)) > blowupLimit
            stoppedEarly = true;
            stopReason = 'state magnitude exceeded blow-up limit';
            last = k + 1;
            break;
        end
    end

    result.time = time(1:last);
    result.state = state(1:last, :);
    result.stoppedEarly = stoppedEarly;
    result.stopReason = stopReason;
end

function dz = rhs(z, cfg)
    thetaTrue = 2.0;
    thetaNominal = 0.5;
    k1 = 4.0; k2 = 3.0;
    l1 = 4.0; l2 = 4.0;

    x1 = z(1); x2 = z(2);
    xhat1 = z(3); xhat2 = z(4);
    thetaHat = z(5);
    outputError = x1 - xhat1;
    control = -k1 * xhat1 - k2 * xhat2 + thetaHat * xhat1;

    thetaHatDot = cfg.gamma * xhat1 * outputError ...
        - cfg.sigma * (thetaHat - thetaNominal);

    if cfg.projection
        if thetaHat <= cfg.thetaMin && thetaHatDot < 0.0
            thetaHatDot = 0.0;
        end
        if thetaHat >= cfg.thetaMax && thetaHatDot > 0.0
            thetaHatDot = 0.0;
        end
    end

    dz = [ ...
        x2; ...
        -thetaTrue * x1 + control; ...
        xhat2 + l1 * outputError; ...
        -thetaHat * xhat1 + control + l2 * outputError; ...
        thetaHatDot];
end

function printSummary(cfg, result)
    y = result.state(:, 1);
    outputError = result.state(:, 1) - result.state(:, 3);
    rmsY = sqrt(mean(y .^ 2));
    rmsError = sqrt(mean(outputError .^ 2));
    maxState = max(abs(result.state), [], 'all');
    finalTheta = result.state(end, 5);

    fprintf('\nScenario: %s\n', cfg.name);
    fprintf('  stopped early : %d (%s)\n', result.stoppedEarly, result.stopReason);
    fprintf('  final time    : %.6f s\n', result.time(end));
    fprintf('  RMS output    : %.6e\n', rmsY);
    fprintf('  RMS obs. error: %.6e\n', rmsError);
    fprintf('  max |state|   : %.6e\n', maxState);
    fprintf('  final thetaHat: %.6f\n', finalTheta);
end

function writeResult(cfg, result)
    data = array2table([result.time, result.state], ...
        'VariableNames', {'time', 'x1', 'x2', 'xhat1', 'xhat2', 'theta_hat'});
    writetable(data, ['Chapter15_Lesson3_', cfg.name, '.csv']);
end

The Wolfram Language version implements the same RK4 recursion, exports CSV data, prints metrics, and creates output and parameter plots.

Chapter15_Lesson3.nb

(* Chapter 15, Lesson 3: nonlinear coupling in adaptive output feedback. *)
(* Save this text as Chapter15_Lesson3.nb and evaluate it in Wolfram Mathematica. *)

ClearAll["Global`*"];

thetaTrue = 2.0;
thetaNominal = 0.5;
k1 = 4.0; k2 = 3.0;
l1 = 4.0; l2 = 4.0;
dt = 1.0*^-4;
tFinal = 10.0;
blowupLimit = 1.0*^6;
initialState = {1.0, 0.0, 0.0, 0.0, thetaNominal};

scenarios = {
   <|"Name" -> "moderate", "Gamma" -> 10.0, "Sigma" -> 0.0,
     "Projection" -> False, "ThetaMin" -> 0.0, "ThetaMax" -> 4.0|>,
   <|"Name" -> "aggressive", "Gamma" -> 100.0, "Sigma" -> 0.0,
     "Projection" -> False, "ThetaMin" -> 0.0, "ThetaMax" -> 4.0|>,
   <|"Name" -> "robustified", "Gamma" -> 100.0, "Sigma" -> 1.0,
     "Projection" -> True, "ThetaMin" -> 0.0, "ThetaMax" -> 4.0|>
   };

rhs[z_List, cfg_Association] := Module[
  {x1, x2, xhat1, xhat2, thetaHat, outputError, control, thetaHatDot},
  {x1, x2, xhat1, xhat2, thetaHat} = z;
  outputError = x1 - xhat1;
  control = -k1*xhat1 - k2*xhat2 + thetaHat*xhat1;
  thetaHatDot = cfg["Gamma"]*xhat1*outputError
    - cfg["Sigma"]*(thetaHat - thetaNominal);

  If[cfg["Projection"],
   If[thetaHat <= cfg["ThetaMin"] && thetaHatDot < 0.0, thetaHatDot = 0.0];
   If[thetaHat >= cfg["ThetaMax"] && thetaHatDot > 0.0, thetaHatDot = 0.0];
   ];

  {x2,
   -thetaTrue*x1 + control,
   xhat2 + l1*outputError,
   -thetaHat*xhat1 + control + l2*outputError,
   thetaHatDot}
  ];

rk4Step[z_List, cfg_Association] := Module[{a, b, c, d, next},
  a = rhs[z, cfg];
  b = rhs[z + 0.5*dt*a, cfg];
  c = rhs[z + 0.5*dt*b, cfg];
  d = rhs[z + dt*c, cfg];
  next = z + (dt/6.0)*(a + 2.0*b + 2.0*c + d);
  If[cfg["Projection"],
   next[[5]] = Clip[next[[5]], {cfg["ThetaMin"], cfg["ThetaMax"]}]
   ];
  next
  ];

simulate[cfg_Association] := Module[
  {steps, z, rows, stoppedEarly, stopReason, i, maximum},
  steps = Round[tFinal/dt];
  z = initialState;
  rows = { {0.0, Sequence @@ z} };
  stoppedEarly = False;
  stopReason = "completed";

  For[i = 1, i <= steps, i++,
   z = rk4Step[z, cfg];
   AppendTo[rows, {i*dt, Sequence @@ z}];
   maximum = Max[Abs[z]];
   If[! VectorQ[z, NumericQ],
    stoppedEarly = True; stopReason = "non-numeric state"; Break[]];
   If[maximum > blowupLimit,
    stoppedEarly = True;
    stopReason = "state magnitude exceeded blow-up limit";
    Break[]];
   ];

  <|"Data" -> rows, "StoppedEarly" -> stoppedEarly,
    "StopReason" -> stopReason|>
  ];

results = Association@Table[cfg["Name"] -> simulate[cfg], {cfg, scenarios}];

Do[
 data = results[cfg["Name"]]["Data"];
 y = data[[All, 2]];
 outputError = data[[All, 2]] - data[[All, 4]];
 Print["\nScenario: ", cfg["Name"]];
 Print["  stopped early : ", results[cfg["Name"]]["StoppedEarly"],
   " (", results[cfg["Name"]]["StopReason"], ")"];
 Print["  final time    : ", Last[data][[1]], " s"];
 Print["  RMS output    : ", Sqrt[Mean[y^2]]];
 Print["  RMS obs. error: ", Sqrt[Mean[outputError^2]]];
 Print["  max |state|   : ", Max[Abs[data[[All, 2 ;; 6]]]]];
 Print["  final thetaHat: ", Last[data][[6]]];
 Export["Chapter15_Lesson3_" <> cfg["Name"] <> ".csv",
  Prepend[data, {"time", "x1", "x2", "xhat1", "xhat2", "theta_hat"}]];
 , {cfg, scenarios}];

outputPlot = ListLinePlot[
   Table[results[cfg["Name"]]["Data"][[All, {1, 2}]], {cfg, scenarios}],
   PlotLegends -> scenarios[[All, "Name"]],
   AxesLabel -> {"Time (s)", "y = x1"},
   PlotLabel -> "Adaptive output-feedback coupling benchmark",
   PlotRange -> All, ImageSize -> Large];

parameterPlot = ListLinePlot[
   Table[results[cfg["Name"]]["Data"][[All, {1, 6}]], {cfg, scenarios}],
   PlotLegends -> scenarios[[All, "Name"]],
   AxesLabel -> {"Time (s)", "thetaHat"},
   PlotLabel -> "Parameter evolution under observer-controller coupling",
   Epilog -> {Dashed, Line[{ {0, thetaTrue}, {tFinal, thetaTrue} }]},
   PlotRange -> All, ImageSize -> Large];

Column[{outputPlot, parameterPlot}]

Problem 1 (Classical factorization): For the known LTI plant and observer in Section 1, prove that the spectrum of the combined system is the union of the spectra of \( A-BK \) and \( A-LC \).

Solution:

In coordinates \( [x^T,e_o^T]^T \), the closed-loop matrix is block upper triangular. For any block upper-triangular matrix, the determinant is the product of the determinants of its diagonal blocks. Hence

\[ \det(sI-\mathcal A_{cl}) =\det(sI-(A-BK))\det(sI-(A-LC)). \]

The roots are therefore the controller poles together with the observer poles. The result relies on the zero lower-left block.

Problem 2 (Adaptive mismatch): Starting from the plant, observer, and certainty-equivalent controller of Section 2, derive the common nonlinear mismatch term that appears in both \( \dot x \) and \( \dot e_o \).

Solution:

Add and subtract \( \hat\theta^T\phi(x) \):

\[ \theta^{*T}\phi(x)-\hat\theta^T\phi(\hat x) =\tilde\theta^T\phi(x) +\hat\theta^T(\phi(x)-\phi(\hat x)). \]

The first term is caused by parameter error; the second is caused by observer error inside the nonlinear regressor. Both enter the plant and observer-error equations through \( B \), producing bidirectional coupling.

Problem 3 (Failure of naïve cancellation): For the second-order example, compute \( \dot V_o \) and show why an update based only on \( e_1 \) cannot cancel all parameter-error terms.

Solution:

Differentiating the candidate from Section 6 and using \( \dot{\tilde\theta}=-\gamma\hat x_1e_1 \) yields

\[ \dot V_o=-q_1l_1e_1^2 +(q_1-q_2(\theta^*+l_2))e_1e_2 -q_2e_2\tilde\theta\hat x_1 -e_1\tilde\theta\hat x_1. \]

Choosing \( q_1=q_2(\theta^*+l_2) \) removes only the state cross term. Cancellation of the term containing \( e_2 \) would require access to an unmeasured observer-error component or an equivalent filtered signal. Therefore the output-error gradient law is not a complete Lyapunov design.

Problem 4 (Coupled quadratic condition): Determine when the quadratic form

\[ -\alpha_x\xi^2-\alpha_e\eta^2+2\beta\xi\eta \]

is negative definite.

Solution:

Write the negative of the form as \( [\xi,\eta]M[\xi,\eta]^T \), where

\[ M=\begin{bmatrix}\alpha_x&-\beta\\-\beta&\alpha_e\end{bmatrix}. \]

By Sylvester’s criterion, \( M>0 \) if and only if \( \alpha_x>0 \) and \( \det M=\alpha_x\alpha_e-\beta^2>0 \). Thus the required condition is \( \alpha_x\alpha_e>\beta^2 \) together with positive diagonal coefficients.

Problem 5 (Interpreting the benchmark): Explain why the robustified simulation may regulate the state while \( \hat\theta \) approaches the nominal value rather than \( \theta^* \).

Solution:

Once the state approaches zero, both \( \hat x_1 \) and the output error become small, so the informative gradient term loses excitation. The leakage term \( -\sigma(\hat\theta-\theta_0) \) then dominates and attracts the estimate toward \( \theta_0 \). Projection prevents escape but does not create persistent excitation. Consequently, bounded regulation and parameter convergence remain distinct objectives.

Classical linear separation follows from a block-triangular closed-loop matrix. Adaptive output feedback destroys that structure because state estimation errors, parameter errors, regressor mismatch, and adaptation dynamics feed one another. Stability must therefore be established for a composite nonlinear system. Lyapunov inequalities, small-gain bounds, measurable-error constructions, projection, normalization, and robust modifications are the appropriate tools. The numerical benchmark shows the practical consequence: increasing the adaptation gain can destabilize an observer and controller that are each stable when considered in isolation.