After completing this lesson, students should be able to:

• formulate an MRAC tracking problem using only measured input and output signals;
• construct stable filtered regressors that replace unavailable plant states;
• state the matching, minimum-phase, relative-degree, and gain-sign assumptions;
• derive a direct output-feedback controller of the form \(u=\boldsymbol{\theta}^{T}\boldsymbol{\omega}\);
• prove stability of an SPR error model with a quadratic Lyapunov function;
• derive the augmented-error identity and a normalized gradient law;
• distinguish output tracking from convergence of the adaptive parameters; and
• implement the conceptual design in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The lesson does not construct a full adaptive observer. Instead, it uses stable filters driven by \(u(t)\) and \(y(t)\). This is consistent with the preceding lessons: output feedback creates additional dynamics, and the usual linear separation principle cannot simply be invoked for the nonlinear adaptive closed loop.

Consider the unknown SISO plant

\[ y_p(s)=W_p(s)u(s),\qquad W_p(s)=k_p\frac{Z_p(s)}{R_p(s)}, \]

where \(R_p(s)\) is monic of degree \(n\), \(Z_p(s)\) is monic of degree \(m\), and the relative degree is \(n^\star=n-m\). The desired response is generated by

\[ y_m(s)=W_m(s)r(s),\qquad W_m(s)=k_m\frac{Z_m(s)}{R_m(s)}, \]

with \(R_m(s)\) Hurwitz. For the basic design in this lesson, assume:

1. the plant order \(n\) and relative degree \(n^\star=1\) are known;
2. \(Z_p(s)\) is Hurwitz, so the plant is minimum phase;
3. the sign of the high-frequency gain \(k_p\) is known;
4. the plant is controllable and observable in a minimal realization;
5. the command \(r(t)\) is bounded and piecewise continuous;
6. there exists a constant ideal parameter vector \(\boldsymbol{\theta}^{\star}\) satisfying the model-matching equations.

These assumptions are not merely technical decoration. Unknown control direction is treated separately in Chapter 16, while unstable plant zeros generally prevent exact stable inversion and therefore invalidate the simple model-matching construction used here.

Let \(\Lambda(s)\) be an arbitrary monic Hurwitz polynomial of degree \(n-1\), and define

\[ \boldsymbol{\alpha}(s)= \begin{bmatrix} s^{n-2} & s^{n-3} & \cdots & s & 1 \end{bmatrix}^{T}. \]

The controller filters are

\[ \boldsymbol{\omega}_1(t)= \frac{\boldsymbol{\alpha}(s)}{\Lambda(s)}u(t), \qquad \boldsymbol{\omega}_2(t)= \frac{\boldsymbol{\alpha}(s)}{\Lambda(s)}y_p(t). \]

Every component is the output of a strictly proper stable filter. The measurable regressor is

\[ \boldsymbol{\omega}(t)= \begin{bmatrix} \boldsymbol{\omega}_1^{T}(t) & \boldsymbol{\omega}_2^{T}(t) & y_p(t) & r(t) \end{bmatrix}^{T} \in\mathbb{R}^{2n}. \]

The direct adaptive controller is parameterized as

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

Unlike indirect adaptive control, the entries of \(\boldsymbol{\theta}\) are controller coefficients. They need not be interpreted as estimates of the physical plant coefficients.

3.1 Filter-state realization

If \((A_f,b_f,c_f)\) realizes \(1/\Lambda(s)\), a scalar filtered signal \(\omega_q=1/\Lambda(s)[q]\) can be implemented as

\[ \dot{x}_f=A_fx_f+b_fq,\qquad \omega_q=c_f^{T}x_f. \]

Since \(A_f\) is Hurwitz, bounded input implies bounded filter state. Therefore, once boundedness of \(u\) and \(y_p\) is established, all components of \(\boldsymbol{\omega}\) are bounded.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> YM["Desired output ym"]
  R --> REG["Measured-signal regressor"]
  Y["Measured plant output y"] --> REG
  U["Applied input u"] --> REG
  REG --> CTRL["Control law u = thetaT omega"]
  CTRL --> U
  U --> PLANT["Unknown minimum-phase plant"]
  PLANT --> Y
  Y --> ERR["Tracking and augmented-error generator"]
  YM --> ERR
  REG --> ERR
  U --> ERR
  ERR --> ADAPT["Normalized adaptive law"]
  ADAPT --> CTRL
        

There are two dynamic mechanisms besides the plant: stable regressor filters and the adaptive law. Their states interact nonlinearly because the current parameter vector multiplies filtered signals. Consequently, the closed loop is not the cascade of an independently designed observer and an independently designed state-feedback controller.

Under the assumptions of Section 2, polynomial model-matching equations admit a constant vector \(\boldsymbol{\theta}^{\star}\) for which the ideal controller

\[ u^{\star}(t)= {\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\omega}(t) \]

produces the reference transfer function. Define the parameter error

\[ \widetilde{\boldsymbol{\theta}}(t)= \boldsymbol{\theta}(t)-\boldsymbol{\theta}^{\star}. \]

For the relative-degree-one case, the tracking error can be represented conceptually as

\[ e(t)=y_p(t)-y_m(t) = k_pW_e(s) \left[ \widetilde{\boldsymbol{\theta}}^{T}(t) \boldsymbol{\omega}(t) \right], \]

where \(W_e(s)\) is a stable transfer function determined by the reference model and the selected filter polynomial. A clean Lyapunov proof is obtained when \(W_e(s)\) is strictly positive real.

A stable proper transfer function \(W_e(s)\) is strictly positive real (SPR) when, in addition to the standard pole and relative-degree conditions,

\[ \operatorname{Re}\{W_e(j\omega)\}>0 \qquad \text{for every real }\omega. \]

Let a minimal realization of the error model be

\[ \dot{x}_e=A_ex_e+b_ek_p \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}, \qquad e=c_e^{T}x_e. \]

By the Kalman-Yakubovich-Popov property for an SPR realization, for any suitable \(Q=Q^{T}>0\) there exists \(P=P^{T}>0\) such that

\[ A_e^{T}P+PA_e=-Q,\qquad Pb_e=c_e. \]

Choose the direct adaptive law

\[ \dot{\boldsymbol{\theta}} = -\Gamma\operatorname{sgn}(k_p) \boldsymbol{\omega}e, \qquad \Gamma=\Gamma^{T}>0. \]

6.1 Lyapunov proof

Consider

\[ V=x_e^{T}Px_e+ |k_p| \widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1} \widetilde{\boldsymbol{\theta}}. \]

Because \(\boldsymbol{\theta}^{\star}\) is constant,

\[ \begin{aligned} \dot{V} &= x_e^{T}(A_e^{T}P+PA_e)x_e +2k_pe\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} +2|k_p| \widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\boldsymbol{\theta}}\\ &= -x_e^{T}Qx_e +2k_pe\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} -2|k_p|\operatorname{sgn}(k_p)e\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}\\ &=-x_e^{T}Qx_e\le 0. \end{aligned} \]

Hence \(x_e\) and \(\widetilde{\boldsymbol{\theta}}\) are bounded and \(x_e\in\mathcal{L}_2\). Under the usual bounded-signal continuation argument, \(\dot{x}_e\) is bounded. Barbalat's lemma then gives

\[ \lim_{t\to\infty}x_e(t)=0, \qquad \lim_{t\to\infty}e(t)=0. \]

This result guarantees tracking, but not necessarily \(\widetilde{\boldsymbol{\theta}}\to 0\). Parameter convergence additionally requires sufficient excitation.

The compact expression \(W_m(s)[\boldsymbol{\theta}^{T}\boldsymbol{\omega}]\) cannot generally be rewritten as \(\boldsymbol{\theta}^{T}W_m(s)[\boldsymbol{\omega}]\) because \(\boldsymbol{\theta}(t)\) varies with time. A dynamic linear operator commutes with constant coefficients, not with arbitrary time-varying coefficients.

Define the filtered regressor and filtered control signal

\[ \boldsymbol{\zeta}(t)=W_m(s)[\boldsymbol{\omega}(t)], \qquad v(t)=W_m(s)[u(t)] =W_m(s)[\boldsymbol{\theta}^{T}(t)\boldsymbol{\omega}(t)]. \]

Now define the augmented error

\[ e_a(t)=e(t)+ \boldsymbol{\theta}^{T}(t)\boldsymbol{\zeta}(t)-v(t). \]

If the ideal matching relation yields \(e=W_m(s)[\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}]\), then the constancy of \(\boldsymbol{\theta}^{\star}\) gives

\[ \begin{aligned} e_a &= W_m(s)[\boldsymbol{\theta}^{T}\boldsymbol{\omega}] -W_m(s)[{\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\omega}] +\boldsymbol{\theta}^{T}W_m(s)[\boldsymbol{\omega}] -W_m(s)[\boldsymbol{\theta}^{T}\boldsymbol{\omega}]\\ &= \boldsymbol{\theta}^{T}\boldsymbol{\zeta} -{\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\zeta}\\ &= \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\zeta}. \end{aligned} \]

The augmented error therefore converts a dynamic error relation into an algebraic regression error in the unknown parameter mismatch.

Choose the normalization signal

\[ m^2(t)=1+ \boldsymbol{\omega}^{T}\boldsymbol{\omega} +\boldsymbol{\zeta}^{T}\boldsymbol{\zeta}, \qquad m^2(t)\ge 1, \]

and the adaptive law

\[ \dot{\boldsymbol{\theta}} = -\Gamma \frac{\boldsymbol{\zeta}e_a}{m^2}. \]

For the parameter-error Lyapunov function

\[ V_\theta= \frac{1}{2} \widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1} \widetilde{\boldsymbol{\theta}}, \]

the augmented-error identity gives

\[ \dot{V}_\theta = \widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\boldsymbol{\theta}} = -\frac{ \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\zeta}e_a }{m^2} = -\frac{e_a^2}{m^2}\le 0. \]

Thus the parameter vector is bounded and \(e_a/m\in\mathcal{L}_2\). To conclude \(e(t)\to0\), the remaining closed-loop signals must be shown bounded and the filtered swapping term must vanish. This is one of the central technical differences between state-feedback MRAC and output-feedback MRAC.

8.1 Projection and robust modifications

If a compact admissible set \(\Omega_\theta\) is known, replace the raw update by

\[ \dot{\boldsymbol{\theta}} = \operatorname{Proj}_{\Omega_\theta} \left( \boldsymbol{\theta}, -\Gamma\frac{\boldsymbol{\zeta}e_a}{m^2} \right). \]

The projection inequality preserves the sign needed in the Lyapunov derivative. Under noise or neglected high-frequency dynamics, a leakage term may be added:

\[ \dot{\boldsymbol{\theta}} = -\Gamma\frac{\boldsymbol{\zeta}e_a}{m^2} -\sigma\Gamma\boldsymbol{\theta}, \qquad \sigma>0. \]

Leakage generally changes asymptotic convergence into convergence to a residual set, but it reduces parameter drift.

flowchart TD
  A["Specify plant order, relative degree, gain sign, and zero assumptions"] --> B["Choose stable reference model"]
  B --> C["Choose Hurwitz filter polynomial Lambda"]
  C --> D["Generate filters from measured u and y"]
  D --> E["Form regressor omega and controller u = thetaT omega"]
  E --> F["Solve ideal model-matching equations symbolically"]
  F --> G["Check SPR condition or construct augmented error"]
  G --> H["Select Gamma, normalization, and optional projection"]
  H --> I["Prove boundedness before claiming tracking convergence"]
  I --> J["Simulate with command changes, noise tests, and saturation monitoring"]
        

1. Verify structural feasibility: minimum phase, known relative degree, known gain sign.
2. Select \(W_m(s)\): do not demand bandwidth beyond neglected plant dynamics.
3. Select \(\Lambda(s)\): its poles should be stable and normally faster than the reference dynamics, but not so fast that measurement noise is strongly amplified.
4. Build \(\boldsymbol{\omega}\): use only measured \(u\), measured \(y\), and known \(r\).
5. Establish matching: solve for \(\boldsymbol{\theta}^{\star}\) theoretically; the implementation does not use it.
6. Choose the error construction: direct SPR error when available; otherwise augmented-error filtering.
7. Normalize: prevent large regressors from causing excessively large parameter rates.
8. Protect the implementation: use projection, monitoring, and actuator limits when required, while recognizing that hard saturation changes the nominal proof.

Consider

\[ W_p(s)= k_p\frac{s+z_p}{s^2+a_1s+a_0}, \qquad k_p>0,\quad z_p>0, \]

choose \(\Lambda(s)=s+\lambda\), and define

\[ \omega_u=\frac{1}{s+\lambda}u,\qquad \omega_y=\frac{1}{s+\lambda}y. \]

The parameterized controller is

\[ u=\theta_1\omega_u+\theta_2\omega_y+\theta_3y+\theta_4r. \]

Eliminating the filters gives the closed-loop transfer function

\[ \frac{y}{r} = \frac{ k_p\theta_4(s+z_p)(s+\lambda) }{ D_{\mathrm{cl}}(s) }, \]

with

\[ D_{\mathrm{cl}}(s)= (s+\lambda-\theta_1)(s^2+a_1s+a_0) -k_p(s+z_p)\left[\theta_3(s+\lambda)+\theta_2\right]. \]

To obtain the first-order reference model \(W_m(s)=k_m/(s+a_m)\), impose

\[ D_{\mathrm{cl}}^{\star}(s)= (s+z_p)(s+\lambda)(s+a_m), \qquad \theta_4^{\star}=\frac{k_m}{k_p}. \]

Coefficient matching produces

\[ \begin{aligned} a_1+\lambda-\theta_1^{\star}-k_p\theta_3^{\star} &=z_p+\lambda+a_m,\\ a_0+a_1(\lambda-\theta_1^{\star}) -k_p\left[\theta_3^{\star}(\lambda+z_p)+\theta_2^{\star}\right] &=z_p\lambda+z_pa_m+\lambda a_m,\\ a_0(\lambda-\theta_1^{\star}) -k_pz_p(\lambda\theta_3^{\star}+\theta_2^{\star}) &=z_p\lambda a_m. \end{aligned} \]

For

\[ a_1=3,\quad a_0=2,\quad k_p=2,\quad z_p=1.5, \quad \lambda=4,\quad a_m=k_m=1.2, \]

the ideal controller vector is

\[ \boldsymbol{\theta}^{\star} = \begin{bmatrix} 2.5 & 3.0 & -1.1 & 0.6 \end{bmatrix}^{T}. \]

The simulations below do not initialize the adaptive law at this vector. In addition, a piecewise-constant command is not persistently exciting enough to guarantee recovery of every ideal coefficient. Consequently, accurate tracking can occur even when the final adaptive parameter vector differs from \(\boldsymbol{\theta}^{\star}\).

This implementation uses NumPy for vector operations, SciPy for adaptive closed-loop integration, and Matplotlib for visualization. The controller uses only the simulated measurements \(u\) and \(y\).

Chapter15_Lesson4.py

"""
Chapter15_Lesson4.py
Adaptive output-feedback MRAC using filtered regressors and an augmented error.

Dependencies:
    numpy
    scipy
    matplotlib

The plant is
    Wp(s) = 2 (s + 1.5) / (s^2 + 3 s + 2)

Only the plant input u and output y are used by the adaptive controller.
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


@dataclass(frozen=True)
class Parameters:
    a1: float = 3.0
    a0: float = 2.0
    kp: float = 2.0
    zp: float = 1.5
    filter_pole: float = 4.0
    am: float = 1.2
    km: float = 1.2
    gamma: tuple[float, float, float, float] = (8.0, 8.0, 5.0, 5.0)


P = Parameters()
GAMMA = np.diag(P.gamma)


def reference(t: float) -> float:
    """Piecewise-constant command with enough changes to exercise adaptation."""
    if t < 15.0:
        return 1.0
    if t < 30.0:
        return -0.5
    return 0.8


def dynamics(t: float, state: np.ndarray) -> np.ndarray:
    """
    State layout:
      0:2   plant realization states
      2     reference-model output
      3:5   controller filters omega_u and omega_y
      5:9   zeta = Wm(s)[omega], four filter states
      9     v = Wm(s)[u]
      10:14 adaptive controller parameters
    """
    xp = state[0:2]
    ym = state[2]
    omega_u, omega_y = state[3:5]
    zeta = state[5:9]
    v = state[9]
    theta = state[10:14]

    r = reference(t)
    y = P.kp * (xp[0] + P.zp * xp[1])
    omega = np.array([omega_u, omega_y, y, r], dtype=float)
    u = float(theta @ omega)

    tracking_error = y - ym
    augmented_error = tracking_error + float(theta @ zeta) - v
    normalization = 1.0 + float(omega @ omega) + float(zeta @ zeta)

    xp_dot = np.array(
        [-P.a1 * xp[0] - P.a0 * xp[1] + u, xp[0]],
        dtype=float,
    )
    ym_dot = -P.am * ym + P.km * r
    omega_u_dot = -P.filter_pole * omega_u + u
    omega_y_dot = -P.filter_pole * omega_y + y
    zeta_dot = -P.am * zeta + P.km * omega
    v_dot = -P.am * v + P.km * u
    theta_dot = -(GAMMA @ zeta) * augmented_error / normalization

    return np.concatenate(
        (
            xp_dot,
            np.array([ym_dot, omega_u_dot, omega_y_dot]),
            zeta_dot,
            np.array([v_dot]),
            theta_dot,
        )
    )


def run_simulation(
    final_time: float = 45.0,
    max_step: float = 0.01,
) -> tuple[np.ndarray, np.ndarray]:
    initial_state = np.zeros(14, dtype=float)

    # A modest stabilizing initial guess. The ideal parameters are
    # [2.5, 3.0, -1.1, 0.6], but the controller is not initialized there.
    initial_state[10:14] = np.array([0.0, 0.0, -0.5, 0.5])

    solution = solve_ivp(
        dynamics,
        (0.0, final_time),
        initial_state,
        method="RK45",
        max_step=max_step,
        rtol=1.0e-7,
        atol=1.0e-9,
    )
    if not solution.success:
        raise RuntimeError(solution.message)
    return solution.t, solution.y


def postprocess(
    time: np.ndarray,
    state_history: np.ndarray,
) -> dict[str, np.ndarray]:
    xp = state_history[0:2]
    y = P.kp * (xp[0] + P.zp * xp[1])
    ym = state_history[2]
    theta = state_history[10:14]
    r = np.array([reference(t) for t in time])

    u = np.empty_like(time)
    augmented_error = np.empty_like(time)
    for k in range(time.size):
        omega = np.array(
            [
                state_history[3, k],
                state_history[4, k],
                y[k],
                r[k],
            ]
        )
        zeta = state_history[5:9, k]
        v = state_history[9, k]
        u[k] = theta[:, k] @ omega
        augmented_error[k] = y[k] - ym[k] + theta[:, k] @ zeta - v

    return {
        "r": r,
        "y": y,
        "ym": ym,
        "u": u,
        "e": y - ym,
        "ea": augmented_error,
        "theta": theta,
    }


def save_csv(
    filename: str,
    time: np.ndarray,
    signals: dict[str, np.ndarray],
) -> None:
    matrix = np.column_stack(
        (
            time,
            signals["r"],
            signals["y"],
            signals["ym"],
            signals["u"],
            signals["e"],
            signals["ea"],
            signals["theta"].T,
        )
    )
    header = "t,r,y,ym,u,e,e_aug,theta1,theta2,theta3,theta4"
    np.savetxt(filename, matrix, delimiter=",", header=header, comments="")


def plot_results(
    time: np.ndarray,
    signals: dict[str, np.ndarray],
) -> None:
    plt.figure()
    plt.plot(time, signals["y"], label="plant output y")
    plt.plot(time, signals["ym"], "--", label="reference output ym")
    plt.plot(time, signals["r"], ":", label="command r")
    plt.xlabel("Time [s]")
    plt.ylabel("Output")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(time, signals["e"], label="tracking error")
    plt.plot(time, signals["ea"], "--", label="augmented error")
    plt.xlabel("Time [s]")
    plt.ylabel("Error")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    for index in range(4):
        plt.plot(time, signals["theta"][index], label=f"theta{index + 1}")
    plt.xlabel("Time [s]")
    plt.ylabel("Adaptive parameter")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.show()


def main() -> None:
    time, states = run_simulation()
    signals = postprocess(time, states)
    save_csv("Chapter15_Lesson4_results.csv", time, signals)

    rms_error = float(np.sqrt(np.mean(signals["e"] ** 2)))
    final_error = float(signals["e"][-1])
    print(f"RMS tracking error: {rms_error:.6f}")
    print(f"Final tracking error: {final_error:.6e}")
    print("Final adaptive parameters:", signals["theta"][:, -1])

    plot_results(time, signals)


if __name__ == "__main__":
    main()

The C++ version uses only the C++17 standard library. A fixed-step fourth-order Runge-Kutta method is implemented from scratch, and the results are written to CSV.

Chapter15_Lesson4.cpp

/*
Chapter15_Lesson4.cpp
From-scratch fixed-step RK4 simulation of adaptive output-feedback MRAC.

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

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

namespace {

constexpr std::size_t N = 14;
using State = std::array<double, N>;

constexpr double A1 = 3.0;
constexpr double A0 = 2.0;
constexpr double KP = 2.0;
constexpr double ZP = 1.5;
constexpr double FILTER_POLE = 4.0;
constexpr double AM = 1.2;
constexpr double KM = 1.2;
constexpr std::array<double, 4> GAMMA{8.0, 8.0, 5.0, 5.0};

double reference(double t) {
    if (t < 15.0) {
        return 1.0;
    }
    if (t < 30.0) {
        return -0.5;
    }
    return 0.8;
}

double dot4(
    const std::array<double, 4>& a,
    const std::array<double, 4>& b
) {
    double result = 0.0;
    for (std::size_t i = 0; i < 4; ++i) {
        result += a[i] * b[i];
    }
    return result;
}

State derivative(double t, const State& x) {
    const double xp0 = x[0];
    const double xp1 = x[1];
    const double ym = x[2];
    const double omegaU = x[3];
    const double omegaY = x[4];
    const std::array<double, 4> zeta{x[5], x[6], x[7], x[8]};
    const double v = x[9];
    const std::array<double, 4> theta{x[10], x[11], x[12], x[13]};

    const double r = reference(t);
    const double y = KP * (xp0 + ZP * xp1);
    const std::array<double, 4> omega{omegaU, omegaY, y, r};
    const double u = dot4(theta, omega);

    const double e = y - ym;
    const double eAug = e + dot4(theta, zeta) - v;
    const double normalization = 1.0 + dot4(omega, omega) + dot4(zeta, zeta);

    State dx{};
    dx[0] = -A1 * xp0 - A0 * xp1 + u;
    dx[1] = xp0;
    dx[2] = -AM * ym + KM * r;
    dx[3] = -FILTER_POLE * omegaU + u;
    dx[4] = -FILTER_POLE * omegaY + y;

    for (std::size_t i = 0; i < 4; ++i) {
        dx[5 + i] = -AM * zeta[i] + KM * omega[i];
    }

    dx[9] = -AM * v + KM * u;
    for (std::size_t i = 0; i < 4; ++i) {
        dx[10 + i] = -GAMMA[i] * zeta[i] * eAug / normalization;
    }
    return dx;
}

State addScaled(const State& x, const State& k, double scale) {
    State result{};
    for (std::size_t i = 0; i < N; ++i) {
        result[i] = x[i] + scale * k[i];
    }
    return result;
}

State rk4Step(double t, const State& x, double dt) {
    const State k1 = derivative(t, x);
    const State k2 = derivative(t + 0.5 * dt, addScaled(x, k1, 0.5 * dt));
    const State k3 = derivative(t + 0.5 * dt, addScaled(x, k2, 0.5 * dt));
    const State k4 = derivative(t + dt, addScaled(x, k3, dt));

    State next{};
    for (std::size_t i = 0; i < N; ++i) {
        next[i] = x[i] + (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.");
        }
    }
    return next;
}

void writeRow(std::ofstream& out, double t, const State& x) {
    const double r = reference(t);
    const double y = KP * (x[0] + ZP * x[1]);
    const double ym = x[2];
    const std::array<double, 4> omega{x[3], x[4], y, r};
    const std::array<double, 4> zeta{x[5], x[6], x[7], x[8]};
    const std::array<double, 4> theta{x[10], x[11], x[12], x[13]};
    const double u = dot4(theta, omega);
    const double e = y - ym;
    const double eAug = e + dot4(theta, zeta) - x[9];

    out << std::setprecision(12)
        << t << ',' << r << ',' << y << ',' << ym << ',' << u << ','
        << e << ',' << eAug << ','
        << theta[0] << ',' << theta[1] << ','
        << theta[2] << ',' << theta[3] << '\n';
}

}  // namespace

int main() {
    constexpr double DT = 0.001;
    constexpr double FINAL_TIME = 45.0;
    const std::size_t steps = static_cast<std::size_t>(FINAL_TIME / DT);

    State x{};
    x[10] = 0.0;
    x[11] = 0.0;
    x[12] = -0.5;
    x[13] = 0.5;

    std::ofstream out("Chapter15_Lesson4_results_cpp.csv");
    if (!out) {
        std::cerr << "Cannot open output CSV file.\n";
        return 1;
    }
    out << "t,r,y,ym,u,e,e_aug,theta1,theta2,theta3,theta4\n";

    double squaredErrorSum = 0.0;
    for (std::size_t k = 0; k <= steps; ++k) {
        const double t = static_cast<double>(k) * DT;
        writeRow(out, t, x);

        const double y = KP * (x[0] + ZP * x[1]);
        const double e = y - x[2];
        squaredErrorSum += e * e;

        if (k < steps) {
            x = rk4Step(t, x, DT);
        }
    }

    const double rmsError = std::sqrt(
        squaredErrorSum / static_cast<double>(steps + 1)
    );
    const double finalY = KP * (x[0] + ZP * x[1]);
    std::cout << std::fixed << std::setprecision(8);
    std::cout << "RMS tracking error: " << rmsError << '\n';
    std::cout << "Final tracking error: " << finalY - x[2] << '\n';
    std::cout << "Final adaptive parameters: "
              << x[10] << ", " << x[11] << ", "
              << x[12] << ", " << x[13] << '\n';
    return 0;
}

The Java implementation also uses a from-scratch RK4 integrator and the standard I/O library. It therefore has no third-party dependency.

Chapter15_Lesson4.java

/*
Chapter15_Lesson4.java
From-scratch RK4 simulation of adaptive output-feedback MRAC.

Compile:
    javac Chapter15_Lesson4.java
Run:
    java Chapter15_Lesson4
*/

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

public final class Chapter15_Lesson4 {
    private static final int N = 14;

    private static final double A1 = 3.0;
    private static final double A0 = 2.0;
    private static final double KP = 2.0;
    private static final double ZP = 1.5;
    private static final double FILTER_POLE = 4.0;
    private static final double AM = 1.2;
    private static final double KM = 1.2;
    private static final double[] GAMMA = {8.0, 8.0, 5.0, 5.0};

    private Chapter15_Lesson4() {
    }

    private static double reference(double t) {
        if (t < 15.0) {
            return 1.0;
        }
        if (t < 30.0) {
            return -0.5;
        }
        return 0.8;
    }

    private static double dot4(double[] a, double[] b) {
        double result = 0.0;
        for (int i = 0; i < 4; i++) {
            result += a[i] * b[i];
        }
        return result;
    }

    private static double[] derivative(double t, double[] x) {
        double xp0 = x[0];
        double xp1 = x[1];
        double ym = x[2];
        double omegaU = x[3];
        double omegaY = x[4];

        double[] zeta = {x[5], x[6], x[7], x[8]};
        double v = x[9];
        double[] theta = {x[10], x[11], x[12], x[13]};

        double r = reference(t);
        double y = KP * (xp0 + ZP * xp1);
        double[] omega = {omegaU, omegaY, y, r};
        double u = dot4(theta, omega);

        double trackingError = y - ym;
        double augmentedError = trackingError + dot4(theta, zeta) - v;
        double normalization = 1.0 + dot4(omega, omega) + dot4(zeta, zeta);

        double[] dx = new double[N];
        dx[0] = -A1 * xp0 - A0 * xp1 + u;
        dx[1] = xp0;
        dx[2] = -AM * ym + KM * r;
        dx[3] = -FILTER_POLE * omegaU + u;
        dx[4] = -FILTER_POLE * omegaY + y;

        for (int i = 0; i < 4; i++) {
            dx[5 + i] = -AM * zeta[i] + KM * omega[i];
        }

        dx[9] = -AM * v + KM * u;
        for (int i = 0; i < 4; i++) {
            dx[10 + i] = -GAMMA[i] * zeta[i] * augmentedError / normalization;
        }
        return dx;
    }

    private static double[] addScaled(double[] x, double[] k, double scale) {
        double[] result = new double[N];
        for (int i = 0; i < N; i++) {
            result[i] = x[i] + scale * k[i];
        }
        return result;
    }

    private static double[] rk4Step(double t, double[] x, double dt) {
        double[] k1 = derivative(t, x);
        double[] k2 = derivative(t + 0.5 * dt, addScaled(x, k1, 0.5 * dt));
        double[] k3 = derivative(t + 0.5 * dt, addScaled(x, k2, 0.5 * dt));
        double[] k4 = derivative(t + dt, addScaled(x, k3, dt));

        double[] next = new double[N];
        for (int i = 0; i < N; i++) {
            next[i] = x[i] + (dt / 6.0)
                * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(next[i])) {
                throw new IllegalStateException("Non-finite state encountered.");
            }
        }
        return next;
    }

    private static void writeRow(
        BufferedWriter writer,
        double t,
        double[] x
    ) throws IOException {
        double r = reference(t);
        double y = KP * (x[0] + ZP * x[1]);
        double ym = x[2];

        double[] omega = {x[3], x[4], y, r};
        double[] zeta = {x[5], x[6], x[7], x[8]};
        double[] theta = {x[10], x[11], x[12], x[13]};

        double u = dot4(theta, omega);
        double e = y - ym;
        double eAug = e + dot4(theta, zeta) - x[9];

        writer.write(String.format(
            Locale.US,
            "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,"
                + "%.12f,%.12f,%.12f,%.12f%n",
            t, r, y, ym, u, e, eAug,
            theta[0], theta[1], theta[2], theta[3]
        ));
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);

        final double dt = 0.001;
        final double finalTime = 45.0;
        final int steps = (int) Math.round(finalTime / dt);

        double[] x = new double[N];
        x[10] = 0.0;
        x[11] = 0.0;
        x[12] = -0.5;
        x[13] = 0.5;

        double squaredErrorSum = 0.0;

        try (BufferedWriter writer = new BufferedWriter(
            new FileWriter("Chapter15_Lesson4_results_java.csv")
        )) {
            writer.write("t,r,y,ym,u,e,e_aug,theta1,theta2,theta3,theta4\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * dt;
                writeRow(writer, t, x);

                double y = KP * (x[0] + ZP * x[1]);
                double e = y - x[2];
                squaredErrorSum += e * e;

                if (k < steps) {
                    x = rk4Step(t, x, dt);
                }
            }
        } catch (IOException ex) {
            System.err.println("File error: " + ex.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(squaredErrorSum / (steps + 1.0));
        double finalY = KP * (x[0] + ZP * x[1]);

        System.out.printf("RMS tracking error: %.8f%n", rmsError);
        System.out.printf("Final tracking error: %.8f%n", finalY - x[2]);
        System.out.printf(
            "Final adaptive parameters: %.8f, %.8f, %.8f, %.8f%n",
            x[10], x[11], x[12], x[13]
        );
    }
}

The first file is an executable ode45 simulation. The second file programmatically constructs a Simulink architecture showing the reference model, unknown plant, measured-signal filters, augmented-error generator, normalized adaptation law, and parameterized control law.

Chapter15_Lesson4.m

% Chapter15_Lesson4.m
% Adaptive output-feedback MRAC with filtered regressors and augmented error.
% The simulation uses ODE45 and only u and y in the adaptive controller.

clear; close all; clc;

p.a1 = 3.0;
p.a0 = 2.0;
p.kp = 2.0;
p.zp = 1.5;
p.lambda = 4.0;
p.am = 1.2;
p.km = 1.2;
p.Gamma = diag([8.0, 8.0, 5.0, 5.0]);

x0 = zeros(14, 1);
x0(11:14) = [0.0; 0.0; -0.5; 0.5];

options = odeset('RelTol', 1e-7, 'AbsTol', 1e-9, 'MaxStep', 0.01);
[t, x] = ode45(@(time, state) closedLoopDynamics(time, state, p), ...
               [0.0, 45.0], x0, options);

n = numel(t);
r = zeros(n, 1);
y = zeros(n, 1);
ym = x(:, 3);
u = zeros(n, 1);
eAug = zeros(n, 1);

for k = 1:n
    r(k) = commandSignal(t(k));
    y(k) = p.kp * (x(k, 1) + p.zp * x(k, 2));

    omega = [x(k, 4); x(k, 5); y(k); r(k)];
    zeta = x(k, 6:9).';
    theta = x(k, 11:14).';

    u(k) = theta.' * omega;
    eAug(k) = y(k) - ym(k) + theta.' * zeta - x(k, 10);
end

e = y - ym;
rmsError = sqrt(mean(e.^2));

fprintf('RMS tracking error: %.8f\n', rmsError);
fprintf('Final tracking error: %.8e\n', e(end));
fprintf('Final adaptive parameters:\n');
disp(x(end, 11:14));

results = table(t, r, y, ym, u, e, eAug, ...
    x(:, 11), x(:, 12), x(:, 13), x(:, 14), ...
    'VariableNames', {'t', 'r', 'y', 'ym', 'u', 'e', 'e_aug', ...
                      'theta1', 'theta2', 'theta3', 'theta4'});
writetable(results, 'Chapter15_Lesson4_results_matlab.csv');

figure;
plot(t, y, 'LineWidth', 1.2); hold on;
plot(t, ym, '--', 'LineWidth', 1.2);
plot(t, r, ':', 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Output');
legend('Plant output y', 'Reference output y_m', 'Command r', ...
       'Location', 'best');
title('Adaptive Output-Feedback MRAC');

figure;
plot(t, e, 'LineWidth', 1.2); hold on;
plot(t, eAug, '--', 'LineWidth', 1.2);
grid on;
xlabel('Time [s]');
ylabel('Error');
legend('Tracking error', 'Augmented error', 'Location', 'best');

figure;
plot(t, x(:, 11:14), 'LineWidth', 1.1);
grid on;
xlabel('Time [s]');
ylabel('Adaptive parameter');
legend('\theta_1', '\theta_2', '\theta_3', '\theta_4', ...
       'Location', 'best');

function dx = closedLoopDynamics(t, x, p)
    xp = x(1:2);
    ym = x(3);
    omegaU = x(4);
    omegaY = x(5);
    zeta = x(6:9);
    v = x(10);
    theta = x(11:14);

    r = commandSignal(t);
    y = p.kp * (xp(1) + p.zp * xp(2));
    omega = [omegaU; omegaY; y; r];
    u = theta.' * omega;

    trackingError = y - ym;
    augmentedError = trackingError + theta.' * zeta - v;
    normalization = 1.0 + omega.' * omega + zeta.' * zeta;

    dx = zeros(14, 1);
    dx(1) = -p.a1 * xp(1) - p.a0 * xp(2) + u;
    dx(2) = xp(1);
    dx(3) = -p.am * ym + p.km * r;
    dx(4) = -p.lambda * omegaU + u;
    dx(5) = -p.lambda * omegaY + y;
    dx(6:9) = -p.am * zeta + p.km * omega;
    dx(10) = -p.am * v + p.km * u;
    dx(11:14) = -p.Gamma * zeta * augmentedError / normalization;
end

function r = commandSignal(t)
    if t < 15.0
        r = 1.0;
    elseif t < 30.0
        r = -0.5;
    else
        r = 0.8;
    end
end

Chapter15_Lesson4_Simulink.m

% Chapter15_Lesson4_Simulink.m
% Programmatically builds a conceptual Simulink architecture for the lesson.
% The numerical implementation is in Chapter15_Lesson4.m.

model = 'Chapter15_Lesson4_OutputFeedbackMRAC';

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

new_system(model);
open_system(model);
set_param(model, 'StopTime', '45');

add_block('simulink/Sources/Step', [model '/Reference r'], ...
    'Position', [30 70 60 100], ...
    'Time', '0', 'Before', '0', 'After', '1');

add_block('simulink/Continuous/Transfer Fcn', [model '/Reference Model Wm'], ...
    'Position', [130 45 260 105], ...
    'Numerator', '[1.2]', ...
    'Denominator', '[1 1.2]');

add_block('simulink/Ports & Subsystems/Subsystem', ...
    [model '/Adaptive Controller'], ...
    'Position', [130 170 310 270]);

add_block('simulink/Continuous/State-Space', [model '/Unknown Plant'], ...
    'Position', [390 170 540 270], ...
    'A', '[-3 -2; 1 0]', ...
    'B', '[1; 0]', ...
    'C', '[2 3]', ...
    'D', '0');

add_block('simulink/Math Operations/Sum', [model '/Tracking Error'], ...
    'Position', [610 80 640 130], ...
    'Inputs', '+-');

add_block('simulink/Sinks/Scope', [model '/Outputs Scope'], ...
    'Position', [720 55 760 125], ...
    'NumInputPorts', '3');

add_block('simulink/Sinks/Scope', [model '/Control Scope'], ...
    'Position', [720 205 760 255]);

add_line(model, 'Reference r/1', 'Reference Model Wm/1');
add_line(model, 'Reference r/1', 'Adaptive Controller/1');
add_line(model, 'Adaptive Controller/1', 'Unknown Plant/1');
add_line(model, 'Adaptive Controller/1', 'Control Scope/1');
add_line(model, 'Unknown Plant/1', 'Adaptive Controller/2');
add_line(model, 'Reference Model Wm/1', 'Tracking Error/1');
add_line(model, 'Unknown Plant/1', 'Tracking Error/2');
add_line(model, 'Reference r/1', 'Outputs Scope/1');
add_line(model, 'Reference Model Wm/1', 'Outputs Scope/2');
add_line(model, 'Unknown Plant/1', 'Outputs Scope/3');

controller = [model '/Adaptive Controller'];
Simulink.SubSystem.deleteContents(controller);

add_block('simulink/Ports & Subsystems/In1', [controller '/r'], ...
    'Position', [25 40 55 60]);
add_block('simulink/Ports & Subsystems/In1', [controller '/y'], ...
    'Position', [25 105 55 125], 'Port', '2');
add_block('simulink/Ports & Subsystems/Out1', [controller '/u'], ...
    'Position', [430 70 460 90]);

add_block('simulink/Ports & Subsystems/Subsystem', ...
    [controller '/Stable Regressor Filters'], ...
    'Position', [100 25 245 85]);
add_block('simulink/Ports & Subsystems/Subsystem', ...
    [controller '/Augmented Error Generator'], ...
    'Position', [100 110 265 170]);
add_block('simulink/Ports & Subsystems/Subsystem', ...
    [controller '/Normalized Adaptive Law'], ...
    'Position', [300 105 405 165]);
add_block('simulink/Ports & Subsystems/Subsystem', ...
    [controller '/Parameterised Control Law'], ...
    'Position', [285 25 405 85]);

add_line(controller, 'r/1', 'Stable Regressor Filters/1');
add_line(controller, 'y/1', 'Stable Regressor Filters/2');
add_line(controller, 'Stable Regressor Filters/1', ...
    'Parameterised Control Law/1');
add_line(controller, 'Stable Regressor Filters/1', ...
    'Augmented Error Generator/1');
add_line(controller, 'Augmented Error Generator/1', ...
    'Normalized Adaptive Law/1');
add_line(controller, 'Normalized Adaptive Law/1', ...
    'Parameterised Control Law/2');
add_line(controller, 'Parameterised Control Law/1', 'u/1');

set_param([controller '/Stable Regressor Filters'], ...
    'AttributesFormatString', ...
    'omega_u = 1/(s+lambda) u\nomega_y = 1/(s+lambda) y');
set_param([controller '/Augmented Error Generator'], ...
    'AttributesFormatString', ...
    'e_a = e + theta^T zeta - Wm[u]');
set_param([controller '/Normalized Adaptive Law'], ...
    'AttributesFormatString', ...
    'theta_dot = -Gamma zeta e_a / m^2');
set_param([controller '/Parameterised Control Law'], ...
    'AttributesFormatString', ...
    'u = theta^T omega');

Simulink.BlockDiagram.arrangeSystem(model);
save_system(model);
open_system(model);

fprintf('Created %s.slx\n', model);
fprintf(['Use Chapter15_Lesson4.m for the executable numerical model; ' ...
         'this generated diagram documents the signal architecture.\n']);

The Wolfram Language implementation uses list-based RK4 integration, exports CSV data, and produces tracking, error, and parameter plots. The notebook file loads the accompanying .wl implementation from the same directory.

Chapter15_Lesson4.wl

(* Chapter15_Lesson4.wl
   Adaptive output-feedback MRAC with filtered regressors and augmented error.
   Run with Wolfram Mathematica or wolframscript.
*)

ClearAll["Global`*"];

a1 = 3.0;
a0 = 2.0;
kp = 2.0;
zp = 1.5;
filterPole = 4.0;
am = 1.2;
km = 1.2;
gamma = {8.0, 8.0, 5.0, 5.0};

reference[t_?NumericQ] := Piecewise[{
    {1.0, t < 15.0},
    {-0.5, t < 30.0}
}, 0.8];

dot4[a_List, b_List] := Total[a b];

derivative[t_?NumericQ, x_List] := Module[
    {xp0, xp1, ym, omegaU, omegaY, zeta, v, theta,
     r, y, omega, u, trackingError, augmentedError,
     normalization, dx},

    xp0 = x[[1]];
    xp1 = x[[2]];
    ym = x[[3]];
    omegaU = x[[4]];
    omegaY = x[[5]];
    zeta = x[[6 ;; 9]];
    v = x[[10]];
    theta = x[[11 ;; 14]];

    r = reference[t];
    y = kp (xp0 + zp xp1);
    omega = {omegaU, omegaY, y, r};
    u = dot4[theta, omega];

    trackingError = y - ym;
    augmentedError = trackingError + dot4[theta, zeta] - v;
    normalization = 1.0 + dot4[omega, omega] + dot4[zeta, zeta];

    dx = ConstantArray[0.0, 14];
    dx[[1]] = -a1 xp0 - a0 xp1 + u;
    dx[[2]] = xp0;
    dx[[3]] = -am ym + km r;
    dx[[4]] = -filterPole omegaU + u;
    dx[[5]] = -filterPole omegaY + y;
    dx[[6 ;; 9]] = -am zeta + km omega;
    dx[[10]] = -am v + km u;
    dx[[11 ;; 14]] = -gamma zeta augmentedError/normalization;
    dx
];

rk4Step[t_?NumericQ, x_List, dt_?NumericQ] := Module[
    {k1, k2, k3, k4},
    k1 = derivative[t, x];
    k2 = derivative[t + dt/2.0, x + dt k1/2.0];
    k3 = derivative[t + dt/2.0, x + dt k2/2.0];
    k4 = derivative[t + dt, x + dt k3];
    x + dt (k1 + 2.0 k2 + 2.0 k3 + k4)/6.0
];

dt = 0.002;
finalTime = 45.0;
steps = Round[finalTime/dt];

initialState = ConstantArray[0.0, 14];
initialState[[11 ;; 14]] = {0.0, 0.0, -0.5, 0.5};

trajectory = Reap[
    Module[{state = initialState, t = 0.0},
        Sow[Prepend[state, t]];
        Do[
            state = rk4Step[t, state, dt];
            t = k dt;
            Sow[Prepend[state, t]],
            {k, 1, steps}
        ];
    ]
][[2, 1]];

processed = Map[
    Function[row,
        Module[{t, x, r, y, ym, omega, zeta, theta, u, e, eAug},
            t = row[[1]];
            x = row[[2 ;;]];
            r = reference[t];
            y = kp (x[[1]] + zp x[[2]]);
            ym = x[[3]];
            omega = {x[[4]], x[[5]], y, r};
            zeta = x[[6 ;; 9]];
            theta = x[[11 ;; 14]];
            u = dot4[theta, omega];
            e = y - ym;
            eAug = e + dot4[theta, zeta] - x[[10]];
            {t, r, y, ym, u, e, eAug,
             theta[[1]], theta[[2]], theta[[3]], theta[[4]]}
        ]
    ],
    trajectory
];

header = {"t", "r", "y", "ym", "u", "e", "e_aug",
          "theta1", "theta2", "theta3", "theta4"};

Export[
    "Chapter15_Lesson4_results_wolfram.csv",
    Prepend[processed, header]
];

rmsError = Sqrt[Mean[processed[[All, 6]]^2]];
Print["RMS tracking error: ", NumberForm[rmsError, {10, 8}]];
Print["Final tracking error: ",
      ScientificForm[processed[[-1, 6]], 6]];
Print["Final adaptive parameters: ", processed[[-1, 8 ;; 11]]];

trackingPlot = ListLinePlot[
    {
        processed[[All, {1, 3}]],
        processed[[All, {1, 4}]],
        processed[[All, {1, 2}]]
    },
    PlotLegends -> {"Plant output y", "Reference output ym", "Command r"},
    Frame -> True,
    FrameLabel -> {"Time [s]", "Output"},
    PlotRange -> All,
    ImageSize -> Large
];

errorPlot = ListLinePlot[
    {
        processed[[All, {1, 6}]],
        processed[[All, {1, 7}]]
    },
    PlotLegends -> {"Tracking error", "Augmented error"},
    Frame -> True,
    FrameLabel -> {"Time [s]", "Error"},
    PlotRange -> All,
    ImageSize -> Large
];

parameterPlot = ListLinePlot[
    Table[processed[[All, {1, 7 + i}]], {i, 1, 4}],
    PlotLegends -> {"theta1", "theta2", "theta3", "theta4"},
    Frame -> True,
    FrameLabel -> {"Time [s]", "Adaptive parameter"},
    PlotRange -> All,
    ImageSize -> Large
];

Print[trackingPlot];
Print[errorPlot];
Print[parameterPlot];

Chapter15_Lesson4.nb

Notebook[{
  Cell["Chapter 15, Lesson 4: Adaptive Output-Feedback MRAC", "Title"],
  Cell[
    "This notebook loads and executes the accompanying Wolfram Language implementation.",
    "Text"
  ],
  Cell[
    BoxData[
      RowBox[{
        "Get", "[",
        RowBox[{
          "FileNameJoin", "[",
          RowBox[{"{",
            RowBox[{
              RowBox[{"NotebookDirectory", "[", "]"}],
              ",",
              "\"Chapter15_Lesson4.wl\""
            }],
          "}"}],
          "]"
        }],
        "]"
      }]
    ],
    "Input"
  ]
},
WindowTitle -> "Chapter15_Lesson4",
StyleDefinitions -> "Default.nb"
]

The simulated plant is realized by

\[ \dot{x}_p= \begin{bmatrix} -3 & -2\\ 1 & 0 \end{bmatrix}x_p+ \begin{bmatrix} 1\\0 \end{bmatrix}u, \qquad y= \begin{bmatrix} 2 & 3 \end{bmatrix}x_p, \]

which has transfer function \(2(s+1.5)/(s^2+3s+2)\). The reference model is \(1.2/(s+1.2)\). For each regressor component,

\[ \dot{\boldsymbol{\zeta}} =-1.2\boldsymbol{\zeta} +1.2\boldsymbol{\omega}, \qquad \dot{v}=-1.2v+1.2u. \]

With the chosen initial controller parameters, all three independently executable implementations produce an RMS tracking error of approximately \(0.0766\) over 45 seconds and a final tracking error near \(2.3\times10^{-7}\). The final parameter vector is not the ideal vector because the command does not provide persistent excitation of all regressor directions.

16.1 What should be monitored in practice?

• tracking error and augmented error separately;
• the norm of \(\boldsymbol{\theta}\) and proximity to projection bounds;
• the normalization signal \(m^2\);
• filter-state magnitudes;
• control magnitude, rate, and saturation duration;
• noise sensitivity when filter poles or adaptation gains are increased.

1. Moving time-varying parameters through a filter: writing \(W(s)[\theta^{T}(t)\omega(t)]=\theta^{T}(t)W(s)[\omega(t)]\) is generally false. The augmented error compensates for this swapping term.
2. Ignoring unstable zeros: exact output model matching may demand cancellation of nonminimum-phase zeros, which is not a stable operation.
3. Assuming separation: an observer error system coupled to adaptive parameters is not independent of the controller.
4. Equating tracking with identification: \(e\to0\) does not imply \(\widetilde{\theta}\to0\) without excitation.
5. Using very fast filters: large filter bandwidth reduces phase lag but can amplify measurement noise and numerical stiffness.
6. Ignoring actuator saturation: saturation breaks the nominal relation \(u=\theta^{T}\omega\) and can create continued adaptation against an unattainable command.

Problem 1 (Filtered-regressor dimension): For an order-\(n\), relative-degree-one SISO plant, show that the regressor in Section 3 has dimension \(2n\).

Solution:

Each vector \(\boldsymbol{\omega}_1\) and \(\boldsymbol{\omega}_2\) contains \(n-1\) components because \(\boldsymbol{\alpha}(s)\) contains the powers \(s^{n-2},\ldots,1\). Adding the direct output and command components gives

\[ \dim(\boldsymbol{\omega}) =(n-1)+(n-1)+1+1=2n. \]

Problem 2 (SPR Lyapunov cancellation): Starting from the error realization in Section 6, prove that the update \(\dot{\boldsymbol{\theta}}=-\Gamma\operatorname{sgn}(k_p) \boldsymbol{\omega}e\) cancels the mixed state-parameter term in \(\dot V\).

Solution:

From \(Pb_e=c_e\), the state cross term is

\[ 2x_e^{T}Pb_ek_p \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} = 2k_pe\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}. \]

The parameter-energy derivative contributes

\[ 2|k_p| \widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\boldsymbol{\theta}} = -2|k_p|\operatorname{sgn}(k_p)e\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega} = -2k_pe\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\omega}. \]

The terms cancel exactly, leaving \(\dot V=-x_e^{T}Qx_e\le0\).

Problem 3 (Augmented-error identity): Let \(u=\boldsymbol{\theta}^{T}\boldsymbol{\omega}\), \(\boldsymbol{\zeta}=W_m(s)[\boldsymbol{\omega}]\), and \(v=W_m(s)[u]\). Prove that \(e_a=e+\boldsymbol{\theta}^{T}\boldsymbol{\zeta}-v =\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\zeta}\).

Solution:

Use the model-matching error relation \(e=W_m(s)[\boldsymbol{\theta}^{T}\boldsymbol{\omega}] -W_m(s)[{\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\omega}]\). Since \(\boldsymbol{\theta}^{\star}\) is constant,

\[ W_m(s)[{\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\omega}] = {\boldsymbol{\theta}^{\star}}^{T} W_m(s)[\boldsymbol{\omega}] = {\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\zeta}. \]

Substitution gives

\[ \begin{aligned} e_a &= v-{\boldsymbol{\theta}^{\star}}^{T}\boldsymbol{\zeta} +\boldsymbol{\theta}^{T}\boldsymbol{\zeta}-v\\ &= (\boldsymbol{\theta}-\boldsymbol{\theta}^{\star})^{T} \boldsymbol{\zeta} = \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\zeta}. \end{aligned} \]

Problem 4 (Numerical model matching): For the numerical values in Section 10, solve the three polynomial matching equations and find \(\boldsymbol{\theta}^{\star}\).

Solution:

The target polynomial is

\[ (s+1.5)(s+4)(s+1.2) = s^3+6.7s^2+12.6s+7.2. \]

Expanding the controller-dependent denominator and equating coefficients produces

\[ \begin{aligned} \theta_1+2\theta_3&=0.3,\\ 3\theta_1+2\theta_2+11\theta_3&=1.4,\\ 2\theta_1+3\theta_2+12\theta_3&=0.8. \end{aligned} \]

Solving gives

\[ \theta_1^{\star}=2.5,\qquad \theta_2^{\star}=3.0,\qquad \theta_3^{\star}=-1.1, \qquad \theta_4^{\star}=\frac{1.2}{2}=0.6. \]

Problem 5 (Tracking without parameter convergence): Explain mathematically why the tracking error may converge to zero while the adaptive parameter vector converges to a value different from \(\boldsymbol{\theta}^{\star}\).

Solution:

The augmented regression error is \(e_a=\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\zeta}\). If the regressor explores only a lower-dimensional subspace, there may be a nonzero vector \(\widetilde{\boldsymbol{\theta}}_\infty\) satisfying

\[ \widetilde{\boldsymbol{\theta}}_\infty^{T} \boldsymbol{\zeta}(t)=0 \qquad \text{along the limiting trajectory}. \]

Then the prediction or augmented error can vanish without the full parameter error vanishing. Persistent excitation removes this ambiguity by requiring the regressor Gramian to satisfy

\[ \int_t^{t+T} \boldsymbol{\zeta}(\tau) \boldsymbol{\zeta}^{T}(\tau)\,d\tau \ge \alpha I, \qquad T>0,\quad \alpha>0. \]

Under the additional regularity conditions of adaptive-system convergence theory, this full-rank excitation can imply \(\widetilde{\boldsymbol{\theta}}\to0\).

Problem 6 (Effect of normalization): Suppose all regressor signals are multiplied by a large scalar \(c>1\). Compare the unnormalized update \(\dot{\theta}=-\Gamma\zeta e_a\) with the normalized update of Section 8.

Solution:

Since \(e_a=\widetilde{\theta}^{T}\zeta\), scaling \(\zeta\) by \(c\) scales \(e_a\) approximately by \(c\). The unnormalized parameter rate therefore scales approximately as \(c^2\), which can create stiffness and sensitivity to noise. For large signals, \(m^2\) scales as \(c^2\), so

\[ \frac{\zeta e_a}{m^2} \]

remains of comparable magnitude. Normalization therefore changes the effective adaptation gain according to signal amplitude.

Adaptive output-feedback MRAC can be constructed without estimating the complete plant state. Stable filters of measured input and output signals generate a nonminimal regressor, and a direct controller adapts its gains to satisfy a model-matching objective. When the resulting error transfer function is SPR, a KYP-based Lyapunov proof gives boundedness and asymptotic tracking. More general filtered relations require an augmented error because time-varying parameters do not commute with dynamic filters. Normalization and projection improve implementability, but exact tracking and parameter convergence remain distinct properties. The next lesson studies the stability challenges that arise when these structural assumptions, boundedness arguments, or ideal matching conditions are weakened.