After completing this lesson, students should be able to:

• compute the initial Lyapunov energy associated with state and parameter mismatch;
• derive conservative transient bounds from the quadratic Lyapunov function;
• separate reference-model, tracking-error, parameter, and saturation transients;
• show how saturation changes the nominal MRAC error dynamics and Lyapunov derivative;
• interpret command limiting and adaptation freezing as practical safeguards rather than automatic proofs of stability;
• perform matched, mismatched, and saturated simulation experiments in several programming environments.

The state-feedback second-order MRAC structure, matching equations, Lyapunov equation, and nominal update law are assumed from Lessons 1–4. Formal robust modifications, normalization, and rigorous anti-saturation adaptive designs are intentionally deferred to later chapters.

Consider the uncertain second-order plant \( \mathbf{x}=[x_1\;x_2]^\mathsf{T} \):

\[ \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{b}u_a,\qquad \mathbf{A}= \begin{bmatrix}0&1\\a_1&a_2\end{bmatrix},\qquad \mathbf{b}=\begin{bmatrix}0\\b\end{bmatrix},\qquad b>0. \]

The desired dynamics are generated by the Hurwitz reference model

\[ \dot{\mathbf{x} }_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{b}_m r,\qquad \mathbf{A}_m= \begin{bmatrix}0&1\\-\omega_n^2&-2\zeta\omega_n\end{bmatrix},\qquad \mathbf{b}_m=\begin{bmatrix}0\\\omega_n^2\end{bmatrix}. \]

The commanded adaptive control is linearly parameterized:

\[ u_c=\boldsymbol{\theta}^\mathsf{T}\boldsymbol{\phi},\qquad \boldsymbol{\phi}= \begin{bmatrix}x_1&x_2&r\end{bmatrix}^\mathsf{T},\qquad \boldsymbol{\theta}= \begin{bmatrix}\theta_1&\theta_2&\theta_r\end{bmatrix}^\mathsf{T}. \]

In the nominal unconstrained derivation, \(u_a=u_c\). Ideal parameters satisfy

\[ \mathbf{A}+\mathbf{b}\boldsymbol{\theta}_x^{*\mathsf{T} } =\mathbf{A}_m,\qquad \mathbf{b}\theta_r^*=\mathbf{b}_m. \]

With \(\mathbf{e}=\mathbf{x}-\mathbf{x}_m\) and \(\widetilde{\boldsymbol{\theta} } =\boldsymbol{\theta}-\boldsymbol{\theta}^*\), the nominal error model becomes

\[ \dot{\mathbf{e} } =\mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^\mathsf{T} \boldsymbol{\phi}. \]

Choose \(\mathbf{P}=\mathbf{P}^\mathsf{T}>0\) from

\[ \mathbf{A}_m^\mathsf{T}\mathbf{P} +\mathbf{P}\mathbf{A}_m=-\mathbf{Q},\qquad \mathbf{Q}=\mathbf{Q}^\mathsf{T}>0, \]

and use the update law

\[ \dot{\boldsymbol{\theta} } =-\boldsymbol{\Gamma}\boldsymbol{\phi}\, \mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b},\qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^\mathsf{T}>0. \]

flowchart TD
  R["Reference command r"] --> RM["Second-order reference model"]
  R --> REG["Regressor: x1, x2, r"]
  PL["Plant states x1, x2"] --> REG
  REG --> AC["Adaptive parameter vector"]
  AC --> UC["Commanded control uc"]
  UC --> SAT["Actuator limit"]
  SAT --> UA["Applied control ua"]
  UA --> PL
  RM --> ERR["Tracking error e = x - xm"]
  PL --> ERR
  ERR --> AC
        

For the nominal system, use the combined state-parameter Lyapunov function

\[ V(\mathbf{e},\widetilde{\boldsymbol{\theta} }) =\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{e} +\widetilde{\boldsymbol{\theta} }^\mathsf{T} \boldsymbol{\Gamma}^{-1} \widetilde{\boldsymbol{\theta} }. \]

Its initial value is

\[ V(0)= \bigl(\mathbf{x}(0)-\mathbf{x}_m(0)\bigr)^\mathsf{T} \mathbf{P} \bigl(\mathbf{x}(0)-\mathbf{x}_m(0)\bigr) +\bigl(\boldsymbol{\theta}(0)-\boldsymbol{\theta}^*\bigr)^\mathsf{T} \boldsymbol{\Gamma}^{-1} \bigl(\boldsymbol{\theta}(0)-\boldsymbol{\theta}^*\bigr). \]

Thus two independent mismatches inject transient energy: state mismatch \(\mathbf{e}(0)\) and parameter mismatch \(\widetilde{\boldsymbol{\theta} }(0)\). Under the nominal assumptions, the update law gives

\[ \dot{V}=-\mathbf{e}^\mathsf{T}\mathbf{Q}\mathbf{e}\le 0. \]

Therefore \(V(t)\le V(0)\), which immediately yields conservative pointwise bounds:

\[ \|\mathbf{e}(t)\| \le \sqrt{\frac{V(0)}{\lambda_{\min}(\mathbf{P})} }, \qquad \|\widetilde{\boldsymbol{\theta} }(t)\| \le \sqrt{\lambda_{\max}(\boldsymbol{\Gamma})\,V(0)}. \]

These are energy bounds, not tight overshoot predictions. They explain why aggressive adaptation cannot erase a large initial state mismatch instantaneously: the plant and reference model still obey differential equations with finite rates.

3.1 Matching the Initial Reference-Model State

If the desired experiment begins from the measured plant state, selecting \(\mathbf{x}_m(0)=\mathbf{x}(0)\) removes the initial tracking-error term. This is often appropriate for command-following tests. It is not appropriate when the reference-model initial state itself represents a prescribed desired condition that differs from the plant.

3.2 Parameter Initialization

Zero parameter initialization is mathematically permissible under the nominal proof, but it can demand a large adaptation excursion. A defensible nominal model may be used to initialize \(\boldsymbol{\theta}(0)\). The initial estimate should be documented, physically plausible, and independent of the data later used to claim performance.

3.3 Reference-Model Initial Response

Even perfect adaptive tracking does not eliminate the reference model's own zero-input transient:

\[ \mathbf{x}_m(t) =e^{\mathbf{A}_m t}\mathbf{x}_m(0) +\int_0^t e^{\mathbf{A}_m(t-\tau)}\mathbf{b}_m r(\tau)\,d\tau. \]

Consequently, overshoot observed in \(\mathbf{x}(t)\) may originate from the selected damping ratio and natural frequency rather than from adaptation.

A practical MRAC transient contains at least four interacting components:

1. Reference-model transient: determined by \(\mathbf{A}_m\), \(\mathbf{b}_m\), \(r\), and \(\mathbf{x}_m(0)\).
2. Tracking-error transient: determined by \(\mathbf{e}(0)\) and the adaptive coupling.
3. Parameter transient: determined by \(\widetilde{\boldsymbol{\theta} }(0)\), \(\boldsymbol{\Gamma}\), and regressor magnitude.
4. Constraint transient: produced when \(u_c\) differs from the applied actuator signal \(u_a\).

4.1 Initial Parameter Rate

The instantaneous update rate is

\[ \dot{\boldsymbol{\theta} }(0) =-\boldsymbol{\Gamma}\boldsymbol{\phi}(0) \mathbf{e}^\mathsf{T}(0)\mathbf{P}\mathbf{b}. \]

Large initial error, large regressor components, or large adaptation gains can therefore create a sharp parameter-rate transient. In norm form,

\[ \|\dot{\boldsymbol{\theta} }(0)\| \le \|\boldsymbol{\Gamma}\|\, \|\boldsymbol{\phi}(0)\|\, \|\mathbf{P}\mathbf{b}\|\, \|\mathbf{e}(0)\|. \]

4.2 Adaptation Gain Is Not a Conventional Feedback Bandwidth

Increasing \(\boldsymbol{\Gamma}\) accelerates parameter motion, but does not simply move fixed closed-loop poles. Excessive values can amplify measurement noise, worsen numerical stiffness, and produce a large commanded input before the plant has responded. Tuning should therefore be based on state scaling, expected regressor magnitude, actuator authority, and integration step size.

4.3 Integral Error Measures

The nominal derivative identity implies

\[ \int_0^T \mathbf{e}^\mathsf{T}\mathbf{Q}\mathbf{e}\,dt =V(0)-V(T)\le V(0). \]

This relation provides a useful implementation check: in a high-accuracy unconstrained simulation, numerical evaluation of the left side should not substantially exceed the computed initial Lyapunov value.

flowchart TD
  S["Observed large transient"] --> A["Check reference-model response alone"]
  A --> B["Check initial state mismatch"]
  B --> C["Check initial parameter estimate and regressor scale"]
  C --> D["Compare commanded and applied control"]
  D --> E["Measure saturation duration"]
  E --> F["Reduce command rate or model aggressiveness"]
  F --> G["Retune adaptation gain and integration step"]
  G --> H["Repeat matched, mismatched, and saturated tests"]
        

Let the physical actuator apply

\[ u_a=\operatorname{sat}(u_c,u_{\max}) = \begin{cases} u_{\max},&u_c>u_{\max},\\ u_c,&|u_c|\le u_{\max},\\ -u_{\max},&u_c<-u_{\max}. \end{cases} \]

Define the saturation discrepancy

\[ \Delta u=u_a-u_c. \]

The nominal equality \(u_a=u_c\) is now false, and the error dynamics become

\[ \dot{\mathbf{e} } =\mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^\mathsf{T} \boldsymbol{\phi} +\mathbf{b}\Delta u. \]

If the nominal update law is retained, differentiation of the same Lyapunov function gives

\[ \dot{V} =-\mathbf{e}^\mathsf{T}\mathbf{Q}\mathbf{e} +2\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b}\,\Delta u. \]

The additional term has no fixed sign. Therefore the nominal proof does not establish monotonic decay of \(V\) during saturation.

5.1 A Practical Boundedness Inequality

For any \(0<\varepsilon<\lambda_{\min}(\mathbf{Q})\), Young's inequality gives

\[ 2\left|\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b}\Delta u\right| \le \varepsilon\|\mathbf{e}\|^2 +\frac{\|\mathbf{P}\mathbf{b}\|^2}{\varepsilon}|\Delta u|^2. \]

Hence

\[ \dot{V} \le -\bigl(\lambda_{\min}(\mathbf{Q})-\varepsilon\bigr) \|\mathbf{e}\|^2 +\frac{\|\mathbf{P}\mathbf{b}\|^2}{\varepsilon}|\Delta u|^2. \]

This inequality explains practical ultimate-bounded behavior when \(\Delta u\) is bounded, but it is not by itself a complete theorem for the adaptive parameters. Persistent saturation can maintain tracking error while the nominal law continues changing the parameter estimates.

5.2 Parameter Drift Mechanism

Under sustained saturation, the plant cannot generate the acceleration requested by the reference model. The error remains nonzero, so

\[ \dot{\boldsymbol{\theta} } =-\boldsymbol{\Gamma}\boldsymbol{\phi}\, \mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b} \]

may continue integrating even though additional commanded control cannot reach the plant. This is the adaptive analogue of windup. The code examples compare plain adaptation with a simple freeze-while-saturated safeguard. That safeguard can reduce parameter growth, but it changes the adaptive law and is not claimed here as a general stability proof.

6.1 Make the Reference Model Achievable

For the reference model, the requested acceleration is

\[ \dot{x}_{m2} =-\omega_n^2x_{m1}-2\zeta\omega_n x_{m2}+\omega_n^2r. \]

A high natural frequency can demand more actuator authority than the plant possesses. Before adaptation is enabled, compare the required nominal input over representative commands with the physical actuator limit.

6.2 Shape Commands Before They Reach the Adaptive Loop

Step commands are valuable for analysis but can be unnecessarily severe in hardware. A first-order command filter

\[ \dot{r}_f=-\omega_f r_f+\omega_f r \]

or a rate limiter can reduce the initial control demand. This is command shaping, not a change to the MRAC stability proof.

6.3 Scale States and Regressors

If state magnitudes differ greatly, use dimensionally meaningful scaled variables. Otherwise a single diagonal adaptation gain can cause one parameter to update far faster than the others.

6.4 Record Both Commanded and Applied Control

Logging only \(u_c\) hides the physical nonlinearity. Always record \(u_c\), \(u_a\), \(\Delta u\), saturation duration, tracking error, and parameter norm.

6.5 Use a Test Matrix

1. matched plant/reference initial states with a loose actuator limit;
2. mismatched initial states with a loose actuator limit;
3. the same mismatch with the true actuator limit;
4. the saturated case with the selected practical safeguard;
5. parameter and state perturbations near the expected operating envelope.

The implementations use

\[ \mathbf{A}= \begin{bmatrix}0&1\\-1&-0.8\end{bmatrix},\qquad \mathbf{b}=\begin{bmatrix}0\\1.2\end{bmatrix}, \]

\[ \mathbf{A}_m= \begin{bmatrix}0&1\\-4&-3.2\end{bmatrix},\qquad \mathbf{b}_m=\begin{bmatrix}0\\4\end{bmatrix}. \]

The matching equations give

\[ \theta_1^*=\frac{-4-(-1)}{1.2}=-2.5,\qquad \theta_2^*=\frac{-3.2-(-0.8)}{1.2}=-2,\qquad \theta_r^*=\frac{4}{1.2}=\frac{10}{3}. \]

With \(\mathbf{Q}=\mathbf{I}\), the Lyapunov equation yields

\[ \mathbf{P}= \begin{bmatrix} 1.18125&0.125\\ 0.125&0.1953125 \end{bmatrix}. \]

Four cases are simulated. The saturation limit \(u_{\max}=1\) is deliberately restrictive so that commanded/applied-input mismatch and parameter growth are visible. Parameter convergence to the ideal values is not expected from a single step command because that command generally does not provide persistent excitation.

This implementation uses NumPy for vector operations and Matplotlib for plots. A fixed-step fourth-order Runge-Kutta method integrates the combined plant, reference-model, and adaptive-parameter states.

Chapter6_Lesson5.py

# Chapter6_Lesson5.py
# Second-order MRAC: initial-condition transients and actuator saturation.
# Dependencies: numpy, matplotlib

from dataclasses import dataclass
from pathlib import Path
import csv
import numpy as np
import matplotlib.pyplot as plt


@dataclass(frozen=True)
class Case:
    name: str
    x0: tuple[float, float]
    xm0: tuple[float, float]
    u_max: float
    freeze_when_saturated: bool


A_P1 = -1.0
A_P2 = -0.8
B = 1.2

A_M1 = -4.0
A_M2 = -3.2
B_M = 4.0

P = np.array([[1.18125, 0.125],
              [0.125,   0.1953125]], dtype=float)
GAMMA = np.array([4.0, 2.0, 3.0], dtype=float)

DT = 0.002
T_FINAL = 15.0
OUTPUT_DIR = Path("Chapter6_Lesson5_Output")


def reference(t: float) -> float:
    return 0.0 if t < 0.5 else 1.5


def saturation(value: float, limit: float) -> float:
    return float(np.clip(value, -limit, limit))


def derivative(t: float, z: np.ndarray, case: Case) -> tuple[np.ndarray, dict[str, float]]:
    x = z[0:2]
    xm = z[2:4]
    theta = z[4:7]

    r = reference(t)
    phi = np.array([x[0], x[1], r], dtype=float)
    u_cmd = float(theta @ phi)
    u_act = saturation(u_cmd, case.u_max)

    e = x - xm
    pb = P @ np.array([0.0, B], dtype=float)
    sigma = float(e @ pb)

    saturated = abs(u_cmd) > case.u_max + 1.0e-12
    adaptation_gate = 0.0 if case.freeze_when_saturated and saturated else 1.0
    theta_dot = -adaptation_gate * GAMMA * phi * sigma

    x_dot = np.array([
        x[1],
        A_P1 * x[0] + A_P2 * x[1] + B * u_act
    ])
    xm_dot = np.array([
        xm[1],
        A_M1 * xm[0] + A_M2 * xm[1] + B_M * r
    ])

    diagnostics = {
        "r": r,
        "u_cmd": u_cmd,
        "u_act": u_act,
        "e1": e[0],
        "e2": e[1],
        "theta_norm": float(np.linalg.norm(theta)),
        "saturated": float(saturated)
    }
    return np.concatenate((x_dot, xm_dot, theta_dot)), diagnostics


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


def simulate(case: Case) -> dict[str, np.ndarray]:
    times = np.arange(0.0, T_FINAL + 0.5 * DT, DT)
    z = np.array([*case.x0, *case.xm0, 0.0, 0.0, 0.0], dtype=float)

    rows: dict[str, list[float]] = {
        key: [] for key in [
            "t", "x1", "x2", "xm1", "xm2",
            "theta1", "theta2", "theta_r",
            "r", "u_cmd", "u_act", "e1", "e2",
            "theta_norm", "saturated"
        ]
    }

    for t in times:
        _, d = derivative(t, z, case)
        values = {
            "t": t,
            "x1": z[0],
            "x2": z[1],
            "xm1": z[2],
            "xm2": z[3],
            "theta1": z[4],
            "theta2": z[5],
            "theta_r": z[6],
            **d
        }
        for key, value in values.items():
            rows[key].append(float(value))
        z = rk4_step(t, z, DT, case)

    return {key: np.asarray(value) for key, value in rows.items()}


def compute_metrics(data: dict[str, np.ndarray]) -> dict[str, float]:
    error_norm = np.sqrt(data["e1"] ** 2 + data["e2"] ** 2)
    tail = data["t"] >= 10.0
    return {
        "peak_error": float(np.max(error_norm)),
        "tail_rms_error": float(np.sqrt(np.mean(error_norm[tail] ** 2))),
        "saturation_fraction": float(np.mean(data["saturated"])),
        "maximum_theta_norm": float(np.max(data["theta_norm"]))
    }


def save_csv(results: dict[str, dict[str, np.ndarray]]) -> None:
    OUTPUT_DIR.mkdir(exist_ok=True)
    path = OUTPUT_DIR / "Chapter6_Lesson5_results.csv"
    fields = ["case"] + list(next(iter(results.values())).keys())
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.DictWriter(stream, fieldnames=fields)
        writer.writeheader()
        for case_name, data in results.items():
            for i in range(len(data["t"])):
                row = {"case": case_name}
                row.update({key: float(values[i]) for key, values in data.items()})
                writer.writerow(row)


def save_plots(results: dict[str, dict[str, np.ndarray]]) -> None:
    OUTPUT_DIR.mkdir(exist_ok=True)

    plt.figure(figsize=(10, 6))
    for name, data in results.items():
        plt.plot(data["t"], data["x1"], label=f"{name}: x1")
    nominal = results["matched_nominal"]
    plt.plot(nominal["t"], nominal["xm1"], "k--", linewidth=2.0, label="reference model xm1")
    plt.xlabel("Time (s)")
    plt.ylabel("Position state")
    plt.grid(True)
    plt.legend(fontsize=8)
    plt.tight_layout()
    plt.savefig(OUTPUT_DIR / "Chapter6_Lesson5_tracking.png", dpi=180)
    plt.close()

    plt.figure(figsize=(10, 6))
    for name, data in results.items():
        plt.plot(data["t"], data["u_cmd"], label=f"{name}: commanded")
        if "saturated" in name:
            plt.plot(data["t"], data["u_act"], "--", label=f"{name}: applied")
    plt.xlabel("Time (s)")
    plt.ylabel("Control input")
    plt.grid(True)
    plt.legend(fontsize=8)
    plt.tight_layout()
    plt.savefig(OUTPUT_DIR / "Chapter6_Lesson5_control.png", dpi=180)
    plt.close()

    plt.figure(figsize=(10, 6))
    for name, data in results.items():
        plt.plot(data["t"], data["theta_norm"], label=name)
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter-vector norm")
    plt.grid(True)
    plt.legend(fontsize=8)
    plt.tight_layout()
    plt.savefig(OUTPUT_DIR / "Chapter6_Lesson5_parameters.png", dpi=180)
    plt.close()


def main() -> None:
    cases = [
        Case("matched_nominal", (0.0, 0.0), (0.0, 0.0), 50.0, False),
        Case("mismatch_nominal", (1.5, -0.5), (0.0, 0.0), 50.0, False),
        Case("mismatch_saturated", (1.5, -0.5), (0.0, 0.0), 1.0, False),
        Case("mismatch_saturated_freeze", (1.5, -0.5), (0.0, 0.0), 1.0, True)
    ]

    results = {case.name: simulate(case) for case in cases}
    save_csv(results)
    save_plots(results)

    print("Case metrics")
    print("-" * 88)
    for case in cases:
        metrics = compute_metrics(results[case.name])
        print(
            f"{case.name:30s} "
            f"peak_error={metrics['peak_error']:.4f}  "
            f"tail_rms={metrics['tail_rms_error']:.4f}  "
            f"sat_fraction={metrics['saturation_fraction']:.3f}  "
            f"max_theta_norm={metrics['maximum_theta_norm']:.4f}"
        )


if __name__ == "__main__":
    main()

      

The C++17 version uses only the standard library and writes all trajectories to CSV for later visualization.

Chapter6_Lesson5.cpp

// Chapter6_Lesson5.cpp
// Second-order MRAC: initial-condition transients and actuator saturation.
// Build: g++ -std=c++17 -O2 Chapter6_Lesson5.cpp -o Chapter6_Lesson5

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

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

struct Case {
    std::string name;
    std::array<double, 2> x0;
    std::array<double, 2> xm0;
    double uMax;
    bool freezeWhenSaturated;
};

struct Sample {
    double t;
    State z;
    double r;
    double uCmd;
    double uAct;
    double e1;
    double e2;
    double thetaNorm;
    bool saturated;
};

constexpr double AP1 = -1.0;
constexpr double AP2 = -0.8;
constexpr double B = 1.2;
constexpr double AM1 = -4.0;
constexpr double AM2 = -3.2;
constexpr double BM = 4.0;
constexpr double DT = 0.002;
constexpr double TFINAL = 15.0;

constexpr std::array<std::array<double, 2>, 2> P{ {
    { {1.18125, 0.125} },
    { {0.125, 0.1953125} }
} };
constexpr std::array<double, 3> GAMMA{ {4.0, 2.0, 3.0} };

double reference(double t) {
    return t < 0.5 ? 0.0 : 1.5;
}

double saturate(double value, double limit) {
    return std::clamp(value, -limit, limit);
}

State addScaled(const State& a, const State& b, double scale) {
    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = a[i] + scale * b[i];
    }
    return result;
}

State derivative(double t, const State& z, const Case& c, Sample* diagnostic = nullptr) {
    const double r = reference(t);
    const double uCmd = z[4] * z[0] + z[5] * z[1] + z[6] * r;
    const double uAct = saturate(uCmd, c.uMax);

    const double e1 = z[0] - z[2];
    const double e2 = z[1] - z[3];

    const double pb1 = P[0][1] * B;
    const double pb2 = P[1][1] * B;
    const double sigma = e1 * pb1 + e2 * pb2;

    const bool saturated = std::abs(uCmd) > c.uMax + 1.0e-12;
    const double gate = (c.freezeWhenSaturated && saturated) ? 0.0 : 1.0;

    State dz{};
    dz[0] = z[1];
    dz[1] = AP1 * z[0] + AP2 * z[1] + B * uAct;
    dz[2] = z[3];
    dz[3] = AM1 * z[2] + AM2 * z[3] + BM * r;
    dz[4] = -gate * GAMMA[0] * z[0] * sigma;
    dz[5] = -gate * GAMMA[1] * z[1] * sigma;
    dz[6] = -gate * GAMMA[2] * r * sigma;

    if (diagnostic != nullptr) {
        diagnostic->t = t;
        diagnostic->z = z;
        diagnostic->r = r;
        diagnostic->uCmd = uCmd;
        diagnostic->uAct = uAct;
        diagnostic->e1 = e1;
        diagnostic->e2 = e2;
        diagnostic->thetaNorm =
            std::sqrt(z[4] * z[4] + z[5] * z[5] + z[6] * z[6]);
        diagnostic->saturated = saturated;
    }
    return dz;
}

State rk4Step(double t, const State& z, const Case& c) {
    const State k1 = derivative(t, z, c);
    const State k2 = derivative(t + 0.5 * DT, addScaled(z, k1, 0.5 * DT), c);
    const State k3 = derivative(t + 0.5 * DT, addScaled(z, k2, 0.5 * DT), c);
    const State k4 = derivative(t + DT, addScaled(z, k3, DT), c);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (DT / 6.0) *
            (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return next;
}

std::vector<Sample> simulate(const Case& c) {
    State z{ {c.x0[0], c.x0[1], c.xm0[0], c.xm0[1], 0.0, 0.0, 0.0} };
    const int steps = static_cast<int>(std::round(TFINAL / DT));
    std::vector<Sample> samples;
    samples.reserve(static_cast<std::size_t>(steps + 1));

    for (int k = 0; k <= steps; ++k) {
        const double t = k * DT;
        Sample sample{};
        derivative(t, z, c, &sample);
        samples.push_back(sample);
        z = rk4Step(t, z, c);
    }
    return samples;
}

void writeCsv(const std::vector<Case>& cases,
              const std::vector<std::vector<Sample>>& results) {
    std::ofstream file("Chapter6_Lesson5_results_cpp.csv");
    if (!file) {
        throw std::runtime_error("Cannot create CSV output.");
    }
    file << "case,t,x1,x2,xm1,xm2,theta1,theta2,theta_r,r,u_cmd,u_act,e1,e2,"
            "theta_norm,saturated\n";
    file << std::setprecision(12);

    for (std::size_t c = 0; c < cases.size(); ++c) {
        for (const Sample& s : results[c]) {
            file << cases[c].name << ',' << s.t;
            for (double value : s.z) {
                file << ',' << value;
            }
            file << ',' << s.r << ',' << s.uCmd << ',' << s.uAct
                 << ',' << s.e1 << ',' << s.e2 << ',' << s.thetaNorm
                 << ',' << (s.saturated ? 1 : 0) << '\n';
        }
    }
}

void printMetrics(const Case& c, const std::vector<Sample>& samples) {
    double peakError = 0.0;
    double tailSquareSum = 0.0;
    std::size_t tailCount = 0;
    std::size_t saturationCount = 0;
    double maxThetaNorm = 0.0;

    for (const Sample& s : samples) {
        const double errorNorm = std::hypot(s.e1, s.e2);
        peakError = std::max(peakError, errorNorm);
        maxThetaNorm = std::max(maxThetaNorm, s.thetaNorm);
        saturationCount += s.saturated ? 1U : 0U;
        if (s.t >= 10.0) {
            tailSquareSum += errorNorm * errorNorm;
            ++tailCount;
        }
    }

    const double tailRms = std::sqrt(tailSquareSum / static_cast<double>(tailCount));
    const double satFraction =
        static_cast<double>(saturationCount) / static_cast<double>(samples.size());

    std::cout << std::left << std::setw(30) << c.name
              << " peak_error=" << std::setw(9) << peakError
              << " tail_rms=" << std::setw(9) << tailRms
              << " sat_fraction=" << std::setw(9) << satFraction
              << " max_theta_norm=" << maxThetaNorm << '\n';
}

int main() {
    try {
        const std::vector<Case> cases{
            {"matched_nominal", {0.0, 0.0}, {0.0, 0.0}, 50.0, false},
            {"mismatch_nominal", {1.5, -0.5}, {0.0, 0.0}, 50.0, false},
            {"mismatch_saturated", {1.5, -0.5}, {0.0, 0.0}, 1.0, false},
            {"mismatch_saturated_freeze", {1.5, -0.5}, {0.0, 0.0}, 1.0, true}
        };

        std::vector<std::vector<Sample>> results;
        for (const Case& c : cases) {
            results.push_back(simulate(c));
        }

        writeCsv(cases, results);
        std::cout << std::fixed << std::setprecision(5);
        for (std::size_t i = 0; i < cases.size(); ++i) {
            printMetrics(cases[i], results[i]);
        }
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

      

The Java version uses records for immutable case and sample data and requires a modern JDK with record support.

Chapter6_Lesson5.java

// Chapter6_Lesson5.java
// Second-order MRAC: initial-condition transients and actuator saturation.
// Build: javac Chapter6_Lesson5.java
// Run:   java Chapter6_Lesson5

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;

public final class Chapter6_Lesson5 {
    private static final double AP1 = -1.0;
    private static final double AP2 = -0.8;
    private static final double B = 1.2;
    private static final double AM1 = -4.0;
    private static final double AM2 = -3.2;
    private static final double BM = 4.0;
    private static final double DT = 0.002;
    private static final double T_FINAL = 15.0;

    private static final double[][] P = {
        {1.18125, 0.125},
        {0.125, 0.1953125}
    };
    private static final double[] GAMMA = {4.0, 2.0, 3.0};

    private record CaseConfig(
        String name,
        double[] x0,
        double[] xm0,
        double uMax,
        boolean freezeWhenSaturated
    ) {}

    private record Sample(
        double t,
        double[] z,
        double r,
        double uCmd,
        double uAct,
        double e1,
        double e2,
        double thetaNorm,
        boolean saturated
    ) {}

    private static double reference(double t) {
        return t < 0.5 ? 0.0 : 1.5;
    }

    private static double saturate(double value, double limit) {
        return Math.max(-limit, Math.min(limit, value));
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] + scale * b[i];
        }
        return result;
    }

    private static double[] derivative(
        double t,
        double[] z,
        CaseConfig config,
        List<Sample> diagnostics
    ) {
        double r = reference(t);
        double uCmd = z[4] * z[0] + z[5] * z[1] + z[6] * r;
        double uAct = saturate(uCmd, config.uMax());

        double e1 = z[0] - z[2];
        double e2 = z[1] - z[3];

        double pb1 = P[0][1] * B;
        double pb2 = P[1][1] * B;
        double sigma = e1 * pb1 + e2 * pb2;

        boolean saturated = Math.abs(uCmd) > config.uMax() + 1.0e-12;
        double gate = config.freezeWhenSaturated() && saturated ? 0.0 : 1.0;

        double[] dz = new double[7];
        dz[0] = z[1];
        dz[1] = AP1 * z[0] + AP2 * z[1] + B * uAct;
        dz[2] = z[3];
        dz[3] = AM1 * z[2] + AM2 * z[3] + BM * r;
        dz[4] = -gate * GAMMA[0] * z[0] * sigma;
        dz[5] = -gate * GAMMA[1] * z[1] * sigma;
        dz[6] = -gate * GAMMA[2] * r * sigma;

        if (diagnostics != null) {
            double thetaNorm = Math.sqrt(
                z[4] * z[4] + z[5] * z[5] + z[6] * z[6]
            );
            diagnostics.add(new Sample(
                t, z.clone(), r, uCmd, uAct, e1, e2, thetaNorm, saturated
            ));
        }
        return dz;
    }

    private static double[] rk4Step(double t, double[] z, CaseConfig config) {
        double[] k1 = derivative(t, z, config, null);
        double[] k2 = derivative(
            t + 0.5 * DT, addScaled(z, k1, 0.5 * DT), config, null
        );
        double[] k3 = derivative(
            t + 0.5 * DT, addScaled(z, k2, 0.5 * DT), config, null
        );
        double[] k4 = derivative(
            t + DT, addScaled(z, k3, DT), config, null
        );

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 6.0) *
                (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        return next;
    }

    private static List<Sample> simulate(CaseConfig config) {
        double[] z = {
            config.x0()[0], config.x0()[1],
            config.xm0()[0], config.xm0()[1],
            0.0, 0.0, 0.0
        };
        int steps = (int) Math.round(T_FINAL / DT);
        List<Sample> samples = new ArrayList<>(steps + 1);

        for (int k = 0; k <= steps; k++) {
            double t = k * DT;
            derivative(t, z, config, samples);
            z = rk4Step(t, z, config);
        }
        return samples;
    }

    private static void writeCsv(
        List<CaseConfig> cases,
        List<List<Sample>> results
    ) throws IOException {
        Path output = Path.of("Chapter6_Lesson5_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output)) {
            writer.write(
                "case,t,x1,x2,xm1,xm2,theta1,theta2,theta_r,r,u_cmd,u_act,"
                + "e1,e2,theta_norm,saturated\n"
            );
            for (int c = 0; c < cases.size(); c++) {
                for (Sample s : results.get(c)) {
                    writer.write(cases.get(c).name());
                    writer.write(String.format(
                        Locale.US,
                        ",%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,"
                        + "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%d%n",
                        s.t(),
                        s.z()[0], s.z()[1], s.z()[2], s.z()[3],
                        s.z()[4], s.z()[5], s.z()[6],
                        s.r(), s.uCmd(), s.uAct(), s.e1(), s.e2(),
                        s.thetaNorm(), s.saturated() ? 1 : 0
                    ));
                }
            }
        }
    }

    private static void printMetrics(CaseConfig config, List<Sample> samples) {
        double peakError = 0.0;
        double tailSquareSum = 0.0;
        int tailCount = 0;
        int saturationCount = 0;
        double maxThetaNorm = 0.0;

        for (Sample s : samples) {
            double errorNorm = Math.hypot(s.e1(), s.e2());
            peakError = Math.max(peakError, errorNorm);
            maxThetaNorm = Math.max(maxThetaNorm, s.thetaNorm());
            if (s.saturated()) {
                saturationCount++;
            }
            if (s.t() >= 10.0) {
                tailSquareSum += errorNorm * errorNorm;
                tailCount++;
            }
        }

        double tailRms = Math.sqrt(tailSquareSum / tailCount);
        double satFraction = (double) saturationCount / samples.size();

        System.out.printf(
            Locale.US,
            "%-30s peak_error=%8.5f tail_rms=%8.5f "
            + "sat_fraction=%8.5f max_theta_norm=%8.5f%n",
            config.name(), peakError, tailRms, satFraction, maxThetaNorm
        );
    }

    public static void main(String[] args) {
        List<CaseConfig> cases = List.of(
            new CaseConfig(
                "matched_nominal",
                new double[]{0.0, 0.0},
                new double[]{0.0, 0.0},
                50.0,
                false
            ),
            new CaseConfig(
                "mismatch_nominal",
                new double[]{1.5, -0.5},
                new double[]{0.0, 0.0},
                50.0,
                false
            ),
            new CaseConfig(
                "mismatch_saturated",
                new double[]{1.5, -0.5},
                new double[]{0.0, 0.0},
                1.0,
                false
            ),
            new CaseConfig(
                "mismatch_saturated_freeze",
                new double[]{1.5, -0.5},
                new double[]{0.0, 0.0},
                1.0,
                true
            )
        );

        List<List<Sample>> results = new ArrayList<>();
        for (CaseConfig config : cases) {
            results.add(simulate(config));
        }

        try {
            writeCsv(cases, results);
            for (int i = 0; i < cases.size(); i++) {
                printMetrics(cases.get(i), results.get(i));
            }
        } catch (IOException ex) {
            System.err.println("I/O error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

      

The MATLAB script reproduces all four cases, computes diagnostic metrics, saves a CSV file, and creates tracking, control, and parameter plots.

Chapter6_Lesson5.m

% Chapter6_Lesson5.m
% Second-order MRAC: initial-condition transients and actuator saturation.
% Requires MATLAB R2016b or newer for local functions in scripts.

clear; clc; close all;

cases = struct( ...
    'name', { ...
        'matched_nominal', ...
        'mismatch_nominal', ...
        'mismatch_saturated', ...
        'mismatch_saturated_freeze'}, ...
    'x0', { ...
        [0; 0], ...
        [1.5; -0.5], ...
        [1.5; -0.5], ...
        [1.5; -0.5]}, ...
    'xm0', { ...
        [0; 0], ...
        [0; 0], ...
        [0; 0], ...
        [0; 0]}, ...
    'uMax', {50, 50, 1, 1}, ...
    'freezeWhenSaturated', {false, false, false, true});

dt = 0.002;
tFinal = 15;
results = cell(size(cases));
summaryRows = cell(numel(cases), 5);

for k = 1:numel(cases)
    results{k} = simulateCase(cases(k), dt, tFinal);
    d = results{k};
    errorNorm = hypot(d.e1, d.e2);
    tail = d.t >= 10;

    summaryRows(k, :) = { ...
        cases(k).name, ...
        max(errorNorm), ...
        sqrt(mean(errorNorm(tail).^2)), ...
        mean(d.saturated), ...
        max(d.thetaNorm)};
end

summaryTable = cell2table(summaryRows, ...
    'VariableNames', {'Case', 'PeakError', 'TailRMSError', ...
                      'SaturationFraction', 'MaximumThetaNorm'});
disp(summaryTable);

allData = table();
for k = 1:numel(cases)
    d = results{k};
    caseName = repmat(string(cases(k).name), numel(d.t), 1);
    current = table(caseName, d.t, d.x1, d.x2, d.xm1, d.xm2, ...
        d.theta1, d.theta2, d.thetaR, d.r, d.uCmd, d.uAct, ...
        d.e1, d.e2, d.thetaNorm, d.saturated, ...
        'VariableNames', {'Case', 't', 'x1', 'x2', 'xm1', 'xm2', ...
        'theta1', 'theta2', 'thetaR', 'r', 'uCmd', 'uAct', ...
        'e1', 'e2', 'thetaNorm', 'saturated'});
    allData = [allData; current]; %#ok<AGROW>
end
writetable(allData, 'Chapter6_Lesson5_results_matlab.csv');

figure;
hold on;
for k = 1:numel(cases)
    plot(results{k}.t, results{k}.x1, ...
        'DisplayName', cases(k).name);
end
plot(results{1}.t, results{1}.xm1, 'k--', 'LineWidth', 1.8, ...
    'DisplayName', 'reference model xm1');
xlabel('Time (s)');
ylabel('Position state');
grid on;
legend('Location', 'best');
title('Initial-condition and saturation comparison');

figure;
hold on;
for k = 1:numel(cases)
    plot(results{k}.t, results{k}.uCmd, ...
        'DisplayName', [cases(k).name ' commanded']);
    if contains(cases(k).name, 'saturated')
        plot(results{k}.t, results{k}.uAct, '--', ...
            'DisplayName', [cases(k).name ' applied']);
    end
end
xlabel('Time (s)');
ylabel('Control input');
grid on;
legend('Location', 'best');
title('Commanded and applied control');

figure;
hold on;
for k = 1:numel(cases)
    plot(results{k}.t, results{k}.thetaNorm, ...
        'DisplayName', cases(k).name);
end
xlabel('Time (s)');
ylabel('Parameter-vector norm');
grid on;
legend('Location', 'best');
title('Adaptive-parameter transient');

function data = simulateCase(caseData, dt, tFinal)
    t = (0:dt:tFinal).';
    n = numel(t);
    z = zeros(n, 7);
    z(1, :) = [caseData.x0(:).', caseData.xm0(:).', 0, 0, 0];

    r = zeros(n, 1);
    uCmd = zeros(n, 1);
    uAct = zeros(n, 1);
    e1 = zeros(n, 1);
    e2 = zeros(n, 1);
    thetaNorm = zeros(n, 1);
    saturated = false(n, 1);

    for i = 1:n
        [~, diag] = rhs(t(i), z(i, :).', caseData);
        r(i) = diag.r;
        uCmd(i) = diag.uCmd;
        uAct(i) = diag.uAct;
        e1(i) = diag.e1;
        e2(i) = diag.e2;
        thetaNorm(i) = diag.thetaNorm;
        saturated(i) = diag.saturated;

        if i < n
            zi = z(i, :).';
            k1 = rhs(t(i), zi, caseData);
            k2 = rhs(t(i) + 0.5 * dt, zi + 0.5 * dt * k1, caseData);
            k3 = rhs(t(i) + 0.5 * dt, zi + 0.5 * dt * k2, caseData);
            k4 = rhs(t(i) + dt, zi + dt * k3, caseData);
            z(i + 1, :) = (zi + (dt / 6) * ...
                (k1 + 2 * k2 + 2 * k3 + k4)).';
        end
    end

    data = struct( ...
        't', t, ...
        'x1', z(:, 1), ...
        'x2', z(:, 2), ...
        'xm1', z(:, 3), ...
        'xm2', z(:, 4), ...
        'theta1', z(:, 5), ...
        'theta2', z(:, 6), ...
        'thetaR', z(:, 7), ...
        'r', r, ...
        'uCmd', uCmd, ...
        'uAct', uAct, ...
        'e1', e1, ...
        'e2', e2, ...
        'thetaNorm', thetaNorm, ...
        'saturated', saturated);
end

function [dz, diagnostic] = rhs(t, z, caseData)
    ap1 = -1.0;
    ap2 = -0.8;
    b = 1.2;

    am1 = -4.0;
    am2 = -3.2;
    bm = 4.0;

    P = [1.18125, 0.125; 0.125, 0.1953125];
    Gamma = [4.0; 2.0; 3.0];

    if t < 0.5
        reference = 0.0;
    else
        reference = 1.5;
    end

    x = z(1:2);
    xm = z(3:4);
    theta = z(5:7);
    phi = [x(1); x(2); reference];

    commanded = theta.' * phi;
    applied = min(caseData.uMax, max(-caseData.uMax, commanded));

    error = x - xm;
    sigma = error.' * P * [0; b];

    isSaturated = abs(commanded) > caseData.uMax + 1e-12;
    gate = 1.0;
    if caseData.freezeWhenSaturated && isSaturated
        gate = 0.0;
    end

    xDot = [ ...
        x(2); ...
        ap1 * x(1) + ap2 * x(2) + b * applied];

    xmDot = [ ...
        xm(2); ...
        am1 * xm(1) + am2 * xm(2) + bm * reference];

    thetaDot = -gate * Gamma .* phi * sigma;
    dz = [xDot; xmDot; thetaDot];

    diagnostic = struct( ...
        'r', reference, ...
        'uCmd', commanded, ...
        'uAct', applied, ...
        'e1', error(1), ...
        'e2', error(2), ...
        'thetaNorm', norm(theta), ...
        'saturated', isSaturated);
end

      

The next script programmatically creates a Simulink model. It uses the accompanying Level-2 MATLAB S-function as the continuous adaptive plant and controller block.

Chapter6_Lesson5_Simulink.m

% Chapter6_Lesson5_Simulink.m
% Programmatically builds and runs a Simulink model that uses the
% Level-2 MATLAB S-function in Chapter6_Lesson5_SFunction.m.
% Requires Simulink.

clear; clc;

model = 'Chapter6_Lesson5_MRAC';
if bdIsLoaded(model)
    close_system(model, 0);
end
if exist([model '.slx'], 'file')
    delete([model '.slx']);
end

new_system(model);
open_system(model);

add_block('simulink/Sources/Step', [model '/Reference'], ...
    'Time', '0.5', ...
    'Before', '0', ...
    'After', '1.5', ...
    'Position', [40 95 80 125]);

add_block('simulink/User-Defined Functions/S-Function', ...
    [model '/SecondOrderMRAC'], ...
    'FunctionName', 'Chapter6_Lesson5_SFunction', ...
    'Position', [145 70 295 150]);

add_block('simulink/Sinks/Scope', [model '/Scope'], ...
    'NumInputPorts', '1', ...
    'Position', [380 65 430 115]);

add_block('simulink/Sinks/To Workspace', [model '/LoggedSignals'], ...
    'VariableName', 'mracSignals', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [360 135 455 165]);

add_line(model, 'Reference/1', 'SecondOrderMRAC/1');
add_line(model, 'SecondOrderMRAC/1', 'Scope/1');
add_line(model, 'SecondOrderMRAC/1', 'LoggedSignals/1');

set_param(model, ...
    'StopTime', '15', ...
    'Solver', 'ode45', ...
    'MaxStep', '0.002');

save_system(model);
simulationOutput = sim(model);

signals = simulationOutput.mracSignals;
time = signals.time;
values = signals.signals.values;

figure;
plot(time, values(:, 1), time, values(:, 2), '--');
grid on;
xlabel('Time (s)');
ylabel('Position');
legend('x1', 'xm1');
title('Simulink second-order MRAC');

figure;
plot(time, values(:, 5), time, values(:, 6), '--');
grid on;
xlabel('Time (s)');
ylabel('Control input');
legend('u commanded', 'u applied');
title('Actuator saturation');

disp('Signal order:');
disp('[x1, xm1, x2, xm2, uCmd, uAct, theta1, theta2, thetaR]');

      

Chapter6_Lesson5_SFunction.m

function Chapter6_Lesson5_SFunction(block)
% Chapter6_Lesson5_SFunction.m
% Level-2 MATLAB S-function for the Chapter 6, Lesson 5 Simulink model.
% Continuous states:
% [x1, x2, xm1, xm2, theta1, theta2, thetaR]

setup(block);

function setup(block)
    block.NumDialogPrms = 0;
    block.NumInputPorts = 1;
    block.NumOutputPorts = 1;

    block.SetPreCompInpPortInfoToDynamic;
    block.SetPreCompOutPortInfoToDynamic;

    block.InputPort(1).Dimensions = 1;
    block.InputPort(1).DirectFeedthrough = true;

    block.OutputPort(1).Dimensions = 9;

    block.NumContStates = 7;
    block.SampleTimes = [0 0];
    block.SimStateCompliance = 'DefaultSimState';

    block.RegBlockMethod('InitializeConditions', @initializeConditions);
    block.RegBlockMethod('Outputs', @outputs);
    block.RegBlockMethod('Derivatives', @derivatives);
end

function initializeConditions(block)
    block.ContStates.Data = [1.5; -0.5; 0; 0; 0; 0; 0];
end

function outputs(block)
    z = block.ContStates.Data;
    r = block.InputPort(1).Data;

    uMax = 1.0;
    phi = [z(1); z(2); r];
    uCmd = z(5:7).' * phi;
    uAct = min(uMax, max(-uMax, uCmd));

    block.OutputPort(1).Data = [ ...
        z(1); z(3); z(2); z(4); ...
        uCmd; uAct; z(5); z(6); z(7)];
end

function derivatives(block)
    ap1 = -1.0;
    ap2 = -0.8;
    b = 1.2;

    am1 = -4.0;
    am2 = -3.2;
    bm = 4.0;

    P = [1.18125, 0.125; 0.125, 0.1953125];
    Gamma = [4.0; 2.0; 3.0];
    uMax = 1.0;

    z = block.ContStates.Data;
    r = block.InputPort(1).Data;

    x = z(1:2);
    xm = z(3:4);
    theta = z(5:7);
    phi = [x(1); x(2); r];

    uCmd = theta.' * phi;
    uAct = min(uMax, max(-uMax, uCmd));

    error = x - xm;
    sigma = error.' * P * [0; b];

    xDot = [ ...
        x(2); ...
        ap1 * x(1) + ap2 * x(2) + b * uAct];

    xmDot = [ ...
        xm(2); ...
        am1 * xm(1) + am2 * xm(2) + bm * r];

    % Plain adaptation is retained here so that saturation-induced
    % parameter growth can be observed directly.
    thetaDot = -Gamma .* phi * sigma;

    block.Derivatives.Data = [xDot; xmDot; thetaDot];
end
end

      

The notebook uses a from-scratch RK4 integrator, exports the combined result table, and generates three diagnostic plots.

Chapter6_Lesson5.nb

(* Chapter6_Lesson5.nb
   Wolfram Language implementation of second-order MRAC with
   initial-condition mismatch and actuator saturation. *)

ClearAll["Global`*"];

ap1 = -1.0; ap2 = -0.8; b = 1.2;
am1 = -4.0; am2 = -3.2; bm = 4.0;
p = { {1.18125, 0.125}, {0.125, 0.1953125} };
gamma = {4.0, 2.0, 3.0};
dt = 0.002; tFinal = 15.0;

reference[t_] := Piecewise[{ {0.0, t < 0.5} }, 1.5];
sat[value_, limit_] := Clip[value, {-limit, limit}];

rk4Step[f_, t_, z_, h_] := Module[{k1, k2, k3, k4},
  k1 = f[t, z];
  k2 = f[t + h/2, z + h k1/2];
  k3 = f[t + h/2, z + h k2/2];
  k4 = f[t + h, z + h k3];
  z + h (k1 + 2 k2 + 2 k3 + k4)/6
];

simulateCase[name_, x0_, xm0_, uMax_, freeze_] := Module[
  {rhs, times, states, diagnostics, z0},

  rhs[t_, z_] := Module[
    {x, xm, theta, r, phi, uCmd, uAct, error, sigma,
     saturated, gate, xDot, xmDot, thetaDot},

    x = z[[1 ;; 2]];
    xm = z[[3 ;; 4]];
    theta = z[[5 ;; 7]];
    r = reference[t];
    phi = {x[[1]], x[[2]], r};

    uCmd = theta.phi;
    uAct = sat[uCmd, uMax];
    error = x - xm;
    sigma = error.p.{0.0, b};

    saturated = Abs[uCmd] > uMax + 10^-12;
    gate = If[TrueQ[freeze] && saturated, 0.0, 1.0];

    xDot = {x[[2]], ap1 x[[1]] + ap2 x[[2]] + b uAct};
    xmDot = {xm[[2]], am1 xm[[1]] + am2 xm[[2]] + bm r};
    thetaDot = -gate gamma phi sigma;

    Join[xDot, xmDot, thetaDot]
  ];

  times = Range[0.0, tFinal, dt];
  z0 = Join[x0, xm0, {0.0, 0.0, 0.0}];
  states = Rest@FoldList[
      rk4Step[rhs, #2[[1]], #1, dt] &,
      z0,
      Transpose[{Most[times], Rest[times]}]
    ];

  states = Prepend[states, z0];

  diagnostics = MapThread[
    Function[{t, z},
      Module[{x, xm, theta, r, phi, uCmd, uAct, error, saturated},
        x = z[[1 ;; 2]];
        xm = z[[3 ;; 4]];
        theta = z[[5 ;; 7]];
        r = reference[t];
        phi = {x[[1]], x[[2]], r};
        uCmd = theta.phi;
        uAct = sat[uCmd, uMax];
        error = x - xm;
        saturated = Boole[Abs[uCmd] > uMax + 10^-12];
        Join[
          {name, t},
          z,
          {r, uCmd, uAct, error[[1]], error[[2]],
           Norm[theta], saturated}
        ]
      ]
    ],
    {times, states}
  ];

  diagnostics
];

cases = {
  {"matched_nominal", {0.0, 0.0}, {0.0, 0.0}, 50.0, False},
  {"mismatch_nominal", {1.5, -0.5}, {0.0, 0.0}, 50.0, False},
  {"mismatch_saturated", {1.5, -0.5}, {0.0, 0.0}, 1.0, False},
  {"mismatch_saturated_freeze", {1.5, -0.5}, {0.0, 0.0}, 1.0, True}
};

results = simulateCase @@@ cases;
headers = {
  "case", "t", "x1", "x2", "xm1", "xm2",
  "theta1", "theta2", "thetaR", "r", "uCmd", "uAct",
  "e1", "e2", "thetaNorm", "saturated"
};

Export[
  "Chapter6_Lesson5_results_mathematica.csv",
  Prepend[Join @@ results, headers]
];

trackingPlot = ListLinePlot[
  Table[
    results[[k, All, {2, 3}]],
    {k, Length[results]}
  ],
  PlotLegends -> cases[[All, 1]],
  Frame -> True,
  FrameLabel -> {"Time (s)", "x1"},
  PlotRange -> All,
  ImageSize -> Large
];

controlPlot = ListLinePlot[
  {
    results[[3, All, {2, 11}]],
    results[[3, All, {2, 12}]]
  },
  PlotLegends -> {"u commanded", "u applied"},
  Frame -> True,
  FrameLabel -> {"Time (s)", "Control"},
  PlotRange -> All,
  ImageSize -> Large
];

parameterPlot = ListLinePlot[
  Table[
    results[[k, All, {2, 15}]],
    {k, Length[results]}
  ],
  PlotLegends -> cases[[All, 1]],
  Frame -> True,
  FrameLabel -> {"Time (s)", "Parameter-vector norm"},
  PlotRange -> All,
  ImageSize -> Large
];

Export["Chapter6_Lesson5_tracking_mathematica.png", trackingPlot];
Export["Chapter6_Lesson5_control_mathematica.png", controlPlot];
Export["Chapter6_Lesson5_parameters_mathematica.png", parameterPlot];

{trackingPlot, controlPlot, parameterPlot}

      

1. Reference-model check: simulate the model alone and verify the expected damping and settling behavior.
2. Lyapunov-matrix check: numerically evaluate \(\mathbf{A}_m^\mathsf{T}\mathbf{P}+ \mathbf{P}\mathbf{A}_m+\mathbf{Q}\) and confirm that its norm is near machine precision.
3. Nominal identity check: with a loose actuator limit, compare numerical \(\dot{V}\) with \(-\mathbf{e}^\mathsf{T}\mathbf{Q}\mathbf{e}\).
4. Saturation check: verify that \(\Delta u=u_a-u_c\) is zero outside saturation and nonzero inside it.
5. Step-size check: repeat with half the integration step. Important metrics should change only slightly.
6. Parameter interpretation: do not infer identification success merely because tracking error becomes small.

13.1 What the Four Cases Demonstrate

• Matched nominal: isolates the command and adaptation transient without initial state error.
• Mismatched nominal: adds stored initial tracking-error energy while retaining the nominal actuator assumption.
• Mismatched saturated: exposes the discrepancy between the control signal used by the adaptive law and the signal delivered to the plant.
• Mismatched saturated with freeze: demonstrates that a simple implementation safeguard can reduce parameter growth while not necessarily improving the state response when the actuator remains pinned.

Problem 1 (Initial Lyapunov Bound): Let \(\mathbf{P}=\begin{bmatrix}1.18125&0.125\\0.125&0.1953125\end{bmatrix}\), \(\boldsymbol{\Gamma}=\operatorname{diag}(4,2,3)\), \(\mathbf{e}(0)=[1.5\;-0.5]^\mathsf{T}\), and \(\widetilde{\boldsymbol{\theta} }(0)= [2.5\;2\;-10/3]^\mathsf{T}\). Compute \(V(0)\) and a bound on \(\|\mathbf{e}(t)\|\).

Solution:

\[ \mathbf{e}^\mathsf{T}(0)\mathbf{P}\mathbf{e}(0) = \begin{bmatrix}1.5&-0.5\end{bmatrix} \mathbf{P} \begin{bmatrix}1.5\\-0.5\end{bmatrix} =2.519140625. \]

\[ \widetilde{\boldsymbol{\theta} }^\mathsf{T}(0) \boldsymbol{\Gamma}^{-1} \widetilde{\boldsymbol{\theta} }(0) = \frac{2.5^2}{4}+\frac{2^2}{2} +\frac{(10/3)^2}{3} \approx 7.266204. \]

\[ V(0)\approx 9.785344. \]

Since \(\lambda_{\min}(\mathbf{P})\approx0.1797115\),

\[ \|\mathbf{e}(t)\| \le\sqrt{\frac{9.785344}{0.1797115} } \approx7.38. \]

The bound is conservative because it permits all Lyapunov energy to appear in the state-error term at one instant.

Problem 2 (Saturation Perturbation): Starting from the saturated error model, derive the Lyapunov derivative when the nominal adaptive law is retained.

Solution:

\[ \dot{\mathbf{e} } =\mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^\mathsf{T} \boldsymbol{\phi} +\mathbf{b}\Delta u. \]

Differentiating \(V=\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{e}+ \widetilde{\boldsymbol{\theta} }^\mathsf{T} \boldsymbol{\Gamma}^{-1} \widetilde{\boldsymbol{\theta} }\), using the Lyapunov equation, and cancelling the nominal cross term with the update law gives

\[ \dot{V} =-\mathbf{e}^\mathsf{T}\mathbf{Q}\mathbf{e} +2\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b}\Delta u. \]

The second term prevents the conclusion \(\dot{V}\le0\) unless additional conditions or a redesigned adaptive architecture are supplied.

Problem 3 (Maximum Initial Parameter Rate): Show that the initial update rate satisfies a norm inequality and explain how state scaling affects it.

Solution:

\[ \|\dot{\boldsymbol{\theta} }(0)\| = \|\boldsymbol{\Gamma}\boldsymbol{\phi}(0) \mathbf{e}^\mathsf{T}(0)\mathbf{P}\mathbf{b}\| \le \|\boldsymbol{\Gamma}\|\, \|\boldsymbol{\phi}(0)\|\, \|\mathbf{e}(0)\|\, \|\mathbf{P}\mathbf{b}\|. \]

If one state is expressed in units that produce a numerically much larger component, that component dominates \(\|\boldsymbol{\phi}\|\) and its associated parameter update. Scaling the regressor makes adaptation gains interpretable and reduces artificial numerical imbalance.

Problem 4 (Reference-Model Feasibility): For zero initial reference-model state and a step command of magnitude \(r_0\), compute the initial desired acceleration.

Solution:

\[ \dot{x}_{m2}(0^+) =-\omega_n^2x_{m1}(0)-2\zeta\omega_n x_{m2}(0) +\omega_n^2r_0 =\omega_n^2r_0. \]

Doubling \(\omega_n\) multiplies this initial acceleration demand by four. Thus a reference model can be stable yet physically infeasible.

Problem 5 (Freeze-on-Saturation Interpretation): Suppose the update law is multiplied by \(g(t)\), where \(g(t)=0\) during saturation and \(g(t)=1\) otherwise. Does the nominal cancellation remain exact?

Solution:

No. The parameter part of \(\dot{V}\) cancels only a fraction \(g(t)\) of the state-parameter cross term. The remaining derivative contains

\[ 2\bigl(1-g(t)\bigr) \mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b}\, \widetilde{\boldsymbol{\theta} }^\mathsf{T}\boldsymbol{\phi} +2\mathbf{e}^\mathsf{T}\mathbf{P}\mathbf{b}\Delta u. \]

Therefore freezing is a practical parameter-growth safeguard, not a direct consequence of the nominal Lyapunov proof. A formal constrained adaptive design requires additional analysis.

Initial state and parameter mismatch determine the initial Lyapunov energy and therefore influence the magnitude of adaptive transients. The reference model contributes its own transient and must be selected with actuator feasibility in mind. When the actuator saturates, the applied input differs from the commanded adaptive input; the error dynamics acquire the term \(\mathbf{b}\Delta u\), and the nominal Lyapunov derivative acquires an indefinite cross term. Consequently, unconstrained MRAC proofs cannot be transferred unchanged to saturated hardware. Reliable implementation requires explicit logging, achievable commands, sensible initialization, scaled regressors, numerical convergence tests, and separate matched, mismatched, and saturated experiments.