After completing this lesson, students should be able to:

• write a linearly parameterized control law for a second-order SISO plant;
• derive the ideal matching gains from coefficient matching;
• construct the scalar adaptation signal \( s=e^{T}PB \);
• derive vector and componentwise Lyapunov update laws;
• prove boundedness of the tracking and parameter errors under the stated assumptions;
• implement the continuous-time adaptive law without confusing numerical integration with a discrete-time MRAC design.

Required prior material is limited to linear state-space models, controllable canonical form, Hurwitz matrices, Lyapunov equations, and the combined state-parameter Lyapunov candidate introduced in the preceding lesson.

Consider the uncertain plant

\[ \dot{x}_1=x_2, \qquad \dot{x}_2=-a_1x_1-a_2x_2+b u, \]

where \(a_1\), \(a_2\), and the magnitude \(|b|\) are unknown. The sign of \(b\) is assumed known and nonzero. The desired model is

\[ \dot{x}_{m1}=x_{m2}, \qquad \dot{x}_{m2}=-a_{m1}x_{m1}-a_{m2}x_{m2}+b_m r, \]

with \(a_{m1}>0\) and \(a_{m2}>0\), so \(A_m\) is Hurwitz. Use the direct adaptive controller

\[ u=\hat{k}_1x_1+\hat{k}_2x_2+\hat{k}_r r =\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}, \]

\[ \hat{\boldsymbol{\theta}} =\begin{bmatrix}\hat{k}_1&\hat{k}_2&\hat{k}_r\end{bmatrix}^{T}, \qquad \boldsymbol{\phi} =\begin{bmatrix}x_1&x_2&r\end{bmatrix}^{T}. \]

flowchart TD
  R["Command r"] --> RM["Reference model"]
  R --> REG["Regressor: x1, x2, r"]
  P["Plant states: x1, x2"] --> REG
  REG --> C["Adaptive control: u = khat dot phi"]
  C --> P
  RM --> E["Tracking error: e = x - xm"]
  P --> E
  E --> S["Adaptation signal: s = eT P B"]
  S --> U["Parameter update laws"]
  REG --> U
  U --> C
        

If constant ideal gains exist, substituting \(u^{*}=k_1^{*}x_1+k_2^{*}x_2+k_r^{*}r\) into the plant and matching coefficients with the reference model gives

\[ -a_1+b k_1^{*}=-a_{m1}, \qquad -a_2+b k_2^{*}=-a_{m2}, \qquad b k_r^{*}=b_m. \]

\[ k_1^{*}=\frac{a_1-a_{m1}}{b}, \qquad k_2^{*}=\frac{a_2-a_{m2}}{b}, \qquad k_r^{*}=\frac{b_m}{b}. \]

Define tracking and parameter errors by

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

With \(B=\begin{bmatrix}0&1\end{bmatrix}^{T}\), the closed-loop tracking-error dynamics become

\[ \dot{\mathbf{e}} =A_m\mathbf{e}+B b\widetilde{\boldsymbol{\theta}}^{T} \boldsymbol{\phi}, \qquad A_m=\begin{bmatrix}0&1\\-a_{m1}&-a_{m2}\end{bmatrix}. \]

This form isolates all uncertainty in the scalar product \(\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}\), which is precisely what enables a three-parameter update law.

Choose any symmetric matrix \(Q=Q^{T}\succ0\). Since \(A_m\) is Hurwitz, the Lyapunov equation

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

has a unique symmetric positive-definite solution \(P=P^{T}\succ0\). Define the scalar signal

\[ s=\mathbf{e}^{T}PB=B^{T}P\mathbf{e}. \]

For the second-order canonical form and \(P=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}\),

\[ s=p_{12}e_1+p_{22}e_2. \]

Thus, the adaptive law does not use an arbitrary weighted error. Its weighting is generated by the reference-model dynamics through the Lyapunov equation.

Closed-form solution for \(Q=I_2\)

Let

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

Equating entries of \(A_m^{T}P+PA_m=-I_2\) yields

\[ -2a_{m1}p_{12}=-1, \qquad 2p_{12}-2a_{m2}p_{22}=-1, \]

\[ p_{11}-a_{m2}p_{12}-a_{m1}p_{22}=0. \]

Therefore,

\[ p_{12}=\frac{1}{2a_{m1}}, \qquad p_{22}=\frac{1+a_{m1}}{2a_{m1}a_{m2}}, \]

\[ p_{11}=\frac{a_{m2}}{2a_{m1}} +\frac{1+a_{m1}}{2a_{m2}}. \]

Let \(\Gamma=\Gamma^{T}\succ0\) be the adaptation-gain matrix. Consider the Lyapunov candidate

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

Since the ideal parameter vector is constant, \(\dot{\widetilde{\boldsymbol{\theta}}} =\dot{\hat{\boldsymbol{\theta}}}\). Differentiate \(V\) along the error dynamics:

\[ \begin{aligned} \dot{V} &=\dot{\mathbf{e}}^{T}P\mathbf{e} +\mathbf{e}^{T}P\dot{\mathbf{e}} +2|b|\widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\hat{\boldsymbol{\theta}}}\\ &=\mathbf{e}^{T}(A_m^{T}P+PA_m)\mathbf{e} +2b\,\mathbf{e}^{T}PB\, \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +2|b|\widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\hat{\boldsymbol{\theta}}}. \end{aligned} \]

Using the Lyapunov equation and \(s=\mathbf{e}^{T}PB\),

\[ \dot{V}=-\mathbf{e}^{T}Q\mathbf{e} +2b s\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +2|b|\widetilde{\boldsymbol{\theta}}^{T} \Gamma^{-1}\dot{\hat{\boldsymbol{\theta}}}. \]

Choose the update law so that the two terms containing the parameter error cancel:

\[ \boxed{ \dot{\hat{\boldsymbol{\theta}}} =-\Gamma\boldsymbol{\phi}\,s\,\operatorname{sgn}(b) } \]

Indeed,

\[ \begin{aligned} 2|b|\widetilde{\boldsymbol{\theta}}^{T}\Gamma^{-1} \dot{\hat{\boldsymbol{\theta}}} &=-2|b|\operatorname{sgn}(b) \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}s\\ &=-2b\widetilde{\boldsymbol{\theta}}^{T} \boldsymbol{\phi}s. \end{aligned} \]

Consequently,

\[ \boxed{\dot{V}=-\mathbf{e}^{T}Q\mathbf{e}\le0}. \]

For diagonal adaptation gains

\[ \Gamma=\operatorname{diag}(\gamma_1,\gamma_2,\gamma_r), \qquad \gamma_1,\gamma_2,\gamma_r>0, \]

the vector update law expands to

\[ \boxed{\dot{\hat{k}}_1 =-\gamma_1 x_1 s\operatorname{sgn}(b)}, \]

\[ \boxed{\dot{\hat{k}}_2 =-\gamma_2 x_2 s\operatorname{sgn}(b)}, \]

\[ \boxed{\dot{\hat{k}}_r =-\gamma_r r s\operatorname{sgn}(b)}. \]

Every update has the same structure: adaptation gain multiplied by a regressor component and by the common Lyapunov error signal. When \(b>0\), the sign factor is one. When \(b<0\) but its sign is known, all three update directions reverse.

Non-diagonal \(\Gamma\)

A full symmetric positive-definite matrix is also permitted:

\[ \dot{\hat{\boldsymbol{\theta}}} =-\Gamma\boldsymbol{\phi}s\operatorname{sgn}(b). \]

Off-diagonal entries couple parameter updates. For this introductory second-order design, a diagonal matrix is preferable because each gain has a transparent tuning coefficient and the stability proof remains identical.

Because \(V(t)\le V(0)\) and both terms in \(V\) are positive definite, the tracking error and parameter error are bounded:

\[ \mathbf{e}\in\mathcal{L}_{\infty}, \qquad \widetilde{\boldsymbol{\theta}}\in\mathcal{L}_{\infty}. \]

Furthermore,

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

so \(\mathbf{e}\in\mathcal{L}_{2}\). Under bounded reference and closed-loop signals, the error derivative is bounded. Barbalat's lemma then gives

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

The result establishes asymptotic tracking, not automatic recovery of the physical plant coefficients or convergence of \(\hat{\boldsymbol{\theta}}\) to one unique vector. The update law stops changing whenever the products \(\boldsymbol{\phi}s\) vanish. Conditions for parameter convergence are treated later in the course.

The positive constants \(\gamma_1\), \(\gamma_2\), and \(\gamma_r\) determine how strongly each parameter reacts to \(\phi_i s\). Scaling all gains by \(c>1\) makes the differential update faster for the same instantaneous signals, but it does not change the algebraic cancellation in the Lyapunov proof.

Since the three regressors may have different physical units, the numerical values of the adaptation gains generally have different units. A dimensionally correct implementation checks that each right-hand side has the units of its corresponding parameter derivative:

\[ [\gamma_i][\phi_i][s] =[\dot{\hat{\theta}}_i]. \]

The present lesson uses moderate diagonal gains and exact state measurements. Noise-sensitive tuning and additional practical modifications are intentionally deferred to later chapters.

flowchart TD
  A["Choose stable reference model Am, Bm"] --> B["Select positive definite Q"]
  B --> C["Solve AmT P + P Am = -Q"]
  C --> D["Measure x and compute model state xm"]
  D --> E["Compute e = x - xm"]
  E --> F["Compute s = eT P B"]
  F --> G["Form phi = x1, x2, r"]
  G --> H["Update khatDot = -Gamma phi s sign(b)"]
  H --> I["Compute u = khat dot phi"]
  I --> J["Integrate plant, model, and gain states"]
  J --> D
        

In software, the estimated gains are dynamic states. Therefore, the complete simulation state is

\[ \mathbf{z} =\begin{bmatrix}x_1&x_2&x_{m1}&x_{m2}& \hat{k}_1&\hat{k}_2&\hat{k}_r\end{bmatrix}^{T}. \]

A continuous ODE solver integrates all seven differential equations simultaneously. A small numerical step does not convert this controller into a formally derived discrete-time adaptive controller; it only approximates the continuous-time equations.

Use the uncertain plant parameters

\[ a_1=1, \qquad a_2=0.8, \qquad b=1.5, \]

and the reference model

\[ a_{m1}=2.25, \qquad a_{m2}=2.10, \qquad b_m=2.25. \]

For \(Q=I_2\),

\[ P=\begin{bmatrix} 1.24047619&0.22222222\\ 0.22222222&0.34391534 \end{bmatrix}. \]

The ideal matching gains are

\[ \boldsymbol{\theta}^{*} =\begin{bmatrix}-0.83333333&-0.86666667&1.5\end{bmatrix}^{T}. \]

The implementations use \(\Gamma=\operatorname{diag}(4,4,2)\), zero initial states and gains, and a unit step applied at \(t=0.5\) seconds. The expected outcome is bounded gain evolution and a tracking error that approaches zero. The estimated gains need not numerically equal the ideal vector for this command.

Software libraries and alternatives

• Python: NumPy for arrays, Matplotlib for plots, and optionally SciPy's continuous Lyapunov solver. The provided file uses the closed-form \(2\times2\) solution and a from-scratch RK4 integrator.
• C++: the standard library is sufficient; Eigen can be substituted for larger matrix models.
• Java: the standard library is sufficient; Apache Commons Math can be used for general matrix and ODE operations.
• MATLAB: ode45 performs integration, while lyap can solve a general continuous Lyapunov equation.
• Simulink: integrator and MATLAB Function blocks can realize the equations; supported installations may also provide a dedicated Model Reference Adaptive Control block.
• Wolfram Mathematica: NDSolveValue integrates the equations and LyapunovSolve handles general Lyapunov equations.

Chapter6_Lesson4.py

"""Chapter 6, Lesson 4: Parameter update laws for second-order MRAC.

Continuous-time direct MRAC simulation for a second-order SISO plant.
The differential equations are integrated with a fixed-step RK4 method.
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    # Unknown plant coefficients used only by the simulator.
    a1: float = 1.0
    a2: float = 0.8
    b: float = 1.5

    # Desired model: x_m'' + am2*x_m' + am1*x_m = bm*r.
    am1: float = 2.25
    am2: float = 2.10
    bm: float = 2.25

    # Diagonal adaptation-gain matrix Gamma.
    gamma1: float = 4.0
    gamma2: float = 4.0
    gammar: float = 2.0

    dt: float = 0.002
    final_time: float = 20.0


def reference(t: float) -> float:
    """Bounded command used in the numerical experiment."""
    return 0.0 if t < 0.5 else 1.0


def lyapunov_matrix_for_q_identity(am1: float, am2: float) -> np.ndarray:
    """Return P solving A_m.T P + P A_m = -I for A_m=[[0,1],[-am1,-am2]]."""
    if am1 <= 0.0 or am2 <= 0.0:
        raise ValueError("am1 and am2 must be positive for a Hurwitz second-order model")

    p12 = 1.0 / (2.0 * am1)
    p22 = (1.0 + am1) / (2.0 * am1 * am2)
    p11 = am2 / (2.0 * am1) + (1.0 + am1) / (2.0 * am2)
    return np.array([[p11, p12], [p12, p22]], dtype=float)


def derivative(t: float, z: np.ndarray, p: Parameters, P: np.ndarray) -> np.ndarray:
    """Closed-loop plant, model, and parameter-update differential equations."""
    x1, x2, xm1, xm2, khat1, khat2, khatr = z
    r = reference(t)

    phi = np.array([x1, x2, r], dtype=float)
    khat = np.array([khat1, khat2, khatr], dtype=float)
    u = float(khat @ phi)

    e = np.array([x1 - xm1, x2 - xm2], dtype=float)
    B = np.array([0.0, 1.0], dtype=float)
    s = float(e @ P @ B)
    sign_b = 1.0 if p.b > 0.0 else -1.0

    gamma = np.diag([p.gamma1, p.gamma2, p.gammar])
    khat_dot = -gamma @ phi * s * sign_b

    return np.array(
        [
            x2,
            -p.a1 * x1 - p.a2 * x2 + p.b * u,
            xm2,
            -p.am1 * xm1 - p.am2 * xm2 + p.bm * r,
            khat_dot[0],
            khat_dot[1],
            khat_dot[2],
        ],
        dtype=float,
    )


def rk4_step(
    rhs: Callable[[float, np.ndarray, Parameters, np.ndarray], np.ndarray],
    t: float,
    z: np.ndarray,
    dt: float,
    p: Parameters,
    P: np.ndarray,
) -> np.ndarray:
    k1 = rhs(t, z, p, P)
    k2 = rhs(t + 0.5 * dt, z + 0.5 * dt * k1, p, P)
    k3 = rhs(t + 0.5 * dt, z + 0.5 * dt * k2, p, P)
    k4 = rhs(t + dt, z + dt * k3, p, P)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(p: Parameters) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    P = lyapunov_matrix_for_q_identity(p.am1, p.am2)
    n_steps = int(round(p.final_time / p.dt)) + 1
    time = np.linspace(0.0, p.final_time, n_steps)

    # z = [x1, x2, xm1, xm2, khat1, khat2, khatr]
    z = np.zeros(7, dtype=float)
    history = np.zeros((n_steps, z.size), dtype=float)

    for i, t in enumerate(time):
        history[i] = z
        if i + 1 < n_steps:
            z = rk4_step(derivative, t, z, p.dt, p, P)
            if not np.all(np.isfinite(z)):
                raise FloatingPointError("Simulation became non-finite; reduce dt or adaptation gains")

    return time, history, P


def save_results(time: np.ndarray, history: np.ndarray, p: Parameters, P: np.ndarray) -> None:
    output_dir = Path(__file__).resolve().parent
    x1, x2, xm1, xm2 = history[:, 0], history[:, 1], history[:, 2], history[:, 3]
    gains = history[:, 4:7]
    error = history[:, 0:2] - history[:, 2:4]

    csv_data = np.column_stack((time, history, error))
    header = "t,x1,x2,xm1,xm2,khat1,khat2,khatr,e1,e2"
    np.savetxt(output_dir / "Chapter6_Lesson4_results.csv", csv_data, delimiter=",", header=header, comments="")

    fig1, ax1 = plt.subplots()
    ax1.plot(time, x1, label="plant x1")
    ax1.plot(time, xm1, "--", label="model xm1")
    ax1.set_xlabel("Time (s)")
    ax1.set_ylabel("Position-like state")
    ax1.grid(True)
    ax1.legend()
    fig1.tight_layout()
    fig1.savefig(output_dir / "Chapter6_Lesson4_tracking.png", dpi=180)

    fig2, ax2 = plt.subplots()
    ax2.plot(time, gains[:, 0], label="khat1")
    ax2.plot(time, gains[:, 1], label="khat2")
    ax2.plot(time, gains[:, 2], label="khatr")
    ax2.set_xlabel("Time (s)")
    ax2.set_ylabel("Adaptive gains")
    ax2.grid(True)
    ax2.legend()
    fig2.tight_layout()
    fig2.savefig(output_dir / "Chapter6_Lesson4_parameters.png", dpi=180)

    ideal = np.array(
        [
            (p.a1 - p.am1) / p.b,
            (p.a2 - p.am2) / p.b,
            p.bm / p.b,
        ]
    )
    print("P =\n", P)
    print("eigenvalues(P) =", np.linalg.eigvals(P))
    print("final tracking error =", error[-1])
    print("final estimated gains =", gains[-1])
    print("ideal matching gains (not generally guaranteed to be identified) =", ideal)


if __name__ == "__main__":
    parameters = Parameters()
    t_values, state_history, P_matrix = simulate(parameters)
    save_results(t_values, state_history, parameters, P_matrix)

Chapter6_Lesson4.cpp

// Chapter 6, Lesson 4: Parameter update laws for second-order MRAC.
// Compile: g++ -std=c++17 -O2 Chapter6_Lesson4.cpp -o Chapter6_Lesson4

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

struct Parameters {
    double a1{1.0};
    double a2{0.8};
    double b{1.5};
    double am1{2.25};
    double am2{2.10};
    double bm{2.25};
    double gamma1{4.0};
    double gamma2{4.0};
    double gammar{2.0};
    double dt{0.002};
    double finalTime{20.0};
};

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

struct LyapunovMatrix {
    double p11{};
    double p12{};
    double p22{};
};

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

LyapunovMatrix lyapunovMatrixForQIdentity(double am1, double am2) {
    if (am1 <= 0.0 || am2 <= 0.0) {
        throw std::invalid_argument("am1 and am2 must be positive");
    }
    LyapunovMatrix P;
    P.p12 = 1.0 / (2.0 * am1);
    P.p22 = (1.0 + am1) / (2.0 * am1 * am2);
    P.p11 = am2 / (2.0 * am1) + (1.0 + am1) / (2.0 * am2);
    return P;
}

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

State derivative(double t, const State& z, const Parameters& p, const LyapunovMatrix& P) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double khat1 = z[4];
    const double khat2 = z[5];
    const double khatr = z[6];
    const double r = reference(t);

    const double u = khat1 * x1 + khat2 * x2 + khatr * r;
    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double s = P.p12 * e1 + P.p22 * e2;
    const double signB = p.b > 0.0 ? 1.0 : -1.0;

    return State{
        x2,
        -p.a1 * x1 - p.a2 * x2 + p.b * u,
        xm2,
        -p.am1 * xm1 - p.am2 * xm2 + p.bm * r,
        -p.gamma1 * x1 * s * signB,
        -p.gamma2 * x2 * s * signB,
        -p.gammar * r * s * signB,
    };
}

State rk4Step(double t, const State& z, const Parameters& p, const LyapunovMatrix& P) {
    const State k1 = derivative(t, z, p, P);
    const State k2 = derivative(t + 0.5 * p.dt, addScaled(z, k1, 0.5 * p.dt), p, P);
    const State k3 = derivative(t + 0.5 * p.dt, addScaled(z, k2, 0.5 * p.dt), p, P);
    const State k4 = derivative(t + p.dt, addScaled(z, k3, p.dt), p, P);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (p.dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Simulation became non-finite");
        }
    }
    return next;
}

int main() {
    try {
        const Parameters p;
        const LyapunovMatrix P = lyapunovMatrixForQIdentity(p.am1, p.am2);
        const int steps = static_cast<int>(std::llround(p.finalTime / p.dt));
        State z{};

        std::ofstream out("Chapter6_Lesson4_results_cpp.csv");
        if (!out) {
            throw std::runtime_error("Cannot open CSV output file");
        }
        out << "t,x1,x2,xm1,xm2,khat1,khat2,khatr,e1,e2\n";
        out << std::setprecision(12);

        for (int i = 0; i <= steps; ++i) {
            const double t = i * p.dt;
            out << t;
            for (double value : z) {
                out << ',' << value;
            }
            out << ',' << z[0] - z[2] << ',' << z[1] - z[3] << '\n';
            if (i < steps) {
                z = rk4Step(t, z, p, P);
            }
        }

        const double ideal1 = (p.a1 - p.am1) / p.b;
        const double ideal2 = (p.a2 - p.am2) / p.b;
        const double idealR = p.bm / p.b;

        std::cout << "P = [[" << P.p11 << ", " << P.p12 << "], ["
                  << P.p12 << ", " << P.p22 << "]]\n";
        std::cout << "Final tracking error = [" << z[0] - z[2] << ", "
                  << z[1] - z[3] << "]\n";
        std::cout << "Final estimated gains = [" << z[4] << ", " << z[5]
                  << ", " << z[6] << "]\n";
        std::cout << "Ideal matching gains = [" << ideal1 << ", " << ideal2
                  << ", " << idealR << "]\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter6_Lesson4.java

// Chapter 6, Lesson 4: Parameter update laws for second-order MRAC.
// Compile: javac Chapter6_Lesson4.java
// Run:     java Chapter6_Lesson4

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_Lesson4 {
    private static final class Parameters {
        final double a1 = 1.0;
        final double a2 = 0.8;
        final double b = 1.5;
        final double am1 = 2.25;
        final double am2 = 2.10;
        final double bm = 2.25;
        final double gamma1 = 4.0;
        final double gamma2 = 4.0;
        final double gammar = 2.0;
        final double dt = 0.002;
        final double finalTime = 20.0;
    }

    private record LyapunovMatrix(double p11, double p12, double p22) {}

    private Chapter6_Lesson4() {}

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

    private static LyapunovMatrix lyapunovMatrixForQIdentity(double am1, double am2) {
        if (am1 <= 0.0 || am2 <= 0.0) {
            throw new IllegalArgumentException("am1 and am2 must be positive");
        }
        double p12 = 1.0 / (2.0 * am1);
        double p22 = (1.0 + am1) / (2.0 * am1 * am2);
        double p11 = am2 / (2.0 * am1) + (1.0 + am1) / (2.0 * am2);
        return new LyapunovMatrix(p11, p12, p22);
    }

    private static double[] derivative(
            double t, double[] z, Parameters p, LyapunovMatrix P) {
        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double khat1 = z[4];
        double khat2 = z[5];
        double khatr = z[6];
        double r = reference(t);

        double u = khat1 * x1 + khat2 * x2 + khatr * r;
        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double s = P.p12() * e1 + P.p22() * e2;
        double signB = p.b > 0.0 ? 1.0 : -1.0;

        return new double[] {
            x2,
            -p.a1 * x1 - p.a2 * x2 + p.b * u,
            xm2,
            -p.am1 * xm1 - p.am2 * xm2 + p.bm * r,
            -p.gamma1 * x1 * s * signB,
            -p.gamma2 * x2 * s * signB,
            -p.gammar * r * s * signB
        };
    }

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

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

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; ++i) {
            next[i] = z[i] + (p.dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(next[i])) {
                throw new ArithmeticException("Simulation became non-finite");
            }
        }
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Parameters p = new Parameters();
        LyapunovMatrix P = lyapunovMatrixForQIdentity(p.am1, p.am2);
        int steps = (int) Math.round(p.finalTime / p.dt);
        double[] z = new double[7];
        Path output = Path.of("Chapter6_Lesson4_results_java.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(output)) {
            writer.write("t,x1,x2,xm1,xm2,khat1,khat2,khatr,e1,e2\n");
            for (int i = 0; i <= steps; ++i) {
                double t = i * p.dt;
                writer.write(String.format(
                        Locale.US,
                        "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                        t, z[0], z[1], z[2], z[3], z[4], z[5], z[6],
                        z[0] - z[2], z[1] - z[3]));
                if (i < steps) {
                    z = rk4Step(t, z, p, P);
                }
            }
        } catch (IOException ex) {
            System.err.println("Cannot write output: " + ex.getMessage());
            System.exit(1);
        }

        double ideal1 = (p.a1 - p.am1) / p.b;
        double ideal2 = (p.a2 - p.am2) / p.b;
        double idealR = p.bm / p.b;

        System.out.printf("P = [[%.8f, %.8f], [%.8f, %.8f]]%n",
                P.p11(), P.p12(), P.p12(), P.p22());
        System.out.printf("Final tracking error = [%.8g, %.8g]%n", z[0] - z[2], z[1] - z[3]);
        System.out.printf("Final estimated gains = [%.8g, %.8g, %.8g]%n", z[4], z[5], z[6]);
        System.out.printf("Ideal matching gains = [%.8g, %.8g, %.8g]%n", ideal1, ideal2, idealR);
    }
}

Chapter6_Lesson4.m

%% Chapter6_Lesson4.m
% Parameter update laws for a continuous-time second-order direct MRAC.
% The script uses ode45 to integrate the plant, reference model, and gains.

clear; clc; close all;

p.a1 = 1.0;
p.a2 = 0.8;
p.b = 1.5;
p.am1 = 2.25;
p.am2 = 2.10;
p.bm = 2.25;
p.Gamma = diag([4.0, 4.0, 2.0]);
p.sign_b = sign(p.b);

% For A_m = [0 1; -am1 -am2] and Q = I, solve A_m'P + P*A_m = -I.
p.P = [p.am2/(2*p.am1) + (1+p.am1)/(2*p.am2), 1/(2*p.am1); ...
       1/(2*p.am1), (1+p.am1)/(2*p.am1*p.am2)];

z0 = zeros(7,1); % [x1 x2 xm1 xm2 khat1 khat2 khatr]'
tspan = [0 20];
options = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,z] = ode45(@(t,z) mrac_rhs(t,z,p), tspan, z0, options);

x = z(:,1:2);
xm = z(:,3:4);
khat = z(:,5:7);
e = x - xm;

figure;
plot(t,x(:,1),'LineWidth',1.4); hold on;
plot(t,xm(:,1),'--','LineWidth',1.4);
grid on; xlabel('Time (s)'); ylabel('Position-like state');
legend('Plant x_1','Model x_{m1}','Location','best');
title('Second-Order MRAC Tracking');

figure;
plot(t,khat,'LineWidth',1.3);
grid on; xlabel('Time (s)'); ylabel('Adaptive gains');
legend('\hat{k}_1','\hat{k}_2','\hat{k}_r','Location','best');
title('Parameter Update Laws');

ideal = [(p.a1-p.am1)/p.b, (p.a2-p.am2)/p.b, p.bm/p.b];
fprintf('P =\n'); disp(p.P);
fprintf('eig(P) =\n'); disp(eig(p.P));
fprintf('Final tracking error: [%g, %g]\n', e(end,1), e(end,2));
fprintf('Final estimated gains: [%g, %g, %g]\n', khat(end,1), khat(end,2), khat(end,3));
fprintf('Ideal matching gains (not necessarily identified): [%g, %g, %g]\n', ideal);

results = table(t,x(:,1),x(:,2),xm(:,1),xm(:,2),khat(:,1),khat(:,2),khat(:,3),e(:,1),e(:,2), ...
    'VariableNames',{'t','x1','x2','xm1','xm2','khat1','khat2','khatr','e1','e2'});
writetable(results,'Chapter6_Lesson4_results_matlab.csv');

function dz = mrac_rhs(t,z,p)
    x1 = z(1); x2 = z(2);
    xm1 = z(3); xm2 = z(4);
    khat = z(5:7);

    r = double(t >= 0.5);
    phi = [x1; x2; r];
    u = khat.' * phi;

    e = [x1-xm1; x2-xm2];
    B = [0;1];
    s = e.' * p.P * B;
    khat_dot = -p.Gamma * phi * s * p.sign_b;

    dz = [x2; ...
          -p.a1*x1 - p.a2*x2 + p.b*u; ...
          xm2; ...
          -p.am1*xm1 - p.am2*xm2 + p.bm*r; ...
          khat_dot];
end

The following script programmatically creates and saves Chapter6_Lesson4_MRAC.slx. It requires Simulink.

Chapter6_Lesson4_Simulink.m

%% Chapter6_Lesson4_Simulink.m
% Programmatically constructs a Simulink diagram for the second-order MRAC.
% Requires Simulink. The model uses one MATLAB Function block containing the
% plant, model, and continuous parameter-update differential equations.

clear; clc;
model = 'Chapter6_Lesson4_MRAC';

if bdIsLoaded(model)
    close_system(model,0);
end
new_system(model);
open_system(model);

add_block('simulink/Sources/Step',[model '/Reference'], ...
    'Time','0.5','Before','0','After','1', ...
    'Position',[40 95 70 125]);
add_block('simulink/User-Defined Functions/MATLAB Function',[model '/MRAC Dynamics'], ...
    'Position',[145 55 315 165]);
add_block('simulink/Continuous/Integrator',[model '/Integrator'], ...
    'InitialCondition','zeros(7,1)', ...
    'Position',[360 75 390 145]);
add_block('simulink/Signal Routing/Demux',[model '/Demux'], ...
    'Outputs','7','Position',[450 45 455 180]);
add_block('simulink/Sinks/Scope',[model '/Tracking Scope'], ...
    'NumInputPorts','2','Position',[555 45 585 95]);
add_block('simulink/Sinks/Scope',[model '/Gain Scope'], ...
    'NumInputPorts','3','Position',[555 115 585 180]);

add_line(model,'Reference/1','MRAC Dynamics/1');
add_line(model,'Integrator/1','MRAC Dynamics/2');
add_line(model,'MRAC Dynamics/1','Integrator/1');
add_line(model,'Integrator/1','Demux/1');
add_line(model,'Demux/1','Tracking Scope/1');
add_line(model,'Demux/3','Tracking Scope/2');
add_line(model,'Demux/5','Gain Scope/1');
add_line(model,'Demux/6','Gain Scope/2');
add_line(model,'Demux/7','Gain Scope/3');

blockPath = [model '/MRAC Dynamics'];
code = [ ...
"function dz = fcn(r,z)" newline ...
"%#codegen" newline ...
"a1=1.0; a2=0.8; b=1.5;" newline ...
"am1=2.25; am2=2.10; bm=2.25;" newline ...
"Gamma=diag([4.0 4.0 2.0]);" newline ...
"P=[am2/(2*am1)+(1+am1)/(2*am2), 1/(2*am1);" newline ...
"   1/(2*am1), (1+am1)/(2*am1*am2)];" newline ...
"x1=z(1); x2=z(2); xm1=z(3); xm2=z(4);" newline ...
"khat=z(5:7); phi=[x1;x2;r];" newline ...
"u=khat.'*phi; e=[x1-xm1;x2-xm2];" newline ...
"s=e.'*P*[0;1];" newline ...
"khatDot=-Gamma*phi*s*sign(b);" newline ...
"dz=[x2; -a1*x1-a2*x2+b*u; xm2; -am1*xm1-am2*xm2+bm*r; khatDot];" newline ...
"end" ...
];
rt = sfroot;
chart = rt.find('-isa','Stateflow.EMChart','Path',blockPath);
if isempty(chart)
    error('Could not access the MATLAB Function block through the Stateflow API.');
end
chart.Script = char(code);

set_param(model,'StopTime','20','Solver','ode45','MaxStep','0.002');
save_system(model,[model '.slx']);
open_system(model);

Chapter6_Lesson4.nb

(* Chapter6_Lesson4.nb *)
ClearAll["Global`*"];

a1 = 1.0; a2 = 0.8; b = 1.5;
am1 = 2.25; am2 = 2.10; bm = 2.25;
gamma1 = 4.0; gamma2 = 4.0; gammar = 2.0;

p12 = 1/(2 am1);
p22 = (1 + am1)/(2 am1 am2);
p11 = am2/(2 am1) + (1 + am1)/(2 am2);
P = {{p11, p12}, {p12, p22}};
Gamma = DiagonalMatrix[{gamma1, gamma2, gammar}];

r[t_] := Piecewise[{{0, t < 0.5}}, 1];

sol = NDSolveValue[
  {
    x1'[t] == x2[t],
    x2'[t] == -a1 x1[t] - a2 x2[t] +
      b (k1[t] x1[t] + k2[t] x2[t] + kr[t] r[t]),
    xm1'[t] == xm2[t],
    xm2'[t] == -am1 xm1[t] - am2 xm2[t] + bm r[t],
    k1'[t] == -gamma1 x1[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    k2'[t] == -gamma2 x2[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    kr'[t] == -gammar r[t]
      (p12 (x1[t] - xm1[t]) + p22 (x2[t] - xm2[t])),
    x1[0] == 0, x2[0] == 0, xm1[0] == 0, xm2[0] == 0,
    k1[0] == 0, k2[0] == 0, kr[0] == 0
  },
  {x1, x2, xm1, xm2, k1, k2, kr},
  {t, 0, 20}, MaxStepSize -> 0.002
];

Plot[Evaluate[{sol[[1]][t], sol[[3]][t]}], {t, 0, 20},
 PlotLegends -> {"Plant x1", "Model xm1"}, PlotRange -> All]

Plot[Evaluate[Table[sol[[i]][t], {i, 5, 7}]], {t, 0, 20},
 PlotLegends -> {"khat1", "khat2", "khatr"}, PlotRange -> All]

Eigenvalues[P]

The downloadable notebook contains the same equations as executable notebook cells.

A correct implementation should verify all of the following:

1. Reference-model stability: the roots of \(s^2+a_{m2}s+a_{m1}\) have negative real parts.
2. Lyapunov residual: \(A_m^{T}P+PA_m+Q\) is numerically zero.
3. Positive definiteness: both eigenvalues of \(P\) are positive.
4. Finite integration: all plant, model, and gain states remain finite over the simulated interval.
5. Error trend: the tracking error approaches zero after the transient.

For the numerical values in Section 10, the eigenvalues of \(P\) are approximately \(1.29253\) and \(0.291858\), confirming positive definiteness.

Problem 1 (Coefficient Matching): For the plant and reference model in Section 2, derive the ideal gains \(k_1^{*}\), \(k_2^{*}\), and \(k_r^{*}\).

Solution: With the ideal controller, the plant becomes

\[ \dot{x}_2=(-a_1+b k_1^{*})x_1 +(-a_2+b k_2^{*})x_2+b k_r^{*}r. \]

Matching the three coefficients with the model gives

\[ k_1^{*}=\frac{a_1-a_{m1}}{b}, \qquad k_2^{*}=\frac{a_2-a_{m2}}{b}, \qquad k_r^{*}=\frac{b_m}{b}. \]

Problem 2 (Solving the Lyapunov Equation): Let \(a_{m1}=4\), \(a_{m2}=4\), and \(Q=I_2\). Compute \(P\) using the formulas derived in Section 4.

Solution:

\[ p_{12}=\frac{1}{8}, \qquad p_{22}=\frac{5}{32}, \qquad p_{11}=\frac{4}{8}+\frac{5}{8}=\frac{9}{8}. \]

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

Its leading principal minors are \(9/8>0\) and \(\det(P)=43/256>0\), so \(P\succ0\).

Problem 3 (Cross-Term Cancellation): Starting from

\[ \dot{V}=-\mathbf{e}^{T}Q\mathbf{e} +2b s\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +2|b|\widetilde{\boldsymbol{\theta}}^{T}\Gamma^{-1} \dot{\hat{\boldsymbol{\theta}}}, \]

show that the proposed update law gives a negative-semidefinite derivative.

Solution: Substitute

\[ \dot{\hat{\boldsymbol{\theta}}} =-\Gamma\boldsymbol{\phi}s\operatorname{sgn}(b). \]

The last term becomes

\[ -2|b|\operatorname{sgn}(b) \widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}s =-2b\widetilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}s, \]

which exactly cancels the positive cross term. Hence \(\dot V=-\mathbf{e}^{T}Q\mathbf{e}\le0\).

Problem 4 (Component Update at One Instant): Suppose \(b>0\), \(\Gamma=\operatorname{diag}(3,2,1)\), \(\boldsymbol{\phi}=[2,-1,0.5]^{T}\), and \(s=-0.4\). Compute the gain derivatives.

Solution:

\[ \dot{\hat{\boldsymbol{\theta}}} =-\begin{bmatrix}3&0&0\\0&2&0\\0&0&1\end{bmatrix} \begin{bmatrix}2\\-1\\0.5\end{bmatrix}(-0.4) =\begin{bmatrix}2.4\\-0.8\\0.2\end{bmatrix}. \]

Problem 5 (Known Negative Input Gain): Assume the plant has \(b<0\), but its sign is known. Write the three component update laws.

Solution: Since \(\operatorname{sgn}(b)=-1\),

\[ \dot{\hat{k}}_1=+\gamma_1x_1s, \qquad \dot{\hat{k}}_2=+\gamma_2x_2s, \qquad \dot{\hat{k}}_r=+\gamma_r r s. \]

Reversing the update direction is essential because the error dynamics contain the factor \(b\). Using the positive-gain law without this sign correction would add rather than cancel the cross term.

Problem 6 (Why Parameter Equality Is Not Proved): Explain why \(\dot V=-\mathbf{e}^{T}Q\mathbf{e}\) does not imply \(\widetilde{\boldsymbol{\theta}}\to\mathbf{0}\).

Solution:

The derivative contains no negative-definite term in the parameter error. It forces the tracking error energy to decrease, while only bounding the parameter error through the value of \(V\). Once \(s\) or the relevant regressor products vanish, parameter motion may stop at a bounded nonzero error. Additional signal conditions are required to identify a unique ideal vector.

For a matched second-order SISO plant, the tracking error satisfies a stable reference-model dynamics plus a linearly parameterized uncertainty term. Solving the reference-model Lyapunov equation produces the scalar signal \(s=e^{T}PB\). The update law \(\dot{\hat{\boldsymbol{\theta}}} =-\Gamma\boldsymbol{\phi}s\operatorname{sgn}(b)\) cancels the state-parameter cross term exactly, yielding \(\dot V=-e^{T}Qe\le0\). The proof establishes bounded adaptive gains and asymptotic tracking under the stated assumptions, but it does not by itself guarantee convergence of the estimates to the ideal gains.