All three methods replace a single fixed controller by a family of controllers, but they obtain the active controller in fundamentally different ways. The distinction is not merely terminology: it determines what must be measured, what must be proven, and what can fail during rapid plant variation.

flowchart TD
  P["Changing plant or operating condition"] --> Q["What information is available?"]
  Q -->|"Tracking error and regressors"| M["MRAC"]
  Q -->|"Measured scheduling variable"| G["Gain scheduling"]
  Q -->|"Input-output data for online model"| S["Self-tuning regulator"]
  M --> M1["Update controller \nparameters directly"]
  G --> G1["Select or interpolate \noffline gains"]
  S --> S1["Estimate plant, \nthen redesign controller"]
        

In compact notation, let \(K\) denote all controller parameters. The three update mechanisms can be summarized as

\[ \begin{aligned} \text{MRAC:}\quad & \dot{K}=\mathcal{A}(e,\phi), \\ \text{gain scheduling:}\quad & K=\mathcal{S}(\rho), \\ \text{STR:}\quad & \hat{\theta}_{p,k+1}=\mathcal{I}(\hat{\theta}_{p,k},u_k,y_k), \qquad K_k=\mathcal{D}(\hat{\theta}_{p,k}). \end{aligned} \]

Here \(e\) is a tracking error, \(\phi\) is a known regressor, \(\rho\) is a scheduling variable, \(\hat{\theta}_p\) is an estimated plant-parameter vector, \(\mathcal{A}\) is an adaptive law, \(\mathcal{S}\) is a precomputed schedule, \(\mathcal{I}\) is an online identifier, and \(\mathcal{D}\) is a controller-design map.

To compare architectures rather than unrelated plant examples, consider the first-order plant

\[ \dot{y}(t)=-a(t)y(t)+b u(t), \qquad b>0, \]

and the desired reference model

\[ \dot{y}_m(t)=-a_m y_m(t)+b_m r(t), \qquad a_m>0. \]

The comparison objective is model following:

\[ e(t)=y(t)-y_m(t), \qquad e(t)→0. \]

A controller family with two adjustable coefficients is

\[ u(t)=\theta_r(t)r(t)-\theta_y(t)y(t). \]

If \(a\) and \(b\) were known and constant, exact model matching would be obtained from the ideal parameters

\[ \theta_y^*=\frac{a_m-a}{b}, \qquad \theta_r^*=\frac{b_m}{b}. \]

Each architecture can now be interpreted as a different mechanism for producing, approximating, or learning these ideal coefficients.

Define parameter errors \(\tilde{\theta}_y=\theta_y-\theta_y^*\) and \(\tilde{\theta}_r=\theta_r-\theta_r^*\). Substitution of the controller into the plant gives

\[ \dot{e}=-a_m e-b\tilde{\theta}_y y+b\tilde{\theta}_r r. \]

For known positive control direction \(b>0\), choose the Lyapunov function

\[ V=\frac{1}{2}e^2+ rac{b}{2\gamma_y}\tilde{\theta}_y^2 +\frac{b}{2\gamma_r}\tilde{\theta}_r^2, \qquad \gamma_y>0,\;\gamma_r>0. \]

Differentiating along the error dynamics yields

\[ \begin{aligned} \dot{V} &=-a_m e^2-b e\tilde{\theta}_y y+b e\tilde{\theta}_r r +\frac{b}{\gamma_y}\tilde{\theta}_y\dot{\theta}_y +\frac{b}{\gamma_r}\tilde{\theta}_r\dot{\theta}_r \\ &=-a_m e^2 +b\tilde{\theta}_y\left(-ey+\frac{\dot{\theta}_y}{\gamma_y}\right) +b\tilde{\theta}_r\left(er+\frac{\dot{\theta}_r}{\gamma_r}\right). \end{aligned} \]

The cross terms cancel under the direct adaptive laws

\[ \boxed{\dot{\theta}_y=\gamma_y e y}, \qquad \boxed{\dot{\theta}_r=-\gamma_r e r}, \]

and consequently

\[ \dot{V}=-a_m e^2\leq 0. \]

Thus \(V(t)\), \(e(t)\), and the parameter errors remain bounded, while \(e\in L_2\). Under bounded reference and closed-loop signals, \(\dot e\) is bounded; the Barbalat argument developed in Chapter 3 then gives \(e(t)→0\). This tracking result does not by itself imply that every parameter estimate reaches its unique physical value; the extra signal-richness conditions for parameter convergence are deferred to Chapter 10.

Structural interpretation: MRAC does not first identify \(a\) and \(b\). It moves the controller coefficients in directions chosen to decrease the reference-model error while preserving a Lyapunov inequality.

Suppose a measurable operating variable \(\rho(t)\) is strongly related to the plant parameter, ideally \(\rho(t)=a(t)\). Controllers are designed offline at grid points \(\rho_1,\ldots,\rho_N\). For the scalar plant,

\[ \theta_y(\rho_i)=\frac{a_m-\rho_i}{b}, \qquad \theta_r(\rho_i)=\frac{b_m}{b}. \]

Between grid points, a convex interpolation is commonly used. If \(\rho_i\leq\rho\leq\rho_{i+1}\), define

\[ \alpha(\rho)=\frac{\rho-\rho_i}{\rho_{i+1}-\rho_i}, \qquad K(\rho)=\left(1-\alpha(\rho)\right)K_i+\alpha(\rho)K_{i+1}. \]

No tracking-error adaptation law is integrated. The controller can therefore change immediately when \(\rho\) changes, but only according to the relation anticipated during offline design.

4.1 Frozen-Parameter Stability Is Not Automatically Transition Stability

Let a gain-scheduled state model be \(\dot{x}=A_c(\rho)x\). Assume that every frozen matrix \(A_c(\rho)\) is stable and that a differentiable parameter-dependent matrix \(P(\rho)=P^T(\rho)>0\) satisfies

\[ A_c^T(\rho)P(\rho)+P(\rho)A_c(\rho)\leq-Q, \qquad Q=Q^T>0. \]

For \(V=x^TP(\rho)x\),

\[ \dot{V} =x^T\left(A_c^TP+PA_c\right)x +x^T\frac{\partial P}{\partial\rho}x\,\dot{\rho}. \]

Hence the bound

\[ \dot{V}\leq- \left(\lambda_{\min}(Q)- \left\|\frac{\partial P}{\partial\rho}\right\|\,|\dot{\rho}|\right) \|x\|^2. \]

A sufficient slow-variation condition is therefore

\[ |\dot{\rho}|< \frac{\lambda_{\min}(Q)}{ \sup_{\rho}\left\|\partial P(\rho)/\partial\rho\right\|}. \]

This inequality explains why stable local controllers may still behave poorly during aggressive transitions. A common quadratic Lyapunov matrix independent of \(\rho\) removes the derivative term and gives a stronger guarantee, but finding such a common certificate may be conservative or impossible for a chosen schedule.

A self-tuning regulator is an indirect adaptive architecture. For a sampled first-order plant, use the regression model

\[ y_{k+1}=\alpha y_k+\beta u_k =\phi_k^T\theta_p, \qquad \phi_k=\begin{bmatrix}y_k & u_k\end{bmatrix}^T, \quad \theta_p=\begin{bmatrix}\alpha & \beta\end{bmatrix}^T. \]

An online recursive least-squares update with forgetting factor \(0<\lambda\leq1\) is

\[ \begin{aligned} L_k&=\frac{P_k\phi_k}{\lambda+\phi_k^TP_k\phi_k},\\ \hat{\theta}_{p,k+1} &=\hat{\theta}_{p,k}+L_k \left(y_{k+1}-\phi_k^T\hat{\theta}_{p,k}\right),\\ P_{k+1}&=\frac{1}{\lambda} \left(P_k-L_k\phi_k^TP_k\right). \end{aligned} \]

Detailed estimator properties are developed in Chapters 11 and 12. Here the equations are used only to expose the architectural distinction.

5.1 Certainty-Equivalent Pole Placement

Choose a discrete reference model

\[ y_{m,k+1}=\lambda_m y_{m,k}+(1-\lambda_m)r_k, \qquad 0<\lambda_m<1. \]

For the controller \(u_k=k_y y_k+k_r r_k\),

\[ y_{k+1}=(\alpha+\beta k_y)y_k+\beta k_r r_k. \]

Exact matching requires

\[ k_y^*=\frac{\lambda_m-\alpha}{\beta}, \qquad k_r^*=\frac{1-\lambda_m}{\beta}. \]

The self-tuning implementation substitutes estimates:

\[ k_{y,k}=\frac{\lambda_m-\hat{\alpha}_k}{\hat{\beta}_k}, \qquad k_{r,k}=\frac{1-\lambda_m}{\hat{\beta}_k}. \]

This is called certainty equivalence: the current estimate is treated as if it were the true plant. The denominator reveals an immediate practical issue. If \(\hat\beta_k\) approaches zero, controller gains become excessive. Real implementations therefore impose sign knowledge, lower bounds, projection, saturation, or redesign logic.

The three architectures can all be written as time-varying controller parameters \(K(t)\), but the state that generates \(K\) differs:

\[ \begin{array}{lll} \text{MRAC:} & \dot{\eta}=f_{\mathrm{A} }(\eta,e,\phi), & K=h_{\mathrm{A} }(\eta),\\ \text{gain scheduling:} & \eta=\rho, & K=h_{\mathrm{S} }(\rho),\\ \text{STR:} & \eta=\hat\theta_p,\quad \eta_{k+1}=f_{\mathrm{I} }(\eta_k,u_k,y_k), & K_k=h_{\mathrm{D} }(\eta_k). \end{array} \]

This representation highlights three different sources of model knowledge:

• MRAC needs a controller parameterization and matching structure that permit stable error-based adaptation.
• Gain scheduling needs an offline map from operating condition to suitable controller gains and a reliable scheduling measurement.
• STR needs an identifiable online plant model and a controller-design formula that remains well conditioned as estimates evolve.

Therefore, the phrase “adaptive controller” is broad. Gain scheduling changes gains online but usually does not learn from performance error. MRAC and STR both learn online, but MRAC learns controller parameters directly whereas an indirect STR learns plant parameters first.

flowchart TD
  A["Start with a changing plant"] --> B["Is a reliable \noperating variable measured?"]
  B -->|"yes"| C["Can controllers \nbe designed across \nthe full envelope?"]
  C -->|"yes"| G["Gain scheduling is \na strong baseline"]
  C -->|"no"| D["Is a stable model-following \nparameterization available?"]
  B -->|"no"| D
  D -->|"yes"| M["Consider MRAC"]
  D -->|"no"| E["Can an online \nparametric plant model \nbe estimated?"]
  E -->|"yes"| S["Consider an STR"]
  E -->|"no"| F["Use a fixed robust \ndesign or improve \nsensing/model structure"]
  G --> H["Validate interpolation and \ntransition stability"]
  M --> I["Validate matching, \nsignal bounds, \nand adaptation transients"]
  S --> J["Validate estimator conditioning \nand redesign safety"]
        

8.1 Gain Scheduling Has an Advantage When the Change Is Measured

If airspeed, payload configuration, shaft speed, or another operating variable is accurately measured and the plant family has been mapped offline, gain scheduling can change gains before a large tracking error develops. Its weakness is dependence on the validity and coverage of that map.

8.2 MRAC Has an Advantage When Desired Dynamics Are Clear but Parameters Are Not

MRAC is natural when a reference model is a direct performance specification and uncertain terms enter through a known parameterization. It avoids explicit plant identification, but adaptation needs an error transient and can interact strongly with noise, saturation, delay, and omitted dynamics.

8.3 STR Has an Advantage When an Estimated Model Is Valuable Beyond One Gain

An STR is attractive when an online model can support pole placement, minimum variance control, or another redesign calculation. The explicit estimate is interpretable and reusable, but its reliability depends on data quality and numerical conditioning. A good prediction model is not automatically a safe closed-loop controller unless estimator and redesign safeguards are included.

The downloadable implementations use the same plant and command for all three controllers. The unknown pole magnitude changes as

\[ a(t)= \begin{cases} 0.8, & 0\leq t<12,\\ 1.6, & 12\leq t<24,\\ 0.55, & 24\leq t\leq36, \end{cases} \qquad b=1. \]

The reference model is

\[ \dot y_m=-1.2y_m+1.2r. \]

The gain-scheduling benchmark is intentionally favourable: it is given the exact value \(\rho=a(t)\) and a table containing the correct matching gains. It therefore establishes what can be achieved when the scheduling variable perfectly exposes the plant change. MRAC receives no \(a(t)\) measurement and adapts from tracking error. The STR receives input-output samples and estimates a discrete model by RLS.

These numbers are not a universal ranking. Changing measurement noise, scheduling bias, model mismatch, rate of variation, initial estimates, or actuator limits can reverse the outcome. The experiment is designed to make architectural mechanisms visible, not to declare one method globally superior.

The Python program uses NumPy for numerical arrays and Matplotlib for plots. It also writes all trajectories to CSV.

Chapter4_Lesson5.py

"""Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators.

Simulates three controllers on the same first-order plant
    y_dot = -a(t) y + b u
with a piecewise-changing parameter a(t).
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass
class Result:
    t: np.ndarray
    r: np.ndarray
    ym: np.ndarray
    y_mrac: np.ndarray
    y_gs: np.ndarray
    y_str: np.ndarray
    u_mrac: np.ndarray
    u_gs: np.ndarray
    u_str: np.ndarray
    theta_y: np.ndarray
    theta_r: np.ndarray
    alpha_hat: np.ndarray
    beta_hat: np.ndarray


def plant_parameter(t: float) -> float:
    """Unknown/variable plant pole magnitude used in the experiment."""
    if t < 12.0:
        return 0.8
    if t < 24.0:
        return 1.6
    return 0.55


def reference_signal(t: float) -> float:
    if t < 2.0:
        return 0.0
    if t < 18.0:
        return 1.0
    if t < 28.0:
        return -0.5
    return 0.75


def interpolate_gain(rho: float, rho_grid: np.ndarray, gain_grid: np.ndarray) -> float:
    return float(np.interp(rho, rho_grid, gain_grid))


def simulate(dt: float = 0.01, t_final: float = 36.0) -> Result:
    n = int(round(t_final / dt)) + 1
    t = np.linspace(0.0, t_final, n)
    r = np.array([reference_signal(ti) for ti in t])

    # Reference model: y_m_dot = -a_m y_m + b_m r
    a_m, b_m = 1.2, 1.2
    b_true = 1.0

    ym = np.zeros(n)
    y_mrac = np.zeros(n)
    y_gs = np.zeros(n)
    y_str = np.zeros(n)
    u_mrac = np.zeros(n)
    u_gs = np.zeros(n)
    u_str = np.zeros(n)

    # Direct MRAC states. The sign of b is assumed known and positive.
    theta_y = np.zeros(n)
    theta_r = np.zeros(n)
    theta_y[0] = 0.0
    theta_r[0] = 0.5
    gamma_y, gamma_r = 3.0, 2.0

    # Gain-scheduled table designed offline for measured scheduling variable rho=a(t).
    rho_grid = np.array([0.50, 0.80, 1.20, 1.60, 1.80])
    theta_y_grid = (a_m - rho_grid) / b_true
    theta_r_grid = np.full_like(rho_grid, b_m / b_true)

    # STR estimator for y[k+1] = alpha*y[k] + beta*u[k].
    alpha_hat = np.zeros(n)
    beta_hat = np.zeros(n)
    alpha_hat[0] = np.exp(-1.0 * dt)
    beta_hat[0] = 0.8 * dt
    p = np.eye(2) * 200.0
    forgetting = 0.998
    lambda_m = np.exp(-a_m * dt)

    for k in range(n - 1):
        tk = t[k]
        a_true = plant_parameter(tk)

        # Shared reference model.
        ym[k + 1] = ym[k] + dt * (-a_m * ym[k] + b_m * r[k])

        # 1) MRAC: control parameters move directly from tracking error.
        e = y_mrac[k] - ym[k]
        u_mrac[k] = theta_r[k] * r[k] - theta_y[k] * y_mrac[k]
        theta_y_dot = gamma_y * e * y_mrac[k]
        theta_r_dot = -gamma_r * e * r[k]
        theta_y[k + 1] = np.clip(theta_y[k] + dt * theta_y_dot, -5.0, 5.0)
        theta_r[k + 1] = np.clip(theta_r[k] + dt * theta_r_dot, -5.0, 5.0)
        y_mrac[k + 1] = y_mrac[k] + dt * (-a_true * y_mrac[k] + b_true * u_mrac[k])

        # 2) Gain scheduling: use measured rho and interpolate precomputed gains.
        rho_measured = a_true  # ideal scheduling sensor for this comparison
        ky_gs = interpolate_gain(rho_measured, rho_grid, theta_y_grid)
        kr_gs = interpolate_gain(rho_measured, rho_grid, theta_r_grid)
        u_gs[k] = kr_gs * r[k] - ky_gs * y_gs[k]
        y_gs[k + 1] = y_gs[k] + dt * (-a_true * y_gs[k] + b_true * u_gs[k])

        # 3) STR: identify a discrete plant model, then redesign controller gains.
        beta_safe = beta_hat[k]
        if abs(beta_safe) < 0.15 * dt:
            beta_safe = np.copysign(0.15 * dt, beta_safe if beta_safe != 0.0 else 1.0)
        ky_str = (lambda_m - alpha_hat[k]) / beta_safe
        kr_str = (1.0 - lambda_m) / beta_safe
        ky_str = float(np.clip(ky_str, -12.0, 12.0))
        kr_str = float(np.clip(kr_str, -12.0, 12.0))
        u_str[k] = ky_str * y_str[k] + kr_str * r[k]
        y_str[k + 1] = y_str[k] + dt * (-a_true * y_str[k] + b_true * u_str[k])

        # RLS update uses the newly observed y[k+1].
        phi = np.array([y_str[k], u_str[k]])
        denominator = forgetting + phi @ p @ phi
        gain = (p @ phi) / denominator
        prediction_error = y_str[k + 1] - phi @ np.array([alpha_hat[k], beta_hat[k]])
        estimate = np.array([alpha_hat[k], beta_hat[k]]) + gain * prediction_error
        p = (p - np.outer(gain, phi) @ p) / forgetting
        alpha_hat[k + 1] = float(np.clip(estimate[0], 0.85, 1.02))
        beta_hat[k + 1] = float(np.clip(estimate[1], 0.15 * dt, 2.5 * dt))

    u_mrac[-1] = u_mrac[-2]
    u_gs[-1] = u_gs[-2]
    u_str[-1] = u_str[-2]

    return Result(
        t=t,
        r=r,
        ym=ym,
        y_mrac=y_mrac,
        y_gs=y_gs,
        y_str=y_str,
        u_mrac=u_mrac,
        u_gs=u_gs,
        u_str=u_str,
        theta_y=theta_y,
        theta_r=theta_r,
        alpha_hat=alpha_hat,
        beta_hat=beta_hat,
    )


def rms(x: np.ndarray) -> float:
    return float(np.sqrt(np.mean(np.square(x))))


def save_csv(result: Result, path: Path) -> None:
    with path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow(
            [
                "t",
                "r",
                "ym",
                "y_mrac",
                "y_gain_scheduling",
                "y_str",
                "u_mrac",
                "u_gain_scheduling",
                "u_str",
                "theta_y_mrac",
                "theta_r_mrac",
                "alpha_hat_str",
                "beta_hat_str",
            ]
        )
        for row in zip(
            result.t,
            result.r,
            result.ym,
            result.y_mrac,
            result.y_gs,
            result.y_str,
            result.u_mrac,
            result.u_gs,
            result.u_str,
            result.theta_y,
            result.theta_r,
            result.alpha_hat,
            result.beta_hat,
        ):
            writer.writerow(row)


def main() -> None:
    result = simulate()
    print(f"MRAC RMS model-following error: {rms(result.y_mrac - result.ym):.6f}")
    print(f"Gain-scheduling RMS model-following error: {rms(result.y_gs - result.ym):.6f}")
    print(f"STR RMS model-following error: {rms(result.y_str - result.ym):.6f}")

    output_csv = Path(__file__).with_name("Chapter4_Lesson5_results.csv")
    save_csv(result, output_csv)
    print(f"Saved {output_csv}")

    plt.figure(figsize=(10, 6))
    plt.plot(result.t, result.r, "k:", label="command r")
    plt.plot(result.t, result.ym, "k--", label="reference model")
    plt.plot(result.t, result.y_mrac, label="MRAC")
    plt.plot(result.t, result.y_gs, label="gain scheduling")
    plt.plot(result.t, result.y_str, label="STR")
    plt.xlabel("Time [s]")
    plt.ylabel("Output")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(10, 5))
    plt.plot(result.t, result.u_mrac, label="u MRAC")
    plt.plot(result.t, result.u_gs, label="u gain scheduling")
    plt.plot(result.t, result.u_str, label="u STR")
    plt.xlabel("Time [s]")
    plt.ylabel("Control input")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library and writes a CSV file suitable for plotting in MATLAB, Python, a spreadsheet, or another tool.

Chapter4_Lesson5.cpp

// Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators.
// Build: g++ -std=c++17 -O2 Chapter4_Lesson5.cpp -o Chapter4_Lesson5

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

struct Sample {
    double t, r, ym, yMrac, yGs, yStr;
    double uMrac, uGs, uStr;
    double thetaY, thetaR, alphaHat, betaHat;
};

double plantParameter(double t) {
    if (t < 12.0) return 0.8;
    if (t < 24.0) return 1.6;
    return 0.55;
}

double referenceSignal(double t) {
    if (t < 2.0) return 0.0;
    if (t < 18.0) return 1.0;
    if (t < 28.0) return -0.5;
    return 0.75;
}

double clampValue(double x, double lo, double hi) {
    return std::max(lo, std::min(x, hi));
}

double linearInterpolate(double x, const std::vector<double>& grid,
                         const std::vector<double>& values) {
    if (x <= grid.front()) return values.front();
    if (x >= grid.back()) return values.back();
    for (std::size_t i = 0; i + 1 < grid.size(); ++i) {
        if (x <= grid[i + 1]) {
            const double w = (x - grid[i]) / (grid[i + 1] - grid[i]);
            return (1.0 - w) * values[i] + w * values[i + 1];
        }
    }
    return values.back();
}

double rms(const std::vector<double>& errors) {
    double sum = 0.0;
    for (double e : errors) sum += e * e;
    return std::sqrt(sum / static_cast<double>(errors.size()));
}

int main() {
    constexpr double dt = 0.01;
    constexpr double tFinal = 36.0;
    constexpr double aModel = 1.2;
    constexpr double bModel = 1.2;
    constexpr double bTrue = 1.0;
    constexpr double gammaY = 3.0;
    constexpr double gammaR = 2.0;
    constexpr double forgetting = 0.998;

    const int steps = static_cast<int>(std::round(tFinal / dt)) + 1;
    const double lambdaModel = std::exp(-aModel * dt);

    const std::vector<double> rhoGrid{0.50, 0.80, 1.20, 1.60, 1.80};
    std::vector<double> thetaYGrid;
    std::vector<double> thetaRGrid;
    for (double rho : rhoGrid) {
        thetaYGrid.push_back((aModel - rho) / bTrue);
        thetaRGrid.push_back(bModel / bTrue);
    }

    double ym = 0.0, yMrac = 0.0, yGs = 0.0, yStr = 0.0;
    double thetaY = 0.0, thetaR = 0.5;
    double alphaHat = std::exp(-1.0 * dt), betaHat = 0.8 * dt;
    std::array<std::array<double, 2>, 2> P{ { {200.0, 0.0}, {0.0, 200.0} } };

    std::vector<Sample> data;
    data.reserve(steps);
    std::vector<double> eMrac, eGs, eStr;
    eMrac.reserve(steps); eGs.reserve(steps); eStr.reserve(steps);

    for (int k = 0; k < steps; ++k) {
        const double t = k * dt;
        const double r = referenceSignal(t);
        const double aTrue = plantParameter(t);

        const double e = yMrac - ym;
        const double uMrac = thetaR * r - thetaY * yMrac;

        const double thetaYGs = linearInterpolate(aTrue, rhoGrid, thetaYGrid);
        const double thetaRGs = linearInterpolate(aTrue, rhoGrid, thetaRGrid);
        const double uGs = thetaRGs * r - thetaYGs * yGs;

        double betaSafe = betaHat;
        if (std::abs(betaSafe) < 0.15 * dt) betaSafe = 0.15 * dt;
        const double kyStr = clampValue((lambdaModel - alphaHat) / betaSafe, -12.0, 12.0);
        const double krStr = clampValue((1.0 - lambdaModel) / betaSafe, -12.0, 12.0);
        const double uStr = kyStr * yStr + krStr * r;

        data.push_back({t, r, ym, yMrac, yGs, yStr, uMrac, uGs, uStr,
                        thetaY, thetaR, alphaHat, betaHat});
        eMrac.push_back(yMrac - ym);
        eGs.push_back(yGs - ym);
        eStr.push_back(yStr - ym);

        if (k == steps - 1) break;

        const double ymNext = ym + dt * (-aModel * ym + bModel * r);
        const double yMracNext = yMrac + dt * (-aTrue * yMrac + bTrue * uMrac);
        const double yGsNext = yGs + dt * (-aTrue * yGs + bTrue * uGs);
        const double yStrNext = yStr + dt * (-aTrue * yStr + bTrue * uStr);

        thetaY = clampValue(thetaY + dt * gammaY * e * yMrac, -5.0, 5.0);
        thetaR = clampValue(thetaR - dt * gammaR * e * r, -5.0, 5.0);

        const std::array<double, 2> phi{yStr, uStr};
        const std::array<double, 2> Pphi{
            P[0][0] * phi[0] + P[0][1] * phi[1],
            P[1][0] * phi[0] + P[1][1] * phi[1]
        };
        const double denominator = forgetting + phi[0] * Pphi[0] + phi[1] * Pphi[1];
        const std::array<double, 2> gain{Pphi[0] / denominator, Pphi[1] / denominator};
        const double prediction = alphaHat * phi[0] + betaHat * phi[1];
        const double predictionError = yStrNext - prediction;
        alphaHat = clampValue(alphaHat + gain[0] * predictionError, 0.85, 1.02);
        betaHat = clampValue(betaHat + gain[1] * predictionError, 0.15 * dt, 2.5 * dt);

        std::array<std::array<double, 2>, 2> Pnew{};
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                const double correction = gain[i] * (phi[0] * P[0][j] + phi[1] * P[1][j]);
                Pnew[i][j] = (P[i][j] - correction) / forgetting;
            }
        }
        P = Pnew;
        ym = ymNext; yMrac = yMracNext; yGs = yGsNext; yStr = yStrNext;
    }

    std::ofstream file("Chapter4_Lesson5_results_cpp.csv");
    file << "t,r,ym,y_mrac,y_gain_scheduling,y_str,u_mrac,u_gain_scheduling,u_str,"
            "theta_y_mrac,theta_r_mrac,alpha_hat_str,beta_hat_str\n";
    file << std::setprecision(12);
    for (const Sample& s : data) {
        file << s.t << ',' << s.r << ',' << s.ym << ',' << s.yMrac << ',' << s.yGs << ','
             << s.yStr << ',' << s.uMrac << ',' << s.uGs << ',' << s.uStr << ','
             << s.thetaY << ',' << s.thetaR << ',' << s.alphaHat << ',' << s.betaHat << '\n';
    }

    std::cout << std::fixed << std::setprecision(6)
              << "MRAC RMS model-following error: " << rms(eMrac) << '\n'
              << "Gain-scheduling RMS model-following error: " << rms(eGs) << '\n'
              << "STR RMS model-following error: " << rms(eStr) << '\n'
              << "Saved Chapter4_Lesson5_results_cpp.csv\n";
    return 0;
}

The Java implementation mirrors the same equations and writes the simulation history to CSV without external dependencies.

Chapter4_Lesson5.java

// Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators.
// Build and run: javac Chapter4_Lesson5.java && java Chapter4_Lesson5

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

public final class Chapter4_Lesson5 {
    private static final double DT = 0.01;
    private static final double T_FINAL = 36.0;
    private static final double A_MODEL = 1.2;
    private static final double B_MODEL = 1.2;
    private static final double B_TRUE = 1.0;

    private record Sample(double t, double r, double ym, double yMrac, double yGs, double yStr,
                          double uMrac, double uGs, double uStr, double thetaY, double thetaR,
                          double alphaHat, double betaHat) {}

    private static double plantParameter(double t) {
        if (t < 12.0) return 0.8;
        if (t < 24.0) return 1.6;
        return 0.55;
    }

    private static double referenceSignal(double t) {
        if (t < 2.0) return 0.0;
        if (t < 18.0) return 1.0;
        if (t < 28.0) return -0.5;
        return 0.75;
    }

    private static double clamp(double x, double lo, double hi) {
        return Math.max(lo, Math.min(x, hi));
    }

    private static double interpolate(double x, double[] grid, double[] values) {
        if (x <= grid[0]) return values[0];
        if (x >= grid[grid.length - 1]) return values[values.length - 1];
        for (int i = 0; i < grid.length - 1; i++) {
            if (x <= grid[i + 1]) {
                double w = (x - grid[i]) / (grid[i + 1] - grid[i]);
                return (1.0 - w) * values[i] + w * values[i + 1];
            }
        }
        return values[values.length - 1];
    }

    private static double rms(List<Double> errors) {
        double sum = 0.0;
        for (double e : errors) sum += e * e;
        return Math.sqrt(sum / errors.size());
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        int steps = (int) Math.round(T_FINAL / DT) + 1;
        double lambdaModel = Math.exp(-A_MODEL * DT);
        double[] rhoGrid = {0.50, 0.80, 1.20, 1.60, 1.80};
        double[] thetaYGrid = new double[rhoGrid.length];
        double[] thetaRGrid = new double[rhoGrid.length];
        for (int i = 0; i < rhoGrid.length; i++) {
            thetaYGrid[i] = (A_MODEL - rhoGrid[i]) / B_TRUE;
            thetaRGrid[i] = B_MODEL / B_TRUE;
        }

        double ym = 0.0, yMrac = 0.0, yGs = 0.0, yStr = 0.0;
        double thetaY = 0.0, thetaR = 0.5;
        double alphaHat = Math.exp(-1.0 * DT), betaHat = 0.8 * DT;
        double[][] p = { {200.0, 0.0}, {0.0, 200.0} };
        double gammaY = 3.0, gammaR = 2.0, forgetting = 0.998;

        List<Sample> samples = new ArrayList<>(steps);
        List<Double> eMrac = new ArrayList<>(steps);
        List<Double> eGs = new ArrayList<>(steps);
        List<Double> eStr = new ArrayList<>(steps);

        for (int k = 0; k < steps; k++) {
            double t = k * DT;
            double r = referenceSignal(t);
            double aTrue = plantParameter(t);

            double e = yMrac - ym;
            double uMrac = thetaR * r - thetaY * yMrac;

            double thetaYGs = interpolate(aTrue, rhoGrid, thetaYGrid);
            double thetaRGs = interpolate(aTrue, rhoGrid, thetaRGrid);
            double uGs = thetaRGs * r - thetaYGs * yGs;

            double betaSafe = Math.abs(betaHat) < 0.15 * DT ? 0.15 * DT : betaHat;
            double kyStr = clamp((lambdaModel - alphaHat) / betaSafe, -12.0, 12.0);
            double krStr = clamp((1.0 - lambdaModel) / betaSafe, -12.0, 12.0);
            double uStr = kyStr * yStr + krStr * r;

            samples.add(new Sample(t, r, ym, yMrac, yGs, yStr, uMrac, uGs, uStr,
                                   thetaY, thetaR, alphaHat, betaHat));
            eMrac.add(yMrac - ym);
            eGs.add(yGs - ym);
            eStr.add(yStr - ym);

            if (k == steps - 1) break;

            double ymNext = ym + DT * (-A_MODEL * ym + B_MODEL * r);
            double yMracNext = yMrac + DT * (-aTrue * yMrac + B_TRUE * uMrac);
            double yGsNext = yGs + DT * (-aTrue * yGs + B_TRUE * uGs);
            double yStrNext = yStr + DT * (-aTrue * yStr + B_TRUE * uStr);

            thetaY = clamp(thetaY + DT * gammaY * e * yMrac, -5.0, 5.0);
            thetaR = clamp(thetaR - DT * gammaR * e * r, -5.0, 5.0);

            double[] phi = {yStr, uStr};
            double[] pPhi = {
                p[0][0] * phi[0] + p[0][1] * phi[1],
                p[1][0] * phi[0] + p[1][1] * phi[1]
            };
            double denominator = forgetting + phi[0] * pPhi[0] + phi[1] * pPhi[1];
            double[] gain = {pPhi[0] / denominator, pPhi[1] / denominator};
            double predictionError = yStrNext - (alphaHat * phi[0] + betaHat * phi[1]);
            alphaHat = clamp(alphaHat + gain[0] * predictionError, 0.85, 1.02);
            betaHat = clamp(betaHat + gain[1] * predictionError, 0.15 * DT, 2.5 * DT);

            double[][] pNew = new double[2][2];
            for (int i = 0; i < 2; i++) {
                for (int j = 0; j < 2; j++) {
                    double correction = gain[i] * (phi[0] * p[0][j] + phi[1] * p[1][j]);
                    pNew[i][j] = (p[i][j] - correction) / forgetting;
                }
            }
            p = pNew;
            ym = ymNext;
            yMrac = yMracNext;
            yGs = yGsNext;
            yStr = yStrNext;
        }

        try (BufferedWriter writer = new BufferedWriter(new FileWriter("Chapter4_Lesson5_results_java.csv"))) {
            writer.write("t,r,ym,y_mrac,y_gain_scheduling,y_str,u_mrac,u_gain_scheduling,u_str," +
                         "theta_y_mrac,theta_r_mrac,alpha_hat_str,beta_hat_str\n");
            for (Sample s : samples) {
                writer.write(String.format(Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    s.t(), s.r(), s.ym(), s.yMrac(), s.yGs(), s.yStr(), s.uMrac(), s.uGs(), s.uStr(),
                    s.thetaY(), s.thetaR(), s.alphaHat(), s.betaHat()));
            }
        }

        System.out.printf("MRAC RMS model-following error: %.6f%n", rms(eMrac));
        System.out.printf("Gain-scheduling RMS model-following error: %.6f%n", rms(eGs));
        System.out.printf("STR RMS model-following error: %.6f%n", rms(eStr));
        System.out.println("Saved Chapter4_Lesson5_results_java.csv");
    }
}

The MATLAB script executes the common comparison, reports RMS tracking errors, writes a table to CSV, and plots outputs and control inputs.

Chapter4_Lesson5.m

% Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators
% The same uncertain first-order plant is controlled by all three methods.

clear; clc; close all;

dt = 0.01;
tFinal = 36;
t = (0:dt:tFinal)';
N = numel(t);

aModel = 1.2;
bModel = 1.2;
bTrue = 1.0;
lambdaModel = exp(-aModel * dt);

r = zeros(N,1);
r(t >= 2 & t < 18) = 1.0;
r(t >= 18 & t < 28) = -0.5;
r(t >= 28) = 0.75;

aTrue = 0.8 * ones(N,1);
aTrue(t >= 12 & t < 24) = 1.6;
aTrue(t >= 24) = 0.55;

ym = zeros(N,1);
yMrac = zeros(N,1);
yGs = zeros(N,1);
yStr = zeros(N,1);
uMrac = zeros(N,1);
uGs = zeros(N,1);
uStr = zeros(N,1);

% Direct MRAC parameters.
thetaY = zeros(N,1);
thetaR = zeros(N,1);
thetaR(1) = 0.5;
gammaY = 3.0;
gammaR = 2.0;

% Offline gain schedule, indexed by the measurable variable rho=a(t).
rhoGrid = [0.50 0.80 1.20 1.60 1.80];
thetaYGrid = (aModel - rhoGrid) / bTrue;
thetaRGrid = (bModel / bTrue) * ones(size(rhoGrid));

% STR's online discrete-time parameter estimator.
alphaHat = zeros(N,1);
betaHat = zeros(N,1);
alphaHat(1) = exp(-1.0 * dt);
betaHat(1) = 0.8 * dt;
P = 200 * eye(2);
forgetting = 0.998;

for k = 1:N-1
    ym(k+1) = ym(k) + dt * (-aModel * ym(k) + bModel * r(k));

    % MRAC.
    e = yMrac(k) - ym(k);
    uMrac(k) = thetaR(k) * r(k) - thetaY(k) * yMrac(k);
    thetaY(k+1) = min(max(thetaY(k) + dt * gammaY * e * yMrac(k), -5), 5);
    thetaR(k+1) = min(max(thetaR(k) - dt * gammaR * e * r(k), -5), 5);
    yMrac(k+1) = yMrac(k) + dt * (-aTrue(k) * yMrac(k) + bTrue * uMrac(k));

    % Gain scheduling.
    thetaYGs = interp1(rhoGrid, thetaYGrid, aTrue(k), 'linear', 'extrap');
    thetaRGs = interp1(rhoGrid, thetaRGrid, aTrue(k), 'linear', 'extrap');
    uGs(k) = thetaRGs * r(k) - thetaYGs * yGs(k);
    yGs(k+1) = yGs(k) + dt * (-aTrue(k) * yGs(k) + bTrue * uGs(k));

    % STR: estimate [alpha beta], then update pole-placement gains.
    betaSafe = betaHat(k);
    if abs(betaSafe) < 0.15 * dt
        betaSafe = 0.15 * dt;
    end
    kyStr = min(max((lambdaModel - alphaHat(k)) / betaSafe, -12), 12);
    krStr = min(max((1 - lambdaModel) / betaSafe, -12), 12);
    uStr(k) = kyStr * yStr(k) + krStr * r(k);
    yStr(k+1) = yStr(k) + dt * (-aTrue(k) * yStr(k) + bTrue * uStr(k));

    phi = [yStr(k); uStr(k)];
    K = P * phi / (forgetting + phi' * P * phi);
    predictionError = yStr(k+1) - phi' * [alphaHat(k); betaHat(k)];
    estimate = [alphaHat(k); betaHat(k)] + K * predictionError;
    P = (P - K * phi' * P) / forgetting;
    alphaHat(k+1) = min(max(estimate(1), 0.85), 1.02);
    betaHat(k+1) = min(max(estimate(2), 0.15 * dt), 2.5 * dt);
end

uMrac(end) = uMrac(end-1);
uGs(end) = uGs(end-1);
uStr(end) = uStr(end-1);

rmsMrac = sqrt(mean((yMrac - ym).^2));
rmsGs = sqrt(mean((yGs - ym).^2));
rmsStr = sqrt(mean((yStr - ym).^2));
fprintf('MRAC RMS model-following error: %.6f\n', rmsMrac);
fprintf('Gain-scheduling RMS model-following error: %.6f\n', rmsGs);
fprintf('STR RMS model-following error: %.6f\n', rmsStr);

results = table(t,r,ym,yMrac,yGs,yStr,uMrac,uGs,uStr,thetaY,thetaR,alphaHat,betaHat);
writetable(results, 'Chapter4_Lesson5_results_matlab.csv');

figure;
plot(t,r,':',t,ym,'--',t,yMrac,t,yGs,t,yStr,'LineWidth',1.2);
grid on; xlabel('Time [s]'); ylabel('Output');
legend('command r','reference model','MRAC','gain scheduling','STR','Location','best');
title('MRAC vs Gain Scheduling vs STR');

figure;
plot(t,uMrac,t,uGs,t,uStr,'LineWidth',1.2);
grid on; xlabel('Time [s]'); ylabel('Control input');
legend('MRAC','gain scheduling','STR','Location','best');
title('Control effort comparison');

The following MATLAB script creates a fixed-step Simulink model with one discrete MATLAB Function block for each architecture. Each block contains its controller, plant state, and adaptive or estimator state. The model is saved as Chapter4_Lesson5_Simulink_Model.slx and simulated for 36 seconds.

Chapter4_Lesson5_Simulink.m

% Chapter 4, Lesson 5: programmatically construct a discrete Simulink comparison.
% Requires Simulink. The three MATLAB Function blocks contain controller,
% plant, and adaptation/estimation states, all updated at Ts = 0.01 s.

clear; clc;
load_system('simulink');
model = 'Chapter4_Lesson5_Simulink_Model';
Ts = 0.01;

if bdIsLoaded(model)
    close_system(model, 0);
end
new_system(model);
open_system(model);
set_param(model, 'SolverType', 'Fixed-step', 'Solver', 'FixedStepDiscrete', ...
    'FixedStep', num2str(Ts), 'StopTime', '36');

add_block('simulink/Sources/Step', [model '/Command'], ...
    'Time', '2', 'Before', '0', 'After', '1', 'SampleTime', num2str(Ts), ...
    'Position', [35 105 65 135]);

blocks = {'MRAC','GainScheduling','STR'};
yPositions = [45 155 265];
for i = 1:numel(blocks)
    add_block('simulink/User-Defined Functions/MATLAB Function', ...
        [model '/' blocks{i}], 'Position', [145 yPositions(i) 330 yPositions(i)+70]);
    add_block('simulink/Sinks/To Workspace', [model '/y_' blocks{i}], ...
        'VariableName', ['y' blocks{i}], 'SaveFormat', 'Structure With Time', ...
        'Position', [390 yPositions(i)+5 485 yPositions(i)+30]);
    add_line(model, 'Command/1', [blocks{i} '/1'], 'autorouting', 'on');
    add_line(model, [blocks{i} '/1'], ['y_' blocks{i} '/1'], 'autorouting', 'on');
end

mracCode = sprintf([ ...
'function y = fcn(r)\n' ...
'%% Direct MRAC plus uncertain first-order plant\n' ...
'persistent x ym thY thR k\n' ...
'if isempty(k), x=0; ym=0; thY=0; thR=0.5; k=0; end\n' ...
'Ts=0.01; am=1.2; bm=1.2; gy=3; gr=2; t=k*Ts;\n' ...
'if t>=18 && t<28, r=-0.5; elseif t>=28, r=0.75; end\n' ...
'if t<12, a=0.8; elseif t<24, a=1.6; else, a=0.55; end\n' ...
'e=x-ym; u=thR*r-thY*x;\n' ...
'thY=min(max(thY+Ts*gy*e*x,-5),5);\n' ...
'thR=min(max(thR-Ts*gr*e*r,-5),5);\n' ...
'x=x+Ts*(-a*x+u); ym=ym+Ts*(-am*ym+bm*r);\n' ...
'k=k+1; y=x;\n']);

gsCode = sprintf([ ...
'function y = fcn(r)\n' ...
'%% Gain scheduling with exact measured scheduling variable rho=a(t)\n' ...
'persistent x k\n' ...
'if isempty(k), x=0; k=0; end\n' ...
'Ts=0.01; am=1.2; bm=1.2; t=k*Ts;\n' ...
'if t>=18 && t<28, r=-0.5; elseif t>=28, r=0.75; end\n' ...
'if t<12, rho=0.8; elseif t<24, rho=1.6; else, rho=0.55; end\n' ...
'ky=am-rho; kr=bm; u=kr*r-ky*x;\n' ...
'x=x+Ts*(-rho*x+u); k=k+1; y=x;\n']);

strCode = sprintf([ ...
'function y = fcn(r)\n' ...
'%% STR with scalar RLS identifier and certainty-equivalent pole placement\n' ...
'persistent x ah bh P k\n' ...
'if isempty(k), x=0; ah=exp(-0.01); bh=0.008; P=200*eye(2); k=0; end\n' ...
'Ts=0.01; am=1.2; lm=exp(-am*Ts); ff=0.998; t=k*Ts;\n' ...
'if t>=18 && t<28, r=-0.5; elseif t>=28, r=0.75; end\n' ...
'if t<12, a=0.8; elseif t<24, a=1.6; else, a=0.55; end\n' ...
'bs=max(bh,0.15*Ts); ky=min(max((lm-ah)/bs,-12),12);\n' ...
'kr=min(max((1-lm)/bs,-12),12); u=ky*x+kr*r;\n' ...
'xOld=x; x=x+Ts*(-a*x+u); phi=[xOld;u];\n' ...
'K=P*phi/(ff+phi''*P*phi); est=[ah;bh]+K*(x-phi''*[ah;bh]);\n' ...
'P=(P-K*phi''*P)/ff; ah=min(max(est(1),0.85),1.02);\n' ...
'bh=min(max(est(2),0.15*Ts),2.5*Ts); k=k+1; y=x;\n']);

root = sfroot;
mracChart = root.find('-isa', 'Stateflow.EMChart', 'Path', [model '/MRAC']);
gsChart = root.find('-isa', 'Stateflow.EMChart', 'Path', [model '/GainScheduling']);
strChart = root.find('-isa', 'Stateflow.EMChart', 'Path', [model '/STR']);
mracChart.Script = mracCode;
gsChart.Script = gsCode;
strChart.Script = strCode;

save_system(model, [model '.slx']);
sim(model);
disp(['Created and simulated ' model '.slx']);

The Wolfram Language file contains the complete simulation. The notebook file is a lightweight launcher that loads the adjacent implementation, making the code maintainable while still providing a native .nb entry point.

Chapter4_Lesson5.wl

(* Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators *)
ClearAll["Global`*"];

dt = 0.01; tFinal = 36.; n = Round[tFinal/dt] + 1;
aModel = 1.2; bModel = 1.2; bTrue = 1.; lambdaModel = Exp[-aModel dt];
plantA[t_] := Piecewise[{ {0.8, t < 12.}, {1.6, t < 24.} }, 0.55];
reference[t_] := Piecewise[{ {0., t < 2.}, {1., t < 18.}, {-0.5, t < 28.} }, 0.75];
clip[x_, lo_, hi_] := Min[Max[x, lo], hi];

rhoGrid = {0.50, 0.80, 1.20, 1.60, 1.80};
thetaYGrid = (aModel - #)/bTrue & /@ rhoGrid;
thetaRGrid = ConstantArray[bModel/bTrue, Length[rhoGrid]];
thetaYInterpolation = Interpolation[Transpose[{rhoGrid, thetaYGrid}], InterpolationOrder -> 1];
thetaRInterpolation = Interpolation[Transpose[{rhoGrid, thetaRGrid}], InterpolationOrder -> 1];

state = <|"ym" -> 0., "yMrac" -> 0., "yGs" -> 0., "yStr" -> 0.,
  "thetaY" -> 0., "thetaR" -> 0.5, "alphaHat" -> Exp[-dt],
  "betaHat" -> 0.8 dt, "P" -> 200. IdentityMatrix[2]|>;
forgetting = 0.998; gammaY = 3.; gammaR = 2.;

rows = Reap[
   Do[
    t = (k - 1) dt; r = reference[t]; aTrue = plantA[t];
    e = state["yMrac"] - state["ym"];
    uMrac = state["thetaR"] r - state["thetaY"] state["yMrac"];
    kyGs = thetaYInterpolation[aTrue]; krGs = thetaRInterpolation[aTrue];
    uGs = krGs r - kyGs state["yGs"];
    betaSafe = Max[state["betaHat"], 0.15 dt];
    kyStr = clip[(lambdaModel - state["alphaHat"])/betaSafe, -12., 12.];
    krStr = clip[(1. - lambdaModel)/betaSafe, -12., 12.];
    uStr = kyStr state["yStr"] + krStr r;
    Sow[{t, r, state["ym"], state["yMrac"], state["yGs"], state["yStr"],
      uMrac, uGs, uStr, state["thetaY"], state["thetaR"],
      state["alphaHat"], state["betaHat"]}];
    If[k < n,
     ymNext = state["ym"] + dt (-aModel state["ym"] + bModel r);
     yMracNext = state["yMrac"] + dt (-aTrue state["yMrac"] + bTrue uMrac);
     yGsNext = state["yGs"] + dt (-aTrue state["yGs"] + bTrue uGs);
     yStrNext = state["yStr"] + dt (-aTrue state["yStr"] + bTrue uStr);
     thetaYNext = clip[state["thetaY"] + dt gammaY e state["yMrac"], -5., 5.];
     thetaRNext = clip[state["thetaR"] - dt gammaR e r, -5., 5.];
     phi = {state["yStr"], uStr};
     gain = state["P"].phi/(forgetting + phi.state["P"].phi);
     predictionError = yStrNext - phi.{state["alphaHat"], state["betaHat"]};
     estimate = {state["alphaHat"], state["betaHat"]} + gain predictionError;
     pNext = (state["P"] - Outer[Times, gain, phi].state["P"])/forgetting;
     state = <|"ym" -> ymNext, "yMrac" -> yMracNext, "yGs" -> yGsNext,
       "yStr" -> yStrNext, "thetaY" -> thetaYNext, "thetaR" -> thetaRNext,
       "alphaHat" -> clip[estimate[[1]], 0.85, 1.02],
       "betaHat" -> clip[estimate[[2]], 0.15 dt, 2.5 dt], "P" -> pNext|>;
     ],
    {k, 1, n}
    ]][[2, 1]];

header = {"t", "r", "ym", "y_mrac", "y_gain_scheduling", "y_str", "u_mrac",
   "u_gain_scheduling", "u_str", "theta_y_mrac", "theta_r_mrac",
   "alpha_hat_str", "beta_hat_str"};
Export[FileNameJoin[{DirectoryName[$InputFileName], "Chapter4_Lesson5_results_wolfram.csv"}],
  Prepend[rows, header]];

rms[v_] := Sqrt[Mean[v^2]];
Print["MRAC RMS model-following error: ", N[rms[rows[[All, 4]] - rows[[All, 3]]], 7]];
Print["Gain-scheduling RMS model-following error: ", N[rms[rows[[All, 5]] - rows[[All, 3]]], 7]];
Print["STR RMS model-following error: ", N[rms[rows[[All, 6]] - rows[[All, 3]]], 7]];

outputPlot = ListLinePlot[
  {rows[[All, {1, 2}]], rows[[All, {1, 3}]], rows[[All, {1, 4}]],
   rows[[All, {1, 5}]], rows[[All, {1, 6}]]},
  PlotLegends -> {"command r", "reference model", "MRAC", "gain scheduling", "STR"},
  Frame -> True, FrameLabel -> {"Time [s]", "Output"}, PlotRange -> All,
  ImageSize -> Large];
controlPlot = ListLinePlot[
  {rows[[All, {1, 7}]], rows[[All, {1, 8}]], rows[[All, {1, 9}]]},
  PlotLegends -> {"MRAC", "gain scheduling", "STR"}, Frame -> True,
  FrameLabel -> {"Time [s]", "Control input"}, PlotRange -> All, ImageSize -> Large];
Column[{outputPlot, controlPlot}]

Chapter4_Lesson5.nb

Notebook[{
 Cell["Chapter 4, Lesson 5: MRAC vs Gain Scheduling vs Self-Tuning Regulators", "Title"],
 Cell["This notebook launches the complete Wolfram Language implementation stored beside it in Chapter4_Lesson5.wl.", "Text"],
 Cell[BoxData[
  RowBox[{"Get", "[", 
   RowBox[{"FileNameJoin", "[", 
    RowBox[{"{", 
     RowBox[{"NotebookDirectory", "[", "]", ",", "\"Chapter4_Lesson5.wl\""}], 
     "}"}], "]"}], "]"}]], "Input"]
}, WindowTitle -> "Chapter4_Lesson5", StyleDefinitions -> "Default.nb"]

Problem 1 (MRAC Lyapunov redesign): For \(\dot y=-ay+bu\), \(b>0\), the reference model \(\dot y_m=-a_my_m+b_mr\), and \(u=\theta_r r-\theta_y y\), derive adaptive laws that make \(\dot V\leq0\) for a quadratic state-parameter Lyapunov function.

Solution: Matching gives

\[ \theta_y^*=\frac{a_m-a}{b},\qquad \theta_r^*=\frac{b_m}{b}. \]

With \(e=y-y_m\),

\[ \dot e=-a_me-b\tilde\theta_y y+b\tilde\theta_r r. \]

Choose

\[ V=\frac12e^2+\frac{b}{2\gamma_y}\tilde\theta_y^2 +\frac{b}{2\gamma_r}\tilde\theta_r^2. \]

The derivative contains the products \(b\tilde\theta_y(-ey+\dot\theta_y/\gamma_y)\) and \(b\tilde\theta_r(er+\dot\theta_r/\gamma_r)\). Setting

\[ \dot\theta_y=\gamma_yey,\qquad \dot\theta_r=-\gamma_rer \]

cancels both products and yields \(\dot V=-a_me^2\leq0\).

Problem 2 (Interpolated gain schedule): Let \(a_m=1.2\), \(b=1\), and use grid values \(\rho_1=0.8\) and \(\rho_2=1.6\). Compute the interpolated feedback coefficient \(\theta_y(1.1)\).

Solution: The endpoint gains are

\[ \theta_y(0.8)=1.2-0.8=0.4,\qquad \theta_y(1.6)=1.2-1.6=-0.4. \]

The interpolation weight is

\[ \alpha=\frac{1.1-0.8}{1.6-0.8}=0.375. \]

Therefore

\[ \theta_y(1.1)=(1-0.375)(0.4)+0.375(-0.4)=0.1. \]

Since the exact matching formula is \(1.2-1.1=0.1\), linear interpolation is exact for this affine scalar schedule.

Problem 3 (STR certainty-equivalent gains): A sampled model is \(y_{k+1}=\alpha y_k+\beta u_k\). At one instant, \(\hat\alpha=0.93\), \(\hat\beta=0.04\), and the desired pole is \(\lambda_m=0.90\). Determine \(k_y\) and \(k_r\) for \(u_k=k_yy_k+k_rr_k\).

Solution:

\[ k_y=\frac{0.90-0.93}{0.04}=-0.75,\qquad k_r=\frac{1-0.90}{0.04}=2.5. \]

The estimated closed-loop equation is then \(y_{k+1}=0.90y_k+0.10r_k\). If the estimate is inaccurate, the actual coefficients differ from these desired values.

Problem 4 (Rate bound for a schedule): Suppose \(\lambda_{\min}(Q)=2\) and \(\sup_\rho\|\partial P/\partial\rho\|=5\). Give a strict sufficient bound on \(|\dot\rho|\) for the Lyapunov estimate in Section 4 to remain negative definite.

Solution:

\[ |\dot\rho|<\frac{2}{5}=0.4. \]

At equality the simple estimate is only semidefinite, so a strict inequality is required for this sufficient argument.

Problem 5 (Architecture selection): Classify the most natural starting architecture for each case: (a) aircraft gains are accurately indexed by measured airspeed and altitude; (b) the plant has a known matching structure and the main requirement is to follow a specified reference model despite an unknown coefficient; (c) a process model must be estimated online and reused to recompute a pole-placement controller.

Solution:

• (a) Gain scheduling, because the operating condition is measured and can index an offline controller family.
• (b) MRAC, because the desired model and matching parameterization directly support error-driven adaptation.
• (c) STR, because explicit online plant estimates are required before controller redesign.

These are starting choices, not unconditional guarantees. Each design must still be checked for its own stability and implementation assumptions.

MRAC, gain scheduling, and self-tuning regulation all produce time-varying controllers, but they differ in the source of that variation. MRAC changes controller coefficients directly from model-following error and admits a combined error-parameter Lyapunov proof. Gain scheduling changes gains from a measured operating variable and can respond immediately, but stability of frozen controllers does not alone prove safe transitions. An STR estimates a plant model online and applies a certainty-equivalent design map, gaining an explicit model at the cost of estimator-controller coupling and conditioning risks. The correct choice depends on whether plant changes are measured, whether a stable matching parameterization is available, and whether explicit online identification is required.