After completing this lesson, students should be able to:

• define an operating condition as a point on an equilibrium or trajectory family;
• select scheduling variables that are measurable, dynamically relevant, bounded, and sufficiently regular;
• derive frozen local models and identify the terms lost when the operating point moves;
• construct an admissible scheduling set containing both value and rate constraints;
• distinguish scheduling information from online parameter identification used in MRAC or self-tuning regulators.

A gain-scheduled controller changes its coefficients according to a measured or estimated signal that indicates the current operating condition. The scheduling signal is not normally an online estimate of an unknown constant parameter. Instead, it is a known or observable variable such as speed, altitude, payload, temperature, pressure, configuration, or commanded equilibrium.

\[ \rho(t) \in \mathcal{P} \subset \mathbb{R}^{n_\rho}, \qquad K(t)=K\!\left(\rho(t)\right). \]

The notation \( \rho \) will denote the scheduling vector, while \( \mathcal{P} \) denotes its admissible operating envelope.

Consider a nonlinear plant with state \( x\in\mathbb{R}^{n} \), control input \( u\in\mathbb{R}^{m} \), measured output \( y\in\mathbb{R}^{p} \), and an exogenous operating variable \( w\in\mathbb{R}^{q} \):

\[ \dot{x}=f(x,u,w), \qquad y=h(x,u,w). \]

A steady operating condition is a triple \( (x_e,u_e,w_e) \) satisfying

\[ f(x_e,u_e,w_e)=0. \]

A collection of such points forms an equilibrium set. For scheduling, it is convenient to parameterize the relevant portion of that set by \( \rho \):

\[ x_e=X_e(\rho),\qquad u_e=U_e(\rho),\qquad w_e=W_e(\rho), \qquad \rho\in\mathcal{P}. \]

The functions \( X_e,U_e,W_e \) define an equilibrium manifold. They must satisfy the identity

\[ f\!\left(X_e(\rho),U_e(\rho),W_e(\rho)\right)=0 \quad \text{for every } \rho\in\mathcal{P}. \]

Differentiating this identity with respect to the scalar scheduling variable gives the equilibrium sensitivity equation

\[ A_e(\rho)X_e'(\rho) +B_e(\rho)U_e'(\rho) +E_e(\rho)W_e'(\rho)=0, \]

where

\[ A_e(\rho)= \left.\frac{\partial f}{\partial x}\right|_e,\qquad B_e(\rho)= \left.\frac{\partial f}{\partial u}\right|_e,\qquad E_e(\rho)= \left.\frac{\partial f}{\partial w}\right|_e. \]

This equation is useful when constructing trim tables. It relates the change in equilibrium state and input to a change in operating condition.

A scheduling variable should reveal how the plant dynamics change. It may originate outside the plant, from the plant state, or from a command generator.

flowchart TD
  A["Candidate operating information"] --> B["Exogenous variable"]
  A --> C["Endogenous measured \nvariable"]
  A --> D["Estimated or \nreconstructed variable"]
  A --> E["Commanded operating point"]
  B --> B1["Examples: altitude, \nambient temperature, \npayload"]
  C --> C1["Examples: speed, \npressure, joint angle"]
  D --> D1["Examples: filtered speed, \nestimated mass class"]
  E --> E1["Examples: commanded velocity \nor trim condition"]
  B1 --> F["Check relevance, range, rate, \nand measurement quality"]
  C1 --> F
  D1 --> F
  E1 --> F
        

3.1 Exogenous scheduling variables

An exogenous variable is generated outside the controlled state dynamics or is treated as externally specified. Examples include ambient density, known payload configuration, or a commanded operating point. It may vary with time, but its value is not determined by the controller state equation being scheduled.

3.2 Endogenous scheduling variables

An endogenous variable is computed from the plant state or output:

\[ \rho=\psi(x,y,u). \]

Because the controller affects \( x,y,u \), an endogenous schedule creates an additional feedback path. This does not make the method invalid, but it means that \( \rho(t) \) cannot be treated as an arbitrary independent signal during analysis.

3.3 Measured versus estimated scheduling variables

Direct measurement is preferable when the sensor bandwidth and accuracy are adequate. When reconstruction is required, write

\[ \hat{\rho}(t)=\rho(t)+e_\rho(t). \]

The scheduling error \( e_\rho \) changes the controller coefficients. A first-order filter,

\[ \tau_f\dot{\hat{\rho}}+\hat{\rho}=\rho_m, \]

reduces high-frequency noise but introduces lag. If \( |\dot{\rho}(t)|\leq\nu \) and measurement noise is neglected, the asymptotic filtering error satisfies the conservative bound

\[ \limsup_{t\to\infty} |\hat{\rho}(t)-\rho(t)| \leq \tau_f\nu. \]

Therefore, filtering and schedule-rate assumptions must be designed together.

Define perturbation coordinates relative to the scheduled equilibrium:

\[ \delta x=x-X_e(\rho),\qquad \delta u=u-U_e(\rho),\qquad \delta w=w-W_e(\rho). \]

Differentiating the state perturbation gives

\[ \delta\dot{x} =\dot{x}-X_\rho(\rho)\dot{\rho}, \qquad X_\rho(\rho)=\frac{\partial X_e}{\partial\rho}. \]

A first-order Taylor expansion about the equilibrium indexed by \( \rho \) yields

\[ \delta\dot{x} = A(\rho)\delta x +B(\rho)\delta u +E(\rho)\delta w -X_\rho(\rho)\dot{\rho} +r_2, \]

with

\[ A(\rho)= \left.\frac{\partial f}{\partial x}\right|_e,\qquad B(\rho)= \left.\frac{\partial f}{\partial u}\right|_e,\qquad E(\rho)= \left.\frac{\partial f}{\partial w}\right|_e. \]

The remainder \( r_2 \) contains second- and higher-order perturbation terms. If the Jacobian of \( f \) is locally Lipschitz with constant \( L \), then in a sufficiently small neighborhood,

\[ \|r_2\| \leq \frac{L}{2} \left\| \begin{bmatrix} \delta x\\ \delta u\\ \delta w \end{bmatrix} \right\|^2. \]

4.1 Why the coupling term matters

The term \( -X_\rho(\rho)\dot{\rho} \) is created solely because the coordinate origin is moving. A frozen linearization assumes \( \dot{\rho}=0 \) and therefore omits this term. Slow scheduling makes the omission smaller, but does not make it identically zero.

4.2 Vector scheduling variables

For \( \rho\in\mathbb{R}^{n_\rho} \),

\[ X_\rho(\rho)\dot{\rho} = \sum_{i=1}^{n_\rho} \frac{\partial X_e}{\partial\rho_i}\dot{\rho}_i. \]

Each scheduling channel can therefore contribute a different moving-equilibrium forcing term.

At a fixed value \( \rho=\rho_0 \), the local frozen model is

\[ \delta\dot{x} = A(\rho_0)\delta x+B(\rho_0)\delta u, \qquad \delta y= C(\rho_0)\delta x+D(\rho_0)\delta u. \]

Replacing \( \rho_0 \) by the current scheduling signal produces the parameter-varying representation

\[ \delta\dot{x} = A(\rho(t))\delta x+B(\rho(t))\delta u, \qquad \delta y= C(\rho(t))\delta x+D(\rho(t))\delta u. \]

This model is linear in the state and input for a prescribed scheduling trajectory, but its matrices vary with time through \( \rho(t) \). It is therefore linear parameter-varying rather than linear time-invariant.

5.1 Affine parameter dependence

A particularly useful form is affine dependence:

\[ A(\rho)=A_0+\sum_{i=1}^{n_\rho}\rho_iA_i, \qquad B(\rho)=B_0+\sum_{i=1}^{n_\rho}\rho_iB_i. \]

Not every nonlinear plant naturally produces this form. A nonlinear dependence can be retained directly, approximated over a grid, or embedded using additional scheduling coordinates. The choice must not hide algebraic dependence or create scheduling variables that cannot be measured.

5.2 Frozen behavior is local information

Eigenvalues of \( A(\rho_0) \) describe the local response when the operating point is held fixed. They do not by themselves establish global or time-varying closed-loop stability. Detailed stability analysis is developed later in this chapter.

A candidate scheduling variable should be evaluated against six engineering criteria.

1. Dynamic relevance: changes in the variable should explain meaningful changes in plant or controller coefficients.
2. Availability: the variable must be measured, commanded, or estimated with known latency and uncertainty.
3. Single-valuedness: the same scheduling value should not correspond to drastically different dynamics unless additional scheduling coordinates are included.
4. Regularity: the scheduling map and coefficient maps should be continuous or piecewise smooth over their intended domain.
5. Envelope coverage: the schedule must cover all operating conditions expected in service, including margins.
6. Rate compatibility: the expected rate of variation must be compatible with the assumptions used in analysis and implementation.

6.1 A normalized relevance measure

For a scalar candidate \( \rho_i \), one simple local sensitivity score is

\[ S_i(\rho)= \frac{\Delta\rho_i} {1+\|A(\rho)\|_F+\|B(\rho)\|_F} \left( \left\|\frac{\partial A}{\partial\rho_i}\right\|_F + \left\|\frac{\partial B}{\partial\rho_i}\right\|_F \right). \]

A very small score indicates that the candidate contributes little to explaining local dynamic variation. A large score indicates relevance, but does not guarantee that the variable is measurable or that the schedule is well conditioned.

6.2 Avoid redundant scheduling coordinates

If two variables are nearly functionally dependent throughout the operating envelope, using both can create a sparse or ill-conditioned grid. Redundancy can be assessed from sampled data using the covariance matrix

\[ \Sigma_\rho= \frac{1}{N-1} \sum_{k=1}^{N} (\rho_k-\bar{\rho})(\rho_k-\bar{\rho})^\mathsf{T}. \]

Very small eigenvalues of \( \Sigma_\rho \) indicate that the observed schedule occupies a lower-dimensional subset of the proposed scheduling space.

A useful operating envelope includes bounds on both the scheduling values and their rates:

\[ \mathcal{P}= \left\{ \rho\in\mathbb{R}^{n_\rho}: \underline{\rho}_i \leq \rho_i \leq \overline{\rho}_i \right\}, \]

\[ \mathcal{V}= \left\{ \dot{\rho}\in\mathbb{R}^{n_\rho}: |\dot{\rho}_i| \leq \nu_i \right\}. \]

The implemented scheduler should continuously check whether \( (\rho,\dot{\rho})\in\mathcal{P}\times\mathcal{V} \). Outside this set, the controller is operating beyond its validated model and scheduling assumptions.

7.1 Normalization

Mapping each schedule to approximately \( [-1,1] \) improves numerical scaling:

\[ \bar{\rho}_i = 2 \frac{\rho_i-\underline{\rho}_i} {\overline{\rho}_i-\underline{\rho}_i} -1. \]

7.2 A bound associated with slow scheduling

Consider the scalar moving-equilibrium local model

\[ \dot{\delta}=a(\rho)\delta-\dot{\rho}, \qquad a(\rho)\leq-\alpha<0, \qquad |\dot{\rho}|\leq\nu. \]

Variation of constants gives

\[ \delta(t) = \Phi(t,0)\delta(0) - \int_0^t \Phi(t,\tau)\dot{\rho}(\tau)\,d\tau, \]

and the transition factor satisfies \( |\Phi(t,\tau)|\leq e^{-\alpha(t-\tau)} \). Therefore,

\[ |\delta(t)| \leq e^{-\alpha t}|\delta(0)| + \frac{\nu}{\alpha} \left(1-e^{-\alpha t}\right). \]

The ratio \( \nu/\alpha \) quantifies the tracking displacement caused by a moving operating point. Slow scheduling means that this ratio is small relative to the permitted perturbation.

Consider the nonlinear scalar plant

\[ \dot{x}=-a x-bx^3+g u, \qquad a>0,\quad b>0,\quad g\neq0. \]

Choose the equilibrium state itself as the operating coordinate: \( x_e(\rho)=\rho \). The equilibrium input follows from \( \dot{x}=0 \):

\[ u_e(\rho)=\frac{a\rho+b\rho^3}{g}. \]

The frozen local state coefficient is

\[ A(\rho) = \left. \frac{\partial}{\partial x} (-ax-bx^3+gu) \right|_{x=\rho} = -a-3b\rho^2, \]

while \( B(\rho)=g \). The coefficient sensitivity is

\[ \frac{dA}{d\rho}=-6b\rho. \]

Hence the local dynamics change weakly near \( \rho=0 \) and more strongly at larger \( |\rho| \).

8.1 Exact perturbation dynamics

Let \( x=\rho+\delta \) and apply the equilibrium feedforward input \( u=u_e(\rho) \). Direct expansion gives

\[ \dot{\delta} = (-a-3b\rho^2)\delta -3b\rho\delta^2 -b\delta^3 -\dot{\rho}. \]

The rate-aware local approximation retains the first-order terms:

\[ \dot{\delta}_{ra} = A(\rho)\delta_{ra}-\dot{\rho}. \]

The frozen approximation is

\[ \dot{\delta}_{fr} = A(\rho)\delta_{fr}. \]

The implementations below compare these two approximations against the nonlinear plant. With the selected test trajectory, the rate-aware model has a much smaller perturbation error because it includes the moving-origin term.

flowchart TD
  A["Define required operating envelope"] --> B["Compute trim or equilibrium points"]
  B --> C["List measurable candidate scheduling variables"]
  C --> D["Evaluate dynamic relevance and single-valuedness"]
  D --> E["Specify value, rate, noise, and latency limits"]
  E --> F["Linearize or identify local models"]
  F --> G["Validate local models between grid points"]
  G --> H["Create runtime envelope checks and fallback behavior"]
  H --> I["Proceed to controller scheduling design"]
        

The runtime implementation should record the schedule value, estimated rate, envelope status, selected model region, and any fallback action. This diagnostic information is essential when the controller is tested across a wide operating range.

All five implementations use the same nonlinear plant, RK4 integration, equilibrium map, scheduling trajectory, and diagnostics. They generate an operating-point table and compare the nonlinear perturbation with rate-aware and frozen local models.

The Python version uses NumPy and Matplotlib. The C++ version uses the C++17 standard library, and the Java version uses the Java standard library. The MATLAB version uses built-in array, table, plotting, and CSV functions; the same equations can be placed in Simulink using MATLAB Function, Integrator, and Lookup Table blocks. The Wolfram Mathematica notebook uses native list processing, graphics, and export functions.

Chapter20_Lesson1.py

# Chapter20_Lesson1.py
# Scheduling variables and operating conditions for a nonlinear plant.
#
# Plant:
#   x_dot = -a*x - b*x^3 + g*u
#
# Equilibrium family parameterized by rho:
#   x_e(rho) = rho
#   u_e(rho) = (a*rho + b*rho^3)/g
#
# Perturbation delta = x - rho obeys:
#   delta_dot = A(rho)*delta - 3*b*rho*delta^2 - b*delta^3 - rho_dot
#   A(rho) = -a - 3*b*rho^2
#
# The term -rho_dot is the moving-operating-point coupling term. A frozen
# linearization omits it; the rate-aware local model retains it.

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class PlantParameters:
    a: float = 1.2
    b: float = 0.35
    g: float = 1.0
    rho_min: float = -2.0
    rho_max: float = 2.0
    rho_rate_max: float = 0.8


def equilibrium_input(rho: np.ndarray | float, p: PlantParameters) -> np.ndarray | float:
    return (p.a * rho + p.b * np.asarray(rho) ** 3) / p.g


def local_a(rho: np.ndarray | float, p: PlantParameters) -> np.ndarray | float:
    return -p.a - 3.0 * p.b * np.asarray(rho) ** 2


def normalized_schedule(rho: np.ndarray, p: PlantParameters) -> np.ndarray:
    return 2.0 * (rho - p.rho_min) / (p.rho_max - p.rho_min) - 1.0


def scheduling_signal(t: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    rho = 1.40 * np.sin(0.25 * t) + 0.35 * np.sin(1.10 * t)
    rho_dot = 0.35 * np.cos(0.25 * t) + 0.385 * np.cos(1.10 * t)
    return rho, rho_dot


def rk4_step(fun, t: float, x: float, dt: float) -> float:
    k1 = fun(t, x)
    k2 = fun(t + 0.5 * dt, x + 0.5 * dt * k1)
    k3 = fun(t + 0.5 * dt, x + 0.5 * dt * k2)
    k4 = fun(t + dt, x + dt * k3)
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(p: PlantParameters, dt: float = 0.01, tf: float = 30.0):
    t = np.arange(0.0, tf + 0.5 * dt, dt)
    rho, rho_dot = scheduling_signal(t)
    u_e = equilibrium_input(rho, p)

    x_nonlinear = np.empty_like(t)
    delta_rate_aware = np.empty_like(t)
    delta_frozen = np.empty_like(t)

    x_nonlinear[0] = rho[0] + 0.40
    delta_rate_aware[0] = x_nonlinear[0] - rho[0]
    delta_frozen[0] = delta_rate_aware[0]

    def rho_at(time: float) -> float:
        return 1.40 * np.sin(0.25 * time) + 0.35 * np.sin(1.10 * time)

    def rho_dot_at(time: float) -> float:
        return 0.35 * np.cos(0.25 * time) + 0.385 * np.cos(1.10 * time)

    for k in range(t.size - 1):
        tk = t[k]

        def nonlinear_rhs(time: float, x: float) -> float:
            r = rho_at(time)
            ue = equilibrium_input(r, p)
            return -p.a * x - p.b * x**3 + p.g * ue

        def rate_aware_rhs(time: float, delta: float) -> float:
            r = rho_at(time)
            return float(local_a(r, p)) * delta - rho_dot_at(time)

        def frozen_rhs(time: float, delta: float) -> float:
            r = rho_at(time)
            return float(local_a(r, p)) * delta

        x_nonlinear[k + 1] = rk4_step(nonlinear_rhs, tk, x_nonlinear[k], dt)
        delta_rate_aware[k + 1] = rk4_step(
            rate_aware_rhs, tk, delta_rate_aware[k], dt
        )
        delta_frozen[k + 1] = rk4_step(frozen_rhs, tk, delta_frozen[k], dt)

    delta_exact = x_nonlinear - rho
    return t, rho, rho_dot, u_e, x_nonlinear, delta_exact, delta_rate_aware, delta_frozen


def write_operating_table(p: PlantParameters, path: Path) -> None:
    rho_grid = np.linspace(p.rho_min, p.rho_max, 17)
    a_grid = local_a(rho_grid, p)
    tau_grid = -1.0 / a_grid
    u_grid = equilibrium_input(rho_grid, p)
    rho_n = normalized_schedule(rho_grid, p)

    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(
            ["rho", "rho_normalized", "equilibrium_input", "local_A", "time_constant"]
        )
        writer.writerows(zip(rho_grid, rho_n, u_grid, a_grid, tau_grid))


def main() -> None:
    p = PlantParameters()
    t, rho, rho_dot, u_e, x, delta_exact, delta_rate_aware, delta_frozen = simulate(p)

    inside_envelope = np.logical_and(rho >= p.rho_min, rho <= p.rho_max)
    inside_rate = np.abs(rho_dot) <= p.rho_rate_max

    rmse_rate_aware = float(np.sqrt(np.mean((delta_exact - delta_rate_aware) ** 2)))
    rmse_frozen = float(np.sqrt(np.mean((delta_exact - delta_frozen) ** 2)))

    print("Scheduling-variable diagnostics")
    print(f"  rho range: [{rho.min():.4f}, {rho.max():.4f}]")
    print(f"  max |rho_dot|: {np.max(np.abs(rho_dot)):.4f}")
    print(f"  samples inside value envelope: {100.0 * inside_envelope.mean():.2f}%")
    print(f"  samples inside rate envelope:  {100.0 * inside_rate.mean():.2f}%")
    print(f"  RMSE, rate-aware local model:  {rmse_rate_aware:.6f}")
    print(f"  RMSE, frozen local model:      {rmse_frozen:.6f}")

    output_dir = Path(__file__).resolve().parent
    write_operating_table(p, output_dir / "Chapter20_Lesson1_operating_table.csv")

    fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
    axes[0].plot(t, rho, label="rho(t)")
    axes[0].axhline(p.rho_min, linestyle="--", label="value envelope")
    axes[0].axhline(p.rho_max, linestyle="--")
    axes[0].set_ylabel("schedule")
    axes[0].legend()

    axes[1].plot(t, rho_dot, label="rho_dot(t)")
    axes[1].axhline(-p.rho_rate_max, linestyle="--", label="rate envelope")
    axes[1].axhline(p.rho_rate_max, linestyle="--")
    axes[1].set_ylabel("schedule rate")
    axes[1].legend()

    axes[2].plot(t, delta_exact, label="nonlinear perturbation")
    axes[2].plot(t, delta_rate_aware, linestyle="--", label="rate-aware local model")
    axes[2].plot(t, delta_frozen, linestyle=":", label="frozen local model")
    axes[2].set_xlabel("time [s]")
    axes[2].set_ylabel("delta")
    axes[2].legend()

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


if __name__ == "__main__":
    main()

Chapter20_Lesson1.cpp

// Chapter20_Lesson1.cpp
// Scheduling variables and operating conditions for a nonlinear plant.
//
// Plant:
//   x_dot = -a*x - b*x^3 + g*u
//
// Equilibrium family:
//   x_e(rho) = rho
//   u_e(rho) = (a*rho + b*rho^3)/g
//
// Perturbation delta = x - rho:
//   delta_dot = A(rho)*delta - 3*b*rho*delta^2 - b*delta^3 - rho_dot
//   A(rho) = -a - 3*b*rho^2
//
// The program compares a rate-aware local model with the frozen local model.

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

struct PlantParameters {
    double a = 1.2;
    double b = 0.35;
    double g = 1.0;
    double rho_min = -2.0;
    double rho_max = 2.0;
    double rho_rate_max = 0.8;
};

double equilibrium_input(double rho, const PlantParameters& p) {
    return (p.a * rho + p.b * rho * rho * rho) / p.g;
}

double local_a(double rho, const PlantParameters& p) {
    return -p.a - 3.0 * p.b * rho * rho;
}

double normalized_schedule(double rho, const PlantParameters& p) {
    return 2.0 * (rho - p.rho_min) / (p.rho_max - p.rho_min) - 1.0;
}

double schedule(double t) {
    return 1.40 * std::sin(0.25 * t) + 0.35 * std::sin(1.10 * t);
}

double schedule_rate(double t) {
    return 0.35 * std::cos(0.25 * t) + 0.385 * std::cos(1.10 * t);
}

double rk4_step(
    const std::function<double(double, double)>& rhs,
    double t,
    double x,
    double dt
) {
    const double k1 = rhs(t, x);
    const double k2 = rhs(t + 0.5 * dt, x + 0.5 * dt * k1);
    const double k3 = rhs(t + 0.5 * dt, x + 0.5 * dt * k2);
    const double k4 = rhs(t + dt, x + dt * k3);
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

void write_operating_table(const PlantParameters& p, const std::string& filename) {
    std::ofstream out(filename);
    if (!out) {
        throw std::runtime_error("Cannot open " + filename);
    }

    out << "rho,rho_normalized,equilibrium_input,local_A,time_constant\n";
    constexpr int points = 17;
    for (int i = 0; i < points; ++i) {
        const double alpha = static_cast<double>(i) / static_cast<double>(points - 1);
        const double rho = p.rho_min + alpha * (p.rho_max - p.rho_min);
        const double A = local_a(rho, p);
        out << std::setprecision(12)
            << rho << ','
            << normalized_schedule(rho, p) << ','
            << equilibrium_input(rho, p) << ','
            << A << ','
            << -1.0 / A << '\n';
    }
}

int main() {
    try {
        const PlantParameters p;
        const double dt = 0.01;
        const double tf = 30.0;
        const std::size_t n = static_cast<std::size_t>(std::llround(tf / dt)) + 1U;

        std::vector<double> t(n), rho(n), rho_dot(n), x(n);
        std::vector<double> delta_exact(n), delta_rate_aware(n), delta_frozen(n);

        for (std::size_t k = 0; k < n; ++k) {
            t[k] = static_cast<double>(k) * dt;
            rho[k] = schedule(t[k]);
            rho_dot[k] = schedule_rate(t[k]);
        }

        x[0] = rho[0] + 0.40;
        delta_rate_aware[0] = x[0] - rho[0];
        delta_frozen[0] = delta_rate_aware[0];

        for (std::size_t k = 0; k + 1U < n; ++k) {
            const double tk = t[k];

            const auto nonlinear_rhs = [&p](double time, double state) {
                const double r = schedule(time);
                const double ue = equilibrium_input(r, p);
                return -p.a * state - p.b * state * state * state + p.g * ue;
            };

            const auto rate_aware_rhs = [&p](double time, double delta) {
                const double r = schedule(time);
                return local_a(r, p) * delta - schedule_rate(time);
            };

            const auto frozen_rhs = [&p](double time, double delta) {
                return local_a(schedule(time), p) * delta;
            };

            x[k + 1U] = rk4_step(nonlinear_rhs, tk, x[k], dt);
            delta_rate_aware[k + 1U] =
                rk4_step(rate_aware_rhs, tk, delta_rate_aware[k], dt);
            delta_frozen[k + 1U] =
                rk4_step(frozen_rhs, tk, delta_frozen[k], dt);
        }

        double squared_error_rate_aware = 0.0;
        double squared_error_frozen = 0.0;
        std::size_t inside_value = 0U;
        std::size_t inside_rate = 0U;
        double rho_min_seen = rho.front();
        double rho_max_seen = rho.front();
        double max_abs_rate = 0.0;

        std::ofstream series("Chapter20_Lesson1_results_cpp.csv");
        if (!series) {
            throw std::runtime_error("Cannot create results CSV");
        }
        series << "time,rho,rho_dot,x,delta_exact,delta_rate_aware,delta_frozen\n";

        for (std::size_t k = 0; k < n; ++k) {
            delta_exact[k] = x[k] - rho[k];
            const double e_rate = delta_exact[k] - delta_rate_aware[k];
            const double e_frozen = delta_exact[k] - delta_frozen[k];
            squared_error_rate_aware += e_rate * e_rate;
            squared_error_frozen += e_frozen * e_frozen;

            if (rho[k] >= p.rho_min && rho[k] <= p.rho_max) {
                ++inside_value;
            }
            if (std::abs(rho_dot[k]) <= p.rho_rate_max) {
                ++inside_rate;
            }

            rho_min_seen = std::min(rho_min_seen, rho[k]);
            rho_max_seen = std::max(rho_max_seen, rho[k]);
            max_abs_rate = std::max(max_abs_rate, std::abs(rho_dot[k]));

            series << std::setprecision(12)
                   << t[k] << ',' << rho[k] << ',' << rho_dot[k] << ','
                   << x[k] << ',' << delta_exact[k] << ','
                   << delta_rate_aware[k] << ',' << delta_frozen[k] << '\n';
        }

        const double rmse_rate = std::sqrt(squared_error_rate_aware / n);
        const double rmse_frozen = std::sqrt(squared_error_frozen / n);

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Scheduling-variable diagnostics\n";
        std::cout << "  rho range: [" << rho_min_seen << ", " << rho_max_seen << "]\n";
        std::cout << "  max |rho_dot|: " << max_abs_rate << '\n';
        std::cout << "  samples inside value envelope: "
                  << 100.0 * static_cast<double>(inside_value) / n << "%\n";
        std::cout << "  samples inside rate envelope:  "
                  << 100.0 * static_cast<double>(inside_rate) / n << "%\n";
        std::cout << "  RMSE, rate-aware local model:  " << rmse_rate << '\n';
        std::cout << "  RMSE, frozen local model:      " << rmse_frozen << '\n';

        write_operating_table(p, "Chapter20_Lesson1_operating_table_cpp.csv");
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter20_Lesson1.java

// Chapter20_Lesson1.java
// Scheduling variables and operating conditions for a nonlinear plant.
//
// Plant:
//   x_dot = -a*x - b*x^3 + g*u
//
// Equilibrium family:
//   x_e(rho) = rho
//   u_e(rho) = (a*rho + b*rho^3)/g
//
// The program compares a moving-operating-point local model with a frozen
// local model that omits the coupling term -rho_dot.

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;
import java.util.function.DoubleBinaryOperator;

public final class Chapter20_Lesson1 {
    private Chapter20_Lesson1() {
    }

    private static final class PlantParameters {
        final double a = 1.2;
        final double b = 0.35;
        final double g = 1.0;
        final double rhoMin = -2.0;
        final double rhoMax = 2.0;
        final double rhoRateMax = 0.8;
    }

    private static double equilibriumInput(double rho, PlantParameters p) {
        return (p.a * rho + p.b * rho * rho * rho) / p.g;
    }

    private static double localA(double rho, PlantParameters p) {
        return -p.a - 3.0 * p.b * rho * rho;
    }

    private static double normalizedSchedule(double rho, PlantParameters p) {
        return 2.0 * (rho - p.rhoMin) / (p.rhoMax - p.rhoMin) - 1.0;
    }

    private static double schedule(double t) {
        return 1.40 * Math.sin(0.25 * t) + 0.35 * Math.sin(1.10 * t);
    }

    private static double scheduleRate(double t) {
        return 0.35 * Math.cos(0.25 * t) + 0.385 * Math.cos(1.10 * t);
    }

    private static double rk4Step(
            DoubleBinaryOperator rhs,
            double t,
            double x,
            double dt
    ) {
        double k1 = rhs.applyAsDouble(t, x);
        double k2 = rhs.applyAsDouble(t + 0.5 * dt, x + 0.5 * dt * k1);
        double k3 = rhs.applyAsDouble(t + 0.5 * dt, x + 0.5 * dt * k2);
        double k4 = rhs.applyAsDouble(t + dt, x + dt * k3);
        return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    private static void writeOperatingTable(PlantParameters p) throws IOException {
        Path path = Path.of("Chapter20_Lesson1_operating_table_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                path, StandardCharsets.UTF_8)) {
            writer.write(
                    "rho,rho_normalized,equilibrium_input,local_A,time_constant\n");
            final int points = 17;
            for (int i = 0; i < points; ++i) {
                double alpha = (double) i / (double) (points - 1);
                double rho = p.rhoMin + alpha * (p.rhoMax - p.rhoMin);
                double aLocal = localA(rho, p);
                writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        rho,
                        normalizedSchedule(rho, p),
                        equilibriumInput(rho, p),
                        aLocal,
                        -1.0 / aLocal
                ));
            }
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        PlantParameters p = new PlantParameters();
        double dt = 0.01;
        double tf = 30.0;
        int n = (int) Math.round(tf / dt) + 1;

        double[] t = new double[n];
        double[] rho = new double[n];
        double[] rhoDot = new double[n];
        double[] x = new double[n];
        double[] deltaExact = new double[n];
        double[] deltaRateAware = new double[n];
        double[] deltaFrozen = new double[n];

        for (int k = 0; k < n; ++k) {
            t[k] = k * dt;
            rho[k] = schedule(t[k]);
            rhoDot[k] = scheduleRate(t[k]);
        }

        x[0] = rho[0] + 0.40;
        deltaRateAware[0] = x[0] - rho[0];
        deltaFrozen[0] = deltaRateAware[0];

        DoubleBinaryOperator nonlinearRhs = (time, state) -> {
            double r = schedule(time);
            double ue = equilibriumInput(r, p);
            return -p.a * state - p.b * state * state * state + p.g * ue;
        };

        DoubleBinaryOperator rateAwareRhs = (time, delta) -> {
            double r = schedule(time);
            return localA(r, p) * delta - scheduleRate(time);
        };

        DoubleBinaryOperator frozenRhs = (time, delta) ->
                localA(schedule(time), p) * delta;

        for (int k = 0; k < n - 1; ++k) {
            x[k + 1] = rk4Step(nonlinearRhs, t[k], x[k], dt);
            deltaRateAware[k + 1] =
                    rk4Step(rateAwareRhs, t[k], deltaRateAware[k], dt);
            deltaFrozen[k + 1] =
                    rk4Step(frozenRhs, t[k], deltaFrozen[k], dt);
        }

        double squaredErrorRateAware = 0.0;
        double squaredErrorFrozen = 0.0;
        int insideValue = 0;
        int insideRate = 0;
        double rhoMinSeen = rho[0];
        double rhoMaxSeen = rho[0];
        double maxAbsRate = 0.0;

        try (BufferedWriter writer = Files.newBufferedWriter(
                Path.of("Chapter20_Lesson1_results_java.csv"),
                StandardCharsets.UTF_8)) {
            writer.write(
                    "time,rho,rho_dot,x,delta_exact,delta_rate_aware,delta_frozen\n");

            for (int k = 0; k < n; ++k) {
                deltaExact[k] = x[k] - rho[k];
                double eRate = deltaExact[k] - deltaRateAware[k];
                double eFrozen = deltaExact[k] - deltaFrozen[k];
                squaredErrorRateAware += eRate * eRate;
                squaredErrorFrozen += eFrozen * eFrozen;

                if (rho[k] >= p.rhoMin && rho[k] <= p.rhoMax) {
                    ++insideValue;
                }
                if (Math.abs(rhoDot[k]) <= p.rhoRateMax) {
                    ++insideRate;
                }

                rhoMinSeen = Math.min(rhoMinSeen, rho[k]);
                rhoMaxSeen = Math.max(rhoMaxSeen, rho[k]);
                maxAbsRate = Math.max(maxAbsRate, Math.abs(rhoDot[k]));

                writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        t[k],
                        rho[k],
                        rhoDot[k],
                        x[k],
                        deltaExact[k],
                        deltaRateAware[k],
                        deltaFrozen[k]
                ));
            }

            double rmseRate = Math.sqrt(squaredErrorRateAware / n);
            double rmseFrozen = Math.sqrt(squaredErrorFrozen / n);

            System.out.println("Scheduling-variable diagnostics");
            System.out.printf("  rho range: [%.6f, %.6f]%n", rhoMinSeen, rhoMaxSeen);
            System.out.printf("  max |rho_dot|: %.6f%n", maxAbsRate);
            System.out.printf(
                    "  samples inside value envelope: %.2f%%%n",
                    100.0 * insideValue / n
            );
            System.out.printf(
                    "  samples inside rate envelope:  %.2f%%%n",
                    100.0 * insideRate / n
            );
            System.out.printf(
                    "  RMSE, rate-aware local model:  %.6f%n",
                    rmseRate
            );
            System.out.printf(
                    "  RMSE, frozen local model:      %.6f%n",
                    rmseFrozen
            );

            writeOperatingTable(p);
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

Chapter20_Lesson1.m

% Chapter20_Lesson1.m
% Scheduling variables and operating conditions for a nonlinear plant.
%
% Plant:
%   x_dot = -a*x - b*x^3 + g*u
%
% Equilibrium family:
%   x_e(rho) = rho
%   u_e(rho) = (a*rho + b*rho^3)/g
%
% Perturbation delta = x - rho:
%   delta_dot = A(rho)*delta - 3*b*rho*delta^2 - b*delta^3 - rho_dot
%   A(rho) = -a - 3*b*rho^2
%
% The script compares a rate-aware local model with a frozen local model.

clear;
clc;

p.a = 1.2;
p.b = 0.35;
p.g = 1.0;
p.rhoMin = -2.0;
p.rhoMax = 2.0;
p.rhoRateMax = 0.8;

dt = 0.01;
tf = 30.0;
t = (0:dt:tf).';

rho = 1.40*sin(0.25*t) + 0.35*sin(1.10*t);
rhoDot = 0.35*cos(0.25*t) + 0.385*cos(1.10*t);
uEq = equilibriumInput(rho, p);

x = zeros(size(t));
deltaRateAware = zeros(size(t));
deltaFrozen = zeros(size(t));

x(1) = rho(1) + 0.40;
deltaRateAware(1) = x(1) - rho(1);
deltaFrozen(1) = deltaRateAware(1);

for k = 1:numel(t)-1
    tk = t(k);

    nonlinearRhs = @(time, state) ...
        -p.a*state - p.b*state.^3 + ...
        p.g*equilibriumInput(schedule(time), p);

    rateAwareRhs = @(time, delta) ...
        localA(schedule(time), p).*delta - scheduleRate(time);

    frozenRhs = @(time, delta) ...
        localA(schedule(time), p).*delta;

    x(k+1) = rk4Step(nonlinearRhs, tk, x(k), dt);
    deltaRateAware(k+1) = ...
        rk4Step(rateAwareRhs, tk, deltaRateAware(k), dt);
    deltaFrozen(k+1) = ...
        rk4Step(frozenRhs, tk, deltaFrozen(k), dt);
end

deltaExact = x - rho;

insideValue = rho >= p.rhoMin & rho <= p.rhoMax;
insideRate = abs(rhoDot) <= p.rhoRateMax;
rmseRateAware = sqrt(mean((deltaExact - deltaRateAware).^2));
rmseFrozen = sqrt(mean((deltaExact - deltaFrozen).^2));

fprintf('Scheduling-variable diagnostics\n');
fprintf('  rho range: [%.6f, %.6f]\n', min(rho), max(rho));
fprintf('  max |rho_dot|: %.6f\n', max(abs(rhoDot)));
fprintf('  samples inside value envelope: %.2f%%\n', ...
    100*mean(insideValue));
fprintf('  samples inside rate envelope:  %.2f%%\n', ...
    100*mean(insideRate));
fprintf('  RMSE, rate-aware local model:  %.6f\n', rmseRateAware);
fprintf('  RMSE, frozen local model:      %.6f\n', rmseFrozen);

rhoGrid = linspace(p.rhoMin, p.rhoMax, 17).';
AGrid = localA(rhoGrid, p);
tauGrid = -1./AGrid;
uGrid = equilibriumInput(rhoGrid, p);
rhoNormalized = normalizedSchedule(rhoGrid, p);

operatingTable = table( ...
    rhoGrid, rhoNormalized, uGrid, AGrid, tauGrid, ...
    'VariableNames', { ...
        'rho', ...
        'rho_normalized', ...
        'equilibrium_input', ...
        'local_A', ...
        'time_constant' ...
    } ...
);
writetable(operatingTable, 'Chapter20_Lesson1_operating_table_matlab.csv');

resultsTable = table( ...
    t, rho, rhoDot, uEq, x, deltaExact, deltaRateAware, deltaFrozen, ...
    'VariableNames', { ...
        'time', ...
        'rho', ...
        'rho_dot', ...
        'equilibrium_input', ...
        'x', ...
        'delta_exact', ...
        'delta_rate_aware', ...
        'delta_frozen' ...
    } ...
);
writetable(resultsTable, 'Chapter20_Lesson1_results_matlab.csv');

figure('Color', 'w');

subplot(3,1,1);
plot(t, rho, 'LineWidth', 1.2);
hold on;
yline(p.rhoMin, '--');
yline(p.rhoMax, '--');
grid on;
ylabel('\rho');
title('Scheduling signal and admissible value envelope');

subplot(3,1,2);
plot(t, rhoDot, 'LineWidth', 1.2);
hold on;
yline(-p.rhoRateMax, '--');
yline(p.rhoRateMax, '--');
grid on;
ylabel('d\rho/dt');
title('Scheduling-rate envelope');

subplot(3,1,3);
plot(t, deltaExact, 'LineWidth', 1.2);
hold on;
plot(t, deltaRateAware, '--', 'LineWidth', 1.2);
plot(t, deltaFrozen, ':', 'LineWidth', 1.4);
grid on;
xlabel('Time [s]');
ylabel('\delta');
legend( ...
    'Nonlinear perturbation', ...
    'Rate-aware local model', ...
    'Frozen local model', ...
    'Location', ...
    'best' ...
);
title('Effect of the moving-operating-point coupling term');

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

function value = equilibriumInput(rho, p)
    value = (p.a.*rho + p.b.*rho.^3)./p.g;
end

function value = localA(rho, p)
    value = -p.a - 3.0*p.b.*rho.^2;
end

function value = normalizedSchedule(rho, p)
    value = 2.0*(rho - p.rhoMin)/(p.rhoMax - p.rhoMin) - 1.0;
end

function value = schedule(time)
    value = 1.40*sin(0.25*time) + 0.35*sin(1.10*time);
end

function value = scheduleRate(time)
    value = 0.35*cos(0.25*time) + 0.385*cos(1.10*time);
end

function xNext = rk4Step(rhs, t, x, dt)
    k1 = rhs(t, x);
    k2 = rhs(t + 0.5*dt, x + 0.5*dt*k1);
    k3 = rhs(t + 0.5*dt, x + 0.5*dt*k2);
    k4 = rhs(t + dt, x + dt*k3);
    xNext = x + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
end

Chapter20_Lesson1.nb

Notebook[{
  Cell["Chapter 20, Lesson 1: Scheduling Variables and Operating Conditions", "Title"],
  Cell[
    BoxData @ ToBoxes @ Unevaluated[
      Module[
        {
          p, equilibriumInput, localA, normalizedSchedule, schedule,
          scheduleRate, rk4Step, dt, tf, time, n, rho, rhoDot, uEq,
          x, deltaExact, deltaRateAware, deltaFrozen, nonlinearRhs,
          rateAwareRhs, frozenRhs, insideValue, insideRate,
          rmseRateAware, rmseFrozen, rhoGrid, operatingTable,
          resultsTable, plot
        },

        p = <|
          "a" -> 1.2,
          "b" -> 0.35,
          "g" -> 1.0,
          "rhoMin" -> -2.0,
          "rhoMax" -> 2.0,
          "rhoRateMax" -> 0.8
        |>;

        equilibriumInput[r_] :=
          (p["a"] r + p["b"] r^3)/p["g"];

        localA[r_] :=
          -p["a"] - 3.0 p["b"] r^2;

        normalizedSchedule[r_] :=
          2.0 (r - p["rhoMin"])/(p["rhoMax"] - p["rhoMin"]) - 1.0;

        schedule[t_] :=
          1.40 Sin[0.25 t] + 0.35 Sin[1.10 t];

        scheduleRate[t_] :=
          0.35 Cos[0.25 t] + 0.385 Cos[1.10 t];

        rk4Step[rhs_, t_, state_, h_] := Module[
          {k1, k2, k3, k4},
          k1 = rhs[t, state];
          k2 = rhs[t + 0.5 h, state + 0.5 h k1];
          k3 = rhs[t + 0.5 h, state + 0.5 h k2];
          k4 = rhs[t + h, state + h k3];
          state + (h/6.0) (k1 + 2.0 k2 + 2.0 k3 + k4)
        ];

        dt = 0.01;
        tf = 30.0;
        time = N @ Range[0.0, tf, dt];
        n = Length[time];

        rho = schedule /@ time;
        rhoDot = scheduleRate /@ time;
        uEq = equilibriumInput /@ rho;

        x = ConstantArray[0.0, n];
        deltaRateAware = ConstantArray[0.0, n];
        deltaFrozen = ConstantArray[0.0, n];

        x[[1]] = rho[[1]] + 0.40;
        deltaRateAware[[1]] = x[[1]] - rho[[1]];
        deltaFrozen[[1]] = deltaRateAware[[1]];

        nonlinearRhs[t_, state_] := Module[
          {r, ue},
          r = schedule[t];
          ue = equilibriumInput[r];
          -p["a"] state - p["b"] state^3 + p["g"] ue
        ];

        rateAwareRhs[t_, delta_] :=
          localA[schedule[t]] delta - scheduleRate[t];

        frozenRhs[t_, delta_] :=
          localA[schedule[t]] delta;

        Do[
          x[[k + 1]] =
            rk4Step[nonlinearRhs, time[[k]], x[[k]], dt];

          deltaRateAware[[k + 1]] =
            rk4Step[
              rateAwareRhs,
              time[[k]],
              deltaRateAware[[k]],
              dt
            ];

          deltaFrozen[[k + 1]] =
            rk4Step[
              frozenRhs,
              time[[k]],
              deltaFrozen[[k]],
              dt
            ],
          {k, 1, n - 1}
        ];

        deltaExact = x - rho;

        insideValue =
          Map[
            p["rhoMin"] <= # <= p["rhoMax"] &,
            rho
          ];

        insideRate =
          Map[
            Abs[#] <= p["rhoRateMax"] &,
            rhoDot
          ];

        rmseRateAware =
          Sqrt[Mean[(deltaExact - deltaRateAware)^2]];

        rmseFrozen =
          Sqrt[Mean[(deltaExact - deltaFrozen)^2]];

        Print["Scheduling-variable diagnostics"];
        Print[
          "  rho range: ",
          NumberForm[{Min[rho], Max[rho]}, {8, 6}]
        ];
        Print[
          "  max |rho_dot|: ",
          NumberForm[Max[Abs[rhoDot]], {8, 6}]
        ];
        Print[
          "  samples inside value envelope: ",
          NumberForm[100.0 Mean[Boole /@ insideValue], {6, 2}],
          "%"
        ];
        Print[
          "  samples inside rate envelope:  ",
          NumberForm[100.0 Mean[Boole /@ insideRate], {6, 2}],
          "%"
        ];
        Print[
          "  RMSE, rate-aware local model:  ",
          NumberForm[rmseRateAware, {8, 6}]
        ];
        Print[
          "  RMSE, frozen local model:      ",
          NumberForm[rmseFrozen, {8, 6}]
        ];

        rhoGrid = N @ Subdivide[p["rhoMin"], p["rhoMax"], 16];

        operatingTable = Prepend[
          Table[
            {
              r,
              normalizedSchedule[r],
              equilibriumInput[r],
              localA[r],
              -1.0/localA[r]
            },
            {r, rhoGrid}
          ],
          {
            "rho",
            "rho_normalized",
            "equilibrium_input",
            "local_A",
            "time_constant"
          }
        ];

        Export[
          "Chapter20_Lesson1_operating_table_mathematica.csv",
          operatingTable
        ];

        resultsTable = Prepend[
          Transpose[
            {
              time,
              rho,
              rhoDot,
              uEq,
              x,
              deltaExact,
              deltaRateAware,
              deltaFrozen
            }
          ],
          {
            "time",
            "rho",
            "rho_dot",
            "equilibrium_input",
            "x",
            "delta_exact",
            "delta_rate_aware",
            "delta_frozen"
          }
        ];

        Export[
          "Chapter20_Lesson1_results_mathematica.csv",
          resultsTable
        ];

        plot = GraphicsGrid[
          {
            {
              ListLinePlot[
                {
                  Transpose[{time, rho}],
                  Transpose[
                    {
                      time,
                      ConstantArray[p["rhoMin"], n]
                    }
                  ],
                  Transpose[
                    {
                      time,
                      ConstantArray[p["rhoMax"], n]
                    }
                  ]
                },
                Frame -> True,
                FrameLabel -> {"time [s]", "rho"},
                PlotLegends -> {
                  "rho(t)",
                  "lower value bound",
                  "upper value bound"
                },
                ImageSize -> Large
              ]
            },
            {
              ListLinePlot[
                {
                  Transpose[{time, rhoDot}],
                  Transpose[
                    {
                      time,
                      ConstantArray[-p["rhoRateMax"], n]
                    }
                  ],
                  Transpose[
                    {
                      time,
                      ConstantArray[p["rhoRateMax"], n]
                    }
                  ]
                },
                Frame -> True,
                FrameLabel -> {"time [s]", "rho rate"},
                PlotLegends -> {
                  "rho_dot(t)",
                  "lower rate bound",
                  "upper rate bound"
                },
                ImageSize -> Large
              ]
            },
            {
              ListLinePlot[
                {
                  Transpose[{time, deltaExact}],
                  Transpose[{time, deltaRateAware}],
                  Transpose[{time, deltaFrozen}]
                },
                Frame -> True,
                FrameLabel -> {"time [s]", "delta"},
                PlotLegends -> {
                  "nonlinear perturbation",
                  "rate-aware local model",
                  "frozen local model"
                },
                ImageSize -> Large
              ]
            }
          }
        ];

        Export[
          "Chapter20_Lesson1_results_mathematica.png",
          plot,
          ImageResolution -> 180
        ];

        plot
      ]
    ],
    "Input"
  ]
},
WindowTitle -> "Chapter20_Lesson1.nb",
StyleDefinitions -> "Default.nb"
]

Problem 1 (Equilibrium Sensitivity): For \( \dot{x}=f(x,u,\rho) \), suppose an equilibrium family satisfies \( f(X_e(\rho),U_e(\rho),\rho)=0 \). Derive the sensitivity relation for \( X_e'(\rho) \) and \( U_e'(\rho) \).

Solution:

Differentiate the equilibrium identity with respect to \( \rho \). Applying the chain rule gives

\[ f_x X_e' + f_u U_e' + f_\rho=0. \]

Evaluated on the equilibrium manifold, this becomes

\[ A_e(\rho)X_e'(\rho) +B_e(\rho)U_e'(\rho) +F_\rho(\rho)=0. \]

When one of the equilibrium maps is prescribed and the corresponding Jacobian has suitable rank, the relation determines the derivative of the other map locally.

Problem 2 (Moving-Origin Term): Starting with \( \delta x=x-X_e(\rho) \), derive the first-order perturbation model around the moving equilibrium.

Solution:

Differentiation gives

\[ \delta\dot{x}=\dot{x}-X_\rho(\rho)\dot{\rho}. \]

Substitute \( x=X_e+\delta x \), \( u=U_e+\delta u \), and expand \( f \) to first order. Because the zeroth-order term is zero on the equilibrium manifold,

\[ \delta\dot{x} = A(\rho)\delta x +B(\rho)\delta u -X_\rho(\rho)\dot{\rho} +r_2. \]

The last first-order term is absent only when the equilibrium is fixed or when \( X_e \) does not depend on the schedule.

Problem 3 (Schedule-Rate Bound): For \( \dot{z}=-\alpha z-\dot{\rho} \), where \( \alpha>0 \) and \( |\dot{\rho}|\leq\nu \), prove that

\[ |z(t)| \leq e^{-\alpha t}|z(0)| + \frac{\nu}{\alpha} (1-e^{-\alpha t}). \]

Solution:

The exact solution is

\[ z(t)=e^{-\alpha t}z(0) - \int_0^t e^{-\alpha(t-\tau)} \dot{\rho}(\tau)\,d\tau. \]

Apply the triangle inequality and the bound on \( \dot{\rho} \):

\[ |z(t)| \leq e^{-\alpha t}|z(0)| + \nu\int_0^t e^{-\alpha(t-\tau)}d\tau. \]

Evaluation of the integral yields the required result. The ultimate bound is \( \nu/\alpha \).

Problem 4 (Nonlinear Example): For \( \dot{x}=-ax-bx^3+gu \), derive \( u_e(\rho) \), \( A(\rho) \), and the exact perturbation dynamics under \( x=\rho+\delta \).

Solution:

The equilibrium condition gives \( 0=-a\rho-b\rho^3+gu_e \), so

\[ u_e(\rho)=\frac{a\rho+b\rho^3}{g}. \]

Linearization with respect to \( x \) gives

\[ A(\rho)=-a-3b\rho^2. \]

Substitution of \( x=\rho+\delta \) and \( u=u_e(\rho) \), followed by cancellation of the equilibrium terms, gives

\[ \dot{\delta} = (-a-3b\rho^2)\delta -3b\rho\delta^2 -b\delta^3 -\dot{\rho}. \]

Problem 5 (Filtering Trade-Off): A measured schedule is filtered using \( \tau_f\dot{\hat{\rho}}+\hat{\rho}=\rho \). Let \( e=\hat{\rho}-\rho \) and assume \( |\dot{\rho}|\leq\nu \). Derive an ultimate bound on \( e \).

Solution:

The error dynamics are

\[ \tau_f\dot{e}+e=-\tau_f\dot{\rho}. \]

This is a stable first-order system driven by an input whose magnitude is bounded by \( \tau_f\nu \). Its solution satisfies

\[ |e(t)| \leq e^{-t/\tau_f}|e(0)| + \tau_f\nu (1-e^{-t/\tau_f}). \]

Consequently, \( \limsup_{t\to\infty}|e(t)|\leq\tau_f\nu \). A larger filter time constant attenuates more noise but permits more schedule lag.

Problem 6 (Two-Dimensional Scheduling Grid): Suppose the dynamics depend on speed \( v \) and payload \( m_p \). Explain why independently sampling a full rectangular grid may be wasteful when payload and achievable speed are physically coupled.

Solution:

A rectangular grid includes combinations that may be unreachable, such as maximum speed at maximum payload. Local models at unreachable points have no physical validation data and can distort later interpolation. The schedule should instead be defined on the feasible subset \( \mathcal{P}_{feas} \) determined by plant limits, with boundary margins and enough points to resolve coefficient variation. If a rectangular representation is required, unreachable vertices must be handled by a justified embedding rather than by invented data.

Scheduling variables parameterize known changes in operating condition; they are conceptually different from unknown-parameter estimates in adaptive control. A rigorous schedule begins with an equilibrium manifold, local sensitivity analysis, and a validated admissible set for both schedule values and rates. Linearization around a moving equilibrium produces the additional term \( -X_\rho(\rho)\dot{\rho} \), which explains why frozen models may be inaccurate during rapid transitions. The next lesson compares this scheduled architecture directly with true adaptive control.