Earlier chapters introduced recursive least squares, projection, persistent excitation, gain scheduling, LPV ideas, constraints, and implementation safeguards. This lesson assembles those ingredients around an optimization-based controller. Students are assumed to know linear state-space control and the basic meaning of a constrained optimization problem, but no complete MPC theory is assumed.

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

• distinguish certainty-equivalent Adaptive MPC, robust Adaptive MPC, LPV MPC, and robust LPV MPC;
• write an LPV prediction model with measured or predicted scheduling variables;
• couple a projected RLS estimator to a receding-horizon controller;
• state recursive-feasibility and stability conditions at a conceptual but mathematically precise level;
• explain why parameter convergence, constraint satisfaction, and tracking performance are separate objectives;
• implement a compact adaptive LPV-MPC simulation in five programming environments.

Consider the uncertain discrete-time plant

\[ x_{k+1}=A(\theta^\star)x_k+B(\theta^\star)u_k+w_k, \qquad x_k\in\mathcal{X},\quad u_k\in\mathcal{U}, \]

where \( \theta^\star \) is an unknown constant parameter vector and \( w_k \) is a bounded disturbance. A fixed-model MPC controller predicts with one nominal parameter \( \bar{\theta} \). An Adaptive MPC controller updates a point estimate \( \hat{\theta}_k \), an uncertainty set \( \Theta_k \), or both.

The certainty-equivalent finite-horizon problem has the generic form

\[ \begin{aligned} V_N(x_k,\hat{\theta}_k) = \min_{\mathbf{u}_k} &\sum_{i=0}^{N-1} \ell(x_{i|k},u_{i|k}) +V_f(x_{N|k}) \\ \text{subject to}\quad &x_{0|k}=x_k,\\ &x_{i+1|k} = A(\hat{\theta}_k)x_{i|k} +B(\hat{\theta}_k)u_{i|k},\\ &x_{i|k}\in\mathcal{X},\quad u_{i|k}\in\mathcal{U},\\ &x_{N|k}\in\mathcal{X}_f . \end{aligned} \]

Only \(u_{0|k}^\star\) is applied. At sample \(k+1\), new measurements update the model and the optimization is solved again. This is the receding-horizon principle.

An LPV model is linear in the state and input for a fixed scheduling trajectory, while its matrices vary with a measurable scheduling variable \(\rho_k\):

\[ x_{k+1}=A(\rho_k)x_k+B(\rho_k)u_k. \]

A common affine representation is

\[ A(\rho)=A_0+\sum_{j=1}^{n_\rho}\rho_jA_j, \qquad B(\rho)=B_0+\sum_{j=1}^{n_\rho}\rho_jB_j. \]

If \(\rho\) lies in a polytope, the same model can be embedded as a convex combination of vertex systems:

\[ \begin{aligned} A(\rho)&=\sum_{v=1}^{N_v}\alpha_v(\rho)A_v,\\ B(\rho)&=\sum_{v=1}^{N_v}\alpha_v(\rho)B_v,\\ \alpha_v(\rho)&\geq 0,\qquad \sum_{v=1}^{N_v}\alpha_v(\rho)=1. \end{aligned} \]

LPV MPC requires an assumption about future scheduling values. Three standard choices are:

1. Frozen scheduling: \(\rho_{i|k}=\rho_k\) for the entire horizon.
2. Anticipative scheduling: future \(\rho_{i|k}\) values are measured, commanded, or forecast.
3. Set-valued scheduling: \(\rho_{i|k}\in\mathcal{P}_{i|k}\), often with a bounded rate \(\|\rho_{k+1}-\rho_k\|\leq\Delta_\rho\).

Gain scheduling uses \(\rho_k\) to select or interpolate a predesigned controller. LPV MPC instead places the scheduled model inside the online prediction and optimization problem. Neither construction is automatically adaptive: adaptation requires online learning or uncertainty-set updating.

A model may be simultaneously LPV and adaptive. For example,

\[ x_{k+1} = \left(A_0+\sum_{j=1}^{n_\rho}\rho_{j,k}A_j\right)x_k +Bu_k+w_k, \]

where some entries of \(A_j\) and \(B\) are unknown. Vectorization produces a regression

\[ y_k=\phi_k^\mathsf{T}\theta^\star+w_k. \]

The scalar implementation later in this lesson uses

\[ x_{k+1} = (\theta_0^\star+\theta_1^\star\rho_k)x_k +\theta_2^\star u_k+w_k, \]

\[ \phi_k= \begin{bmatrix} x_k & \rho_kx_k & u_k \end{bmatrix}^{\mathsf{T}}, \qquad x_{k+1}=\phi_k^\mathsf{T}\theta^\star+w_k. \]

Reusing Chapter 12, projected RLS is

\[ \begin{aligned} K_k&= \frac{P_k\phi_k} {\lambda+\phi_k^\mathsf{T}P_k\phi_k},\\ \tilde{y}_k&= x_{k+1}-\phi_k^\mathsf{T}\hat{\theta}_k,\\ \hat{\theta}_{k+1}&= \operatorname{Proj}_{\Theta} \left(\hat{\theta}_k+K_k\tilde{y}_k\right),\\ P_{k+1}&= \lambda^{-1} \left(P_k-K_k\phi_k^\mathsf{T}P_k\right). \end{aligned} \]

Projection encodes prior physical bounds and prevents an optimizer from receiving implausible matrices. The forgetting factor \(0<\lambda\leq 1\) determines how quickly old data are discounted.

Parameter adaptation and LPV scheduling play different roles: \(\rho_k\) is assumed measurable and explains known operating-point variation; \(\hat{\theta}_k\) estimates uncertain coefficients that remain unknown after scheduling is included.

With a scheduling prediction \(\boldsymbol{\rho}_k= \{\rho_{0|k},\ldots,\rho_{N-1|k}\}\), define

\[ A_{i|k}=A(\rho_{i|k},\hat{\theta}_k), \qquad B_{i|k}=B(\rho_{i|k},\hat{\theta}_k). \]

The quadratic tracking problem is

\[ \begin{aligned} \min_{\mathbf{u}_k,\boldsymbol{\epsilon}_k}\quad &\sum_{i=0}^{N-1} \left( \|x_{i|k}-r_{i|k}\|_Q^2+ \|u_{i|k}\|_R^2+ \mu\|\epsilon_{i|k}\|_1 \right)\\ &+\|x_{N|k}-r_{N|k}\|_P^2\\ \text{subject to}\quad &x_{i+1|k}=A_{i|k}x_{i|k}+B_{i|k}u_{i|k},\\ &x_{i|k}\in\mathcal{X}\oplus\{\epsilon_{i|k}\},\\ &u_{i|k}\in\mathcal{U},\qquad \epsilon_{i|k}\geq 0. \end{aligned} \]

Slack variables can preserve numerical feasibility, but they convert hard state constraints into soft constraints. Safety-critical constraints should not be softened unless a separate certified protection layer is present.

For a fixed scheduling sequence and affine dynamics, the condensed prediction is

\[ \mathbf{x}_k= \mathcal{A}_k x_k+\mathcal{B}_k\mathbf{u}_k, \]

where the block rows contain products such as

\[ x_{2|k} = A_{1|k}A_{0|k}x_k+ A_{1|k}B_{0|k}u_{0|k}+ B_{1|k}u_{1|k}. \]

Substitution yields a quadratic program when the constraints are linear. Nonlinear dependence on scheduling or parameters may instead produce a nonlinear program.

6.1 Certainty-equivalent Adaptive MPC

The controller predicts only with \(\hat{\theta}_k\). This is computationally attractive and often performs well when estimation error is small, but constraint guarantees generally depend on additional robustness assumptions.

6.2 Robust Adaptive MPC

Let \(\Theta_k\) contain all parameters consistent with the prior set, measurements, and a disturbance bound:

\[ \Theta_{k+1} = \Theta_k\cap \left\{ \theta: \|x_{k+1}-\Phi_k\theta\|\leq\bar{w} \right\}. \]

Therefore,

\[ \theta^\star\in\Theta_{k+1}\subseteq\Theta_k. \]

A tube or min-max MPC formulation can enforce constraints for every admissible parameter and disturbance. The nominal trajectory \(z_{i|k}\) is surrounded by an error tube \(\mathcal{S}_{i|k}\):

\[ x_{i|k}\in z_{i|k}\oplus\mathcal{S}_{i|k}. \]

The nominal constraints are tightened:

\[ z_{i|k}\in\mathcal{X}\ominus\mathcal{S}_{i|k}, \qquad v_{i|k}\in\mathcal{U}\ominus K\mathcal{S}_{i|k}. \]

As \(\Theta_k\) contracts, the tube may shrink and the controller becomes less conservative.

6.3 Dual control effect

An input influences both the physical state and the information available to the estimator. A purely tracking-oriented input may fail to excite the uncertain dynamics. A dual or exploration-aware cost may add

\[ J_{\mathrm{dual}} = J_{\mathrm{track}} +\gamma\,\Psi(P_{k+N|k}), \]

where \(\Psi\) penalizes predicted uncertainty, for example \(\operatorname{tr}(P)\) or \(\log\det(P)\). This improves identification only if excitation can be introduced without violating constraints.

7.1 Fixed-model nominal MPC argument

Suppose the terminal set \(\mathcal{X}_f\) is invariant under a local controller \(\kappa_f(x)\) and satisfies

\[ V_f(f(x,\kappa_f(x)))-V_f(x) \leq -\ell(x,\kappa_f(x)). \]

Shift the optimal sequence at time \(k\) and append the terminal action:

\[ \tilde{\mathbf{u}}_{k+1} = \left\{ u_{1|k}^\star,\ldots, u_{N-1|k}^\star, \kappa_f(x_{N|k}^\star) \right\}. \]

Feasibility of the shifted sequence gives

\[ V_N(x_{k+1})-V_N(x_k) \leq -\ell(x_k,u_k). \]

Thus the optimal value acts as a Lyapunov function.

7.2 What adaptation changes

When the prediction model changes from \(\hat{\theta}_k\) to \(\hat{\theta}_{k+1}\), the shifted sequence may no longer generate the same predicted trajectory. Under a local Lipschitz bound on the value function,

\[ \begin{aligned} &V_N(x_{k+1},\hat{\theta}_{k+1}) -V_N(x_k,\hat{\theta}_k)\\ &\qquad\leq -\ell(x_k,u_k) +L_\theta \|\hat{\theta}_{k+1}-\hat{\theta}_k\| +c_w\|w_k\|. \end{aligned} \]

Consequently, exact asymptotic stability is not obtained merely because the nominal MPC problem is stabilizing. Typical sufficient routes are:

• summable or convergent model updates, \(\sum_k\|\Delta\hat{\theta}_k\|<\infty\);
• robust tube bounds that absorb model mismatch and disturbances;
• a common Lyapunov function over all admissible LPV vertices;
• terminal ingredients parameterized consistently with \(\rho\) and \(\Theta_k\);
• input-to-state stability or practical stability instead of an unjustified asymptotic claim.

7.3 Why nested uncertainty sets help feasibility

Suppose a robust MPC problem at time \(k\) enforces its constraints for every \(\theta\in\Theta_k\). If the estimator guarantees \(\Theta_{k+1}\subseteq\Theta_k\), every parameter considered at the next step was already included in the previous robust design. Together with a robustly invariant terminal construction, this monotonicity is a central ingredient in recursive-feasibility proofs.

flowchart TD
  R["Reference and constraints"] --> O["Finite-horizon optimizer"]
  X["Measured state x(k)"] --> O
  P["Measured or predicted schedule rho"] --> O
  E["Projected estimator: theta_hat and uncertainty set"] --> O
  O --> U["Apply first input u(k)"]
  U --> G["Uncertain LPV plant"]
  G --> XN["New state x(k+1)"]
  XN --> E
  XN --> X
  P --> E
  U --> E
        

The estimator and optimizer form coupled dynamics. The estimator changes the prediction model; the optimizer changes the input and therefore the informativeness of future data.

The implementations use the scalar model

\[ x_{k+1} = (\theta_0+\theta_1\rho_k)x_k+\theta_2u_k+w_k. \]

At every sample:

1. obtain the measured state and a horizon of scheduling values;
2. predict with the current projected RLS estimate;
3. minimize a quadratic tracking cost by projected gradient descent;
4. project each candidate input onto \([-u_{\max},u_{\max}]\);
5. apply the first input;
6. update the parameter estimate with the new transition;
7. shift the previous optimizer solution to warm-start the next problem.

The demonstration uses a quadratic penalty for state-constraint violation. Therefore, the state limit is soft; this code is not a certificate of robust constraint satisfaction. Its purpose is to expose the estimator-predictor-optimizer coupling without hiding the mechanism behind a large software framework.

For the scalar prediction, define \(a_i=\hat{\theta}_0+\hat{\theta}_1\rho_{i|k}\) and \(b=\hat{\theta}_2\). The adjoint recursion used to compute the cost gradient is

\[ \begin{aligned} \lambda_N &= 2P(x_{N|k}-r_N)+ \nabla p_{\mathcal{X}}(x_{N|k}),\\ \frac{\partial J}{\partial u_{i|k}} &= 2Ru_{i|k}+b\lambda_{i+1},\\ \lambda_i &= 2Q(x_{i|k}-r_i)+ \nabla p_{\mathcal{X}}(x_{i|k})+ a_i\lambda_{i+1}. \end{aligned} \]

The projected-gradient iteration is

\[ \mathbf{u}^{(j+1)} = \Pi_{\mathcal{U}^N} \left( \mathbf{u}^{(j)} -\alpha_j\nabla_{\mathbf{u}}J \right). \]

The example requires NumPy. For larger constrained problems, common choices include CVXPY with OSQP, CasADi, do-mpc, GEKKO, or direct solver interfaces.

Chapter28_Lesson4.py

"""
Chapter28_Lesson4.py
Adaptive LPV model predictive control for a scalar constrained plant.

Model:
    x[k+1] = (theta0 + theta1*rho[k]) * x[k] + theta2*u[k] + w[k]

The controller:
1. estimates theta online with projected recursive least squares;
2. predicts the measured scheduling variable rho over the MPC horizon;
3. solves a box-constrained finite-horizon problem by projected gradient descent;
4. applies only the first optimized input (receding-horizon control).

Dependencies:
    numpy
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Tuple

import numpy as np


@dataclass(frozen=True)
class MPCConfig:
    horizon: int = 12
    q: float = 5.0
    r: float = 0.15
    p_terminal: float = 8.0
    u_max: float = 1.5
    x_max: float = 3.0
    state_penalty: float = 80.0
    step_size: float = 0.012
    iterations: int = 90


class ProjectedRLS:
    """Recursive least squares with forgetting and componentwise projection."""

    def __init__(
        self,
        theta0: np.ndarray,
        covariance_scale: float = 80.0,
        forgetting: float = 0.995,
        lower: np.ndarray | None = None,
        upper: np.ndarray | None = None,
    ) -> None:
        self.theta = np.asarray(theta0, dtype=float).copy()
        self.P = covariance_scale * np.eye(self.theta.size)
        self.forgetting = float(forgetting)
        self.lower = (
            np.full_like(self.theta, -np.inf)
            if lower is None
            else np.asarray(lower, dtype=float)
        )
        self.upper = (
            np.full_like(self.theta, np.inf)
            if upper is None
            else np.asarray(upper, dtype=float)
        )

    def update(self, phi: np.ndarray, y: float) -> Tuple[np.ndarray, float]:
        phi = np.asarray(phi, dtype=float)
        denominator = self.forgetting + phi @ self.P @ phi
        gain = (self.P @ phi) / denominator
        prediction_error = float(y - phi @ self.theta)

        self.theta = self.theta + gain * prediction_error
        self.theta = np.clip(self.theta, self.lower, self.upper)

        self.P = (self.P - np.outer(gain, phi) @ self.P) / self.forgetting
        self.P = 0.5 * (self.P + self.P.T)
        return self.theta.copy(), prediction_error


def scheduling_signal(k: int) -> float:
    """Measured exogenous scheduling variable in [-1, 1]."""
    return float(0.75 * np.sin(0.045 * k) + 0.20 * np.sin(0.11 * k))


def reference_signal(k: int) -> float:
    if k < 25:
        return 0.0
    if k < 95:
        return 1.0
    if k < 145:
        return -0.7
    return 0.4


def state_penalty_derivative(x: float, x_max: float, weight: float) -> float:
    excess = abs(x) - x_max
    if excess <= 0.0:
        return 0.0
    return 2.0 * weight * excess * np.sign(x)


def predict_states(
    x0: float,
    u_sequence: np.ndarray,
    rho_sequence: np.ndarray,
    theta: np.ndarray,
) -> np.ndarray:
    horizon = u_sequence.size
    x = np.empty(horizon + 1, dtype=float)
    x[0] = x0
    b_hat = float(theta[2])

    for i in range(horizon):
        a_hat = float(theta[0] + theta[1] * rho_sequence[i])
        x[i + 1] = a_hat * x[i] + b_hat * u_sequence[i]
    return x


def objective_and_gradient(
    x0: float,
    u_sequence: np.ndarray,
    rho_sequence: np.ndarray,
    reference: float,
    theta: np.ndarray,
    cfg: MPCConfig,
) -> Tuple[float, np.ndarray, np.ndarray]:
    """Return finite-horizon cost, its input gradient, and predicted states."""
    x = predict_states(x0, u_sequence, rho_sequence, theta)
    horizon = cfg.horizon

    cost = 0.0
    for i in range(horizon):
        state_error = x[i] - reference
        excess = max(abs(x[i]) - cfg.x_max, 0.0)
        cost += (
            cfg.q * state_error**2
            + cfg.r * u_sequence[i] ** 2
            + cfg.state_penalty * excess**2
        )

    terminal_error = x[horizon] - reference
    terminal_excess = max(abs(x[horizon]) - cfg.x_max, 0.0)
    cost += (
        cfg.p_terminal * terminal_error**2
        + cfg.state_penalty * terminal_excess**2
    )

    gradient = np.zeros(horizon, dtype=float)
    adjoint = (
        2.0 * cfg.p_terminal * terminal_error
        + state_penalty_derivative(
            x[horizon], cfg.x_max, cfg.state_penalty
        )
    )

    b_hat = float(theta[2])
    for i in range(horizon - 1, -1, -1):
        gradient[i] = 2.0 * cfg.r * u_sequence[i] + b_hat * adjoint

        a_hat = float(theta[0] + theta[1] * rho_sequence[i])
        stage_derivative = (
            2.0 * cfg.q * (x[i] - reference)
            + state_penalty_derivative(
                x[i], cfg.x_max, cfg.state_penalty
            )
        )
        adjoint = stage_derivative + a_hat * adjoint

    return float(cost), gradient, x


def solve_adaptive_lpv_mpc(
    x0: float,
    rho_sequence: np.ndarray,
    reference: float,
    theta: np.ndarray,
    warm_start: np.ndarray,
    cfg: MPCConfig,
) -> Tuple[np.ndarray, np.ndarray, float]:
    """Projected-gradient MPC with warm starting and box input constraints."""
    u = np.clip(warm_start.copy(), -cfg.u_max, cfg.u_max)
    best_u = u.copy()
    best_cost = np.inf
    best_x = predict_states(x0, u, rho_sequence, theta)

    for iteration in range(cfg.iterations):
        cost, gradient, x_prediction = objective_and_gradient(
            x0, u, rho_sequence, reference, theta, cfg
        )

        if cost < best_cost:
            best_cost = cost
            best_u = u.copy()
            best_x = x_prediction.copy()

        diminishing_step = cfg.step_size / np.sqrt(1.0 + 0.03 * iteration)
        u = np.clip(
            u - diminishing_step * gradient,
            -cfg.u_max,
            cfg.u_max,
        )

    return best_u, best_x, float(best_cost)


def main() -> None:
    rng = np.random.default_rng(7)
    cfg = MPCConfig()

    theta_true = np.array([0.72, 0.18, 0.55], dtype=float)
    estimator = ProjectedRLS(
        theta0=np.array([0.55, 0.02, 0.35], dtype=float),
        covariance_scale=60.0,
        forgetting=0.995,
        lower=np.array([0.35, -0.05, 0.15]),
        upper=np.array([1.05, 0.40, 0.90]),
    )

    steps = 190
    x = 0.0
    warm_start = np.zeros(cfg.horizon, dtype=float)

    x_log = np.zeros(steps + 1)
    u_log = np.zeros(steps)
    rho_log = np.zeros(steps)
    reference_log = np.zeros(steps)
    theta_log = np.zeros((steps, 3))
    prediction_error_log = np.zeros(steps)

    for k in range(steps):
        rho_horizon = np.array(
            [scheduling_signal(k + i) for i in range(cfg.horizon)],
            dtype=float,
        )
        reference = reference_signal(k)

        optimal_u, _, _ = solve_adaptive_lpv_mpc(
            x0=x,
            rho_sequence=rho_horizon,
            reference=reference,
            theta=estimator.theta,
            warm_start=warm_start,
            cfg=cfg,
        )
        u = float(optimal_u[0])

        rho = scheduling_signal(k)
        a_true = theta_true[0] + theta_true[1] * rho
        process_noise = 0.006 * rng.standard_normal()
        x_next = a_true * x + theta_true[2] * u + process_noise

        phi = np.array([x, rho * x, u], dtype=float)
        theta_hat, prediction_error = estimator.update(phi, x_next)

        x_log[k] = x
        u_log[k] = u
        rho_log[k] = rho
        reference_log[k] = reference
        theta_log[k] = theta_hat
        prediction_error_log[k] = prediction_error

        warm_start[:-1] = optimal_u[1:]
        warm_start[-1] = optimal_u[-1]
        x = float(x_next)

    x_log[-1] = x

    tracking_rmse = np.sqrt(
        np.mean((x_log[:-1] - reference_log) ** 2)
    )
    maximum_input = np.max(np.abs(u_log))
    maximum_state = np.max(np.abs(x_log))
    final_parameter_error = np.linalg.norm(estimator.theta - theta_true)

    print("Adaptive LPV-MPC simulation completed.")
    print(f"Tracking RMSE: {tracking_rmse:.5f}")
    print(f"Maximum |u|:  {maximum_input:.5f}")
    print(f"Maximum |x|:  {maximum_state:.5f}")
    print("True theta:  ", theta_true)
    print("Final theta: ", estimator.theta)
    print(f"Final parameter error norm: {final_parameter_error:.5f}")
    print(
        "Mean absolute one-step prediction error: "
        f"{np.mean(np.abs(prediction_error_log)):.5f}"
    )


if __name__ == "__main__":
    main()

The code is self-contained C++17. Production implementations commonly use Eigen for linear algebra and OSQP, qpOASES, HPIPM, acados, or CasADi generated code for optimization.

Chapter28_Lesson4.cpp

/*
Chapter28_Lesson4.cpp
Adaptive LPV model predictive control for a scalar constrained plant.

Build:
    g++ -std=c++17 -O2 Chapter28_Lesson4.cpp -o Chapter28_Lesson4
Run:
    ./Chapter28_Lesson4
*/

#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <random>
#include <vector>

struct MPCConfig {
    int horizon = 12;
    double q = 5.0;
    double r = 0.15;
    double p_terminal = 8.0;
    double u_max = 1.5;
    double x_max = 3.0;
    double state_penalty = 80.0;
    double step_size = 0.012;
    int iterations = 90;
};

class ProjectedRLS {
public:
    ProjectedRLS(
        const std::array<double, 3>& theta0,
        double covariance_scale,
        double forgetting,
        const std::array<double, 3>& lower,
        const std::array<double, 3>& upper)
        : theta_(theta0),
          forgetting_(forgetting),
          lower_(lower),
          upper_(upper) {
        for (auto& row : P_) {
            row.fill(0.0);
        }
        for (int i = 0; i < 3; ++i) {
            P_[i][i] = covariance_scale;
        }
    }

    std::pair<std::array<double, 3>, double> update(
        const std::array<double, 3>& phi,
        double y) {
        std::array<double, 3> Pphi{};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                Pphi[i] += P_[i][j] * phi[j];
            }
        }

        double denominator = forgetting_;
        for (int i = 0; i < 3; ++i) {
            denominator += phi[i] * Pphi[i];
        }

        std::array<double, 3> gain{};
        for (int i = 0; i < 3; ++i) {
            gain[i] = Pphi[i] / denominator;
        }

        double prediction = 0.0;
        for (int i = 0; i < 3; ++i) {
            prediction += phi[i] * theta_[i];
        }
        const double error = y - prediction;

        for (int i = 0; i < 3; ++i) {
            theta_[i] += gain[i] * error;
            theta_[i] = std::clamp(theta_[i], lower_[i], upper_[i]);
        }

        std::array<std::array<double, 3>, 3> newP{};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                double correction = 0.0;
                for (int m = 0; m < 3; ++m) {
                    correction += gain[i] * phi[m] * P_[m][j];
                }
                newP[i][j] = (P_[i][j] - correction) / forgetting_;
            }
        }

        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                P_[i][j] = 0.5 * (newP[i][j] + newP[j][i]);
            }
        }

        return {theta_, error};
    }

    const std::array<double, 3>& theta() const {
        return theta_;
    }

private:
    std::array<double, 3> theta_{};
    std::array<std::array<double, 3>, 3> P_{};
    double forgetting_;
    std::array<double, 3> lower_{};
    std::array<double, 3> upper_{};
};

double schedulingSignal(int k) {
    return 0.75 * std::sin(0.045 * k) + 0.20 * std::sin(0.11 * k);
}

double referenceSignal(int k) {
    if (k < 25) return 0.0;
    if (k < 95) return 1.0;
    if (k < 145) return -0.7;
    return 0.4;
}

double statePenaltyDerivative(double x, const MPCConfig& cfg) {
    const double excess = std::abs(x) - cfg.x_max;
    if (excess <= 0.0) {
        return 0.0;
    }
    const double sign = (x >= 0.0) ? 1.0 : -1.0;
    return 2.0 * cfg.state_penalty * excess * sign;
}

std::vector<double> predictStates(
    double x0,
    const std::vector<double>& u,
    const std::vector<double>& rho,
    const std::array<double, 3>& theta) {
    std::vector<double> x(u.size() + 1, 0.0);
    x[0] = x0;
    const double bHat = theta[2];

    for (std::size_t i = 0; i < u.size(); ++i) {
        const double aHat = theta[0] + theta[1] * rho[i];
        x[i + 1] = aHat * x[i] + bHat * u[i];
    }
    return x;
}

struct ObjectiveResult {
    double cost;
    std::vector<double> gradient;
    std::vector<double> states;
};

ObjectiveResult objectiveAndGradient(
    double x0,
    const std::vector<double>& u,
    const std::vector<double>& rho,
    double reference,
    const std::array<double, 3>& theta,
    const MPCConfig& cfg) {
    const int N = cfg.horizon;
    auto x = predictStates(x0, u, rho, theta);
    double cost = 0.0;

    for (int i = 0; i < N; ++i) {
        const double error = x[i] - reference;
        const double excess = std::max(std::abs(x[i]) - cfg.x_max, 0.0);
        cost += cfg.q * error * error
              + cfg.r * u[i] * u[i]
              + cfg.state_penalty * excess * excess;
    }

    const double terminalError = x[N] - reference;
    const double terminalExcess =
        std::max(std::abs(x[N]) - cfg.x_max, 0.0);
    cost += cfg.p_terminal * terminalError * terminalError
          + cfg.state_penalty * terminalExcess * terminalExcess;

    std::vector<double> gradient(N, 0.0);
    double adjoint =
        2.0 * cfg.p_terminal * terminalError
        + statePenaltyDerivative(x[N], cfg);

    const double bHat = theta[2];
    for (int i = N - 1; i >= 0; --i) {
        gradient[i] = 2.0 * cfg.r * u[i] + bHat * adjoint;

        const double aHat = theta[0] + theta[1] * rho[i];
        const double stageDerivative =
            2.0 * cfg.q * (x[i] - reference)
            + statePenaltyDerivative(x[i], cfg);
        adjoint = stageDerivative + aHat * adjoint;
    }

    return {cost, gradient, x};
}

struct MPCResult {
    std::vector<double> input;
    std::vector<double> statePrediction;
    double cost;
};

MPCResult solveAdaptiveLPVMPC(
    double x0,
    const std::vector<double>& rho,
    double reference,
    const std::array<double, 3>& theta,
    const std::vector<double>& warmStart,
    const MPCConfig& cfg) {
    std::vector<double> u = warmStart;
    for (double& value : u) {
        value = std::clamp(value, -cfg.u_max, cfg.u_max);
    }

    std::vector<double> bestU = u;
    std::vector<double> bestX = predictStates(x0, u, rho, theta);
    double bestCost = std::numeric_limits<double>::infinity();

    for (int iteration = 0; iteration < cfg.iterations; ++iteration) {
        const auto result =
            objectiveAndGradient(x0, u, rho, reference, theta, cfg);

        if (result.cost < bestCost) {
            bestCost = result.cost;
            bestU = u;
            bestX = result.states;
        }

        const double step =
            cfg.step_size / std::sqrt(1.0 + 0.03 * iteration);
        for (int i = 0; i < cfg.horizon; ++i) {
            u[i] = std::clamp(
                u[i] - step * result.gradient[i],
                -cfg.u_max,
                cfg.u_max);
        }
    }

    return {bestU, bestX, bestCost};
}

int main() {
    const MPCConfig cfg;
    const std::array<double, 3> thetaTrue{0.72, 0.18, 0.55};

    ProjectedRLS estimator(
        {0.55, 0.02, 0.35},
        60.0,
        0.995,
        {0.35, -0.05, 0.15},
        {1.05, 0.40, 0.90});

    std::mt19937 generator(7);
    std::normal_distribution<double> normal(0.0, 1.0);

    const int steps = 190;
    double x = 0.0;
    std::vector<double> warmStart(cfg.horizon, 0.0);

    double squaredTrackingError = 0.0;
    double maxInput = 0.0;
    double maxState = 0.0;
    double absolutePredictionError = 0.0;

    for (int k = 0; k < steps; ++k) {
        std::vector<double> rhoHorizon(cfg.horizon, 0.0);
        for (int i = 0; i < cfg.horizon; ++i) {
            rhoHorizon[i] = schedulingSignal(k + i);
        }

        const double reference = referenceSignal(k);
        const auto result = solveAdaptiveLPVMPC(
            x,
            rhoHorizon,
            reference,
            estimator.theta(),
            warmStart,
            cfg);

        const double u = result.input.front();
        const double rho = schedulingSignal(k);
        const double aTrue = thetaTrue[0] + thetaTrue[1] * rho;
        const double processNoise = 0.006 * normal(generator);
        const double xNext =
            aTrue * x + thetaTrue[2] * u + processNoise;

        const std::array<double, 3> phi{x, rho * x, u};
        const auto update = estimator.update(phi, xNext);

        squaredTrackingError += (x - reference) * (x - reference);
        maxInput = std::max(maxInput, std::abs(u));
        maxState = std::max(maxState, std::abs(x));
        absolutePredictionError += std::abs(update.second);

        for (int i = 0; i < cfg.horizon - 1; ++i) {
            warmStart[i] = result.input[i + 1];
        }
        warmStart.back() = result.input.back();
        x = xNext;
    }

    maxState = std::max(maxState, std::abs(x));
    const auto thetaFinal = estimator.theta();
    double parameterErrorSquared = 0.0;
    for (int i = 0; i < 3; ++i) {
        const double difference = thetaFinal[i] - thetaTrue[i];
        parameterErrorSquared += difference * difference;
    }

    std::cout << std::fixed << std::setprecision(5);
    std::cout << "Adaptive LPV-MPC simulation completed.\n";
    std::cout << "Tracking RMSE: "
              << std::sqrt(squaredTrackingError / steps) << "\n";
    std::cout << "Maximum |u|:  " << maxInput << "\n";
    std::cout << "Maximum |x|:  " << maxState << "\n";
    std::cout << "True theta:   ["
              << thetaTrue[0] << ", "
              << thetaTrue[1] << ", "
              << thetaTrue[2] << "]\n";
    std::cout << "Final theta:  ["
              << thetaFinal[0] << ", "
              << thetaFinal[1] << ", "
              << thetaFinal[2] << "]\n";
    std::cout << "Final parameter error norm: "
              << std::sqrt(parameterErrorSquared) << "\n";
    std::cout << "Mean absolute one-step prediction error: "
              << absolutePredictionError / steps << "\n";

    return 0;
}

The Java version is dependency-free. EJML can replace manual matrix operations, while ojAlgo provides quadratic and convex optimization facilities for larger MPC formulations.

Chapter28_Lesson4.java

/*
Chapter28_Lesson4.java
Adaptive LPV model predictive control for a scalar constrained plant.

Build:
    javac Chapter28_Lesson4.java
Run:
    java Chapter28_Lesson4
*/

import java.util.Arrays;
import java.util.Random;

public final class Chapter28_Lesson4 {
    private Chapter28_Lesson4() {
    }

    private static final class MPCConfig {
        final int horizon = 12;
        final double q = 5.0;
        final double r = 0.15;
        final double pTerminal = 8.0;
        final double uMax = 1.5;
        final double xMax = 3.0;
        final double statePenalty = 80.0;
        final double stepSize = 0.012;
        final int iterations = 90;
    }

    private static final class RLSUpdate {
        final double[] theta;
        final double predictionError;

        RLSUpdate(double[] theta, double predictionError) {
            this.theta = theta;
            this.predictionError = predictionError;
        }
    }

    private static final class ProjectedRLS {
        private final double[] theta;
        private final double[][] covariance;
        private final double forgetting;
        private final double[] lower;
        private final double[] upper;

        ProjectedRLS(
                double[] theta0,
                double covarianceScale,
                double forgetting,
                double[] lower,
                double[] upper) {
            this.theta = theta0.clone();
            this.covariance = new double[3][3];
            this.forgetting = forgetting;
            this.lower = lower.clone();
            this.upper = upper.clone();

            for (int i = 0; i < 3; i++) {
                covariance[i][i] = covarianceScale;
            }
        }

        RLSUpdate update(double[] phi, double y) {
            double[] pPhi = new double[3];
            for (int i = 0; i < 3; i++) {
                for (int j = 0; j < 3; j++) {
                    pPhi[i] += covariance[i][j] * phi[j];
                }
            }

            double denominator = forgetting;
            for (int i = 0; i < 3; i++) {
                denominator += phi[i] * pPhi[i];
            }

            double[] gain = new double[3];
            for (int i = 0; i < 3; i++) {
                gain[i] = pPhi[i] / denominator;
            }

            double prediction = 0.0;
            for (int i = 0; i < 3; i++) {
                prediction += phi[i] * theta[i];
            }
            double error = y - prediction;

            for (int i = 0; i < 3; i++) {
                theta[i] = clamp(
                        theta[i] + gain[i] * error,
                        lower[i],
                        upper[i]);
            }

            double[][] updated = new double[3][3];
            for (int i = 0; i < 3; i++) {
                for (int j = 0; j < 3; j++) {
                    double correction = 0.0;
                    for (int m = 0; m < 3; m++) {
                        correction +=
                                gain[i] * phi[m] * covariance[m][j];
                    }
                    updated[i][j] =
                            (covariance[i][j] - correction) / forgetting;
                }
            }

            for (int i = 0; i < 3; i++) {
                for (int j = 0; j < 3; j++) {
                    covariance[i][j] =
                            0.5 * (updated[i][j] + updated[j][i]);
                }
            }

            return new RLSUpdate(theta.clone(), error);
        }

        double[] theta() {
            return theta.clone();
        }
    }

    private static final class ObjectiveResult {
        final double cost;
        final double[] gradient;
        final double[] states;

        ObjectiveResult(double cost, double[] gradient, double[] states) {
            this.cost = cost;
            this.gradient = gradient;
            this.states = states;
        }
    }

    private static final class MPCResult {
        final double[] input;
        final double[] statePrediction;
        final double cost;

        MPCResult(double[] input, double[] statePrediction, double cost) {
            this.input = input;
            this.statePrediction = statePrediction;
            this.cost = cost;
        }
    }

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

    private static double schedulingSignal(int k) {
        return 0.75 * Math.sin(0.045 * k)
                + 0.20 * Math.sin(0.11 * k);
    }

    private static double referenceSignal(int k) {
        if (k < 25) {
            return 0.0;
        }
        if (k < 95) {
            return 1.0;
        }
        if (k < 145) {
            return -0.7;
        }
        return 0.4;
    }

    private static double statePenaltyDerivative(
            double x,
            MPCConfig cfg) {
        double excess = Math.abs(x) - cfg.xMax;
        if (excess <= 0.0) {
            return 0.0;
        }
        double sign = x >= 0.0 ? 1.0 : -1.0;
        return 2.0 * cfg.statePenalty * excess * sign;
    }

    private static double[] predictStates(
            double x0,
            double[] input,
            double[] rho,
            double[] theta) {
        double[] states = new double[input.length + 1];
        states[0] = x0;
        double bHat = theta[2];

        for (int i = 0; i < input.length; i++) {
            double aHat = theta[0] + theta[1] * rho[i];
            states[i + 1] =
                    aHat * states[i] + bHat * input[i];
        }
        return states;
    }

    private static ObjectiveResult objectiveAndGradient(
            double x0,
            double[] input,
            double[] rho,
            double reference,
            double[] theta,
            MPCConfig cfg) {
        int horizon = cfg.horizon;
        double[] states =
                predictStates(x0, input, rho, theta);
        double cost = 0.0;

        for (int i = 0; i < horizon; i++) {
            double error = states[i] - reference;
            double excess =
                    Math.max(Math.abs(states[i]) - cfg.xMax, 0.0);
            cost += cfg.q * error * error
                    + cfg.r * input[i] * input[i]
                    + cfg.statePenalty * excess * excess;
        }

        double terminalError = states[horizon] - reference;
        double terminalExcess =
                Math.max(
                        Math.abs(states[horizon]) - cfg.xMax,
                        0.0);
        cost += cfg.pTerminal * terminalError * terminalError
                + cfg.statePenalty
                * terminalExcess * terminalExcess;

        double[] gradient = new double[horizon];
        double adjoint =
                2.0 * cfg.pTerminal * terminalError
                + statePenaltyDerivative(
                        states[horizon],
                        cfg);

        double bHat = theta[2];
        for (int i = horizon - 1; i >= 0; i--) {
            gradient[i] =
                    2.0 * cfg.r * input[i] + bHat * adjoint;

            double aHat = theta[0] + theta[1] * rho[i];
            double stageDerivative =
                    2.0 * cfg.q * (states[i] - reference)
                    + statePenaltyDerivative(states[i], cfg);
            adjoint = stageDerivative + aHat * adjoint;
        }

        return new ObjectiveResult(cost, gradient, states);
    }

    private static MPCResult solveAdaptiveLPVMPC(
            double x0,
            double[] rho,
            double reference,
            double[] theta,
            double[] warmStart,
            MPCConfig cfg) {
        double[] input = warmStart.clone();
        for (int i = 0; i < input.length; i++) {
            input[i] =
                    clamp(input[i], -cfg.uMax, cfg.uMax);
        }

        double[] bestInput = input.clone();
        double[] bestStates =
                predictStates(x0, input, rho, theta);
        double bestCost = Double.POSITIVE_INFINITY;

        for (int iteration = 0;
                iteration < cfg.iterations;
                iteration++) {
            ObjectiveResult result =
                    objectiveAndGradient(
                            x0,
                            input,
                            rho,
                            reference,
                            theta,
                            cfg);

            if (result.cost < bestCost) {
                bestCost = result.cost;
                bestInput = input.clone();
                bestStates = result.states.clone();
            }

            double step =
                    cfg.stepSize
                    / Math.sqrt(1.0 + 0.03 * iteration);
            for (int i = 0; i < cfg.horizon; i++) {
                input[i] = clamp(
                        input[i] - step * result.gradient[i],
                        -cfg.uMax,
                        cfg.uMax);
            }
        }

        return new MPCResult(bestInput, bestStates, bestCost);
    }

    public static void main(String[] args) {
        MPCConfig cfg = new MPCConfig();
        double[] thetaTrue = {0.72, 0.18, 0.55};

        ProjectedRLS estimator = new ProjectedRLS(
                new double[] {0.55, 0.02, 0.35},
                60.0,
                0.995,
                new double[] {0.35, -0.05, 0.15},
                new double[] {1.05, 0.40, 0.90});

        Random random = new Random(7L);
        int steps = 190;
        double x = 0.0;
        double[] warmStart = new double[cfg.horizon];

        double squaredTrackingError = 0.0;
        double maxInput = 0.0;
        double maxState = 0.0;
        double absolutePredictionError = 0.0;

        for (int k = 0; k < steps; k++) {
            double[] rhoHorizon = new double[cfg.horizon];
            for (int i = 0; i < cfg.horizon; i++) {
                rhoHorizon[i] = schedulingSignal(k + i);
            }

            double reference = referenceSignal(k);
            MPCResult result = solveAdaptiveLPVMPC(
                    x,
                    rhoHorizon,
                    reference,
                    estimator.theta(),
                    warmStart,
                    cfg);

            double input = result.input[0];
            double rho = schedulingSignal(k);
            double aTrue =
                    thetaTrue[0] + thetaTrue[1] * rho;
            double processNoise =
                    0.006 * random.nextGaussian();
            double xNext =
                    aTrue * x
                    + thetaTrue[2] * input
                    + processNoise;

            double[] phi = {x, rho * x, input};
            RLSUpdate update = estimator.update(phi, xNext);

            squaredTrackingError +=
                    (x - reference) * (x - reference);
            maxInput = Math.max(maxInput, Math.abs(input));
            maxState = Math.max(maxState, Math.abs(x));
            absolutePredictionError +=
                    Math.abs(update.predictionError);

            for (int i = 0; i < cfg.horizon - 1; i++) {
                warmStart[i] = result.input[i + 1];
            }
            warmStart[cfg.horizon - 1] =
                    result.input[cfg.horizon - 1];
            x = xNext;
        }

        maxState = Math.max(maxState, Math.abs(x));
        double[] thetaFinal = estimator.theta();
        double parameterErrorSquared = 0.0;
        for (int i = 0; i < 3; i++) {
            double difference =
                    thetaFinal[i] - thetaTrue[i];
            parameterErrorSquared += difference * difference;
        }

        System.out.println(
                "Adaptive LPV-MPC simulation completed.");
        System.out.printf(
                "Tracking RMSE: %.5f%n",
                Math.sqrt(squaredTrackingError / steps));
        System.out.printf("Maximum |u|:  %.5f%n", maxInput);
        System.out.printf("Maximum |x|:  %.5f%n", maxState);
        System.out.println(
                "True theta:   " + Arrays.toString(thetaTrue));
        System.out.println(
                "Final theta:  " + Arrays.toString(thetaFinal));
        System.out.printf(
                "Final parameter error norm: %.5f%n",
                Math.sqrt(parameterErrorSquared));
        System.out.printf(
                "Mean absolute one-step prediction error: %.5f%n",
                absolutePredictionError / steps);
    }
}

The script needs no toolbox. In MATLAB projects, the same structure can use quadprog, fmincon, the Model Predictive Control Toolbox, YALMIP, or CasADi. In Simulink, place the estimator in a discrete subsystem and execute the optimizer in a MATLAB Function block, an S-Function, or generated solver code.

Chapter28_Lesson4.m

% Chapter28_Lesson4.m
% Adaptive LPV model predictive control for a scalar constrained plant.
%
% Model:
%   x(k+1) = (theta(1) + theta(2)*rho(k))*x(k) + theta(3)*u(k) + w(k)
%
% The script combines:
%   1) projected recursive least squares,
%   2) measured/predicted LPV scheduling,
%   3) constrained finite-horizon optimization by projected gradient descent,
%   4) receding-horizon implementation.
%
% No toolbox is required. For production MPC, consider Model Predictive
% Control Toolbox, Optimization Toolbox (quadprog/fmincon), or YALMIP.

clear;
clc;
rng(7);

cfg.horizon = 12;
cfg.q = 5.0;
cfg.r = 0.15;
cfg.pTerminal = 8.0;
cfg.uMax = 1.5;
cfg.xMax = 3.0;
cfg.statePenalty = 80.0;
cfg.stepSize = 0.012;
cfg.iterations = 90;

thetaTrue = [0.72; 0.18; 0.55];

thetaHat = [0.55; 0.02; 0.35];
covariance = 60.0 * eye(3);
forgetting = 0.995;
lowerTheta = [0.35; -0.05; 0.15];
upperTheta = [1.05; 0.40; 0.90];

steps = 190;
x = 0.0;
warmStart = zeros(cfg.horizon, 1);

xLog = zeros(steps + 1, 1);
uLog = zeros(steps, 1);
rhoLog = zeros(steps, 1);
referenceLog = zeros(steps, 1);
thetaLog = zeros(steps, 3);
predictionErrorLog = zeros(steps, 1);

for k = 0:(steps - 1)
    rhoHorizon = zeros(cfg.horizon, 1);
    for i = 1:cfg.horizon
        rhoHorizon(i) = schedulingSignal(k + i - 1);
    end

    reference = referenceSignal(k);

    [optimalInput, ~, ~] = solveAdaptiveLPVMPC( ...
        x, ...
        rhoHorizon, ...
        reference, ...
        thetaHat, ...
        warmStart, ...
        cfg);

    u = optimalInput(1);
    rho = schedulingSignal(k);
    aTrue = thetaTrue(1) + thetaTrue(2) * rho;
    processNoise = 0.006 * randn();
    xNext = aTrue * x + thetaTrue(3) * u + processNoise;

    phi = [x; rho * x; u];
    [thetaHat, covariance, predictionError] = projectedRLSUpdate( ...
        thetaHat, ...
        covariance, ...
        phi, ...
        xNext, ...
        forgetting, ...
        lowerTheta, ...
        upperTheta);

    xLog(k + 1) = x;
    uLog(k + 1) = u;
    rhoLog(k + 1) = rho;
    referenceLog(k + 1) = reference;
    thetaLog(k + 1, :) = thetaHat.';
    predictionErrorLog(k + 1) = predictionError;

    warmStart(1:end-1) = optimalInput(2:end);
    warmStart(end) = optimalInput(end);
    x = xNext;
end

xLog(end) = x;

trackingRMSE = sqrt(mean((xLog(1:end-1) - referenceLog).^2));
maximumInput = max(abs(uLog));
maximumState = max(abs(xLog));
finalParameterError = norm(thetaHat - thetaTrue);
meanAbsolutePredictionError = mean(abs(predictionErrorLog));

fprintf('Adaptive LPV-MPC simulation completed.\n');
fprintf('Tracking RMSE: %.5f\n', trackingRMSE);
fprintf('Maximum |u|:  %.5f\n', maximumInput);
fprintf('Maximum |x|:  %.5f\n', maximumState);
fprintf('True theta:   [%.5f, %.5f, %.5f]\n', thetaTrue);
fprintf('Final theta:  [%.5f, %.5f, %.5f]\n', thetaHat);
fprintf('Final parameter error norm: %.5f\n', finalParameterError);
fprintf('Mean absolute one-step prediction error: %.5f\n', ...
    meanAbsolutePredictionError);

figure('Name', 'Adaptive LPV-MPC');
subplot(3, 1, 1);
plot(0:(steps - 1), xLog(1:end-1), 'LineWidth', 1.2);
hold on;
plot(0:(steps - 1), referenceLog, '--', 'LineWidth', 1.2);
grid on;
ylabel('State');
legend('x', 'reference', 'Location', 'best');
title('Adaptive LPV-MPC tracking');

subplot(3, 1, 2);
plot(0:(steps - 1), uLog, 'LineWidth', 1.2);
hold on;
yline(cfg.uMax, '--');
yline(-cfg.uMax, '--');
grid on;
ylabel('Input');

subplot(3, 1, 3);
plot(0:(steps - 1), thetaLog, 'LineWidth', 1.1);
grid on;
xlabel('Sample');
ylabel('Parameter estimates');
legend('\theta_0', '\theta_1', '\theta_2', 'Location', 'best');


function rho = schedulingSignal(k)
    rho = 0.75 * sin(0.045 * k) + 0.20 * sin(0.11 * k);
end


function reference = referenceSignal(k)
    if k < 25
        reference = 0.0;
    elseif k < 95
        reference = 1.0;
    elseif k < 145
        reference = -0.7;
    else
        reference = 0.4;
    end
end


function derivative = statePenaltyDerivative(x, cfg)
    excess = abs(x) - cfg.xMax;
    if excess <= 0.0
        derivative = 0.0;
        return;
    end

    derivative = 2.0 * cfg.statePenalty * excess * sign(x);
end


function states = predictStates(x0, input, rho, theta)
    horizon = numel(input);
    states = zeros(horizon + 1, 1);
    states(1) = x0;
    bHat = theta(3);

    for i = 1:horizon
        aHat = theta(1) + theta(2) * rho(i);
        states(i + 1) = aHat * states(i) + bHat * input(i);
    end
end


function [cost, gradient, states] = objectiveAndGradient( ...
    x0, input, rho, reference, theta, cfg)

    horizon = cfg.horizon;
    states = predictStates(x0, input, rho, theta);
    cost = 0.0;

    for i = 1:horizon
        stateError = states(i) - reference;
        excess = max(abs(states(i)) - cfg.xMax, 0.0);
        cost = cost ...
            + cfg.q * stateError^2 ...
            + cfg.r * input(i)^2 ...
            + cfg.statePenalty * excess^2;
    end

    terminalError = states(horizon + 1) - reference;
    terminalExcess = max( ...
        abs(states(horizon + 1)) - cfg.xMax, ...
        0.0);

    cost = cost ...
        + cfg.pTerminal * terminalError^2 ...
        + cfg.statePenalty * terminalExcess^2;

    gradient = zeros(horizon, 1);
    adjoint = 2.0 * cfg.pTerminal * terminalError ...
        + statePenaltyDerivative(states(horizon + 1), cfg);

    bHat = theta(3);
    for i = horizon:-1:1
        gradient(i) = 2.0 * cfg.r * input(i) + bHat * adjoint;

        aHat = theta(1) + theta(2) * rho(i);
        stageDerivative = 2.0 * cfg.q * (states(i) - reference) ...
            + statePenaltyDerivative(states(i), cfg);
        adjoint = stageDerivative + aHat * adjoint;
    end
end


function [bestInput, bestStates, bestCost] = solveAdaptiveLPVMPC( ...
    x0, rho, reference, theta, warmStart, cfg)

    input = min(max(warmStart, -cfg.uMax), cfg.uMax);
    bestInput = input;
    bestStates = predictStates(x0, input, rho, theta);
    bestCost = inf;

    for iteration = 0:(cfg.iterations - 1)
        [cost, gradient, states] = objectiveAndGradient( ...
            x0, ...
            input, ...
            rho, ...
            reference, ...
            theta, ...
            cfg);

        if cost < bestCost
            bestCost = cost;
            bestInput = input;
            bestStates = states;
        end

        step = cfg.stepSize / sqrt(1.0 + 0.03 * iteration);
        input = input - step * gradient;
        input = min(max(input, -cfg.uMax), cfg.uMax);
    end
end


function [thetaNew, covarianceNew, predictionError] = projectedRLSUpdate( ...
    theta, covariance, phi, y, forgetting, lowerTheta, upperTheta)

    denominator = forgetting + phi.' * covariance * phi;
    gain = covariance * phi / denominator;
    predictionError = y - phi.' * theta;

    thetaNew = theta + gain * predictionError;
    thetaNew = min(max(thetaNew, lowerTheta), upperTheta);

    covarianceNew = ( ...
        covariance - gain * phi.' * covariance ...
        ) / forgetting;
    covarianceNew = 0.5 * (covarianceNew + covarianceNew.');
end

The notebook uses core Wolfram Language operations. For exact quadratic programs, QuadraticOptimization is appropriate; for nonlinear adaptive prediction models, NMinimize can be used.

Chapter28_Lesson4.nb

(* Chapter28_Lesson4.nb
   Adaptive LPV model predictive control for a scalar constrained plant.

   The notebook combines projected recursive least squares, an LPV prediction
   model, projected-gradient MPC, and receding-horizon implementation.
*)

ClearAll["Global`*"];
SeedRandom[7];

cfg = <|
   "Horizon" -> 12,
   "Q" -> 5.0,
   "R" -> 0.15,
   "PTerminal" -> 8.0,
   "UMax" -> 1.5,
   "XMax" -> 3.0,
   "StatePenalty" -> 80.0,
   "StepSize" -> 0.012,
   "Iterations" -> 90
   |>;

schedulingSignal[k_Integer] :=
  0.75 Sin[0.045 k] + 0.20 Sin[0.11 k];

referenceSignal[k_Integer] := Piecewise[
  {
   {0.0, k < 25},
   {1.0, k < 95},
   {-0.7, k < 145}
   },
  0.4
  ];

statePenaltyDerivative[x_?NumericQ] := Module[
  {excess, sign},
  excess = Abs[x] - cfg["XMax"];
  If[excess <= 0.0,
   0.0,
   sign = If[x >= 0.0, 1.0, -1.0];
   2.0 cfg["StatePenalty"] excess sign
   ]
  ];

predictStates[
   x0_?NumericQ,
   input_List,
   rho_List,
   theta_List
   ] := Module[
  {horizon, states, bHat, i, aHat},
  horizon = Length[input];
  states = ConstantArray[0.0, horizon + 1];
  states[[1]] = x0;
  bHat = theta[[3]];

  For[i = 1, i <= horizon, i++,
   aHat = theta[[1]] + theta[[2]] rho[[i]];
   states[[i + 1]] =
    aHat states[[i]] + bHat input[[i]];
   ];

  states
  ];

objectiveAndGradient[
   x0_?NumericQ,
   input_List,
   rho_List,
   reference_?NumericQ,
   theta_List
   ] := Module[
  {
   horizon, states, cost, i, stateError, excess,
   terminalError, terminalExcess, gradient,
   adjoint, bHat, aHat, stageDerivative
   },

  horizon = cfg["Horizon"];
  states = predictStates[x0, input, rho, theta];
  cost = 0.0;

  For[i = 1, i <= horizon, i++,
   stateError = states[[i]] - reference;
   excess = Max[Abs[states[[i]]] - cfg["XMax"], 0.0];
   cost += cfg["Q"] stateError^2
     + cfg["R"] input[[i]]^2
     + cfg["StatePenalty"] excess^2;
   ];

  terminalError = states[[horizon + 1]] - reference;
  terminalExcess =
   Max[Abs[states[[horizon + 1]]] - cfg["XMax"], 0.0];

  cost += cfg["PTerminal"] terminalError^2
    + cfg["StatePenalty"] terminalExcess^2;

  gradient = ConstantArray[0.0, horizon];
  adjoint =
   2.0 cfg["PTerminal"] terminalError
    + statePenaltyDerivative[states[[horizon + 1]]];

  bHat = theta[[3]];
  For[i = horizon, i >= 1, i--,
   gradient[[i]] =
    2.0 cfg["R"] input[[i]] + bHat adjoint;

   aHat = theta[[1]] + theta[[2]] rho[[i]];
   stageDerivative =
    2.0 cfg["Q"] (states[[i]] - reference)
     + statePenaltyDerivative[states[[i]]];

   adjoint = stageDerivative + aHat adjoint;
   ];

  <|
   "Cost" -> cost,
   "Gradient" -> gradient,
   "States" -> states
   |>
  ];

solveAdaptiveLPVMPC[
   x0_?NumericQ,
   rho_List,
   reference_?NumericQ,
   theta_List,
   warmStart_List
   ] := Module[
  {
   input, bestInput, bestStates, bestCost,
   iteration, result, step
   },

  input = Clip[warmStart, {-cfg["UMax"], cfg["UMax"]}];
  bestInput = input;
  bestStates = predictStates[x0, input, rho, theta];
  bestCost = Infinity;

  For[iteration = 0, iteration < cfg["Iterations"], iteration++,
   result =
    objectiveAndGradient[x0, input, rho, reference, theta];

   If[result["Cost"] < bestCost,
    bestCost = result["Cost"];
    bestInput = input;
    bestStates = result["States"];
    ];

   step =
    cfg["StepSize"]/Sqrt[1.0 + 0.03 iteration];

   input = Clip[
     input - step result["Gradient"],
     {-cfg["UMax"], cfg["UMax"]}
     ];
   ];

  <|
   "Input" -> bestInput,
   "States" -> bestStates,
   "Cost" -> bestCost
   |>
  ];

projectedRLSUpdate[
   theta_List,
   covariance_List,
   phi_List,
   y_?NumericQ,
   forgetting_?NumericQ,
   lower_List,
   upper_List
   ] := Module[
  {
   denominator, gain, predictionError,
   thetaNew, covarianceNew
   },

  denominator =
   forgetting + phi . covariance . phi;
  gain = covariance . phi/denominator;
  predictionError = y - phi . theta;

  thetaNew =
   Clip[theta + gain predictionError, Transpose[{lower, upper}]];

  covarianceNew =
   (covariance - Outer[Times, gain, phi] . covariance)/
    forgetting;
  covarianceNew =
   0.5 (covarianceNew + Transpose[covarianceNew]);

  <|
   "Theta" -> thetaNew,
   "Covariance" -> covarianceNew,
   "PredictionError" -> predictionError
   |>
  ];

thetaTrue = {0.72, 0.18, 0.55};
thetaHat = {0.55, 0.02, 0.35};
covariance = 60.0 IdentityMatrix[3];
forgetting = 0.995;
lowerTheta = {0.35, -0.05, 0.15};
upperTheta = {1.05, 0.40, 0.90};

steps = 190;
x = 0.0;
warmStart = ConstantArray[0.0, cfg["Horizon"]];

xLog = ConstantArray[0.0, steps + 1];
uLog = ConstantArray[0.0, steps];
rhoLog = ConstantArray[0.0, steps];
referenceLog = ConstantArray[0.0, steps];
thetaLog = ConstantArray[0.0, {steps, 3}];
predictionErrorLog = ConstantArray[0.0, steps];

For[k = 0, k < steps, k++,
 rhoHorizon =
  Table[schedulingSignal[k + i], {i, 0, cfg["Horizon"] - 1}];
 reference = referenceSignal[k];

 mpcResult =
  solveAdaptiveLPVMPC[
   x,
   rhoHorizon,
   reference,
   thetaHat,
   warmStart
   ];

 u = mpcResult["Input"][[1]];
 rho = schedulingSignal[k];
 aTrue = thetaTrue[[1]] + thetaTrue[[2]] rho;
 processNoise =
  0.006 RandomVariate[NormalDistribution[0.0, 1.0]];
 xNext = aTrue x + thetaTrue[[3]] u + processNoise;

 phi = {x, rho x, u};
 rlsResult =
  projectedRLSUpdate[
   thetaHat,
   covariance,
   phi,
   xNext,
   forgetting,
   lowerTheta,
   upperTheta
   ];

 thetaHat = rlsResult["Theta"];
 covariance = rlsResult["Covariance"];

 xLog[[k + 1]] = x;
 uLog[[k + 1]] = u;
 rhoLog[[k + 1]] = rho;
 referenceLog[[k + 1]] = reference;
 thetaLog[[k + 1]] = thetaHat;
 predictionErrorLog[[k + 1]] =
  rlsResult["PredictionError"];

 warmStart =
  Join[
   Rest[mpcResult["Input"]],
   {Last[mpcResult["Input"]]}
   ];

 x = xNext;
 ];

xLog[[-1]] = x;

trackingRMSE =
 Sqrt[Mean[(Most[xLog] - referenceLog)^2]];
maximumInput = Max[Abs[uLog]];
maximumState = Max[Abs[xLog]];
finalParameterError = Norm[thetaHat - thetaTrue];
meanAbsolutePredictionError =
 Mean[Abs[predictionErrorLog]];

Print["Adaptive LPV-MPC simulation completed."];
Print["Tracking RMSE: ", NumberForm[trackingRMSE, {8, 5}]];
Print["Maximum |u|:  ", NumberForm[maximumInput, {8, 5}]];
Print["Maximum |x|:  ", NumberForm[maximumState, {8, 5}]];
Print["True theta:   ", thetaTrue];
Print["Final theta:  ", thetaHat];
Print[
 "Final parameter error norm: ",
 NumberForm[finalParameterError, {8, 5}]
 ];
Print[
 "Mean absolute one-step prediction error: ",
 NumberForm[meanAbsolutePredictionError, {8, 5}]
 ];

trackingPlot = ListLinePlot[
  {
   Transpose[{Range[0, steps - 1], Most[xLog]}],
   Transpose[{Range[0, steps - 1], referenceLog}]
   },
  PlotLegends -> {"x", "reference"},
  Frame -> True,
  FrameLabel -> {"Sample", "State"},
  PlotLabel -> "Adaptive LPV-MPC Tracking",
  ImageSize -> Large
  ];

parameterPlot = ListLinePlot[
  Transpose[thetaLog],
  PlotLegends -> {"theta0", "theta1", "theta2"},
  Frame -> True,
  FrameLabel -> {"Sample", "Estimate"},
  PlotLabel -> "Online Parameter Estimates",
  ImageSize -> Large
  ];

Column[{trackingPlot, parameterPlot}]

15.1 Numerical safeguards

• scale states, inputs, regressors, and cost weights;
• symmetrize and regularize covariance matrices;
• project parameter estimates onto physically valid bounds;
• warm-start the optimizer by shifting the prior solution;
• impose iteration and wall-clock limits;
• retain a verified fallback controller when optimization fails;
• log feasibility, solver residuals, active constraints, and estimator covariance.

15.2 Sampling and time-scale separation

Adaptation, optimization, and plant dynamics should not all be tuned as if they were independent. Excessively rapid parameter updates can make the prediction model vary faster than the optimizer and terminal design can tolerate. A practical implementation may filter estimates, update the model at a slower supervisory rate, or limit \(\|\hat{\theta}_{k+1}-\hat{\theta}_k\|\).

15.3 Validation hierarchy

1. offline Monte Carlo simulation over the initial uncertainty set;
2. constraint stress tests and infeasibility injection;
3. processor-in-the-loop timing tests;
4. hardware-in-the-loop testing with sensor and actuator faults;
5. limited-envelope experiments with an independent safety monitor.

flowchart TD
  S["Start with constrained varying dynamics"] --> M["Is a reliable scheduling \nvariable measured?"]
  M -->|yes| L["Use an LPV prediction model"]
  M -->|no| U["Use an uncertain LTI \nor nonlinear model"]
  L --> K["Are LPV coefficients \naccurately known?"]
  K -->|yes| LPV["LPV MPC"]
  K -->|no| ALPV["Adaptive LPV MPC"]
  U --> C["Is point-estimate \nconstraint risk acceptable?"]
  C -->|yes| CE["Certainty-equivalent \nAdaptive MPC"]
  C -->|no| RAMPC["Robust Adaptive \nMPC with uncertainty \nsets or tubes"]
  ALPV --> G["Need guaranteed constraints?"]
  G -->|yes| RLPV["Robust Adaptive LPV MPC"]
  G -->|no| CELPV["Certainty-equivalent \nAdaptive LPV MPC"]
        

The strongest architecture is not always the best engineering choice. Robust adaptive LPV-MPC may be too expensive for a fast embedded loop, while a simpler scheduled controller with supervisory adaptation and independent constraint protection may be more verifiable.

Problem 1: LPV Condensed Prediction

For \(x_{i+1}=A_ix_i+B_iu_i\), derive \(x_1,x_2,x_3\) as affine functions of \(x_0,u_0,u_1,u_2\).

Solution:

\[ x_1=A_0x_0+B_0u_0. \]

\[ \begin{aligned} x_2 &=A_1x_1+B_1u_1\\ &=A_1A_0x_0+A_1B_0u_0+B_1u_1. \end{aligned} \]

\[ \begin{aligned} x_3 &=A_2x_2+B_2u_2\\ &=A_2A_1A_0x_0 +A_2A_1B_0u_0 +A_2B_1u_1 +B_2u_2. \end{aligned} \]

These block products form the matrices \(\mathcal{A}_k\) and \(\mathcal{B}_k\). Unlike LTI MPC, the products depend on the predicted scheduling sequence.

Problem 2: Regression for the Adaptive LPV Model

Show that \(x_{k+1}=(\theta_0+\theta_1\rho_k)x_k+\theta_2u_k+w_k\) is linear in the unknown parameter vector.

Solution:

\[ x_{k+1} = \theta_0x_k+ \theta_1(\rho_kx_k)+ \theta_2u_k+w_k. \]

Define

\[ \theta= \begin{bmatrix} \theta_0 & \theta_1 & \theta_2 \end{bmatrix}^{\mathsf{T}}, \qquad \phi_k= \begin{bmatrix} x_k & \rho_kx_k & u_k \end{bmatrix}^{\mathsf{T}}. \]

Then

\[ x_{k+1}=\phi_k^\mathsf{T}\theta+w_k, \]

which is suitable for RLS even though the physical dynamics vary with \(\rho_k\).

Problem 3: Numerical RLS Update

Let \(\hat{\theta}_k=[0.6,\;0.1,\;0.4]^\mathsf{T}\), \(P_k=10I\), \(\lambda=1\), \(x_k=1\), \(\rho_k=0.5\), \(u_k=0.2\), and \(x_{k+1}=0.82\). Compute the unprojected RLS update.

Solution:

\[ \phi_k= \begin{bmatrix} 1 & 0.5 & 0.2 \end{bmatrix}^{\mathsf{T}}. \]

\[ \hat{x}_{k+1|k} = \phi_k^\mathsf{T}\hat{\theta}_k = 0.6+0.05+0.08=0.73. \]

\[ \tilde{y}_k=0.82-0.73=0.09. \]

\[ \phi_k^\mathsf{T}P_k\phi_k = 10(1^2+0.5^2+0.2^2)=12.9. \]

\[ K_k= \frac{10\phi_k}{1+12.9} = \frac{1}{13.9} \begin{bmatrix} 10 & 5 & 2 \end{bmatrix}^{\mathsf{T}}. \]

\[ \hat{\theta}_{k+1} = \hat{\theta}_k+K_k(0.09) \approx \begin{bmatrix} 0.66475 & 0.13237 & 0.41295 \end{bmatrix}^{\mathsf{T}}. \]

Problem 4: Nested Sets and Robust Feasibility

Explain why \(\Theta_{k+1}\subseteq\Theta_k\) is useful but not by itself sufficient for recursive feasibility.

Solution:

The inclusion means that uncertainty has not expanded, so any constraint enforced for all models in \(\Theta_k\) also covers all models in \(\Theta_{k+1}\). However, recursive feasibility also requires the realized disturbed successor state to remain inside the robust feasible region. This normally needs a robustly invariant tube or terminal set, a bounded disturbance model, and a valid shifted candidate sequence. Nested parameter sets remove one source of expansion but do not replace the invariant-set argument.

Problem 5: Stability Under Model Updates

Suppose

\[ V_{k+1}-V_k \leq -c\|x_k\|^2+ L\|\Delta\hat{\theta}_k\|. \]

State one condition under which asymptotic convergence remains plausible and one condition that yields only practical stability.

Solution:

If \(\sum_{k=0}^{\infty}\|\Delta\hat{\theta}_k\|<\infty\), the cumulative perturbation of the Lyapunov decrease is finite. With standard boundedness and detectability conditions, this can support an asymptotic convergence argument. If instead \(\|\Delta\hat{\theta}_k\|\leq\bar{d}_\theta\) with a nonzero constant bound, then \(V_k\) is guaranteed to decrease only outside a neighborhood whose size depends on \(L\bar{d}_\theta/c\). The appropriate conclusion is ultimate boundedness or practical stability.

Problem 6: Frozen Versus Anticipative Scheduling

A plant has \(A(\rho)=A_0+\rho A_1\). Compare frozen scheduling with an exact known future scheduling trajectory.

Solution:

Frozen scheduling sets every prediction matrix to \(A(\rho_k)\). It is simple but introduces prediction mismatch when \(\rho\) changes appreciably over the horizon. Anticipative scheduling uses \(A(\rho_{i|k})\) at each prediction stage. When the trajectory is known accurately, this removes that source of mismatch and can reduce conservatism. When the schedule forecast is uncertain, a set-valued or robust treatment is needed; treating a poor forecast as exact can be less safe than freezing the current value.

Problem 7: Why Parameter Convergence Is Not Guaranteed by MPC Tracking

Explain why small tracking error does not imply \(\hat{\theta}_k→\theta^\star\).

Solution:

MPC may regulate the state near a constant reference using an input-state trajectory that does not persistently excite every regressor direction. Several parameter vectors can then predict nearly identical local behavior. Tracking error can converge while the information matrix remains rank deficient. Parameter convergence still requires an excitation or identifiability condition, and deliberate excitation must respect the MPC constraints.

Adaptive MPC places an online-updated model inside a receding-horizon optimization. LPV MPC uses measured or predicted scheduling variables to describe known operating-condition dependence. The two ideas can be combined by estimating uncertain LPV coefficients online.

Certainty-equivalent implementations are relatively simple but do not automatically protect constraints against estimation error. Robust Adaptive MPC augments the point estimate with parameter sets, tubes, tightened constraints, or min-max predictions. Nested uncertainty sets, invariant terminal ingredients, and bounded model updates are central to recursive-feasibility and stability arguments.

Full derivations of terminal sets, robust invariant tubes, multiparametric quadratic programming, nonlinear programming, and real-time solver design belong in dedicated MPC and Robust Control courses. The adaptive-control contribution is the estimator-controller coupling, the information effect of control inputs, and the need to reconcile learning with stability and constraints.