A practical gain-scheduled controller is a nonlinear controller whose coefficients depend on a measurable operating-condition vector \( \boldsymbol{\rho}(t) \in \mathcal{P} \subset \mathbb{R}^{n_\rho} \). Typical aircraft scheduling variables include airspeed, Mach number, dynamic pressure, altitude, and configuration state. Typical robotic scheduling variables include payload estimate, arm extension, joint configuration, battery voltage, and contact mode.

The complete implementation contains more than a gain table. It includes signal validation, normalization, filtering, rate limiting, interpolation, trim/feedforward scheduling, actuator limiting, monitoring, and fallback logic.

flowchart TD
  R["Reference command"] --> C["Scheduled controller"]
  Y["Measured states and outputs"] --> C
  S["Raw scheduling measurements"] --> V["Validation and plausibility checks"]
  V --> F["Filtering and rate limiting"]
  F --> I["Controller and trim interpolation"]
  I --> C
  C --> L["Actuator limits and safety logic"]
  L --> P["Aircraft or robot plant"]
  P --> Y
  M["Envelope monitor and fallback"] --> C
  F --> M
        

The central distinction from true adaptive control is that the map from operating condition to controller coefficients is designed offline. During operation, the controller evaluates this predesigned map; it does not infer unknown parameters through an online adaptation law.

Begin with a nonlinear plant model parameterized by the measurable operating condition \( \boldsymbol{\rho} \):

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

At each design point \( \boldsymbol{\rho}_i \), compute a trim pair \( (\bar{\mathbf{x}}_i,\bar{\mathbf{u}}_i) \) satisfying

\[ \mathbf{f}(\bar{\mathbf{x}}_i,\bar{\mathbf{u}}_i,\boldsymbol{\rho}_i)=\mathbf{0}. \]

For trajectory-dependent scheduling, the nominal point may be time-varying; the local perturbation coordinates are still formed relative to the nominal family:

\[ \delta\mathbf{x}=\mathbf{x}-\bar{\mathbf{x}}(\boldsymbol{\rho}), \qquad \delta\mathbf{u}=\mathbf{u}-\bar{\mathbf{u}}(\boldsymbol{\rho}). \]

Linearization at a fixed design point gives

\[ \delta\dot{\mathbf{x}} = \mathbf{A}_i\delta\mathbf{x} +\mathbf{B}_i\delta\mathbf{u}, \qquad \delta\mathbf{y}=\mathbf{C}_i\delta\mathbf{x} +\mathbf{D}_i\delta\mathbf{u}, \]

\[ \mathbf{A}_i=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right|_i, \quad \mathbf{B}_i=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\right|_i, \quad \mathbf{C}_i=\left.\frac{\partial\mathbf{h}}{\partial\mathbf{x}}\right|_i, \quad \mathbf{D}_i=\left.\frac{\partial\mathbf{h}}{\partial\mathbf{u}}\right|_i. \]

All local models must use the same state ordering, input ordering, output definitions, units, signs, actuator conventions, and controller-state realization. Interpolating incompatible coordinate systems produces a mathematically smooth table but a physically meaningless controller.

A useful scheduling variable should satisfy five engineering conditions:

1. It is measured or estimated causally at the controller update rate.
2. It explains a substantial part of the plant-dynamics variation.
3. Its admissible range and rate of variation can be bounded.
4. Its noise and latency are compatible with the desired scheduling bandwidth.
5. Its failure can be detected and handled by a defined fallback mode.

Normalize each scalar scheduling coordinate to reduce numerical scaling problems:

\[ \eta_j = \frac{\rho_j-\rho_{j,\min}}{\rho_{j,\max}-\rho_{j,\min}}, \qquad 0\le \eta_j\le 1. \]

A measured scheduling signal is normally filtered and rate-limited before entering the gain table. A first-order filter with a symmetric rate limit is

\[ \dot{\rho}_{f,j}= \operatorname{sat}_{[-\nu_j,\nu_j]} \left(\frac{\rho_j-\rho_{f,j}}{\tau_j}\right), \]

where \( \tau_j \) is the filter time constant and \( \nu_j \) is the maximum scheduling rate admitted by the implementation. The filtered value must also be clamped to the certified envelope. Filtering reduces gain chatter but introduces lag; therefore the nonlinear verification campaign must include rapid operating-condition changes.

Endogenous scheduling variables, such as angle of attack, joint angle, or a plant output, deserve special attention because the scheduling map then forms an additional nonlinear feedback path. Frozen-point stability alone does not certify this closed loop.

Let the certified envelope be partitioned by design nodes \( \{\boldsymbol{\rho}_1,\ldots,\boldsymbol{\rho}_N\} \). Node placement should be dense where plant matrices, trim quantities, actuator effectiveness, or desired bandwidth change rapidly. A rectangular uniform grid is convenient but often wasteful; a physics-informed nonuniform grid is usually more efficient.

At each node, design a controller with a common architecture. For a state-feedback controller with scheduled feedforward,

\[ \mathbf{u}=\bar{\mathbf{u}}(\boldsymbol{\rho}) -\mathbf{K}(\boldsymbol{\rho}) \left[\mathbf{x}-\bar{\mathbf{x}}(\boldsymbol{\rho})\right] +\mathbf{K}_r(\boldsymbol{\rho})\mathbf{r}. \]

Local controllers should be coordinated rather than tuned independently. Useful coordination constraints include consistent closed-loop bandwidth, damping ratios, integral-action conventions, estimator poles, actuator usage, and sign/phase behavior. Abruptly different local designs are difficult to interpolate safely even when every node is stable.

4.1 Closed-Form Local Design for the Lesson Example

The code laboratories use a second-order axis that can represent an aircraft attitude channel or a robotic joint/position channel:

\[ \ddot{x}+a(\rho)\dot{x}+b(\rho)x=g(\rho)u+d_{nl}(x), \]

with the scheduled control law

\[ u=-k_1(\rho)x-k_2(\rho)\dot{x}+k_r(\rho)r. \]

Ignoring the bounded nonlinear residual during local design, the closed-loop characteristic polynomial is

\[ s^2+\left[a(\rho)+g(\rho)k_2(\rho)\right]s +\left[b(\rho)+g(\rho)k_1(\rho)\right]. \]

Matching it to \( s^2+2\zeta(\rho)\omega_n(\rho)s+\omega_n^2(\rho) \) gives

\[ k_1(\rho)=\frac{\omega_n^2(\rho)-b(\rho)}{g(\rho)}, \qquad k_2(\rho)=\frac{2\zeta(\rho)\omega_n(\rho)-a(\rho)}{g(\rho)}. \]

Requiring unit steady-state gain for a constant command yields

\[ k_r(\rho)=\frac{\omega_n^2(\rho)}{g(\rho)}. \]

These equations make the node designs transparent and allow the same algorithm to be implemented from scratch in every programming language.

For a scalar schedule between nodes \( \rho_i \) and \( \rho_{i+1} \), define

\[ \alpha(\rho)=\frac{\rho-\rho_i}{\rho_{i+1}-\rho_i}, \qquad 0\le\alpha\le 1, \]

\[ \mathbf{K}(\rho)=(1-\alpha)\mathbf{K}_i+\alpha\mathbf{K}_{i+1}. \]

For a polytopic envelope with barycentric coordinates \( \alpha_i(\boldsymbol{\rho}) \),

\[ \alpha_i\ge 0, \qquad \sum_{i=1}^{N_v}\alpha_i=1, \qquad \mathbf{K}(\boldsymbol{\rho})= \sum_{i=1}^{N_v}\alpha_i(\boldsymbol{\rho})\mathbf{K}_i. \]

The trim state, trim input, feedforward map, estimator matrices, and anti-windup parameters should be scheduled consistently with the feedback gains. Otherwise, a controller may be locally well designed yet generate command jumps during transitions.

5.1 Why Interpolating Stable Gains Is Not Automatically Safe

Suppose the plant and controller are both interpolated between two nodes:

\[ \mathbf{A}(\alpha)=(1-\alpha)\mathbf{A}_1+\alpha\mathbf{A}_2, \quad \mathbf{B}(\alpha)=(1-\alpha)\mathbf{B}_1+\alpha\mathbf{B}_2, \quad \mathbf{K}(\alpha)=(1-\alpha)\mathbf{K}_1+\alpha\mathbf{K}_2. \]

The actual closed-loop matrix is

\[ \mathbf{A}_{cl}(\alpha)=\mathbf{A}(\alpha)-\mathbf{B}(\alpha)\mathbf{K}(\alpha). \]

Expanding the product gives

\[ \begin{aligned} \mathbf{B}(\alpha)\mathbf{K}(\alpha) ={}&(1-\alpha)^2\mathbf{B}_1\mathbf{K}_1 +\alpha^2\mathbf{B}_2\mathbf{K}_2 \\ &+\alpha(1-\alpha) \left(\mathbf{B}_1\mathbf{K}_2+\mathbf{B}_2\mathbf{K}_1\right). \end{aligned} \]

The cross terms show that \( \mathbf{A}_{cl}(\alpha) \) is generally not the convex interpolation of the endpoint closed-loop matrices. Therefore endpoint stability is insufficient. One must either use a synthesis method that preserves the required convex structure or verify the true closed-loop model throughout the envelope.

5.2 Dynamic Controllers and State-Realization Compatibility

For a dynamic controller

\[ \dot{\mathbf{x}}_c=\mathbf{A}_c(\rho)\mathbf{x}_c +\mathbf{B}_c(\rho)\mathbf{e}, \qquad \mathbf{u}=\mathbf{C}_c(\rho)\mathbf{x}_c +\mathbf{D}_c(\rho)\mathbf{e}, \]

two equivalent local transfer functions can have incompatible controller-state coordinates. Direct matrix interpolation can then create unintended internal dynamics. A practical design must impose a common realization, use an interpolation method with stability guarantees, or reset/map controller states through a bumpless-transfer mechanism.

6.1 Common Quadratic Lyapunov Certificate

Assume the closed-loop matrix is genuinely polytopic:

\[ \mathbf{A}_{cl}(\boldsymbol{\rho})= \sum_{i=1}^{N_v}\alpha_i(\boldsymbol{\rho})\mathbf{A}_{cl,i}, \qquad \alpha_i\ge 0, \qquad \sum_i\alpha_i=1. \]

If there exists a common matrix \( \mathbf{P}=\mathbf{P}^{T}\succ 0 \) and \( \mathbf{Q}=\mathbf{Q}^{T}\succ 0 \) such that

\[ \mathbf{A}_{cl,i}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{cl,i} \preceq -\mathbf{Q}, \qquad i=1,\ldots,N_v, \]

then the scheduled system is quadratically stable for arbitrary scheduling speed. To prove this, choose \( V=\mathbf{x}^{T}\mathbf{P}\mathbf{x} \). By convexity,

\[ \begin{aligned} \dot{V} &=\mathbf{x}^{T} \left(\mathbf{A}_{cl}^{T}\mathbf{P}+ \mathbf{P}\mathbf{A}_{cl}\right)\mathbf{x} \\ &=\sum_i\alpha_i\mathbf{x}^{T} \left(\mathbf{A}_{cl,i}^{T}\mathbf{P}+ \mathbf{P}\mathbf{A}_{cl,i}\right)\mathbf{x} \\ &\le -\mathbf{x}^{T}\mathbf{Q}\mathbf{x}. \end{aligned} \]

Thus \( \dot{V} \) is negative definite uniformly over the parameter polytope. The certificate can be conservative because a single ellipsoidal Lyapunov function must work for every operating condition.

6.2 Parameter-Dependent Lyapunov Function and Rate Bounds

A less conservative analysis may use \( V=\mathbf{x}^{T}\mathbf{P}(\boldsymbol{\rho})\mathbf{x} \). Its derivative contains an explicit scheduling-rate term:

\[ \dot{V}=\mathbf{x}^{T} \left[ \mathbf{A}_{cl}^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_{cl} +\sum_{j=1}^{n_\rho} \frac{\partial\mathbf{P}}{\partial\rho_j}\dot{\rho}_j \right]\mathbf{x}. \]

If the implementation guarantees \( |\dot{\rho}_j|\le\nu_j \), sufficient stability conditions require the bracketed matrix to remain negative definite for all admissible parameter values and rate vertices. This formula explains why scheduling-rate limits are part of the mathematical design, not merely a software convenience.

6.3 Frozen-Time Eigenvalues Are Necessary but Not Sufficient

Checking that every frozen matrix \( \mathbf{A}_{cl}(\boldsymbol{\rho}) \) is Hurwitz is an important verification step. However, a time-varying system can be unstable even when every frozen system is stable. A rigorous workflow therefore combines frozen-point analysis, Lyapunov/rate analysis where feasible, dense-grid checks, and nonlinear transition simulations.

7.1 Deterministic Table Evaluation

Store grid coordinates and controller data in immutable, versioned tables. At each control cycle: validate the schedule, clamp it to the certified envelope, identify the enclosing cell, compute interpolation weights, evaluate gains and trim values, apply control, enforce actuator limits, and record health flags. The table-search execution time must be bounded.

7.2 Bumpless Scheduling

For static state feedback, continuous interpolation usually gives continuous gains. For controllers with integrators, observers, filters, or notch states, bumpless transfer requires controller-state consistency. If a controller state is reparameterized by a schedule-dependent transformation \( \mathbf{x}_c=\mathbf{T}(\rho)\mathbf{z}_c \), then

\[ \dot{\mathbf{x}}_c=\mathbf{T}(\rho)\dot{\mathbf{z}}_c +\frac{\partial\mathbf{T}}{\partial\rho}\dot{\rho}\mathbf{z}_c. \]

Neglecting the second term creates a hidden coupling proportional to the scheduling rate. This is one reason to preserve a common controller realization across all nodes.

7.3 Saturation and Anti-Windup

Actuator position, rate, current, torque, and thermal limits vary across many aircraft and robotic envelopes. Limit values may themselves be scheduled. Because saturation invalidates the local linear design, integral states should use anti-windup or conditional integration, and the verification campaign must include prolonged saturation and recovery.

7.4 Fault Handling

Define explicit behavior for stale, impossible, discontinuous, or mutually inconsistent scheduling measurements. Typical responses are hold-last-valid, transition to a conservative fixed controller, command limiting, or controlled shutdown. Silent extrapolation beyond the certified envelope should be avoided.

8.1 Aircraft Example

Consider longitudinal control across speed and altitude. A practical schedule may use \( \boldsymbol{\rho}=[M,\bar{q},h,c_f]^{T} \), where \( M \) is Mach number, \( \bar{q} \) is dynamic pressure, \( h \) is altitude, and \( c_f \) is a discrete configuration indicator. The design points are trimmed flight conditions. Local controllers should coordinate short-period damping, command bandwidth, control-surface rate, structural-mode separation, and load-factor limits.

Configuration changes such as flap or landing-gear transitions are often better represented by explicit modes with guarded transitions than by naively treating every discrete variable as continuous. Within each mode, continuous gain scheduling can be applied over speed and altitude.

8.2 Robotics Example

For a manipulator mounted on a mobile base, a useful schedule may contain payload estimate, arm extension, selected joint coordinates, and battery-bus voltage. The local models capture inertia variation, gravity sensitivity, actuator effectiveness, and flexible-mode movement. Controller bandwidth is usually reduced at large extension or high payload to preserve torque margin and avoid exciting structural modes.

Scheduling directly on every joint angle can create a high-dimensional table. A practical reduction uses physically meaningful aggregate variables, such as effective reflected inertia or horizontal reach, while retaining nonlinear simulation to detect unmodeled configuration dependence.

8.3 Common Certification-Oriented Data

For either domain, preserve the model version, trim solver settings, grid, controller-design settings, table checksum, test cases, expected metrics, and software version. Reproducibility is essential because a gain table is an executable engineering artifact, not merely a tuning worksheet.

flowchart TD
  E["Define certified operating envelope"] --> T["Compute trims and compatible local models"]
  T --> D["Design coordinated local controllers"]
  D --> N["Verify every design node"]
  N --> G["Check dense off-grid frozen points"]
  G --> R["Analyze parameter-rate and transition effects"]
  R --> S["Run nonlinear nominal simulations"]
  S --> U["Run uncertainty and Monte Carlo cases"]
  U --> H["Run sensor fault, latency, saturation, and fallback tests"]
  H --> C["Code-level and hardware-in-the-loop validation"]
  C --> A["Approve table, limits, and evidence package"]
        

9.1 Frozen-Point Tests

Evaluate closed-loop poles, damping ratios, bandwidth, sensitivity, control effort, estimator dynamics, and stability margins at all nodes and at a dense set of off-grid points. Off-grid points are mandatory because interpolation can produce worse behavior between locally acceptable designs.

9.2 Transition Tests

Exercise ramps, steps, reversals, and noisy trajectories of the scheduling variables at their maximum admissible rates. Include simultaneous reference commands and disturbances so that transitions occur away from equilibrium.

9.3 Nonlinear and Uncertainty Tests

Use the highest-fidelity nonlinear model available. Vary aerodynamic coefficients, mass properties, payload, friction, sensor bias, delay, actuator dynamics, structural modes, and environmental disturbances. For sampled cases, useful aggregate metrics include

\[ e_{\mathrm{RMS}}= \sqrt{\frac{1}{T}\int_0^T e^2(t)\,dt}, \qquad e_{\max}=\max_{t\in[0,T]}|e(t)|, \qquad u_{\max}=\max_{t\in[0,T]}|u(t)|. \]

Pass/fail criteria should also include envelope violations, saturation duration, recovery time, gain discontinuity, schedule-clamping duration, and fallback activation.

9.4 Software and Hardware Tests

Unit-test interpolation at nodes, cell boundaries, and out-of-range values. Compare generated code against the design environment. Then perform processor-in-the-loop and hardware-in-the-loop tests with realistic timing, quantization, sensor transport, actuator interfaces, and fault injection.

The standalone implementation below uses only the Python standard library. For larger projects, useful libraries include NumPy and SciPy for numerical work, python-control for state-space and frequency-domain analysis, and CVXPY for LMI-based offline design. The script creates the local gain table, filters and rate-limits the schedule, performs piecewise-linear interpolation, simulates the nonlinear plant by RK4, and writes validation data to CSV.

Chapter20_Lesson5.py

"""Chapter20_Lesson5.py

Practical gain-scheduled control of a second-order aircraft/robotics axis.
The scheduling variable rho is filtered and rate-limited before interpolating
precomputed controller gains. Results are written to CSV for inspection.
"""

from __future__ import annotations

import csv
import math
from dataclasses import dataclass
from typing import Iterable, Sequence


@dataclass(frozen=True)
class GainNode:
    rho: float
    k1: float
    k2: float
    kr: float


def clamp(value: float, lower: float, upper: float) -> float:
    return max(lower, min(upper, value))


def plant_coefficients(rho: float) -> tuple[float, float, float]:
    """Return damping a(rho), stiffness b(rho), and input gain g(rho)."""
    a = 0.8 + 0.7 * rho
    b = 1.4 - 0.5 * rho
    g = 1.0 / (1.0 + 0.8 * rho)
    return a, b, g


def desired_dynamics(rho: float) -> tuple[float, float]:
    """Desired damping ratio and natural frequency at a design node."""
    zeta = 0.85
    omega_n = 2.2 - 0.4 * rho
    return zeta, omega_n


def design_node(rho: float) -> GainNode:
    """Pole-placement formulas for u = -k1*x1 - k2*x2 + kr*r."""
    a, b, g = plant_coefficients(rho)
    zeta, omega_n = desired_dynamics(rho)
    k1 = (omega_n**2 - b) / g
    k2 = (2.0 * zeta * omega_n - a) / g
    kr = omega_n**2 / g
    return GainNode(rho, k1, k2, kr)


def interpolate(nodes: Sequence[GainNode], rho: float) -> GainNode:
    """Piecewise-linear interpolation with endpoint clamping."""
    if len(nodes) < 2:
        raise ValueError("At least two gain nodes are required.")

    rho_c = clamp(rho, nodes[0].rho, nodes[-1].rho)
    for left, right in zip(nodes[:-1], nodes[1:]):
        if rho_c <= right.rho:
            span = right.rho - left.rho
            if span <= 0.0:
                raise ValueError("Gain nodes must have strictly increasing rho.")
            alpha = (rho_c - left.rho) / span
            lerp = lambda x0, x1: (1.0 - alpha) * x0 + alpha * x1
            return GainNode(
                rho=rho_c,
                k1=lerp(left.k1, right.k1),
                k2=lerp(left.k2, right.k2),
                kr=lerp(left.kr, right.kr),
            )
    return nodes[-1]


def reference_command(t: float) -> float:
    step = 0.20 if t >= 8.0 else 0.0
    return 0.40 * math.sin(0.30 * t) + step


def raw_schedule(t: float) -> float:
    value = 0.50 + 0.40 * math.sin(0.12 * t) + 0.08 * math.sin(0.90 * t)
    return clamp(value, 0.0, 1.0)


def schedule_rate(rho_raw: float, rho_filtered: float, tau: float, rate_limit: float) -> float:
    requested_rate = (rho_raw - rho_filtered) / tau
    return clamp(requested_rate, -rate_limit, rate_limit)


def state_derivative(x1: float, x2: float, control: float, rho: float) -> tuple[float, float]:
    a, b, g = plant_coefficients(rho)
    nonlinear_term = 0.08 * x1**3
    return x2, -a * x2 - b * x1 - nonlinear_term + g * control


def rk4_step(x1: float, x2: float, control: float, rho: float, dt: float) -> tuple[float, float]:
    def f(s1: float, s2: float) -> tuple[float, float]:
        return state_derivative(s1, s2, control, rho)

    k11, k12 = f(x1, x2)
    k21, k22 = f(x1 + 0.5 * dt * k11, x2 + 0.5 * dt * k12)
    k31, k32 = f(x1 + 0.5 * dt * k21, x2 + 0.5 * dt * k22)
    k41, k42 = f(x1 + dt * k31, x2 + dt * k32)

    next_x1 = x1 + dt * (k11 + 2.0 * k21 + 2.0 * k31 + k41) / 6.0
    next_x2 = x2 + dt * (k12 + 2.0 * k22 + 2.0 * k32 + k42) / 6.0
    return next_x1, next_x2


def rms(values: Iterable[float]) -> float:
    values = list(values)
    return math.sqrt(sum(value * value for value in values) / len(values))


def main() -> None:
    nodes = [design_node(rho) for rho in (0.0, 0.5, 1.0)]
    for node in nodes:
        print(
            f"node rho={node.rho:.2f}: k1={node.k1:.6f}, "
            f"k2={node.k2:.6f}, kr={node.kr:.6f}"
        )

    dt = 0.002
    final_time = 40.0
    tau_schedule = 0.25
    schedule_rate_limit = 0.90
    control_limit = 8.0

    x1 = 0.0
    x2 = 0.0
    rho_filtered = raw_schedule(0.0)
    rows: list[tuple[float, ...]] = []
    errors: list[float] = []
    controls: list[float] = []

    steps = int(round(final_time / dt))
    for index in range(steps + 1):
        t = index * dt
        rho_raw = raw_schedule(t)
        rho_dot = schedule_rate(rho_raw, rho_filtered, tau_schedule, schedule_rate_limit)
        rho_filtered = clamp(rho_filtered + dt * rho_dot, 0.0, 1.0)

        gains = interpolate(nodes, rho_filtered)
        reference = reference_command(t)
        control_unsaturated = -gains.k1 * x1 - gains.k2 * x2 + gains.kr * reference
        control = clamp(control_unsaturated, -control_limit, control_limit)

        error = reference - x1
        rows.append(
            (
                t,
                reference,
                x1,
                x2,
                error,
                control,
                rho_raw,
                rho_filtered,
                gains.k1,
                gains.k2,
                gains.kr,
            )
        )
        errors.append(error)
        controls.append(control)

        x1, x2 = rk4_step(x1, x2, control, rho_raw, dt)

    output_name = "Chapter20_Lesson5_results.csv"
    with open(output_name, "w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(
            [
                "time",
                "reference",
                "position",
                "rate",
                "tracking_error",
                "control",
                "rho_raw",
                "rho_filtered",
                "k1",
                "k2",
                "kr",
            ]
        )
        writer.writerows(rows)

    print(f"RMS tracking error: {rms(errors):.6f}")
    print(f"Peak absolute error: {max(abs(value) for value in errors):.6f}")
    print(f"Peak absolute control: {max(abs(value) for value in controls):.6f}")
    print(f"Wrote {output_name}")


if __name__ == "__main__":
    main()

The C++17 version uses the standard library and deterministic table search. In production control software, Eigen is commonly used for matrix operations; offline optimization can be connected to solvers such as OSQP or commercial LMI tools. Avoid dynamic allocation inside the real-time loop when targeting an embedded controller.

Chapter20_Lesson5.cpp

// Chapter20_Lesson5.cpp
// Practical gain-scheduled control of a second-order aircraft/robotics axis.

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

struct GainNode {
    double rho;
    double k1;
    double k2;
    double kr;
};

struct State {
    double x1;
    double x2;
};

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

std::array<double, 3> plantCoefficients(double rho) {
    const double a = 0.8 + 0.7 * rho;
    const double b = 1.4 - 0.5 * rho;
    const double g = 1.0 / (1.0 + 0.8 * rho);
    return {a, b, g};
}

std::array<double, 2> desiredDynamics(double rho) {
    const double zeta = 0.85;
    const double omegaN = 2.2 - 0.4 * rho;
    return {zeta, omegaN};
}

GainNode designNode(double rho) {
    const auto [a, b, g] = plantCoefficients(rho);
    const auto [zeta, omegaN] = desiredDynamics(rho);
    const double k1 = (omegaN * omegaN - b) / g;
    const double k2 = (2.0 * zeta * omegaN - a) / g;
    const double kr = omegaN * omegaN / g;
    return {rho, k1, k2, kr};
}

GainNode interpolate(const std::vector<GainNode>& nodes, double rho) {
    if (nodes.size() < 2) {
        throw std::invalid_argument("At least two gain nodes are required.");
    }

    const double rhoClamped = clampValue(rho, nodes.front().rho, nodes.back().rho);
    for (std::size_t i = 0; i + 1 < nodes.size(); ++i) {
        const GainNode& left = nodes[i];
        const GainNode& right = nodes[i + 1];
        if (rhoClamped <= right.rho) {
            const double span = right.rho - left.rho;
            if (span <= 0.0) {
                throw std::invalid_argument("Gain nodes must have strictly increasing rho.");
            }
            const double alpha = (rhoClamped - left.rho) / span;
            const auto lerp = [alpha](double x0, double x1) {
                return (1.0 - alpha) * x0 + alpha * x1;
            };
            return {rhoClamped,
                    lerp(left.k1, right.k1),
                    lerp(left.k2, right.k2),
                    lerp(left.kr, right.kr)};
        }
    }
    return nodes.back();
}

double referenceCommand(double t) {
    const double step = (t >= 8.0) ? 0.20 : 0.0;
    return 0.40 * std::sin(0.30 * t) + step;
}

double rawSchedule(double t) {
    const double value = 0.50 + 0.40 * std::sin(0.12 * t) + 0.08 * std::sin(0.90 * t);
    return clampValue(value, 0.0, 1.0);
}

double scheduleRate(double rhoRaw, double rhoFiltered, double tau, double rateLimit) {
    const double requestedRate = (rhoRaw - rhoFiltered) / tau;
    return clampValue(requestedRate, -rateLimit, rateLimit);
}

State derivative(const State& state, double control, double rho) {
    const auto [a, b, g] = plantCoefficients(rho);
    const double nonlinearTerm = 0.08 * state.x1 * state.x1 * state.x1;
    return {state.x2, -a * state.x2 - b * state.x1 - nonlinearTerm + g * control};
}

State addScaled(const State& state, const State& increment, double scale) {
    return {state.x1 + scale * increment.x1, state.x2 + scale * increment.x2};
}

State rk4Step(const State& state, double control, double rho, double dt) {
    const State k1 = derivative(state, control, rho);
    const State k2 = derivative(addScaled(state, k1, 0.5 * dt), control, rho);
    const State k3 = derivative(addScaled(state, k2, 0.5 * dt), control, rho);
    const State k4 = derivative(addScaled(state, k3, dt), control, rho);

    return {
        state.x1 + dt * (k1.x1 + 2.0 * k2.x1 + 2.0 * k3.x1 + k4.x1) / 6.0,
        state.x2 + dt * (k1.x2 + 2.0 * k2.x2 + 2.0 * k3.x2 + k4.x2) / 6.0
    };
}

int main() {
    try {
        std::vector<GainNode> nodes{designNode(0.0), designNode(0.5), designNode(1.0)};
        std::cout << std::fixed << std::setprecision(6);
        for (const auto& node : nodes) {
            std::cout << "node rho=" << node.rho << ": k1=" << node.k1
                      << ", k2=" << node.k2 << ", kr=" << node.kr << '\n';
        }

        constexpr double dt = 0.002;
        constexpr double finalTime = 40.0;
        constexpr double tauSchedule = 0.25;
        constexpr double scheduleRateLimit = 0.90;
        constexpr double controlLimit = 8.0;

        State state{0.0, 0.0};
        double rhoFiltered = rawSchedule(0.0);
        double sumSquaredError = 0.0;
        double peakError = 0.0;
        double peakControl = 0.0;
        std::size_t sampleCount = 0;

        std::ofstream csv("Chapter20_Lesson5_results.csv");
        if (!csv) {
            throw std::runtime_error("Could not open output CSV file.");
        }
        csv << "time,reference,position,rate,tracking_error,control,rho_raw,rho_filtered,k1,k2,kr\n";
        csv << std::setprecision(12);

        const int steps = static_cast<int>(std::llround(finalTime / dt));
        for (int index = 0; index <= steps; ++index) {
            const double t = index * dt;
            const double rhoRaw = rawSchedule(t);
            const double rhoDot = scheduleRate(rhoRaw, rhoFiltered, tauSchedule, scheduleRateLimit);
            rhoFiltered = clampValue(rhoFiltered + dt * rhoDot, 0.0, 1.0);

            const GainNode gains = interpolate(nodes, rhoFiltered);
            const double reference = referenceCommand(t);
            const double controlUnsaturated =
                -gains.k1 * state.x1 - gains.k2 * state.x2 + gains.kr * reference;
            const double control = clampValue(controlUnsaturated, -controlLimit, controlLimit);
            const double error = reference - state.x1;

            csv << t << ',' << reference << ',' << state.x1 << ',' << state.x2 << ','
                << error << ',' << control << ',' << rhoRaw << ',' << rhoFiltered << ','
                << gains.k1 << ',' << gains.k2 << ',' << gains.kr << '\n';

            sumSquaredError += error * error;
            peakError = std::max(peakError, std::abs(error));
            peakControl = std::max(peakControl, std::abs(control));
            ++sampleCount;

            state = rk4Step(state, control, rhoRaw, dt);
        }

        const double rmsError = std::sqrt(sumSquaredError / static_cast<double>(sampleCount));
        std::cout << "RMS tracking error: " << rmsError << '\n';
        std::cout << "Peak absolute error: " << peakError << '\n';
        std::cout << "Peak absolute control: " << peakControl << '\n';
        std::cout << "Wrote Chapter20_Lesson5_results.csv\n";
        return 0;
    } catch (const std::exception& exception) {
        std::cerr << "Error: " << exception.what() << '\n';
        return 1;
    }
}

The Java implementation uses records for immutable gain nodes and states. EJML is suitable for dense and sparse matrix computations, while Hipparchus provides numerical algorithms and ODE utilities. For hard real-time use, object creation, garbage collection, and file I/O should be removed from the control task and confined to initialization or telemetry threads.

Chapter20_Lesson5.java

// Chapter20_Lesson5.java
// Practical gain-scheduled control of a second-order aircraft/robotics axis.

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.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter20_Lesson5 {
    private Chapter20_Lesson5() {}

    private record GainNode(double rho, double k1, double k2, double kr) {}
    private record State(double x1, double x2) {}

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

    private static double[] plantCoefficients(double rho) {
        double a = 0.8 + 0.7 * rho;
        double b = 1.4 - 0.5 * rho;
        double g = 1.0 / (1.0 + 0.8 * rho);
        return new double[] {a, b, g};
    }

    private static double[] desiredDynamics(double rho) {
        double zeta = 0.85;
        double omegaN = 2.2 - 0.4 * rho;
        return new double[] {zeta, omegaN};
    }

    private static GainNode designNode(double rho) {
        double[] plant = plantCoefficients(rho);
        double[] desired = desiredDynamics(rho);
        double a = plant[0];
        double b = plant[1];
        double g = plant[2];
        double zeta = desired[0];
        double omegaN = desired[1];

        double k1 = (omegaN * omegaN - b) / g;
        double k2 = (2.0 * zeta * omegaN - a) / g;
        double kr = omegaN * omegaN / g;
        return new GainNode(rho, k1, k2, kr);
    }

    private static GainNode interpolate(List<GainNode> nodes, double rho) {
        if (nodes.size() < 2) {
            throw new IllegalArgumentException("At least two gain nodes are required.");
        }
        double rhoClamped = clamp(rho, nodes.get(0).rho(), nodes.get(nodes.size() - 1).rho());

        for (int i = 0; i + 1 < nodes.size(); i++) {
            GainNode left = nodes.get(i);
            GainNode right = nodes.get(i + 1);
            if (rhoClamped <= right.rho()) {
                double span = right.rho() - left.rho();
                if (span <= 0.0) {
                    throw new IllegalArgumentException("Gain nodes must have strictly increasing rho.");
                }
                double alpha = (rhoClamped - left.rho()) / span;
                return new GainNode(
                    rhoClamped,
                    lerp(left.k1(), right.k1(), alpha),
                    lerp(left.k2(), right.k2(), alpha),
                    lerp(left.kr(), right.kr(), alpha)
                );
            }
        }
        return nodes.get(nodes.size() - 1);
    }

    private static double lerp(double x0, double x1, double alpha) {
        return (1.0 - alpha) * x0 + alpha * x1;
    }

    private static double referenceCommand(double t) {
        double step = t >= 8.0 ? 0.20 : 0.0;
        return 0.40 * Math.sin(0.30 * t) + step;
    }

    private static double rawSchedule(double t) {
        double value = 0.50 + 0.40 * Math.sin(0.12 * t) + 0.08 * Math.sin(0.90 * t);
        return clamp(value, 0.0, 1.0);
    }

    private static double scheduleRate(
        double rhoRaw,
        double rhoFiltered,
        double tau,
        double rateLimit
    ) {
        double requestedRate = (rhoRaw - rhoFiltered) / tau;
        return clamp(requestedRate, -rateLimit, rateLimit);
    }

    private static State derivative(State state, double control, double rho) {
        double[] coefficients = plantCoefficients(rho);
        double a = coefficients[0];
        double b = coefficients[1];
        double g = coefficients[2];
        double nonlinearTerm = 0.08 * state.x1() * state.x1() * state.x1();
        return new State(
            state.x2(),
            -a * state.x2() - b * state.x1() - nonlinearTerm + g * control
        );
    }

    private static State addScaled(State state, State increment, double scale) {
        return new State(
            state.x1() + scale * increment.x1(),
            state.x2() + scale * increment.x2()
        );
    }

    private static State rk4Step(State state, double control, double rho, double dt) {
        State k1 = derivative(state, control, rho);
        State k2 = derivative(addScaled(state, k1, 0.5 * dt), control, rho);
        State k3 = derivative(addScaled(state, k2, 0.5 * dt), control, rho);
        State k4 = derivative(addScaled(state, k3, dt), control, rho);

        return new State(
            state.x1() + dt * (k1.x1() + 2.0 * k2.x1() + 2.0 * k3.x1() + k4.x1()) / 6.0,
            state.x2() + dt * (k1.x2() + 2.0 * k2.x2() + 2.0 * k3.x2() + k4.x2()) / 6.0
        );
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        List<GainNode> nodes = new ArrayList<>();
        nodes.add(designNode(0.0));
        nodes.add(designNode(0.5));
        nodes.add(designNode(1.0));

        for (GainNode node : nodes) {
            System.out.printf(
                "node rho=%.2f: k1=%.6f, k2=%.6f, kr=%.6f%n",
                node.rho(), node.k1(), node.k2(), node.kr()
            );
        }

        final double dt = 0.002;
        final double finalTime = 40.0;
        final double tauSchedule = 0.25;
        final double scheduleRateLimit = 0.90;
        final double controlLimit = 8.0;

        State state = new State(0.0, 0.0);
        double rhoFiltered = rawSchedule(0.0);
        double sumSquaredError = 0.0;
        double peakError = 0.0;
        double peakControl = 0.0;
        long sampleCount = 0;

        Path output = Path.of("Chapter20_Lesson5_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("time,reference,position,rate,tracking_error,control,rho_raw,rho_filtered,k1,k2,kr");
            writer.newLine();

            int steps = (int) Math.round(finalTime / dt);
            for (int index = 0; index <= steps; index++) {
                double t = index * dt;
                double rhoRaw = rawSchedule(t);
                double rhoDot = scheduleRate(rhoRaw, rhoFiltered, tauSchedule, scheduleRateLimit);
                rhoFiltered = clamp(rhoFiltered + dt * rhoDot, 0.0, 1.0);

                GainNode gains = interpolate(nodes, rhoFiltered);
                double reference = referenceCommand(t);
                double controlUnsaturated =
                    -gains.k1() * state.x1() - gains.k2() * state.x2() + gains.kr() * reference;
                double control = clamp(controlUnsaturated, -controlLimit, controlLimit);
                double error = reference - state.x1();

                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                    t, reference, state.x1(), state.x2(), error, control, rhoRaw,
                    rhoFiltered, gains.k1(), gains.k2(), gains.kr()
                ));
                writer.newLine();

                sumSquaredError += error * error;
                peakError = Math.max(peakError, Math.abs(error));
                peakControl = Math.max(peakControl, Math.abs(control));
                sampleCount++;

                state = rk4Step(state, control, rhoRaw, dt);
            }
        } catch (IOException exception) {
            System.err.println("Could not write CSV: " + exception.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(sumSquaredError / sampleCount);
        System.out.printf("RMS tracking error: %.6f%n", rmsError);
        System.out.printf("Peak absolute error: %.6f%n", peakError);
        System.out.printf("Peak absolute control: %.6f%n", peakControl);
        System.out.println("Wrote " + output);
    }
}

MATLAB supports offline node design and analysis through Control System Toolbox. Robust Control Toolbox provides LPV and gain-scheduled synthesis and analysis functions in supported releases. Simulink implements the online map with Prelookup and Interpolation Using Prelookup blocks, Lookup Table blocks, Rate Limiter, low-pass filtering, Saturation, and Stateflow fallback logic. Aerospace Blockset and Robotics System Toolbox can supply higher-fidelity plant and sensor models.

In a Simulink realization, keep the gain table and trim table in data dictionaries, explicitly define extrapolation behavior, and place schedule validation before lookup execution. Log both raw and conditioned scheduling signals so that transition behavior can be audited.

Chapter20_Lesson5.m

% Chapter20_Lesson5.m
% Practical gain-scheduled control of a second-order aircraft/robotics axis.

clear; clc; close all;

rhoNodes = [0.0, 0.5, 1.0];
gainTable = zeros(numel(rhoNodes), 3);
for i = 1:numel(rhoNodes)
    gainTable(i, :) = designNode(rhoNodes(i));
    fprintf('node rho=%.2f: k1=%.6f, k2=%.6f, kr=%.6f\n', ...
        rhoNodes(i), gainTable(i, 1), gainTable(i, 2), gainTable(i, 3));
end

dt = 0.002;
finalTime = 40.0;
t = (0:dt:finalTime)';
numberOfSamples = numel(t);

tauSchedule = 0.25;
scheduleRateLimit = 0.90;
controlLimit = 8.0;

x = zeros(numberOfSamples, 2);
r = zeros(numberOfSamples, 1);
u = zeros(numberOfSamples, 1);
rhoRaw = zeros(numberOfSamples, 1);
rhoFiltered = zeros(numberOfSamples, 1);
gains = zeros(numberOfSamples, 3);

rhoFiltered(1) = rawSchedule(0.0);

for k = 1:numberOfSamples
    rhoRaw(k) = rawSchedule(t(k));
    if k > 1
        requestedRate = (rhoRaw(k) - rhoFiltered(k - 1)) / tauSchedule;
        rhoDot = clampValue(requestedRate, -scheduleRateLimit, scheduleRateLimit);
        rhoFiltered(k) = clampValue(rhoFiltered(k - 1) + dt * rhoDot, 0.0, 1.0);
    end

    gains(k, :) = interpolateGains(rhoNodes, gainTable, rhoFiltered(k));
    r(k) = referenceCommand(t(k));
    uUnsaturated = -gains(k, 1) * x(k, 1) ...
                   -gains(k, 2) * x(k, 2) ...
                   +gains(k, 3) * r(k);
    u(k) = clampValue(uUnsaturated, -controlLimit, controlLimit);

    if k < numberOfSamples
        x(k + 1, :) = rk4Step(x(k, :), u(k), rhoRaw(k), dt);
    end
end

trackingError = r - x(:, 1);
rmsError = sqrt(mean(trackingError .^ 2));
peakError = max(abs(trackingError));
peakControl = max(abs(u));

fprintf('RMS tracking error: %.6f\n', rmsError);
fprintf('Peak absolute error: %.6f\n', peakError);
fprintf('Peak absolute control: %.6f\n', peakControl);

results = table(t, r, x(:, 1), x(:, 2), trackingError, u, rhoRaw, rhoFiltered, ...
    gains(:, 1), gains(:, 2), gains(:, 3), ...
    'VariableNames', {'time', 'reference', 'position', 'rate', 'tracking_error', ...
    'control', 'rho_raw', 'rho_filtered', 'k1', 'k2', 'kr'});
writetable(results, 'Chapter20_Lesson5_results.csv');

figure;
plot(t, r, '--', t, x(:, 1), '-');
grid on;
xlabel('Time (s)');
ylabel('Command and response');
legend('Reference', 'Response', 'Location', 'best');
title('Gain-Scheduled Tracking');

figure;
plot(t, rhoRaw, '--', t, rhoFiltered, '-');
grid on;
xlabel('Time (s)');
ylabel('Scheduling variable');
legend('Raw schedule', 'Filtered schedule', 'Location', 'best');
title('Schedule Filtering and Rate Limiting');

figure;
plot(t, gains(:, 1), t, gains(:, 2), t, gains(:, 3));
grid on;
xlabel('Time (s)');
ylabel('Interpolated gains');
legend('k_1', 'k_2', 'k_r', 'Location', 'best');
title('Scheduled Controller Gains');

function gains = designNode(rho)
    [a, b, g] = plantCoefficients(rho);
    zeta = 0.85;
    omegaN = 2.2 - 0.4 * rho;
    k1 = (omegaN^2 - b) / g;
    k2 = (2.0 * zeta * omegaN - a) / g;
    kr = omegaN^2 / g;
    gains = [k1, k2, kr];
end

function gains = interpolateGains(rhoNodes, gainTable, rho)
    rhoClamped = clampValue(rho, rhoNodes(1), rhoNodes(end));
    gains = zeros(1, 3);
    for column = 1:3
        gains(column) = interp1(rhoNodes, gainTable(:, column), rhoClamped, 'linear');
    end
end

function command = referenceCommand(t)
    step = 0.0;
    if t >= 8.0
        step = 0.20;
    end
    command = 0.40 * sin(0.30 * t) + step;
end

function rho = rawSchedule(t)
    value = 0.50 + 0.40 * sin(0.12 * t) + 0.08 * sin(0.90 * t);
    rho = clampValue(value, 0.0, 1.0);
end

function [a, b, g] = plantCoefficients(rho)
    a = 0.8 + 0.7 * rho;
    b = 1.4 - 0.5 * rho;
    g = 1.0 / (1.0 + 0.8 * rho);
end

function dx = stateDerivative(x, control, rho)
    [a, b, g] = plantCoefficients(rho);
    nonlinearTerm = 0.08 * x(1)^3;
    dx = [x(2), -a * x(2) - b * x(1) - nonlinearTerm + g * control];
end

function nextState = rk4Step(state, control, rho, dt)
    k1 = stateDerivative(state, control, rho);
    k2 = stateDerivative(state + 0.5 * dt * k1, control, rho);
    k3 = stateDerivative(state + 0.5 * dt * k2, control, rho);
    k4 = stateDerivative(state + dt * k3, control, rho);
    nextState = state + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
end

function value = clampValue(value, lower, upper)
    value = max(lower, min(upper, value));
end

Wolfram Language functions such as StateSpaceModel, ControllableModelQ, LQRegulatorGains, NDSolveValue, and Linearize are useful for local-model analysis and simulation. The notebook source below implements the same gain table and sampled RK4 simulation used in the other languages.

Chapter20_Lesson5.nb

Notebook[{
 Cell["Chapter20_Lesson5.nb", "Title"],
 Cell["Practical gain-scheduled control of a second-order aircraft/robotics axis", "Subtitle"],
 Cell[BoxData["ClearAll[\"Global`*\"];"] , "Input"],
 Cell[BoxData["clamp[value_, lower_, upper_] := Max[lower, Min[upper, value]];"] , "Input"],
 Cell[BoxData["plantCoefficients[rho_] := {0.8 + 0.7 rho, 1.4 - 0.5 rho, 1/(1 + 0.8 rho)};"] , "Input"],
 Cell[BoxData["desiredDynamics[rho_] := {0.85, 2.2 - 0.4 rho};"] , "Input"],
 Cell[BoxData["designNode[rho_] := Module[{a, b, g, zeta, omegaN}, {a, b, g} = plantCoefficients[rho]; {zeta, omegaN} = desiredDynamics[rho]; {rho, (omegaN^2 - b)/g, (2 zeta omegaN - a)/g, omegaN^2/g}];"] , "Input"],
 Cell[BoxData["nodes = designNode /@ {0.0, 0.5, 1.0}; nodes // TableForm"] , "Input"],
 Cell[BoxData["interpolate[nodes_, rho_] := Module[{rhoC, pair, alpha, left, right}, rhoC = clamp[rho, nodes[[1, 1]], nodes[[-1, 1]]]; pair = SelectFirst[Partition[nodes, 2, 1], rhoC <= #[[2, 1]] &]; left = pair[[1]]; right = pair[[2]]; alpha = (rhoC - left[[1]])/(right[[1]] - left[[1]]); {rhoC, Sequence @@ ((1 - alpha) left[[2 ;; 4]] + alpha right[[2 ;; 4]])}];"] , "Input"],
 Cell[BoxData["referenceCommand[t_] := 0.40 Sin[0.30 t] + If[t >= 8.0, 0.20, 0.0];"] , "Input"],
 Cell[BoxData["rawSchedule[t_] := clamp[0.50 + 0.40 Sin[0.12 t] + 0.08 Sin[0.90 t], 0.0, 1.0];"] , "Input"],
 Cell[BoxData["stateDerivative[{x1_, x2_}, control_, rho_] := Module[{a, b, g}, {a, b, g} = plantCoefficients[rho]; {x2, -a x2 - b x1 - 0.08 x1^3 + g control}];"] , "Input"],
 Cell[BoxData["rk4Step[state_, control_, rho_, dt_] := Module[{k1, k2, k3, k4}, k1 = stateDerivative[state, control, rho]; k2 = stateDerivative[state + 0.5 dt k1, control, rho]; k3 = stateDerivative[state + 0.5 dt k2, control, rho]; k4 = stateDerivative[state + dt k3, control, rho]; state + dt (k1 + 2 k2 + 2 k3 + k4)/6];"] , "Input"],
 Cell[BoxData["dt = 0.002; finalTime = 40.0; tauSchedule = 0.25; scheduleRateLimit = 0.90; controlLimit = 8.0;"] , "Input"],
 Cell[BoxData["simulation = Reap[Module[{state = {0.0, 0.0}, rhoFiltered = rawSchedule[0.0], rhoRaw, rhoDot, gains, reference, control, error}, Do[rhoRaw = rawSchedule[t]; rhoDot = clamp[(rhoRaw - rhoFiltered)/tauSchedule, -scheduleRateLimit, scheduleRateLimit]; rhoFiltered = clamp[rhoFiltered + dt rhoDot, 0.0, 1.0]; gains = interpolate[nodes, rhoFiltered]; reference = referenceCommand[t]; control = clamp[-gains[[2]] state[[1]] - gains[[3]] state[[2]] + gains[[4]] reference, -controlLimit, controlLimit]; error = reference - state[[1]]; Sow[{t, reference, state[[1]], state[[2]], error, control, rhoRaw, rhoFiltered, gains[[2]], gains[[3]], gains[[4]]}]; state = rk4Step[state, control, rhoRaw, dt], {t, 0.0, finalTime, dt}]]][[2, 1]];"] , "Input"],
 Cell[BoxData["headers = {\"time\", \"reference\", \"position\", \"rate\", \"tracking_error\", \"control\", \"rho_raw\", \"rho_filtered\", \"k1\", \"k2\", \"kr\"}; Export[\"Chapter20_Lesson5_results.csv\", Prepend[simulation, headers]];"] , "Input"],
 Cell[BoxData["rmsError = Sqrt[Mean[simulation[[All, 5]]^2]]; peakError = Max[Abs[simulation[[All, 5]]]]; peakControl = Max[Abs[simulation[[All, 6]]]]; {rmsError, peakError, peakControl}"] , "Input"],
 Cell[BoxData["ListLinePlot[{simulation[[All, {1, 2}]], simulation[[All, {1, 3}]]}, PlotLegends -> {\"Reference\", \"Response\"}, Frame -> True, FrameLabel -> {\"Time (s)\", \"Command and response\"}, PlotLabel -> \"Gain-Scheduled Tracking\"]"] , "Input"],
 Cell[BoxData["ListLinePlot[{simulation[[All, {1, 7}]], simulation[[All, {1, 8}]]}, PlotLegends -> {\"Raw schedule\", \"Filtered schedule\"}, Frame -> True, FrameLabel -> {\"Time (s)\", \"Scheduling variable\"}, PlotLabel -> \"Schedule Filtering and Rate Limiting\"]"] , "Input"]
}, WindowTitle -> "Chapter20_Lesson5"]

Problem 1 (Local Gain Derivation): For \( \ddot{x}+a\dot{x}+bx=gu \) with \( g\ne 0 \), derive the scheduled gains in \( u=-k_1x-k_2\dot{x}+k_rr \) that produce the desired characteristic polynomial \( s^2+2\zeta\omega_ns+\omega_n^2 \) and unit steady-state reference gain.

Solution: Substitution of the control law gives

\[ \ddot{x}+(a+gk_2)\dot{x}+(b+gk_1)x=gk_rr. \]

Coefficient matching requires

\[ a+gk_2=2\zeta\omega_n, \qquad b+gk_1=\omega_n^2. \]

Therefore

\[ k_1=\frac{\omega_n^2-b}{g}, \qquad k_2=\frac{2\zeta\omega_n-a}{g}. \]

For a constant reference, set derivatives to zero and require \( x_{ss}=r \):

\[ (b+gk_1)r=gk_rr \quad\Longrightarrow\quad k_r=\frac{b+gk_1}{g}=\frac{\omega_n^2}{g}. \]

Problem 2 (Common-Lyapunov Proof): Let \( \mathbf{A}_{cl}(\rho)=\sum_i\alpha_i(\rho)\mathbf{A}_{cl,i} \) with nonnegative weights summing to one. Prove uniform exponential stability if a common \( \mathbf{P}\succ 0 \) satisfies \( \mathbf{A}_{cl,i}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{cl,i} \preceq-\mathbf{Q} \) for all vertices and some \( \mathbf{Q}\succ 0 \).

Solution: With \( V=\mathbf{x}^{T}\mathbf{P}\mathbf{x} \),

\[ \dot{V}\le-\lambda_{\min}(\mathbf{Q})\|\mathbf{x}\|^2. \]

Also,

\[ \lambda_{\min}(\mathbf{P})\|\mathbf{x}\|^2 \le V \le\lambda_{\max}(\mathbf{P})\|\mathbf{x}\|^2. \]

Hence

\[ \dot{V}\le- \frac{\lambda_{\min}(\mathbf{Q})}{\lambda_{\max}(\mathbf{P})}V. \]

Integration yields

\[ V(t)\le V(0) \exp\left[-\frac{\lambda_{\min}(\mathbf{Q})} {\lambda_{\max}(\mathbf{P})}t\right], \]

which proves a uniform exponential bound independent of the scheduling trajectory.

Problem 3 (Interpolation Cross Terms): For two-node affine interpolation of \( \mathbf{B} \) and \( \mathbf{K} \), compute the difference between \( \mathbf{B}(\alpha)\mathbf{K}(\alpha) \) and the endpoint interpolation \( (1-\alpha)\mathbf{B}_1\mathbf{K}_1+ \alpha\mathbf{B}_2\mathbf{K}_2 \).

Solution: Direct expansion gives

\[ \begin{aligned} \Delta(\alpha) ={}&\mathbf{B}(\alpha)\mathbf{K}(\alpha) -\left[(1-\alpha)\mathbf{B}_1\mathbf{K}_1 +\alpha\mathbf{B}_2\mathbf{K}_2\right] \\ ={}&\alpha(1-\alpha) \left(\mathbf{B}_1-\mathbf{B}_2\right) \left(\mathbf{K}_2-\mathbf{K}_1\right). \end{aligned} \]

The difference vanishes at the endpoints, but can be largest near the cell midpoint. It also vanishes if either \( \mathbf{B} \) or \( \mathbf{K} \) is constant across the cell.

Problem 4 (Rate-Limited Schedule Filter): Consider \( \dot{\rho}_f=\operatorname{sat}_{[-\nu,\nu]} ((\rho-\rho_f)/\tau) \). Show that the implementation guarantees \( |\dot{\rho}_f|\le\nu \), and determine the unsaturated small-signal transfer function from \( \rho \) to \( \rho_f \).

Solution: The saturation definition directly gives

\[ -\nu\le\dot{\rho}_f\le\nu. \]

When the limiter is inactive,

\[ \tau\dot{\rho}_f+\rho_f=\rho. \]

Taking Laplace transforms with zero initial condition gives

\[ \frac{\rho_f(s)}{\rho(s)}=\frac{1}{\tau s+1}. \]

Thus the filter attenuates high-frequency schedule noise while the limiter enforces the scheduling-rate assumption used in analysis.

Problem 5 (Grid-Spacing Error Bound): Let a scalar gain \( k(\rho) \) be twice continuously differentiable on a grid cell of width \( h \), and suppose \( |k''(\rho)|\le M \). Bound the error of linear interpolation.

Solution: The interpolation remainder is

\[ k(\rho)-k_{lin}(\rho)= \frac{k''(\xi)}{2}(\rho-\rho_i)(\rho-\rho_{i+1}) \]

for some \( \xi \) in the cell. The product magnitude is maximized at the midpoint and equals \( h^2/4 \). Therefore

\[ |k(\rho)-k_{lin}(\rho)|\le\frac{Mh^2}{8}. \]

To enforce an interpolation tolerance \( \varepsilon_k \), a sufficient grid-width condition is

\[ h\le\sqrt{\frac{8\varepsilon_k}{M}}. \]

Practical gain scheduling requires a certified envelope, compatible trim and linearization data, coordinated local controllers, carefully conditioned scheduling signals, stability-aware interpolation, and explicit safety logic. Stable local controllers do not by themselves guarantee a stable scheduled system. Common or parameter-dependent Lyapunov analysis, dense off-grid checks, scheduling-rate tests, nonlinear uncertainty simulation, and code/hardware verification together form the required evidence. The five implementations demonstrate the same portable architecture for an aircraft or robotic axis.