After completing this lesson, students should be able to:

• relate aerodynamic coefficient changes to variations in dimensional stability derivatives;
• separate a certified or otherwise trusted baseline controller from an adaptive augmentation;
• state and test the matched-uncertainty condition for a flight-dynamics model;
• derive an MRAC update law from a composite state-and-parameter Lyapunov function;
• explain why tracking can converge even when aerodynamic parameter estimates do not;
• apply projection, normalization, and sigma modification to a practical adaptive flight loop;
• evaluate transient recovery after abrupt or gradual aerodynamic changes.

Students are assumed to know linear state-space control, pole placement or LQR, Lyapunov equations, aircraft small-disturbance modeling, and the adaptive-control concepts developed in Chapters 1–24. The lesson does not attempt to replace a full flight-dynamics or flight-certification course.

An aircraft controller is usually designed around a family of linearizations, not one immutable plant. Aerodynamic coefficients depend on Mach number, Reynolds number, dynamic pressure, angle of attack, control-surface deflection, configuration, and aeroelastic state. Mass, inertia, and center of gravity also vary with fuel usage and payload motion. Icing, battle damage, surface jams, and structural deformation can produce changes that were not included in the nominal gain schedule.

For a longitudinal model, the pitching-moment coefficient can be represented locally as

\[ C_m = C_{m_0} + C_{m_\alpha}\alpha + C_{m_q}\frac{\bar c}{2V}q + C_{m_{\delta_e}}\delta_e + \Delta C_m(\chi,t), \]

where \( \chi \) collects flight condition and configuration. The dimensional derivatives scale with dynamic pressure and geometry. For example,

\[ M_\alpha = \frac{\bar q S\bar c}{I_y}C_{m_\alpha}, \qquad M_q = \frac{\bar q S\bar c^2}{2I_yV}C_{m_q}, \qquad M_{\delta_e} = \frac{\bar q S\bar c}{I_y}C_{m_{\delta_e}}. \]

Consequently, even unchanged nondimensional coefficient slopes can produce different state matrices when \( \bar q \), \(V\), or \(I_y\) changes. An adaptive augmentation is therefore designed to compensate the residual uncertainty left after scheduling and baseline feedback, rather than to discard all prior modeling.

Consider the short-period states \( \mathbf{x}=[\alpha\;q]^\mathsf{T} \) and elevator input \(u=\delta_e\). Around a trimmed operating point,

\[ \dot{\mathbf{x}} = \mathbf{A}(\chi,t)\mathbf{x} + \mathbf{B}(\chi,t)u + \mathbf{B}_r r + \mathbf{d}(t), \]

\[ \mathbf{A}(\chi,t)= \begin{bmatrix} Z_\alpha(\chi,t) & 1+Z_q(\chi,t)\\ M_\alpha(\chi,t) & M_q(\chi,t) \end{bmatrix}, \qquad \mathbf{B}(\chi,t)= \begin{bmatrix}Z_{\delta_e}(\chi,t)\\M_{\delta_e}(\chi,t)\end{bmatrix}. \]

Decompose the plant into nominal and uncertain parts:

\[ \mathbf{A}(\chi,t)=\mathbf{A}_0+\Delta\mathbf{A}(\chi,t), \qquad \mathbf{B}(\chi,t)=\mathbf{B}_0\Lambda(\chi,t) +\Delta\mathbf{B}_\perp(\chi,t). \]

The scalar or diagonal matrix \(\Lambda\) represents control effectiveness. The term \(\Delta\mathbf{B}_\perp\) represents uncertainty outside the nominal input direction. Classical direct MRAC cancels matched uncertainty; unmatched uncertainty requires robust baseline authority, additional dynamics, control allocation, or a different adaptive architecture.

Let the baseline controller render the nominal closed loop \(\mathbf{A}_m\) Hurwitz. A useful adaptive representation is

\[ \dot{\mathbf{x}} = \mathbf{A}_m\mathbf{x}+\mathbf{B}_r r +\mathbf{B}\left(\mathbf{W}^{*\mathsf{T}}(t)\boldsymbol{\phi}(\mathbf{x},r) +u_{ad}\right)+\mathbf{d}_u(t), \]

where \(\mathbf{W}^*(t)\) is the ideal uncertainty parameter, \(\boldsymbol{\phi}\) is a known regressor, and \(\mathbf{d}_u\) contains unmatched and unmodeled effects. A common flight regressor is

\[ \boldsymbol{\phi}=\begin{bmatrix} \alpha & q & r & 1 \end{bmatrix}^{\mathsf{T}}, \]

possibly augmented by filtered signals, angle-of-attack nonlinearities, or scheduled basis functions. The exact matching condition for matrix uncertainty is

\[ \Delta\mathbf{A}=\mathbf{B}\mathbf{W}_x^{*\mathsf{T}}, \qquad \Delta\mathbf{B}_r=\mathbf{B}\mathbf{W}_r^{*\mathsf{T}}. \]

Equivalently, every column of the uncertainty to be cancelled must lie in the column space of \(\mathbf{B}\). For a single elevator channel, this is restrictive. In practice, the baseline loop and control allocation are designed so that the dominant uncertain moment enters through the available control direction, while residual unmatched effects remain bounded.

The total command is partitioned as \(u=u_{bl}+u_{ad}\). The baseline component may be an LQR, dynamic-inversion, pole-placement, or scheduled flight-control law. The adaptive term is

\[ u_{ad}=-\widehat{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi}. \]

The desired response is generated by

\[ \dot{\mathbf{x}}_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{B}_r r, \qquad \mathbf{A}_m\text{ Hurwitz}. \]

Define tracking and parameter errors as \(\mathbf{e}=\mathbf{x}-\mathbf{x}_m\) and \(\widetilde{\mathbf{W}}=\widehat{\mathbf{W}}-\mathbf{W}^*\). Under exact matching and constant ideal parameters,

\[ \dot{\mathbf{e}}=\mathbf{A}_m\mathbf{e} -\mathbf{B}\widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi}. \]

flowchart TD
  R["Pilot or guidance command r"] --> RM["Reference model"]
  R --> BL["Baseline flight controller"]
  X["Measured aircraft states"] --> BL
  X --> REG["Regressor phi"]
  RM --> ERR["Tracking error e"]
  X --> ERR
  ERR --> LAW["Adaptive update law"]
  REG --> LAW
  LAW --> UAD["Adaptive command u_ad"]
  BL --> SUM["Command summation"]
  UAD --> SUM
  SUM --> LIM["Rate and position limits"]
  LIM --> ACT["Actuator and aircraft"]
  ACT --> X
        

The augmentation should not be allowed to erase actuator limits, structural-load limits, or protection logic. Adaptation operates inside a supervised control architecture.

Choose any \(\mathbf{Q}=\mathbf{Q}^{\mathsf{T}}\succ0\) and solve the Lyapunov equation

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

Use the composite candidate

\[ V=\mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{e} +\operatorname{tr}\!\left( \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1} \widetilde{\mathbf{W}}\right), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\mathsf{T}}\succ0. \]

For constant \(\mathbf{W}^*\), differentiation gives

\[ \begin{aligned} \dot V &=\mathbf{e}^{\mathsf{T}} (\mathbf{A}_m^{\mathsf{T}}\mathbf{P}+\mathbf{P}\mathbf{A}_m)\mathbf{e} -2\mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B} \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi}\\ &\quad+2\operatorname{tr}\!\left( \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1} \dot{\widehat{\mathbf{W}}}\right). \end{aligned} \]

Select the update law

\[ \boxed{\dot{\widehat{\mathbf{W}}} =\boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B}}. \]

The trace identity \(\operatorname{tr}(\widetilde{\mathbf{W}}^{\mathsf{T}} \boldsymbol{\phi}\mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B}) =\mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B} \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi}\) cancels the cross term, yielding

\[ \dot V=-\mathbf{e}^{\mathsf{T}}\mathbf{Q}\mathbf{e}\leq0. \]

Therefore \(\mathbf{e}\) and \(\widetilde{\mathbf{W}}\) are bounded, and \(\mathbf{e}\in L_2\cap L_\infty\). If the closed-loop signals make \(\dot{\mathbf{e}}\) bounded, Barbalat's lemma gives \(\mathbf{e}(t)\) converging to zero. This theorem proves tracking, not parameter convergence. Convergence of \(\widehat{\mathbf{W}}\) to \(\mathbf{W}^*\) requires sufficient excitation and identifiability.

During fuel burn, icing growth, aeroelastic deformation, or configuration transition, \(\mathbf{W}^*(t)\) is not constant. Then \(\dot{\widetilde{\mathbf{W}}} =\dot{\widehat{\mathbf{W}}}-\dot{\mathbf{W}}^*\), and the ideal cancellation proof acquires the disturbance term

\[ -2\operatorname{tr}\!\left( \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1} \dot{\mathbf{W}}^*\right). \]

If \(\|\dot{\mathbf{W}}^*\|_F\leq\bar w_d\), the system is generally shown to be uniformly ultimately bounded rather than asymptotically convergent. A standard sigma-modified law is

\[ \dot{\widehat{\mathbf{W}}} =\boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B} -\sigma\boldsymbol{\Gamma}\widehat{\mathbf{W}}, \qquad \sigma>0. \]

The leakage term suppresses drift under noise and unmodeled dynamics but creates bias. A representative bound has the structure

\[ \dot V\leq-\lambda_{\min}(\mathbf{Q})\|\mathbf{e}\|^2 -\sigma\|\widetilde{\mathbf{W}}\|_F^2+c_1\bar w_d^2+c_2\sigma\|\mathbf{W}^*\|_F^2. \]

Other practical modifications are:

• projection: keeps estimates inside a physically justified compact set;
• normalization: divides the update by a signal-energy term to avoid excessive learning under large regressors;
• dead zone: freezes adaptation when error is consistent with the noise floor;
• e-modification: scales leakage with tracking-error magnitude;
• filtered regressors: reduce sensitivity to sensor noise and neglected actuator dynamics.

\[ \dot{\widehat{\mathbf{W}}} =\operatorname{Proj}\!\left( \widehat{\mathbf{W}}, \frac{\boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B}} {1+\nu\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\phi}} -\sigma\boldsymbol{\Gamma}\widehat{\mathbf{W}}\right). \]

Elevator position and rate limits break the ideal input channel assumed in the Lyapunov proof. Let the commanded and realized actuator signals be \(u_c\) and \(u_a\). The deficiency is

\[ \Delta u=u_a-u_c. \]

If adaptation interprets saturation-induced tracking error as aerodynamic uncertainty, parameter estimates can grow in the wrong direction. A hedged reference model can subtract the unrealized pseudo-control:

\[ \dot{\mathbf{x}}_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{B}_r r +\mathbf{B}\Delta u. \]

This does not create actuator authority; it prevents the adaptive loop from trying to learn an impossible response. A deployable architecture also includes command limiting, anti-windup, load-factor and angle-of-attack protection, estimator health checks, adaptation-rate limits, and a monitored reversion path to the baseline controller.

Unmatched uncertainty enters outside the span of \(\mathbf{B}\). With \(\dot{\mathbf{e}}=\mathbf{A}_m\mathbf{e} -\mathbf{B}\widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi} +\mathbf{d}_u\), the best classical result is normally an input-to-state or ultimate bound proportional to \(\|\mathbf{d}_u\|_\infty\). The design goal is therefore bounded performance degradation, not exact cancellation of every damaged-aircraft mode.

For a multiaxis aircraft with state \(\mathbf{x}\in\mathbb{R}^n\), control \(\mathbf{u}\in\mathbb{R}^m\), and regressor \(\boldsymbol{\phi}\in\mathbb{R}^p\), use

\[ \mathbf{u}_{ad}=-\widehat{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\phi}, \qquad \widehat{\mathbf{W}}\in\mathbb{R}^{p\times m}, \]

\[ \dot{\widehat{\mathbf{W}}} =\boldsymbol{\Gamma}\boldsymbol{\phi} \mathbf{e}^{\mathsf{T}}\mathbf{P}\mathbf{B}. \]

Control allocation maps the requested moment vector to redundant effectors. If damage changes the allocation matrix, the adaptive inner loop and allocator must be coordinated; otherwise, the adaptive controller may request moments that the degraded effector set cannot produce.

flowchart TD
  A["Define flight envelope and failure set"] --> B["Build nominal scheduled plant family"]
  B --> C["Design baseline controller and protection logic"]
  C --> D["Identify matched uncertainty channels"]
  D --> E["Choose regressors and reference model"]
  E --> F["Derive Lyapunov update law"]
  F --> G["Add projection normalization and leakage"]
  G --> H["Test actuator limits delays noise and unmatched modes"]
  H --> I["Monte Carlo and worst-case nonlinear simulation"]
  I --> J["Hardware-in-the-loop and piloted evaluation"]
  J --> K["Define monitors reversion and evidence"]
        

The mathematical proof covers an idealized model class. Verification must then quantify how far the implemented system departs from that class. Important sweeps include aerodynamic derivatives, inertia, center of gravity, actuator bandwidth, delay, sensor noise, sampling, turbulence, flexible modes, saturation, and command aggressiveness.

The implementations below use a normalized two-state longitudinal example. The nominal reference dynamics are

\[ \mathbf{A}_m=\begin{bmatrix}-1.2&1.0\\-2.5&-1.8\end{bmatrix}, \qquad \mathbf{B}=\mathbf{B}_r=\begin{bmatrix}0\\2\end{bmatrix}. \]

With \(\boldsymbol{\phi}=[x_1\;x_2\;1]^{\mathsf{T}}\), the ideal uncertainty changes at 10 and 20 seconds:

\[ \mathbf{W}^*(t)= \begin{cases} [0.25\;-0.15\;0.02]^{\mathsf{T}}, & 0\leq t<10,\\ [0.80\;-0.55\;0.12]^{\mathsf{T}}, & 10\leq t<20,\\ [0.45\;-0.25\;0.05]^{\mathsf{T}}, & t\geq20. \end{cases} \]

The code uses \(\Gamma=18\), \(\sigma=0.08\), a componentwise projection bound of 3, and the Lyapunov solution

\[ \mathbf{P}=\begin{bmatrix} 0.50608011&-0.04291845\\-0.04291845&0.25393419 \end{bmatrix}. \]

This is a pedagogical normalized model, not a flight-qualified aircraft model. Its purpose is to expose the adaptation mechanism with identical equations across languages.

Dependencies: numpy and matplotlib. The simulation uses a fixed-step RK4 integrator, writes CSV data, and saves plots.

Chapter26_Lesson1.py

"""Chapter 26, Lesson 1: Adaptive flight control for changing aerodynamic parameters.

A compact two-state longitudinal MRAC example. The matched uncertainty vector is
changed twice to represent aerodynamic changes caused by a new flight condition or
loss of control effectiveness. The adaptive augmentation attempts to preserve the
reference-model response.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 0.002
    final_time: float = 30.0
    gamma: float = 18.0
    sigma: float = 0.08
    parameter_limit: float = 3.0


A_M = np.array([[-1.2, 1.0], [-2.5, -1.8]], dtype=float)
B = np.array([0.0, 2.0], dtype=float)
B_R = np.array([0.0, 2.0], dtype=float)
P = np.array([[0.50608011, -0.04291845], [-0.04291845, 0.25393419]])


def command(t: float) -> float:
    """Bounded pitch-rate command in rad/s."""
    step = 0.08 if t >= 4.0 else 0.0
    return step + 0.12 * np.sin(0.7 * t)


def true_uncertainty(t: float) -> np.ndarray:
    """Piecewise aerodynamic uncertainty [alpha, q, bias] coefficients."""
    if t < 10.0:
        return np.array([0.25, -0.15, 0.02])
    if t < 20.0:
        return np.array([0.80, -0.55, 0.12])
    return np.array([0.45, -0.25, 0.05])


def project_box(theta: np.ndarray, limit: float) -> np.ndarray:
    """Simple componentwise projection used only as a numerical safety bound."""
    return np.clip(theta, -limit, limit)


def derivative(t: float, z: np.ndarray, cfg: Config) -> np.ndarray:
    """Continuous-time plant, reference model, and adaptive-law dynamics."""
    x = z[0:2]
    x_m = z[2:4]
    theta_hat = z[4:7]

    phi = np.array([x[0], x[1], 1.0])
    theta_star = true_uncertainty(t)
    r = command(t)

    # u_ad = -theta_hat^T phi. The true matched uncertainty is theta_star^T phi.
    residual = float((theta_star - theta_hat) @ phi)
    x_dot = A_M @ x + B_R * r + B * residual
    x_m_dot = A_M @ x_m + B_R * r

    e = x - x_m
    s = float(e @ P @ B)
    theta_dot = cfg.gamma * phi * s - cfg.sigma * cfg.gamma * theta_hat

    # Stop an estimate from integrating farther outward at a projection boundary.
    for i, value in enumerate(theta_hat):
        if value >= cfg.parameter_limit and theta_dot[i] > 0.0:
            theta_dot[i] = 0.0
        elif value <= -cfg.parameter_limit and theta_dot[i] < 0.0:
            theta_dot[i] = 0.0

    return np.concatenate((x_dot, x_m_dot, theta_dot))


def rk4_step(t: float, z: np.ndarray, dt: float, cfg: Config) -> np.ndarray:
    k1 = derivative(t, z, cfg)
    k2 = derivative(t + 0.5 * dt, z + 0.5 * dt * k1, cfg)
    k3 = derivative(t + 0.5 * dt, z + 0.5 * dt * k2, cfg)
    k4 = derivative(t + dt, z + dt * k3, cfg)
    z_next = z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    z_next[4:7] = project_box(z_next[4:7], cfg.parameter_limit)
    return z_next


def simulate(cfg: Config) -> dict[str, np.ndarray]:
    times = np.arange(0.0, cfg.final_time + cfg.dt, cfg.dt)
    state = np.zeros(7)

    log = np.zeros((times.size, 15))
    for k, t in enumerate(times):
        theta_star = true_uncertainty(t)
        r = command(t)
        x = state[0:2]
        x_m = state[2:4]
        theta_hat = state[4:7]
        phi = np.array([x[0], x[1], 1.0])
        u_ad = -float(theta_hat @ phi)
        error = x - x_m

        log[k] = np.array(
            [
                t,
                r,
                x[0],
                x[1],
                x_m[0],
                x_m[1],
                error[0],
                error[1],
                u_ad,
                *theta_hat,
                *theta_star,
            ]
        )
        if k + 1 < times.size:
            state = rk4_step(t, state, cfg.dt, cfg)

    return {
        "time": log[:, 0],
        "command": log[:, 1],
        "alpha": log[:, 2],
        "q": log[:, 3],
        "alpha_m": log[:, 4],
        "q_m": log[:, 5],
        "e_alpha": log[:, 6],
        "e_q": log[:, 7],
        "u_ad": log[:, 8],
        "theta_hat": log[:, 9:12],
        "theta_star": log[:, 12:15],
        "raw": log,
    }


def save_results(results: dict[str, np.ndarray], output_dir: Path) -> None:
    output_dir.mkdir(parents=True, exist_ok=True)
    header = (
        "time,command,alpha,q,alpha_reference,q_reference,e_alpha,e_q,"
        "u_ad,theta_hat_alpha,theta_hat_q,theta_hat_bias,"
        "theta_true_alpha,theta_true_q,theta_true_bias"
    )
    np.savetxt(output_dir / "Chapter26_Lesson1_results.csv", results["raw"], delimiter=",", header=header, comments="")

    t = results["time"]
    fig = plt.figure(figsize=(10, 8))
    ax1 = fig.add_subplot(3, 1, 1)
    ax1.plot(t, results["q"], label="aircraft q")
    ax1.plot(t, results["q_m"], "--", label="reference q_m")
    ax1.set_ylabel("pitch rate (rad/s)")
    ax1.grid(True)
    ax1.legend()

    ax2 = fig.add_subplot(3, 1, 2)
    ax2.plot(t, results["e_alpha"], label="e_alpha")
    ax2.plot(t, results["e_q"], label="e_q")
    ax2.set_ylabel("tracking error")
    ax2.grid(True)
    ax2.legend()

    ax3 = fig.add_subplot(3, 1, 3)
    labels = ["alpha", "q", "bias"]
    for i, label in enumerate(labels):
        ax3.plot(t, results["theta_hat"][:, i], label=f"hat theta {label}")
        ax3.plot(t, results["theta_star"][:, i], "--", label=f"true theta {label}")
    ax3.set_xlabel("time (s)")
    ax3.set_ylabel("parameter")
    ax3.grid(True)
    ax3.legend(ncol=2)

    fig.tight_layout()
    fig.savefig(output_dir / "Chapter26_Lesson1_results.png", dpi=180)
    plt.close(fig)


def main() -> None:
    cfg = Config()
    results = simulate(cfg)
    output_dir = Path(__file__).resolve().parent / "output"
    save_results(results, output_dir)

    e_norm = np.sqrt(results["e_alpha"] ** 2 + results["e_q"] ** 2)
    rms = float(np.sqrt(np.mean(e_norm**2)))
    peak_after_change = float(np.max(e_norm[results["time"] >= 10.0]))
    print(f"RMS state-tracking error: {rms:.6f}")
    print(f"Peak error after aerodynamic changes: {peak_after_change:.6f}")
    print(f"Results written to: {output_dir}")


if __name__ == "__main__":
    main()

This version requires only a C++17 compiler and writes a CSV file suitable for plotting in Python, MATLAB, GNU Octave, or a spreadsheet.

Chapter26_Lesson1.cpp

// Chapter 26, Lesson 1: Adaptive flight control for changing aerodynamic parameters.
// Standard-library implementation of the same two-state MRAC simulation used in
// the Python, Java, MATLAB, and Mathematica examples.

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

namespace {
constexpr double DT = 0.002;
constexpr double FINAL_TIME = 30.0;
constexpr double GAMMA = 18.0;
constexpr double SIGMA = 0.08;
constexpr double PARAMETER_LIMIT = 3.0;

using State = std::array<double, 7>;  // x(2), xm(2), theta_hat(3)
using Vec3 = std::array<double, 3>;

constexpr std::array<std::array<double, 2>, 2> A_M{{
    {{-1.2, 1.0}},
    {{-2.5, -1.8}},
}};
constexpr std::array<double, 2> B{{0.0, 2.0}};
constexpr std::array<double, 2> B_R{{0.0, 2.0}};
constexpr std::array<std::array<double, 2>, 2> P{{
    {{0.50608011, -0.04291845}},
    {{-0.04291845, 0.25393419}},
}};

double command(double t) {
    const double step = (t >= 4.0) ? 0.08 : 0.0;
    return step + 0.12 * std::sin(0.7 * t);
}

Vec3 trueUncertainty(double t) {
    if (t < 10.0) return {{0.25, -0.15, 0.02}};
    if (t < 20.0) return {{0.80, -0.55, 0.12}};
    return {{0.45, -0.25, 0.05}};
}

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(double t, const State& z) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const Vec3 thetaHat{{z[4], z[5], z[6]}};
    const Vec3 phi{{x1, x2, 1.0}};
    const Vec3 thetaStar = trueUncertainty(t);
    const double r = command(t);

    double residual = 0.0;
    for (std::size_t i = 0; i < 3; ++i) residual += (thetaStar[i] - thetaHat[i]) * phi[i];

    State dz{};
    dz[0] = A_M[0][0] * x1 + A_M[0][1] * x2 + B_R[0] * r + B[0] * residual;
    dz[1] = A_M[1][0] * x1 + A_M[1][1] * x2 + B_R[1] * r + B[1] * residual;
    dz[2] = A_M[0][0] * xm1 + A_M[0][1] * xm2 + B_R[0] * r;
    dz[3] = A_M[1][0] * xm1 + A_M[1][1] * xm2 + B_R[1] * r;

    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double pb1 = P[0][0] * B[0] + P[0][1] * B[1];
    const double pb2 = P[1][0] * B[0] + P[1][1] * B[1];
    const double s = e1 * pb1 + e2 * pb2;

    for (std::size_t i = 0; i < 3; ++i) {
        dz[4 + i] = GAMMA * phi[i] * s - SIGMA * GAMMA * thetaHat[i];
        if (thetaHat[i] >= PARAMETER_LIMIT && dz[4 + i] > 0.0) dz[4 + i] = 0.0;
        if (thetaHat[i] <= -PARAMETER_LIMIT && dz[4 + i] < 0.0) dz[4 + i] = 0.0;
    }
    return dz;
}

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

    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]);
    }
    for (std::size_t i = 4; i < 7; ++i) {
        next[i] = std::clamp(next[i], -PARAMETER_LIMIT, PARAMETER_LIMIT);
    }
    return next;
}
}  // namespace

int main(int argc, char* argv[]) {
    try {
        const std::string outputPath = (argc > 1) ? argv[1] : "Chapter26_Lesson1_results_cpp.csv";
        std::ofstream csv(outputPath);
        if (!csv) throw std::runtime_error("Cannot open output file: " + outputPath);

        csv << "time,command,alpha,q,alpha_reference,q_reference,e_alpha,e_q,u_ad,"
               "theta_hat_alpha,theta_hat_q,theta_hat_bias,theta_true_alpha,theta_true_q,theta_true_bias\n";
        csv << std::setprecision(12);

        State state{};
        double squaredErrorSum = 0.0;
        double peakAfterChange = 0.0;
        std::size_t sampleCount = 0;

        const std::size_t steps = static_cast<std::size_t>(std::llround(FINAL_TIME / DT));
        for (std::size_t k = 0; k <= steps; ++k) {
            const double t = static_cast<double>(k) * DT;
            const Vec3 thetaStar = trueUncertainty(t);
            const Vec3 phi{{state[0], state[1], 1.0}};
            const double uAd = -(state[4] * phi[0] + state[5] * phi[1] + state[6]);
            const double e1 = state[0] - state[2];
            const double e2 = state[1] - state[3];
            const double eNorm = std::sqrt(e1 * e1 + e2 * e2);

            squaredErrorSum += eNorm * eNorm;
            if (t >= 10.0) peakAfterChange = std::max(peakAfterChange, eNorm);
            ++sampleCount;

            csv << t << ',' << command(t) << ',' << state[0] << ',' << state[1] << ','
                << state[2] << ',' << state[3] << ',' << e1 << ',' << e2 << ',' << uAd << ','
                << state[4] << ',' << state[5] << ',' << state[6] << ','
                << thetaStar[0] << ',' << thetaStar[1] << ',' << thetaStar[2] << '\n';

            if (k < steps) state = rk4Step(t, state);
        }

        const double rms = std::sqrt(squaredErrorSum / static_cast<double>(sampleCount));
        std::cout << "RMS state-tracking error: " << rms << '\n';
        std::cout << "Peak error after aerodynamic changes: " << peakAfterChange << '\n';
        std::cout << "Results written to: " << outputPath << '\n';
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

The Java implementation is dependency-free and uses the same RK4 equations and test schedule.

Chapter26_Lesson1.java

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

/**
 * Chapter 26, Lesson 1: Adaptive flight control for changing aerodynamic parameters.
 *
 * A dependency-free Java implementation of a two-state MRAC simulation. The true
 * matched aerodynamic uncertainty changes at t = 10 s and t = 20 s.
 */
public final class Chapter26_Lesson1 {
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 30.0;
    private static final double GAMMA = 18.0;
    private static final double SIGMA = 0.08;
    private static final double PARAMETER_LIMIT = 3.0;

    private static final double[][] A_M = {
        {-1.2, 1.0},
        {-2.5, -1.8}
    };
    private static final double[] B = {0.0, 2.0};
    private static final double[] B_R = {0.0, 2.0};
    private static final double[][] P = {
        {0.50608011, -0.04291845},
        {-0.04291845, 0.25393419}
    };

    private Chapter26_Lesson1() {
    }

    private static double command(double t) {
        double step = t >= 4.0 ? 0.08 : 0.0;
        return step + 0.12 * Math.sin(0.7 * t);
    }

    private static double[] trueUncertainty(double t) {
        if (t < 10.0) {
            return new double[] {0.25, -0.15, 0.02};
        }
        if (t < 20.0) {
            return new double[] {0.80, -0.55, 0.12};
        }
        return new double[] {0.45, -0.25, 0.05};
    }

    private static double[] derivative(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double[] thetaHat = {z[4], z[5], z[6]};
        double[] phi = {x1, x2, 1.0};
        double[] thetaStar = trueUncertainty(t);
        double r = command(t);

        double residual = 0.0;
        for (int i = 0; i < 3; i++) {
            residual += (thetaStar[i] - thetaHat[i]) * phi[i];
        }

        double[] dz = new double[7];
        dz[0] = A_M[0][0] * x1 + A_M[0][1] * x2 + B_R[0] * r + B[0] * residual;
        dz[1] = A_M[1][0] * x1 + A_M[1][1] * x2 + B_R[1] * r + B[1] * residual;
        dz[2] = A_M[0][0] * xm1 + A_M[0][1] * xm2 + B_R[0] * r;
        dz[3] = A_M[1][0] * xm1 + A_M[1][1] * xm2 + B_R[1] * r;

        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double pb1 = P[0][0] * B[0] + P[0][1] * B[1];
        double pb2 = P[1][0] * B[0] + P[1][1] * B[1];
        double s = e1 * pb1 + e2 * pb2;

        for (int i = 0; i < 3; i++) {
            dz[4 + i] = GAMMA * phi[i] * s - SIGMA * GAMMA * thetaHat[i];
            if (thetaHat[i] >= PARAMETER_LIMIT && dz[4 + i] > 0.0) {
                dz[4 + i] = 0.0;
            } else if (thetaHat[i] <= -PARAMETER_LIMIT && dz[4 + i] < 0.0) {
                dz[4 + i] = 0.0;
            }
        }
        return dz;
    }

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

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

        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]);
        }
        for (int i = 4; i < 7; i++) {
            next[i] = Math.max(-PARAMETER_LIMIT, Math.min(PARAMETER_LIMIT, next[i]));
        }
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Path outputPath = args.length > 0
            ? Paths.get(args[0])
            : Paths.get("Chapter26_Lesson1_results_java.csv");

        try {
            Path parent = outputPath.toAbsolutePath().getParent();
            if (parent != null) {
                Files.createDirectories(parent);
            }

            try (BufferedWriter writer = Files.newBufferedWriter(outputPath)) {
                writer.write("time,command,alpha,q,alpha_reference,q_reference,e_alpha,e_q,u_ad,"
                    + "theta_hat_alpha,theta_hat_q,theta_hat_bias,theta_true_alpha,theta_true_q,theta_true_bias\n");

                double[] state = new double[7];
                double squaredErrorSum = 0.0;
                double peakAfterChange = 0.0;
                long samples = 0;
                int steps = (int) Math.round(FINAL_TIME / DT);

                for (int k = 0; k <= steps; k++) {
                    double t = k * DT;
                    double[] thetaStar = trueUncertainty(t);
                    double[] phi = {state[0], state[1], 1.0};
                    double uAd = -(state[4] * phi[0] + state[5] * phi[1] + state[6]);
                    double e1 = state[0] - state[2];
                    double e2 = state[1] - state[3];
                    double eNorm = Math.hypot(e1, e2);

                    squaredErrorSum += eNorm * eNorm;
                    if (t >= 10.0) {
                        peakAfterChange = Math.max(peakAfterChange, eNorm);
                    }
                    samples++;

                    writer.write(String.format(Locale.ROOT,
                        "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                        t, command(t), state[0], state[1], state[2], state[3], e1, e2, uAd,
                        state[4], state[5], state[6], thetaStar[0], thetaStar[1], thetaStar[2]));

                    if (k < steps) {
                        state = rk4Step(t, state);
                    }
                }

                double rms = Math.sqrt(squaredErrorSum / samples);
                System.out.printf(Locale.ROOT, "RMS state-tracking error: %.7f%n", rms);
                System.out.printf(Locale.ROOT, "Peak error after aerodynamic changes: %.7f%n", peakAfterChange);
                System.out.println("Results written to: " + outputPath.toAbsolutePath());
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

The script performs the simulation directly, exports a table, and plots tracking errors and parameter estimates.

Chapter26_Lesson1.m

%% Chapter26_Lesson1.m
% Adaptive flight control for changing aerodynamic parameters.
% Two-state longitudinal model-reference adaptive-control simulation.

clear; clc; close all;

cfg.dt = 0.002;
cfg.finalTime = 30.0;
cfg.gamma = 18.0;
cfg.sigma = 0.08;
cfg.parameterLimit = 3.0;

A_m = [-1.2, 1.0; -2.5, -1.8];
B = [0.0; 2.0];
B_r = [0.0; 2.0];
P = [0.50608011, -0.04291845; -0.04291845, 0.25393419];

time = (0:cfg.dt:cfg.finalTime).';
state = zeros(7, 1); % [x(2); x_m(2); theta_hat(3)]
logData = zeros(numel(time), 15);

for k = 1:numel(time)
    t = time(k);
    thetaStar = trueUncertainty(t);
    r = commandSignal(t);
    x = state(1:2);
    xm = state(3:4);
    thetaHat = state(5:7);
    phi = [x(1); x(2); 1.0];
    uAd = -thetaHat.' * phi;
    e = x - xm;

    logData(k, :) = [t, r, x.', xm.', e.', uAd, thetaHat.', thetaStar.'];

    if k < numel(time)
        state = rk4Step(t, state, cfg, A_m, B, B_r, P);
    end
end

names = {'time','command','alpha','q','alpha_reference','q_reference', ...
    'e_alpha','e_q','u_ad','theta_hat_alpha','theta_hat_q','theta_hat_bias', ...
    'theta_true_alpha','theta_true_q','theta_true_bias'};
results = array2table(logData, 'VariableNames', names);
writetable(results, 'Chapter26_Lesson1_results_matlab.csv');

errorNorm = hypot(results.e_alpha, results.e_q);
rmsError = sqrt(mean(errorNorm.^2));
peakAfterChange = max(errorNorm(results.time >= 10.0));
fprintf('RMS state-tracking error: %.7f\n', rmsError);
fprintf('Peak error after aerodynamic changes: %.7f\n', peakAfterChange);

figure('Name', 'Chapter 26 Lesson 1: Adaptive Flight Control');
tiledlayout(3,1);

nexttile;
plot(time, results.q, 'LineWidth', 1.2); hold on;
plot(time, results.q_reference, '--', 'LineWidth', 1.2);
grid on; ylabel('pitch rate (rad/s)');
legend('aircraft q','reference q_m','Location','best');

nexttile;
plot(time, results.e_alpha, 'LineWidth', 1.1); hold on;
plot(time, results.e_q, 'LineWidth', 1.1);
grid on; ylabel('tracking error');
legend('e_alpha','e_q','Location','best');

nexttile;
plot(time, results.theta_hat_alpha, 'LineWidth', 1.0); hold on;
plot(time, results.theta_true_alpha, '--', 'LineWidth', 1.0);
plot(time, results.theta_hat_q, 'LineWidth', 1.0);
plot(time, results.theta_true_q, '--', 'LineWidth', 1.0);
plot(time, results.theta_hat_bias, 'LineWidth', 1.0);
plot(time, results.theta_true_bias, '--', 'LineWidth', 1.0);
grid on; xlabel('time (s)'); ylabel('parameter');
legend('hat theta alpha','true theta alpha','hat theta q','true theta q', ...
    'hat theta bias','true theta bias','Location','best');

exportgraphics(gcf, 'Chapter26_Lesson1_results_matlab.png', 'Resolution', 180);

%% Local functions
function r = commandSignal(t)
    step = 0.08 * double(t >= 4.0);
    r = step + 0.12 * sin(0.7 * t);
end

function theta = trueUncertainty(t)
    if t < 10.0
        theta = [0.25; -0.15; 0.02];
    elseif t < 20.0
        theta = [0.80; -0.55; 0.12];
    else
        theta = [0.45; -0.25; 0.05];
    end
end

function dz = derivative(t, z, cfg, A_m, B, B_r, P)
    x = z(1:2);
    xm = z(3:4);
    thetaHat = z(5:7);
    phi = [x(1); x(2); 1.0];
    thetaStar = trueUncertainty(t);
    r = commandSignal(t);

    residual = (thetaStar - thetaHat).' * phi;
    xDot = A_m * x + B_r * r + B * residual;
    xmDot = A_m * xm + B_r * r;

    e = x - xm;
    s = e.' * P * B;
    thetaDot = cfg.gamma * phi * s - cfg.sigma * cfg.gamma * thetaHat;

    outwardUpper = thetaHat >= cfg.parameterLimit & thetaDot > 0;
    outwardLower = thetaHat <= -cfg.parameterLimit & thetaDot < 0;
    thetaDot(outwardUpper | outwardLower) = 0;

    dz = [xDot; xmDot; thetaDot];
end

function next = rk4Step(t, z, cfg, A_m, B, B_r, P)
    dt = cfg.dt;
    k1 = derivative(t, z, cfg, A_m, B, B_r, P);
    k2 = derivative(t + dt/2, z + dt*k1/2, cfg, A_m, B, B_r, P);
    k3 = derivative(t + dt/2, z + dt*k2/2, cfg, A_m, B, B_r, P);
    k4 = derivative(t + dt, z + dt*k3, cfg, A_m, B, B_r, P);
    next = z + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
    next(5:7) = max(-cfg.parameterLimit, min(cfg.parameterLimit, next(5:7)));
end

Running this MATLAB script creates and simulates Chapter26_Lesson1_Simulink.slx. It uses a fixed-step discrete MATLAB Function block so that the complete example can be generated without a manually prepared model file.

Chapter26_Lesson1_Simulink.m

%% Chapter26_Lesson1_Simulink.m
% Programmatically builds and runs a fixed-step Simulink demonstration of the
% Chapter 26 Lesson 1 adaptive flight-control example.
%
% Requirements: MATLAB and Simulink. The script creates
% Chapter26_Lesson1_Simulink.slx in the current directory.

clear; clc;
model = 'Chapter26_Lesson1_Simulink';

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

new_system(model);
open_system(model);
set_param(model, ...
    'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '0.002', ...
    'StopTime', '30', ...
    'ReturnWorkspaceOutputs', 'on');

add_block('simulink/Sources/Clock', [model '/Clock'], ...
    'Position', [45 85 75 115]);
add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/Adaptive Flight MRAC'], ...
    'Position', [150 55 355 145]);
add_block('simulink/Sinks/To Workspace', [model '/Log'], ...
    'VariableName', 'adaptiveFlightLog', ...
    'SaveFormat', 'Array', ...
    'MaxDataPoints', 'inf', ...
    'Position', [440 85 540 115]);

add_line(model, 'Clock/1', 'Adaptive Flight MRAC/1', 'autorouting', 'on');
add_line(model, 'Adaptive Flight MRAC/1', 'Log/1', 'autorouting', 'on');

root = sfroot;
chart = root.find('-isa', 'Stateflow.EMChart', 'Path', [model '/Adaptive Flight MRAC']);
if isempty(chart)
    error('Could not access the MATLAB Function block through the Stateflow API.');
end

chart.Script = sprintf([ ...
'function y = fcn(t)\n' ...
'%%#codegen\n' ...
'%% y = [t r alpha q alpha_m q_m e_alpha e_q u_ad thetaHat(1:3) thetaStar(1:3)]\n' ...
'persistent z lastT\n' ...
'if isempty(z)\n' ...
'    z = zeros(7,1);\n' ...
'    lastT = t;\n' ...
'end\n' ...
'dt = t - lastT;\n' ...
'if dt > 0\n' ...
'    dz = localDerivative(t, z);\n' ...
'    z = z + dt*dz;\n' ...
'    z(5:7) = min(3.0, max(-3.0, z(5:7)));\n' ...
'    lastT = t;\n' ...
'end\n' ...
'x = z(1:2); xm = z(3:4); thetaHat = z(5:7);\n' ...
'phi = [x(1); x(2); 1.0];\n' ...
'r = localCommand(t); thetaStar = localTheta(t);\n' ...
'uAd = -thetaHat.''*phi; e = x-xm;\n' ...
'y = [t; r; x; xm; e; uAd; thetaHat; thetaStar];\n' ...
'end\n' ...
'\n' ...
'function dz = localDerivative(t, z)\n' ...
'A = [-1.2 1.0; -2.5 -1.8]; B = [0;2]; Br = [0;2];\n' ...
'P = [0.50608011 -0.04291845; -0.04291845 0.25393419];\n' ...
'gamma = 18.0; sigma = 0.08;\n' ...
'x = z(1:2); xm = z(3:4); thetaHat = z(5:7);\n' ...
'phi = [x(1); x(2); 1.0]; thetaStar = localTheta(t); r = localCommand(t);\n' ...
'residual = (thetaStar-thetaHat).''*phi;\n' ...
'xDot = A*x + Br*r + B*residual; xmDot = A*xm + Br*r;\n' ...
'e = x-xm; s = e.''*P*B;\n' ...
'parameterDot = gamma*phi*s - sigma*gamma*thetaHat;\n' ...
'for i=1:3\n' ...
'    if thetaHat(i) >= 3.0 && parameterDot(i) > 0\n' ...
'        parameterDot(i) = 0;\n' ...
'    elseif thetaHat(i) <= -3.0 && parameterDot(i) < 0\n' ...
'        parameterDot(i) = 0;\n' ...
'    end\n' ...
'end\n' ...
'dz = [xDot; xmDot; parameterDot];\n' ...
'end\n' ...
'\n' ...
'function r = localCommand(t)\n' ...
'r = 0.12*sin(0.7*t);\n' ...
'if t >= 4.0, r = r + 0.08; end\n' ...
'end\n' ...
'\n' ...
'function theta = localTheta(t)\n' ...
'if t < 10.0\n' ...
'    theta = [0.25; -0.15; 0.02];\n' ...
'elseif t < 20.0\n' ...
'    theta = [0.80; -0.55; 0.12];\n' ...
'else\n' ...
'    theta = [0.45; -0.25; 0.05];\n' ...
'end\n' ...
'end\n']);

save_system(model);
simOut = sim(model);
logData = simOut.adaptiveFlightLog;

headers = {'time','command','alpha','q','alpha_reference','q_reference', ...
    'e_alpha','e_q','u_ad','theta_hat_alpha','theta_hat_q','theta_hat_bias', ...
    'theta_true_alpha','theta_true_q','theta_true_bias'};
if size(logData, 2) == numel(headers)
    writetable(array2table(logData, 'VariableNames', headers), ...
        'Chapter26_Lesson1_results_simulink.csv');
end

figure('Name','Simulink Adaptive Flight MRAC');
plot(logData(:,1), logData(:,4), 'LineWidth', 1.2); hold on;
plot(logData(:,1), logData(:,6), '--', 'LineWidth', 1.2);
grid on; xlabel('time (s)'); ylabel('pitch rate (rad/s)');
legend('aircraft q','reference q_m','Location','best');
title('Adaptive response under two aerodynamic-parameter changes');

fprintf('Created and simulated %s.slx\n', model);

The notebook solves the coupled plant, reference-model, and adaptation equations with NDSolveValue, exports results, and produces tracking and parameter plots.

Chapter26_Lesson1.nb

(* Chapter26_Lesson1.wl
   Adaptive flight control for changing aerodynamic parameters. *)

ClearAll["Global`*"];

finalTime = 30.0;
gamma = 18.0;
sigma = 0.08;

am = {{-1.2, 1.0}, {-2.5, -1.8}};
b = {0.0, 2.0};
br = {0.0, 2.0};
p = {{0.50608011, -0.04291845}, {-0.04291845, 0.25393419}};
pb = p.b;

command[t_?NumericQ] := 0.12 Sin[0.7 t] + Piecewise[{{0.08, t >= 4.0}}, 0.0];
theta1[t_?NumericQ] := Piecewise[{{0.25, t < 10.0}, {0.80, t < 20.0}}, 0.45];
theta2[t_?NumericQ] := Piecewise[{{-0.15, t < 10.0}, {-0.55, t < 20.0}}, -0.25];
theta3[t_?NumericQ] := Piecewise[{{0.02, t < 10.0}, {0.12, t < 20.0}}, 0.05];

trackingScalar[t_] := (x1[t] - xm1[t]) pb[[1]] + (x2[t] - xm2[t]) pb[[2]];
residual[t_] := (theta1[t] - wh1[t]) x1[t] + (theta2[t] - wh2[t]) x2[t] + theta3[t] - wh3[t];

eqns = {
   x1'[t] == am[[1, 1]] x1[t] + am[[1, 2]] x2[t] + br[[1]] command[t] + b[[1]] residual[t],
   x2'[t] == am[[2, 1]] x1[t] + am[[2, 2]] x2[t] + br[[2]] command[t] + b[[2]] residual[t],
   xm1'[t] == am[[1, 1]] xm1[t] + am[[1, 2]] xm2[t] + br[[1]] command[t],
   xm2'[t] == am[[2, 1]] xm1[t] + am[[2, 2]] xm2[t] + br[[2]] command[t],
   wh1'[t] == gamma x1[t] trackingScalar[t] - sigma gamma wh1[t],
   wh2'[t] == gamma x2[t] trackingScalar[t] - sigma gamma wh2[t],
   wh3'[t] == gamma trackingScalar[t] - sigma gamma wh3[t],
   x1[0] == 0, x2[0] == 0, xm1[0] == 0, xm2[0] == 0,
   wh1[0] == 0, wh2[0] == 0, wh3[0] == 0
   };

sol = NDSolveValue[
   eqns,
   {x1, x2, xm1, xm2, wh1, wh2, wh3},
   {t, 0, finalTime},
   Method -> {"EquationSimplification" -> "Residual"},
   MaxStepSize -> 0.01
   ];

{x1s, x2s, xm1s, xm2s, wh1s, wh2s, wh3s} = sol;

sampleTimes = Range[0.0, finalTime, 0.01];
data = Table[
   With[{a = x1s[t], q = x2s[t], amRef = xm1s[t], qmRef = xm2s[t],
     h1 = wh1s[t], h2 = wh2s[t], h3 = wh3s[t]},
    {t, command[t], a, q, amRef, qmRef, a - amRef, q - qmRef,
     -(h1 a + h2 q + h3), h1, h2, h3, theta1[t], theta2[t], theta3[t]}],
   {t, sampleTimes}
   ];

header = {"time", "command", "alpha", "q", "alpha_reference", "q_reference",
   "e_alpha", "e_q", "u_ad", "theta_hat_alpha", "theta_hat_q",
   "theta_hat_bias", "theta_true_alpha", "theta_true_q", "theta_true_bias"};
Export["Chapter26_Lesson1_results_mathematica.csv", Prepend[data, header]];

rmsError = Sqrt[Mean[(data[[All, 7]]^2 + data[[All, 8]]^2)]];
peakAfterChange = Max[Map[Sqrt[#[[7]]^2 + #[[8]]^2] &, Select[data, #[[1]] >= 10.0 &]]];
Print["RMS state-tracking error: ", N[rmsError, 8]];
Print["Peak error after aerodynamic changes: ", N[peakAfterChange, 8]];

trackingPlot = Plot[
   Evaluate[{x2s[t], xm2s[t]}], {t, 0, finalTime},
   PlotLegends -> {"aircraft q", "reference q_m"},
   Frame -> True, FrameLabel -> {"time (s)", "pitch rate (rad/s)"},
   PlotRange -> All, ImageSize -> Large
   ];

parameterPlot = Plot[
   Evaluate[{wh1s[t], theta1[t], wh2s[t], theta2[t], wh3s[t], theta3[t]}],
   {t, 0, finalTime},
   PlotLegends -> {"hat theta alpha", "true theta alpha", "hat theta q",
     "true theta q", "hat theta bias", "true theta bias"},
   Frame -> True, FrameLabel -> {"time (s)", "parameter"},
   PlotRange -> All, ImageSize -> Large
   ];

Export["Chapter26_Lesson1_tracking_mathematica.png", trackingPlot, ImageResolution -> 180];
Export["Chapter26_Lesson1_parameters_mathematica.png", parameterPlot, ImageResolution -> 180];
Column[{trackingPlot, parameterPlot}]

A single tracking plot is insufficient for adaptive flight-control assessment. Useful quantitative measures include

\[ J_e=\left(\frac{1}{T}\int_0^T\mathbf{e}^{\mathsf{T}}\mathbf{e}\,dt\right)^{1/2}, \qquad J_u=\left(\frac{1}{T}\int_0^T\mathbf{u}_{ad}^{\mathsf{T}} \mathbf{u}_{ad}\,dt\right)^{1/2}, \]

\[ J_{\dot u}=\max_{t\in[0,T]}\|\dot{\mathbf{u}}(t)\|, \qquad J_\theta=\max_{t\in[0,T]}\|\widehat{\mathbf{W}}(t)\|_F. \]

Also record peak angle of attack, load factor, control position and rate margins, recovery time after each plant change, delay margin, sensitivity to flexible modes, and the number and duration of protection or reversion events. For piloted aircraft, handling-quality and pilot-induced-oscillation assessments are separate from pure reference-model tracking.

In the supplied test case, the independently compiled Python, C++, and Java versions produce an RMS state-tracking error of approximately \(0.0144\) and a peak post-change error of approximately \(0.0317\). Agreement across implementations is a useful software-verification check, not evidence of flight safety.

Problem 1 — Dynamic-pressure dependence: Suppose \(C_{m_\alpha}\), geometry, inertia, and airspeed remain fixed, but dynamic pressure rises from \(\bar q_1\) to \(\bar q_2=1.4\bar q_1\). Determine the change in \(M_\alpha\).

Solution:

\[ M_\alpha=\frac{\bar qS\bar c}{I_y}C_{m_\alpha}. \]

Every factor except \(\bar q\) is fixed, so \(M_{\alpha,2}=1.4M_{\alpha,1}\). If \(C_{m_\alpha}<0\), the derivative becomes 40 percent more negative. The short-period frequency and damping generally change, so a fixed controller no longer produces exactly the original closed-loop poles.

Problem 2 — Testing the matching condition: Let \(\mathbf{B}=[0\;2]^{\mathsf{T}}\) and \(\Delta\mathbf{A}=\begin{bmatrix}0&0\\0.8&-0.5\end{bmatrix}\). Find \(\mathbf{W}_x^*\) satisfying \(\Delta\mathbf{A}=\mathbf{B}\mathbf{W}_x^{*\mathsf{T}}\).

Solution:

Write \(\mathbf{W}_x^*=[w_1\;w_2]^{\mathsf{T}}\). Then

\[ \mathbf{B}\mathbf{W}_x^{*\mathsf{T}} =\begin{bmatrix}0\\2\end{bmatrix} \begin{bmatrix}w_1&w_2\end{bmatrix} =\begin{bmatrix}0&0\\2w_1&2w_2\end{bmatrix}. \]

Hence \(w_1=0.4\) and \(w_2=-0.25\). The uncertainty is matched because its first row is zero and both uncertain columns lie in the span of \(\mathbf{B}\).

Problem 3 — Slowly varying ideal parameters: Repeat the Lyapunov derivative when \(\mathbf{W}^*=\mathbf{W}^*(t)\) and show why exact asymptotic tracking is no longer obtained directly.

Solution:

Since \(\dot{\widetilde{\mathbf{W}}} =\dot{\widehat{\mathbf{W}}}-\dot{\mathbf{W}}^*\), the nominal update law cancels only the state-parameter cross term. Therefore

\[ \dot V=-\mathbf{e}^{\mathsf{T}}\mathbf{Q}\mathbf{e} -2\operatorname{tr}\!\left( \widetilde{\mathbf{W}}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1} \dot{\mathbf{W}}^*\right). \]

The second term has no fixed sign. Applying Cauchy–Schwarz and Young's inequality bounds it by a quadratic parameter-error term plus a constant proportional to \(\|\dot{\mathbf{W}}^*\|_F^2\). With leakage or projection, one obtains uniform ultimate boundedness whose radius increases with parameter variation rate and unmatched disturbance magnitude.

Problem 4 — Leakage bias: At zero tracking error, determine the direction of parameter motion under sigma modification.

Solution:

When \(\mathbf{e}=\mathbf{0}\), the gradient term vanishes and

\[ \dot{\widehat{\mathbf{W}}} =-\sigma\boldsymbol{\Gamma}\widehat{\mathbf{W}}. \]

Thus the estimate decays toward zero. This prevents indefinite drift, but if the true ideal parameter is nonzero, the estimate is biased away from it. The resulting residual uncertainty causes a nonzero ultimate tracking-error bound.

Problem 5 — Saturation-induced false adaptation: Explain why an elevator position limit can drive parameter estimates even when the aerodynamic model is exact, and state one remedy.

Solution:

When \(u_c\) exceeds the position limit, the aircraft receives \(u_a\neq u_c\). The reference model still predicts the response associated with the unattainable command, so tracking error appears. The update law treats this error as evidence of uncertainty and changes \(\widehat{\mathbf{W}}\), even though the plant parameters are correct. A pseudo-control-hedged reference model adds \(\mathbf{B}(u_a-u_c)\) to the reference dynamics, or adaptation can be frozen or reduced while saturation persists.

Adaptive flight control is most effective as an augmentation to a stable baseline architecture. Aerodynamic and inertial changes are represented through a regressor and ideal uncertainty parameters, while a reference model defines the desired closed-loop behavior. Under exact matching and constant parameters, a Lyapunov update law guarantees bounded signals and asymptotic tracking. Real aircraft add slowly varying parameters, unmatched dynamics, actuator limits, delay, noise, and flexible modes; projection, normalization, leakage, hedging, protection logic, and systematic verification are therefore essential. The supplied multi-language implementations demonstrate recovery from two abrupt matched aerodynamic changes using one common mathematical model.