After completing this lesson, the student should be able to:

• construct normalized interpolation coordinates from measured scheduling variables;
• interpolate state-feedback gains, feedforward gains, and compatible controller realizations;
• distinguish a practical LPV-like schedule from a formally synthesized LPV controller;
• derive an interpolation-error bound for a smooth ideal gain map;
• state and prove a common-quadratic-Lyapunov stability certificate;
• recognize why stable endpoint controllers do not automatically imply a stable interpolation;
• implement clamping, breakpoint selection, simulation, and frozen-time validation.

The phrase LPV-like is intentional. We use an online measured parameter \( \rho(t) \) to blend controllers, but we do not claim that arbitrary interpolation supplies the global guarantees of a complete LPV synthesis. The method is best understood as a disciplined gain-scheduling construction whose dependence on \( \rho \) is made explicit.

Suppose a nonlinear or parameter-dependent plant is represented locally by a family of linear models

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

At grid points \( \rho_0,\rho_1,\ldots,\rho_N \), linear-control methods produce local controllers. For state feedback these may be

\[ \mathbf{u}=-\mathbf{K}_i\mathbf{x}+\mathbf{K}_{r,i}\mathbf{r}, \qquad i=0,1,\ldots,N. \]

Online interpolation replaces the finite table by a continuous map \( \rho \mapsto \mathbf{K}(\rho) \). The implementation therefore contains a measurement path, breakpoint logic, interpolation weights, and a controller evaluation path.

flowchart TD
  A["Plant measurements"] --> B["Compute scheduling variable rho"]
  B --> C["Clamp rho to validated envelope"]
  C --> D["Locate neighboring grid points"]
  D --> E["Compute interpolation weights"]
  E --> F["Blend gains and feedforward terms"]
  F --> G["Evaluate scheduled control law"]
  G --> H["Actuator command"]
  H --> A
        

The scheduling variable should describe the changing plant dynamics and should be measurable or estimable with acceptable delay and noise. Examples include speed, altitude, dynamic pressure, payload mass, temperature, or joint configuration.

Let \( \rho_i \le \rho \le \rho_{i+1} \). Define the normalized local coordinate

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

The two nonnegative convex weights are

\[ \lambda_i(\rho)=1-\alpha(\rho), \qquad \lambda_{i+1}(\rho)=\alpha(\rho), \qquad \lambda_i+\lambda_{i+1}=1. \]

A matrix-valued controller gain is interpolated elementwise:

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

The same weights should normally be used for all mutually designed quantities: \( \mathbf{K} \), reference feedforward, observer gains, command filters, and anti-windup parameters. Mixing unrelated schedules can destroy the relationships imposed by the local design.

flowchart TD
  C0["Controller at \nleft grid point"] --> W0["Multiply by 1-alpha"]
  C1["Controller at \nright grid point"] --> W1["Multiply by alpha"]
  W0 --> S["Add weighted controller data"]
  W1 --> S
  R["Measured rho"] --> A["Compute alpha"]
  A --> W0
  A --> W1
  S --> U["Scheduled control command"]
        

Outside the validated range, the safest basic policy is usually endpoint clamping:

\[ \rho_c=\min\!\left(\rho_N,\max(\rho_0,\rho)\right). \]

Unrestricted extrapolation can create gains that were never analyzed and can amplify sensor faults or envelope excursions.

4.1 Convex-hull property

Because the weights are nonnegative and sum to one, every interpolated gain lies in the line segment joining its two endpoint gains. For any convex matrix norm,

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

This bounds the interpolated gain norm by a convex combination of endpoint norms, but it does not prove closed-loop stability because eigenvalues and feedback interconnections are nonlinear functions of controller coefficients.

4.2 Continuity and differentiability

Piecewise-linear interpolation is continuous when adjacent intervals share the same stored breakpoint values. Its derivative is generally discontinuous at breakpoints. Consequently, a smoothly varying \( \rho(t) \) gives a continuous controller gain but may produce slope changes in the command. Rate limiting or a higher-order interpolant may be appropriate when actuator jerk is important.

4.3 Error theorem for a smooth ideal gain map

Let a scalar ideal gain \( k^{\star}(\rho) \) be twice continuously differentiable on \( [a,b] \), and let \( k_L(\rho) \) be the straight line through the endpoint values. The classical interpolation remainder gives, for some \( \xi \in (a,b) \),

\[ k^{\star}(\rho)-k_L(\rho) =\frac{k^{\star\prime\prime}(\xi)}{2}(\rho-a)(\rho-b). \]

Since the maximum of \( |(\rho-a)(\rho-b)| \) on the interval is \( (b-a)^2/4 \),

\[ |k^{\star}(\rho)-k_L(\rho)| \le \frac{(b-a)^2}{8} \max_{\eta\in[a,b]}|k^{\star\prime\prime}(\eta)|. \]

Therefore, halving the grid spacing reduces the worst-case interpolation error by a factor of approximately four when the curvature bound remains unchanged. The same estimate can be applied entrywise to a matrix gain; a compatible matrix-norm version follows by bounding \( \|\mathbf{K}^{\star\prime\prime}(\rho)\| \).

Proof sketch: define the interpolation error, subtract a scaled quadratic that vanishes at \( a \), \( \rho \), and \( b \), and apply Rolle's theorem twice. The second derivative must vanish at an intermediate point, which yields the remainder formula above.

For two normalized scheduling coordinates \( \theta_1,\theta_2\in[0,1] \), bilinear interpolation over a rectangular cell uses four vertex weights:

\[ \begin{aligned} \lambda_{00}&=(1-\theta_1)(1-\theta_2),\\ \lambda_{10}&=\theta_1(1-\theta_2),\\ \lambda_{01}&=(1-\theta_1)\theta_2,\\ \lambda_{11}&=\theta_1\theta_2. \end{aligned} \]

These satisfy

\[ \lambda_{jk}\ge 0, \qquad \sum_{j=0}^{1}\sum_{k=0}^{1}\lambda_{jk}=1, \qquad \mathbf{K}(\boldsymbol{\theta})= \sum_{j=0}^{1}\sum_{k=0}^{1}\lambda_{jk}\mathbf{K}_{jk}. \]

With \( p \) scheduling variables, a hyper-rectangular cell has \( 2^p \) vertices. This exponential growth motivates sparse grids, triangulation, tensor-product approximations, or direct parameter-dependent synthesis. On a simplex, barycentric coordinates replace tensor-product weights and use only \( p+1 \) vertices per cell.

6.1 Static gains

State-feedback, output-feedback, feedforward, and PID gains are the simplest objects to interpolate. Even here, endpoint stability alone is insufficient: the closed-loop matrix between endpoints must be analyzed.

6.2 Dynamic controller realizations

Consider local dynamic controllers

\[ \dot{\mathbf{x} }_K=\mathbf{A}_{K,i}\mathbf{x}_K +\mathbf{B}_{K,i}\mathbf{e}, \qquad \mathbf{u}=\mathbf{C}_{K,i}\mathbf{x}_K+\mathbf{D}_{K,i}\mathbf{e}. \]

Direct interpolation of the four matrices is meaningful only when every local controller uses the same state dimension and compatible state coordinates. Two controllers can have the same transfer function but realizations related by a similarity transformation. Interpolating unaligned realizations can create an entirely different and possibly unstable controller.

6.3 Parallel blending and bumpless behavior

An alternative is to keep neighboring controllers active in parallel and blend their outputs:

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

If both controller states are continuously updated and the weights are continuous, the output blend is naturally continuous. This costs extra computation but avoids abrupt state initialization. Integrator states, saturations, and anti-windup feedback must still be coordinated.

6.4 Scheduling on endogenous variables

If \( \rho=h(\mathbf{x}) \), the schedule is part of a nonlinear feedback loop. Frozen-time analysis treats \( \rho \) as fixed and therefore omits this coupling. The resulting controller should be validated on the complete nonlinear model, not only on the local linear models.

Assume the interpolated closed-loop matrix is an affine convex combination

\[ \mathbf{A}_{cl}(\rho)=\sum_{j=1}^{M}\lambda_j(\rho)\mathbf{A}_{cl,j}, \qquad \lambda_j\ge 0, \qquad \sum_{j=1}^{M}\lambda_j=1. \]

Suppose there exists one symmetric matrix \( \mathbf{P}=\mathbf{P}^{T}>0 \) and one symmetric \( \mathbf{Q}=\mathbf{Q}^{T}>0 \) such that every vertex satisfies

\[ \mathbf{A}_{cl,j}^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_{cl,j}\le -\mathbf{Q}, \qquad j=1,\ldots,M. \]

Choose \( V(\mathbf{x})=\mathbf{x}^{T}\mathbf{P}\mathbf{x} \). Then

\[ \begin{aligned} \dot V &=\mathbf{x}^{T}\!\left(\mathbf{A}_{cl}(\rho)^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_{cl}(\rho)\right)\!\mathbf{x}\\ &=\sum_{j=1}^{M}\lambda_j \mathbf{x}^{T}\!\left(\mathbf{A}_{cl,j}^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_{cl,j}\right)\!\mathbf{x}\\ &\le -\mathbf{x}^{T}\mathbf{Q}\mathbf{x}. \end{aligned} \]

Thus the entire convex family is uniformly exponentially stable, even when the scheduling variable changes with time, because the Lyapunov matrix is constant and no \( \dot\rho \) term appears.

This certificate applies only when the actual closed-loop matrix is affine in the weights. If both \( \mathbf{B}(\rho) \) and \( \mathbf{K}(\rho) \) are interpolated, then

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

generally contains products of weights and is not simply the convex combination of endpoint closed-loop matrices. A parameter-dependent Lyapunov function \( \mathbf{P}(\rho) \) can reduce conservatism, but then

\[ \dot V=\mathbf{x}^{T}\!\left( \mathbf{A}_{cl}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{cl} +\frac{\partial\mathbf{P} }{\partial\rho}\dot\rho \right)\!\mathbf{x}, \]

so the admissible rate of parameter variation becomes part of the certificate.

Consider the LPV-like second-order model

\[ m(\rho)\ddot q+c\dot q+kq=u, \qquad m(\rho)=1+2\rho, \qquad \rho\in[0,1], \]

with \( c=0.5 \) and \( k=2 \). Define \( \mathbf{x}=[q\;\;\dot q]^{T} \) and use

\[ u=-k_p(\rho)q-k_v(\rho)\dot q+k_r(\rho)r. \]

The frozen closed-loop characteristic polynomial is

\[ s^2+\frac{c+k_v(\rho)}{m(\rho)}s +\frac{k+k_p(\rho)}{m(\rho)}. \]

Select desired poles at \( -2 \) and \( -3 \), so the desired polynomial is \( s^2+5s+6 \). Coefficient matching gives

\[ k_p(\rho)=6m(\rho)-k, \qquad k_v(\rho)=5m(\rho)-c. \]

At the scheduling endpoints:

\[ \begin{aligned} \rho=0,\;m=1:&\quad \mathbf{K}_0=[4\;\;4.5],\quad k_{r,0}=6,\\ \rho=1,\;m=3:&\quad \mathbf{K}_1=[16\;\;14.5],\quad k_{r,1}=18. \end{aligned} \]

Because \( m(\rho) \) is affine and the required gains are affine in mass, ordinary interpolation with \( \alpha=\rho \) reproduces the exact analytical schedule:

\[ \begin{aligned} \mathbf{K}(\rho)&=(1-\rho)\mathbf{K}_0+\rho\mathbf{K}_1 =[4+12\rho\;\;4.5+10\rho],\\ k_r(\rho)&=(1-\rho)6+18\rho=6+12\rho=6m(\rho). \end{aligned} \]

Substitution produces

\[ \ddot q+5\dot q+6q=6r, \]

so every frozen model has the same poles and unit steady-state gain. In this particular simplified model the cancellation also removes explicit dependence on a time-varying \( \rho(t) \). This should not be generalized: a physical variable-mass momentum model can contain additional terms proportional to the mass rate, and general LPV plants do not admit such exact cancellation.

1. Acquire or estimate \( \rho \).
2. Apply plausibility checks, filtering if justified, and envelope clamping.
3. Locate the active interval using a monotone breakpoint table.
4. Compute normalized weights with protection against duplicate breakpoints.
5. Blend all coordinated controller quantities with the same weights.
6. Apply actuator limits and anti-windup logic.
7. Log \( \rho \), interval index, weights, gains, command, and saturation status.

At minimum, validation should include:

• endpoint equality and continuity at every breakpoint;
• frozen-time eigenvalue or frequency-response sweeps over a dense parameter grid;
• time-varying scheduling trajectories at nominal and maximum expected rates;
• noise, delay, quantization, and out-of-range scheduling tests;
• actuator saturation and bumplessness tests;
• nonlinear-plant simulation, not only local-model simulation.

This implementation uses NumPy and SciPy. It interpolates the endpoint gains, checks the frozen-time eigenvalues, integrates the time-varying closed loop, and writes a CSV file for plotting or regression testing.

Chapter20_Lesson4.py

"""Chapter 20, Lesson 4: Simple LPV-Like Interpolation of Linear Controllers.

Simulates a mass-spring-damper plant with a measured scheduling variable rho.
Two endpoint state-feedback controllers are linearly interpolated online.
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import numpy as np
from scipy.integrate import solve_ivp


@dataclass(frozen=True)
class ScheduledController:
    rho_min: float = 0.0
    rho_max: float = 1.0
    k_spring: float = 2.0
    c_damping: float = 0.5

    # Endpoint gains for desired poles {-2, -3}.
    k0: tuple[float, float] = (4.0, 4.5)
    k1: tuple[float, float] = (16.0, 14.5)
    kr0: float = 6.0
    kr1: float = 18.0

    def alpha(self, rho: float) -> float:
        raw = (rho - self.rho_min) / (self.rho_max - self.rho_min)
        return float(np.clip(raw, 0.0, 1.0))

    def gains(self, rho: float) -> tuple[np.ndarray, float]:
        a = self.alpha(rho)
        k0 = np.asarray(self.k0, dtype=float)
        k1 = np.asarray(self.k1, dtype=float)
        k = (1.0 - a) * k0 + a * k1
        kr = (1.0 - a) * self.kr0 + a * self.kr1
        return k, kr

    @staticmethod
    def mass(rho: float) -> float:
        return 1.0 + 2.0 * rho


def scheduling_variable(t: float) -> float:
    """Measured operating condition in [0, 1]."""
    return 0.5 * (1.0 + math.sin(0.35 * t))


def reference(t: float) -> float:
    return 1.0 if t >= 1.0 else 0.0


def closed_loop_rhs(
    t: float, x: np.ndarray, controller: ScheduledController
) -> np.ndarray:
    rho = scheduling_variable(t)
    mass = controller.mass(rho)
    k_gain, kr = controller.gains(rho)
    r = reference(t)

    u = -float(k_gain @ x) + kr * r
    position, velocity = x
    acceleration = (
        u - controller.c_damping * velocity - controller.k_spring * position
    ) / mass
    return np.array([velocity, acceleration], dtype=float)


def frozen_closed_loop_eigenvalues(
    rho: float, controller: ScheduledController
) -> np.ndarray:
    mass = controller.mass(rho)
    k_gain, _ = controller.gains(rho)
    a_cl = np.array(
        [
            [0.0, 1.0],
            [-(controller.k_spring + k_gain[0]) / mass,
             -(controller.c_damping + k_gain[1]) / mass],
        ]
    )
    return np.linalg.eigvals(a_cl)


def main() -> None:
    controller = ScheduledController()
    t_eval = np.linspace(0.0, 25.0, 2501)
    solution = solve_ivp(
        fun=lambda t, x: closed_loop_rhs(t, x, controller),
        t_span=(t_eval[0], t_eval[-1]),
        y0=np.array([0.0, 0.0]),
        t_eval=t_eval,
        rtol=1e-9,
        atol=1e-11,
    )
    if not solution.success:
        raise RuntimeError(solution.message)

    sample_rho = np.linspace(0.0, 1.0, 6)
    print("Frozen-time eigenvalue check:")
    for rho in sample_rho:
        eig = np.sort_complex(frozen_closed_loop_eigenvalues(rho, controller))
        print(f"rho={rho:4.2f}, eigenvalues={eig}")

    output = np.column_stack(
        (
            solution.t,
            solution.y[0],
            solution.y[1],
            [scheduling_variable(t) for t in solution.t],
            [reference(t) for t in solution.t],
        )
    )
    np.savetxt(
        "Chapter20_Lesson4_python.csv",
        output,
        delimiter=",",
        header="time,position,velocity,rho,reference",
        comments="",
    )
    print("Saved Chapter20_Lesson4_python.csv")


if __name__ == "__main__":
    main()

The C++17 version has no external numerical dependency. A fixed-step fourth-order Runge-Kutta integrator is included from scratch, making the scheduling and integration logic explicit.

Chapter20_Lesson4.cpp

// Chapter 20, Lesson 4: Simple LPV-Like Interpolation of Linear Controllers
// Build: g++ -std=c++17 -O2 Chapter20_Lesson4.cpp -o Chapter20_Lesson4

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

struct Controller {
    static constexpr double kSpring = 2.0;
    static constexpr double cDamping = 0.5;
    static constexpr std::array<double, 2> K0{4.0, 4.5};
    static constexpr std::array<double, 2> K1{16.0, 14.5};
    static constexpr double kr0 = 6.0;
    static constexpr double kr1 = 18.0;

    static double clamp01(double value) {
        return std::clamp(value, 0.0, 1.0);
    }

    static double mass(double rho) {
        return 1.0 + 2.0 * rho;
    }

    static std::pair<std::array<double, 2>, double> gains(double rho) {
        const double alpha = clamp01(rho);
        std::array<double, 2> K{
            (1.0 - alpha) * K0[0] + alpha * K1[0],
            (1.0 - alpha) * K0[1] + alpha * K1[1]
        };
        const double kr = (1.0 - alpha) * kr0 + alpha * kr1;
        return {K, kr};
    }
};

using State = std::array<double, 2>;

static double rhoOfTime(double t) {
    return 0.5 * (1.0 + std::sin(0.35 * t));
}

static double reference(double t) {
    return t >= 1.0 ? 1.0 : 0.0;
}

static State rhs(double t, const State& x) {
    const double rho = rhoOfTime(t);
    const double m = Controller::mass(rho);
    const auto [K, kr] = Controller::gains(rho);
    const double u = -(K[0] * x[0] + K[1] * x[1]) + kr * reference(t);
    const double acceleration =
        (u - Controller::cDamping * x[1] - Controller::kSpring * x[0]) / m;
    return {x[1], acceleration};
}

static State addScaled(const State& x, const State& dx, double scale) {
    return {x[0] + scale * dx[0], x[1] + scale * dx[1]};
}

static State rk4Step(double t, const State& x, double h) {
    const State k1 = rhs(t, x);
    const State k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h));
    const State k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h));
    const State k4 = rhs(t + h, addScaled(x, k3, h));
    return {
        x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
        x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
    };
}

static std::array<double, 2> frozenEigenvalues(double rho) {
    // For this worked design the characteristic polynomial is s^2 + 5 s + 6.
    const auto [K, unusedKr] = Controller::gains(rho);
    (void)unusedKr;
    const double m = Controller::mass(rho);
    const double a1 = (Controller::cDamping + K[1]) / m;
    const double a0 = (Controller::kSpring + K[0]) / m;
    const double discriminant = a1 * a1 - 4.0 * a0;
    if (discriminant < 0.0) {
        throw std::runtime_error("This demonstration expects real frozen poles.");
    }
    return {(-a1 - std::sqrt(discriminant)) / 2.0,
            (-a1 + std::sqrt(discriminant)) / 2.0};
}

int main() {
    std::cout << "Frozen-time eigenvalue check:\n";
    for (int i = 0; i <= 5; ++i) {
        const double rho = i / 5.0;
        const auto eig = frozenEigenvalues(rho);
        std::cout << "rho=" << std::fixed << std::setprecision(2) << rho
                  << ", eigenvalues=[" << eig[0] << ", " << eig[1] << "]\n";
    }

    std::ofstream csv("Chapter20_Lesson4_cpp.csv");
    if (!csv) {
        std::cerr << "Cannot create output CSV.\n";
        return 1;
    }
    csv << "time,position,velocity,rho,reference\n";

    const double h = 0.01;
    const int steps = static_cast<int>(25.0 / h);
    State x{0.0, 0.0};
    double t = 0.0;
    for (int step = 0; step <= steps; ++step) {
        csv << std::setprecision(12) << t << ',' << x[0] << ',' << x[1] << ','
            << rhoOfTime(t) << ',' << reference(t) << '\n';
        x = rk4Step(t, x, h);
        t += h;
    }
    std::cout << "Saved Chapter20_Lesson4_cpp.csv\n";
    return 0;
}

The Java version mirrors the C++ structure, uses a hand-coded RK4 step, and writes a deterministic CSV record for verification.

Chapter20_Lesson4.java

// Chapter 20, Lesson 4: Simple LPV-Like Interpolation of Linear Controllers
// Build: javac Chapter20_Lesson4.java
// Run:   java Chapter20_Lesson4

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

public final class Chapter20_Lesson4 {
    private static final double K_SPRING = 2.0;
    private static final double C_DAMPING = 0.5;
    private static final double[] K0 = {4.0, 4.5};
    private static final double[] K1 = {16.0, 14.5};
    private static final double KR0 = 6.0;
    private static final double KR1 = 18.0;

    private Chapter20_Lesson4() {}

    private static double clamp01(double value) {
        return Math.max(0.0, Math.min(1.0, value));
    }

    private static double mass(double rho) {
        return 1.0 + 2.0 * rho;
    }

    private static double[] gains(double rho) {
        double alpha = clamp01(rho);
        return new double[] {
            (1.0 - alpha) * K0[0] + alpha * K1[0],
            (1.0 - alpha) * K0[1] + alpha * K1[1],
            (1.0 - alpha) * KR0 + alpha * KR1
        };
    }

    private static double rhoOfTime(double t) {
        return 0.5 * (1.0 + Math.sin(0.35 * t));
    }

    private static double reference(double t) {
        return t >= 1.0 ? 1.0 : 0.0;
    }

    private static double[] rhs(double t, double[] x) {
        double rho = rhoOfTime(t);
        double[] g = gains(rho);
        double u = -(g[0] * x[0] + g[1] * x[1]) + g[2] * reference(t);
        double acceleration =
            (u - C_DAMPING * x[1] - K_SPRING * x[0]) / mass(rho);
        return new double[] {x[1], acceleration};
    }

    private static double[] addScaled(double[] x, double[] dx, double scale) {
        return new double[] {x[0] + scale * dx[0], x[1] + scale * dx[1]};
    }

    private static double[] rk4Step(double t, double[] x, double h) {
        double[] k1 = rhs(t, x);
        double[] k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h));
        double[] k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h));
        double[] k4 = rhs(t + h, addScaled(x, k3, h));
        return new double[] {
            x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        };
    }

    private static double[] frozenEigenvalues(double rho) {
        double[] g = gains(rho);
        double a1 = (C_DAMPING + g[1]) / mass(rho);
        double a0 = (K_SPRING + g[0]) / mass(rho);
        double discriminant = a1 * a1 - 4.0 * a0;
        if (discriminant < 0.0) {
            throw new IllegalStateException("Expected real frozen poles.");
        }
        return new double[] {
            (-a1 - Math.sqrt(discriminant)) / 2.0,
            (-a1 + Math.sqrt(discriminant)) / 2.0
        };
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        System.out.println("Frozen-time eigenvalue check:");
        for (int i = 0; i <= 5; i++) {
            double rho = i / 5.0;
            double[] eig = frozenEigenvalues(rho);
            System.out.printf("rho=%.2f, eigenvalues=[%.6f, %.6f]%n",
                              rho, eig[0], eig[1]);
        }

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter20_Lesson4_java.csv"))) {
            writer.write("time,position,velocity,rho,reference\n");
            double h = 0.01;
            int steps = (int) Math.round(25.0 / h);
            double[] x = {0.0, 0.0};
            double t = 0.0;
            for (int step = 0; step <= steps; step++) {
                writer.write(String.format(Locale.US, "%.12f,%.12f,%.12f,%.12f,%.1f%n",
                    t, x[0], x[1], rhoOfTime(t), reference(t)));
                x = rk4Step(t, x, h);
                t += h;
            }
        }
        System.out.println("Saved Chapter20_Lesson4_java.csv");
    }
}

The MATLAB script uses ode45, performs the frozen-pole sweep, exports data, and plots the response. For Simulink, the functions scheduledGains, scheduledMass, and the control equation can be placed in a MATLAB Function block. Connect rho, x1, x2, and reference as inputs, and return u. A 1-D Lookup Table block is an equivalent choice for each stored gain table; enable clipping at the endpoints rather than extrapolation.

Chapter20_Lesson4.m

% Chapter 20, Lesson 4: Simple LPV-Like Interpolation of Linear Controllers
% Requires only base MATLAB. The script simulates with ode45 and verifies
% the frozen-time poles across the scheduling interval.

clear; clc; close all;

p.kSpring = 2.0;
p.cDamping = 0.5;
p.K0 = [4.0, 4.5];
p.K1 = [16.0, 14.5];
p.kr0 = 6.0;
p.kr1 = 18.0;

rhoGrid = linspace(0, 1, 6);
fprintf('Frozen-time eigenvalue check:\n');
for rho = rhoGrid
    [K, ~] = scheduledGains(rho, p);
    m = scheduledMass(rho);
    Acl = [0, 1; -(p.kSpring + K(1))/m, -(p.cDamping + K(2))/m];
    lambda = sort(eig(Acl));
    fprintf('rho=%.2f, eigenvalues=[%.6f, %.6f]\n', rho, lambda(1), lambda(2));
end

opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, x] = ode45(@(t, x) closedLoopRhs(t, x, p), [0, 25], [0; 0], opts);
rho = arrayfun(@schedulingVariable, t);
r = arrayfun(@referenceSignal, t);

results = table(t, x(:,1), x(:,2), rho, r, ...
    'VariableNames', {'time','position','velocity','rho','reference'});
writetable(results, 'Chapter20_Lesson4_matlab.csv');

figure;
plot(t, x(:,1), 'LineWidth', 1.2); hold on;
plot(t, r, '--', 'LineWidth', 1.0);
grid on;
xlabel('Time (s)'); ylabel('Position');
legend('x_1', 'reference', 'Location', 'best');
title('LPV-like interpolated state-feedback response');

figure;
plot(t, rho, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)'); ylabel('\rho(t)');
title('Scheduling variable');

fprintf('Saved Chapter20_Lesson4_matlab.csv\n');

function dx = closedLoopRhs(t, x, p)
    rho = schedulingVariable(t);
    m = scheduledMass(rho);
    [K, kr] = scheduledGains(rho, p);
    u = -K*x + kr*referenceSignal(t);
    acceleration = (u - p.cDamping*x(2) - p.kSpring*x(1))/m;
    dx = [x(2); acceleration];
end

function rho = schedulingVariable(t)
    rho = 0.5*(1 + sin(0.35*t));
end

function r = referenceSignal(t)
    r = double(t >= 1.0);
end

function m = scheduledMass(rho)
    m = 1 + 2*rho;
end

function [K, kr] = scheduledGains(rho, p)
    alpha = min(max(rho, 0), 1);
    K = (1-alpha)*p.K0 + alpha*p.K1;
    kr = (1-alpha)*p.kr0 + alpha*p.kr1;
end

The Mathematica source defines the schedule symbolically, verifies the local eigenvalues, integrates the closed-loop equations with NDSolveValue, exports CSV data, and produces response plots.

Chapter20_Lesson4.nb

(* Chapter 20, Lesson 4: Simple LPV-Like Interpolation of Linear Controllers *)

ClearAll["Global`*"];

kSpring = 2.0;
cDamping = 0.5;
k0 = {4.0, 4.5};
k1 = {16.0, 14.5};
kr0 = 6.0;
kr1 = 18.0;

clip01[z_] := Clip[z, {0.0, 1.0}];
mass[rho_] := 1.0 + 2.0 rho;
rhoOfTime[t_] := 0.5 (1.0 + Sin[0.35 t]);
reference[t_] := Piecewise[{ {1.0, t >= 1.0} }, 0.0];

gains[rho_] := Module[{alpha = clip01[rho]},
  {(1.0 - alpha) k0 + alpha k1,
   (1.0 - alpha) kr0 + alpha kr1}
];

frozenClosedLoop[rho_] := Module[{k, kr, m},
  {k, kr} = gains[rho];
  m = mass[rho];
  { {0.0, 1.0},
   {-(kSpring + k[[1]])/m, -(cDamping + k[[2]])/m} }
];

Print["Frozen-time eigenvalue check:"];
Table[
  Print["rho=", NumberForm[rho, {3, 2}], ", eigenvalues=",
    Eigenvalues[frozenClosedLoop[rho]]],
  {rho, 0.0, 1.0, 0.2}
];

solution = NDSolveValue[
  {
    x1'[t] == x2[t],
    x2'[t] == Module[{rho, k, kr, u, m},
      rho = rhoOfTime[t];
      {k, kr} = gains[rho];
      m = mass[rho];
      u = -(k[[1]] x1[t] + k[[2]] x2[t]) + kr reference[t];
      (u - cDamping x2[t] - kSpring x1[t])/m
    ],
    x1[0] == 0.0,
    x2[0] == 0.0
  },
  {x1, x2},
  {t, 0.0, 25.0},
  Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4} }
];

samples = Table[
  {tt, solution[[1]][tt], solution[[2]][tt], rhoOfTime[tt], reference[tt]},
  {tt, 0.0, 25.0, 0.01}
];
Export["Chapter20_Lesson4_mathematica.csv",
  Prepend[samples, {"time", "position", "velocity", "rho", "reference"}]
];

responsePlot = Plot[
  {solution[[1]][t], reference[t]},
  {t, 0.0, 25.0},
  PlotLegends -> {"x1", "reference"},
  AxesLabel -> {"time", "position"},
  PlotRange -> All,
  GridLines -> Automatic
];

schedulePlot = Plot[
  rhoOfTime[t],
  {t, 0.0, 25.0},
  AxesLabel -> {"time", "rho"},
  PlotRange -> {0.0, 1.0},
  GridLines -> Automatic
];

Print[responsePlot];
Print[schedulePlot];
Print["Saved Chapter20_Lesson4_mathematica.csv"];

Problem 1 (Convex weights): For two breakpoints \( \rho_i<\rho_{i+1} \), prove that the linear interpolation weights are nonnegative and sum to one for every \( \rho\in[\rho_i,\rho_{i+1}] \).

Solution: Since

\[ 0\le \rho-\rho_i\le \rho_{i+1}-\rho_i, \]

division by the positive interval length gives \( 0\le\alpha\le1 \). Therefore \( \lambda_i=1-\alpha\ge0 \) and \( \lambda_{i+1}=\alpha\ge0 \). Their sum is identically one. Thus the interpolated controller belongs to the convex hull of the two stored controllers.

Problem 2 (Numerical scheduled gains): For the mechanical example, compute the mass, feedback gains, and reference feedforward at \( \rho=0.35 \).

Solution:

\[ m=1+2(0.35)=1.7. \]

\[ \begin{aligned} k_p&=4+12(0.35)=8.2,\\ k_v&=4.5+10(0.35)=8.0,\\ k_r&=6+12(0.35)=10.2. \end{aligned} \]

Substituting these values yields coefficients \( (c+k_v)/m=5 \) and \( (k+k_p)/m=6 \), hence the frozen poles remain \( -2 \) and \( -3 \).

Problem 3 (Interpolation-error bound): Let \( g(\rho)=2+3\rho+4\rho^2 \) on \( [0,0.25] \). Find the linear-interpolation error bound and verify it at the midpoint.

Solution: Since \( g^{\prime\prime}=8 \) and \( h=0.25 \),

\[ |g-g_L|\le\frac{h^2}{8}\max|g^{\prime\prime}| =\frac{0.25^2}{8}(8)=0.0625. \]

At \( \rho=0.125 \), the true value is \( 2.4375 \), whereas the line joining \( g(0)=2 \) and \( g(0.25)=3 \) gives \( 2.5 \). The absolute error is exactly \( 0.0625 \), so the bound is attained.

Problem 4 (Common Lyapunov proof): Assume \( \mathbf{A}_{cl}(\rho)=\sum_j\lambda_j\mathbf{A}_{cl,j} \) and that one \( \mathbf{P}>0 \) satisfies \( \mathbf{A}_{cl,j}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}_{cl,j} \le-\mathbf{Q} \) for all vertices. Prove uniform exponential stability.

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

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

By the comparison lemma,

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

Bounds between \( V \) and \( \|\mathbf{x}\|^2 \) then establish uniform exponential decay of the state.

Problem 5 (Why endpoint interpolation can be non-affine): Let \( \mathbf{B}(\alpha)=(1-\alpha)\mathbf{B}_0+\alpha\mathbf{B}_1 \) and \( \mathbf{K}(\alpha)=(1-\alpha)\mathbf{K}_0+\alpha\mathbf{K}_1 \). Expand their product and identify the cross terms.

Solution:

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

The mixed products are not contained in a simple interpolation of \( \mathbf{B}_0\mathbf{K}_0 \) and \( \mathbf{B}_1\mathbf{K}_1 \). Therefore, endpoint closed-loop matrices alone may not describe the interior family.

Problem 6 (Breakpoint continuity): A three-point schedule uses breakpoints \( \rho_0<\rho_1<\rho_2 \). Show that piecewise-linear interpolation is continuous at \( \rho_1 \) when both adjacent intervals use the same stored controller \( \mathbf{K}_1 \).

Solution: On the left interval, taking \( \rho\rightarrow\rho_1 \) gives \( \alpha\rightarrow1 \) and hence \( \mathbf{K}\rightarrow\mathbf{K}_1 \). On the right interval, taking \( \rho\rightarrow\rho_1 \) gives the new local coordinate \( \alpha\rightarrow0 \), again yielding \( \mathbf{K}\rightarrow\mathbf{K}_1 \). The one-sided limits and stored value coincide. The slopes need not coincide, so differentiability is not guaranteed.

Simple LPV-like controller interpolation converts a table of local linear designs into an online parameter-dependent controller. Linear and barycentric weights provide continuous convex blending, while the interpolation remainder quantifies the effect of grid spacing and gain-map curvature. The principal theoretical warning is that stable local controllers do not by themselves guarantee a stable scheduled loop. A common quadratic Lyapunov matrix supplies a clean sufficient certificate when the closed-loop family is genuinely affine in the interpolation weights. Controller-state alignment, scheduling noise and delay, rate effects, saturation, and nonlinear-model validation remain essential implementation concerns.