An autopilot attempts to make measured variables follow commanded variables. In a longitudinal channel these variables may represent a normalized pitch response, angle-of-attack response, normal-acceleration response, or a carefully selected combination of states. Adaptation can recover tracking performance when aerodynamic derivatives, mass properties, actuator effectiveness, or environmental conditions differ from their nominal values. Adaptation alone, however, does not imply constraint satisfaction.

The conceptual architecture therefore assigns two distinct responsibilities:

1. The adaptive autopilot stabilizes the uncertain plant and reduces the tracking error relative to a reference model.
2. The envelope-protection layer replaces an inadmissible raw command by the closest command predicted to preserve a tightened safe set.

flowchart TD
  R["Raw pilot or guidance command"] --> G["Envelope protection and command governor"]
  X["Measured normalized flight variables"] --> G
  G --> RP["Protected command"]
  RP --> M["Reference model"]
  RP --> C["Baseline plus adaptive autopilot"]
  M --> E["Tracking error"]
  P["Uncertain normalized plant"] --> X
  C --> P
  X --> E
  E --> A["Parameter update with projection and leakage"]
  A --> C
        

The separation is important. A fast parameter estimator is not a safety certificate, and a safety governor is not a replacement for stabilization. The protection calculation must be based on explicit constraints and on a defensible bound for the difference between the actual response and the protected reference-model response.

To expose the mathematics without introducing aircraft-specific data, consider a normalized scalar channel \( y(t) \). Its uncertain dynamics are written in a form already familiar from linearly parameterized adaptive control:

\[ \dot y = -a_0 y + b u + \boldsymbol{\theta}^{*\top} \boldsymbol{\phi}(t,y) + d(t), \qquad b > 0. \]

Here \( a_0 \) is a nominal damping coefficient, \( b \) is the known-sign control effectiveness, \( \boldsymbol{\theta}^{*} \) is an unknown constant or slowly varying matched-uncertainty vector, \( \boldsymbol{\phi}(t,y) \) is a known regressor, and \( d(t) \) collects bounded unmatched effects and numerical imperfections.

The protected command \( r_p \) drives a stable first-order reference model:

\[ \dot y_m = -a_m y_m + b_m r_p, \qquad a_m > 0, \qquad b_m > 0. \]

A direct adaptive control law is

\[ u = \frac{1}{b}\left[(a_0-a_m)y+b_m r_p- \hat{\boldsymbol{\theta} }^{\top}\boldsymbol{\phi}(t,y)\right]. \]

Define tracking and parameter errors by \( e=y-y_m \) and \( \tilde{\boldsymbol{\theta} }= \boldsymbol{\theta}^{*}-\hat{\boldsymbol{\theta} } \). Substitution gives the error dynamics

\[ \dot e=-a_m e+\tilde{\boldsymbol{\theta} }^{\top} \boldsymbol{\phi}(t,y)+d(t). \]

This equation is the bridge between adaptation and protection. The governor constrains \( y_m \); the adaptive law must make \( e \) sufficiently small that the actual output \( y=y_m+e \) remains safe.

First set \( d(t)=0 \) and use the ideal gradient update

\[ \dot{\hat{\boldsymbol{\theta} } }=\mathbf{\Gamma} \boldsymbol{\phi}(t,y)e, \qquad \mathbf{\Gamma}=\mathbf{\Gamma}^{\top} > 0. \]

Consider the composite Lyapunov function

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

Since \( \dot{\tilde{\boldsymbol{\theta} } }=-\mathbf{\Gamma} \boldsymbol{\phi}e \), differentiation yields

\[ \begin{aligned} \dot V &=e\left(-a_m e+\tilde{\boldsymbol{\theta} }^{\top} \boldsymbol{\phi}\right) +\tilde{\boldsymbol{\theta} }^{\top}\mathbf{\Gamma}^{-1} \left(-\mathbf{\Gamma}\boldsymbol{\phi}e\right)\\ &=-a_m e^2. \end{aligned} \]

Therefore all signals represented in \( V \) are bounded. Under the standard bounded-regressor assumptions, Barbalat's lemma gives \( e(t)\rightarrow 0 \). As emphasized in Chapter 10, this does not by itself imply convergence of \( \hat{\boldsymbol{\theta} } \) to the physical parameter vector; persistent excitation is required for that stronger result.

Flight-envelope protection cannot rely on the disturbance-free asymptotic statement alone. Measurement noise, unmodeled dynamics, actuator saturation, and sampling errors motivate projection and leakage:

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

Assume the true parameter belongs to the projection set \( \Theta \), and temporarily analyze the interior where projection is inactive. The Lyapunov derivative becomes

\[ \dot V=-a_m e^2+e d+ \sigma\tilde{\boldsymbol{\theta} }^{\top} \mathbf{\Gamma}^{-1}\hat{\boldsymbol{\theta} }. \]

Using \( \hat{\boldsymbol{\theta} }=\boldsymbol{\theta}^{*}- \tilde{\boldsymbol{\theta} } \), Young's inequality, and \( |d(t)|\leq d_{\max} \), one obtains

\[ \dot V\leq -\frac{a_m}{2}e^2- \frac{\sigma}{2}\left\|\tilde{\boldsymbol{\theta} }\right\|_{ \mathbf{\Gamma}^{-1} }^2+c, \]

\[ c=\frac{d_{\max}^2}{2a_m}+ \frac{\sigma}{2}\left\|\boldsymbol{\theta}^{*}\right\|_{ \mathbf{\Gamma}^{-1} }^2. \]

With \( \eta=\min(a_m,\sigma) \), the comparison inequality

\[ \dot V\leq -\eta V+c \]

gives

\[ V(t)\leq e^{-\eta t}V(0)+\frac{c}{\eta} \left(1-e^{-\eta t}\right). \]

Consequently, an asymptotic error estimate is

\[ \limsup_{t\rightarrow\infty}|e(t)| \leq \sqrt{\frac{2c}{\eta} }. \]

In an engineering assurance process, the adopted error budget \( \varepsilon_e \) must be no smaller than the largest bound supported by analysis, Monte Carlo testing, numerical-error analysis, sensor tolerances, actuator limits, and model-validation evidence. It is not legitimate to select \( \varepsilon_e \) only because a nominal simulation happens to track well.

Let the normalized protected output satisfy the hard envelope

\[ \mathcal{Y}=\left\{y\in\mathbb{R}: |y|\leq y_{\max}\right\}. \]

Because \( y=y_m+e \), the reference model must use the tightened set

\[ \mathcal{Y}_m=\left\{y_m\in\mathbb{R}: |y_m|\leq y_{\mathrm{tight} }\right\}, \qquad y_{\mathrm{tight} }=y_{\max}-\varepsilon_e. \]

A meaningful tightened set requires \( 0<\varepsilon_e<y_{\max} \). The difference \( y_{\max}-y_{\mathrm{tight} } \) is not unused performance margin; it is explicitly reserved for adaptive tracking error and other quantified discrepancies.

In a multivariable design, the envelope may be represented by inequalities such as \( \mathbf{H}x\leq\mathbf{h} \), ellipsoids, lookup-table boundaries, or reachable-set approximations. The scalar construction is used here because it admits a transparent invariance proof and exact equivalent implementations in five programming languages.

Sample the reference model with period \( T_s \). Its exact zero-order-hold representation is

\[ y_{m,k+1}=F y_{m,k}+G r_{p,k}, \qquad F=e^{-a_m T_s}, \qquad G=\frac{b_m}{a_m}(1-F). \]

Since \( a_m>0 \), \( b_m>0 \), and \( T_s>0 \), one has \( 0<F<1 \) and \( G>0 \). Requiring the next reference-model output to remain in the tightened set gives

\[ -y_{\mathrm{tight} }\leq F y_{m,k}+G r_{p,k} \leq y_{\mathrm{tight} }. \]

Solving for the admissible protected-command interval yields

\[ \frac{-y_{\mathrm{tight} }-F y_{m,k} }{G} \leq r_{p,k}\leq \frac{y_{\mathrm{tight} }-F y_{m,k} }{G}. \]

The closest admissible scalar command is therefore

\[ r_{p,k}=\operatorname{clip}\!\left( r_k, \frac{-y_{\mathrm{tight} }-F y_{m,k} }{G}, \frac{y_{\mathrm{tight} }-F y_{m,k} }{G} \right). \]

flowchart TD
  A["Read raw command and current model output"] --> B["Compute tightened output limit"]
  B --> C["Compute exact sampled model coefficients F and G"]
  C --> D["Compute lower and upper admissible commands"]
  D --> E["Clip raw command to admissible interval"]
  E --> F["Advance protected reference model"]
  F --> G["Adaptive autopilot tracks protected model"]
  G --> H["Check error budget, actuator limits, and monitors"]
        

For vector commands, the clipping operation becomes a small constrained optimization problem, commonly a quadratic program that minimizes command modification subject to predicted state, output, and actuator constraints. The scalar governor here is the closed-form special case.

Theorem. Suppose:

1. \( |y_{m,0}|\leq y_{\mathrm{tight} } \);
2. the reference governor enforces \( |y_{m,k}|\leq y_{\mathrm{tight} } \) for every sample;
3. the adaptive closed loop satisfies the verified bound \( |e_k|\leq\varepsilon_e \) for every sample; and
4. \( y_{\mathrm{tight} }=y_{\max}-\varepsilon_e \).

Then \( |y_k|\leq y_{\max} \) for every sample.

Proof of reference-model invariance. The governor selects

\[ r_{p,k}\in\left[ \frac{-y_{\mathrm{tight} }-F y_{m,k} }{G}, \frac{y_{\mathrm{tight} }-F y_{m,k} }{G} \right]. \]

Multiplication by \( G>0 \) and addition of \( F y_{m,k} \) give

\[ -y_{\mathrm{tight} }\leq F y_{m,k}+G r_{p,k} \leq y_{\mathrm{tight} }. \]

Since the center expression equals \( y_{m,k+1} \), the tightened set is positively invariant under the sampled reference-model dynamics.

Proof of actual-output safety. By the triangle inequality,

\[ |y_k|=|y_{m,k}+e_k| \leq |y_{m,k}|+|e_k| \leq y_{\mathrm{tight} }+\varepsilon_e =y_{\max}. \]

The theorem is conditional, not absolute. Any failure of the assumed error budget, sampling model, actuator authority, state measurement, or constraint model invalidates the conclusion. A practical system therefore needs monitors, fallback modes, sensor validation, saturation logic, and independent verification in addition to the mathematical governor.

The downloadable programs implement the same normalized experiment. The raw command intentionally requests values outside the feasible steady region. The reference governor tightens the normalized hard envelope from \( y_{\max}=1 \) to \( y_{\mathrm{tight} }=0.75 \) using the design budget \( \varepsilon_e=0.25 \). The uncertain plant contains two matched basis functions and a small sinusoidal disturbance.

The controller uses projection, sigma modification, actuator saturation, and adaptation freezing during saturation. The latter is an implementation safeguard rather than part of the ideal Lyapunov proof. The provided test run gives a maximum absolute plant output of approximately \( 0.931 \) and a maximum absolute tracking error of approximately \( 0.181 \). These numbers demonstrate the example but do not establish a general certified error bound.

The important observation is structural: the raw command can be much larger than the envelope, yet the protected reference remains invariant and the adaptive loop follows it with bounded error. Increasing the adaptation gain may reduce some low-frequency error but can amplify noise and sampling sensitivity; it cannot replace constraint tightening.

The Python implementation uses NumPy for vector calculations, Matplotlib for visualization, and the standard csv module for reproducible data export.

Chapter26_Lesson2.py

"""Chapter 26, Lesson 2: adaptive autopilot with conceptual envelope protection.

Educational normalized simulation only. It is not flight-certified software and
must not be used to control a real vehicle.
"""
from __future__ import annotations

import csv
from dataclasses import dataclass
from math import exp, sin
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 0.01
    duration: float = 30.0
    a_nominal: float = 1.0
    b_known: float = 1.0
    a_model: float = 1.4
    b_model: float = 1.4
    gamma: float = 4.0
    sigma: float = 0.03
    theta_limit: float = 2.5
    output_limit: float = 1.0
    tracking_error_budget: float = 0.25
    actuator_limit: float = 3.0


def raw_command(t: float) -> float:
    if t < 1.0:
        return 0.0
    if t < 10.0:
        return 1.40
    if t < 18.0:
        return -1.30
    if t < 25.0:
        return 0.90
    return 0.0


def project(vector: np.ndarray, limit: float) -> np.ndarray:
    norm = float(np.linalg.norm(vector))
    return vector if norm <= limit else vector * (limit / norm)


def protected_command(
    command: float, model_output: float, cfg: Config
) -> tuple[float, float, float]:
    """One-step reference governor for the exact sampled reference model."""
    tightened_limit = cfg.output_limit - cfg.tracking_error_budget
    f = exp(-cfg.a_model * cfg.dt)
    g = (cfg.b_model / cfg.a_model) * (1.0 - f)
    lower = (-tightened_limit - f * model_output) / g
    upper = (tightened_limit - f * model_output) / g
    return float(np.clip(command, lower, upper)), lower, upper


def simulate(cfg: Config) -> dict[str, np.ndarray]:
    steps = int(round(cfg.duration / cfg.dt)) + 1
    time = np.linspace(0.0, cfg.duration, steps)

    y = np.zeros(steps)
    y_model = np.zeros(steps)
    theta_hat = np.zeros((steps, 2))
    command = np.zeros(steps)
    command_protected = np.zeros(steps)
    control = np.zeros(steps)

    theta_true = np.array([0.55, 0.18])
    f = exp(-cfg.a_model * cfg.dt)
    g = (cfg.b_model / cfg.a_model) * (1.0 - f)

    for k in range(steps - 1):
        t = time[k]
        command[k] = raw_command(t)
        rp, _, _ = protected_command(command[k], y_model[k], cfg)
        command_protected[k] = rp

        phi = np.array([y[k], sin(0.5 * t)])
        u_unsat = (
            (cfg.a_nominal - cfg.a_model) * y[k]
            + cfg.b_model * rp
            - float(theta_hat[k] @ phi)
        ) / cfg.b_known
        control[k] = float(np.clip(u_unsat, -cfg.actuator_limit, cfg.actuator_limit))

        disturbance = 0.04 * sin(1.7 * t)
        y_dot = (
            -cfg.a_nominal * y[k]
            + cfg.b_known * control[k]
            + float(theta_true @ phi)
            + disturbance
        )
        y[k + 1] = y[k] + cfg.dt * y_dot
        y_model[k + 1] = f * y_model[k] + g * rp

        error = y[k] - y_model[k]
        adaptation_enabled = abs(u_unsat) <= cfg.actuator_limit
        theta_dot = (
            cfg.gamma * phi * error - cfg.sigma * theta_hat[k]
            if adaptation_enabled
            else -cfg.sigma * theta_hat[k]
        )
        theta_hat[k + 1] = project(
            theta_hat[k] + cfg.dt * theta_dot, cfg.theta_limit
        )

    command[-1] = raw_command(time[-1])
    command_protected[-1] = protected_command(command[-1], y_model[-1], cfg)[0]
    control[-1] = control[-2]

    return {
        "time": time,
        "raw_command": command,
        "protected_command": command_protected,
        "model_output": y_model,
        "plant_output": y,
        "tracking_error": y - y_model,
        "control": control,
        "theta_hat_1": theta_hat[:, 0],
        "theta_hat_2": theta_hat[:, 1],
    }


def save_results(data: dict[str, np.ndarray], directory: Path, cfg: Config) -> None:
    directory.mkdir(parents=True, exist_ok=True)
    csv_path = directory / "Chapter26_Lesson2_results.csv"
    names = list(data.keys())
    with csv_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(names)
        writer.writerows(zip(*(data[name] for name in names)))

    plt.figure(figsize=(10, 6))
    plt.plot(data["time"], data["raw_command"], "--", label="raw command")
    plt.plot(data["time"], data["protected_command"], label="protected command")
    plt.plot(data["time"], data["model_output"], label="reference-model output")
    plt.plot(data["time"], data["plant_output"], label="plant output")
    plt.axhline(cfg.output_limit, linestyle=":", label="envelope")
    plt.axhline(-cfg.output_limit, linestyle=":")
    plt.xlabel("time")
    plt.ylabel("normalized response")
    plt.title("Adaptive autopilot with tightened reference governor")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(directory / "Chapter26_Lesson2_response.png", dpi=160)
    plt.close()


if __name__ == "__main__":
    configuration = Config()
    results = simulate(configuration)
    save_results(results, Path.cwd(), configuration)
    max_output = float(np.max(np.abs(results["plant_output"])))
    max_error = float(np.max(np.abs(results["tracking_error"])))
    print(f"Maximum absolute plant output: {max_output:.4f}")
    print(f"Maximum absolute tracking error: {max_error:.4f}")
    print(f"Envelope respected in this run: {max_output <= configuration.output_limit}")

Both implementations depend only on their standard libraries and write CSV files. This keeps the adaptive law and governor transparent and makes them convenient for unit testing.

Chapter26_Lesson2.cpp

// Chapter 26, Lesson 2: conceptual adaptive autopilot and envelope protection.
// Educational normalized simulation only; not flight-certified software.
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>

struct Config {
    double dt = 0.01;
    double duration = 30.0;
    double aNominal = 1.0;
    double bKnown = 1.0;
    double aModel = 1.4;
    double bModel = 1.4;
    double gamma = 4.0;
    double sigma = 0.03;
    double thetaLimit = 2.5;
    double outputLimit = 1.0;
    double trackingErrorBudget = 0.25;
    double actuatorLimit = 3.0;
};

double rawCommand(double t) {
    if (t < 1.0) return 0.0;
    if (t < 10.0) return 1.40;
    if (t < 18.0) return -1.30;
    if (t < 25.0) return 0.90;
    return 0.0;
}

double clip(double value, double lower, double upper) {
    return std::max(lower, std::min(value, upper));
}

std::array<double, 2> project(std::array<double, 2> value, double limit) {
    const double norm = std::hypot(value[0], value[1]);
    if (norm > limit && norm > 0.0) {
        value[0] *= limit / norm;
        value[1] *= limit / norm;
    }
    return value;
}

double protectedCommand(double command, double ym, const Config& cfg) {
    const double tightened = cfg.outputLimit - cfg.trackingErrorBudget;
    const double f = std::exp(-cfg.aModel * cfg.dt);
    const double g = (cfg.bModel / cfg.aModel) * (1.0 - f);
    const double lower = (-tightened - f * ym) / g;
    const double upper = ( tightened - f * ym) / g;
    return clip(command, lower, upper);
}

int main() {
    try {
        const Config cfg;
        const int steps = static_cast<int>(std::llround(cfg.duration / cfg.dt)) + 1;
        const std::array<double, 2> thetaTrue{0.55, 0.18};
        std::array<double, 2> thetaHat{0.0, 0.0};
        double y = 0.0;
        double ym = 0.0;
        double maxAbsY = 0.0;
        double maxAbsError = 0.0;
        const double f = std::exp(-cfg.aModel * cfg.dt);
        const double g = (cfg.bModel / cfg.aModel) * (1.0 - f);

        std::ofstream csv("Chapter26_Lesson2_results_cpp.csv");
        if (!csv) throw std::runtime_error("Cannot open output CSV file.");
        csv << "time,raw_command,protected_command,model_output,plant_output,tracking_error,control,theta_hat_1,theta_hat_2\n";
        csv << std::setprecision(12);

        for (int k = 0; k < steps; ++k) {
            const double t = k * cfg.dt;
            const double command = rawCommand(t);
            const double rp = protectedCommand(command, ym, cfg);
            const std::array<double, 2> phi{y, std::sin(0.5 * t)};
            const double estimate = thetaHat[0] * phi[0] + thetaHat[1] * phi[1];
            const double uUnsat = ((cfg.aNominal - cfg.aModel) * y + cfg.bModel * rp - estimate) / cfg.bKnown;
            const double u = clip(uUnsat, -cfg.actuatorLimit, cfg.actuatorLimit);
            const double error = y - ym;

            csv << t << ',' << command << ',' << rp << ',' << ym << ',' << y << ','
                << error << ',' << u << ',' << thetaHat[0] << ',' << thetaHat[1] << '\n';
            maxAbsY = std::max(maxAbsY, std::abs(y));
            maxAbsError = std::max(maxAbsError, std::abs(error));
            if (k == steps - 1) break;

            const double disturbance = 0.04 * std::sin(1.7 * t);
            const double yDot = -cfg.aNominal * y + cfg.bKnown * u
                + thetaTrue[0] * phi[0] + thetaTrue[1] * phi[1] + disturbance;
            y += cfg.dt * yDot;
            ym = f * ym + g * rp;

            const bool adaptationEnabled = std::abs(uUnsat) <= cfg.actuatorLimit;
            std::array<double, 2> thetaDot{
                -cfg.sigma * thetaHat[0],
                -cfg.sigma * thetaHat[1]
            };
            if (adaptationEnabled) {
                thetaDot[0] += cfg.gamma * phi[0] * error;
                thetaDot[1] += cfg.gamma * phi[1] * error;
            }
            thetaHat[0] += cfg.dt * thetaDot[0];
            thetaHat[1] += cfg.dt * thetaDot[1];
            thetaHat = project(thetaHat, cfg.thetaLimit);
        }

        std::cout << "Maximum absolute plant output: " << maxAbsY << '\n';
        std::cout << "Maximum absolute tracking error: " << maxAbsError << '\n';
        std::cout << "Envelope respected in this run: "
                  << (maxAbsY <= cfg.outputLimit ? "true" : "false") << '\n';
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter26_Lesson2.java

// Chapter 26, Lesson 2: conceptual adaptive autopilot and envelope protection.
// Educational normalized simulation only; not flight-certified software.
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;

public final class Chapter26_Lesson2 {
    private static final double DT = 0.01;
    private static final double DURATION = 30.0;
    private static final double A_NOMINAL = 1.0;
    private static final double B_KNOWN = 1.0;
    private static final double A_MODEL = 1.4;
    private static final double B_MODEL = 1.4;
    private static final double GAMMA = 4.0;
    private static final double SIGMA = 0.03;
    private static final double THETA_LIMIT = 2.5;
    private static final double OUTPUT_LIMIT = 1.0;
    private static final double ERROR_BUDGET = 0.25;
    private static final double ACTUATOR_LIMIT = 3.0;

    private Chapter26_Lesson2() {}

    private static double rawCommand(double t) {
        if (t < 1.0) return 0.0;
        if (t < 10.0) return 1.40;
        if (t < 18.0) return -1.30;
        if (t < 25.0) return 0.90;
        return 0.0;
    }

    private static double clip(double value, double lower, double upper) {
        return Math.max(lower, Math.min(value, upper));
    }

    private static double protectedCommand(double command, double ym) {
        double tightened = OUTPUT_LIMIT - ERROR_BUDGET;
        double f = Math.exp(-A_MODEL * DT);
        double g = (B_MODEL / A_MODEL) * (1.0 - f);
        double lower = (-tightened - f * ym) / g;
        double upper = ( tightened - f * ym) / g;
        return clip(command, lower, upper);
    }

    private static void project(double[] theta) {
        double norm = Math.hypot(theta[0], theta[1]);
        if (norm > THETA_LIMIT && norm > 0.0) {
            theta[0] *= THETA_LIMIT / norm;
            theta[1] *= THETA_LIMIT / norm;
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        int steps = (int) Math.round(DURATION / DT) + 1;
        double[] thetaTrue = {0.55, 0.18};
        double[] thetaHat = {0.0, 0.0};
        double y = 0.0;
        double ym = 0.0;
        double maxAbsY = 0.0;
        double maxAbsError = 0.0;
        double f = Math.exp(-A_MODEL * DT);
        double g = (B_MODEL / A_MODEL) * (1.0 - f);

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter26_Lesson2_results_java.csv"))) {
            writer.write("time,raw_command,protected_command,model_output,plant_output,tracking_error,control,theta_hat_1,theta_hat_2\n");

            for (int k = 0; k < steps; k++) {
                double t = k * DT;
                double command = rawCommand(t);
                double rp = protectedCommand(command, ym);
                double[] phi = {y, Math.sin(0.5 * t)};
                double estimate = thetaHat[0] * phi[0] + thetaHat[1] * phi[1];
                double uUnsat = ((A_NOMINAL - A_MODEL) * y + B_MODEL * rp - estimate) / B_KNOWN;
                double u = clip(uUnsat, -ACTUATOR_LIMIT, ACTUATOR_LIMIT);
                double error = y - ym;

                writer.write(String.format(Locale.US,
                    "%.8f,%.8f,%.8f,%.8f,%.8f,%.8f,%.8f,%.8f,%.8f%n",
                    t, command, rp, ym, y, error, u, thetaHat[0], thetaHat[1]));
                maxAbsY = Math.max(maxAbsY, Math.abs(y));
                maxAbsError = Math.max(maxAbsError, Math.abs(error));
                if (k == steps - 1) break;

                double disturbance = 0.04 * Math.sin(1.7 * t);
                double yDot = -A_NOMINAL * y + B_KNOWN * u
                    + thetaTrue[0] * phi[0] + thetaTrue[1] * phi[1] + disturbance;
                y += DT * yDot;
                ym = f * ym + g * rp;

                boolean adaptationEnabled = Math.abs(uUnsat) <= ACTUATOR_LIMIT;
                double thetaDot0 = -SIGMA * thetaHat[0];
                double thetaDot1 = -SIGMA * thetaHat[1];
                if (adaptationEnabled) {
                    thetaDot0 += GAMMA * phi[0] * error;
                    thetaDot1 += GAMMA * phi[1] * error;
                }
                thetaHat[0] += DT * thetaDot0;
                thetaHat[1] += DT * thetaDot1;
                project(thetaHat);
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }

        System.out.printf(Locale.US, "Maximum absolute plant output: %.4f%n", maxAbsY);
        System.out.printf(Locale.US, "Maximum absolute tracking error: %.4f%n", maxAbsError);
        System.out.println("Envelope respected in this run: " + (maxAbsY <= OUTPUT_LIMIT));
    }
}

The MATLAB script exports a table and figure. A Simulink realization can map the plant and reference model to Integrator blocks, the controller and command governor to MATLAB Function blocks, the parameter vector to an Integrator with projection logic, and actuator authority to a Saturation block. In a serious model, discrete rates, data types, reset behavior, and saturation flags must be explicit rather than inherited from solver defaults.

Chapter26_Lesson2.m

% Chapter 26, Lesson 2: conceptual adaptive autopilot and envelope protection.
% Educational normalized simulation only; not flight-certified software.
clear; clc; close all;

dt = 0.01;
T = 30.0;
t = (0:dt:T)';
N = numel(t);

aNominal = 1.0;
bKnown = 1.0;
aModel = 1.4;
bModel = 1.4;
gamma = 4.0;
sigma = 0.03;
thetaLimit = 2.5;
outputLimit = 1.0;
errorBudget = 0.25;
actuatorLimit = 3.0;
thetaTrue = [0.55; 0.18];

y = zeros(N,1);
ym = zeros(N,1);
thetaHat = zeros(N,2);
raw = zeros(N,1);
protected = zeros(N,1);
u = zeros(N,1);

F = exp(-aModel*dt);
G = (bModel/aModel)*(1-F);
tightened = outputLimit-errorBudget;

for k = 1:N-1
    raw(k) = commandProfile(t(k));
    lower = (-tightened-F*ym(k))/G;
    upper = ( tightened-F*ym(k))/G;
    protected(k) = min(max(raw(k),lower),upper);

    phi = [y(k); sin(0.5*t(k))];
    uUnsat = ((aNominal-aModel)*y(k)+bModel*protected(k)-thetaHat(k,:)*phi)/bKnown;
    u(k) = min(max(uUnsat,-actuatorLimit),actuatorLimit);

    disturbance = 0.04*sin(1.7*t(k));
    yDot = -aNominal*y(k)+bKnown*u(k)+thetaTrue.'*phi+disturbance;
    y(k+1) = y(k)+dt*yDot;
    ym(k+1) = F*ym(k)+G*protected(k);

    error = y(k)-ym(k);
    if abs(uUnsat) <= actuatorLimit
        thetaDot = gamma*phi*error-sigma*thetaHat(k,:).';
    else
        thetaDot = -sigma*thetaHat(k,:).';
    end
    candidate = thetaHat(k,:).'+dt*thetaDot;
    if norm(candidate) > thetaLimit
        candidate = thetaLimit*candidate/norm(candidate);
    end
    thetaHat(k+1,:) = candidate.';
end

raw(N) = commandProfile(t(N));
lower = (-tightened-F*ym(N))/G;
upper = ( tightened-F*ym(N))/G;
protected(N) = min(max(raw(N),lower),upper);
u(N) = u(N-1);
error = y-ym;

results = table(t,raw,protected,ym,y,error,u,thetaHat(:,1),thetaHat(:,2), ...
    'VariableNames',{'time','raw_command','protected_command','model_output', ...
    'plant_output','tracking_error','control','theta_hat_1','theta_hat_2'});
writetable(results,'Chapter26_Lesson2_results_matlab.csv');

figure;
plot(t,raw,'--',t,protected,t,ym,t,y,'LineWidth',1.1); hold on;
yline(outputLimit,':'); yline(-outputLimit,':');
grid on; xlabel('time'); ylabel('normalized response');
title('Adaptive autopilot with tightened reference governor');
legend('raw command','protected command','reference-model output', ...
    'plant output','envelope','Location','best');
exportgraphics(gcf,'Chapter26_Lesson2_response_matlab.png','Resolution',160);

fprintf('Maximum absolute plant output: %.4f\n',max(abs(y)));
fprintf('Maximum absolute tracking error: %.4f\n',max(abs(error)));
fprintf('Envelope respected in this run: %d\n',max(abs(y)) <= outputLimit);

function r = commandProfile(time)
if time < 1.0
    r = 0.0;
elseif time < 10.0
    r = 1.40;
elseif time < 18.0
    r = -1.30;
elseif time < 25.0
    r = 0.90;
else
    r = 0.0;
end
end

The Mathematica file is a plain-text Wolfram Language notebook script. It can be evaluated with Get["Chapter26_Lesson2.nb"] and exports matching CSV and PNG results.

Chapter26_Lesson2.nb

(* Chapter 26, Lesson 2: conceptual adaptive autopilot and envelope protection. *)
(* Educational normalized simulation only; not flight-certified software. *)
ClearAll["Global`*"];

dt = 0.01; duration = 30.;
aNominal = 1.; bKnown = 1.; aModel = 1.4; bModel = 1.4;
gamma = 4.; sigma = 0.03; thetaLimit = 2.5;
outputLimit = 1.; errorBudget = 0.25; actuatorLimit = 3.;
thetaTrue = {0.55, 0.18};
f = Exp[-aModel dt]; g = (bModel/aModel) (1 - f);
tightened = outputLimit - errorBudget;

rawCommand[t_] := Piecewise[{ {0., t < 1.}, {1.40, t < 10.},
   {-1.30, t < 18.}, {0.90, t < 25.} }, 0.];
clip[x_, lo_, hi_] := Min[Max[x, lo], hi];
project[v_] := If[Norm[v] <= thetaLimit, v, thetaLimit v/Norm[v]];
protectedCommand[r_, ym_] := Module[{lo, hi},
  lo = (-tightened - f ym)/g;
  hi = ( tightened - f ym)/g;
  clip[r, lo, hi]
];

steps = Round[duration/dt] + 1;
y = 0.; ym = 0.; thetaHat = {0., 0.};
records = Reap[
   Do[
    t = (k - 1) dt;
    r = rawCommand[t]; rp = protectedCommand[r, ym];
    phi = {y, Sin[0.5 t]};
    uUnsat = ((aNominal - aModel) y + bModel rp - thetaHat.phi)/bKnown;
    u = clip[uUnsat, -actuatorLimit, actuatorLimit];
    e = y - ym;
    Sow[{t, r, rp, ym, y, e, u, thetaHat[[1]], thetaHat[[2]]}];
    If[k < steps,
      disturbance = 0.04 Sin[1.7 t];
      yDot = -aNominal y + bKnown u + thetaTrue.phi + disturbance;
      y = y + dt yDot;
      ym = f ym + g rp;
      thetaDot = If[Abs[uUnsat] <= actuatorLimit,
        gamma phi e - sigma thetaHat, -sigma thetaHat];
      thetaHat = project[thetaHat + dt thetaDot];
    ],
    {k, 1, steps}
   ]
][[2, 1]];

header = {"time", "raw_command", "protected_command", "model_output",
   "plant_output", "tracking_error", "control", "theta_hat_1", "theta_hat_2"};
Export["Chapter26_Lesson2_results_mathematica.csv", Prepend[records, header]];

responsePlot = ListLinePlot[
  {records[[All, {1, 2}]], records[[All, {1, 3}]],
   records[[All, {1, 4}]], records[[All, {1, 5}]]},
  PlotLegends -> {"raw command", "protected command", "reference-model output", "plant output"},
  Frame -> True, FrameLabel -> {"time", "normalized response"},
  PlotRange -> All, GridLines -> {None, {-outputLimit, outputLimit} },
  PlotLabel -> "Adaptive autopilot with tightened reference governor"
];
Export["Chapter26_Lesson2_response_mathematica.png", responsePlot, ImageResolution -> 160];

maxOutput = Max[Abs[records[[All, 5]]]];
maxError = Max[Abs[records[[All, 6]]]];
Print["Maximum absolute plant output: ", NumberForm[maxOutput, {5, 4}]];
Print["Maximum absolute tracking error: ", NumberForm[maxError, {5, 4}]];
Print["Envelope respected in this run: ", maxOutput <= outputLimit];

The conceptual theorem must be embedded in a much broader assurance process. At minimum, a real design distinguishes command limits, state limits, structural limits, aerodynamic limits, actuator limits, sensor-validity regions, and model-validity regions. A single scalar limiter cannot represent these coupled constraints.

Separation of nominal and adaptive authority. The baseline controller should provide acceptable stability and handling over a declared nominal region. Adaptive augmentation is then bounded so that an estimator fault, poor excitation, or transient parameter excursion cannot command arbitrary control effort.

Constraint hierarchy. Hard constraints require conservative margins and deterministic fallback behavior. Soft performance constraints may be relaxed with penalties. Confusing the two categories can produce smooth simulations but unsafe edge cases.

Sampling and delay. The proof uses an exact sampled reference model with a known period. Jitter, computation delay, sensor transport delay, and asynchronous updates change the prediction. The admissible interval must be recomputed for the implemented timing model.

Saturation interaction. When the actuator saturates, the matched uncertainty cancellation assumed in the ideal error dynamics is no longer exact. Common safeguards include command prefilters, pseudo-control hedging, adaptation freezing, anti-windup states, and explicit actuator models. Each safeguard changes the proof obligations.

Validation evidence. Useful evidence includes linear and nonlinear analysis, reachability calculations, Monte Carlo campaigns, worst-case optimization, processor-in-the-loop testing, hardware-in-the-loop testing, fault injection, independent code review, and traceability from every envelope boundary to a verified requirement.

The downloadable code is intentionally excluded from any claim of airworthiness, certification compliance, fault tolerance, or real-time determinism.

Problem 1 (Lyapunov cancellation): For the disturbance-free error dynamics \( \dot e=-a_m e+\tilde{\boldsymbol{\theta} }^{\top} \boldsymbol{\phi} \), derive an adaptive law that cancels the mixed term in \( \dot V \) for the Lyapunov function in Section 3.

Solution: Differentiation gives

\[ \dot V=-a_m e^2+e\tilde{\boldsymbol{\theta} }^{\top} \boldsymbol{\phi}+ ilde{\boldsymbol{\theta} }^{\top} \mathbf{\Gamma}^{-1}\dot{\tilde{\boldsymbol{\theta} } }. \]

Since \( \dot{\tilde{\boldsymbol{\theta} } }=- \dot{\hat{\boldsymbol{\theta} } } \), choose \( \dot{\hat{\boldsymbol{\theta} } }= \mathbf{\Gamma}\boldsymbol{\phi}e \). Then the two mixed terms are equal and opposite, leaving \( \dot V=-a_m e^2\leq 0 \).

Problem 2 (Admissible command interval): Let \( a_m=1.4 \), \( b_m=1.4 \), \( T_s=0.01 \), \( y_{\mathrm{tight} }=0.75 \), and \( y_{m,k}=0.74 \). Compute the admissible interval for \( r_{p,k} \).

Solution:

\[ F=e^{-1.4(0.01)}\approx 0.98610, \qquad G=1-F\approx 0.01390. \]

\[ r_{p,k}^{\min}=\frac{-0.75-(0.98610)(0.74)}{0.01390} \approx -106.44, \]

\[ r_{p,k}^{\max}=\frac{0.75-(0.98610)(0.74)}{0.01390} \approx 1.46. \]

Thus a large negative command remains one-step feasible because it moves the model inward, whereas the positive command is tightly limited near the upper boundary. Additional actuator and rate constraints would narrow this interval.

Problem 3 (Constraint tightening): A verified campaign gives \( |e(t)|\leq 0.08 \) and the hard normalized limit is \( y_{\max}=1.2 \). Determine the largest symmetric reference-model envelope that supports the triangle-inequality guarantee.

Solution:

\[ y_{\mathrm{tight} }=y_{\max}-\varepsilon_e=1.2-0.08=1.12. \]

Therefore the governor must enforce \( |y_m|\leq 1.12 \). If the reference model reaches \( 1.12 \) and the error reaches its worst positive value, the actual output reaches but does not exceed \( 1.20 \).

Problem 4 (Why simulation is not a proof): A nominal simulation produces \( \max |e|=0.03 \). Explain why choosing \( \varepsilon_e=0.03 \) may be invalid.

Solution: The observed maximum applies only to the tested initial state, parameter realization, disturbance, command sequence, solver, and sampling pattern. A defensible budget must also cover the declared ranges of parameter uncertainty, disturbances, sensor errors, delays, saturation, numerical error, mode transitions, and faults. The budget must be supported by analysis or systematic verification over those ranges, normally with an additional assurance margin.

Problem 5 (Vector generalization): Let the sampled protected model be \( x_{m,k+1}=\mathbf{F}x_{m,k}+\mathbf{G}v_k \) and the tightened polyhedral constraint be \( \mathbf{H}x_{m,k+1}\leq\mathbf{h}_{\mathrm{tight} } \). Formulate the closest-command governor for raw command \( r_k \).

Solution: With positive-definite weighting matrix \( \mathbf{W} \), solve

\[ \begin{aligned} v_k^{*}=\arg\min_{v}\quad & \frac{1}{2}(v-r_k)^{\top}\mathbf{W}(v-r_k)\\ \text{subject to}\quad & \mathbf{H}(\mathbf{F}x_{m,k}+\mathbf{G}v) \leq\mathbf{h}_{\mathrm{tight} }. \end{aligned} \]

This is a convex quadratic program when the constraints are linear and \( \mathbf{W} \) is positive definite. Recursive feasibility still requires an invariant admissible set or a suitable terminal condition; one-step feasibility alone is not automatically sufficient for every multivariable system.

An adaptive autopilot and a flight-envelope protection system solve different problems. The adaptive controller reduces the uncertain plant's tracking error, while the command governor constrains the protected reference. A scalar exact sampled reference model gives a closed-form admissible command interval. If the tightened reference set is invariant and the adaptive error remains below a verified budget, the triangle inequality guarantees that the actual normalized output remains inside the original envelope. Projection, leakage, saturation handling, timing analysis, fallback logic, and systematic validation are necessary to translate this conceptual result into a credible safety case.