After completing this lesson, students should be able to:

• compute the ideal feedback and feedforward parameters from the canonical-form matching equations;
• construct the Lyapunov equation, filtered tracking-error signal, and direct adaptive law for an arbitrary order;
• distinguish guaranteed tracking convergence from the stronger and excitation-dependent claim of parameter convergence;
• diagnose structural mismatch caused by incorrect order, actuator dynamics, unavailable states, saturation, and an unknown control direction; and
• implement and numerically test the controller without relying on a specialized adaptive-control library.

The lesson remains within the plain continuous-time, full-state, known-control-direction MRAC framework of this chapter. Robust modifications are only motivated here; their design is developed in Chapter 8.

Consider an n-th order controllable SISO plant whose complete state is available:

\[ \dot{\mathbf{x} } = \mathbf{A}_p\mathbf{x}+\mathbf{b}_p u, \qquad y=\mathbf{c}^{T}\mathbf{x}, \qquad \mathbf{b}_p=b\mathbf{e}_n, \qquad b\neq 0. \]

In controllable companion form, \( \mathbf{e}_n=[0\;\cdots\;0\;1]^T \) and

\[ \mathbf{A}_p= \begin{bmatrix} 0&1&0&\cdots&0\\ 0&0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1\\ -a_{p0}&-a_{p1}&-a_{p2}&\cdots&-a_{p,n-1} \end{bmatrix}. \]

The desired bounded response is generated by a stable reference model with the same canonical structure:

\[ \dot{\mathbf{x} }_m=\mathbf{A}_m\mathbf{x}_m+ \mathbf{b}_m r, \qquad \mathbf{b}_m=b_m\mathbf{e}_n, \qquad \mathbf{A}_m\text{ Hurwitz}. \]

The direct state-feedback controller is parameterized as

\[ u=\hat{\mathbf{k} }_x^{T}\mathbf{x}+\hat{k}_r r =\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\omega}, \quad \hat{\boldsymbol{\theta} }= \begin{bmatrix}\hat{\mathbf{k} }_x\\\hat{k}_r\end{bmatrix}, \quad \boldsymbol{\omega}= \begin{bmatrix}\mathbf{x}\\r\end{bmatrix}. \]

Exact model matching requires constant ideal parameters satisfying

\[ \mathbf{A}_p+\mathbf{b}_p\mathbf{k}_x^{*T}=\mathbf{A}_m, \qquad \mathbf{b}_p k_r^*=\mathbf{b}_m. \]

Because only the final row differs in the two companion matrices, these equations reduce to

\[ \boxed{ \mathbf{k}_x^*=\frac{\mathbf{a}_p-\mathbf{a}_m}{b}, \qquad k_r^*=\frac{b_m}{b} } \]

where \( \mathbf{a}_p=[a_{p0},\ldots,a_{p,n-1}]^T \) and \( \mathbf{a}_m=[a_{m0},\ldots,a_{m,n-1}]^T \). The controller does not know these ideal values, but their existence is a structural assumption used in the stability proof.

Define tracking and parameter errors by

\[ \mathbf{e}=\mathbf{x}-\mathbf{x}_m, \qquad \tilde{\boldsymbol{\theta} }= \hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^*. \]

Substitution of the matching equations gives the compact error model

\[ \dot{\mathbf{e} } =\mathbf{A}_m\mathbf{e} +b\mathbf{e}_n\tilde{\boldsymbol{\theta} }^T \boldsymbol{\omega}. \]

Select any symmetric matrix \( \mathbf{Q}=\mathbf{Q}^T\succ 0 \). Since \( \mathbf{A}_m \) is Hurwitz, the Lyapunov equation

\[ \mathbf{A}_m^T\mathbf{P}+\mathbf{P}\mathbf{A}_m=-\mathbf{Q} \]

has a unique symmetric solution \( \mathbf{P}=\mathbf{P}^T\succ 0 \). Let \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^T\succ 0 \) be the adaptation-gain matrix and choose

\[ \boxed{ \dot{\hat{\boldsymbol{\theta} } } =-\operatorname{sgn}(b)\boldsymbol{\Gamma} \boldsymbol{\omega}\,\mathbf{e}^T\mathbf{P}\mathbf{e}_n } \]

The scalar \( s_e=\mathbf{e}^T\mathbf{P}\mathbf{e}_n \) is a filtered combination of the state-tracking errors. It is the signal that couples the tracking dynamics to the parameter update.

Theorem. Suppose the plant is exactly represented by the stated canonical model, the sign of \( b \) is known, the matching equations admit constant ideal parameters, the reference input is bounded, and all plant states are measured. Then the adaptive law in Section 3 makes every closed-loop signal bounded and yields asymptotic state tracking.

Proof. Consider

\[ V=\mathbf{e}^T\mathbf{P}\mathbf{e} +|b|\tilde{\boldsymbol{\theta} }^T \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta} }. \]

Since \( \boldsymbol{\theta}^* \) is constant, \( \dot{\tilde{\boldsymbol{\theta} } }= \dot{\hat{\boldsymbol{\theta} } } \). Differentiating gives

\[ \begin{aligned} \dot V &=\mathbf{e}^T (\mathbf{A}_m^T\mathbf{P}+\mathbf{P}\mathbf{A}_m) \mathbf{e}\\ &\quad+2b\,\mathbf{e}^T\mathbf{P}\mathbf{e}_n \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}\\ &\quad+2|b|\tilde{\boldsymbol{\theta} }^T \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \end{aligned} \]

Using the Lyapunov equation and adaptive law,

\[ \begin{aligned} \dot V &=-\mathbf{e}^T\mathbf{Q}\mathbf{e} +2\operatorname{sgn}(b)|b|s_e \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}\\ &\quad-2|b|\operatorname{sgn}(b) \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}s_e\\ &=-\mathbf{e}^T\mathbf{Q}\mathbf{e}\le 0. \end{aligned} \]

Hence \( V(t) \) is nonincreasing, so \( \mathbf{e} \) and \( \tilde{\boldsymbol{\theta} } \) are bounded. A bounded reference passed through a Hurwitz model gives bounded \( \mathbf{x}_m \); therefore \( \mathbf{x}=\mathbf{e}+\mathbf{x}_m \), the regressor, control parameters, and control input are bounded. Also,

\[ \int_0^{\infty}\mathbf{e}^T\mathbf{Q}\mathbf{e}\,dt \le V(0), \]

so \( \mathbf{e}\in\mathcal{L}_2 \). The bounded error dynamics imply bounded \( \dot{\mathbf{e} } \). Barbalat's lemma then gives \( \lim_{t\to\infty}\mathbf{e}(t)=\mathbf{0} \). This proves asymptotic tracking. Notice that the proof establishes bounded parameter estimates, not parameter convergence. ■

flowchart TD
  A["Specify plant order n and measured canonical state x"] --> B["Choose Hurwitz reference-model polynomial"]
  B --> C["Check matching: same order, input in final canonical equation, known sign of b"]
  C --> D["Solve Am^T P + P Am = -Q"]
  D --> E["Form regressor omega = [x; r]"]
  E --> F["Control u = theta_hat^T omega"]
  F --> G["Compute error e = x - xm and scalar se = e^T P en"]
  G --> H["Update theta_hat_dot = -sign(b) Gamma omega se"]
  H --> I["Check tracking, control effort, parameter motion, and assumption violations"]
        

A correct Lyapunov derivation does not remove the need for engineering checks. The final step is essential: the theorem applies to a precise model class, while a simulation or physical plant may contain saturation, delays, hidden states, disturbances, or dynamics omitted from the canonical model.

Consider the uncertain plant

\[ \ddot y+0.8\dot y+1.2y=1.5u, \qquad \mathbf{x}=\begin{bmatrix}y\\\dot y\end{bmatrix}. \]

Choose the critically damped reference model

\[ \ddot y_m+4\dot y_m+4y_m=4r, \]

so that

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

The ideal matching gains are

\[ \begin{aligned} k_{x1}^*&=\frac{1.2-4}{1.5}=-1.866\overline{6},\\ k_{x2}^*&=\frac{0.8-4}{1.5}=-2.133\overline{3},\\ k_r^*&=\frac{4}{1.5}=2.666\overline{6}. \end{aligned} \]

With \( \mathbf{Q}=\mathbf{I}_2 \), let \( \mathbf{P}=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22} \end{bmatrix} \). Expanding the Lyapunov equation gives

\[ \begin{cases} -8p_{12}=-1,\\ p_{11}-4p_{12}-4p_{22}=0,\\ 2p_{12}-8p_{22}=-1, \end{cases} \qquad \mathbf{P}=\begin{bmatrix} \frac{9}{8}&\frac{1}{8}\\ \frac{1}{8}&\frac{5}{32} \end{bmatrix}. \]

The filtered error is therefore

\[ s_e=\mathbf{e}^T\mathbf{P}\mathbf{e}_2 =\frac{1}{8}e_1+\frac{5}{32}e_2. \]

For diagonal \( \boldsymbol{\Gamma}=\operatorname{diag}(20,20,10) \) and positive \( b \), the complete adaptive law is

\[ \begin{aligned} u&=\hat k_{x1}x_1+\hat k_{x2}x_2+\hat k_r r,\\ \dot{\hat k}_{x1}&=-20x_1s_e,\\ \dot{\hat k}_{x2}&=-20x_2s_e,\\ \dot{\hat k}_r&=-10rs_e. \end{aligned} \]

The selected gains are simulation values, not universal tuning rules. The Lyapunov proof remains valid for any positive-definite adaptation gain, but transient behavior and numerical sensitivity change substantially with its magnitude.

Consider

\[ y^{(3)}+1.1\ddot y+2.0\dot y+1.5y=0.8u. \]

Select a reference model with characteristic polynomial \( (s+2)^3=s^3+6s^2+12s+8 \) and unit steady-state gain:

\[ y_m^{(3)}+6\ddot y_m+12\dot y_m+8y_m=8r. \]

The ideal controller parameters are

\[ \mathbf{k}_x^*=\frac{1}{0.8} \begin{bmatrix}1.5-8\\2.0-12\\1.1-6\end{bmatrix} =\begin{bmatrix}-8.125\\-12.5\\-6.125\end{bmatrix}, \qquad k_r^*=10. \]

This example is structurally valid, yet the required gains are much larger than in the second-order example. The reason is visible from

\[ \|\mathbf{k}_x^*\|_2= \frac{\|\mathbf{a}_p-\mathbf{a}_m\|_2}{|b|}, \qquad |k_r^*|=\frac{|b_m|}{|b|}. \]

Fast reference poles enlarge the coefficient difference, while a small plant high-frequency gain enlarges both feedback and feedforward gains. Thus a mathematically admissible model may still demand excessive control effort or exhibit undesirable transients.

The Lyapunov result gives \( \mathbf{e}(t) \) convergence, but the parameter error appears in the error dynamics only through the scalar product \( \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega} \). If the regressor does not explore all parameter directions, different parameter vectors can generate the same asymptotic input-output behavior.

A standard sufficient condition for identifying all constant parameters is persistent excitation: there must exist constants \( T > 0 \) and \( \alpha > 0 \) such that

\[ \int_t^{t+T}\boldsymbol{\omega}(\tau) \boldsymbol{\omega}^T(\tau)\,d\tau \succeq \alpha\mathbf{I} \qquad \text{for every }t\ge 0. \]

Under suitable boundedness and excitation assumptions, this condition permits \( \tilde{\boldsymbol{\theta} }(t) \) convergence. Without it, the estimates can settle at non-ideal constants while the plant tracks the reference model. The supplied programs make this distinction visible by comparing:

• rich command: \( r(t)=\sin(0.7t)+0.5\sin(1.3t)+0.25\sin(2.1t) \);
• weak command: a constant unit step.

In the step case, the simulated tracking error becomes essentially zero even though the final estimates differ markedly from the ideal parameter vector. This is not a failure of the tracking theorem; it is a reminder that control and identification have different objectives.

flowchart LR
  A["Basic canonical-form \nMRAC assumptions"] --> B["Correct known order and \nexact canonical structure"]
  A --> C["All states measured"]
  A --> D["Known sign of control gain b"]
  A --> E["No omitted actuator or \nhigh-frequency dynamics"]
  A --> F["Command and regressor sufficiently \nexciting for identification"]
  A --> G["Control input not clipped and \nmeasurements reasonably clean"]
  B --> H["Order or structure mismatch \nbreaks exact matching"]
  C --> I["Output-only plant needs \nadditional estimation structure"]
  D --> J["Wrong sign reverses \nLyapunov cancellation"]
  E --> K["Residual dynamics enter \nthe error equation"]
  F --> L["Tracking may converge \nwhile parameters do not"]
  G --> M["Saturation or noise \nadds uncancelled terms \nand poor transients"]
        

9.1 Exact order and matching are indispensable

Suppose the physical actuator has an unmodeled first-order lag:

\[ \dot{\mathbf{x} }=\mathbf{A}_p\mathbf{x}+b\mathbf{e}_n u_a, \qquad \tau_a\dot u_a=-u_a+u_c, \qquad \tau_a > 0. \]

The true closed-loop plant now has order \( n+1 \), while the reference model has order \( n \). No constant n-dimensional feedback vector can make the augmented plant matrix equal to \( \mathbf{A}_m \). The ideal parameter vector assumed in the proof therefore does not exist for this model pair.

9.2 Full-state measurement is a real implementation requirement

The regressor contains every state. For the second-order example, measuring position alone does not provide the velocity term required by the control and update laws. Numerical differentiation generally amplifies noise. An observer- or filter-based adaptive architecture is a different design problem and is treated later in the course.

9.3 The control direction must be correct

The cancellation in \( \dot V \) uses \( \operatorname{sgn}(b) \). If the implemented sign is wrong, the cross terms reinforce rather than cancel each other. Knowing only that \( b\neq 0 \) is insufficient for this basic adaptive law.

9.4 Saturation destroys the nominal cancellation

If the actuator applies \( u_a=\operatorname{sat}(u) \), define \( \Delta_u=u_a-u \). The error model becomes

\[ \dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +b\mathbf{e}_n\tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega} +b\mathbf{e}_n\Delta_u. \]

After the nominal adaptive cancellation,

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

The final term has no fixed sign. Consequently, the nominal theorem cannot be cited unchanged when the commanded control is clipped.

9.5 Adaptation gain is not a free performance knob

Increasing \( \boldsymbol{\Gamma} \) makes the estimates react faster to the filtered error, but also increases sensitivity to measurement noise, discretization error, and abrupt transients. Large gains can make the combined ODE numerically stiff even though the ideal continuous-time proof remains valid.

9.6 A more aggressive reference model can be counterproductive

Faster model poles increase ideal feedback magnitude and often demand larger transient input. Reference-model selection must therefore account for actuator limits, sampling rate, sensor bandwidth, and the uncertainty range rather than being chosen only from a desired settling time.

All implementations integrate the combined seven-state vector

\[ \mathbf{z}= \begin{bmatrix} x_1&x_2&x_{m1}&x_{m2}&\hat k_{x1}&\hat k_{x2}&\hat k_r \end{bmatrix}^T \]

using a fourth-order Runge-Kutta method with \( \Delta t=0.002\text{ s} \). The plant starts at \( [x_1(0),x_2(0)]^T=[0.8,-0.4]^T \); the model and adaptive parameters start at zero. Each program reports the root-mean- square state error over the last 20 percent of the simulation:

\[ e_{\mathrm{RMS,tail} }= \sqrt{\frac{1}{N_t}\sum_{k\in\mathcal{I}_{\mathrm{tail} } } \|\mathbf{e}(t_k)\|_2^2}. \]

A fixed integration step is used so that the Python, C++, Java, and MATLAB trajectories are directly comparable. A smaller step should be used when adaptation gains or model bandwidth are increased.

Required libraries are numpy and matplotlib. The controller and RK4 integrator are implemented from scratch.

Chapter7_Lesson5.py

"""Chapter 7, Lesson 5: Canonical-form MRAC examples and limitations.

This script simulates a second-order SISO plant in controllable canonical form
with a direct Lyapunov MRAC law. Two reference inputs are compared:

1. A multi-sine reference that supplies richer excitation.
2. A constant step reference that generally permits tracking but does not
   guarantee convergence of every adaptive parameter to its ideal value.

Outputs:
    Chapter7_Lesson5_rich.csv
    Chapter7_Lesson5_step.csv
    Chapter7_Lesson5_tracking.png
    Chapter7_Lesson5_parameters.png
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class MRACConfig:
    # True plant: y_ddot + a1*y_dot + a0*y = b*u
    a0: float = 1.2
    a1: float = 0.8
    b: float = 1.5

    # Reference model: y_m_ddot + am1*y_m_dot + am0*y_m = bm*r
    am0: float = 4.0
    am1: float = 4.0
    bm: float = 4.0

    # Q = I gives P = [[9/8, 1/8], [1/8, 5/32]]
    p12: float = 1.0 / 8.0
    p22: float = 5.0 / 32.0

    # Diagonal adaptation gain Gamma
    gamma_x1: float = 20.0
    gamma_x2: float = 20.0
    gamma_r: float = 10.0

    dt: float = 0.002
    final_time: float = 30.0

    @property
    def ideal_parameters(self) -> np.ndarray:
        """Return [k_x1*, k_x2*, k_r*] from the matching conditions."""
        return np.array(
            [
                (self.a0 - self.am0) / self.b,
                (self.a1 - self.am1) / self.b,
                self.bm / self.b,
            ],
            dtype=float,
        )


def rich_reference(t: float) -> float:
    """Bounded multi-sine command used to provide richer excitation."""
    return float(np.sin(0.7 * t) + 0.5 * np.sin(1.3 * t) + 0.25 * np.sin(2.1 * t))


def step_reference(t: float) -> float:
    """Constant command; useful for demonstrating weak excitation."""
    del t
    return 1.0


def closed_loop_derivative(
    t: float,
    z: np.ndarray,
    cfg: MRACConfig,
    reference: Callable[[float], float],
) -> np.ndarray:
    """Evaluate the seven-state plant/model/adaptation dynamics.

    z = [x1, x2, xm1, xm2, kx1_hat, kx2_hat, kr_hat]
    """
    if z.shape != (7,):
        raise ValueError(f"Expected a seven-state vector, received shape {z.shape}.")

    x1, x2, xm1, xm2, kx1_hat, kx2_hat, kr_hat = z
    r = reference(t)

    u = kx1_hat * x1 + kx2_hat * x2 + kr_hat * r

    e1 = x1 - xm1
    e2 = x2 - xm2
    filtered_error = cfg.p12 * e1 + cfg.p22 * e2

    sign_b = 1.0 if cfg.b > 0.0 else -1.0

    dz = np.empty(7, dtype=float)
    dz[0] = x2
    dz[1] = -cfg.a0 * x1 - cfg.a1 * x2 + cfg.b * u
    dz[2] = xm2
    dz[3] = -cfg.am0 * xm1 - cfg.am1 * xm2 + cfg.bm * r
    dz[4] = -sign_b * cfg.gamma_x1 * x1 * filtered_error
    dz[5] = -sign_b * cfg.gamma_x2 * x2 * filtered_error
    dz[6] = -sign_b * cfg.gamma_r * r * filtered_error
    return dz


def rk4_step(
    t: float,
    z: np.ndarray,
    dt: float,
    derivative: Callable[[float, np.ndarray], np.ndarray],
) -> np.ndarray:
    """Advance one fixed-step classical fourth-order Runge-Kutta update."""
    k1 = derivative(t, z)
    k2 = derivative(t + 0.5 * dt, z + 0.5 * dt * k1)
    k3 = derivative(t + 0.5 * dt, z + 0.5 * dt * k2)
    k4 = derivative(t + dt, z + dt * k3)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(
    cfg: MRACConfig,
    reference: Callable[[float], float],
    label: str,
) -> dict[str, np.ndarray | float | str]:
    """Simulate one reference case and return trajectories and diagnostics."""
    if cfg.dt <= 0.0 or cfg.final_time <= 0.0:
        raise ValueError("dt and final_time must both be positive.")

    time = np.arange(0.0, cfg.final_time + 0.5 * cfg.dt, cfg.dt)
    states = np.zeros((time.size, 7), dtype=float)
    control = np.zeros(time.size, dtype=float)
    reference_values = np.zeros(time.size, dtype=float)

    # Nonzero plant initial condition tests the transient response.
    states[0] = np.array([0.8, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0], dtype=float)

    derivative = lambda current_t, current_z: closed_loop_derivative(
        current_t, current_z, cfg, reference
    )

    for index, current_t in enumerate(time):
        x1, x2, _, _, kx1_hat, kx2_hat, kr_hat = states[index]
        r = reference(float(current_t))
        reference_values[index] = r
        control[index] = kx1_hat * x1 + kx2_hat * x2 + kr_hat * r

        if index + 1 < time.size:
            states[index + 1] = rk4_step(
                float(current_t), states[index], cfg.dt, derivative
            )
            if not np.all(np.isfinite(states[index + 1])):
                raise FloatingPointError(
                    f"Non-finite state encountered in the {label!r} simulation "
                    f"at t={current_t:.6f} s."
                )

    error = states[:, 0:2] - states[:, 2:4]
    tail = time >= 0.8 * cfg.final_time
    tail_rms_error = float(np.sqrt(np.mean(np.sum(error[tail] ** 2, axis=1))))

    return {
        "label": label,
        "time": time,
        "states": states,
        "reference": reference_values,
        "control": control,
        "error": error,
        "tail_rms_error": tail_rms_error,
    }


def save_csv(result: dict[str, np.ndarray | float | str], path: Path) -> None:
    """Save a simulation trajectory using a portable comma-separated format."""
    time = np.asarray(result["time"])
    states = np.asarray(result["states"])
    reference_values = np.asarray(result["reference"])
    control = np.asarray(result["control"])
    error = np.asarray(result["error"])

    table = np.column_stack(
        [
            time,
            reference_values,
            states[:, 0],
            states[:, 1],
            states[:, 2],
            states[:, 3],
            error[:, 0],
            error[:, 1],
            states[:, 4],
            states[:, 5],
            states[:, 6],
            control,
        ]
    )
    header = (
        "time,r,x1,x2,xm1,xm2,e1,e2,kx1_hat,kx2_hat,kr_hat,u"
    )
    np.savetxt(path, table, delimiter=",", header=header, comments="")


def create_plots(
    rich: dict[str, np.ndarray | float | str],
    step: dict[str, np.ndarray | float | str],
    cfg: MRACConfig,
    output_directory: Path,
) -> None:
    """Create tracking and parameter-comparison figures."""
    rich_time = np.asarray(rich["time"])
    rich_states = np.asarray(rich["states"])
    rich_reference_values = np.asarray(rich["reference"])
    rich_error = np.asarray(rich["error"])

    step_time = np.asarray(step["time"])
    step_states = np.asarray(step["states"])
    step_error = np.asarray(step["error"])

    figure = plt.figure(figsize=(10.0, 8.0))
    axis1 = figure.add_subplot(3, 1, 1)
    axis1.plot(rich_time, rich_reference_values, label="r")
    axis1.plot(rich_time, rich_states[:, 0], label="y")
    axis1.plot(rich_time, rich_states[:, 2], "--", label="y_m")
    axis1.set_ylabel("Output")
    axis1.set_title("Canonical-form MRAC with rich excitation")
    axis1.grid(True)
    axis1.legend()

    axis2 = figure.add_subplot(3, 1, 2)
    axis2.plot(rich_time, rich_error[:, 0], label="e1")
    axis2.plot(rich_time, rich_error[:, 1], label="e2")
    axis2.set_ylabel("Tracking error")
    axis2.grid(True)
    axis2.legend()

    axis3 = figure.add_subplot(3, 1, 3)
    axis3.plot(step_time, step_error[:, 0], label="step e1")
    axis3.plot(step_time, step_error[:, 1], label="step e2")
    axis3.set_xlabel("Time (s)")
    axis3.set_ylabel("Tracking error")
    axis3.set_title("Weak excitation can still produce good tracking")
    axis3.grid(True)
    axis3.legend()

    figure.tight_layout()
    figure.savefig(output_directory / "Chapter7_Lesson5_tracking.png", dpi=180)
    plt.close(figure)

    ideal = cfg.ideal_parameters
    figure = plt.figure(figsize=(10.0, 7.0))
    axis1 = figure.add_subplot(2, 1, 1)
    for index, name in enumerate(("kx1_hat", "kx2_hat", "kr_hat")):
        axis1.plot(rich_time, rich_states[:, 4 + index], label=name)
        axis1.axhline(ideal[index], linestyle="--", label=f"{name} ideal")
    axis1.set_ylabel("Parameter value")
    axis1.set_title("Parameter evolution under richer excitation")
    axis1.grid(True)
    axis1.legend(ncol=2)

    axis2 = figure.add_subplot(2, 1, 2)
    for index, name in enumerate(("kx1_hat", "kx2_hat", "kr_hat")):
        axis2.plot(step_time, step_states[:, 4 + index], label=name)
        axis2.axhline(ideal[index], linestyle="--", label=f"{name} ideal")
    axis2.set_xlabel("Time (s)")
    axis2.set_ylabel("Parameter value")
    axis2.set_title("A step command need not identify all ideal parameters")
    axis2.grid(True)
    axis2.legend(ncol=2)

    figure.tight_layout()
    figure.savefig(output_directory / "Chapter7_Lesson5_parameters.png", dpi=180)
    plt.close(figure)


def main() -> None:
    cfg = MRACConfig()
    output_directory = Path(__file__).resolve().parent

    rich = simulate(cfg, rich_reference, "rich")
    step = simulate(cfg, step_reference, "step")

    save_csv(rich, output_directory / "Chapter7_Lesson5_rich.csv")
    save_csv(step, output_directory / "Chapter7_Lesson5_step.csv")
    create_plots(rich, step, cfg, output_directory)

    ideal = cfg.ideal_parameters
    rich_final = np.asarray(rich["states"])[-1, 4:7]
    step_final = np.asarray(step["states"])[-1, 4:7]

    print("Ideal parameters [kx1*, kx2*, kr*]:", ideal)
    print("Rich-reference final estimates:", rich_final)
    print("Step-reference final estimates:", step_final)
    print(f"Rich-reference tail RMS state error: {rich['tail_rms_error']:.6f}")
    print(f"Step-reference tail RMS state error: {step['tail_rms_error']:.6f}")
    print("Files written to:", output_directory)


if __name__ == "__main__":
    main()

This C++17 implementation uses only the standard library. It writes CSV files for subsequent plotting in any preferred tool.

Chapter7_Lesson5.cpp

// Chapter 7, Lesson 5: Canonical-form MRAC examples and limitations.
//
// This self-contained C++17 program simulates a second-order SISO plant and
// reference model with a Lyapunov direct-MRAC law. It compares a multi-sine
// command with a constant step command and writes two CSV files.

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

namespace {

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

enum class ReferenceCase { Rich, Step };

struct Config {
    double a0 = 1.2;
    double a1 = 0.8;
    double b = 1.5;

    double am0 = 4.0;
    double am1 = 4.0;
    double bm = 4.0;

    double p12 = 1.0 / 8.0;
    double p22 = 5.0 / 32.0;

    double gamma_x1 = 20.0;
    double gamma_x2 = 20.0;
    double gamma_r = 10.0;

    double dt = 0.002;
    double final_time = 30.0;
};

double reference_value(double t, ReferenceCase reference_case) {
    if (reference_case == ReferenceCase::Step) {
        return 1.0;
    }
    return std::sin(0.7 * t) + 0.5 * std::sin(1.3 * t) +
           0.25 * std::sin(2.1 * t);
}

State add_scaled(const State& first, const State& second, double scale) {
    State result{};
    for (std::size_t i = 0; i < kStateSize; ++i) {
        result[i] = first[i] + scale * second[i];
    }
    return result;
}

State derivative(double t, const State& z, const Config& cfg,
                 ReferenceCase reference_case) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double kx1_hat = z[4];
    const double kx2_hat = z[5];
    const double kr_hat = z[6];

    const double r = reference_value(t, reference_case);
    const double u = kx1_hat * x1 + kx2_hat * x2 + kr_hat * r;

    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double filtered_error = cfg.p12 * e1 + cfg.p22 * e2;
    const double sign_b = cfg.b > 0.0 ? 1.0 : -1.0;

    return State{
        x2,
        -cfg.a0 * x1 - cfg.a1 * x2 + cfg.b * u,
        xm2,
        -cfg.am0 * xm1 - cfg.am1 * xm2 + cfg.bm * r,
        -sign_b * cfg.gamma_x1 * x1 * filtered_error,
        -sign_b * cfg.gamma_x2 * x2 * filtered_error,
        -sign_b * cfg.gamma_r * r * filtered_error,
    };
}

State rk4_step(double t, const State& z, const Config& cfg,
               ReferenceCase reference_case) {
    const State k1 = derivative(t, z, cfg, reference_case);
    const State k2 = derivative(t + 0.5 * cfg.dt,
                                add_scaled(z, k1, 0.5 * cfg.dt), cfg,
                                reference_case);
    const State k3 = derivative(t + 0.5 * cfg.dt,
                                add_scaled(z, k2, 0.5 * cfg.dt), cfg,
                                reference_case);
    const State k4 = derivative(t + cfg.dt, add_scaled(z, k3, cfg.dt), cfg,
                                reference_case);

    State next{};
    for (std::size_t i = 0; i < kStateSize; ++i) {
        next[i] = z[i] + (cfg.dt / 6.0) *
                             (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Non-finite state encountered during RK4 integration.");
        }
    }
    return next;
}

struct SimulationResult {
    State final_state{};
    double tail_rms_error = 0.0;
};

SimulationResult simulate(const Config& cfg, ReferenceCase reference_case,
                          const std::string& output_path) {
    if (cfg.dt <= 0.0 || cfg.final_time <= 0.0) {
        throw std::invalid_argument("dt and final_time must be positive.");
    }

    std::ofstream csv(output_path);
    if (!csv) {
        throw std::runtime_error("Cannot open output file: " + output_path);
    }
    csv << "time,r,x1,x2,xm1,xm2,e1,e2,kx1_hat,kx2_hat,kr_hat,u\n";
    csv << std::setprecision(12);

    State z{0.8, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0};
    const long long number_of_steps =
        static_cast<long long>(std::llround(cfg.final_time / cfg.dt));
    const double tail_start = 0.8 * cfg.final_time;

    double tail_squared_error_sum = 0.0;
    long long tail_count = 0;

    for (long long step = 0; step <= number_of_steps; ++step) {
        const double t = static_cast<double>(step) * cfg.dt;
        const double r = reference_value(t, reference_case);
        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        const double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;

        csv << t << ',' << r << ',' << z[0] << ',' << z[1] << ',' << z[2]
            << ',' << z[3] << ',' << e1 << ',' << e2 << ',' << z[4] << ','
            << z[5] << ',' << z[6] << ',' << u << '\n';

        if (t >= tail_start) {
            tail_squared_error_sum += e1 * e1 + e2 * e2;
            ++tail_count;
        }

        if (step < number_of_steps) {
            z = rk4_step(t, z, cfg, reference_case);
        }
    }

    SimulationResult result;
    result.final_state = z;
    result.tail_rms_error =
        std::sqrt(tail_squared_error_sum / static_cast<double>(tail_count));
    return result;
}

void print_result(const std::string& name, const SimulationResult& result) {
    std::cout << name << " final estimates [kx1_hat, kx2_hat, kr_hat] = ["
              << result.final_state[4] << ", " << result.final_state[5] << ", "
              << result.final_state[6] << "]\n";
    std::cout << name << " tail RMS state error = " << result.tail_rms_error
              << "\n";
}

}  // namespace

int main() {
    try {
        const Config cfg;
        const SimulationResult rich =
            simulate(cfg, ReferenceCase::Rich, "Chapter7_Lesson5_cpp_rich.csv");
        const SimulationResult step =
            simulate(cfg, ReferenceCase::Step, "Chapter7_Lesson5_cpp_step.csv");

        const std::array<double, 3> ideal{
            (cfg.a0 - cfg.am0) / cfg.b,
            (cfg.a1 - cfg.am1) / cfg.b,
            cfg.bm / cfg.b,
        };

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Ideal parameters [kx1*, kx2*, kr*] = [" << ideal[0] << ", "
                  << ideal[1] << ", " << ideal[2] << "]\n";
        print_result("Rich reference", rich);
        print_result("Step reference", step);
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Simulation failed: " << error.what() << '\n';
        return 1;
    }
}

This Java 17 implementation uses only the standard library and writes trajectories to CSV.

Chapter7_Lesson5.java

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

/**
 * Chapter 7, Lesson 5: Canonical-form MRAC examples and limitations.
 *
 * <p>This Java 17 program simulates a second-order plant and reference model
 * with a direct Lyapunov MRAC law. It compares a multi-sine command with a
 * constant step command and writes portable CSV output.</p>
 */
public final class Chapter7_Lesson5 {
    private static final int STATE_SIZE = 7;

    private enum ReferenceCase {
        RICH,
        STEP
    }

    private static final class Config {
        final double a0 = 1.2;
        final double a1 = 0.8;
        final double b = 1.5;

        final double am0 = 4.0;
        final double am1 = 4.0;
        final double bm = 4.0;

        final double p12 = 1.0 / 8.0;
        final double p22 = 5.0 / 32.0;

        final double gammaX1 = 20.0;
        final double gammaX2 = 20.0;
        final double gammaR = 10.0;

        final double dt = 0.002;
        final double finalTime = 30.0;
    }

    private record SimulationResult(double[] finalState, double tailRmsError) {
    }

    private Chapter7_Lesson5() {
        // Utility class.
    }

    private static double referenceValue(double time, ReferenceCase referenceCase) {
        if (referenceCase == ReferenceCase.STEP) {
            return 1.0;
        }
        return Math.sin(0.7 * time)
                + 0.5 * Math.sin(1.3 * time)
                + 0.25 * Math.sin(2.1 * time);
    }

    private static double[] derivative(
            double time,
            double[] state,
            Config config,
            ReferenceCase referenceCase) {
        if (state.length != STATE_SIZE) {
            throw new IllegalArgumentException("Expected a seven-state vector.");
        }

        double x1 = state[0];
        double x2 = state[1];
        double xm1 = state[2];
        double xm2 = state[3];
        double kx1Hat = state[4];
        double kx2Hat = state[5];
        double krHat = state[6];

        double reference = referenceValue(time, referenceCase);
        double control = kx1Hat * x1 + kx2Hat * x2 + krHat * reference;

        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double filteredError = config.p12 * e1 + config.p22 * e2;
        double signB = config.b > 0.0 ? 1.0 : -1.0;

        return new double[] {
            x2,
            -config.a0 * x1 - config.a1 * x2 + config.b * control,
            xm2,
            -config.am0 * xm1 - config.am1 * xm2 + config.bm * reference,
            -signB * config.gammaX1 * x1 * filteredError,
            -signB * config.gammaX2 * x2 * filteredError,
            -signB * config.gammaR * reference * filteredError
        };
    }

    private static double[] addScaled(double[] first, double[] second, double scale) {
        double[] result = new double[STATE_SIZE];
        for (int index = 0; index < STATE_SIZE; index++) {
            result[index] = first[index] + scale * second[index];
        }
        return result;
    }

    private static double[] rk4Step(
            double time,
            double[] state,
            Config config,
            ReferenceCase referenceCase) {
        double[] k1 = derivative(time, state, config, referenceCase);
        double[] k2 = derivative(
                time + 0.5 * config.dt,
                addScaled(state, k1, 0.5 * config.dt),
                config,
                referenceCase);
        double[] k3 = derivative(
                time + 0.5 * config.dt,
                addScaled(state, k2, 0.5 * config.dt),
                config,
                referenceCase);
        double[] k4 = derivative(
                time + config.dt,
                addScaled(state, k3, config.dt),
                config,
                referenceCase);

        double[] next = new double[STATE_SIZE];
        for (int index = 0; index < STATE_SIZE; index++) {
            next[index] = state[index]
                    + (config.dt / 6.0)
                    * (k1[index] + 2.0 * k2[index] + 2.0 * k3[index] + k4[index]);
            if (!Double.isFinite(next[index])) {
                throw new ArithmeticException("Non-finite state encountered during RK4 integration.");
            }
        }
        return next;
    }

    private static SimulationResult simulate(
            Config config,
            ReferenceCase referenceCase,
            Path outputPath) throws IOException {
        if (config.dt <= 0.0 || config.finalTime <= 0.0) {
            throw new IllegalArgumentException("dt and finalTime must be positive.");
        }

        long numberOfSteps = Math.round(config.finalTime / config.dt);
        double tailStart = 0.8 * config.finalTime;
        double[] state = {0.8, -0.4, 0.0, 0.0, 0.0, 0.0, 0.0};

        double tailSquaredErrorSum = 0.0;
        long tailCount = 0L;

        try (BufferedWriter writer = Files.newBufferedWriter(
                outputPath, StandardCharsets.UTF_8)) {
            writer.write("time,r,x1,x2,xm1,xm2,e1,e2,kx1_hat,kx2_hat,kr_hat,u");
            writer.newLine();

            for (long step = 0L; step <= numberOfSteps; step++) {
                double time = step * config.dt;
                double reference = referenceValue(time, referenceCase);
                double e1 = state[0] - state[2];
                double e2 = state[1] - state[3];
                double control = state[4] * state[0]
                        + state[5] * state[1]
                        + state[6] * reference;

                writer.write(String.format(
                        Locale.ROOT,
                        "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                        time,
                        reference,
                        state[0],
                        state[1],
                        state[2],
                        state[3],
                        e1,
                        e2,
                        state[4],
                        state[5],
                        state[6],
                        control));

                if (time >= tailStart) {
                    tailSquaredErrorSum += e1 * e1 + e2 * e2;
                    tailCount++;
                }

                if (step < numberOfSteps) {
                    state = rk4Step(time, state, config, referenceCase);
                }
            }
        }

        double tailRmsError = Math.sqrt(tailSquaredErrorSum / tailCount);
        return new SimulationResult(state.clone(), tailRmsError);
    }

    private static void printResult(String name, SimulationResult result) {
        double[] finalState = result.finalState();
        System.out.printf(
                Locale.ROOT,
                "%s final estimates [kx1_hat, kx2_hat, kr_hat] = [%.6f, %.6f, %.6f]%n",
                name,
                finalState[4],
                finalState[5],
                finalState[6]);
        System.out.printf(
                Locale.ROOT,
                "%s tail RMS state error = %.6f%n",
                name,
                result.tailRmsError());
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Config config = new Config();

        try {
            SimulationResult rich = simulate(
                    config,
                    ReferenceCase.RICH,
                    Path.of("Chapter7_Lesson5_java_rich.csv"));
            SimulationResult step = simulate(
                    config,
                    ReferenceCase.STEP,
                    Path.of("Chapter7_Lesson5_java_step.csv"));

            double idealKx1 = (config.a0 - config.am0) / config.b;
            double idealKx2 = (config.a1 - config.am1) / config.b;
            double idealKr = config.bm / config.b;

            System.out.printf(
                    Locale.ROOT,
                    "Ideal parameters [kx1*, kx2*, kr*] = [%.6f, %.6f, %.6f]%n",
                    idealKx1,
                    idealKx2,
                    idealKr);
            printResult("Rich reference", rich);
            printResult("Step reference", step);
        } catch (IOException | RuntimeException error) {
            System.err.println("Simulation failed: " + error.getMessage());
            System.exit(1);
        }
    }
}

The MATLAB script contains the RK4 solver, two reference experiments, CSV export, diagnostics, and plots.

Chapter7_Lesson5.m

%% Chapter 7, Lesson 5: Canonical-form MRAC examples and limitations
% This script simulates a second-order SISO plant in controllable canonical
% form using a direct Lyapunov MRAC law. It compares a multi-sine reference
% with a constant step reference. The step case illustrates that asymptotic
% tracking does not imply convergence of every adaptive parameter.

clear;
close all;
clc;

cfg.a0 = 1.2;
cfg.a1 = 0.8;
cfg.b = 1.5;

cfg.am0 = 4.0;
cfg.am1 = 4.0;
cfg.bm = 4.0;

% For A_m = [0 1; -4 -4] and Q = I:
% A_m' P + P A_m = -Q gives P = [9/8 1/8; 1/8 5/32].
cfg.p12 = 1/8;
cfg.p22 = 5/32;

cfg.gamma = [20; 20; 10];
cfg.dt = 0.002;
cfg.finalTime = 30.0;

idealParameters = [
    (cfg.a0 - cfg.am0) / cfg.b;
    (cfg.a1 - cfg.am1) / cfg.b;
    cfg.bm / cfg.b
];

rich = simulateCanonicalMRAC(cfg, @richReference, 'rich');
step = simulateCanonicalMRAC(cfg, @stepReference, 'step');

writeResultCsv(rich, 'Chapter7_Lesson5_matlab_rich.csv');
writeResultCsv(step, 'Chapter7_Lesson5_matlab_step.csv');

fprintf('Ideal parameters [kx1*, kx2*, kr*] = [%.6f, %.6f, %.6f]\n', ...
    idealParameters(1), idealParameters(2), idealParameters(3));
fprintf('Rich-reference final estimates = [%.6f, %.6f, %.6f]\n', ...
    rich.z(end, 5), rich.z(end, 6), rich.z(end, 7));
fprintf('Step-reference final estimates = [%.6f, %.6f, %.6f]\n', ...
    step.z(end, 5), step.z(end, 6), step.z(end, 7));
fprintf('Rich-reference tail RMS state error = %.6f\n', rich.tailRmsError);
fprintf('Step-reference tail RMS state error = %.6f\n', step.tailRmsError);

figure('Name', 'Chapter07 Lesson5 Tracking');
tiledlayout(3, 1);

nexttile;
plot(rich.t, rich.r, 'DisplayName', 'r');
hold on;
plot(rich.t, rich.z(:, 1), 'DisplayName', 'y');
plot(rich.t, rich.z(:, 3), '--', 'DisplayName', 'y_m');
grid on;
ylabel('Output');
title('Canonical-form MRAC with rich excitation');
legend('Location', 'best');

nexttile;
plot(rich.t, rich.error(:, 1), 'DisplayName', 'e_1');
hold on;
plot(rich.t, rich.error(:, 2), 'DisplayName', 'e_2');
grid on;
ylabel('Tracking error');
legend('Location', 'best');

nexttile;
plot(step.t, step.error(:, 1), 'DisplayName', 'step e_1');
hold on;
plot(step.t, step.error(:, 2), 'DisplayName', 'step e_2');
grid on;
xlabel('Time (s)');
ylabel('Tracking error');
title('Weak excitation can still produce good tracking');
legend('Location', 'best');
exportgraphics(gcf, 'Chapter7_Lesson5_matlab_tracking.png', 'Resolution', 180);

figure('Name', 'Chapter07 Lesson5 Parameters');
tiledlayout(2, 1);

nexttile;
plotParameterTrajectories(rich.t, rich.z(:, 5:7), idealParameters);
title('Parameter evolution under richer excitation');

nexttile;
plotParameterTrajectories(step.t, step.z(:, 5:7), idealParameters);
title('A step command need not identify all ideal parameters');
xlabel('Time (s)');
exportgraphics(gcf, 'Chapter7_Lesson5_matlab_parameters.png', 'Resolution', 180);

%% Local functions
function result = simulateCanonicalMRAC(cfg, referenceFunction, label)
    if cfg.dt <= 0 || cfg.finalTime <= 0
        error('dt and finalTime must be positive.');
    end

    t = (0:cfg.dt:cfg.finalTime).';
    z = zeros(numel(t), 7);
    z(1, :) = [0.8, -0.4, 0, 0, 0, 0, 0];

    r = zeros(numel(t), 1);
    u = zeros(numel(t), 1);

    for index = 1:numel(t)
        currentTime = t(index);
        r(index) = referenceFunction(currentTime);
        u(index) = z(index, 5) * z(index, 1) ...
            + z(index, 6) * z(index, 2) ...
            + z(index, 7) * r(index);

        if index < numel(t)
            derivative = @(time, state) closedLoopDerivative( ...
                time, state, cfg, referenceFunction);
            z(index + 1, :) = rk4Step( ...
                currentTime, z(index, :).', cfg.dt, derivative).';

            if any(~isfinite(z(index + 1, :)))
                error('Non-finite state in %s simulation at t = %.6f s.', ...
                    label, currentTime);
            end
        end
    end

    trackingError = z(:, 1:2) - z(:, 3:4);
    tailMask = t >= 0.8 * cfg.finalTime;
    tailRmsError = sqrt(mean(sum(trackingError(tailMask, :) .^ 2, 2)));

    result.label = label;
    result.t = t;
    result.z = z;
    result.r = r;
    result.u = u;
    result.error = trackingError;
    result.tailRmsError = tailRmsError;
end

function dz = closedLoopDerivative(time, z, cfg, referenceFunction)
    x1 = z(1);
    x2 = z(2);
    xm1 = z(3);
    xm2 = z(4);
    kx1Hat = z(5);
    kx2Hat = z(6);
    krHat = z(7);

    reference = referenceFunction(time);
    control = kx1Hat * x1 + kx2Hat * x2 + krHat * reference;

    e1 = x1 - xm1;
    e2 = x2 - xm2;
    filteredError = cfg.p12 * e1 + cfg.p22 * e2;
    signB = sign(cfg.b);

    dz = zeros(7, 1);
    dz(1) = x2;
    dz(2) = -cfg.a0 * x1 - cfg.a1 * x2 + cfg.b * control;
    dz(3) = xm2;
    dz(4) = -cfg.am0 * xm1 - cfg.am1 * xm2 + cfg.bm * reference;
    dz(5) = -signB * cfg.gamma(1) * x1 * filteredError;
    dz(6) = -signB * cfg.gamma(2) * x2 * filteredError;
    dz(7) = -signB * cfg.gamma(3) * reference * filteredError;
end

function nextState = rk4Step(time, state, dt, derivative)
    k1 = derivative(time, state);
    k2 = derivative(time + 0.5 * dt, state + 0.5 * dt * k1);
    k3 = derivative(time + 0.5 * dt, state + 0.5 * dt * k2);
    k4 = derivative(time + dt, state + dt * k3);
    nextState = state + (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
end

function value = richReference(time)
    value = sin(0.7 * time) + 0.5 * sin(1.3 * time) ...
        + 0.25 * sin(2.1 * time);
end

function value = stepReference(~)
    value = 1.0;
end

function writeResultCsv(result, fileName)
    outputTable = table( ...
        result.t, ...
        result.r, ...
        result.z(:, 1), ...
        result.z(:, 2), ...
        result.z(:, 3), ...
        result.z(:, 4), ...
        result.error(:, 1), ...
        result.error(:, 2), ...
        result.z(:, 5), ...
        result.z(:, 6), ...
        result.z(:, 7), ...
        result.u, ...
        'VariableNames', { ...
            'time', 'r', 'x1', 'x2', 'xm1', 'xm2', ...
            'e1', 'e2', 'kx1_hat', 'kx2_hat', 'kr_hat', 'u'});
    writetable(outputTable, fileName);
end

function plotParameterTrajectories(time, parameters, idealParameters)
    names = {'k_{x1}', 'k_{x2}', 'k_r'};
    hold on;
    for index = 1:3
        plot(time, parameters(:, index), ...
            'DisplayName', [names{index}, ' estimate']);
        yline(idealParameters(index), '--', ...
            'DisplayName', [names{index}, ' ideal']);
    end
    grid on;
    ylabel('Parameter value');
    legend('Location', 'best', 'NumColumns', 2);
end

The following MATLAB script creates a Simulink model from standard blocks. It implements the plant, reference model, regressor, control dot product, filtered error, adaptive integrator, and data logging. A Simulink license is required.

Chapter7_Lesson5_Simulink.m

%% Programmatic Simulink model for Chapter 7, Lesson 5
% This script creates Chapter7_Lesson5_MRAC.slx using standard Simulink
% blocks. The model implements the same second-order canonical-form MRAC as
% Chapter7_Lesson5.m. Simulink is required.

clear;
clc;

modelName = 'Chapter7_Lesson5_MRAC';

if ~license('test', 'Simulink')
    error('This script requires a Simulink license.');
end

if bdIsLoaded(modelName)
    close_system(modelName, 0);
end

new_system(modelName);
open_system(modelName);

% Model and plant constants.
Aplant = [0, 1; -1.2, -0.8];
Bplant = [0; 1.5];
Amodel = [0, 1; -4, -4];
Bmodel = [0; 4];
Cstate = eye(2);
Dstate = zeros(2, 1);
Prow = [1/8, 5/32];
Gamma = diag([20, 20, 10]);

% Rich reference r(t) = sin(0.7t) + 0.5 sin(1.3t) + 0.25 sin(2.1t).
add_block('simulink/Sources/Sine Wave', [modelName, '/Sine1'], ...
    'Amplitude', '1', 'Frequency', '0.7', ...
    'Position', [30, 40, 80, 70]);
add_block('simulink/Sources/Sine Wave', [modelName, '/Sine2'], ...
    'Amplitude', '0.5', 'Frequency', '1.3', ...
    'Position', [30, 90, 80, 120]);
add_block('simulink/Sources/Sine Wave', [modelName, '/Sine3'], ...
    'Amplitude', '0.25', 'Frequency', '2.1', ...
    'Position', [30, 140, 80, 170]);
add_block('simulink/Math Operations/Add', [modelName, '/ReferenceSum'], ...
    'Inputs', '+++', 'Position', [120, 75, 150, 145]);

% Reference model and true plant, both exposing the complete state vector.
add_block('simulink/Continuous/State-Space', [modelName, '/ReferenceModel'], ...
    'A', 'Amodel', 'B', 'Bmodel', 'C', 'Cstate', 'D', 'Dstate', ...
    'InitialCondition', '[0;0]', 'Position', [210, 40, 350, 100]);
add_block('simulink/Continuous/State-Space', [modelName, '/Plant'], ...
    'A', 'Aplant', 'B', 'Bplant', 'C', 'Cstate', 'D', 'Dstate', ...
    'InitialCondition', '[0.8;-0.4]', 'Position', [690, 230, 830, 290]);

% Tracking error e = x - x_m and scalar filtered error s = e' P e_n.
add_block('simulink/Math Operations/Add', [modelName, '/TrackingError'], ...
    'Inputs', '+-', 'Position', [470, 85, 500, 145]);
add_block('simulink/Math Operations/Gain', [modelName, '/FilteredError'], ...
    'Gain', 'Prow', 'Multiplication', 'Matrix(K*u)', ...
    'Position', [540, 95, 640, 135]);

% Regressor omega = [x1; x2; r].
add_block('simulink/Signal Routing/Mux', [modelName, '/RegressorMux'], ...
    'Inputs', '2', 'Position', [480, 230, 485, 315]);

% Control u = theta_hat' omega.
add_block('simulink/Math Operations/Dot Product', [modelName, '/ControlDot'], ...
    'Position', [580, 235, 630, 285]);

% theta_dot = -Gamma omega s, theta_hat(0) = 0.
add_block('simulink/Math Operations/Product', [modelName, '/OmegaTimesS'], ...
    'Inputs', '**', 'Position', [700, 80, 740, 125]);
add_block('simulink/Math Operations/Gain', [modelName, '/AdaptationGain'], ...
    'Gain', '-Gamma', 'Multiplication', 'Matrix(K*u)', ...
    'Position', [780, 80, 880, 125]);
add_block('simulink/Continuous/Integrator', [modelName, '/ThetaIntegrator'], ...
    'InitialCondition', '[0;0;0]', 'Position', [920, 80, 960, 125]);

% Data logging.
add_block('simulink/Sinks/To Workspace', [modelName, '/LogReference'], ...
    'VariableName', 'r_out', 'SaveFormat', 'Timeseries', ...
    'Position', [210, 145, 310, 175]);
add_block('simulink/Sinks/To Workspace', [modelName, '/LogModelState'], ...
    'VariableName', 'xm_out', 'SaveFormat', 'Timeseries', ...
    'Position', [380, 35, 480, 65]);
add_block('simulink/Sinks/To Workspace', [modelName, '/LogPlantState'], ...
    'VariableName', 'x_out', 'SaveFormat', 'Timeseries', ...
    'Position', [870, 245, 970, 275]);
add_block('simulink/Sinks/To Workspace', [modelName, '/LogTheta'], ...
    'VariableName', 'theta_out', 'SaveFormat', 'Timeseries', ...
    'Position', [1000, 80, 1100, 110]);
add_block('simulink/Sinks/To Workspace', [modelName, '/LogControl'], ...
    'VariableName', 'u_out', 'SaveFormat', 'Timeseries', ...
    'Position', [650, 315, 750, 345]);

% Source and reference-model connections.
add_line(modelName, 'Sine1/1', 'ReferenceSum/1', 'autorouting', 'on');
add_line(modelName, 'Sine2/1', 'ReferenceSum/2', 'autorouting', 'on');
add_line(modelName, 'Sine3/1', 'ReferenceSum/3', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'ReferenceModel/1', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'LogReference/1', 'autorouting', 'on');
add_line(modelName, 'ReferenceSum/1', 'RegressorMux/2', 'autorouting', 'on');
add_line(modelName, 'ReferenceModel/1', 'TrackingError/2', 'autorouting', 'on');
add_line(modelName, 'ReferenceModel/1', 'LogModelState/1', 'autorouting', 'on');

% Plant, error, and regressor connections.
add_line(modelName, 'Plant/1', 'TrackingError/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'RegressorMux/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'LogPlantState/1', 'autorouting', 'on');
add_line(modelName, 'TrackingError/1', 'FilteredError/1', 'autorouting', 'on');

% Adaptive law and control connections.
add_line(modelName, 'RegressorMux/1', 'ControlDot/1', 'autorouting', 'on');
add_line(modelName, 'ThetaIntegrator/1', 'ControlDot/2', 'autorouting', 'on');
add_line(modelName, 'ControlDot/1', 'Plant/1', 'autorouting', 'on');
add_line(modelName, 'ControlDot/1', 'LogControl/1', 'autorouting', 'on');
add_line(modelName, 'RegressorMux/1', 'OmegaTimesS/1', 'autorouting', 'on');
add_line(modelName, 'FilteredError/1', 'OmegaTimesS/2', 'autorouting', 'on');
add_line(modelName, 'OmegaTimesS/1', 'AdaptationGain/1', 'autorouting', 'on');
add_line(modelName, 'AdaptationGain/1', 'ThetaIntegrator/1', 'autorouting', 'on');
add_line(modelName, 'ThetaIntegrator/1', 'LogTheta/1', 'autorouting', 'on');

set_param(modelName, ...
    'Solver', 'ode45', ...
    'StopTime', '30', ...
    'SaveTime', 'on', ...
    'TimeSaveName', 'tout');

save_system(modelName, [modelName, '.slx']);
set_param(modelName, 'SimulationCommand', 'update');

fprintf('Created %s.slx. Run sim(''%s'') to execute it.\n', ...
    modelName, modelName);

The notebook file loads the accompanying Wolfram Language source from the same directory. The source uses NDSolveValue, exports CSV files and figures, and compares the two excitation cases.

Chapter7_Lesson5.nb

Notebook[{
 Cell["Chapter 7, Lesson 5: Example Designs and Limitations of Canonical-Form MRAC", "Title"],
 Cell["This notebook loads and evaluates the accompanying Chapter7_Lesson5.wl source file. Keep both files in the same directory.", "Text"],
 Cell[BoxData[
  RowBox[{"Get", "[", 
   RowBox[{"FileNameJoin", "[", 
    RowBox[{"{", 
     RowBox[{"NotebookDirectory", "[", "]", ",", "\"Chapter7_Lesson5.wl\""}], 
     "}"}], "]"}], "]"}]], "Input"]
},
WindowSize->{1100, 800},
WindowMargins->{ {Automatic, 80}, {Automatic, 40} },
FrontEndVersion->"14.0 for Linux x86 (64-bit)",
StyleDefinitions->"Default.nb"
]

Chapter7_Lesson5.wl

(* Chapter 7, Lesson 5: Canonical-form MRAC examples and limitations. *)

ClearAll["Global`*"];

(* True plant and reference-model coefficients. *)
a0 = 1.2;
a1 = 0.8;
b = 1.5;
am0 = 4.0;
am1 = 4.0;
bm = 4.0;

(* P solves Transpose[Am].P + P.Am == -IdentityMatrix[2]. *)
p12 = 1/8;
p22 = 5/32;

gammaX1 = 20.0;
gammaX2 = 20.0;
gammaR = 10.0;
finalTime = 30.0;

idealParameters = {(a0 - am0)/b, (a1 - am1)/b, bm/b};

richReference[t_] := Sin[0.7 t] + 0.5 Sin[1.3 t] + 0.25 Sin[2.1 t];
stepReference[t_] := 1.0;

ClearAll[simulateCanonicalMRAC];
simulateCanonicalMRAC[referenceFunction_, label_String] := Module[
  {
    x1, x2, xm1, xm2, kx1, kx2, kr, t, equations, solution,
    sampleTimes, sampleTable, tailRows, tailRmsError
  },

  equations = {
    x1'[t] == x2[t],
    x2'[t] == -a0 x1[t] - a1 x2[t]
      + b (kx1[t] x1[t] + kx2[t] x2[t] + kr[t] referenceFunction[t]),
    xm1'[t] == xm2[t],
    xm2'[t] == -am0 xm1[t] - am1 xm2[t] + bm referenceFunction[t],
    kx1'[t] == -gammaX1 x1[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    kx2'[t] == -gammaX2 x2[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    kr'[t] == -gammaR referenceFunction[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    x1[0] == 0.8,
    x2[0] == -0.4,
    xm1[0] == 0.0,
    xm2[0] == 0.0,
    kx1[0] == 0.0,
    kx2[0] == 0.0,
    kr[0] == 0.0
  };

  solution = NDSolveValue[
    equations,
    {x1, x2, xm1, xm2, kx1, kx2, kr},
    {t, 0, finalTime},
    Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4} },
    MaxSteps -> Infinity
  ];

  sampleTimes = Range[0.0, finalTime, 0.01];
  sampleTable = Table[
    With[
      {
        plant1 = solution[[1]][time],
        plant2 = solution[[2]][time],
        model1 = solution[[3]][time],
        model2 = solution[[4]][time],
        theta1 = solution[[5]][time],
        theta2 = solution[[6]][time],
        thetaR = solution[[7]][time],
        reference = referenceFunction[time]
      },
      {
        time,
        reference,
        plant1,
        plant2,
        model1,
        model2,
        plant1 - model1,
        plant2 - model2,
        theta1,
        theta2,
        thetaR,
        theta1 plant1 + theta2 plant2 + thetaR reference
      }
    ],
    {time, sampleTimes}
  ];

  tailRows = Select[sampleTable, First[#] >= 0.8 finalTime &];
  tailRmsError = Sqrt[Mean[(#[[7]]^2 + #[[8]]^2) & /@ tailRows]];

  <|
    "Label" -> label,
    "Solution" -> solution,
    "Samples" -> sampleTable,
    "TailRMSError" -> tailRmsError
  |>
];

richResult = simulateCanonicalMRAC[richReference, "rich"];
stepResult = simulateCanonicalMRAC[stepReference, "step"];

header = {
  "time", "r", "x1", "x2", "xm1", "xm2", "e1", "e2",
  "kx1_hat", "kx2_hat", "kr_hat", "u"
};

Export[
  FileNameJoin[{Directory[], "Chapter7_Lesson5_mathematica_rich.csv"}],
  Prepend[richResult["Samples"], header]
];
Export[
  FileNameJoin[{Directory[], "Chapter7_Lesson5_mathematica_step.csv"}],
  Prepend[stepResult["Samples"], header]
];

Print["Ideal parameters [kx1*, kx2*, kr*] = ", N[idealParameters, 8]];
Print[
  "Rich-reference final estimates = ",
  N[Through[richResult["Solution"][[5 ;; 7]][finalTime]], 8]
];
Print[
  "Step-reference final estimates = ",
  N[Through[stepResult["Solution"][[5 ;; 7]][finalTime]], 8]
];
Print["Rich-reference tail RMS state error = ", richResult["TailRMSError"]];
Print["Step-reference tail RMS state error = ", stepResult["TailRMSError"]];

richTrackingPlot = Plot[
  Evaluate[{
    richResult["Solution"][[1]][t],
    richResult["Solution"][[3]][t],
    richReference[t]
  }],
  {t, 0, finalTime},
  PlotLegends -> {"y", "y_m", "r"},
  PlotRange -> All,
  Frame -> True,
  FrameLabel -> {"Time (s)", "Output"},
  PlotLabel -> "Canonical-form MRAC with rich excitation",
  ImageSize -> Large
];

parameterPlot = Plot[
  Evaluate[Join[
    Through[richResult["Solution"][[5 ;; 7]][t]],
    idealParameters
  ]],
  {t, 0, finalTime},
  PlotLegends -> {
    "kx1 estimate", "kx2 estimate", "kr estimate",
    "kx1 ideal", "kx2 ideal", "kr ideal"
  },
  PlotRange -> All,
  Frame -> True,
  FrameLabel -> {"Time (s)", "Parameter value"},
  PlotLabel -> "Parameter evolution under richer excitation",
  ImageSize -> Large
];

Export[
  FileNameJoin[{Directory[], "Chapter7_Lesson5_mathematica_tracking.png"}],
  richTrackingPlot,
  ImageResolution -> 180
];
Export[
  FileNameJoin[{Directory[], "Chapter7_Lesson5_mathematica_parameters.png"}],
  parameterPlot,
  ImageResolution -> 180
];

Column[{richTrackingPlot, parameterPlot}]

Problem 1 (General Matching Parameters): For the n-th order canonical plant and model in Section 2, derive \( \mathbf{k}_x^* \) and \( k_r^* \) directly from the final rows of the matching equations.

Solution: The final row of

\[ \mathbf{A}_p+b\mathbf{e}_n\mathbf{k}_x^{*T}=\mathbf{A}_m \]

is

\[ -\mathbf{a}_p^T+b\mathbf{k}_x^{*T}=-\mathbf{a}_m^T. \]

Therefore,

\[ b\mathbf{k}_x^{*T}=\mathbf{a}_p^T-\mathbf{a}_m^T, \qquad \mathbf{k}_x^*=\frac{\mathbf{a}_p-\mathbf{a}_m}{b}. \]

Similarly, \( b\mathbf{e}_n k_r^*=b_m\mathbf{e}_n \) gives \( k_r^*=b_m/b \).

Problem 2 (Lyapunov Cross-Term Cancellation): Starting from the error equation, show explicitly why the sign of \( b \) is sufficient for the adaptive law even when its magnitude is unknown.

Solution: The plant-dependent cross term is

\[ 2b s_e\tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega} =2\operatorname{sgn}(b)|b|s_e \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}. \]

The parameter part of the Lyapunov derivative is

\[ 2|b|\tilde{\boldsymbol{\theta} }^T \boldsymbol{\Gamma}^{-1} \left[-\operatorname{sgn}(b)\boldsymbol{\Gamma} \boldsymbol{\omega}s_e\right] =-2\operatorname{sgn}(b)|b|s_e \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}. \]

The terms cancel without using the numerical value of \( |b| \). The correct sign is nevertheless essential.

Problem 3 (Second-Order Lyapunov Matrix): Verify that the matrix found in Section 6 is positive definite.

Solution: For

\[ \mathbf{P}=\begin{bmatrix}9/8&1/8\\1/8&5/32\end{bmatrix}, \]

the leading principal minors are

\[ \frac{9}{8} > 0, \qquad \det(\mathbf{P})= \frac{9}{8}\frac{5}{32}-\left(\frac{1}{8}\right)^2 =\frac{41}{256} > 0. \]

Sylvester's criterion therefore gives \( \mathbf{P}\succ 0 \).

Problem 4 (Perfect Tracking with Wrong Parameters): Explain mathematically how a step command can yield zero tracking error while the three parameter estimates do not equal their ideal values.

Solution: When tracking is exact, the error dynamics require

\[ \tilde{\boldsymbol{\theta} }^T\boldsymbol{\omega}(t)=0. \]

This is one scalar constraint on a three-dimensional parameter error. Under a constant command, the asymptotic regressor may span only a lower-dimensional subspace. Every parameter error in the orthogonal complement of that subspace is invisible to the tracking equation. Persistent excitation is needed to eliminate all such directions.

Problem 5 (Actuator-Lag Mismatch): A designer applies the second-order controller to a plant whose actuator satisfies \( \tau_a\dot u_a=-u_a+u_c \). Identify the failed assumption and explain why the original ideal parameter vector cannot exactly match the reference model.

Solution: The actual state is

\[ \mathbf{x}_a= \begin{bmatrix}x_1&x_2&u_a\end{bmatrix}^T, \]

so the physical plant is third order, while the chosen model remains second order. The original matching equation equates matrices of the wrong state dimension and assumes the commanded input acts directly on the second plant equation. Both assumptions are false. No constant three-parameter vector in the original controller can remove the actuator pole.

Problem 6 (Reference-Model Trade-Off): For a family of second-order models \( s^2+2\zeta\omega_n s+\omega_n^2 \), determine how the ideal gains scale as \( \omega_n \) increases.

Solution: With plant coefficients \( a_0,a_1 \),

\[ k_{x1}^*=\frac{a_0-\omega_n^2}{b}, \qquad k_{x2}^*=\frac{a_1-2\zeta\omega_n}{b}. \]

Thus the position-feedback magnitude grows quadratically with \( \omega_n \), while the velocity-feedback magnitude grows approximately linearly. A faster model can therefore require rapidly increasing control authority.

Canonical-form MRAC converts coefficient uncertainty into a linear parameter-error term and uses a Lyapunov-designed update law to cancel that term in the stability derivative. Under exact matching, bounded reference signals, full-state measurement, and a known control-gain sign, all closed-loop signals remain bounded and the plant state tracks the reference-model state asymptotically. The result does not by itself identify the true or ideal controller parameters; that stronger claim requires sufficient excitation. The same derivation also exposes the method's boundaries: incorrect order, actuator dynamics, unavailable states, saturation, noise, excessive adaptation gain, and an overly aggressive reference model can invalidate assumptions or degrade practical behavior. These limitations motivate the robust modifications introduced in the next chapter.