After completing this lesson, a student should be able to:

• define state-tracking, output-tracking, and parameter-estimation errors;
• derive an MRAC error model from plant, controller, and reference-model equations;
• identify the homogeneous stable dynamics and the parameter-mismatch forcing term;
• derive the parameter-error dynamics for constant ideal parameters;
• construct a state–parameter Lyapunov candidate and perform the cross-term cancellation proof;
• explain rigorously why \( e(t) → 0 \) does not generally imply \( \widetilde{\boldsymbol{\theta} }(t) → \mathbf{0} \).

The analysis uses full-state notation and a known input direction so that the central error-model geometry remains visible. Direct versus indirect MRAC, scalar first-order design, higher-order canonical forms, robustness modifications, output feedback, and formal excitation conditions are treated in later chapters.

Let \( \mathbf{x}(t)\in\mathbb{R}^{n} \) be the plant state and \( \mathbf{x}_m(t)\in\mathbb{R}^{n} \) the state of the reference model introduced in Lesson 1. The state-tracking error is defined as

\[ \mathbf{e}(t)=\mathbf{x}(t)-\mathbf{x}_m(t). \]

If the measured or regulated output is \( y=\mathbf{c}^{T}\mathbf{x} \) and the reference output is \( y_m=\mathbf{c}^{T}\mathbf{x}_m \), then

\[ e_y(t)=y(t)-y_m(t)=\mathbf{c}^{T}\mathbf{e}(t). \]

Tracking error is a difference between trajectories. Parameter error is a difference between an adjustable parameter and the ideal parameter that would make exact model matching possible:

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

The sign convention above is used throughout this lesson. Some books define the reverse sign; either convention is valid, but the adaptive-law sign must be changed consistently.

Differentiation gives the general parameter-error dynamics

\[ \dot{\widetilde{\boldsymbol{\theta} } } =\dot{\widehat{\boldsymbol{\theta} } } -\dot{\boldsymbol{\theta} }^{*}. \]

Under the standard introductory MRAC assumption that the unknown plant parameters and therefore the ideal controller parameters are constant,

\[ \dot{\boldsymbol{\theta} }^{*}=\mathbf{0}, \qquad \dot{\widetilde{\boldsymbol{\theta} } } =\dot{\widehat{\boldsymbol{\theta} } }. \]

Thus an adaptive law is simultaneously a differential equation for the controller estimate and, under constant ideal parameters, for the parameter error.

Consider a single-input linear plant

\[ \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{b}u, \qquad \mathbf{x}\in\mathbb{R}^{n}, \quad u\in\mathbb{R}, \]

and the desired reference model

\[ \dot{\mathbf{x} }_m =\mathbf{A}_m\mathbf{x}_m+\mathbf{b}_m r, \]

where \( \mathbf{A}_m \) is Hurwitz and \( r(t) \) is a bounded command. Use an adjustable state-feedback and command-feedforward controller:

\[ u =\widehat{\mathbf{k} }_x^{T}\mathbf{x} +\widehat{k}_r r. \]

The matching-condition hypothesis states that constant ideal parameters \( \mathbf{k}_x^{*} \) and \( k_r^{*} \) exist such that

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

These equalities do not say that the ideal parameters are known. They say that the selected controller structure is capable of reproducing the desired model dynamics if those parameters were available.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> XM["Desired state x_m"]
  R --> REG["Regressor omega = [x ; r]"]
  X["Plant state x"] --> REG
  REG --> CTRL["Adjustable controller theta_hat"]
  CTRL --> U["Control input u"]
  U --> PLANT["Unknown plant"]
  PLANT --> X
  X --> ERR["Tracking error e = x - x_m"]
  XM --> ERR
  ERR --> ADAPT["Parameter update"]
  REG --> ADAPT
  ADAPT --> CTRL
        

Substitute the adjustable controller into the plant:

\[ \dot{\mathbf{x} } =\left(\mathbf{A} +\mathbf{b}\widehat{\mathbf{k} }_x^{T}\right)\mathbf{x} +\mathbf{b}\widehat{k}_r r. \]

Add and subtract the ideal gains:

\[ \begin{aligned} \dot{\mathbf{x} } &= \left(\mathbf{A} +\mathbf{b}\mathbf{k}_x^{*T}\right)\mathbf{x} +\mathbf{b}k_r^{*}r \\ &\quad +\mathbf{b} \left[ \left(\widehat{\mathbf{k} }_x-\mathbf{k}_x^{*}\right)^{T}\mathbf{x} +\left(\widehat{k}_r-k_r^{*}\right)r \right]. \end{aligned} \]

Apply the matching conditions:

\[ \dot{\mathbf{x} } =\mathbf{A}_m\mathbf{x}+\mathbf{b}_m r +\mathbf{b} \left( \widetilde{\mathbf{k} }_x^{T}\mathbf{x} +\widetilde{k}_r r \right). \]

Since \( \mathbf{e}=\mathbf{x}-\mathbf{x}_m \), subtract the reference-model equation:

\[ \begin{aligned} \dot{\mathbf{e} } &=\dot{\mathbf{x} }-\dot{\mathbf{x} }_m \\ &=\mathbf{A}_m(\mathbf{x}-\mathbf{x}_m) +\mathbf{b} \left( \widetilde{\mathbf{k} }_x^{T}\mathbf{x} +\widetilde{k}_r r \right). \end{aligned} \]

Define the stacked regressor and parameter vectors

\[ \boldsymbol{\omega} = \begin{bmatrix} \mathbf{x}\\ r \end{bmatrix}, \qquad \widehat{\boldsymbol{\theta} } = \begin{bmatrix} \widehat{\mathbf{k} }_x\\ \widehat{k}_r \end{bmatrix}, \qquad \boldsymbol{\theta}^{*} = \begin{bmatrix} \mathbf{k}_x^{*}\\ k_r^{*} \end{bmatrix}. \]

The compact MRAC error model is therefore

\[ \boxed{ \dot{\mathbf{e} } =\mathbf{A}_m\mathbf{e} +\mathbf{b}\, \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} } \]

This equation has two parts: the desired stable dynamics \( \mathbf{A}_m\mathbf{e} \) and a forcing term generated by the instantaneous parameter mismatch projected onto the regressor. Adaptation is designed to control the coupled dynamics of both \( \mathbf{e} \) and \( \widetilde{\boldsymbol{\theta} } \).

Variation of constants gives the exact representation

\[ \mathbf{e}(t) = \mathrm{e}^{\mathbf{A}_m t}\mathbf{e}(0) + \int_{0}^{t} \mathrm{e}^{\mathbf{A}_m(t-\tau)} \mathbf{b}\, \widetilde{\boldsymbol{\theta} }^{T}(\tau) \boldsymbol{\omega}(\tau) \,d\tau. \]

Because \( \mathbf{A}_m \) is Hurwitz, constants \( M\ge 1 \) and \( \alpha>0 \) exist such that

\[ \left\| \mathrm{e}^{\mathbf{A}_m t} \right\| \le M\mathrm{e}^{-\alpha t}. \]

Consequently,

\[ \begin{aligned} \|\mathbf{e}(t)\| &\le M\mathrm{e}^{-\alpha t}\|\mathbf{e}(0)\|\\ &\quad+ M\|\mathbf{b}\| \int_{0}^{t} \mathrm{e}^{-\alpha(t-\tau)} \left| \widetilde{\boldsymbol{\theta} }^{T}(\tau) \boldsymbol{\omega}(\tau) \right| d\tau. \end{aligned} \]

A stable reference model removes the initial tracking error, but it does not automatically remove the forcing caused by parameter mismatch. The adaptive law must organize the parameter motion so that the coupled system remains stable and the forcing becomes compatible with tracking convergence.

Let \( \mathbf{Q}=\mathbf{Q}^{T}\succ\mathbf{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 positive-definite solution \( \mathbf{P}=\mathbf{P}^{T}\succ\mathbf{0} \). Choose \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}\succ\mathbf{0} \) and define

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

The first term measures tracking-error energy in a metric compatible with the reference dynamics. The second term measures parameter-error energy, weighted by the inverse adaptation-gain matrix.

Differentiate along the coupled error dynamics:

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

A Lyapunov-compatible conceptual update law is

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } = -\boldsymbol{\Gamma} \boldsymbol{\omega} \mathbf{e}^{T}\mathbf{P}\mathbf{b} } \]

because \( \dot{\widetilde{\boldsymbol{\theta} } } =\dot{\widehat{\boldsymbol{\theta} } } \) for constant ideal parameters. Substitution gives

\[ \begin{aligned} \dot V &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} + 2\mathbf{e}^{T}\mathbf{P}\mathbf{b}\, \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega}\\ &\quad- 2\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} \mathbf{e}^{T}\mathbf{P}\mathbf{b}\\ &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \le 0. \end{aligned} \]

The two mixed state–parameter terms are identical scalars with opposite signs. This cancellation is the central design mechanism of basic Lyapunov MRAC.

Since \( V(t)\le V(0) \), both \( \mathbf{e}(t) \) and \( \widetilde{\boldsymbol{\theta} }(t) \) are bounded. Integration of \( \dot V \) gives

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

so \( \mathbf{e}\in L_2 \). If the command, plant state, and therefore the regressor are bounded, then the error model and adaptive law imply bounded \( \dot{\mathbf{e} } \). Barbalat’s lemma, reviewed in Chapter 3, then yields

\[ \lim_{t → \infty}\mathbf{e}(t)=\mathbf{0}. \]

However, the derivative contains no negative-definite term of the form \( -\|\widetilde{\boldsymbol{\theta} }\|^2 \). Therefore this proof establishes bounded parameter error, not necessarily \( \widetilde{\boldsymbol{\theta} }(t) → \mathbf{0} \).

The result is a tracking theorem. It is not automatically an identification theorem.

The parameter error enters the state dynamics only through the scalar projection \( \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} \). If the regressor does not expose a direction in parameter space, error in that direction cannot influence the tracking dynamics and cannot be corrected by the update law.

For example, let

\[ \boldsymbol{\omega}(t) = \begin{bmatrix} \omega_1(t)\\ 0 \end{bmatrix}, \qquad \widetilde{\boldsymbol{\theta} } = \begin{bmatrix} \widetilde{\theta}_1\\ \widetilde{\theta}_2 \end{bmatrix}. \]

Then

\[ \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} = \widetilde{\theta}_1\omega_1, \qquad \dot{\widehat{\theta} }_2=0. \]

The second parameter may remain wrong forever while the tracking error still converges to zero. More generally, the accumulated regressor matrix

\[ \mathbf{G}(t_0,t_1) = \int_{t_0}^{t_1} \boldsymbol{\omega}(\tau) \boldsymbol{\omega}^{T}(\tau) \, d\tau \]

reveals which parameter directions have been activated over an interval. If this matrix is rank deficient, at least one nonzero parameter direction is invisible to the collected data. Formal conditions for parameter convergence are postponed to Chapter 10.

flowchart TD
  PE["Parameter error theta_tilde"] --> PROJ["Projection theta_tilde^T omega"]
  REG["Regressor omega"] --> PROJ
  PROJ --> FORCE["Forcing of stable error dynamics"]
  FORCE --> E["Tracking error e"]
  E --> UPDATE["Adaptive update"]
  REG --> UPDATE
  UPDATE --> PE
  E --> ZERO{"Does e approach zero?"}
  ZERO -->|yes| TRACK["Tracking objective achieved"]
  TRACK --> RICH{"Did omega activate \nevery parameter direction?"}
  RICH -->|not necessarily| BOUNDED["Parameter error may \nremain nonzero"]
  RICH -->|sufficiently rich| IDENT["Parameter convergence \nmay be established"]
        

A scalar version makes the coupling transparent:

\[ \dot e=-a_m e+\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}, \qquad a_m>0. \]

Choose

\[ V=\frac{1}{2}e^2 +\frac{1}{2} \widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \widetilde{\boldsymbol{\theta} }, \qquad \dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\phi}e. \]

Then

\[ \begin{aligned} \dot V &= e\left(-a_m e +\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\right) + \widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \left(-\boldsymbol{\Gamma}\boldsymbol{\phi}e\right)\\ &=-a_m e^2. \end{aligned} \]

Adaptation transfers energy between the tracking-error and parameter-error components of \( V \), while the stable reference-model term dissipates tracking-error energy. The parameter component need not vanish because \( \dot V \) does not directly dissipate it.

All five implementations below simulate the same abstract error system:

\[ \dot{\mathbf{e} } = \begin{bmatrix} 0 & 1\\ -2 & -3 \end{bmatrix}\mathbf{e} + \begin{bmatrix} 0\\ 1 \end{bmatrix} \widetilde{\boldsymbol{\theta} }^{T} \begin{bmatrix} \sin(0.7t)\\ \cos(1.1t) \end{bmatrix}. \]

With \( \mathbf{Q}=\mathbf{I} \),

\[ \mathbf{P} = \begin{bmatrix} 1.25 & 0.25\\ 0.25 & 0.25 \end{bmatrix} \]

solves \( \mathbf{A}_m^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_m=-\mathbf{I} \). The experiment is not a new plant-specific MRAC design; it isolates the coupled tracking-error and parameter-error equations derived in this lesson.

Chapter4_Lesson2.py

"""
Chapter4_Lesson2.py
Numerical experiment for the MRAC conceptual error model

    e_dot = A_m e + b * (tilde_theta^T phi)
    hat_theta_dot = -Gamma * phi * (b^T P e)

where tilde_theta = hat_theta - theta_star and
A_m^T P + P A_m = -Q.

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    am: np.ndarray
    b: np.ndarray
    p: np.ndarray
    gamma: np.ndarray
    theta_star: np.ndarray


def regressor(t: float) -> np.ndarray:
    """A bounded, signal-rich regressor used only for the conceptual experiment."""
    return np.array([np.sin(0.7 * t), np.cos(1.1 * t)], dtype=float)


def derivative(t: float, z: np.ndarray, prm: Parameters) -> np.ndarray:
    """Return [e_dot; hat_theta_dot]."""
    e = z[:2]
    theta_hat = z[2:]
    phi = regressor(t)
    theta_tilde = theta_hat - prm.theta_star

    forcing = float(theta_tilde @ phi)
    e_dot = prm.am @ e + prm.b * forcing

    filtered_error = float(prm.b @ prm.p @ e)
    theta_hat_dot = -prm.gamma @ phi * filtered_error

    return np.concatenate((e_dot, theta_hat_dot))


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


def lyapunov_value(z: np.ndarray, prm: Parameters) -> float:
    e = z[:2]
    theta_tilde = z[2:] - prm.theta_star
    gamma_inv = np.linalg.inv(prm.gamma)
    return float(e @ prm.p @ e + theta_tilde @ gamma_inv @ theta_tilde)


def main() -> None:
    am = np.array([[0.0, 1.0], [-2.0, -3.0]])
    b = np.array([0.0, 1.0])
    q = np.eye(2)

    # Exact solution of A_m^T P + P A_m = -I for this A_m.
    p = np.array([[1.25, 0.25], [0.25, 0.25]])
    gamma = np.diag([4.0, 3.0])
    theta_star = np.array([1.2, -0.7])
    prm = Parameters(am=am, b=b, p=p, gamma=gamma, theta_star=theta_star)

    residual = am.T @ p + p @ am + q
    if np.linalg.norm(residual, ord=np.inf) > 1.0e-12:
        raise RuntimeError("The supplied P does not solve the Lyapunov equation.")

    dt = 0.002
    final_time = 30.0
    steps = int(round(final_time / dt))

    time = np.linspace(0.0, final_time, steps + 1)
    history = np.zeros((steps + 1, 4))
    history[0] = np.array([1.0, -0.5, 0.0, 0.0])

    for k in range(steps):
        history[k + 1] = rk4_step(time[k], history[k], dt, prm)

    e_history = history[:, :2]
    theta_hat_history = history[:, 2:]
    theta_tilde_history = theta_hat_history - theta_star
    v_history = np.array([lyapunov_value(z, prm) for z in history])

    csv_path = Path("Chapter4_Lesson2_results_python.csv")
    with csv_path.open("w", newline="", encoding="utf-8") as file:
        writer = csv.writer(file)
        writer.writerow(
            ["time", "e1", "e2", "theta_hat_1", "theta_hat_2",
             "theta_tilde_1", "theta_tilde_2", "V"]
        )
        for k, t in enumerate(time):
            writer.writerow([t, *history[k], *theta_tilde_history[k], v_history[k]])

    print(f"Final tracking-error norm: {np.linalg.norm(e_history[-1]):.6e}")
    print(f"Final parameter error: {theta_tilde_history[-1]}")
    print(f"Initial V: {v_history[0]:.6f}")
    print(f"Final V:   {v_history[-1]:.6f}")
    print(f"Saved: {csv_path.resolve()}")

    plt.figure()
    plt.plot(time, e_history[:, 0], label="e1")
    plt.plot(time, e_history[:, 1], label="e2")
    plt.xlabel("Time (s)")
    plt.ylabel("Tracking error")
    plt.title("Tracking-error dynamics")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter4_Lesson2_tracking_error_python.png", dpi=180)

    plt.figure()
    plt.plot(time, theta_hat_history[:, 0], label="theta_hat_1")
    plt.plot(time, theta_hat_history[:, 1], label="theta_hat_2")
    plt.axhline(theta_star[0], linestyle="--", label="theta_star_1")
    plt.axhline(theta_star[1], linestyle="--", label="theta_star_2")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.title("Parameter estimates")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter4_Lesson2_parameters_python.png", dpi=180)

    plt.figure()
    plt.plot(time, v_history)
    plt.xlabel("Time (s)")
    plt.ylabel("V")
    plt.title("Composite Lyapunov function")
    plt.grid(True)
    plt.tight_layout()
    plt.savefig("Chapter4_Lesson2_lyapunov_python.png", dpi=180)

    plt.show()


if __name__ == "__main__":
    main()

Chapter4_Lesson2.cpp

/*
Chapter4_Lesson2.cpp
From-scratch C++17 RK4 simulation of the conceptual MRAC error model.

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

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

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

struct Parameters {
    std::array<std::array<double, 2>, 2> am{};
    std::array<double, 2> b{};
    std::array<std::array<double, 2>, 2> p{};
    std::array<double, 2> gamma{};
    std::array<double, 2> theta_star{};
};

std::array<double, 2> regressor(double t) {
    return {std::sin(0.7 * t), std::cos(1.1 * t)};
}

State derivative(double t, const State& z, const Parameters& prm) {
    const std::array<double, 2> e{z[0], z[1]};
    const std::array<double, 2> theta_hat{z[2], z[3]};
    const auto phi = regressor(t);

    const std::array<double, 2> theta_tilde{
        theta_hat[0] - prm.theta_star[0],
        theta_hat[1] - prm.theta_star[1]
    };

    const double forcing =
        theta_tilde[0] * phi[0] + theta_tilde[1] * phi[1];

    std::array<double, 2> e_dot{};
    for (std::size_t i = 0; i < 2; ++i) {
        e_dot[i] =
            prm.am[i][0] * e[0] +
            prm.am[i][1] * e[1] +
            prm.b[i] * forcing;
    }

    const std::array<double, 2> pe{
        prm.p[0][0] * e[0] + prm.p[0][1] * e[1],
        prm.p[1][0] * e[0] + prm.p[1][1] * e[1]
    };
    const double filtered_error = prm.b[0] * pe[0] + prm.b[1] * pe[1];

    const std::array<double, 2> theta_hat_dot{
        -prm.gamma[0] * phi[0] * filtered_error,
        -prm.gamma[1] * phi[1] * filtered_error
    };

    return {e_dot[0], e_dot[1], theta_hat_dot[0], theta_hat_dot[1]};
}

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

State rk4_step(double t, const State& z, double dt, const Parameters& prm) {
    const State k1 = derivative(t, z, prm);
    const State k2 = derivative(t + 0.5 * dt, add_scaled(z, k1, 0.5 * dt), prm);
    const State k3 = derivative(t + 0.5 * dt, add_scaled(z, k2, 0.5 * dt), prm);
    const State k4 = derivative(t + dt, add_scaled(z, k3, dt), prm);

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

double lyapunov_value(const State& z, const Parameters& prm) {
    const std::array<double, 2> e{z[0], z[1]};
    const std::array<double, 2> theta_tilde{
        z[2] - prm.theta_star[0],
        z[3] - prm.theta_star[1]
    };

    const double state_energy =
        e[0] * (prm.p[0][0] * e[0] + prm.p[0][1] * e[1]) +
        e[1] * (prm.p[1][0] * e[0] + prm.p[1][1] * e[1]);

    const double parameter_energy =
        theta_tilde[0] * theta_tilde[0] / prm.gamma[0] +
        theta_tilde[1] * theta_tilde[1] / prm.gamma[1];

    return state_energy + parameter_energy;
}

int main() {
    try {
        const Parameters prm{
            { { {0.0, 1.0}, {-2.0, -3.0} } },
            { {0.0, 1.0} },
            { { {1.25, 0.25}, {0.25, 0.25} } },
            { {4.0, 3.0} },
            { {1.2, -0.7} }
        };

        const double dt = 0.002;
        const double final_time = 30.0;
        const int steps = static_cast<int>(std::llround(final_time / dt));

        State z{1.0, -0.5, 0.0, 0.0};
        const double initial_v = lyapunov_value(z, prm);

        std::ofstream csv("Chapter4_Lesson2_results_cpp.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open output CSV file.");
        }

        csv << "time,e1,e2,theta_hat_1,theta_hat_2,"
               "theta_tilde_1,theta_tilde_2,V\n";
        csv << std::setprecision(15);

        for (int k = 0; k <= steps; ++k) {
            const double t = k * dt;
            const double tilde_1 = z[2] - prm.theta_star[0];
            const double tilde_2 = z[3] - prm.theta_star[1];
            csv << t << ',' << z[0] << ',' << z[1] << ','
                << z[2] << ',' << z[3] << ','
                << tilde_1 << ',' << tilde_2 << ','
                << lyapunov_value(z, prm) << '\n';

            if (k < steps) {
                z = rk4_step(t, z, dt, prm);
            }
        }

        const double final_error_norm = std::hypot(z[0], z[1]);
        std::cout << std::scientific;
        std::cout << "Final tracking-error norm: " << final_error_norm << '\n';
        std::cout << "Final parameter error: ["
                  << z[2] - prm.theta_star[0] << ", "
                  << z[3] - prm.theta_star[1] << "]\n";
        std::cout << "Initial V: " << initial_v << '\n';
        std::cout << "Final V:   " << lyapunov_value(z, prm) << '\n';
        std::cout << "Saved Chapter4_Lesson2_results_cpp.csv\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter4_Lesson2.java

/*
Chapter4_Lesson2.java
From-scratch Java simulation of the conceptual MRAC error model.

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

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;

public final class Chapter4_Lesson2 {
    private static final double[][] AM = {
        {0.0, 1.0},
        {-2.0, -3.0}
    };
    private static final double[] B = {0.0, 1.0};
    private static final double[][] P = {
        {1.25, 0.25},
        {0.25, 0.25}
    };
    private static final double[] GAMMA = {4.0, 3.0};
    private static final double[] THETA_STAR = {1.2, -0.7};

    private Chapter4_Lesson2() {
    }

    private static double[] regressor(double t) {
        return new double[] {Math.sin(0.7 * t), Math.cos(1.1 * t)};
    }

    private static double[] derivative(double t, double[] z) {
        double[] phi = regressor(t);
        double tilde0 = z[2] - THETA_STAR[0];
        double tilde1 = z[3] - THETA_STAR[1];
        double forcing = tilde0 * phi[0] + tilde1 * phi[1];

        double eDot0 = AM[0][0] * z[0] + AM[0][1] * z[1] + B[0] * forcing;
        double eDot1 = AM[1][0] * z[0] + AM[1][1] * z[1] + B[1] * forcing;

        double pe0 = P[0][0] * z[0] + P[0][1] * z[1];
        double pe1 = P[1][0] * z[0] + P[1][1] * z[1];
        double filteredError = B[0] * pe0 + B[1] * pe1;

        double thetaDot0 = -GAMMA[0] * phi[0] * filteredError;
        double thetaDot1 = -GAMMA[1] * phi[1] * filteredError;

        return new double[] {eDot0, eDot1, thetaDot0, thetaDot1};
    }

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

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

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

    private static double lyapunovValue(double[] z) {
        double e0 = z[0];
        double e1 = z[1];
        double tilde0 = z[2] - THETA_STAR[0];
        double tilde1 = z[3] - THETA_STAR[1];

        double stateEnergy =
            e0 * (P[0][0] * e0 + P[0][1] * e1)
            + e1 * (P[1][0] * e0 + P[1][1] * e1);

        double parameterEnergy =
            tilde0 * tilde0 / GAMMA[0]
            + tilde1 * tilde1 / GAMMA[1];

        return stateEnergy + parameterEnergy;
    }

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

        double dt = 0.002;
        double finalTime = 30.0;
        int steps = (int) Math.round(finalTime / dt);
        double[] z = {1.0, -0.5, 0.0, 0.0};
        double initialV = lyapunovValue(z);

        Path csvPath = Path.of("Chapter4_Lesson2_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                csvPath, StandardCharsets.UTF_8)) {
            writer.write(
                "time,e1,e2,theta_hat_1,theta_hat_2,"
                + "theta_tilde_1,theta_tilde_2,V\n"
            );

            for (int k = 0; k <= steps; k++) {
                double t = k * dt;
                writer.write(String.format(
                    Locale.ROOT,
                    "%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g%n",
                    t, z[0], z[1], z[2], z[3],
                    z[2] - THETA_STAR[0],
                    z[3] - THETA_STAR[1],
                    lyapunovValue(z)
                ));

                if (k < steps) {
                    z = rk4Step(t, z, dt);
                }
            }
        } catch (IOException ex) {
            System.err.println("Cannot write CSV: " + ex.getMessage());
            System.exit(1);
        }

        System.out.printf(
            Locale.ROOT,
            "Final tracking-error norm: %.6e%n",
            Math.hypot(z[0], z[1])
        );
        System.out.printf(
            Locale.ROOT,
            "Final parameter error: [%.6e, %.6e]%n",
            z[2] - THETA_STAR[0],
            z[3] - THETA_STAR[1]
        );
        System.out.printf(Locale.ROOT, "Initial V: %.6f%n", initialV);
        System.out.printf(Locale.ROOT, "Final V:   %.6f%n", lyapunovValue(z));
        System.out.println("Saved " + csvPath.toAbsolutePath());
    }
}

Chapter4_Lesson2.m

% Chapter4_Lesson2.m
% Conceptual MRAC error-model simulation:
%   e_dot = A_m e + b*(tilde_theta' * phi)
%   theta_hat_dot = -Gamma*phi*(b' * P * e)
%
% Requires MATLAB. The Control System Toolbox is used for lyap when present.

clear; clc; close all;

Am = [0 1; -2 -3];
b = [0; 1];
Q = eye(2);
Gamma = diag([4 3]);
thetaStar = [1.2; -0.7];

if exist('lyap', 'file') == 2
    P = lyap(Am', Q);
else
    % Exact solution for this numerical A_m and Q.
    P = [1.25 0.25; 0.25 0.25];
end

residual = Am' * P + P * Am + Q;
assert(norm(residual, inf) < 1e-10, ...
    'P does not satisfy the continuous Lyapunov equation.');

z0 = [1.0; -0.5; 0.0; 0.0];
tspan = [0 30];

rhs = @(t, z) errorModelRhs(t, z, Am, b, P, Gamma, thetaStar);
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, z] = ode45(rhs, tspan, z0, options);

e = z(:, 1:2);
thetaHat = z(:, 3:4);
thetaTilde = thetaHat - thetaStar.';
V = zeros(size(t));

GammaInv = inv(Gamma);
for k = 1:numel(t)
    ek = e(k, :).';
    tk = thetaTilde(k, :).';
    V(k) = ek' * P * ek + tk' * GammaInv * tk;
end

results = table( ...
    t, e(:, 1), e(:, 2), ...
    thetaHat(:, 1), thetaHat(:, 2), ...
    thetaTilde(:, 1), thetaTilde(:, 2), V, ...
    'VariableNames', { ...
        'time', 'e1', 'e2', ...
        'theta_hat_1', 'theta_hat_2', ...
        'theta_tilde_1', 'theta_tilde_2', 'V' ...
    } ...
);
writetable(results, 'Chapter4_Lesson2_results_matlab.csv');

fprintf('Final tracking-error norm: %.6e\n', norm(e(end, :)));
fprintf('Final parameter error: [%.6e, %.6e]\n', ...
    thetaTilde(end, 1), thetaTilde(end, 2));
fprintf('Initial V: %.6f\n', V(1));
fprintf('Final V:   %.6f\n', V(end));

figure;
plot(t, e(:, 1), t, e(:, 2), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Tracking error');
title('Tracking-error dynamics');
legend('e_1', 'e_2', 'Location', 'best');

figure;
plot(t, thetaHat(:, 1), t, thetaHat(:, 2), 'LineWidth', 1.2);
hold on;
yline(thetaStar(1), '--');
yline(thetaStar(2), '--');
grid on;
xlabel('Time (s)');
ylabel('Parameter estimate');
title('Parameter estimates');
legend('\hat{\theta}_1', '\hat{\theta}_2', ...
    '\theta_1^*', '\theta_2^*', 'Location', 'best');

figure;
plot(t, V, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('V');
title('Composite Lyapunov function');

% Simulink implementation note:
% Use one State-Space block for A_m and b, a MATLAB Function block for
% phi(t), Product blocks for tilde_theta' * phi and b' * P * e, and an
% Integrator block for theta_hat_dot. The ODE implemented below is the
% continuous-time equation that those blocks realize.

function dz = errorModelRhs(t, z, Am, b, P, Gamma, thetaStar)
    e = z(1:2);
    thetaHat = z(3:4);
    phi = [sin(0.7 * t); cos(1.1 * t)];
    thetaTilde = thetaHat - thetaStar;

    eDot = Am * e + b * (thetaTilde' * phi);
    filteredError = b' * P * e;
    thetaHatDot = -Gamma * phi * filteredError;

    dz = [eDot; thetaHatDot];
end

The MATLAB script states the corresponding Simulink realization. A State-Space block implements \( \mathbf{A}_m \) and \( \mathbf{b} \); a Clock and MATLAB Function block generate the regressor; Product blocks form \( \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \) and \( \mathbf{b}^{T}\mathbf{P}\mathbf{e} \); and an Integrator block produces \( \widehat{\boldsymbol{\theta} } \).

Chapter4_Lesson2.nb

(* Chapter4_Lesson2.nb
   Wolfram Language implementation of the conceptual MRAC error model. *)

ClearAll["Global`*"];

am = { {0., 1.}, {-2., -3.} };
b = {0., 1.};
q = IdentityMatrix[2];
gamma = DiagonalMatrix[{4., 3.}];
thetaStar = {1.2, -0.7};

p = LyapunovSolve[Transpose[am], q];
residual = Transpose[am].p + p.am + q;

If[Norm[residual, Infinity] > 10^-10,
  Print["Lyapunov-equation residual is too large."];
  Abort[];
];

phi[t_] := {Sin[0.7 t], Cos[1.1 t]};

equations = {
  e1'[t] == (am.{e1[t], e2[t]} +
       b.(({th1[t], th2[t]} - thetaStar).phi[t]))[[1]],
  e2'[t] == (am.{e1[t], e2[t]} +
       b.(({th1[t], th2[t]} - thetaStar).phi[t]))[[2]],
  th1'[t] == (-gamma.phi[t]*
       (b.p.{e1[t], e2[t]}))[[1]],
  th2'[t] == (-gamma.phi[t]*
       (b.p.{e1[t], e2[t]}))[[2]],
  e1[0] == 1.,
  e2[0] == -0.5,
  th1[0] == 0.,
  th2[0] == 0.
};

solution = NDSolveValue[
  equations,
  {e1, e2, th1, th2},
  {t, 0, 30},
  Method -> {"TimeIntegration" -> {"ExplicitRungeKutta",
      "DifferenceOrder" -> 4} }
];

e[t_] := {solution[[1]][t], solution[[2]][t]};
thetaHat[t_] := {solution[[3]][t], solution[[4]][t]};
thetaTilde[t_] := thetaHat[t] - thetaStar;
v[t_] := e[t].p.e[t] +
  thetaTilde[t].Inverse[gamma].thetaTilde[t];

sampleTimes = Subdivide[0., 30., 3000];
data = Table[
  Join[{tt}, e[tt], thetaHat[tt], thetaTilde[tt], {v[tt]}],
  {tt, sampleTimes}
];

Export[
  "Chapter4_Lesson2_results_mathematica.csv",
  Prepend[data, {"time", "e1", "e2", "theta_hat_1",
    "theta_hat_2", "theta_tilde_1", "theta_tilde_2", "V"}]
];

Print["Final tracking-error norm: ", Norm[e[30.]]];
Print["Final parameter error: ", thetaTilde[30.]];
Print["Initial V: ", v[0.]];
Print["Final V: ", v[30.]];

trackingPlot = Plot[
  Evaluate[e[t]],
  {t, 0, 30},
  PlotLegends -> {"e1", "e2"},
  AxesLabel -> {"Time (s)", "Tracking error"},
  PlotLabel -> "Tracking-error dynamics",
  GridLines -> Automatic
];

parameterPlot = Plot[
  Evaluate[Join[thetaHat[t], thetaStar]],
  {t, 0, 30},
  PlotLegends -> {"thetaHat1", "thetaHat2", "thetaStar1", "thetaStar2"},
  AxesLabel -> {"Time (s)", "Parameter"},
  PlotLabel -> "Parameter estimates",
  GridLines -> Automatic
];

lyapunovPlot = Plot[
  v[t],
  {t, 0, 30},
  AxesLabel -> {"Time (s)", "V"},
  PlotLabel -> "Composite Lyapunov function",
  GridLines -> Automatic
];

Column[{trackingPlot, parameterPlot, lyapunovPlot}]

Each program records the tracking error, parameter estimate, parameter error, and composite Lyapunov function. The expected observations are:

• the tracking-error norm decreases toward zero;
• the Lyapunov function is nonincreasing up to numerical integration error;
• the parameter estimates remain bounded;
• parameter estimates may approach the ideal values in this signal-rich example, but that behavior is stronger than the basic tracking proof guarantees;
• changing the regressor to \( [\sin(0.7t),\,0]^{T} \) freezes the second estimate and demonstrates non-identifiability directly.

Adaptation gains in \( \boldsymbol{\Gamma} \) change the speed of parameter motion. Large gains are not automatically better; noise, sampling, delays, and unmodeled dynamics make aggressive adaptation a practical concern addressed in later chapters.

Problem 1 (Error-model derivation): For \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{b}u \), \( u=\widehat{\mathbf{k} }_x^{T}\mathbf{x}+\widehat{k}_r r \), and \( \dot{\mathbf{x} }_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{b}_m r \), assume ideal gains satisfy the matching conditions. Derive the compact tracking-error equation.

Solution:

\[ \begin{aligned} \dot{\mathbf{e} } &= \dot{\mathbf{x} }-\dot{\mathbf{x} }_m\\ &= \left(\mathbf{A} +\mathbf{b}\widehat{\mathbf{k} }_x^{T}\right)\mathbf{x} +\mathbf{b}\widehat{k}_r r -\mathbf{A}_m\mathbf{x}_m-\mathbf{b}_m r. \end{aligned} \]

Insert \( \widehat{\mathbf{k} }_x=\mathbf{k}_x^{*} +\widetilde{\mathbf{k} }_x \) and \( \widehat{k}_r=k_r^{*}+\widetilde{k}_r \), then use \( \mathbf{A}+\mathbf{b}\mathbf{k}_x^{*T}=\mathbf{A}_m \) and \( \mathbf{b}k_r^{*}=\mathbf{b}_m \). The result is

\[ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} +\mathbf{b} \left( \widetilde{\mathbf{k} }_x^{T}\mathbf{x} +\widetilde{k}_r r \right) = \mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\omega}. \]

Problem 2 (Lyapunov cancellation): Given \( \dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} \), prove that \( \dot{\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\Gamma}\boldsymbol{\omega} \mathbf{e}^{T}\mathbf{P}\mathbf{b} \) produces \( \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \).

Solution:

\[ \begin{aligned} \dot V &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} +2\mathbf{e}^{T}\mathbf{P}\mathbf{b} \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega}\\ &\quad+ 2\widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \left( -\boldsymbol{\Gamma}\boldsymbol{\omega} \mathbf{e}^{T}\mathbf{P}\mathbf{b} \right). \end{aligned} \]

The last term equals \( -2\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} \mathbf{e}^{T}\mathbf{P}\mathbf{b} \). Both mixed terms are scalars and therefore cancel exactly, leaving \( \dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}\le0 \).

Problem 3 (Tracking without identification): Suppose \( \boldsymbol{\omega}(t)=[\omega_1(t),\,0]^{T} \) and \( \widehat{\theta}_2(0)\ne\theta_2^{*} \). Show that the second parameter cannot converge under the gradient update.

Solution:

\[ \dot{\widehat{\theta} }_2 = -\gamma_2\omega_2 \mathbf{e}^{T}\mathbf{P}\mathbf{b} =0. \]

Hence

\[ \widehat{\theta}_2(t)=\widehat{\theta}_2(0), \qquad \widetilde{\theta}_2(t) =\widehat{\theta}_2(0)-\theta_2^{*}. \]

The second parameter error remains constant and nonzero, yet it does not enter the forcing term because \( \omega_2=0 \). Tracking can therefore converge independently of that parameter.

Problem 4 (Solve a Lyapunov equation): For \( \mathbf{A}_m=\begin{bmatrix}0&1\\-2&-3\end{bmatrix} \) and \( \mathbf{Q}=\mathbf{I} \), let \( \mathbf{P}=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix} \). Find \( \mathbf{P} \).

Solution: Expanding \( \mathbf{A}_m^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_m=-\mathbf{I} \) gives

\[ -4p_{12}=-1, \qquad p_{11}-3p_{12}-2p_{22}=0, \qquad 2p_{12}-6p_{22}=-1. \]

Therefore

\[ p_{12}=\frac{1}{4}, \qquad p_{22}=\frac{1}{4}, \qquad p_{11}=\frac{5}{4}, \]

and

\[ \mathbf{P} = \begin{bmatrix} 5/4 & 1/4\\ 1/4 & 1/4 \end{bmatrix}. \]

Its leading principal minors are \( 5/4 \) and \( 1/4 \), so \( \mathbf{P}\succ\mathbf{0} \).

Problem 5 (Output error versus state error): Let \( e_y=\mathbf{c}^{T}\mathbf{e} \). Does \( e_y(t) → 0 \) necessarily imply \( \mathbf{e}(t) → \mathbf{0} \)?

Solution:

Not in general. The equality \( e_y=\mathbf{c}^{T}\mathbf{e}=0 \) constrains the state error to the null space of \( \mathbf{c}^{T} \). For \( n>1 \), this null space usually contains nonzero vectors. Full state convergence requires additional system structure and analysis. This lesson proves state-error convergence because its Lyapunov function contains the complete vector \( \mathbf{e} \).

The MRAC error model separates desired stable dynamics from uncertainty: \( \dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +\mathbf{b}\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\omega} \). For constant ideal parameters, \( \dot{\widetilde{\boldsymbol{\theta} } } =\dot{\widehat{\boldsymbol{\theta} } } \). A composite Lyapunov function combines state-tracking and parameter-error energies, and the adaptive law is selected to cancel their mixed term. The resulting proof establishes bounded estimates and asymptotic tracking under bounded-signal assumptions. It does not by itself establish parameter convergence, because only the regressor-projected parameter error affects the tracking dynamics.

1. Parks, P.C. (1966). Lyapunov redesign of model reference adaptive control systems. IEEE Transactions on Automatic Control, 11(3), 362–367. doi:10.1109/TAC.1966.1098361.
2. Landau, I.D. (1969). A hyperstability criterion for model reference adaptive control systems. IEEE Transactions on Automatic Control, 14(5), 552–555. doi:10.1109/TAC.1969.1099237.
3. Monopoli, R.V. (1974). Model reference adaptive control with an augmented error signal. IEEE Transactions on Automatic Control, 19(5), 474–484. doi:10.1109/TAC.1974.1100670.
4. Narendra, K.S., & Valavani, L.S. (1979). Direct and indirect model reference adaptive control. Automatica, 15(6), 653–664. doi:10.1016/0005-1098(79)90033-5.
5. Narendra, K.S., & Valavani, L.S. (1980). A comparison of Lyapunov and hyperstability approaches to adaptive control of continuous systems. IEEE Transactions on Automatic Control, 25(2), 243–247.
6. Winsor, C.A., & Roy, R.J. (1968). Design of model reference adaptive control systems by Lyapunov’s second method. IEEE Transactions on Automatic Control, 13(2), 204.