After completing this lesson, students should be able to:

• express a second-order uncertain plant and its reference model in controllable canonical state-space form;
• derive the ideal matching equations and distinguish the unknown ideal controller parameters from their online estimates;
• obtain the compact error equation \( \dot{\mathbf e}=\mathbf A_m\mathbf e+\mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega \);
• solve the continuous-time Lyapunov equation for a positive-definite matrix \( \mathbf P \);
• construct a composite state-parameter Lyapunov candidate and identify the cross term that the next lesson's update law must cancel; and
• verify the derivation numerically in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The parameter-update law is not finalized here. We derive only the algebraic condition it must satisfy. This preserves the course sequence: Lesson 3 develops the analysis model; Lesson 4 converts that model into implementable adaptation laws.

Consider the second-order plant introduced in the preceding lessons:

\[ \ddot y+a_1\dot y+a_0y=b\,u, \qquad b\neq 0. \]

Define \(x_1=y\) and \(x_2=\dot y\). The plant is

\[ \dot{\mathbf x} = \underbrace{\begin{bmatrix} 0 & 1\\ -a_0 & -a_1 \end{bmatrix} }_{\mathbf A_p}\mathbf x + \underbrace{\begin{bmatrix} 0\\ b \end{bmatrix} }_{\mathbf B}u, \qquad \mathbf x=\begin{bmatrix}x_1\\x_2\end{bmatrix}. \]

The desired second-order reference model is selected as

\[ \dot{\mathbf x}_m = \underbrace{\begin{bmatrix} 0 & 1\\ -a_{m0} & -a_{m1} \end{bmatrix} }_{\mathbf A_m}\mathbf x_m + \underbrace{\begin{bmatrix} 0\\ b_m \end{bmatrix} }_{\mathbf B_m}r, \]

where \(a_{m0}\) and \(a_{m1}\) are chosen so that \(\mathbf A_m\) is Hurwitz. For the introductory design in this chapter, the input direction and nonzero gain \(b\) are treated as known. More difficult high-frequency-gain assumptions are deferred to later chapters.

Use the linearly parameterized control law

\[ u = \theta_1x_1+\theta_2x_2+\theta_r r = \boldsymbol\theta^{T}\boldsymbol\omega, \qquad \boldsymbol\omega = \begin{bmatrix}x_1&x_2&r\end{bmatrix}^{T}. \]

The regressor \(\boldsymbol\omega\) is measurable. The adjustable controller vector is \(\boldsymbol\theta=[\theta_1,\theta_2,\theta_r]^T\).

Substitute the control law into the plant:

\[ \dot{\mathbf x} = \left( \mathbf A_p+ \mathbf B \begin{bmatrix}\theta_1&\theta_2\end{bmatrix} \right)\mathbf x + \mathbf B\theta_r r. \]

The ideal parameter vector \(\boldsymbol\theta^{*}\) is the constant vector that would make the plant's closed-loop dynamics identical to the reference model if the uncertain coefficients were known. It satisfies

\[ \mathbf A_p+ \mathbf B \begin{bmatrix}\theta_1^{*}&\theta_2^{*}\end{bmatrix} = \mathbf A_m, \qquad \mathbf B\theta_r^{*}=\mathbf B_m. \]

Comparing the second rows gives the scalar matching equations

\[ -a_0+b\theta_1^{*}=-a_{m0}, \qquad -a_1+b\theta_2^{*}=-a_{m1}, \qquad b\theta_r^{*}=b_m. \]

Therefore,

\[ \theta_1^{*}=\frac{a_0-a_{m0} }{b}, \qquad \theta_2^{*}=\frac{a_1-a_{m1} }{b}, \qquad \theta_r^{*}=\frac{b_m}{b}. \]

These equations establish the matching assumption. They do not imply that the controller knows \(\boldsymbol\theta^{*}\). Instead, adaptation attempts to adjust \(\boldsymbol\theta(t)\) so that the effect of the parameter error on tracking is controlled.

\[ \tilde{\boldsymbol\theta} = \boldsymbol\theta-\boldsymbol\theta^{*}, \qquad \dot{\tilde{\boldsymbol\theta} } = \dot{\boldsymbol\theta}, \]

because the ideal parameters are constant for the time-invariant plant considered here.

Define the state-tracking error as

\[ \mathbf e=\mathbf x-\mathbf x_m. \]

Write the actual parameter vector as \(\boldsymbol\theta=\boldsymbol\theta^{*}+ \tilde{\boldsymbol\theta}\). The plant dynamics become

\[ \begin{aligned} \dot{\mathbf x} &= \mathbf A_p\mathbf x+ \mathbf B \left(\boldsymbol\theta^{*}+ \tilde{\boldsymbol\theta}\right)^{T} \boldsymbol\omega\\ &= \mathbf A_p\mathbf x+ \mathbf B{\boldsymbol\theta^{*} }^{T}\boldsymbol\omega+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \end{aligned} \]

Apply the matching equations to the ideal-parameter part:

\[ \mathbf A_p\mathbf x+ \mathbf B{\boldsymbol\theta^{*} }^{T}\boldsymbol\omega = \mathbf A_m\mathbf x+\mathbf B_m r. \]

Hence,

\[ \dot{\mathbf x} = \mathbf A_m\mathbf x+\mathbf B_m r+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \]

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\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \end{aligned} \]

The central error equation is therefore

\[ \boxed{ \dot{\mathbf e} = \mathbf A_m\mathbf e+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega }. \]

In scalar coordinates, this equation is

\[ \dot e_1=e_2, \qquad \dot e_2= -a_{m0}e_1-a_{m1}e_2+ b\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \]

flowchart TD
  R["Command r"] --> RM["Reference model: \nx_m_dot = A_m x_m + B_m r"]
  R --> REG["Regressor omega = \n[x1, x2, r]"]
  XP["Plant state x"] --> REG
  TH["Parameter error theta_tilde"] --> MM["Mismatch scalar theta_tilde^T omega"]
  REG --> MM
  MM --> ERR["Error dynamics: \ne_dot = A_m e + B mismatch"]
  RM --> SUB["e = x - x_m"]
  XP --> SUB
  SUB --> ERR
        

First consider the ideal case \(\tilde{\boldsymbol\theta}=\mathbf 0\). Then

\[ \dot{\mathbf e}=\mathbf A_m\mathbf e. \]

Since \(\mathbf A_m\) is Hurwitz, for every symmetric matrix \(\mathbf Q\succ 0\), there is a unique symmetric matrix \(\mathbf P\succ 0\) satisfying

\[ \boxed{ \mathbf A_m^{T}\mathbf P+ \mathbf P\mathbf A_m = -\mathbf Q }. \]

Choose the nominal state-error energy

\[ V_e(\mathbf e)=\mathbf e^{T}\mathbf P\mathbf e. \]

Its derivative along the nominal error system is

\[ \begin{aligned} \dot V_e &= \dot{\mathbf e}^{T}\mathbf P\mathbf e+ \mathbf e^{T}\mathbf P\dot{\mathbf e}\\ &= \mathbf e^{T} \left( \mathbf A_m^{T}\mathbf P+ \mathbf P\mathbf A_m \right)\mathbf e\\ &= -\mathbf e^{T}\mathbf Q\mathbf e. \end{aligned} \]

Using eigenvalue bounds,

\[ \lambda_{\min}(\mathbf P)\|\mathbf e\|^2 \le V_e \le \lambda_{\max}(\mathbf P)\|\mathbf e\|^2, \]

\[ \dot V_e \le -\lambda_{\min}(\mathbf Q)\|\mathbf e\|^2 \le -\frac{\lambda_{\min}(\mathbf Q)} {\lambda_{\max}(\mathbf P)}V_e. \]

Thus, with perfect parameter matching, the tracking-error origin is exponentially stable. The adaptive problem is harder because \(\tilde{\boldsymbol\theta}\) is generally nonzero and time varying.

Evaluate the same state-only candidate along the actual error dynamics:

\[ \begin{aligned} \dot V_e &= 2\mathbf e^{T}\mathbf P \left( \mathbf A_m\mathbf e+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega \right)\\ &= -\mathbf e^{T}\mathbf Q\mathbf e+ 2\mathbf e^{T}\mathbf P\mathbf B \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \end{aligned} \]

The first term is negative definite in the tracking error. The second term has no fixed sign. Its factors may reinforce or oppose each other, so \(\dot V_e\) can be positive even though \(\mathbf A_m\) is Hurwitz. Therefore, proving nominal model stability is not enough to prove stability of the adaptive closed loop.

The key design idea is to add a parameter-error energy whose derivative can be shaped to cancel this mixed state-parameter term.

Let \(\boldsymbol\Gamma\succ 0\) be a symmetric adaptation-gain matrix. Define

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

Both terms are nonnegative and vanish simultaneously only at \(\mathbf e=\mathbf 0\) and \(\tilde{\boldsymbol\theta}=\mathbf 0\). If

\[ \boldsymbol\eta = \begin{bmatrix} \mathbf e\\ \tilde{\boldsymbol\theta} \end{bmatrix}, \]

then

\[ \underline{\alpha}\|\boldsymbol\eta\|^2 \le V \le \overline{\alpha}\|\boldsymbol\eta\|^2, \]

\[ \underline{\alpha} = \min\left\{ \lambda_{\min}(\mathbf P), \lambda_{\min}(\boldsymbol\Gamma^{-1}) \right\}, \qquad \overline{\alpha} = \max\left\{ \lambda_{\max}(\mathbf P), \lambda_{\max}(\boldsymbol\Gamma^{-1}) \right\}. \]

Hence the candidate is positive definite and radially unbounded in the combined state-parameter error vector.

Differentiate the composite candidate:

\[ \begin{aligned} \dot V &= 2\mathbf e^{T}\mathbf P\dot{\mathbf e} + 2\tilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1} \dot{\tilde{\boldsymbol\theta} }\\ &= -\mathbf e^{T}\mathbf Q\mathbf e + 2\mathbf e^{T}\mathbf P\mathbf B \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega + 2\tilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1} \dot{\tilde{\boldsymbol\theta} }. \end{aligned} \]

Since each mixed term is scalar, the cross term can also be written as

\[ 2\mathbf e^{T}\mathbf P\mathbf B \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega = 2\tilde{\boldsymbol\theta}^{T} \boldsymbol\omega\mathbf B^{T}\mathbf P\mathbf e. \]

A sufficient algebraic cancellation condition is

\[ \boldsymbol\Gamma^{-1} \dot{\tilde{\boldsymbol\theta} } = -\boldsymbol\omega\mathbf B^{T}\mathbf P\mathbf e. \]

Under that condition,

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

This equation is the bridge to Lesson 4. It is not yet a complete discussion of update-law implementation, gain selection, robustness, normalization, or projection.

flowchart TD
  E["Derive e_dot = A_m e + B theta_tilde^T omega"] --> Q["Choose Q positive definite"]
  Q --> P["Solve A_m^T P + P A_m = -Q"]
  P --> V["Choose V = e^T P e + theta_tilde^T Gamma^-1 theta_tilde"]
  V --> DV["Differentiate V"]
  DV --> C["Identify mixed cross term"]
  C --> R["Require parameter dynamics to cancel cross term"]
  R --> N["Obtain V_dot = -e^T Q e"]
        

If the parameter dynamics satisfy the cancellation condition, then \(V(t)\) is nonincreasing. Consequently, \(\mathbf e(t)\) and \(\tilde{\boldsymbol\theta}(t)\) are bounded, and

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

Therefore \(\mathbf e\in\mathcal L_2\). If the closed-loop signals are bounded so that \(\dot{\mathbf e}\) is also bounded, Barbalat's lemma gives

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

However, \(\dot V\) is only negative semidefinite in the combined vector \((\mathbf e,\tilde{\boldsymbol\theta})\). The result does not by itself imply \(\tilde{\boldsymbol\theta}(t)\to\mathbf 0\). Parameter convergence requires additional excitation conditions developed later in the course.

It is also important to distinguish a Lyapunov candidate from a completed proof. Positive definiteness of \(V\) is only the first step; the closed-loop dynamics and the resulting sign of \(\dot V\) determine the stability conclusion.

Choose the reference-model matrices

\[ \mathbf A_m= \begin{bmatrix} 0&1\\ -4&-3 \end{bmatrix}, \qquad \mathbf Q=\mathbf I_2. \]

The characteristic polynomial is \(s^2+3s+4\), whose roots have negative real parts, so \(\mathbf A_m\) is Hurwitz. Let

\[ \mathbf P= \begin{bmatrix} p_{11}&p_{12}\\ p_{12}&p_{22} \end{bmatrix}. \]

Expanding \(\mathbf A_m^T\mathbf P+\mathbf P\mathbf A_m=-\mathbf I_2\) gives

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

Solving these equations yields

\[ \boxed{ \mathbf P= \begin{bmatrix} \frac{29}{24}&\frac18\\ \frac18&\frac{5}{24} \end{bmatrix} }. \]

Sylvester's criterion confirms positive definiteness:

\[ p_{11}=\frac{29}{24}>0, \qquad \det(\mathbf P) = \frac{17}{72}>0. \]

For the numerical programs below, use

\[ \mathbf B= \begin{bmatrix}0\\1.4\end{bmatrix}, \qquad \tilde{\boldsymbol\theta} = \begin{bmatrix}0.35\\-0.20\\0.15\end{bmatrix}, \qquad \boldsymbol\Gamma=\operatorname{diag}(4,4,2). \]

The parameter error is deliberately frozen. The programs verify the identity

\[ 2\mathbf e^{T}\mathbf P\dot{\mathbf e} = -\mathbf e^{T}\mathbf Q\mathbf e+ 2\mathbf e^{T}\mathbf P\mathbf B \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \]

Because no update law is active, the composite candidate need not decrease monotonically. This is an intentional demonstration of the uncancelled cross term rather than a completed adaptive controller.

• Python: NumPy handles vectors and matrices, SciPy solves the continuous Lyapunov equation, and Matplotlib is optional for plots.
• C++: the example uses only the standard library and implements fixed-size matrix operations and fourth-order Runge-Kutta integration from scratch. Eigen and Boost.Odeint are suitable production alternatives.
• Java: the example uses only the standard library. EJML or Apache Commons Math can replace the hand-written linear algebra in larger systems.
• MATLAB/Simulink: lyap solves the matrix equation, ode45 integrates the combined model, and the script can programmatically create an equivalent Simulink block diagram.
• Wolfram Mathematica: Solve, NDSolveValue, Eigenvalues, and NMaxValue provide symbolic and numerical verification.

All implementations use the same matrices, initial conditions, reference signal, and frozen parameter error so that their numerical results can be compared directly.

Chapter6_Lesson3.py

# Chapter6_Lesson3.py
"""
MRAC for a second-order system:
derive and numerically verify the tracking-error dynamics and a composite
Lyapunov candidate. The controller-parameter error is intentionally held
constant because parameter-update laws are introduced in the next lesson.
"""

from __future__ import annotations

import math
from pathlib import Path

import numpy as np

try:
    from scipy.linalg import solve_continuous_lyapunov
except ImportError as exc:
    raise SystemExit(
        "SciPy is required. Install it with: python -m pip install numpy scipy matplotlib"
    ) from exc


AM = np.array([[0.0, 1.0], [-4.0, -3.0]])
BM = np.array([0.0, 4.0])
B = np.array([0.0, 1.4])
Q = np.eye(2)
GAMMA = np.diag([4.0, 4.0, 2.0])
THETA_TILDE = np.array([0.35, -0.20, 0.15])

T0 = 0.0
TF = 20.0
DT = 0.002
Z0 = np.array([0.0, 0.0, 0.8, -0.4])  # [x_m1, x_m2, e1, e2]


def reference(t: float) -> float:
    """Bounded command used by both the model and the regressor."""
    return math.sin(0.7 * t) + 0.25 * math.sin(1.9 * t)


def dynamics(t: float, z: np.ndarray) -> np.ndarray:
    """Combined reference-model and error dynamics."""
    xm = z[0:2]
    e = z[2:4]
    x = xm + e
    r = reference(t)
    omega = np.array([x[0], x[1], r])
    mismatch = float(THETA_TILDE @ omega)

    xm_dot = AM @ xm + BM * r
    e_dot = AM @ e + B * mismatch
    return np.concatenate((xm_dot, e_dot))


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


def simulate() -> tuple[np.ndarray, np.ndarray]:
    steps = int(round((TF - T0) / DT))
    times = np.linspace(T0, TF, steps + 1)
    states = np.zeros((steps + 1, 4))
    states[0] = Z0

    for k in range(steps):
        states[k + 1] = rk4_step(times[k], states[k], DT)

    return times, states


def main() -> None:
    # SciPy solves A X + X A^T = Q. Therefore A=AM^T and Q=-Q
    # produce AM^T P + P AM = -Q.
    p_matrix = solve_continuous_lyapunov(AM.T, -Q)
    p_matrix = 0.5 * (p_matrix + p_matrix.T)
    gamma_inv = np.linalg.inv(GAMMA)

    lyapunov_residual = AM.T @ p_matrix + p_matrix @ AM + Q
    p_eigenvalues = np.linalg.eigvalsh(p_matrix)

    times, states = simulate()
    candidate = np.zeros_like(times)
    vdot_from_gradient = np.zeros_like(times)
    vdot_from_identity = np.zeros_like(times)
    tracking_norm = np.zeros_like(times)

    parameter_energy = float(THETA_TILDE @ gamma_inv @ THETA_TILDE)

    for k, t in enumerate(times):
        xm = states[k, 0:2]
        e = states[k, 2:4]
        x = xm + e
        r = reference(float(t))
        omega = np.array([x[0], x[1], r])
        mismatch = float(THETA_TILDE @ omega)
        e_dot = AM @ e + B * mismatch

        candidate[k] = float(e @ p_matrix @ e) + parameter_energy
        vdot_from_gradient[k] = 2.0 * float(e @ p_matrix @ e_dot)
        vdot_from_identity[k] = (
            -float(e @ Q @ e)
            + 2.0 * float(e @ p_matrix @ B) * mismatch
        )
        tracking_norm[k] = np.linalg.norm(e)

    identity_error = np.max(np.abs(vdot_from_gradient - vdot_from_identity))

    print("P =")
    print(p_matrix)
    print("eigenvalues(P) =", p_eigenvalues)
    print("||AM^T P + P AM + Q||_F =", np.linalg.norm(lyapunov_residual))
    print("max Lyapunov-derivative identity error =", identity_error)
    print("initial ||e|| =", tracking_norm[0])
    print("final   ||e|| =", tracking_norm[-1])
    print(
        "Important: V need not decrease monotonically here because theta_tilde "
        "is frozen; Lesson 4 chooses parameter dynamics to cancel the cross term."
    )

    csv_path = Path("Chapter6_Lesson3_python_results.csv")
    data = np.column_stack(
        (
            times,
            states[:, 0],
            states[:, 1],
            states[:, 2],
            states[:, 3],
            tracking_norm,
            candidate,
            vdot_from_gradient,
            vdot_from_identity,
        )
    )
    np.savetxt(
        csv_path,
        data,
        delimiter=",",
        header="t,xm1,xm2,e1,e2,error_norm,V,Vdot_gradient,Vdot_identity",
        comments="",
    )
    print("saved:", csv_path.resolve())

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(times, states[:, 2], label="e1")
        plt.plot(times, states[:, 3], label="e2")
        plt.xlabel("Time (s)")
        plt.ylabel("Tracking-error states")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        plt.figure()
        plt.plot(times, candidate, label="V")
        plt.plot(times, vdot_from_identity, label="Vdot")
        plt.xlabel("Time (s)")
        plt.ylabel("Lyapunov quantities")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()
    except ImportError:
        print("Matplotlib not installed; numerical verification completed without plots.")


if __name__ == "__main__":
    main()

Chapter6_Lesson3.cpp

// Chapter6_Lesson3.cpp
// From-scratch RK4 verification of second-order MRAC error dynamics.
// The parameter error is intentionally frozen; update laws belong to Lesson 4.

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

using Vec2 = std::array<double, 2>;
using Vec3 = std::array<double, 3>;
using Vec4 = std::array<double, 4>;
using Mat2 = std::array<std::array<double, 2>, 2>;

constexpr double DT = 0.002;
constexpr double TF = 20.0;

const Mat2 AM{ { {0.0, 1.0}, {-4.0, -3.0} } };
const Mat2 Q{ { {1.0, 0.0}, {0.0, 1.0} } };
const Mat2 P{ { {29.0 / 24.0, 1.0 / 8.0}, {1.0 / 8.0, 5.0 / 24.0} } };
const Vec2 BM{ {0.0, 4.0} };
const Vec2 B{ {0.0, 1.4} };
const Vec3 THETA_TILDE{ {0.35, -0.20, 0.15} };
const Vec3 GAMMA_INV_DIAG{ {1.0 / 4.0, 1.0 / 4.0, 1.0 / 2.0} };

double reference(double t) {
    return std::sin(0.7 * t) + 0.25 * std::sin(1.9 * t);
}

Vec2 matVec(const Mat2& a, const Vec2& x) {
    return { {
        a[0][0] * x[0] + a[0][1] * x[1],
        a[1][0] * x[0] + a[1][1] * x[1]
    } };
}

double dot2(const Vec2& a, const Vec2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

double dot3(const Vec3& a, const Vec3& b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

Vec4 dynamics(double t, const Vec4& z) {
    const Vec2 xm{ {z[0], z[1]} };
    const Vec2 e{ {z[2], z[3]} };
    const Vec2 x{ {xm[0] + e[0], xm[1] + e[1]} };
    const double r = reference(t);
    const Vec3 omega{ {x[0], x[1], r} };
    const double mismatch = dot3(THETA_TILDE, omega);

    const Vec2 amXm = matVec(AM, xm);
    const Vec2 amE = matVec(AM, e);

    return { {
        amXm[0] + BM[0] * r,
        amXm[1] + BM[1] * r,
        amE[0] + B[0] * mismatch,
        amE[1] + B[1] * mismatch
    } };
}

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

Vec4 rk4Step(double t, const Vec4& z, double h) {
    const Vec4 k1 = dynamics(t, z);
    const Vec4 k2 = dynamics(t + 0.5 * h, addScaled(z, k1, 0.5 * h));
    const Vec4 k3 = dynamics(t + 0.5 * h, addScaled(z, k2, 0.5 * h));
    const Vec4 k4 = dynamics(t + h, addScaled(z, k3, h));

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

double quadratic(const Vec2& x, const Mat2& m) {
    return dot2(x, matVec(m, x));
}

double parameterEnergy() {
    double value = 0.0;
    for (std::size_t i = 0; i < 3; ++i) {
        value += THETA_TILDE[i] * GAMMA_INV_DIAG[i] * THETA_TILDE[i];
    }
    return value;
}

double lyapunovResidualNorm() {
    Mat2 residual{};
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 2; ++j) {
            double value = Q[i][j];
            for (int k = 0; k < 2; ++k) {
                value += AM[k][i] * P[k][j] + P[i][k] * AM[k][j];
            }
            residual[i][j] = value;
        }
    }

    double sum = 0.0;
    for (const auto& row : residual) {
        for (double value : row) {
            sum += value * value;
        }
    }
    return std::sqrt(sum);
}

int main() {
    Vec4 z{ {0.0, 0.0, 0.8, -0.4} };
    const int steps = static_cast<int>(std::lround(TF / DT));
    const double paramEnergy = parameterEnergy();
    double maxIdentityError = 0.0;

    std::ofstream csv("Chapter6_Lesson3_cpp_results.csv");
    if (!csv) {
        throw std::runtime_error("Could not open the CSV output file.");
    }
    csv << "t,xm1,xm2,e1,e2,error_norm,V,Vdot_gradient,Vdot_identity\n";
    csv << std::setprecision(12);

    const double initialNorm = std::hypot(z[2], z[3]);

    for (int k = 0; k <= steps; ++k) {
        const double t = k * DT;
        const Vec2 xm{ {z[0], z[1]} };
        const Vec2 e{ {z[2], z[3]} };
        const Vec2 x{ {xm[0] + e[0], xm[1] + e[1]} };
        const double r = reference(t);
        const Vec3 omega{ {x[0], x[1], r} };
        const double mismatch = dot3(THETA_TILDE, omega);
        const Vec2 eDotBase = matVec(AM, e);
        const Vec2 eDot{ {
            eDotBase[0] + B[0] * mismatch,
            eDotBase[1] + B[1] * mismatch
        } };

        const double v = quadratic(e, P) + paramEnergy;
        const double vdotGradient = 2.0 * dot2(e, matVec(P, eDot));
        const double ePb = dot2(e, matVec(P, B));
        const double vdotIdentity = -quadratic(e, Q) + 2.0 * ePb * mismatch;
        maxIdentityError = std::max(
            maxIdentityError, std::abs(vdotGradient - vdotIdentity)
        );

        csv << t << ',' << z[0] << ',' << z[1] << ',' << z[2] << ','
            << z[3] << ',' << std::hypot(z[2], z[3]) << ',' << v << ','
            << vdotGradient << ',' << vdotIdentity << '\n';

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

    std::cout << "P = [[29/24, 1/8], [1/8, 5/24]]\n";
    std::cout << "||AM^T P + P AM + Q||_F = "
              << lyapunovResidualNorm() << '\n';
    std::cout << "max Lyapunov-derivative identity error = "
              << maxIdentityError << '\n';
    std::cout << "initial ||e|| = " << initialNorm << '\n';
    std::cout << "final   ||e|| = " << std::hypot(z[2], z[3]) << '\n';
    std::cout << "V is not guaranteed to decrease because theta_tilde is frozen.\n";
    return 0;
}

Chapter6_Lesson3.java

// Chapter6_Lesson3.java
// From-scratch RK4 verification of second-order MRAC error dynamics.
// The parameter error remains constant because update laws are Lesson 4 material.

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

public final class Chapter6_Lesson3 {
    private static final double DT = 0.002;
    private static final double TF = 20.0;

    private static final double[][] AM = { {0.0, 1.0}, {-4.0, -3.0} };
    private static final double[][] Q = { {1.0, 0.0}, {0.0, 1.0} };
    private static final double[][] P = {
        {29.0 / 24.0, 1.0 / 8.0},
        {1.0 / 8.0, 5.0 / 24.0}
    };
    private static final double[] BM = {0.0, 4.0};
    private static final double[] B = {0.0, 1.4};
    private static final double[] THETA_TILDE = {0.35, -0.20, 0.15};
    private static final double[] GAMMA_INV_DIAG = {0.25, 0.25, 0.5};

    private Chapter6_Lesson3() {
    }

    private static double reference(double t) {
        return Math.sin(0.7 * t) + 0.25 * Math.sin(1.9 * t);
    }

    private static double[] matVec(double[][] a, double[] x) {
        return new double[] {
            a[0][0] * x[0] + a[0][1] * x[1],
            a[1][0] * x[0] + a[1][1] * x[1]
        };
    }

    private static double dot(double[] a, double[] b) {
        double value = 0.0;
        for (int i = 0; i < a.length; i++) {
            value += a[i] * b[i];
        }
        return value;
    }

    private static double[] dynamics(double t, double[] z) {
        double[] xm = {z[0], z[1]};
        double[] e = {z[2], z[3]};
        double[] x = {xm[0] + e[0], xm[1] + e[1]};
        double r = reference(t);
        double[] omega = {x[0], x[1], r};
        double mismatch = dot(THETA_TILDE, omega);

        double[] amXm = matVec(AM, xm);
        double[] amE = matVec(AM, e);

        return new double[] {
            amXm[0] + BM[0] * r,
            amXm[1] + BM[1] * r,
            amE[0] + B[0] * mismatch,
            amE[1] + B[1] * mismatch
        };
    }

    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 h) {
        double[] k1 = dynamics(t, z);
        double[] k2 = dynamics(t + 0.5 * h, addScaled(z, k1, 0.5 * h));
        double[] k3 = dynamics(t + 0.5 * h, addScaled(z, k2, 0.5 * h));
        double[] k4 = dynamics(t + h, addScaled(z, k3, h));

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

    private static double quadratic(double[] x, double[][] m) {
        return dot(x, matVec(m, x));
    }

    private static double parameterEnergy() {
        double value = 0.0;
        for (int i = 0; i < THETA_TILDE.length; i++) {
            value += THETA_TILDE[i] * GAMMA_INV_DIAG[i] * THETA_TILDE[i];
        }
        return value;
    }

    private static double lyapunovResidualNorm() {
        double sumSquares = 0.0;
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                double value = Q[i][j];
                for (int k = 0; k < 2; k++) {
                    value += AM[k][i] * P[k][j] + P[i][k] * AM[k][j];
                }
                sumSquares += value * value;
            }
        }
        return Math.sqrt(sumSquares);
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.ROOT);
        double[] z = {0.0, 0.0, 0.8, -0.4};
        int steps = (int) Math.round(TF / DT);
        double paramEnergy = parameterEnergy();
        double maxIdentityError = 0.0;
        double initialNorm = Math.hypot(z[2], z[3]);

        Path csvPath = Path.of("Chapter6_Lesson3_java_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(csvPath)) {
            writer.write("t,xm1,xm2,e1,e2,error_norm,V,Vdot_gradient,Vdot_identity");
            writer.newLine();

            for (int k = 0; k <= steps; k++) {
                double t = k * DT;
                double[] xm = {z[0], z[1]};
                double[] e = {z[2], z[3]};
                double[] x = {xm[0] + e[0], xm[1] + e[1]};
                double r = reference(t);
                double[] omega = {x[0], x[1], r};
                double mismatch = dot(THETA_TILDE, omega);

                double[] eDotBase = matVec(AM, e);
                double[] eDot = {
                    eDotBase[0] + B[0] * mismatch,
                    eDotBase[1] + B[1] * mismatch
                };

                double v = quadratic(e, P) + paramEnergy;
                double vdotGradient = 2.0 * dot(e, matVec(P, eDot));
                double ePb = dot(e, matVec(P, B));
                double vdotIdentity = -quadratic(e, Q) + 2.0 * ePb * mismatch;
                maxIdentityError = Math.max(
                    maxIdentityError, Math.abs(vdotGradient - vdotIdentity)
                );

                writer.write(String.format(
                    Locale.ROOT,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                    t, z[0], z[1], z[2], z[3], Math.hypot(z[2], z[3]),
                    v, vdotGradient, vdotIdentity
                ));
                writer.newLine();

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

        System.out.println("P = [[29/24, 1/8], [1/8, 5/24]]");
        System.out.println("||AM^T P + P AM + Q||_F = " + lyapunovResidualNorm());
        System.out.println(
            "max Lyapunov-derivative identity error = " + maxIdentityError
        );
        System.out.println("initial ||e|| = " + initialNorm);
        System.out.println("final   ||e|| = " + Math.hypot(z[2], z[3]));
        System.out.println(
            "V is not guaranteed to decrease because theta_tilde is frozen."
        );
        System.out.println("saved: " + csvPath.toAbsolutePath());
    }
}

The script first verifies the equations with ode45. Set buildSimulinkModel = true to create Chapter6_Lesson3_Simulink.slx using State-Space, Sum, Mux, Demux, Gain, and To Workspace blocks.

Chapter6_Lesson3.m

% Chapter6_Lesson3.m
% MRAC second-order error dynamics and Lyapunov-candidate verification.
% Parameter error is deliberately held constant; adaptation laws are Lesson 4.

clear; clc; close all;

Am = [0 1; -4 -3];
Bm = [0; 4];
B = [0; 1.4];
Q = eye(2);
Gamma = diag([4 4 2]);
thetaTilde = [0.35; -0.20; 0.15];

if exist('lyap', 'file') == 2
    % lyap(A,Q) solves A*P + P*A' + Q = 0.
    P = lyap(Am', Q);
else
    warning('Control System Toolbox not found; using the exact P for this example.');
    P = [29/24, 1/8; 1/8, 5/24];
end
P = 0.5*(P + P');

z0 = [0; 0; 0.8; -0.4]; % [xm1; xm2; e1; e2]
tspan = linspace(0, 20, 10001);
[t, z] = ode45(@combinedDynamics, tspan, z0);

xm = z(:, 1:2);
e = z(:, 3:4);
x = xm + e;
r = sin(0.7*t) + 0.25*sin(1.9*t);
omega = [x(:,1), x(:,2), r];
mismatch = omega*thetaTilde;

V = zeros(size(t));
VdotGradient = zeros(size(t));
VdotIdentity = zeros(size(t));
parameterEnergy = thetaTilde'*(Gamma\thetaTilde);

for k = 1:numel(t)
    ek = e(k,:)';
    eDot = Am*ek + B*mismatch(k);
    V(k) = ek'*P*ek + parameterEnergy;
    VdotGradient(k) = 2*ek'*P*eDot;
    VdotIdentity(k) = -ek'*Q*ek + 2*(ek'*P*B)*mismatch(k);
end

fprintf('P =\n');
disp(P);
fprintf('eig(P) =\n');
disp(eig(P));
fprintf('||Am''P + P Am + Q||_F = %.3e\n', norm(Am'*P + P*Am + Q, 'fro'));
fprintf('max derivative identity error = %.3e\n', ...
    max(abs(VdotGradient - VdotIdentity)));
fprintf('initial ||e|| = %.6f\n', norm(e(1,:)));
fprintf('final   ||e|| = %.6f\n', norm(e(end,:)));
fprintf(['V is not guaranteed to decrease because thetaTilde is frozen. ', ...
    'Lesson 4 supplies parameter dynamics that cancel the cross term.\n']);

results = table(t, xm(:,1), xm(:,2), e(:,1), e(:,2), vecnorm(e,2,2), ...
    V, VdotGradient, VdotIdentity, ...
    'VariableNames', {'t','xm1','xm2','e1','e2','errorNorm', ...
    'V','VdotGradient','VdotIdentity'});
writetable(results, 'Chapter6_Lesson3_matlab_results.csv');

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

figure;
plot(t, V, 'LineWidth', 1.2); hold on;
plot(t, VdotIdentity, 'LineWidth', 1.2);
grid on; xlabel('Time (s)'); ylabel('Lyapunov quantities');
legend('V','dV/dt');

% Set true to generate an equivalent Simulink block model.
buildSimulinkModel = false;
if buildSimulinkModel
    buildLesson3SimulinkModel(Am, Bm, B, thetaTilde);
end

function dz = combinedDynamics(t, z)
    Am = [0 1; -4 -3];
    Bm = [0; 4];
    B = [0; 1.4];
    thetaTilde = [0.35; -0.20; 0.15];

    xm = z(1:2);
    e = z(3:4);
    x = xm + e;
    r = sin(0.7*t) + 0.25*sin(1.9*t);
    omega = [x(1); x(2); r];
    mismatch = thetaTilde'*omega;

    dz = [Am*xm + Bm*r; Am*e + B*mismatch];
end

function buildLesson3SimulinkModel(Am, Bm, B, thetaTilde)
    if ~license('test', 'Simulink')
        error('Simulink is not licensed on this MATLAB installation.');
    end

    model = 'Chapter6_Lesson3_Simulink';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    if isfile([model '.slx'])
        delete([model '.slx']);
    end

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

    add_block('simulink/Sources/Sine Wave', [model '/r1'], ...
        'Amplitude', '1', 'Frequency', '0.7', ...
        'Position', [30 45 75 75]);
    add_block('simulink/Sources/Sine Wave', [model '/r2'], ...
        'Amplitude', '0.25', 'Frequency', '1.9', ...
        'Position', [30 105 75 135]);
    add_block('simulink/Math Operations/Sum', [model '/r'], ...
        'Inputs', '++', 'Position', [110 60 135 120]);

    add_block('simulink/Continuous/State-Space', [model '/Reference Model'], ...
        'A', 'Am', 'B', 'Bm', 'C', 'eye(2)', 'D', 'zeros(2,1)', ...
        'X0', '[0;0]', 'Position', [190 45 310 105]);

    add_block('simulink/Continuous/State-Space', [model '/Error Dynamics'], ...
        'A', 'Am', 'B', 'B', 'C', 'eye(2)', 'D', 'zeros(2,1)', ...
        'X0', '[0.8;-0.4]', 'Position', [560 175 680 235]);

    add_block('simulink/Math Operations/Sum', [model '/x=xm+e'], ...
        'Inputs', '++', 'Position', [350 85 380 145]);
    add_block('simulink/Signal Routing/Demux', [model '/Demux x'], ...
        'Outputs', '2', 'Position', [415 80 420 150]);
    add_block('simulink/Signal Routing/Mux', [model '/omega'], ...
        'Inputs', '3', 'Position', [470 75 475 165]);
    add_block('simulink/Math Operations/Gain', [model '/thetaTilde dot omega'], ...
        'Gain', 'thetaTilde''', 'Multiplication', 'Matrix(K*u)', ...
        'Position', [500 95 545 145]);

    add_block('simulink/Sinks/To Workspace', [model '/xm workspace'], ...
        'VariableName', 'xm_sim', 'SaveFormat', 'Structure With Time', ...
        'Position', [350 20 450 50]);
    add_block('simulink/Sinks/To Workspace', [model '/e workspace'], ...
        'VariableName', 'e_sim', 'SaveFormat', 'Structure With Time', ...
        'Position', [720 175 820 205]);

    add_line(model, 'r1/1', 'r/1');
    add_line(model, 'r2/1', 'r/2');
    add_line(model, 'r/1', 'Reference Model/1');
    add_line(model, 'Reference Model/1', 'x=xm+e/1');
    add_line(model, 'Error Dynamics/1', 'x=xm+e/2');
    add_line(model, 'x=xm+e/1', 'Demux x/1');
    add_line(model, 'Demux x/1', 'omega/1');
    add_line(model, 'Demux x/2', 'omega/2');
    add_line(model, 'r/1', 'omega/3');
    add_line(model, 'omega/1', 'thetaTilde dot omega/1');
    add_line(model, 'thetaTilde dot omega/1', 'Error Dynamics/1');
    add_line(model, 'Reference Model/1', 'xm workspace/1');
    add_line(model, 'Error Dynamics/1', 'e workspace/1');

    save_system(model);
    fprintf('Generated %s.slx\n', model);
end

Chapter6_Lesson3.nb

Notebook[{
Cell["Chapter 6, Lesson 3: Error Dynamics and Lyapunov Candidate", "Title"],
Cell["The parameter error is held constant so that the notebook verifies the Lesson 3 identities without introducing the Lesson 4 update law.", "Text"],
Cell[BoxData["ClearAll[\"Global`*\"];

Am = { {0., 1.}, {-4., -3.} };
Bm = {0., 4.};
B = {0., 1.4};
Q = IdentityMatrix[2];
Gamma = DiagonalMatrix[{4., 4., 2.}];
thetaTilde = {0.35, -0.20, 0.15};

pSolution = First@Solve[
   Transpose[Am].{ {p11, p12}, {p12, p22} } +
      { {p11, p12}, {p12, p22} }.Am == -Q,
   {p11, p12, p22}
   ];
P = { {p11, p12}, {p12, p22} } /. pSolution // N;

reference[t_] := Sin[0.7 t] + 0.25 Sin[1.9 t];

solution = NDSolveValue[
   {
    xm1'[t] == xm2[t],
    xm2'[t] == -4 xm1[t] - 3 xm2[t] + 4 reference[t],
    e1'[t] == e2[t],
    e2'[t] == -4 e1[t] - 3 e2[t] +
      1.4 thetaTilde.{
        xm1[t] + e1[t],
        xm2[t] + e2[t],
        reference[t]
        },
    xm1[0] == 0, xm2[0] == 0,
    e1[0] == 0.8, e2[0] == -0.4
    },
   {xm1, xm2, e1, e2},
   {t, 0, 20}
   ];

xm[t_] := {solution[[1]][t], solution[[2]][t]};
e[t_] := {solution[[3]][t], solution[[4]][t]};
x[t_] := xm[t] + e[t];
omega[t_] := {x[t][[1]], x[t][[2]], reference[t]};
mismatch[t_] := thetaTilde.omega[t];
eDot[t_] := Am.e[t] + B mismatch[t];

parameterEnergy = thetaTilde.Inverse[Gamma].thetaTilde;
V[t_] := e[t].P.e[t] + parameterEnergy;
VdotGradient[t_] := 2 e[t].P.eDot[t];
VdotIdentity[t_] := -e[t].Q.e[t] +
   2 (e[t].P.B) mismatch[t];

Print[\"P = \", MatrixForm[P]];
Print[\"Eigenvalues(P) = \", Eigenvalues[P]];
Print[\"Lyapunov residual norm = \",
  Norm[Transpose[Am].P + P.Am + Q, \"Frobenius\"]];
Print[\"Maximum derivative identity error = \",
  NMaxValue[
   {Abs[VdotGradient[t] - VdotIdentity[t]], 0 <= t <= 20},
   t
   ]];

errorPlot = Plot[
   Evaluate[{solution[[3]][t], solution[[4]][t]}],
   {t, 0, 20},
   PlotLegends -> {\"e1\", \"e2\"},
   AxesLabel -> {\"t\", \"tracking error\"},
   PlotRange -> All
   ];

lyapunovPlot = Plot[
   Evaluate[{V[t], VdotIdentity[t]}],
   {t, 0, 20},
   PlotLegends -> {\"V\", \"Vdot\"},
   AxesLabel -> {\"t\", \"Lyapunov quantities\"},
   PlotRange -> All
   ];

Column[{errorPlot, lyapunovPlot}]
"], "Input"]
},
WindowSize->{1200, 800},
StyleDefinitions->"Default.nb"]

Common error 1: subtracting unmatched dynamics.

The compact error equation contains \(\mathbf A_m\mathbf e\) only after the ideal matching equations have been applied. Subtracting the plant and model before introducing the ideal parameters often leaves unexplained terms.

Common error 2: reversing the parameter-error definition.

This lesson uses \(\tilde{\boldsymbol\theta}=\boldsymbol\theta- \boldsymbol\theta^{*}\). If the opposite convention is used, the mismatch term and the later update-law sign must both change.

Common error 3: choosing an arbitrary positive matrix.

The matrix \(\mathbf P\) must satisfy the Lyapunov equation associated with the chosen \(\mathbf A_m\) and \(\mathbf Q\). Positive definiteness alone does not produce the required derivative identity.

Common error 4: claiming parameter convergence from tracking convergence.

The semidefinite derivative establishes boundedness and, under standard regularity conditions, tracking-error convergence. It does not establish convergence of every controller parameter without sufficient excitation.

Common error 5: treating a candidate as a proof.

The expression for \(V\) is a candidate until its derivative is evaluated along the complete closed-loop system and shown to have the required sign.

Problem 1 (Ideal Matching Parameters): Consider \(\ddot y+1.2\dot y+2.5y=1.4u\) and the reference model \(\ddot y_m+3\dot y_m+4y_m=4r\). Determine \(\theta_1^{*}\), \(\theta_2^{*}\), and \(\theta_r^{*}\).

Solution:

\[ \theta_1^{*} = \frac{2.5-4}{1.4} = -\frac{15}{14}, \qquad \theta_2^{*} = \frac{1.2-3}{1.4} = -\frac97, \qquad \theta_r^{*} = \frac4{1.4} = \frac{20}{7}. \]

These values would exactly reproduce the selected reference-model coefficients if the plant parameters were known.

Problem 2 (Coordinate Error Dynamics): Starting from \(\dot{\mathbf e}=\mathbf A_m\mathbf e+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega\), derive the two scalar equations.

Solution:

\[ \begin{bmatrix}\dot e_1\\\dot e_2\end{bmatrix} = \begin{bmatrix} 0&1\\-a_{m0}&-a_{m1} \end{bmatrix} \begin{bmatrix}e_1\\e_2\end{bmatrix} + \begin{bmatrix}0\\b\end{bmatrix} \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \]

Therefore,

\[ \dot e_1=e_2, \qquad \dot e_2=-a_{m0}e_1-a_{m1}e_2+ b\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega. \]

Problem 3 (Solve a Lyapunov Equation): For \(\mathbf A_m=\begin{bmatrix}0&1\\-4&-3\end{bmatrix}\) and \(\mathbf Q=\mathbf I_2\), solve for the symmetric matrix \(\mathbf P\).

Solution:

Substitute \(\mathbf P=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}\) into the Lyapunov equation. Equality of the matrix entries gives

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

Thus \(p_{12}=1/8\), \(p_{22}=5/24\), and \(p_{11}=29/24\). Therefore

\[ \mathbf P= \begin{bmatrix} 29/24&1/8\\ 1/8&5/24 \end{bmatrix}\succ0. \]

Problem 4 (Positive-Definiteness Bounds): Prove that the composite candidate is bounded above and below by quadratic functions of \(\boldsymbol\eta=[\mathbf e^T,\tilde{\boldsymbol\theta}^T]^T\).

Solution:

For every symmetric positive-definite matrix \(\mathbf M\),

\[ \lambda_{\min}(\mathbf M)\|\mathbf z\|^2 \le \mathbf z^T\mathbf M\mathbf z \le \lambda_{\max}(\mathbf M)\|\mathbf z\|^2. \]

Apply this inequality separately to the state and parameter terms, then combine the lower bounds using the smaller eigenvalue and the upper bounds using the larger eigenvalue. This gives

\[ \underline{\alpha}\|\boldsymbol\eta\|^2 \le V\le \overline{\alpha}\|\boldsymbol\eta\|^2, \]

with \(\underline{\alpha}\) and \(\overline{\alpha}\) defined in Section 7.

Problem 5 (Why the Cross Term Matters): Let \(\mathbf e=[0,0.1]^T\), \(\mathbf Q=\mathbf I_2\), \(\mathbf P\) be the matrix from Problem 3, \(\mathbf B=[0,1.4]^T\), \(\tilde{\boldsymbol\theta}=[0.35,-0.20,0.15]^T\), and \(\boldsymbol\omega=[10,10,10]^T\). Evaluate \(\dot V_e\) when the parameter error is frozen.

Solution:

\[ \tilde{\boldsymbol\theta}^{T}\boldsymbol\omega = 3, \qquad \mathbf e^{T}\mathbf P\mathbf B = \frac{7}{240}. \]

Therefore,

\[ \dot V_e = -0.1^2+ 2\left(\frac{7}{240}\right)(3) = -0.01+0.175 = 0.165. \]

The derivative is positive at this point. This explicitly demonstrates that the state-only Lyapunov function does not prove stability in the presence of parameter mismatch.

Problem 6 (Conditional Tracking Convergence): Suppose the parameter dynamics cancel the cross term so that \(\dot V=-\mathbf e^T\mathbf Q\mathbf e\). Show that \(\mathbf e\in\mathcal L_2\), and state the additional condition needed to conclude \(\mathbf e(t)\to\mathbf 0\).

Solution:

Integrating from zero to \(T\) gives

\[ V(T)-V(0) = -\int_0^T \mathbf e^T(t)\mathbf Q\mathbf e(t)\,dt. \]

Since \(V(T)\ge0\),

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

so \(\mathbf e\in\mathcal L_2\). If \(\dot{\mathbf e}\) is bounded, then the tracking error is uniformly continuous. Barbalat's lemma then yields \(\mathbf e(t)\to\mathbf 0\).

For the second-order MRAC problem, the matching equations define a constant ideal controller vector. Subtracting the reference model from the matched plant yields \(\dot{\mathbf e}=\mathbf A_m\mathbf e+ \mathbf B\tilde{\boldsymbol\theta}^{T}\boldsymbol\omega\). A state-only quadratic function reveals a sign-indefinite cross term, motivating the composite candidate \(V=\mathbf e^T\mathbf P\mathbf e+ \tilde{\boldsymbol\theta}^T\boldsymbol\Gamma^{-1} \tilde{\boldsymbol\theta}\). The Lyapunov equation determines \(\mathbf P\), while the derivative of the parameter-error energy provides the mechanism needed to cancel the mixed term. Lesson 4 converts this cancellation requirement into explicit parameter-update laws.