In gain scheduling, a measurable operating variable \( \rho(t) \in \mathcal{P} \subset \mathbb{R}^{n_\rho} \) selects or interpolates controller coefficients. A continuous-time LPV-like plant can be represented as

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

For scheduled state feedback, \( u(t)=-K(\rho(t))x(t) \), the closed-loop state matrix is

\[ A_{\mathrm{cl} }(\rho)=A(\rho)-B(\rho)K(\rho), \qquad \dot{x}=A_{\mathrm{cl} }(\rho(t))x. \]

A local design typically verifies that every frozen matrix \( A_{\mathrm{cl} }(\bar{\rho}) \) is Hurwitz for each fixed \( \bar{\rho}\in\mathcal{P} \). This is a necessary engineering check, but it is not a complete time-domain stability proof because the actual system is linear time-varying: \( \rho=\rho(t) \).

Throughout the lesson, assume that the scheduling signal remains inside a compact admissible set and is measurable:

\[ \rho(t)\in\mathcal{P}, \qquad \|\dot{\rho}(t)\|\leq \nu, \]

whenever a rate bound is required. A common quadratic certificate will later remove the need for the second restriction.

Suppose that every frozen matrix is Hurwitz:

\[ \max_i \operatorname{Re}\lambda_i \bigl(A_{\mathrm{cl} }(\rho)\bigr)<0 \qquad \forall \rho\in\mathcal{P}. \]

This condition only describes systems in which \( \rho(t)\equiv\bar{\rho} \). For a varying schedule, the state-transition matrix depends on the ordered product of infinitesimal transitions. In general,

\[ A_{\mathrm{cl} }(\rho_1)A_{\mathrm{cl} }(\rho_2) \neq A_{\mathrm{cl} }(\rho_2)A_{\mathrm{cl} }(\rho_1). \]

Therefore, the decay directions of one frozen system can be amplified when the schedule moves into another system. Stable eigenvalues at each point do not control the interaction among changing eigenvectors, controller states, interpolation rules, and scheduling dynamics.

The correct question is not merely whether every frozen system is stable, but whether there exists a Lyapunov function that decreases along every admissible scheduling trajectory.

flowchart TD
  A["Local controllers are stable at design points"] --> B["Interpolate or schedule controller coefficients"]
  B --> C["Form the actual time-varying closed loop"]
  C --> D{"Which certificate \nis available?"}
  D --> E["Common Lyapunov matrix P"]
  D --> F["Parameter-dependent P(rho) \nand rate bounds"]
  D --> G["Only frozen poles"]
  E --> H["Uniform stability for \nevery admissible schedule"]
  F --> I["Stability for schedules \nsatisfying rate limits"]
  G --> J["No global guarantee; \nperform redesign or \nstronger analysis"]
        

The strongest elementary certificate used in this lesson is a common quadratic Lyapunov function (CQLF):

\[ V(x)=x^\top P x, \qquad P=P^\top>0. \]

If one matrix \( P \) satisfies

\[ A_{\mathrm{cl} }^\top(\rho)P+ P A_{\mathrm{cl} }(\rho) \leq -Q \qquad \forall \rho\in\mathcal{P}, \]

for some constant \( Q=Q^\top>0 \), then the scheduled system is uniformly exponentially stable for every measurable scheduling trajectory that remains in \( \mathcal{P} \). No bound on \( \dot{\rho} \) is needed.

Theorem 1: Uniform Exponential Stability under a CQLF

Let \( P>0 \) and \( Q>0 \) satisfy the inequality above. Then

\[ \|x(t)\| \leq \sqrt{\frac{\lambda_{\max}(P)}{\lambda_{\min}(P)} } \exp\!\left( -\frac{\lambda_{\min}(Q)} {2\lambda_{\max}(P)}t \right) \|x(0)\|. \]

Proof

Along trajectories,

\[ \dot{V} = x^\top \left( A_{\mathrm{cl} }^\top(\rho)P+ P A_{\mathrm{cl} }(\rho) \right)x \leq -x^\top Qx. \]

Using eigenvalue bounds,

\[ \lambda_{\min}(P)\|x\|^2 \leq V(x) \leq \lambda_{\max}(P)\|x\|^2, \]

\[ \dot{V} \leq -\lambda_{\min}(Q)\|x\|^2 \leq -\frac{\lambda_{\min}(Q)} {\lambda_{\max}(P)}V. \]

Integration gives \( V(t)\leq e^{-ct}V(0) \), where \( c=\lambda_{\min}(Q)/\lambda_{\max}(P) \). Converting the Lyapunov bound back to the Euclidean norm proves the result.

Suppose the closed-loop matrix belongs to the convex hull of \( N \) vertex matrices:

\[ A_{\mathrm{cl} }(\rho) = \sum_{i=1}^{N}\alpha_i(\rho)A_i, \qquad \alpha_i(\rho)\geq 0, \qquad \sum_{i=1}^{N}\alpha_i(\rho)=1. \]

For a common matrix \( P \), the Lyapunov expression is affine in the convex weights:

\[ A_{\mathrm{cl} }^\top(\rho)P+ P A_{\mathrm{cl} }(\rho) = \sum_{i=1}^{N}\alpha_i(\rho) \left(A_i^\top P+P A_i\right). \]

Hence it is sufficient to solve the finite family of LMIs

\[ P>0, \qquad A_i^\top P+P A_i<0, \qquad i=1,\ldots,N. \]

Every convex combination is then negative definite. This vertex property is one reason polytopic LPV embeddings are useful. However, the test can be conservative: a stable family may not possess one common quadratic matrix even though another Lyapunov function exists.

To reduce CQLF conservatism, use

\[ V(x,\rho)=x^\top P(\rho)x, \qquad P(\rho)=P^\top(\rho)>0. \]

Differentiation introduces an additional term:

\[ \dot{V} = x^\top \left[ A_{\mathrm{cl} }^\top(\rho)P(\rho) +P(\rho)A_{\mathrm{cl} }(\rho) +\dot{P}(\rho) \right]x, \]

\[ \dot{P}(\rho) = \sum_{j=1}^{n_\rho} \frac{\partial P}{\partial \rho_j}\dot{\rho}_j. \]

This term is the mathematical reason that a schedule-rate limit can be essential. Assume

\[ A_{\mathrm{cl} }^\top P+ P A_{\mathrm{cl} } \leq -Q(\rho), \qquad \lambda_{\min}(Q(\rho))\geq q>0, \]

\[ \left\| \frac{\partial P}{\partial \rho_j} \right\| \leq \mu_j, \qquad |\dot{\rho}_j|\leq \nu_j. \]

Then

\[ \dot{V} \leq -\left( q-\sum_{j=1}^{n_\rho}\mu_j\nu_j \right)\|x\|^2. \]

A sufficient condition for uniform exponential stability is therefore

\[ \sum_{j=1}^{n_\rho}\mu_j\nu_j<q. \]

Unlike a CQLF certificate, a parameter-dependent certificate can be explicitly tied to the maximum admissible scheduling speed.

Another classical route starts with a uniformly Hurwitz frozen family. For every fixed \( \rho \), define \( P(\rho) \) as the solution of

\[ A_{\mathrm{cl} }^\top(\rho)P(\rho) +P(\rho)A_{\mathrm{cl} }(\rho) =-I. \]

If the family is smooth, uniformly separated from instability, and \( P(\rho) \) has bounded derivatives on the compact set \( \mathcal{P} \), then a sufficiently small scheduling rate guarantees stability. For one scalar schedule,

\[ \dot{V} \leq -\|x\|^2 + \left\| \frac{\partial P}{\partial \rho} \right\| |\dot{\rho}| \|x\|^2. \]

Thus any rate satisfying

\[ |\dot{\rho}|< \left( \sup_{\rho\in\mathcal{P} } \left\| \frac{\partial P}{\partial \rho} \right\| \right)^{-1} \]

is sufficient for this normalized equation. The bound is often conservative, but it exposes the distinction between frozen stability and trajectory stability.

Interpolating gains is not automatically equivalent to interpolating stable closed-loop systems. For static state feedback, one can directly analyze \( A(\rho)-B(\rho)K(\rho) \). For a dynamic controller,

\[ \dot{x}_c=A_c(\rho)x_c+B_c(\rho)e, \qquad u=C_c(\rho)x_c+D_c(\rho)e, \]

changing controller realizations may introduce state-coordinate inconsistencies. If a scheduled coordinate transformation is \( x_c=T(\rho)z_c \), then

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

The second term is a hidden coupling term. Ignoring it can make a realization look correct at every frozen point while producing undesirable transients or instability during scheduling.

Practical safeguards include a consistent controller realization, interpolation methods that preserve a common certificate, bumpless transfer, state initialization rules, rate limiting, saturation handling, and explicit analysis of the augmented plant-controller state.

Consider the parameter-varying plant

\[ A(\rho)= \begin{bmatrix} 0 & 1\\ -(1+0.8\rho) & -(0.4+0.2\rho) \end{bmatrix}, \qquad B= \begin{bmatrix} 0\\1 \end{bmatrix}, \qquad 0\leq\rho\leq 1. \]

Interpolate the state-feedback gains as

\[ K(\rho)= \begin{bmatrix} 3+\rho & 2.5+0.5\rho \end{bmatrix}. \]

The closed-loop matrix is

\[ A_{\mathrm{cl} }(\rho)= \begin{bmatrix} 0 & 1\\ -(4+1.8\rho) & -(2.9+0.7\rho) \end{bmatrix}. \]

At the midpoint, solve

\[ A_{\mathrm{cl} }^\top(0.5)P+ P A_{\mathrm{cl} }(0.5)=-I. \]

This gives approximately

\[ P= \begin{bmatrix} 1.23932496 & 0.10204082\\ 0.10204082 & 0.18524333 \end{bmatrix}, \qquad \lambda_{\min}(P)\approx 0.175456. \]

Defining \( Q(\rho)=- [A_{\mathrm{cl} }^\top(\rho)P+P A_{\mathrm{cl} }(\rho)] \), the minimum eigenvalue over the complete interval is approximately \( 0.639102 \). Since \( A_{\mathrm{cl} }(\rho) \) is affine in \( \rho \), checking the two endpoints is already sufficient for the CQLF inequality; the dense grid used in the code is an additional numerical verification.

The schedule used in simulation is

\[ \rho(t)= \operatorname{sat}_{[0,1]} \left(0.5+0.5\sin(0.8t)\right). \]

Because the certificate uses a constant \( P \), its proof is independent of \( \dot{\rho} \). The scheduling frequency can be changed without invalidating the mathematical certificate, provided the model and schedule remain inside the certified set.

flowchart TD
  R["Measured schedule rho(t)"] --> G["Interpolate gain K(rho)"]
  X["Measured state x(t)"] --> C["Compute u = -K(rho) x"]
  G --> C
  C --> P["Parameter-varying plant"]
  P --> X
  X --> V["Evaluate V = xT P x"]
  R --> Q["Evaluate Q(rho)"]
  Q --> M{"Margins positive?"}
  V --> M
  M -->|yes| S["Continue scheduled operation"]
  M -->|no| F["Limit schedule or enter safe fallback"]
        

The Python program uses NumPy, SciPy, and Matplotlib. It computes \( P \) from the midpoint Lyapunov equation, verifies all frozen poles, checks the Lyapunov margin on a dense grid, integrates the time-varying system using RK4, and exports numerical results.

Chapter20_Lesson3.py

#!/usr/bin/env python3
"""
Chapter20_Lesson3.py
Stability verification and simulation for a gain-scheduled controller.

Dependencies:
    pip install numpy scipy matplotlib
"""

from __future__ import annotations

import csv
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import eigvals, solve_continuous_lyapunov


P_REFERENCE = np.array(
    [
        [1.23932496, 0.10204082],
        [0.10204082, 0.18524333],
    ],
    dtype=float,
)


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


def scheduling_variable(t: float, omega: float = 0.8) -> float:
    """Measured scheduling variable constrained to rho in [0, 1]."""
    return clamp(0.5 + 0.5 * np.sin(omega * t), 0.0, 1.0)


def plant_matrix(rho: float) -> np.ndarray:
    return np.array(
        [
            [0.0, 1.0],
            [-(1.0 + 0.8 * rho), -(0.4 + 0.2 * rho)],
        ],
        dtype=float,
    )


def input_matrix() -> np.ndarray:
    return np.array([[0.0], [1.0]], dtype=float)


def scheduled_gain(rho: float) -> np.ndarray:
    """Linear interpolation between K(0)=[3,2.5] and K(1)=[4,3]."""
    return np.array([[3.0 + rho, 2.5 + 0.5 * rho]], dtype=float)


def closed_loop_matrix(rho: float) -> np.ndarray:
    return plant_matrix(rho) - input_matrix() @ scheduled_gain(rho)


def common_lyapunov_matrix() -> np.ndarray:
    """
    Compute P at rho=0.5 from A_mid^T P + P A_mid = -I.

    The returned P is then checked over the complete scheduling interval.
    """
    a_mid = closed_loop_matrix(0.5)
    return solve_continuous_lyapunov(a_mid.T, -np.eye(2))


def symmetric_min_eigenvalue(matrix: np.ndarray) -> float:
    return float(np.min(np.linalg.eigvalsh(0.5 * (matrix + matrix.T))))


def verify_certificate(samples: int = 1001) -> tuple[np.ndarray, float, float]:
    p = common_lyapunov_matrix()
    p_margin = symmetric_min_eigenvalue(p)
    q_margin = float("inf")

    for rho in np.linspace(0.0, 1.0, samples):
        a_cl = closed_loop_matrix(float(rho))
        q = -(a_cl.T @ p + p @ a_cl)
        q_margin = min(q_margin, symmetric_min_eigenvalue(q))

        poles = eigvals(a_cl)
        if np.max(np.real(poles)) >= 0.0:
            raise RuntimeError(f"Frozen closed loop is unstable at rho={rho:.6f}")

    if p_margin <= 0.0 or q_margin <= 0.0:
        raise RuntimeError("The sampled common quadratic Lyapunov test failed.")

    return p, p_margin, q_margin


def derivative(t: float, x: np.ndarray, omega: float) -> np.ndarray:
    rho = scheduling_variable(t, omega)
    return closed_loop_matrix(rho) @ x


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


def simulate(
    p: np.ndarray,
    omega: float = 0.8,
    final_time: float = 20.0,
    dt: float = 0.005,
) -> dict[str, np.ndarray]:
    steps = int(round(final_time / dt))
    time = np.linspace(0.0, final_time, steps + 1)
    state = np.zeros((steps + 1, 2), dtype=float)
    rho_hist = np.zeros(steps + 1, dtype=float)
    v_hist = np.zeros(steps + 1, dtype=float)
    vdot_hist = np.zeros(steps + 1, dtype=float)

    state[0] = np.array([1.0, -0.5], dtype=float)

    for k, t in enumerate(time):
        rho = scheduling_variable(float(t), omega)
        a_cl = closed_loop_matrix(rho)
        x = state[k]

        rho_hist[k] = rho
        v_hist[k] = float(x.T @ p @ x)
        vdot_hist[k] = float(x.T @ (a_cl.T @ p + p @ a_cl) @ x)

        if k < steps:
            state[k + 1] = rk4_step(float(t), x, dt, omega)

    return {
        "time": time,
        "state": state,
        "rho": rho_hist,
        "V": v_hist,
        "Vdot": vdot_hist,
    }


def save_results(results: dict[str, np.ndarray], output_dir: Path) -> None:
    output_dir.mkdir(parents=True, exist_ok=True)
    csv_path = output_dir / "Chapter20_Lesson3_python_results.csv"

    with csv_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(["time", "rho", "x1", "x2", "V", "Vdot"])
        for row in zip(
            results["time"],
            results["rho"],
            results["state"][:, 0],
            results["state"][:, 1],
            results["V"],
            results["Vdot"],
        ):
            writer.writerow([f"{value:.12g}" for value in row])

    figure, axes = plt.subplots(3, 1, figsize=(9, 9), sharex=True)
    axes[0].plot(results["time"], results["state"][:, 0], label="x1")
    axes[0].plot(results["time"], results["state"][:, 1], label="x2")
    axes[0].set_ylabel("state")
    axes[0].legend()
    axes[0].grid(True)

    axes[1].plot(results["time"], results["rho"])
    axes[1].set_ylabel("rho(t)")
    axes[1].grid(True)

    axes[2].plot(results["time"], results["V"], label="V")
    axes[2].plot(results["time"], results["Vdot"], label="Vdot")
    axes[2].set_xlabel("time [s]")
    axes[2].set_ylabel("Lyapunov quantities")
    axes[2].legend()
    axes[2].grid(True)

    figure.tight_layout()
    figure.savefig(output_dir / "Chapter20_Lesson3_python.png", dpi=180)
    plt.close(figure)


def main() -> None:
    p, p_margin, q_margin = verify_certificate()
    reference_error = float(np.max(np.abs(p - P_REFERENCE)))

    print("Computed common Lyapunov matrix P:")
    print(p)
    print(f"lambda_min(P) = {p_margin:.9f}")
    print(f"minimum sampled lambda_min(Q(rho)) = {q_margin:.9f}")
    print(f"maximum difference from documented P = {reference_error:.3e}")

    for rho in (0.0, 0.5, 1.0):
        poles = eigvals(closed_loop_matrix(rho))
        print(f"rho={rho:.1f}, frozen poles={poles}")

    results = simulate(p=p, omega=0.8)
    save_results(results, Path("."))

    final_norm = float(np.linalg.norm(results["state"][-1]))
    max_vdot = float(np.max(results["Vdot"]))
    print(f"final state norm = {final_norm:.6e}")
    print(f"maximum analytical Vdot along trajectory = {max_vdot:.6e}")

    if max_vdot > 1e-10:
        raise RuntimeError("Vdot became positive; certificate or implementation is inconsistent.")


if __name__ == "__main__":
    main()

The C++ implementation uses only the standard library. For the two-state example, Hurwitz stability is checked through trace and determinant conditions, while positive definiteness is checked with the closed-form eigenvalues of a symmetric two-by-two matrix.

Chapter20_Lesson3.cpp

// Chapter20_Lesson3.cpp
// Self-contained gain-scheduled stability verification and RK4 simulation.
//
// Build:
//   g++ -std=c++17 -O2 Chapter20_Lesson3.cpp -o Chapter20_Lesson3
// Run:
//   ./Chapter20_Lesson3

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

struct Vec2 {
    double x1{};
    double x2{};
};

struct Mat2 {
    double a11{};
    double a12{};
    double a21{};
    double a22{};
};

constexpr Mat2 P{
    1.23932496, 0.10204082,
    0.10204082, 0.18524333
};

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

double rho_of_t(double t, double omega = 0.8) {
    return clamp(0.5 + 0.5 * std::sin(omega * t), 0.0, 1.0);
}

Mat2 closed_loop_matrix(double rho) {
    // A(rho) - B K(rho), where
    // A = [[0,1],[-(1+0.8rho),-(0.4+0.2rho)]],
    // B = [[0],[1]], K = [3+rho, 2.5+0.5rho].
    return Mat2{
        0.0,
        1.0,
        -(4.0 + 1.8 * rho),
        -(2.9 + 0.7 * rho)
    };
}

Mat2 transpose(const Mat2& m) {
    return Mat2{m.a11, m.a21, m.a12, m.a22};
}

Mat2 multiply(const Mat2& lhs, const Mat2& rhs) {
    return Mat2{
        lhs.a11 * rhs.a11 + lhs.a12 * rhs.a21,
        lhs.a11 * rhs.a12 + lhs.a12 * rhs.a22,
        lhs.a21 * rhs.a11 + lhs.a22 * rhs.a21,
        lhs.a21 * rhs.a12 + lhs.a22 * rhs.a22
    };
}

Mat2 add(const Mat2& lhs, const Mat2& rhs) {
    return Mat2{
        lhs.a11 + rhs.a11,
        lhs.a12 + rhs.a12,
        lhs.a21 + rhs.a21,
        lhs.a22 + rhs.a22
    };
}

Mat2 negate(const Mat2& m) {
    return Mat2{-m.a11, -m.a12, -m.a21, -m.a22};
}

Vec2 multiply(const Mat2& m, const Vec2& x) {
    return Vec2{
        m.a11 * x.x1 + m.a12 * x.x2,
        m.a21 * x.x1 + m.a22 * x.x2
    };
}

Vec2 add(const Vec2& lhs, const Vec2& rhs) {
    return Vec2{lhs.x1 + rhs.x1, lhs.x2 + rhs.x2};
}

Vec2 scale(double scalar, const Vec2& x) {
    return Vec2{scalar * x.x1, scalar * x.x2};
}

double dot(const Vec2& lhs, const Vec2& rhs) {
    return lhs.x1 * rhs.x1 + lhs.x2 * rhs.x2;
}

double norm(const Vec2& x) {
    return std::sqrt(dot(x, x));
}

double minimum_eigenvalue_symmetric(const Mat2& m) {
    const double b = 0.5 * (m.a12 + m.a21);
    const double trace = m.a11 + m.a22;
    const double discriminant =
        std::sqrt((m.a11 - m.a22) * (m.a11 - m.a22) + 4.0 * b * b);
    return 0.5 * (trace - discriminant);
}

double quadratic_form(const Vec2& x, const Mat2& m) {
    return dot(x, multiply(m, x));
}

Vec2 derivative(double t, const Vec2& x, double omega) {
    return multiply(closed_loop_matrix(rho_of_t(t, omega)), x);
}

Vec2 rk4_step(double t, const Vec2& x, double dt, double omega) {
    const Vec2 k1 = derivative(t, x, omega);
    const Vec2 k2 = derivative(t + 0.5 * dt, add(x, scale(0.5 * dt, k1)), omega);
    const Vec2 k3 = derivative(t + 0.5 * dt, add(x, scale(0.5 * dt, k2)), omega);
    const Vec2 k4 = derivative(t + dt, add(x, scale(dt, k3)), omega);

    return add(
        x,
        scale(
            dt / 6.0,
            add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))
        )
    );
}

void check_frozen_hurwitz(const Mat2& a, double rho) {
    const double trace = a.a11 + a.a22;
    const double determinant = a.a11 * a.a22 - a.a12 * a.a21;

    // For a real 2x2 matrix, trace < 0 and determinant > 0 imply Hurwitz.
    if (!(trace < 0.0 && determinant > 0.0)) {
        throw std::runtime_error("Frozen matrix is not Hurwitz at rho=" + std::to_string(rho));
    }
}

int main() {
    try {
        const double p_margin = minimum_eigenvalue_symmetric(P);
        double q_margin = 1.0e100;

        for (int i = 0; i <= 1000; ++i) {
            const double rho = static_cast<double>(i) / 1000.0;
            const Mat2 a = closed_loop_matrix(rho);
            check_frozen_hurwitz(a, rho);

            const Mat2 q = negate(add(multiply(transpose(a), P), multiply(P, a)));
            q_margin = std::min(q_margin, minimum_eigenvalue_symmetric(q));
        }

        if (!(p_margin > 0.0 && q_margin > 0.0)) {
            throw std::runtime_error("Common Lyapunov certificate failed.");
        }

        std::cout << std::setprecision(10);
        std::cout << "lambda_min(P) = " << p_margin << '\n';
        std::cout << "minimum sampled lambda_min(Q(rho)) = " << q_margin << '\n';

        const double dt = 0.005;
        const double final_time = 20.0;
        const double omega = 0.8;
        const int steps = static_cast<int>(std::round(final_time / dt));

        Vec2 x{1.0, -0.5};
        double max_vdot = -1.0e100;

        std::ofstream output("Chapter20_Lesson3_cpp_results.csv");
        if (!output) {
            throw std::runtime_error("Could not open CSV output file.");
        }
        output << "time,rho,x1,x2,V,Vdot\n";
        output << std::setprecision(12);

        for (int k = 0; k <= steps; ++k) {
            const double t = k * dt;
            const double rho = rho_of_t(t, omega);
            const Mat2 a = closed_loop_matrix(rho);
            const Mat2 s = add(multiply(transpose(a), P), multiply(P, a));
            const double v = quadratic_form(x, P);
            const double vdot = quadratic_form(x, s);
            max_vdot = std::max(max_vdot, vdot);

            output << t << ',' << rho << ',' << x.x1 << ',' << x.x2
                   << ',' << v << ',' << vdot << '\n';

            if (k < steps) {
                x = rk4_step(t, x, dt, omega);
            }
        }

        std::cout << "final state norm = " << norm(x) << '\n';
        std::cout << "maximum analytical Vdot along trajectory = " << max_vdot << '\n';

        if (max_vdot > 1.0e-10) {
            throw std::runtime_error("Vdot became positive.");
        }
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }

    return 0;
}

The Java version is self-contained and uses immutable records for two-dimensional vectors and matrices. It verifies the same certificate and writes a CSV trajectory for independent plotting or regression testing.

Chapter20_Lesson3.java

// Chapter20_Lesson3.java
// Self-contained gain-scheduled stability verification and RK4 simulation.
//
// Build:
//   javac Chapter20_Lesson3.java
// Run:
//   java Chapter20_Lesson3

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Locale;

public final class Chapter20_Lesson3 {
    private static final Matrix2 P = new Matrix2(
        1.23932496, 0.10204082,
        0.10204082, 0.18524333
    );

    private Chapter20_Lesson3() {
    }

    private record Vector2(double x1, double x2) {
        Vector2 add(Vector2 other) {
            return new Vector2(x1 + other.x1, x2 + other.x2);
        }

        Vector2 scale(double scalar) {
            return new Vector2(scalar * x1, scalar * x2);
        }

        double dot(Vector2 other) {
            return x1 * other.x1 + x2 * other.x2;
        }

        double norm() {
            return Math.sqrt(dot(this));
        }
    }

    private record Matrix2(double a11, double a12, double a21, double a22) {
        Matrix2 transpose() {
            return new Matrix2(a11, a21, a12, a22);
        }

        Matrix2 multiply(Matrix2 other) {
            return new Matrix2(
                a11 * other.a11 + a12 * other.a21,
                a11 * other.a12 + a12 * other.a22,
                a21 * other.a11 + a22 * other.a21,
                a21 * other.a12 + a22 * other.a22
            );
        }

        Vector2 multiply(Vector2 x) {
            return new Vector2(
                a11 * x.x1 + a12 * x.x2,
                a21 * x.x1 + a22 * x.x2
            );
        }

        Matrix2 add(Matrix2 other) {
            return new Matrix2(
                a11 + other.a11,
                a12 + other.a12,
                a21 + other.a21,
                a22 + other.a22
            );
        }

        Matrix2 negate() {
            return new Matrix2(-a11, -a12, -a21, -a22);
        }

        double minimumEigenvalueSymmetric() {
            double b = 0.5 * (a12 + a21);
            double trace = a11 + a22;
            double discriminant = Math.sqrt(
                (a11 - a22) * (a11 - a22) + 4.0 * b * b
            );
            return 0.5 * (trace - discriminant);
        }
    }

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

    private static double rhoOfTime(double t, double omega) {
        return clamp(0.5 + 0.5 * Math.sin(omega * t), 0.0, 1.0);
    }

    private static Matrix2 closedLoopMatrix(double rho) {
        return new Matrix2(
            0.0,
            1.0,
            -(4.0 + 1.8 * rho),
            -(2.9 + 0.7 * rho)
        );
    }

    private static double quadraticForm(Vector2 x, Matrix2 matrix) {
        return x.dot(matrix.multiply(x));
    }

    private static Vector2 derivative(double t, Vector2 x, double omega) {
        return closedLoopMatrix(rhoOfTime(t, omega)).multiply(x);
    }

    private static Vector2 rk4Step(double t, Vector2 x, double dt, double omega) {
        Vector2 k1 = derivative(t, x, omega);
        Vector2 k2 = derivative(t + 0.5 * dt, x.add(k1.scale(0.5 * dt)), omega);
        Vector2 k3 = derivative(t + 0.5 * dt, x.add(k2.scale(0.5 * dt)), omega);
        Vector2 k4 = derivative(t + dt, x.add(k3.scale(dt)), omega);

        return x.add(
            k1.add(k2.scale(2.0)).add(k3.scale(2.0)).add(k4).scale(dt / 6.0)
        );
    }

    private static void checkFrozenHurwitz(Matrix2 matrix, double rho) {
        double trace = matrix.a11 + matrix.a22;
        double determinant = matrix.a11 * matrix.a22 - matrix.a12 * matrix.a21;

        if (!(trace < 0.0 && determinant > 0.0)) {
            throw new IllegalStateException(
                "Frozen matrix is not Hurwitz at rho=" + rho
            );
        }
    }

    private static void run() throws IOException {
        double pMargin = P.minimumEigenvalueSymmetric();
        double qMargin = Double.POSITIVE_INFINITY;

        for (int i = 0; i <= 1000; i++) {
            double rho = i / 1000.0;
            Matrix2 a = closedLoopMatrix(rho);
            checkFrozenHurwitz(a, rho);

            Matrix2 q = a.transpose().multiply(P).add(P.multiply(a)).negate();
            qMargin = Math.min(qMargin, q.minimumEigenvalueSymmetric());
        }

        if (!(pMargin > 0.0 && qMargin > 0.0)) {
            throw new IllegalStateException("Common Lyapunov certificate failed.");
        }

        System.out.printf(Locale.US, "lambda_min(P) = %.10f%n", pMargin);
        System.out.printf(
            Locale.US,
            "minimum sampled lambda_min(Q(rho)) = %.10f%n",
            qMargin
        );

        double dt = 0.005;
        double finalTime = 20.0;
        double omega = 0.8;
        int steps = (int) Math.round(finalTime / dt);
        Vector2 x = new Vector2(1.0, -0.5);
        double maxVdot = -Double.MAX_VALUE;

        try (
            PrintWriter output = new PrintWriter(
                new BufferedWriter(
                    new FileWriter("Chapter20_Lesson3_java_results.csv")
                )
            )
        ) {
            output.println("time,rho,x1,x2,V,Vdot");

            for (int k = 0; k <= steps; k++) {
                double t = k * dt;
                double rho = rhoOfTime(t, omega);
                Matrix2 a = closedLoopMatrix(rho);
                Matrix2 s = a.transpose().multiply(P).add(P.multiply(a));
                double v = quadraticForm(x, P);
                double vdot = quadraticForm(x, s);
                maxVdot = Math.max(maxVdot, vdot);

                output.printf(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t, rho, x.x1, x.x2, v, vdot
                );

                if (k < steps) {
                    x = rk4Step(t, x, dt, omega);
                }
            }
        }

        System.out.printf(Locale.US, "final state norm = %.10e%n", x.norm());
        System.out.printf(
            Locale.US,
            "maximum analytical Vdot along trajectory = %.10e%n",
            maxVdot
        );

        if (maxVdot > 1.0e-10) {
            throw new IllegalStateException("Vdot became positive.");
        }
    }

    public static void main(String[] args) {
        try {
            run();
        } catch (Exception error) {
            System.err.println("Error: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script checks the common Lyapunov margin and integrates the scheduled closed loop with ode45. The local function gs_dynamics is also suitable for a MATLAB Function block in Simulink: connect simulation time and the integrated state to the function, and return the state derivative to an Integrator block.

Chapter20_Lesson3.m

% Chapter20_Lesson3.m
% Gain-scheduled stability verification and simulation.
%
% Requirements:
%   MATLAB base functions. Control System Toolbox is optional.
%
% Simulink use:
%   The local function gs_dynamics(t,x) can be copied into a MATLAB
%   Function block. Feed Clock time t and a two-element state vector x,
%   then integrate the returned derivative with an Integrator block.

clear; clc; close all;

P = [1.23932496, 0.10204082; ...
     0.10204082, 0.18524333];

rhoGrid = linspace(0, 1, 1001);
qMargin = inf;
frozenSpectralAbscissa = -inf;

for k = 1:numel(rhoGrid)
    rho = rhoGrid(k);
    Acl = closedLoopMatrix(rho);
    Q = -(Acl' * P + P * Acl);

    qMargin = min(qMargin, min(eig((Q + Q') / 2)));
    frozenSpectralAbscissa = max( ...
        frozenSpectralAbscissa, max(real(eig(Acl))) ...
    );
end

pMargin = min(eig((P + P') / 2));

fprintf('lambda_min(P) = %.10f\n', pMargin);
fprintf('minimum sampled lambda_min(Q(rho)) = %.10f\n', qMargin);
fprintf('worst frozen spectral abscissa = %.10f\n', frozenSpectralAbscissa);

assert(pMargin > 0, 'P is not positive definite.');
assert(qMargin > 0, 'The common Lyapunov certificate failed.');
assert(frozenSpectralAbscissa < 0, 'A frozen closed loop is unstable.');

x0 = [1.0; -0.5];
timeSpan = [0, 20];
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, x] = ode45(@gs_dynamics, timeSpan, x0, options);

rhoHistory = arrayfun(@rhoOfTime, t);
V = zeros(size(t));
Vdot = zeros(size(t));

for k = 1:numel(t)
    xk = x(k, :)';
    Acl = closedLoopMatrix(rhoHistory(k));
    V(k) = xk' * P * xk;
    Vdot(k) = xk' * (Acl' * P + P * Acl) * xk;
end

fprintf('final state norm = %.10e\n', norm(x(end, :)));
fprintf('maximum analytical Vdot along trajectory = %.10e\n', max(Vdot));
assert(max(Vdot) <= 1e-10, 'Vdot became positive.');

results = table(t, rhoHistory, x(:,1), x(:,2), V, Vdot, ...
    'VariableNames', {'time','rho','x1','x2','V','Vdot'});
writetable(results, 'Chapter20_Lesson3_matlab_results.csv');

figure('Name', 'Chapter20 Lesson 3');
tiledlayout(3,1);

nexttile;
plot(t, x(:,1), t, x(:,2), 'LineWidth', 1.2);
grid on;
ylabel('state');
legend('x_1','x_2');

nexttile;
plot(t, rhoHistory, 'LineWidth', 1.2);
grid on;
ylabel('\rho(t)');

nexttile;
plot(t, V, t, Vdot, 'LineWidth', 1.2);
grid on;
xlabel('time [s]');
ylabel('Lyapunov quantities');
legend('V','dV/dt');

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


function dx = gs_dynamics(t, x)
    rho = rhoOfTime(t);
    dx = closedLoopMatrix(rho) * x;
end


function rho = rhoOfTime(t)
    omega = 0.8;
    rho = min(1.0, max(0.0, 0.5 + 0.5 * sin(omega * t)));
end


function Acl = closedLoopMatrix(rho)
    A = [0.0, 1.0; ...
        -(1.0 + 0.8 * rho), -(0.4 + 0.2 * rho)];
    B = [0.0; 1.0];
    K = [3.0 + rho, 2.5 + 0.5 * rho];
    Acl = A - B * K;
end

The notebook expression performs the same frozen-pole, Lyapunov-margin, trajectory, and export calculations using built-in Wolfram Language functions.

Chapter20_Lesson3.nb

Notebook[{
  Cell["Chapter 20, Lesson 3: Stability Considerations in Gain-Scheduled Controllers", "Title"],
  Cell["This notebook verifies frozen-point stability, checks a common quadratic Lyapunov certificate over rho in [0,1], and simulates the scheduled closed loop.", "Text"],
  Cell[
    BoxData[StringJoin[
      "ClearAll[rhoOfTime, plantMatrix, inputMatrix, scheduledGain,\n",
      "  closedLoopMatrix, p, rhoGrid, qMargins, frozenAbscissa,\n",
      "  solution, timeGrid, stateHistory, rhoHistory, vHistory, vdotHistory];\n\n",
      "rhoOfTime[t_] := Clip[0.5 + 0.5 Sin[0.8 t], {0.0, 1.0}];\n\n",
      "plantMatrix[rho_] := { {0.0, 1.0},\n  {-(1.0 + 0.8 rho), -(0.4 + 0.2 rho)} };\n\n",
      "inputMatrix = { {0.0}, {1.0} };\n",
      "scheduledGain[rho_] := { {3.0 + rho, 2.5 + 0.5 rho} };\n",
      "closedLoopMatrix[rho_] :=\n  plantMatrix[rho] - inputMatrix . scheduledGain[rho];\n\n",
      "p = { {1.23932496, 0.10204082},\n  {0.10204082, 0.18524333} };\n\n",
      "rhoGrid = Subdivide[0.0, 1.0, 1000];\n\n",
      "qMargins = Table[\n  Min[Eigenvalues[-(Transpose[closedLoopMatrix[rho]] . p +\n",
      "    p . closedLoopMatrix[rho])]],\n  {rho, rhoGrid}\n];\n\n",
      "frozenAbscissa = Max[\n  Flatten[Re[Eigenvalues[closedLoopMatrix[#]]] & /@ rhoGrid]\n];\n\n",
      "Print[\"lambda_min(P) = \", Min[Eigenvalues[p]]];\n",
      "Print[\"minimum sampled lambda_min(Q(rho)) = \", Min[qMargins]];\n",
      "Print[\"worst frozen spectral abscissa = \", frozenAbscissa];\n\n",
      "If[Min[Eigenvalues[p]] <= 0 || Min[qMargins] <= 0 || frozenAbscissa >= 0,\n",
      "  Print[\"Certificate failed.\"];\n  Abort[]\n];\n\n",
      "solution = NDSolveValue[\n  {\n",
      "    x1'[t] == (closedLoopMatrix[rhoOfTime[t]] . {x1[t], x2[t]})[[1]],\n",
      "    x2'[t] == (closedLoopMatrix[rhoOfTime[t]] . {x1[t], x2[t]})[[2]],\n",
      "    x1[0] == 1.0,\n    x2[0] == -0.5\n  },\n",
      "  {x1, x2},\n  {t, 0, 20},\n",
      "  Method -> {\"TimeIntegration\" -> {\"ExplicitRungeKutta\"} }\n];\n\n",
      "timeGrid = Subdivide[0.0, 20.0, 4000];\n",
      "stateHistory = ({solution[[1]][#], solution[[2]][#]} &) /@ timeGrid;\n",
      "rhoHistory = rhoOfTime /@ timeGrid;\n",
      "vHistory = Map[# . p . # &, stateHistory];\n\n",
      "vdotHistory = MapThread[\n  Function[{x, rho},\n",
      "    x . (Transpose[closedLoopMatrix[rho]] . p +\n      p . closedLoopMatrix[rho]) . x\n  ],\n",
      "  {stateHistory, rhoHistory}\n];\n\n",
      "Print[\"final state norm = \", Norm[Last[stateHistory]]];\n",
      "Print[\"maximum analytical Vdot = \", Max[vdotHistory]];\n\n",
      "Export[\"Chapter20_Lesson3_mathematica_results.csv\",\n  Prepend[\n",
      "    MapThread[Join[{ #1, #2 }, #3, { #4, #5 }] &,\n",
      "      {timeGrid, rhoHistory, stateHistory, vHistory, vdotHistory}],\n",
      "    {\"time\", \"rho\", \"x1\", \"x2\", \"V\", \"Vdot\"}\n  ]\n];\n\n",
      "Export[\"Chapter20_Lesson3_mathematica.png\",\n  GraphicsGrid[{\n",
      "    {ListLinePlot[Transpose[{timeGrid, stateHistory[[All, 1]]}],\n",
      "       PlotLegends -> {\"x1\"}, Frame -> True, FrameLabel -> {\"time\", \"state\"}],\n",
      "     ListLinePlot[Transpose[{timeGrid, stateHistory[[All, 2]]}],\n",
      "       PlotLegends -> {\"x2\"}, Frame -> True, FrameLabel -> {\"time\", \"state\"}]},\n",
      "    {ListLinePlot[Transpose[{timeGrid, rhoHistory}],\n",
      "       PlotLegends -> {\"rho\"}, Frame -> True, FrameLabel -> {\"time\", \"schedule\"}],\n",
      "     ListLinePlot[{Transpose[{timeGrid, vHistory}],\n",
      "       Transpose[{timeGrid, vdotHistory}]}, PlotLegends -> {\"V\", \"Vdot\"},\n",
      "       Frame -> True, FrameLabel -> {\"time\", \"Lyapunov quantities\"}]}\n",
      "  }],\n  ImageResolution -> 180\n];"
    ]],
    "Input"
  ]}]        

A defensible gain-scheduled design should be checked in the following order.

1. Define the admissible scheduling set, measurement errors, filtering, delay, and rate limits.
2. Form the complete augmented closed-loop model, including controller states, filters, integrators, actuator dynamics, and interpolation.
3. Verify local stabilizability and frozen closed-loop poles at design points and between them.
4. Search first for a common quadratic Lyapunov certificate. For a polytopic affine model, test all vertices.
5. If the CQLF test fails, use a parameter-dependent Lyapunov function and include the scheduling-rate term explicitly.
6. Test scheduling-variable mismatch, noise, delay, rate saturation, actuator saturation, and controller-state initialization.
7. Simulate aggressive trajectories that traverse the scheduling envelope, not only constant operating points.
8. Implement online envelope checks and a safe fallback controller for operation outside the certified region.

Dense simulation is evidence, not proof. The analytical certificate establishes stability; simulation is then used to expose modeling and implementation errors not represented in the proof.

Problem 1 (Frozen Stability versus Scheduled Stability): Explain why the condition \( A_{\mathrm{cl} }(\rho) \) Hurwitz for every fixed \( \rho \) does not by itself prove stability for \( \rho=\rho(t) \).

Solution: Frozen stability studies a family of time-invariant systems separately. The scheduled system is time-varying, and its state-transition matrix depends on the ordered evolution of the matrices. Different matrices generally do not commute, their decay directions can rotate, and dynamic-controller realizations can introduce terms proportional to \( \dot{\rho} \). A common or parameter-dependent Lyapunov inequality is needed to control all admissible trajectories.

Problem 2 (CQLF Decay Bound): Suppose \( P>0 \) and \( A_{\mathrm{cl} }^\top P+P A_{\mathrm{cl} } \leq -Q \) for all schedules, with \( Q>0 \). Derive an exponential norm bound.

Solution: Since \( \dot{V}\leq-\lambda_{\min}(Q)\|x\|^2 \) and \( V\leq\lambda_{\max}(P)\|x\|^2 \),

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

Therefore

\[ V(t) \leq e^{-ct}V(0), \qquad c=\frac{\lambda_{\min}(Q)} {\lambda_{\max}(P)}. \]

Applying lower and upper eigenvalue bounds to \( V \) gives

\[ \|x(t)\| \leq \sqrt{\frac{\lambda_{\max}(P)} {\lambda_{\min}(P)} } e^{-ct/2}\|x(0)\|. \]

Problem 3 (Vertex Sufficiency): Let \( A(\rho)=\sum_i\alpha_i(\rho)A_i \) with convex weights. Show that one common \( P>0 \) satisfying \( A_i^\top P+P A_i<0 \) at every vertex also certifies every interior point.

Solution: Linearity gives

\[ A^\top(\rho)P+P A(\rho) = \sum_i\alpha_i(\rho) \left(A_i^\top P+P A_i\right). \]

A convex combination of negative-definite matrices is negative definite. Hence the same \( P \) proves stability over the entire polytope.

Problem 4 (Rate Bound for a Parameter-Dependent Certificate): Assume \( A^\top P+PA\leq-qI \), \( \|\partial P/\partial\rho\|\leq\mu \), and \( |\dot{\rho}|\leq\nu \). Find a sufficient stability condition.

Solution:

\[ \dot{V} \leq -q\|x\|^2+ \mu\nu\|x\|^2 = -(q-\mu\nu)\|x\|^2. \]

Uniform decay follows if \( \mu\nu<q \), or equivalently \( \nu<q/\mu \) when \( \mu>0 \).

Problem 5 (Worked-Example Certificate): For the matrix

\[ A_{\mathrm{cl} }(\rho)= \begin{bmatrix} 0 & 1\\ -(4+1.8\rho) & -(2.9+0.7\rho) \end{bmatrix}, \]

show that every frozen system is Hurwitz for \( 0\leq\rho\leq 1 \).

Solution: For a real two-by-two matrix, Hurwitz stability is equivalent to negative trace and positive determinant. Here,

\[ \operatorname{tr}A_{\mathrm{cl} } =-(2.9+0.7\rho)<0, \qquad \det A_{\mathrm{cl} } =4+1.8\rho>0. \]

Both inequalities hold throughout the interval. This proves frozen stability. The common matrix in Section 8 provides the stronger trajectory-stability result.

Stable frozen controllers are only the first step in a gain-scheduled design. A common quadratic Lyapunov function proves uniform exponential stability for arbitrary admissible scheduling rates. A parameter-dependent Lyapunov function can reduce conservatism, but its derivative introduces explicit scheduling-rate terms. Polytopic vertex tests convert some infinite-dimensional checks into finite LMIs. Dynamic-controller interpolation must preserve realization consistency and account for hidden coupling. The worked example demonstrates a complete certificate whose validity is independent of scheduling speed.