After completing this lesson, students should be able to:

• formulate a general SISO plant and reference model in controllable canonical form;
• derive the model-matching equations and identify the ideal controller parameters;
• obtain a linearly parameterized tracking-error system;
• construct a Lyapunov function containing both tracking and parameter errors;
• derive a vector adaptive law that cancels the indefinite cross term;
• prove global boundedness and asymptotic tracking using Barbalat's lemma;
• implement and numerically verify the design for a third-order plant.

The lesson assumes that the complete plant state is measurable, the plant order is known, the plant is represented in the canonical structure introduced in Lesson 1, and the sign of the scalar high-frequency gain is known. The unknown-sign case is deliberately postponed to later lessons.

Consider a SISO plant with state \(\mathbf{x}\in\mathbb{R}^{n}\), scalar input \(u\in\mathbb{R}\), and measurable output \(y=x_1\):

\[ \dot{\mathbf{x} }(t) = \mathbf{A}_p\mathbf{x}(t) + \mathbf{B}b_0u(t), \qquad y(t)=\mathbf{c}^{T}\mathbf{x}(t), \]

where

\[ \mathbf{B} = \begin{bmatrix} 0 & 0 & \cdots & 0 & 1 \end{bmatrix}^{T}, \qquad \mathbf{c} = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}^{T}, \]

and the unknown companion matrix is

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

The scalar \(b_0\neq 0\) is the plant input gain. Its magnitude may be unknown; for the present design, its sign is known. A desired stable reference model is selected as

\[ \dot{\mathbf{x} }_m(t) = \mathbf{A}_m\mathbf{x}_m(t) + \mathbf{B}b_m r(t), \qquad y_m(t)=\mathbf{c}^{T}\mathbf{x}_m(t), \]

where \(\mathbf{A}_m\) is Hurwitz, \(b_m\) is selected to obtain the desired command gain, and \(r(t)\) is bounded and piecewise continuous.

\[ \mathbf{A}_m = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ -a_{m0} & -a_{m1} & -a_{m2} & \cdots & -a_{m,n-1} \end{bmatrix}. \]

Use the direct adaptive state-feedback structure

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

Introduce the combined parameter vector and regressor

\[ \widehat{\boldsymbol{\theta} } = \begin{bmatrix} \widehat{\mathbf{k} }_x \\ \widehat{k}_r \end{bmatrix} \in\mathbb{R}^{n+1}, \qquad \boldsymbol{\phi} = \begin{bmatrix} \mathbf{x} \\ r \end{bmatrix} \in\mathbb{R}^{n+1}, \qquad u=\widehat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]

Suppose constant ideal parameters \(\mathbf{k}_x^{*}\) and \(k_r^{*}\) exist such that the plant under ideal feedback exactly reproduces the reference-model dynamics. The matching equations are

\[ \mathbf{A}_p+\mathbf{B}b_0\mathbf{k}_x^{*T} = \mathbf{A}_m, \qquad b_0k_r^{*}=b_m. \]

Because only the final row differs between the two companion matrices, the first equation can be solved componentwise:

\[ k_{x,i}^{*} = \frac{a_{p,i-1}-a_{m,i-1} }{b_0}, \qquad i=1,\ldots,n, \qquad k_r^{*}=\frac{b_m}{b_0}. \]

Define \(\boldsymbol{\theta}^{*} = [\,\mathbf{k}_x^{*T}\;k_r^{*}\,]^{T}\) and the parameter error

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

Since the ideal plant parameters are constant, \(\dot{\widetilde{\boldsymbol{\theta} } } =\dot{\widehat{\boldsymbol{\theta} } }\).

Define the full-state tracking error

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

Substitution of the adaptive controller into the plant yields

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

Add and subtract the ideal parameters and apply the matching equations:

\[ \begin{aligned} \dot{\mathbf{x} } &= \left( \mathbf{A}_p+\mathbf{B}b_0\mathbf{k}_x^{*T} \right)\mathbf{x} + \mathbf{B}b_0k_r^{*}r + \mathbf{B}b_0 \left( \widetilde{\mathbf{k} }_x^{T}\mathbf{x} + \widetilde{k}_r r \right) \\ &= \mathbf{A}_m\mathbf{x} + \mathbf{B}b_mr + \mathbf{B}b_0 \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \end{aligned} \]

Subtracting the reference-model equation gives the central error model:

\[ \boxed{ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} + \mathbf{B}b_0 \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} }. \]

The unknown parameters now appear linearly through the scalar inner product \(\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\). This structure makes a Lyapunov cancellation design possible.

flowchart TD
  R["Bounded command r"] --> RM["Stable reference model"]
  R --> PHI["Regressor phi = [plant state; command]"]
  X["Measured plant state x"] --> PHI
  PHI --> C["Adaptive control u = theta_hat dot phi"]
  C --> P["Unknown canonical-form plant"]
  P --> X
  X --> E["State error e = x - xm"]
  RM --> XM["Reference state xm"]
  XM --> E
  E --> S["Scalar error s = e transpose P B"]
  S --> A["Parameter update law"]
  PHI --> A
  A --> C
        

Since \(\mathbf{A}_m\) is Hurwitz, for every symmetric positive-definite matrix \(\mathbf{Q}\), 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}\). Let \(\boldsymbol{\Gamma} =\boldsymbol{\Gamma}^{T}>0\) be the adaptation-gain matrix. Choose

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

This candidate is positive definite and radially unbounded in the augmented error vector \([\,\mathbf{e}^{T}\; \widetilde{\boldsymbol{\theta} }^{T}\,]^{T}\). Differentiating along the error dynamics gives

\[ \begin{aligned} \dot V &= \dot{\mathbf{e} }^{T}\mathbf{P}\mathbf{e} + \mathbf{e}^{T}\mathbf{P}\dot{\mathbf{e} } + 2|b_0| \widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\widehat{\boldsymbol{\theta} } } \\ &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} + 2b_0 \left(\mathbf{e}^{T}\mathbf{P}\mathbf{B}\right) \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} + 2|b_0| \widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\widehat{\boldsymbol{\theta} } }. \end{aligned} \]

The first term is negative definite in the tracking error, whereas the second term has indefinite sign. The adaptive law is selected to cancel this cross term exactly.

Define the scalar filtered tracking error

\[ s(t)=\mathbf{e}^{T}(t)\mathbf{P}\mathbf{B}. \]

Select the parameter update law

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } = -\operatorname{sgn}(b_0) \boldsymbol{\Gamma} \boldsymbol{\phi} \mathbf{e}^{T}\mathbf{P}\mathbf{B} }. \]

Equivalently, the state-feedback and feedforward parts are

\[ \dot{\widehat{\mathbf{k} } }_x = -\operatorname{sgn}(b_0) \boldsymbol{\Gamma}_x \mathbf{x} \mathbf{e}^{T}\mathbf{P}\mathbf{B}, \qquad \dot{\widehat{k} }_r = -\operatorname{sgn}(b_0) \gamma_r r \mathbf{e}^{T}\mathbf{P}\mathbf{B}. \]

Substitution into \(\dot V\) gives

\[ \begin{aligned} \dot V &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} + 2b_0s \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} - 2|b_0|\operatorname{sgn}(b_0)s \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \\ &= -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \leq 0, \end{aligned} \]

because \(|b_0|\operatorname{sgn}(b_0)=b_0\). Thus the controller and update law are not guessed independently: the update law is the exact algebraic choice required by the composite Lyapunov derivative.

Theorem. Assume that:

1. the plant has the stated controllable canonical form and the full state is measurable;
2. the matching equations admit constant ideal parameters;
3. the sign of \(b_0\) is known;
4. \(\mathbf{A}_m\) is Hurwitz and \(\mathbf{Q}=\mathbf{Q}^{T}>0\);
5. the reference input \(r(t)\) is bounded and piecewise continuous.

Then, under the adaptive controller and update law derived above, all closed-loop signals remain bounded and \(\lim_{t\to\infty}\mathbf{e}(t)=\mathbf{0}\).

Proof, Step 1 — Lyapunov monotonicity.

\[ \dot V = -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} \leq -\lambda_{\min}(\mathbf{Q})\|\mathbf{e}\|^2 \leq 0. \]

Hence \(V(t)\leq V(0)\) for all \(t\geq 0\).

Proof, Step 2 — Bounded error and parameter error.

\[ \lambda_{\min}(\mathbf{P})\|\mathbf{e}\|^2 + |b_0|\lambda_{\min} \left(\boldsymbol{\Gamma}^{-1}\right) \|\widetilde{\boldsymbol{\theta} }\|^2 \leq V(t)\leq V(0). \]

Therefore \(\mathbf{e}\in\mathcal{L}_{\infty}\) and \(\widetilde{\boldsymbol{\theta} } \in\mathcal{L}_{\infty}\). Since \(\boldsymbol{\theta}^{*}\) is constant, \(\widehat{\boldsymbol{\theta} }\) is bounded.

Proof, Step 3 — Square integrability of the tracking error.

Integration of the Lyapunov derivative gives

\[ V(T)-V(0) = -\int_{0}^{T}\mathbf{e}^{T}(\tau) \mathbf{Q}\mathbf{e}(\tau)\,d\tau. \]

Since \(V(T)\geq 0\),

\[ \lambda_{\min}(\mathbf{Q}) \int_{0}^{\infty}\|\mathbf{e}(\tau)\|^2d\tau \leq V(0), \]

so \(\mathbf{e}\in\mathcal{L}_2\).

Proof, Step 4 — Boundedness of all remaining signals.

The stable reference model driven by bounded \(r(t)\) gives \(\mathbf{x}_m\in\mathcal{L}_{\infty}\). Since \(\mathbf{x}=\mathbf{e}+\mathbf{x}_m\), the plant state is bounded. Therefore the regressor \(\boldsymbol{\phi}=[\,\mathbf{x}^{T}\;r\,]^{T}\), adaptive control input \(u=\widehat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\), parameter derivative, and tracking-error derivative are bounded.

Proof, Step 5 — Asymptotic tracking.

The boundedness of \(\dot{\mathbf{e} }\) implies uniform continuity of \(\mathbf{e}\). Together with \(\mathbf{e}\in\mathcal{L}_2\), Barbalat's lemma yields

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

Consequently, \(\lim_{t\to\infty}(y(t)-y_m(t))=0\). The proof establishes tracking, but it does not by itself establish convergence of every adaptive parameter to its ideal value. Formal conditions for parameter convergence are developed later in the course.

flowchart TD
  A["Choose Hurwitz reference matrix Am"] --> B["Solve Am transpose P + P Am = -Q"]
  B --> C["Form V = state-error energy + parameter-error energy"]
  C --> D["Differentiate V along error dynamics"]
  D --> E["Select update law to cancel cross term"]
  E --> F["Obtain Vdot = -e transpose Q e"]
  F --> G["e and parameter error are bounded"]
  G --> H["Stable model plus bounded command gives bounded plant signals"]
  H --> I["e is square integrable and e-dot is bounded"]
  I --> J["Barbalat argument"]
  J --> K["State tracking error tends to zero"]
        

The matrix \(\mathbf{Q}\) does not directly set the reference-model poles; those are selected through \(\mathbf{A}_m\). Instead, \(\mathbf{Q}\) determines \(\mathbf{P}\) and therefore shapes the scalar error \(s=\mathbf{e}^{T}\mathbf{P}\mathbf{B}\) used by the parameter update law.

For a diagonal adaptation matrix

\[ \boldsymbol{\Gamma} = \operatorname{diag} \left( \gamma_{x1},\ldots,\gamma_{xn},\gamma_r \right), \qquad \gamma_{xi}>0,\quad\gamma_r>0, \]

increasing a diagonal entry increases the instantaneous rate of change of the corresponding adaptive parameter. The ideal continuous-time proof remains valid for every positive-definite \(\boldsymbol{\Gamma}\); however, very large gains can produce sharp transients and make a numerical or sampled-data implementation difficult. These practical effects do not contradict the continuous-time theorem.

A useful a priori error-energy bound follows directly from the proof:

\[ \int_{0}^{\infty}\|\mathbf{e}(t)\|^2dt \leq \frac{V(0)}{\lambda_{\min}(\mathbf{Q})}. \]

The bound shows explicitly that the initial state error and initial parameter error both contribute to the total tracking-error energy.

Consider the uncertain third-order plant

\[ \dot{\mathbf{x} } = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & -2 & -1.2 \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} 1.4u. \]

Choose the model polynomial \((s+1)(s+2)(s+3)=s^3+6s^2+11s+6\), giving

\[ \mathbf{A}_m = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -6 & -11 & -6 \end{bmatrix}, \qquad b_m=6. \]

The ideal parameters, used only for analysis and plotting, are

\[ \mathbf{k}_x^{*} = \frac{1}{1.4} \begin{bmatrix} -5 \\ -9 \\ -4.8 \end{bmatrix} = \begin{bmatrix} -3.5714 \\ -6.4286 \\ -3.4286 \end{bmatrix}, \qquad k_r^{*}=\frac{6}{1.4}=4.2857. \]

Set \(\mathbf{Q}=\mathbf{I}_3\). Solving the Lyapunov equation gives

\[ \mathbf{P} = \begin{bmatrix} 1.8166667 & 1.15 & 0.0833333 \\ 1.15 & 2.0083333 & 0.15 \\ 0.0833333 & 0.15 & 0.1083333 \end{bmatrix}. \]

Thus \(\mathbf{P}\mathbf{B} =[\,0.0833333\;0.15\;0.1083333\,]^{T}\) and

\[ s = 0.0833333e_1+0.15e_2+0.1083333e_3. \]

The implementations use \(\boldsymbol{\Gamma} =\operatorname{diag}(8,8,8,4)\), \(\widehat{\boldsymbol{\theta} }(0)=\mathbf{0}\), and \(\mathbf{x}(0)=[\,0.5\;-0.2\;0.1\,]^{T}\). The command is

\[ r(t) = \begin{cases} 1, & 0\leq t < 8, \\ -0.5, & 8\leq t < 16, \\ 0.8, & t\geq 16. \end{cases} \]

The Python version uses NumPy for matrix operations, SciPy for adaptive ODE integration, and Matplotlib for visualization. Install the dependencies with python -m pip install numpy scipy matplotlib.

Chapter7_Lesson3.py

"""
Chapter7_Lesson3.py
Lyapunov-based MRAC for a third-order SISO plant in controllable canonical form.

Dependencies:
    numpy
    scipy
    matplotlib

The controller is
    u = theta_hat^T phi,  phi = [x1, x2, x3, r]^T

and the adaptive law is
    theta_hat_dot = -sign(b0) Gamma phi e^T P B.
"""

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

N = 3
B = np.array([0.0, 0.0, 1.0])
B0 = 1.4
BM = 6.0

AP = np.array(
    [
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0],
        [-1.0, -2.0, -1.2],
    ]
)

AM = np.array(
    [
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0],
        [-6.0, -11.0, -6.0],
    ]
)

# P solves AM.T @ P + P @ AM = -I.
P = np.array(
    [
        [1.8166666666666667, 1.15, 0.0833333333333333],
        [1.15, 2.0083333333333333, 0.15],
        [0.0833333333333333, 0.15, 0.1083333333333333],
    ]
)

GAMMA = np.diag([8.0, 8.0, 8.0, 4.0])


def reference(t: float) -> float:
    """Bounded piecewise-constant command."""
    if t < 8.0:
        return 1.0
    if t < 16.0:
        return -0.5
    return 0.8


def dynamics(t: float, z: np.ndarray) -> np.ndarray:
    """Augmented plant, model, and adaptive-parameter dynamics."""
    x = z[0:3]
    xm = z[3:6]
    theta_hat = z[6:10]

    r = reference(t)
    phi = np.concatenate((x, np.array([r])))
    u = float(theta_hat @ phi)

    x_dot = AP @ x + B * (B0 * u)
    xm_dot = AM @ xm + B * (BM * r)

    e = x - xm
    scalar_error = float(e @ P @ B)
    theta_dot = -np.sign(B0) * (GAMMA @ phi) * scalar_error

    return np.concatenate((x_dot, xm_dot, theta_dot))


def main() -> None:
    z0 = np.zeros(10)
    z0[0:3] = np.array([0.5, -0.2, 0.1])

    t_eval = np.linspace(0.0, 60.0, 12001)
    solution = solve_ivp(
        dynamics,
        (t_eval[0], t_eval[-1]),
        z0,
        t_eval=t_eval,
        method="RK45",
        rtol=1.0e-8,
        atol=1.0e-10,
    )
    if not solution.success:
        raise RuntimeError(solution.message)

    t = solution.t
    x = solution.y[0:3, :]
    xm = solution.y[3:6, :]
    theta_hat = solution.y[6:10, :]

    r = np.array([reference(ti) for ti in t])
    u = np.sum(theta_hat[0:3, :] * x, axis=0) + theta_hat[3, :] * r
    e = x - xm
    error_norm = np.linalg.norm(e, axis=0)

    ideal_theta = np.array(
        [
            (-6.0 - (-1.0)) / B0,
            (-11.0 - (-2.0)) / B0,
            (-6.0 - (-1.2)) / B0,
            BM / B0,
        ]
    )

    print("Ideal matching parameters:", ideal_theta)
    print("Final adaptive parameters:", theta_hat[:, -1])
    print("Final tracking-error norm:", error_norm[-1])

    data = np.column_stack(
        (
            t,
            x.T,
            xm.T,
            theta_hat.T,
            u,
            error_norm,
        )
    )
    header = (
        "t,x1,x2,x3,xm1,xm2,xm3,"
        "theta1,theta2,theta3,theta_r,u,error_norm"
    )
    np.savetxt(
        "Chapter7_Lesson3_python_results.csv",
        data,
        delimiter=",",
        header=header,
        comments="",
    )

    plt.figure()
    plt.plot(t, x[0, :], label="plant x1")
    plt.plot(t, xm[0, :], "--", label="model xm1")
    plt.xlabel("Time (s)")
    plt.ylabel("Output/state")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.semilogy(t, np.maximum(error_norm, 1.0e-12))
    plt.xlabel("Time (s)")
    plt.ylabel("||e||")
    plt.grid(True)
    plt.tight_layout()

    plt.figure()
    for i in range(4):
        plt.plot(t, theta_hat[i, :], label=f"theta_hat[{i + 1}]")
        plt.axhline(ideal_theta[i], linestyle=":")
    plt.xlabel("Time (s)")
    plt.ylabel("Adaptive parameters")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, u)
    plt.xlabel("Time (s)")
    plt.ylabel("Control input u")
    plt.grid(True)
    plt.tight_layout()

    plt.show()


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library. It implements the augmented plant–model–parameter dynamics and a fourth-order Runge–Kutta integrator from scratch.

Chapter7_Lesson3.cpp

/*
Chapter7_Lesson3.cpp
Lyapunov-based MRAC for a third-order SISO canonical-form plant.

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

The program uses a fixed-step fourth-order Runge-Kutta integrator and writes
Chapter7_Lesson3_cpp_results.csv.
*/

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

namespace {

constexpr int STATE_SIZE = 10;
constexpr double B0 = 1.4;
constexpr double BM = 6.0;
constexpr double DT = 0.001;
constexpr double T_FINAL = 60.0;

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

constexpr std::array<std::array<double, 3>, 3> AP{ {
    { {0.0, 1.0, 0.0} },
    { {0.0, 0.0, 1.0} },
    { {-1.0, -2.0, -1.2} }
} };

constexpr std::array<std::array<double, 3>, 3> AM{ {
    { {0.0, 1.0, 0.0} },
    { {0.0, 0.0, 1.0} },
    { {-6.0, -11.0, -6.0} }
} };

constexpr std::array<std::array<double, 3>, 3> P{ {
    { {1.8166666666666667, 1.15, 0.0833333333333333} },
    { {1.15, 2.0083333333333333, 0.15} },
    { {0.0833333333333333, 0.15, 0.1083333333333333} }
} };

constexpr std::array<double, 4> GAMMA{ {8.0, 8.0, 8.0, 4.0} };

double reference(double t) {
    if (t < 8.0) {
        return 1.0;
    }
    if (t < 16.0) {
        return -0.5;
    }
    return 0.8;
}

std::array<double, 3> matVec(
    const std::array<std::array<double, 3>, 3>& matrix,
    const std::array<double, 3>& vector
) {
    std::array<double, 3> result{};
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            result[i] += matrix[i][j] * vector[j];
        }
    }
    return result;
}

State derivative(double t, const State& z) {
    State dz{};
    const std::array<double, 3> x{ {z[0], z[1], z[2]} };
    const std::array<double, 3> xm{ {z[3], z[4], z[5]} };
    const std::array<double, 4> theta{ {z[6], z[7], z[8], z[9]} };
    const double r = reference(t);

    const std::array<double, 4> phi{ {x[0], x[1], x[2], r} };
    double u = 0.0;
    for (int i = 0; i < 4; ++i) {
        u += theta[i] * phi[i];
    }

    auto xDot = matVec(AP, x);
    auto xmDot = matVec(AM, xm);
    xDot[2] += B0 * u;
    xmDot[2] += BM * r;

    const std::array<double, 3> e{ {
        x[0] - xm[0],
        x[1] - xm[1],
        x[2] - xm[2]
    } };

    // B = [0, 0, 1]^T, so P*B is the third column of P.
    const double scalarError =
        e[0] * P[0][2] + e[1] * P[1][2] + e[2] * P[2][2];

    dz[0] = xDot[0];
    dz[1] = xDot[1];
    dz[2] = xDot[2];
    dz[3] = xmDot[0];
    dz[4] = xmDot[1];
    dz[5] = xmDot[2];

    const double gainSign = (B0 >= 0.0) ? 1.0 : -1.0;
    for (int i = 0; i < 4; ++i) {
        dz[6 + i] = -gainSign * GAMMA[i] * phi[i] * scalarError;
    }
    return dz;
}

State addScaled(const State& a, const State& b, double scale) {
    State result{};
    for (int i = 0; i < STATE_SIZE; ++i) {
        result[i] = a[i] + scale * b[i];
    }
    return result;
}

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

    for (int i = 0; i < STATE_SIZE; ++i) {
        z[i] += (dt / 6.0) *
                (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(z[i])) {
            throw std::runtime_error("Non-finite state encountered.");
        }
    }
}

double controlInput(double t, const State& z) {
    const double r = reference(t);
    return z[6] * z[0] + z[7] * z[1] + z[8] * z[2] + z[9] * r;
}

double errorNorm(const State& z) {
    const double e1 = z[0] - z[3];
    const double e2 = z[1] - z[4];
    const double e3 = z[2] - z[5];
    return std::sqrt(e1 * e1 + e2 * e2 + e3 * e3);
}

}  // namespace

int main() {
    State z{};
    z[0] = 0.5;
    z[1] = -0.2;
    z[2] = 0.1;

    std::ofstream csv("Chapter7_Lesson3_cpp_results.csv");
    if (!csv) {
        std::cerr << "Cannot open output CSV file.\n";
        return 1;
    }

    csv << "t,x1,x2,x3,xm1,xm2,xm3,"
           "theta1,theta2,theta3,theta_r,u,error_norm\n";
    csv << std::setprecision(12);

    const int steps = static_cast<int>(std::round(T_FINAL / DT));
    for (int step = 0; step <= steps; ++step) {
        const double t = step * DT;

        if (step % 10 == 0) {
            csv << t;
            for (double value : z) {
                csv << ',' << value;
            }
            csv << ',' << controlInput(t, z)
                << ',' << errorNorm(z) << '\n';
        }

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

    const std::array<double, 4> idealTheta{ {
        (-6.0 - (-1.0)) / B0,
        (-11.0 - (-2.0)) / B0,
        (-6.0 - (-1.2)) / B0,
        BM / B0
    } };

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Ideal matching parameters: ";
    for (double value : idealTheta) {
        std::cout << value << ' ';
    }
    std::cout << "\nFinal adaptive parameters: ";
    for (int i = 6; i < 10; ++i) {
        std::cout << z[i] << ' ';
    }
    std::cout << "\nFinal tracking-error norm: " << errorNorm(z) << '\n';
    return 0;
}

The Java implementation also uses a self-contained RK4 method and the standard I/O library. It writes the same principal states, adaptive parameters, control signal, and error norm to a CSV file.

Chapter7_Lesson3.java

/*
Chapter7_Lesson3.java
Lyapunov-based MRAC for a third-order SISO canonical-form plant.

Build:
    javac Chapter7_Lesson3.java
Run:
    java Chapter7_Lesson3

The program uses fixed-step RK4 and writes
Chapter7_Lesson3_java_results.csv.
*/

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

public final class Chapter7_Lesson3 {
    private static final int STATE_SIZE = 10;
    private static final double B0 = 1.4;
    private static final double BM = 6.0;
    private static final double DT = 0.001;
    private static final double T_FINAL = 60.0;

    private static final double[][] AP = {
        {0.0, 1.0, 0.0},
        {0.0, 0.0, 1.0},
        {-1.0, -2.0, -1.2}
    };

    private static final double[][] AM = {
        {0.0, 1.0, 0.0},
        {0.0, 0.0, 1.0},
        {-6.0, -11.0, -6.0}
    };

    private static final double[][] P = {
        {1.8166666666666667, 1.15, 0.0833333333333333},
        {1.15, 2.0083333333333333, 0.15},
        {0.0833333333333333, 0.15, 0.1083333333333333}
    };

    private static final double[] GAMMA = {8.0, 8.0, 8.0, 4.0};

    private Chapter7_Lesson3() {
    }

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

    private static double[] matVec(double[][] matrix, double[] vector) {
        double[] result = new double[3];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                result[i] += matrix[i][j] * vector[j];
            }
        }
        return result;
    }

    private static double[] derivative(double t, double[] z) {
        double[] dz = new double[STATE_SIZE];
        double[] x = {z[0], z[1], z[2]};
        double[] xm = {z[3], z[4], z[5]};
        double[] theta = {z[6], z[7], z[8], z[9]};
        double r = reference(t);
        double[] phi = {x[0], x[1], x[2], r};

        double u = 0.0;
        for (int i = 0; i < 4; ++i) {
            u += theta[i] * phi[i];
        }

        double[] xDot = matVec(AP, x);
        double[] xmDot = matVec(AM, xm);
        xDot[2] += B0 * u;
        xmDot[2] += BM * r;

        double[] e = {
            x[0] - xm[0],
            x[1] - xm[1],
            x[2] - xm[2]
        };

        // B = [0, 0, 1]^T; therefore P*B is column 3 of P.
        double scalarError =
            e[0] * P[0][2] + e[1] * P[1][2] + e[2] * P[2][2];

        dz[0] = xDot[0];
        dz[1] = xDot[1];
        dz[2] = xDot[2];
        dz[3] = xmDot[0];
        dz[4] = xmDot[1];
        dz[5] = xmDot[2];

        double gainSign = Math.signum(B0);
        for (int i = 0; i < 4; ++i) {
            dz[6 + i] = -gainSign * GAMMA[i] * phi[i] * scalarError;
        }
        return dz;
    }

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

    private static void rk4Step(double t, double dt, double[] z) {
        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));

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

    private static double controlInput(double t, double[] z) {
        double r = reference(t);
        return z[6] * z[0] + z[7] * z[1] + z[8] * z[2] + z[9] * r;
    }

    private static double errorNorm(double[] z) {
        double e1 = z[0] - z[3];
        double e2 = z[1] - z[4];
        double e3 = z[2] - z[5];
        return Math.sqrt(e1 * e1 + e2 * e2 + e3 * e3);
    }

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

        double[] z = new double[STATE_SIZE];
        z[0] = 0.5;
        z[1] = -0.2;
        z[2] = 0.1;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter7_Lesson3_java_results.csv"))) {
            writer.write(
                "t,x1,x2,x3,xm1,xm2,xm3,"
                + "theta1,theta2,theta3,theta_r,u,error_norm\n"
            );

            int steps = (int) Math.round(T_FINAL / DT);
            for (int step = 0; step <= steps; ++step) {
                double t = step * DT;

                if (step % 10 == 0) {
                    writer.write(String.format(Locale.US, "%.12g", t));
                    for (double value : z) {
                        writer.write(String.format(Locale.US, ",%.12g", value));
                    }
                    writer.write(String.format(
                        Locale.US,
                        ",%.12g,%.12g%n",
                        controlInput(t, z),
                        errorNorm(z)
                    ));
                }

                if (step < steps) {
                    rk4Step(t, DT, z);
                }
            }
        } catch (IOException exception) {
            System.err.println("I/O error: " + exception.getMessage());
            System.exit(1);
        }

        double[] idealTheta = {
            (-6.0 - (-1.0)) / B0,
            (-11.0 - (-2.0)) / B0,
            (-6.0 - (-1.2)) / B0,
            BM / B0
        };

        System.out.print("Ideal matching parameters: ");
        for (double value : idealTheta) {
            System.out.printf(Locale.US, "%.6f ", value);
        }
        System.out.print("\nFinal adaptive parameters: ");
        for (int i = 6; i < 10; ++i) {
            System.out.printf(Locale.US, "%.6f ", z[i]);
        }
        System.out.printf(
            Locale.US,
            "%nFinal tracking-error norm: %.9g%n",
            errorNorm(z)
        );
    }
}

The MATLAB script uses ode45 and therefore requires no control-design toolbox for the simulation. When Control System Toolbox is installed, the commented lyap check can be enabled to verify the supplied matrix \(\mathbf{P}\).

Chapter7_Lesson3.m

% Chapter7_Lesson3.m
% Lyapunov-based MRAC for a third-order SISO canonical-form plant.
%
% The script integrates the plant, reference model, and adaptive law with
% ode45, plots the principal signals, and writes
% Chapter7_Lesson3_matlab_results.csv.

clear;
clc;
close all;

p.AP = [0 1 0; 0 0 1; -1 -2 -1.2];
p.AM = [0 1 0; 0 0 1; -6 -11 -6];
p.B = [0; 0; 1];
p.b0 = 1.4;
p.bm = 6.0;
p.P = [1.8166666666666667, 1.15, 0.0833333333333333; ...
       1.15, 2.0083333333333333, 0.15; ...
       0.0833333333333333, 0.15, 0.1083333333333333];
p.Gamma = diag([8, 8, 8, 4]);

% Optional verification when Control System Toolbox is available:
% Q = eye(3);
% Pcheck = lyap(p.AM', Q);
% disp(norm(Pcheck - p.P, 'fro'));

z0 = zeros(10, 1);
z0(1:3) = [0.5; -0.2; 0.1];

tspan = linspace(0, 60, 12001);
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, z] = ode45(@(t, z) augmentedDynamics(t, z, p), tspan, z0, options);

x = z(:, 1:3);
xm = z(:, 4:6);
thetaHat = z(:, 7:10);
r = arrayfun(@referenceSignal, t);
u = sum(thetaHat(:, 1:3) .* x, 2) + thetaHat(:, 4) .* r;
e = x - xm;
errorNorm = vecnorm(e, 2, 2);

idealTheta = [(-6 - (-1)) / p.b0; ...
              (-11 - (-2)) / p.b0; ...
              (-6 - (-1.2)) / p.b0; ...
              p.bm / p.b0];

fprintf('Ideal matching parameters:\n');
disp(idealTheta.');
fprintf('Final adaptive parameters:\n');
disp(thetaHat(end, :));
fprintf('Final tracking-error norm: %.9g\n', errorNorm(end));

results = table(t, x(:,1), x(:,2), x(:,3), ...
    xm(:,1), xm(:,2), xm(:,3), ...
    thetaHat(:,1), thetaHat(:,2), thetaHat(:,3), thetaHat(:,4), ...
    u, errorNorm, ...
    'VariableNames', {'t','x1','x2','x3','xm1','xm2','xm3', ...
    'theta1','theta2','theta3','theta_r','u','error_norm'});
writetable(results, 'Chapter7_Lesson3_matlab_results.csv');

figure;
plot(t, x(:,1), 'LineWidth', 1.2);
hold on;
plot(t, xm(:,1), '--', 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Output/state');
legend('plant x_1', 'model x_{m1}', 'Location', 'best');
title('MRAC tracking');

figure;
semilogy(t, max(errorNorm, 1e-12), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('||e||_2');
title('Tracking-error norm');

figure;
plot(t, thetaHat, 'LineWidth', 1.1);
hold on;
for k = 1:4
    yline(idealTheta(k), ':');
end
grid on;
xlabel('Time (s)');
ylabel('Adaptive parameters');
legend('\hat{k}_{x1}', '\hat{k}_{x2}', '\hat{k}_{x3}', ...
    '\hat{k}_r', 'Location', 'best');
title('Parameter estimates');

figure;
plot(t, u, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Control input u');
title('Adaptive control effort');

function dz = augmentedDynamics(t, z, p)
    x = z(1:3);
    xm = z(4:6);
    thetaHat = z(7:10);

    r = referenceSignal(t);
    phi = [x; r];
    u = thetaHat.' * phi;

    xDot = p.AP * x + p.B * (p.b0 * u);
    xmDot = p.AM * xm + p.B * (p.bm * r);

    e = x - xm;
    scalarError = e.' * p.P * p.B;
    thetaDot = -sign(p.b0) * p.Gamma * phi * scalarError;

    dz = [xDot; xmDot; thetaDot];
end

function r = referenceSignal(t)
    if t < 8
        r = 1.0;
    elseif t < 16
        r = -0.5;
    else
        r = 0.8;
    end
end

The following MATLAB script creates a Simulink model programmatically. The model contains separate plant and reference-model state-space blocks, a regressor multiplexer, a dot-product adaptive controller, a tracking error path, the scalar signal \(s=\mathbf{e}^{T}\mathbf{P}\mathbf{B}\), and a vector integrator for the adaptive parameters. Running the script generates Chapter7_Lesson3_MRAC.slx.

Chapter7_Lesson3_Simulink.m

% Chapter7_Lesson3_Simulink.m
% Programmatically builds and simulates a Simulink implementation of the
% third-order Lyapunov-based MRAC example from Chapter 7, Lesson 3.
%
% Requirements:
%   MATLAB
%   Simulink
%
% Output:
%   Chapter7_Lesson3_MRAC.slx
%   Chapter7_Lesson3_simulink_results.mat

clear;
clc;

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

AP = [0 1 0; 0 0 1; -1 -2 -1.2];
AM = [0 1 0; 0 0 1; -6 -11 -6];
B = [0; 0; 1];
b0 = 1.4;
bm = 6.0;
P = [1.8166666666666667, 1.15, 0.0833333333333333; ...
     1.15, 2.0083333333333333, 0.15; ...
     0.0833333333333333, 0.15, 0.1083333333333333];
Gamma = diag([8, 8, 8, 4]);
PB = P * B;

% Piecewise-constant reference as a timeseries.
tRef = [0; 7.999; 8; 15.999; 16; 60];
rRef = [1; 1; -0.5; -0.5; 0.8; 0.8];
referenceData = timeseries(rRef, tRef);

new_system(model);
open_system(model);
set_param(model, ...
    'Solver', 'ode45', ...
    'StopTime', '60', ...
    'RelTol', '1e-8', ...
    'AbsTol', '1e-10');

add_block('simulink/Sources/From Workspace', [model '/Reference'], ...
    'VariableName', 'referenceData', ...
    'Position', [40 220 150 250]);

add_block('simulink/Continuous/State-Space', [model '/Plant'], ...
    'A', 'AP', ...
    'B', 'B*b0', ...
    'C', 'eye(3)', ...
    'D', 'zeros(3,1)', ...
    'X0', '[0.5; -0.2; 0.1]', ...
    'Position', [590 80 720 140]);

add_block('simulink/Continuous/State-Space', [model '/Reference Model'], ...
    'A', 'AM', ...
    'B', 'B*bm', ...
    'C', 'eye(3)', ...
    'D', 'zeros(3,1)', ...
    'X0', 'zeros(3,1)', ...
    'Position', [300 290 450 350]);

add_block('simulink/Signal Routing/Mux', [model '/Regressor Mux'], ...
    'Inputs', '2', ...
    'DisplayOption', 'bar', ...
    'Position', [300 100 305 180]);

add_block('simulink/Math Operations/Dot Product', [model '/Control Dot Product'], ...
    'Position', [470 100 530 150]);

add_block('simulink/Math Operations/Sum', [model '/Tracking Error'], ...
    'Inputs', '+-', ...
    'Position', [790 170 820 230]);

add_block('simulink/Math Operations/Gain', [model '/PB Transpose'], ...
    'Gain', 'PB''', ...
    'Multiplication', 'Matrix(K*u)', ...
    'Position', [870 170 960 230]);

add_block('simulink/Math Operations/Product', [model '/Scalar Times Regressor'], ...
    'Inputs', '**', ...
    'Multiplication', 'Element-wise(.*)', ...
    'Position', [480 210 550 270]);

add_block('simulink/Math Operations/Gain', [model '/Adaptation Gain'], ...
    'Gain', '-sign(b0)*Gamma', ...
    'Multiplication', 'Matrix(K*u)', ...
    'Position', [600 210 720 270]);

add_block('simulink/Continuous/Integrator', [model '/Parameter Integrator'], ...
    'InitialCondition', 'zeros(4,1)', ...
    'Position', [770 260 810 310]);

add_block('simulink/Sinks/To Workspace', [model '/Plant State Output'], ...
    'VariableName', 'xSim', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [870 60 990 90]);

add_block('simulink/Sinks/To Workspace', [model '/Model State Output'], ...
    'VariableName', 'xmSim', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [510 320 640 350]);

add_block('simulink/Sinks/To Workspace', [model '/Parameter Output'], ...
    'VariableName', 'thetaSim', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [870 280 1000 310]);

add_block('simulink/Sinks/To Workspace', [model '/Control Output'], ...
    'VariableName', 'uSim', ...
    'SaveFormat', 'Structure With Time', ...
    'Position', [590 150 710 180]);

% Reference to model and regressor.
add_line(model, 'Reference/1', 'Reference Model/1', 'autorouting', 'on');
add_line(model, 'Reference/1', 'Regressor Mux/2', 'autorouting', 'on');

% Plant state to regressor, tracking error, and logging.
add_line(model, 'Plant/1', 'Regressor Mux/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Tracking Error/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Plant State Output/1', 'autorouting', 'on');

% Model state to error and logging.
add_line(model, 'Reference Model/1', 'Tracking Error/2', 'autorouting', 'on');
add_line(model, 'Reference Model/1', 'Model State Output/1', 'autorouting', 'on');

% Adaptive control u = theta_hat^T phi.
add_line(model, 'Regressor Mux/1', 'Control Dot Product/1', 'autorouting', 'on');
add_line(model, 'Parameter Integrator/1', 'Control Dot Product/2', 'autorouting', 'on');
add_line(model, 'Control Dot Product/1', 'Plant/1', 'autorouting', 'on');
add_line(model, 'Control Dot Product/1', 'Control Output/1', 'autorouting', 'on');

% s = e^T P B and theta_dot = -sign(b0) Gamma phi s.
add_line(model, 'Tracking Error/1', 'PB Transpose/1', 'autorouting', 'on');
add_line(model, 'PB Transpose/1', 'Scalar Times Regressor/2', 'autorouting', 'on');
add_line(model, 'Regressor Mux/1', 'Scalar Times Regressor/1', 'autorouting', 'on');
add_line(model, 'Scalar Times Regressor/1', 'Adaptation Gain/1', 'autorouting', 'on');
add_line(model, 'Adaptation Gain/1', 'Parameter Integrator/1', 'autorouting', 'on');
add_line(model, 'Parameter Integrator/1', 'Parameter Output/1', 'autorouting', 'on');

save_system(model);
simulationOutput = sim(model);

save('Chapter7_Lesson3_simulink_results.mat', ...
    'simulationOutput', 'AP', 'AM', 'B', 'b0', 'bm', 'P', 'Gamma');

disp(['Created and simulated ' model '.slx']);

The Wolfram Language implementation uses native list and matrix operations, a from-scratch RK4 recursion, CSV export, and ListLinePlot/ListLogPlot visualization.

Chapter7_Lesson3.nb

(* Content-type: application/vnd.wolfram.mathematica *)

Notebook[{
  Cell[
    "Chapter 7, Lesson 3: Lyapunov-Based MRAC for n-th Order Systems",
    "Title"
  ],
  Cell[
    "Keep Chapter7_Lesson3.wl in the same directory as this notebook. Evaluate the input cell below to run the complete simulation, export the CSV results, and generate the plots.",
    "Text"
  ],
  Cell[
    BoxData[
      RowBox[{
        "Get", "[",
        RowBox[{
          "FileNameJoin", "[",
          RowBox[{
            "{",
            RowBox[{
              RowBox[{"NotebookDirectory", "[", "]"}],
              ",",
              "\"Chapter7_Lesson3.wl\""
            }],
            "}"
          }],
          "]"
        }],
        "]"
      }]
    ],
    "Input"
  ]
},
WindowSize -> {1100, 800},
WindowMargins -> { {Automatic, 50}, {Automatic, 50} },
StyleDefinitions -> "Default.nb"
]

The notebook loads the complete Wolfram Language implementation stored beside it:

Chapter7_Lesson3.wl

(* Chapter7_Lesson3.wl
   Wolfram Language notebook source for Lyapunov-based MRAC applied to a
   third-order SISO plant in controllable canonical form.

   Paste this content into a Wolfram notebook and evaluate all cells.
*)

ClearAll["Global`*"];

ap = { {0., 1., 0.}, {0., 0., 1.}, {-1., -2., -1.2} };
am = { {0., 1., 0.}, {0., 0., 1.}, {-6., -11., -6.} };
b = {0., 0., 1.};
b0 = 1.4;
bm = 6.0;

p = {
   {1.8166666666666667, 1.15, 0.0833333333333333},
   {1.15, 2.0083333333333333, 0.15},
   {0.0833333333333333, 0.15, 0.1083333333333333}
   };

gamma = DiagonalMatrix[{8., 8., 8., 4.}];
dt = 0.001;
tFinal = 60.0;

reference[t_] := Piecewise[{ {1.0, t < 8.0}, {-0.5, t < 16.0} }, 0.8];

derivative[t_, z_] := Module[
  {x, xm, thetaHat, r, phi, u, xDot, xmDot, e, scalarError, thetaDot},
  x = z[[1 ;; 3]];
  xm = z[[4 ;; 6]];
  thetaHat = z[[7 ;; 10]];
  r = reference[t];
  phi = Join[x, {r}];
  u = thetaHat . phi;

  xDot = ap . x + b (b0 u);
  xmDot = am . xm + b (bm r);

  e = x - xm;
  scalarError = e . p . b;
  thetaDot = -Sign[b0] gamma . phi scalarError;

  Join[xDot, xmDot, thetaDot]
  ];

rk4Step[t_, z_] := Module[{k1, k2, k3, k4},
  k1 = derivative[t, z];
  k2 = derivative[t + dt/2, z + (dt/2) k1];
  k3 = derivative[t + dt/2, z + (dt/2) k2];
  k4 = derivative[t + dt, z + dt k3];
  z + (dt/6) (k1 + 2 k2 + 2 k3 + k4)
  ];

z0 = ConstantArray[0., 10];
z0[[1 ;; 3]] = {0.5, -0.2, 0.1};

steps = Round[tFinal/dt];
trajectory = NestList[
   Function[pair, {pair[[1]] + dt, rk4Step[pair[[1]], pair[[2]]]}],
   {0., z0},
   steps
   ];

sampled = trajectory[[1 ;; ;; 10]];
times = sampled[[All, 1]];
states = sampled[[All, 2]];
xData = states[[All, 1 ;; 3]];
xmData = states[[All, 4 ;; 6]];
thetaData = states[[All, 7 ;; 10]];
rData = reference /@ times;
uData = MapThread[#1[[1 ;; 3]] . #2 + #1[[4]] #3 &,
   {thetaData, xData, rData}];
errorNorm = Norm /@ (xData - xmData);

idealTheta = {
   (-6.0 - (-1.0))/b0,
   (-11.0 - (-2.0))/b0,
   (-6.0 - (-1.2))/b0,
   bm/b0
   };

Print["Ideal matching parameters: ", idealTheta];
Print["Final adaptive parameters: ", Last[thetaData]];
Print["Final tracking-error norm: ", Last[errorNorm]];

header = {
   "t", "x1", "x2", "x3", "xm1", "xm2", "xm3",
   "theta1", "theta2", "theta3", "theta_r", "u", "error_norm"
   };

rows = MapThread[
   Join[{ #1 }, #2, #3, #4, { #5, #6 }] &,
   {times, xData, xmData, thetaData, uData, errorNorm}
   ];

Export["Chapter7_Lesson3_mathematica_results.csv",
  Prepend[rows, header]];

trackingPlot = ListLinePlot[
   {
    Transpose[{times, xData[[All, 1]]}],
    Transpose[{times, xmData[[All, 1]]}]
    },
   PlotLegends -> {"plant x1", "model xm1"},
   AxesLabel -> {"Time (s)", "Output/state"},
   PlotRange -> All,
   ImageSize -> Large
   ];

errorPlot = ListLogPlot[
   Transpose[{times, Clip[errorNorm, {10^-12, Infinity}]}],
   AxesLabel -> {"Time (s)", "||e||"},
   PlotRange -> All,
   ImageSize -> Large
   ];

parameterPlot = ListLinePlot[
   Table[Transpose[{times, thetaData[[All, k]]}], {k, 1, 4}],
   PlotLegends -> {"theta1", "theta2", "theta3", "theta_r"},
   AxesLabel -> {"Time (s)", "Adaptive parameters"},
   PlotRange -> All,
   ImageSize -> Large
   ];

controlPlot = ListLinePlot[
   Transpose[{times, uData}],
   AxesLabel -> {"Time (s)", "Control input u"},
   PlotRange -> All,
   ImageSize -> Large
   ];

Column[{trackingPlot, errorPlot, parameterPlot, controlPlot}]

A correct implementation should be checked through mathematical invariants rather than through a plot alone:

1. verify numerically that \(\mathbf{P}=\mathbf{P}^{T}>0\);
2. verify that \(\mathbf{A}_m^{T}\mathbf{P} +\mathbf{P}\mathbf{A}_m+\mathbf{Q}\) is numerically zero;
3. confirm that the plant and model use the same command \(r(t)\);
4. compute \(s=\mathbf{e}^{T}\mathbf{P}\mathbf{B}\) with the correct sign and orientation;
5. confirm that the input-gain sign in the update law is \(\operatorname{sgn}(b_0)\);
6. reduce the numerical step or integration tolerances and verify that the trajectories do not materially change;
7. verify that the tracking-error norm decreases toward zero after each command transient.

The adaptive parameter trajectories need not approach the plotted ideal values for this command. That observation is consistent with the theorem, which guarantees asymptotic tracking but does not yet impose a separate condition for exact parameter recovery.

• Using output error instead of state error: this full-state proof requires \(\mathbf{e}=\mathbf{x}-\mathbf{x}_m\), not merely \(y-y_m\).
• Dropping the factor \(|b_0|\): the parameter-energy term must be scaled so that the sign-based adaptive law cancels the error cross term.
• Solving the wrong Lyapunov equation: the equation is formed with the reference matrix, not the uncertain plant matrix.
• Reversing the error definition: changing from \(\mathbf{x}-\mathbf{x}_m\) to \(\mathbf{x}_m-\mathbf{x}\) requires a corresponding change in the adaptive-law sign.
• Assuming parameter convergence from \(\dot V\leq 0\): negative semidefiniteness establishes boundedness and supports the tracking proof, but it does not by itself identify every ideal parameter.
• Ignoring continuous-to-discrete effects: a large fixed numerical step can destroy the behavior predicted by a continuous-time Lyapunov proof.

Problem 1 — Matching Parameters for a Fourth-Order Plant. Consider

\[ \mathbf{A}_p = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -2 & -3 & -4 & -1 \end{bmatrix}, \qquad b_0=2. \]

The desired model has last row \([\, -10\;-35\;-50\;-24\,]\) and \(b_m=10\). Find \(\mathbf{k}_x^{*}\) and \(k_r^{*}\).

Solution. From the matching equation,

\[ 2\mathbf{k}_x^{*T} = \begin{bmatrix} -10 & -35 & -50 & -24 \end{bmatrix} - \begin{bmatrix} -2 & -3 & -4 & -1 \end{bmatrix}. \]

\[ \boxed{ \mathbf{k}_x^{*} = \begin{bmatrix} -4 \\ -16 \\ -23 \\ -11.5 \end{bmatrix}, \qquad k_r^{*}=\frac{10}{2}=5 }. \]

Problem 2 — Derive the Vector Adaptive Law. Starting from

\[ \dot{\mathbf{e} } = \mathbf{A}_m\mathbf{e} + \mathbf{B}b_0 \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}, \]

and \(V=\mathbf{e}^{T}\mathbf{P}\mathbf{e} +|b_0|\widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1}\widetilde{\boldsymbol{\theta} }\), derive the adaptive law.

Solution. Differentiation gives

\[ \dot V = -\mathbf{e}^{T}\mathbf{Q}\mathbf{e} + 2b_0s\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} + 2|b_0| \widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \dot{\widehat{\boldsymbol{\theta} } }, \qquad s=\mathbf{e}^{T}\mathbf{P}\mathbf{B}. \]

To cancel the two terms linear in the parameter error, require

\[ |b_0| \boldsymbol{\Gamma}^{-1} \dot{\widehat{\boldsymbol{\theta} } } = -b_0\boldsymbol{\phi}s. \]

Therefore,

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } = -\operatorname{sgn}(b_0) \boldsymbol{\Gamma}\boldsymbol{\phi}s }. \]

Problem 3 — Explicit Error-Energy Bound. Let \(\lambda_{\min}(\mathbf{Q})=0.8\) and \(V(0)=12\). Find an upper bound on \(\int_{0}^{\infty}\|\mathbf{e}(t)\|^2dt\).

Solution.

\[ \int_{0}^{\infty}\|\mathbf{e}(t)\|^2dt \leq \frac{V(0)}{\lambda_{\min}(\mathbf{Q})} = \frac{12}{0.8} = \boxed{15}. \]

Problem 4 — Effect of Reversing the Error Definition. Suppose a programmer defines \(\bar{\mathbf{e} }=\mathbf{x}_m-\mathbf{x}=-\mathbf{e}\) but leaves the update law unchanged. Determine the correct adaptive law in terms of \(\bar{\mathbf{e} }\).

Solution. Since

\[ \mathbf{e}^{T}\mathbf{P}\mathbf{B} = -\bar{\mathbf{e} }^{T}\mathbf{P}\mathbf{B}, \]

the original law becomes

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } } = +\operatorname{sgn}(b_0) \boldsymbol{\Gamma} \boldsymbol{\phi} \bar{\mathbf{e} }^{T}\mathbf{P}\mathbf{B} }. \]

Leaving the original minus sign while reversing the error definition destroys the cancellation in \(\dot V\).

Problem 5 — Why the Tracking Proof Does Not Prove Parameter Recovery. From \(\dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}\), explain why \(\widetilde{\boldsymbol{\theta} }(t)\to\mathbf{0}\) does not follow immediately.

Solution. The derivative contains no negative-definite term in the parameter error:

\[ \dot V = -\mathbf{e}^{T}\mathbf{Q}\mathbf{e}, \]

so the proof forces tracking-error energy to decrease but only bounds the parameter error through the value of \(V\). When the tracking error becomes zero, the update may stop at a parameter vector that reproduces the observed trajectory without necessarily equaling the unique ideal vector. Additional excitation conditions are required for a separate parameter-convergence theorem.

For an \(n\)-th order SISO plant in controllable canonical form, the matching equations transform the adaptive state-feedback problem into the linearly parameterized error system \(\dot{\mathbf{e} }=\mathbf{A}_m\mathbf{e} +\mathbf{B}b_0\widetilde{\boldsymbol{\theta} }^{T} \boldsymbol{\phi}\). A composite Lyapunov function containing state and parameter errors leads directly to the update law \(\dot{\widehat{\boldsymbol{\theta} } } =-\operatorname{sgn}(b_0)\boldsymbol{\Gamma} \boldsymbol{\phi}\mathbf{e}^{T}\mathbf{P}\mathbf{B}\). The cross term cancels exactly, yielding \(\dot V=-\mathbf{e}^{T}\mathbf{Q}\mathbf{e}\). Lyapunov boundedness, square integrability, and Barbalat's lemma then establish asymptotic state tracking. The same mathematical design was implemented in six computational forms for a third-order example.