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

• simulate a first-order Lyapunov-based MRAC loop as one coupled nonlinear differential equation;
• verify the matching parameters and the closed-loop tracking-error model;
• explain why an adaptive gain moves, slows, reverses direction, or becomes nearly constant;
• distinguish tracking convergence from controller-parameter convergence;
• check whether numerical results are consistent with the continuous-time Lyapunov analysis; and
• reproduce the experiment in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The central question is not merely whether the plant output follows the reference-model output. We also ask:

\[ \begin{gathered} \text{What information about the unknown plant is actually encoded} \\ \text{in the observed evolution of } \hat{\boldsymbol{\theta} }(t)\text{?} \end{gathered} \]

The answer requires simultaneous examination of the tracking error, the regressor signals, the adaptive-law direction, and the command waveform.

Consider the scalar plant used throughout this chapter:

\[ \dot{y}(t)=-a y(t)+b u(t), \qquad a > 0, \quad b > 0, \]

where the numerical plant used only by the simulator is

\[ a=1.2, \qquad b=0.8. \]

The controller is not assumed to know these two values. It knows only that the control direction is positive, namely \( b > 0 \). The desired first-order behavior is specified by

\[ \dot{y}_m(t)=-a_m y_m(t)+b_m r(t), \qquad a_m=2, \quad b_m=2. \]

Therefore the reference model has unit steady-state gain \( b_m/a_m=1 \) and time constant \( 1/a_m=0.5\,\text{s} \). Use the direct adaptive control law

\[ u(t)=\hat{\theta}_y(t)y(t)+\hat{\theta}_r(t)r(t). \]

If constant ideal parameters existed and were known, coefficient matching would require

\[ -a+b\theta_y^*=-a_m, \qquad b\theta_r^*=b_m. \]

Solving these equations gives

\[ \theta_y^*=\frac{a-a_m}{b}=-1, \qquad \theta_r^*=\frac{b_m}{b}=2.5. \]

These ideal values are used only for post-simulation interpretation. The adaptive controller starts from

\[ \hat{\theta}_y(0)=0, \qquad \hat{\theta}_r(0)=0.5, \qquad y(0)=y_m(0)=0. \]

Define the tracking and parameter errors as

\[ e=y-y_m, \qquad \widetilde{\theta}_y=\hat{\theta}_y-\theta_y^*, \qquad \widetilde{\theta}_r=\hat{\theta}_r-\theta_r^*. \]

Substituting the adaptive control law into the plant gives

\[ \dot{y}=(-a+b\hat{\theta}_y)y+b\hat{\theta}_r r. \]

Add and subtract the matched coefficients. Because \( -a+b\theta_y^*=-a_m \) and \( b\theta_r^*=b_m \), the tracking-error equation is

\[ \dot{e}=-a_m e+b\widetilde{\theta}_y y +b\widetilde{\theta}_r r. \]

Introduce the regressor and parameter-error vectors

\[ \boldsymbol{\phi}=\begin{bmatrix}y\\r\end{bmatrix}, \qquad \widetilde{\boldsymbol{\theta} }= \begin{bmatrix}\widetilde{\theta}_y\\ \widetilde{\theta}_r\end{bmatrix}. \]

Then the compact error model is

\[ \dot{e}=-a_m e+b\widetilde{\boldsymbol{\theta} }^T \boldsymbol{\phi}. \]

Select positive adaptation gains \( \gamma_y=\gamma_r=4 \) and use

\[ \dot{\hat{\theta} }_y=-\gamma_y y e, \qquad \dot{\hat{\theta} }_r=-\gamma_r r e. \]

Since the ideal parameters are constant, \( \dot{\widetilde{\theta} }_i= \dot{\hat{\theta} }_i \). Consider

\[ V=\frac{1}{2}e^2+ \frac{b}{2\gamma_y}\widetilde{\theta}_y^2+ \frac{b}{2\gamma_r}\widetilde{\theta}_r^2. \]

Differentiating along the closed-loop trajectories yields

\[ \begin{aligned} \dot{V} &=e\left(-a_m e+b\widetilde{\theta}_y y +b\widetilde{\theta}_r r\right) +\frac{b}{\gamma_y}\widetilde{\theta}_y \dot{\hat{\theta} }_y +\frac{b}{\gamma_r}\widetilde{\theta}_r \dot{\hat{\theta} }_r\\ &=-a_m e^2+b\widetilde{\theta}_y y e +b\widetilde{\theta}_r r e -b\widetilde{\theta}_y y e -b\widetilde{\theta}_r r e\\ &=-a_m e^2\le 0. \end{aligned} \]

Hence \( V(t) \), \( e(t) \), and both parameter errors remain bounded. Moreover,

\[ \int_0^\infty e^2(\tau)\,d\tau \le \frac{V(0)}{a_m} < \infty. \]

For bounded reference commands, all right-hand-side signals are bounded, so \( \dot{e} \) is bounded. Barbalat's lemma, introduced in Chapter 3, then gives

\[ e(t) → 0 \qquad \text{as} \qquad t → \infty. \]

This conclusion guarantees asymptotic tracking for the ideal continuous-time equations. It does not yet prove that each parameter error converges to zero. The simulation will make that distinction visible.

flowchart TD
  R["Command r"] --> RM["Reference model: ym_dot = -am ym + bm r"]
  R --> C["Adaptive controller: u = theta_y y + theta_r r"]
  Y["Plant output y"] --> C
  C --> P["Unknown plant: y_dot = -a y + b u"]
  P --> Y
  Y --> E["Tracking error e = y - ym"]
  RM --> E
  E --> AY["theta_y_dot = -gamma_y y e"]
  E --> AR["theta_r_dot = -gamma_r r e"]
  AY --> C
  AR --> C
        

The command is intentionally divided into three intervals:

\[ r(t)=\begin{cases} 1, & 0\le t < 10,\\ -0.5, & 10\le t < 20,\\ 0.8\sin(1.1t)+0.45\sin(0.37t), & 20\le t\le 60. \end{cases} \]

The first two intervals test step tracking and expose a parameter ambiguity under nearly constant steady-state signals. The final interval contains two distinct time scales and causes the regressor pair \( [y(t),r(t)]^T \) to change direction repeatedly. Later chapters formalize signal conditions for parameter convergence; here we use only direct simulation-based interpretation.

The complete state integrated by the numerical solver is

\[ \mathbf{x}= \begin{bmatrix} y & y_m & \hat{\theta}_y & \hat{\theta}_r & J_e \end{bmatrix}^T, \qquad \dot{J}_e=e^2, \quad J_e(0)=0. \]

Thus \( J_e(t)=\int_0^t e^2(\tau)d\tau \) provides a cumulative tracking-performance index.

The adaptive law can be written as

\[ \dot{\hat{\boldsymbol{\theta} } }=- \begin{bmatrix}\gamma_y&0\\0&\gamma_r\end{bmatrix} \boldsymbol{\phi}e. \]

Therefore parameter motion requires both a nonzero tracking error and a nonzero regressor component. In particular,

\[ e=0 \quad \Longrightarrow \quad \dot{\hat{\theta} }_y=\dot{\hat{\theta} }_r=0. \]

This explains the common observation that parameters become nearly constant as soon as tracking becomes accurate. It does not follow that the frozen values are the unique ideal values.

To see the ambiguity analytically, let a nonzero constant command be applied long enough that \( y\approx y_m \) and both derivatives are approximately zero. Define the reference-model steady-state ratio

\[ \rho=\frac{b_m}{a_m}, \qquad y_m\approx \rho r. \]

Substitution into the steady-state plant equation gives

\[ 0=-a\rho r+b\left(\hat{\theta}_y\rho r+ \hat{\theta}_r r\right). \]

For \( r\ne0 \), divide by \( r \):

\[ \rho\hat{\theta}_y+\hat{\theta}_r= \frac{a\rho}{b}. \]

This is one linear equation in two unknown controller parameters. With the present numerical values, \( \rho=1 \), so every pair on

\[ \hat{\theta}_y+\hat{\theta}_r=1.5 \]

produces the correct constant steady-state output. The ideal pair \( (-1,2.5) \) belongs to this line, but it is not the only pair on it. In the supplied simulation, the estimates near \( t=10\,\text{s} \) are approximately

\[ \hat{\theta}_y\approx-0.067, \qquad \hat{\theta}_r\approx1.565, \qquad \hat{\theta}_y+\hat{\theta}_r\approx1.498. \]

Tracking is already accurate because the required combination is nearly correct, even though the individual parameters are still far from \( -1 \) and \( 2.5 \).

6.1 Tracking output and transient error

At \( t=0 \), the controller has insufficient feedforward gain, so the plant initially responds more slowly than the reference model. The nonzero error drives both updates. Each command change creates a new transient, after which the Lyapunov law again reduces the tracking error.

With a Runge–Kutta step of \( h=0.005\,\text{s} \), the representative root-mean-square errors are

6.2 Direction and speed of parameter motion

The instantaneous rates satisfy

\[ \left|\dot{\hat{\theta} }_y\right|= \gamma_y|y||e|, \qquad \left|\dot{\hat{\theta} }_r\right|= \gamma_r|r||e|. \]

Thus a parameter changes rapidly when its associated regressor is large and the tracking error is large. It changes slowly when either factor is small. The sign is also directly readable:

\[ \operatorname{sgn}\!\left(\dot{\hat{\theta} }_y\right) =-\operatorname{sgn}(ye), \qquad \operatorname{sgn}\!\left(\dot{\hat{\theta} }_r\right) =-\operatorname{sgn}(re). \]

6.3 Parameter values at the end of the experiment

A representative run gives

\[ \hat{\theta}_y(60)\approx-0.944, \qquad \hat{\theta}_r(60)\approx2.438, \qquad e(60)\approx8.13\times10^{-3}. \]

These values are close to the ideal pair but not exactly equal. The finite simulation horizon, finite numerical step, and remaining small tracking error all contribute. More importantly, the continuous-time stability proof guarantees tracking and boundedness; it does not state that an arbitrary command must uniquely determine every parameter.

6.4 Lyapunov-function interpretation

For the selected initial conditions,

\[ V(0)=\frac{0.8}{2(4)} \left[(0-(-1))^2+(0.5-2.5)^2\right]=0.5. \]

The continuous-time theory predicts a nonincreasing \( V(t) \). A very small local increase in a computed sample sequence can occur at command discontinuities or because of finite integration accuracy. Such increases should shrink when the step size and solver tolerances are tightened.

The adaptive closed loop is nonlinear because products such as \( \hat{\theta}_y y \), \( ye \), and \( re \) appear in the differential equations. The supplied Python, C++, and Java programs use a classical fourth-order Runge–Kutta step:

\[ \begin{aligned} \mathbf{k}_1&=f(t_k,\mathbf{x}_k),\\ \mathbf{k}_2&=f\!\left(t_k+\frac{h}{2}, \mathbf{x}_k+\frac{h}{2}\mathbf{k}_1\right),\\ \mathbf{k}_3&=f\!\left(t_k+\frac{h}{2}, \mathbf{x}_k+\frac{h}{2}\mathbf{k}_2\right),\\ \mathbf{k}_4&=f(t_k+h,\mathbf{x}_k+h\mathbf{k}_3),\\ \mathbf{x}_{k+1}&=\mathbf{x}_k+ \frac{h}{6}\left(\mathbf{k}_1+2\mathbf{k}_2+ 2\mathbf{k}_3+\mathbf{k}_4\right). \end{aligned} \]

The continuous-time Lyapunov proof does not automatically certify every discretization. A simulation should therefore pass the following checks:

1. Halving \( h \) should produce nearly the same output, error, and parameter trajectories.
2. The computed solution must remain bounded and finite; no state may become NaN or infinite.
3. The ideal parameters must satisfy the two matching equations numerically.
4. The cumulative index \( J_e(t) \) must be nondecreasing because \( \dot{J}_e=e^2\ge0 \).
5. A reconstructed Lyapunov sequence should be nearly nonincreasing, with numerical deviations reduced by tighter integration settings.

Excessively large adaptation gains can create fast parameter transients, which require a smaller integration step. Therefore adaptation gain and solver step size must be selected together in simulation.

This implementation uses NumPy for arrays and Matplotlib for plotting. The Runge–Kutta integrator is written from scratch so that the coupled adaptive equations are explicit.

Chapter5_Lesson5.py

"""Chapter 5, Lesson 5: first-order Lyapunov MRAC simulation.

Plant:
    y_dot = -a*y + b*u
Reference model:
    ym_dot = -am*ym + bm*r
Controller:
    u = theta_y*y + theta_r*r
Adaptive laws (known b > 0):
    theta_y_dot = -gamma_y*y*e
    theta_r_dot = -gamma_r*r*e
where e = y - ym.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    a: float = 1.2
    b: float = 0.8
    am: float = 2.0
    bm: float = 2.0
    gamma_y: float = 4.0
    gamma_r: float = 4.0
    t_final: float = 60.0
    dt: float = 0.005

    @property
    def theta_y_star(self) -> float:
        return (self.a - self.am) / self.b

    @property
    def theta_r_star(self) -> float:
        return self.bm / self.b


def reference(t: float) -> float:
    """Reference signal with constant and dynamically rich intervals."""
    if t < 10.0:
        return 1.0
    if t < 20.0:
        return -0.5
    return 0.8 * np.sin(1.1 * t) + 0.45 * np.sin(0.37 * t)


def rhs(t: float, state: np.ndarray, p: Parameters) -> np.ndarray:
    """Closed-loop differential equations.

    state = [y, ym, theta_y, theta_r, integral_error_squared]
    """
    y, ym, theta_y, theta_r, _ = state
    r = reference(t)
    error = y - ym
    control = theta_y * y + theta_r * r

    y_dot = -p.a * y + p.b * control
    ym_dot = -p.am * ym + p.bm * r
    theta_y_dot = -p.gamma_y * y * error
    theta_r_dot = -p.gamma_r * r * error
    performance_dot = error * error

    return np.array(
        [y_dot, ym_dot, theta_y_dot, theta_r_dot, performance_dot],
        dtype=float,
    )


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


def simulate(p: Parameters) -> tuple[np.ndarray, np.ndarray]:
    sample_count = int(round(p.t_final / p.dt)) + 1
    time = np.linspace(0.0, p.t_final, sample_count)
    state = np.zeros((sample_count, 5), dtype=float)
    state[0, :] = np.array([0.0, 0.0, 0.0, 0.5, 0.0])

    for k in range(sample_count - 1):
        state[k + 1, :] = rk4_step(time[k], state[k, :], p.dt, p)

    return time, state


def interval_rms(time: np.ndarray, error: np.ndarray, start: float, stop: float) -> float:
    mask = (time >= start) & (time <= stop)
    return float(np.sqrt(np.mean(error[mask] ** 2)))


def save_results(time: np.ndarray, state: np.ndarray, p: Parameters) -> None:
    y = state[:, 0]
    ym = state[:, 1]
    theta_y = state[:, 2]
    theta_r = state[:, 3]
    integral_error_squared = state[:, 4]
    r = np.array([reference(t) for t in time])
    error = y - ym
    control = theta_y * y + theta_r * r

    output_dir = Path(__file__).resolve().parent
    csv_path = output_dir / "Chapter5_Lesson5_results.csv"
    np.savetxt(
        csv_path,
        np.column_stack(
            [time, r, y, ym, error, control, theta_y, theta_r, integral_error_squared]
        ),
        delimiter=",",
        header="time,r,y,ym,error,u,theta_y,theta_r,integral_error_squared",
        comments="",
    )

    fig1, ax1 = plt.subplots(figsize=(10, 5))
    ax1.plot(time, r, label="reference input r")
    ax1.plot(time, ym, label="reference model ym")
    ax1.plot(time, y, "--", label="plant output y")
    ax1.set_xlabel("Time [s]")
    ax1.set_ylabel("Signal")
    ax1.set_title("First-order MRAC tracking")
    ax1.grid(True)
    ax1.legend()
    fig1.tight_layout()
    fig1.savefig(output_dir / "Chapter5_Lesson5_tracking.png", dpi=180)

    fig2, ax2 = plt.subplots(figsize=(10, 5))
    ax2.plot(time, theta_y, label="theta_y")
    ax2.plot(time, theta_r, label="theta_r")
    ax2.axhline(p.theta_y_star, linestyle="--", label="theta_y ideal")
    ax2.axhline(p.theta_r_star, linestyle="--", label="theta_r ideal")
    ax2.set_xlabel("Time [s]")
    ax2.set_ylabel("Adaptive parameter")
    ax2.set_title("Controller-parameter evolution")
    ax2.grid(True)
    ax2.legend()
    fig2.tight_layout()
    fig2.savefig(output_dir / "Chapter5_Lesson5_parameters.png", dpi=180)

    plt.show()


def main() -> None:
    p = Parameters()
    time, state = simulate(p)
    error = state[:, 0] - state[:, 1]

    print(f"Ideal theta_y = {p.theta_y_star:.6f}")
    print(f"Ideal theta_r = {p.theta_r_star:.6f}")
    print(f"Final theta_y = {state[-1, 2]:.6f}")
    print(f"Final theta_r = {state[-1, 3]:.6f}")
    print(f"Final tracking error = {error[-1]:.6e}")
    print(f"Integral of e^2 = {state[-1, 4]:.6f}")
    print(f"RMS error, 0-10 s = {interval_rms(time, error, 0.0, 10.0):.6f}")
    print(f"RMS error, 10-20 s = {interval_rms(time, error, 10.0, 20.0):.6f}")
    print(f"RMS error, 20-60 s = {interval_rms(time, error, 20.0, 60.0):.6f}")

    save_results(time, state, p)


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. It writes all simulated signals to a CSV file that can be plotted with Python, MATLAB, GNUplot, or a spreadsheet program.

Chapter5_Lesson5.cpp

// Chapter 5, Lesson 5: first-order Lyapunov MRAC simulation.
// Compile: g++ -std=c++17 -O2 Chapter5_Lesson5.cpp -o Chapter5_Lesson5

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

struct Parameters {
    double a = 1.2;
    double b = 0.8;
    double am = 2.0;
    double bm = 2.0;
    double gamma_y = 4.0;
    double gamma_r = 4.0;
    double t_final = 60.0;
    double dt = 0.005;

    [[nodiscard]] double thetaYStar() const { return (a - am) / b; }
    [[nodiscard]] double thetaRStar() const { return bm / b; }
};

using State = std::array<double, 5>;  // y, ym, theta_y, theta_r, integral(e^2)

double reference(const double t) {
    if (t < 10.0) {
        return 1.0;
    }
    if (t < 20.0) {
        return -0.5;
    }
    return 0.8 * std::sin(1.1 * t) + 0.45 * std::sin(0.37 * t);
}

State rhs(const double t, const State& x, const Parameters& p) {
    const double y = x[0];
    const double ym = x[1];
    const double theta_y = x[2];
    const double theta_r = x[3];
    const double r = reference(t);
    const double error = y - ym;
    const double control = theta_y * y + theta_r * r;

    return State{
        -p.a * y + p.b * control,
        -p.am * ym + p.bm * r,
        -p.gamma_y * y * error,
        -p.gamma_r * r * error,
        error * error,
    };
}

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

State rk4Step(const double t, const State& x, const double h, const Parameters& p) {
    const State k1 = rhs(t, x, p);
    const State k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h), p);
    const State k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h), p);
    const State k4 = rhs(t + h, addScaled(x, k3, h), p);

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

int main() {
    try {
        const Parameters p{};
        const int steps = static_cast<int>(std::llround(p.t_final / p.dt));
        State x{0.0, 0.0, 0.0, 0.5, 0.0};

        std::ofstream csv("Chapter5_Lesson5_results.csv");
        if (!csv) {
            throw std::runtime_error("Could not create Chapter5_Lesson5_results.csv");
        }
        csv << "time,r,y,ym,error,u,theta_y,theta_r,integral_error_squared\n";
        csv << std::setprecision(12);

        double squared_error_sum_0_10 = 0.0;
        double squared_error_sum_10_20 = 0.0;
        double squared_error_sum_20_60 = 0.0;
        int count_0_10 = 0;
        int count_10_20 = 0;
        int count_20_60 = 0;

        for (int k = 0; k <= steps; ++k) {
            const double t = k * p.dt;
            const double r = reference(t);
            const double error = x[0] - x[1];
            const double control = x[2] * x[0] + x[3] * r;

            csv << t << ',' << r << ',' << x[0] << ',' << x[1] << ',' << error << ','
                << control << ',' << x[2] << ',' << x[3] << ',' << x[4] << '\n';

            if (t <= 10.0) {
                squared_error_sum_0_10 += error * error;
                ++count_0_10;
            } else if (t <= 20.0) {
                squared_error_sum_10_20 += error * error;
                ++count_10_20;
            } else {
                squared_error_sum_20_60 += error * error;
                ++count_20_60;
            }

            if (k < steps) {
                x = rk4Step(t, x, p.dt, p);
            }
        }

        const double final_error = x[0] - x[1];
        std::cout << std::fixed << std::setprecision(6)
                  << "Ideal theta_y = " << p.thetaYStar() << '\n'
                  << "Ideal theta_r = " << p.thetaRStar() << '\n'
                  << "Final theta_y = " << x[2] << '\n'
                  << "Final theta_r = " << x[3] << '\n'
                  << "Final tracking error = " << final_error << '\n'
                  << "Integral of e^2 = " << x[4] << '\n'
                  << "RMS error, 0-10 s = "
                  << std::sqrt(squared_error_sum_0_10 / count_0_10) << '\n'
                  << "RMS error, 10-20 s = "
                  << std::sqrt(squared_error_sum_10_20 / count_10_20) << '\n'
                  << "RMS error, 20-60 s = "
                  << std::sqrt(squared_error_sum_20_60 / count_20_60) << '\n';
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }

    return 0;
}

The Java version uses the standard JDK numerical and file APIs. No external control library is required for this scalar example.

Chapter5_Lesson5.java

// Chapter 5, Lesson 5: first-order Lyapunov MRAC simulation.
// Compile: javac Chapter5_Lesson5.java
// Run:     java Chapter5_Lesson5

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

public final class Chapter5_Lesson5 {
    private static final class Parameters {
        final double a = 1.2;
        final double b = 0.8;
        final double am = 2.0;
        final double bm = 2.0;
        final double gammaY = 4.0;
        final double gammaR = 4.0;
        final double tFinal = 60.0;
        final double dt = 0.005;

        double thetaYStar() {
            return (a - am) / b;
        }

        double thetaRStar() {
            return bm / b;
        }
    }

    private static double reference(double t) {
        if (t < 10.0) {
            return 1.0;
        }
        if (t < 20.0) {
            return -0.5;
        }
        return 0.8 * Math.sin(1.1 * t) + 0.45 * Math.sin(0.37 * t);
    }

    private static double[] rhs(double t, double[] x, Parameters p) {
        double y = x[0];
        double ym = x[1];
        double thetaY = x[2];
        double thetaR = x[3];
        double r = reference(t);
        double error = y - ym;
        double control = thetaY * y + thetaR * r;

        return new double[] {
            -p.a * y + p.b * control,
            -p.am * ym + p.bm * r,
            -p.gammaY * y * error,
            -p.gammaR * r * error,
            error * error
        };
    }

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

    private static double[] rk4Step(double t, double[] x, double h, Parameters p) {
        double[] k1 = rhs(t, x, p);
        double[] k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h), p);
        double[] k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h), p);
        double[] k4 = rhs(t + h, addScaled(x, k3, h), p);

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

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        Parameters p = new Parameters();
        int steps = (int) Math.round(p.tFinal / p.dt);
        double[] x = {0.0, 0.0, 0.0, 0.5, 0.0};

        double sum0To10 = 0.0;
        double sum10To20 = 0.0;
        double sum20To60 = 0.0;
        int count0To10 = 0;
        int count10To20 = 0;
        int count20To60 = 0;

        Path output = Path.of("Chapter5_Lesson5_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write("time,r,y,ym,error,u,theta_y,theta_r,integral_error_squared\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * p.dt;
                double r = reference(t);
                double error = x[0] - x[1];
                double control = x[2] * x[0] + x[3] * r;

                writer.write(String.format(
                        Locale.ROOT,
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        t, r, x[0], x[1], error, control, x[2], x[3], x[4]));

                if (t <= 10.0) {
                    sum0To10 += error * error;
                    count0To10++;
                } else if (t <= 20.0) {
                    sum10To20 += error * error;
                    count10To20++;
                } else {
                    sum20To60 += error * error;
                    count20To60++;
                }

                if (k < steps) {
                    x = rk4Step(t, x, p.dt, p);
                }
            }
        } catch (IOException exception) {
            System.err.println("Could not write results: " + exception.getMessage());
            System.exit(1);
        }

        double finalError = x[0] - x[1];
        System.out.printf("Ideal theta_y = %.6f%n", p.thetaYStar());
        System.out.printf("Ideal theta_r = %.6f%n", p.thetaRStar());
        System.out.printf("Final theta_y = %.6f%n", x[2]);
        System.out.printf("Final theta_r = %.6f%n", x[3]);
        System.out.printf("Final tracking error = %.6e%n", finalError);
        System.out.printf("Integral of e^2 = %.6f%n", x[4]);
        System.out.printf("RMS error, 0-10 s = %.6f%n",
                Math.sqrt(sum0To10 / count0To10));
        System.out.printf("RMS error, 10-20 s = %.6f%n",
                Math.sqrt(sum10To20 / count10To20));
        System.out.printf("RMS error, 20-60 s = %.6f%n",
                Math.sqrt(sum20To60 / count20To60));
    }
}

The MATLAB script uses ode45, tables, and plotting functions. Setting CREATE_SIMULINK_MODEL = true constructs an equivalent Simulink model from standard blocks: Integrator, Gain, Product, Sum, From Workspace, Mux, and To Workspace.

Chapter5_Lesson5.m

%% Chapter5_Lesson5.m
% Chapter 5, Lesson 5: first-order Lyapunov MRAC simulation.
% The script uses ODE45 and can optionally construct an equivalent Simulink model.

clear; clc; close all;

p.a = 1.2;
p.b = 0.8;
p.am = 2.0;
p.bm = 2.0;
p.gamma_y = 4.0;
p.gamma_r = 4.0;
p.t_final = 60.0;
p.output_step = 0.005;

p.theta_y_star = (p.a - p.am) / p.b;
p.theta_r_star = p.bm / p.b;

% State: [y; ym; theta_y; theta_r; integral(e^2)]
x0 = [0; 0; 0; 0.5; 0];
t_eval = (0:p.output_step:p.t_final).';
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, x] = ode45(@(time, state) mracRhs(time, state, p), t_eval, x0, options);

y = x(:, 1);
ym = x(:, 2);
theta_y = x(:, 3);
theta_r = x(:, 4);
integral_error_squared = x(:, 5);
r = arrayfun(@referenceSignal, t);
e = y - ym;
u = theta_y .* y + theta_r .* r;

fprintf('Ideal theta_y = %.6f\n', p.theta_y_star);
fprintf('Ideal theta_r = %.6f\n', p.theta_r_star);
fprintf('Final theta_y = %.6f\n', theta_y(end));
fprintf('Final theta_r = %.6f\n', theta_r(end));
fprintf('Final tracking error = %.6e\n', e(end));
fprintf('Integral of e^2 = %.6f\n', integral_error_squared(end));
fprintf('RMS error, 0-10 s = %.6f\n', intervalRms(t, e, 0, 10));
fprintf('RMS error, 10-20 s = %.6f\n', intervalRms(t, e, 10, 20));
fprintf('RMS error, 20-60 s = %.6f\n', intervalRms(t, e, 20, 60));

results = table(t, r, y, ym, e, u, theta_y, theta_r, ...
    integral_error_squared);
writetable(results, 'Chapter5_Lesson5_results.csv');

figure('Name', 'Chapter 5 Lesson 5 - Tracking');
plot(t, r, 'LineWidth', 1.0); hold on;
plot(t, ym, 'LineWidth', 1.5);
plot(t, y, '--', 'LineWidth', 1.5);
grid on;
xlabel('Time [s]');
ylabel('Signal');
title('First-order MRAC tracking');
legend('reference input r', 'reference model y_m', 'plant output y', ...
    'Location', 'best');

figure('Name', 'Chapter 5 Lesson 5 - Parameters');
plot(t, theta_y, 'LineWidth', 1.5); hold on;
plot(t, theta_r, 'LineWidth', 1.5);
yline(p.theta_y_star, '--', 'theta_y ideal');
yline(p.theta_r_star, '--', 'theta_r ideal');
grid on;
xlabel('Time [s]');
ylabel('Adaptive parameter');
title('Controller-parameter evolution');
legend('theta_y', 'theta_r', 'Location', 'best');

% Set this flag to true to construct an equivalent block-diagram model.
CREATE_SIMULINK_MODEL = false;
if CREATE_SIMULINK_MODEL
    if license('test', 'Simulink')
        r_ts = timeseries(r, t); %#ok<NASGU>
        assignin('base', 'r_ts', r_ts);
        assignin('base', 'p', p);
        buildSimulinkModel();
        fprintf('Created Chapter5_Lesson5_Simulink.slx\n');
    else
        warning('Simulink is not licensed on this MATLAB installation.');
    end
end

function dx = mracRhs(t, x, p)
    y = x(1);
    ym = x(2);
    theta_y = x(3);
    theta_r = x(4);
    r = referenceSignal(t);
    e = y - ym;
    u = theta_y * y + theta_r * r;

    dx = [
        -p.a * y + p.b * u;
        -p.am * ym + p.bm * r;
        -p.gamma_y * y * e;
        -p.gamma_r * r * e;
        e^2
    ];
end

function r = referenceSignal(t)
    if t < 10
        r = 1.0;
    elseif t < 20
        r = -0.5;
    else
        r = 0.8 * sin(1.1 * t) + 0.45 * sin(0.37 * t);
    end
end

function value = intervalRms(t, e, start_time, stop_time)
    mask = (t >= start_time) & (t <= stop_time);
    value = sqrt(mean(e(mask).^2));
end

function buildSimulinkModel()
    model = 'Chapter5_Lesson5_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', 'p.t_final', 'Solver', 'ode45', ...
        'ReturnWorkspaceOutputs', 'on');

    add_block('simulink/Sources/From Workspace', [model '/Reference'], ...
        'VariableName', 'r_ts', 'Position', [35 230 135 260]);

    add_block('simulink/Continuous/Integrator', [model '/Plant y'], ...
        'InitialCondition', '0', 'Position', [650 75 680 105]);
    add_block('simulink/Continuous/Integrator', [model '/Model ym'], ...
        'InitialCondition', '0', 'Position', [650 300 680 330]);
    add_block('simulink/Continuous/Integrator', [model '/Theta y'], ...
        'InitialCondition', '0', 'Position', [650 430 680 460]);
    add_block('simulink/Continuous/Integrator', [model '/Theta r'], ...
        'InitialCondition', '0.5', 'Position', [650 535 680 565]);

    add_block('simulink/Math Operations/Product', [model '/theta_y times y'], ...
        'Position', [230 70 270 105]);
    add_block('simulink/Math Operations/Product', [model '/theta_r times r'], ...
        'Position', [230 145 270 180]);
    add_block('simulink/Math Operations/Sum', [model '/Control sum'], ...
        'Inputs', '++', 'Position', [320 100 350 160]);
    add_block('simulink/Math Operations/Gain', [model '/b'], ...
        'Gain', 'p.b', 'Position', [405 115 455 145]);
    add_block('simulink/Math Operations/Gain', [model '/minus a'], ...
        'Gain', '-p.a', 'Position', [405 50 455 80]);
    add_block('simulink/Math Operations/Sum', [model '/Plant derivative'], ...
        'Inputs', '++', 'Position', [535 65 565 125]);

    add_block('simulink/Math Operations/Gain', [model '/bm'], ...
        'Gain', 'p.bm', 'Position', [405 260 455 290]);
    add_block('simulink/Math Operations/Gain', [model '/minus am'], ...
        'Gain', '-p.am', 'Position', [405 325 455 355]);
    add_block('simulink/Math Operations/Sum', [model '/Model derivative'], ...
        'Inputs', '++', 'Position', [535 280 565 340]);

    add_block('simulink/Math Operations/Sum', [model '/Tracking error'], ...
        'Inputs', '+-', 'Position', [740 175 770 235]);
    add_block('simulink/Math Operations/Product', [model '/y times e'], ...
        'Position', [825 395 865 430]);
    add_block('simulink/Math Operations/Product', [model '/r times e'], ...
        'Position', [825 500 865 535]);
    add_block('simulink/Math Operations/Gain', [model '/minus gamma y'], ...
        'Gain', '-p.gamma_y', 'Position', [915 395 985 425]);
    add_block('simulink/Math Operations/Gain', [model '/minus gamma r'], ...
        'Gain', '-p.gamma_r', 'Position', [915 500 985 530]);

    add_block('simulink/Signal Routing/Mux', [model '/Output mux'], ...
        'Inputs', '6', 'Position', [1060 120 1065 335]);
    add_block('simulink/Sinks/To Workspace', [model '/Simulation data'], ...
        'VariableName', 'simout', 'SaveFormat', 'Timeseries', ...
        'Position', [1135 205 1235 235]);

    add_line(model, 'Theta y/1', 'theta_y times y/1', 'autorouting', 'on');
    add_line(model, 'Plant y/1', 'theta_y times y/2', 'autorouting', 'on');
    add_line(model, 'Theta r/1', 'theta_r times r/1', 'autorouting', 'on');
    add_line(model, 'Reference/1', 'theta_r times r/2', 'autorouting', 'on');
    add_line(model, 'theta_y times y/1', 'Control sum/1', 'autorouting', 'on');
    add_line(model, 'theta_r times r/1', 'Control sum/2', 'autorouting', 'on');
    add_line(model, 'Control sum/1', 'b/1', 'autorouting', 'on');
    add_line(model, 'Plant y/1', 'minus a/1', 'autorouting', 'on');
    add_line(model, 'minus a/1', 'Plant derivative/1', 'autorouting', 'on');
    add_line(model, 'b/1', 'Plant derivative/2', 'autorouting', 'on');
    add_line(model, 'Plant derivative/1', 'Plant y/1', 'autorouting', 'on');

    add_line(model, 'Reference/1', 'bm/1', 'autorouting', 'on');
    add_line(model, 'Model ym/1', 'minus am/1', 'autorouting', 'on');
    add_line(model, 'bm/1', 'Model derivative/1', 'autorouting', 'on');
    add_line(model, 'minus am/1', 'Model derivative/2', 'autorouting', 'on');
    add_line(model, 'Model derivative/1', 'Model ym/1', 'autorouting', 'on');

    add_line(model, 'Plant y/1', 'Tracking error/1', 'autorouting', 'on');
    add_line(model, 'Model ym/1', 'Tracking error/2', 'autorouting', 'on');
    add_line(model, 'Plant y/1', 'y times e/1', 'autorouting', 'on');
    add_line(model, 'Tracking error/1', 'y times e/2', 'autorouting', 'on');
    add_line(model, 'Reference/1', 'r times e/1', 'autorouting', 'on');
    add_line(model, 'Tracking error/1', 'r times e/2', 'autorouting', 'on');
    add_line(model, 'y times e/1', 'minus gamma y/1', 'autorouting', 'on');
    add_line(model, 'r times e/1', 'minus gamma r/1', 'autorouting', 'on');
    add_line(model, 'minus gamma y/1', 'Theta y/1', 'autorouting', 'on');
    add_line(model, 'minus gamma r/1', 'Theta r/1', 'autorouting', 'on');

    add_line(model, 'Reference/1', 'Output mux/1', 'autorouting', 'on');
    add_line(model, 'Plant y/1', 'Output mux/2', 'autorouting', 'on');
    add_line(model, 'Model ym/1', 'Output mux/3', 'autorouting', 'on');
    add_line(model, 'Tracking error/1', 'Output mux/4', 'autorouting', 'on');
    add_line(model, 'Theta y/1', 'Output mux/5', 'autorouting', 'on');
    add_line(model, 'Theta r/1', 'Output mux/6', 'autorouting', 'on');
    add_line(model, 'Output mux/1', 'Simulation data/1', 'autorouting', 'on');

    save_system(model, [model '.slx']);
end

The notebook uses NDSolveValue for the coupled nonlinear differential equations, Plot for trajectory visualization, and Export for CSV generation.

Chapter5_Lesson5.nb

Notebook[{
  Cell["Chapter 5, Lesson 5: First-Order Lyapunov MRAC", "Title"],
  Cell["The notebook simulates the plant, reference model, and two adaptive controller parameters. The reference signal begins with constant commands and then becomes a two-frequency signal so that parameter evolution can be interpreted.", "Text"],
  Cell[BoxData["ClearAll[\"Global`*\"];

a = 1.2; b = 0.8; am = 2.0; bm = 2.0;
gammaY = 4.0; gammaR = 4.0; tFinal = 60.0; dt = 0.005;

thetaYStar = (a - am)/b; thetaRStar = bm/b;

r[t_?NumericQ] := Piecewise[{ {1.0, t < 10.0}, {-0.5, t < 20.0} },
  0.8 Sin[1.1 t] + 0.45 Sin[0.37 t]];

solution = NDSolveValue[{
   y'[t] == -a y[t] + b (thetaY[t] y[t] + thetaR[t] r[t]),
   ym'[t] == -am ym[t] + bm r[t],
   thetaY'[t] == -gammaY y[t] (y[t] - ym[t]),
   thetaR'[t] == -gammaR r[t] (y[t] - ym[t]),
   performance'[t] == (y[t] - ym[t])^2,
   y[0] == 0, ym[0] == 0, thetaY[0] == 0, thetaR[0] == 0.5,
   performance[0] == 0
   }, {y, ym, thetaY, thetaR, performance}, {t, 0, tFinal},
  Method -> {\"TimeIntegration\" -> {\"ExplicitRungeKutta\", \"DifferenceOrder\" -> 4} },
  AccuracyGoal -> 10, PrecisionGoal -> 10];

ySol = solution[[1]]; ymSol = solution[[2]];
thetaYSol = solution[[3]]; thetaRSol = solution[[4]];
performanceSol = solution[[5]];

Print[\"Ideal theta_y = \", N[thetaYStar, 8]];
Print[\"Ideal theta_r = \", N[thetaRStar, 8]];
Print[\"Final theta_y = \", N[thetaYSol[tFinal], 8]];
Print[\"Final theta_r = \", N[thetaRSol[tFinal], 8]];
Print[\"Final tracking error = \", N[ySol[tFinal] - ymSol[tFinal], 8]];
Print[\"Integral of e^2 = \", N[performanceSol[tFinal], 8]];

trackingPlot = Plot[{r[t], ymSol[t], ySol[t]}, {t, 0, tFinal},
  PlotLegends -> {\"reference input r\", \"reference model ym\", \"plant output y\"},
  PlotLabel -> \"First-order MRAC tracking\", AxesLabel -> {\"Time [s]\", \"Signal\"},
  PlotRange -> All, ImageSize -> Large];

parameterPlot = Plot[{thetaYSol[t], thetaRSol[t], thetaYStar, thetaRStar},
  {t, 0, tFinal},
  PlotLegends -> {\"theta_y\", \"theta_r\", \"theta_y ideal\", \"theta_r ideal\"},
  PlotLabel -> \"Controller-parameter evolution\",
  AxesLabel -> {\"Time [s]\", \"Adaptive parameter\"}, PlotRange -> All,
  ImageSize -> Large];

Print[trackingPlot]; Print[parameterPlot];

sampleTimes = Range[0, tFinal, dt];
data = Table[{tt, r[tt], ySol[tt], ymSol[tt], ySol[tt] - ymSol[tt],
    thetaYSol[tt] ySol[tt] + thetaRSol[tt] r[tt], thetaYSol[tt],
    thetaRSol[tt], performanceSol[tt]}, {tt, sampleTimes}];
Export[\"Chapter5_Lesson5_results.csv\",
  Prepend[data, {\"time\", \"r\", \"y\", \"ym\", \"error\", \"u\",
    \"theta_y\", \"theta_r\", \"integral_error_squared\"}]];"], "Input"]
},
WindowSize -> {1200, 800}, StyleDefinitions -> "Default.nb"
]        
      

flowchart TD
  A["Inspect tracking error e"] --> B{"Is e appreciably \nnonzero?"}
  B -->|"No"| C["Adaptation nearly stops; \nparameter values may freeze"]
  B -->|"Yes"| D["Inspect y and r magnitudes"]
  D --> E{"Is the associated \nregressor near zero?"}
  E -->|"Yes"| F["That parameter \nchanges slowly even \nthough error exists"]
  E -->|"No"| G["Use sign of \nregressor times error \nto predict motion"]
  G --> H["Compare tracking, \nparameter combination, \nand ideal values"]
  H --> I{"Only constant \nsteady-state data?"}
  I -->|"Yes"| J["Several parameter \npairs may produce \nthe same tracking"]
  I -->|"No"| K["Changing dynamic \ndata can separate \nparameter effects"]
  C --> L["Check Lyapunov trend \nand numerical step size"]
  J --> L
  K --> L
        

The workflow prevents a frequent interpretation error: declaring a simulation unsuccessful because the parameters do not immediately equal their ideal values, even though the proven objective—stable tracking—is being achieved. Conversely, parameter motion alone is not evidence of good control; the output, error, control effort, and boundedness must all be checked.

Problem 1 (Matching Parameters): For \( \dot{y}=-1.2y+0.8u \), \( \dot{y}_m=-2y_m+2r \), and \( u=\theta_y y+\theta_r r \), compute the ideal controller parameters.

Solution: Matching the coefficients gives

\[ -1.2+0.8\theta_y^*=-2, \qquad 0.8\theta_r^*=2. \]

Therefore

\[ \theta_y^*=\frac{1.2-2}{0.8}=-1, \qquad \theta_r^*=\frac{2}{0.8}=2.5. \]

Problem 2 (Initial Adaptive-Law Direction): Suppose \( r(0)=1 \), \( y(0)=y_m(0)=0 \), and the initial parameters are those used in the lesson. Find \( \dot{\hat{\theta} }_y(0) \) and \( \dot{\hat{\theta} }_r(0) \). Why can the parameters nevertheless start changing immediately after \( t=0 \)?

Solution: Since \( e(0)=0 \),

\[ \dot{\hat{\theta} }_y(0)=-4y(0)e(0)=0, \qquad \dot{\hat{\theta} }_r(0)=-4r(0)e(0)=0. \]

However, the initial model derivative is \( \dot{y}_m(0)=2 \), whereas the initial plant control is \( u(0)=0.5 \), so \( \dot{y}(0)=0.8(0.5)=0.4 \). Hence

\[ \dot{e}(0)=0.4-2=-1.6. \]

The error becomes negative immediately, causing the parameter update to begin for \( t > 0 \).

Problem 3 (Lyapunov Cancellation): Starting from \( \dot{e}=-a_m e+b\widetilde{\theta}_y y+ b\widetilde{\theta}_r r \), show that the selected update laws cancel the parameter cross terms in \( \dot{V} \).

Solution: Differentiate

\[ V=\frac{1}{2}e^2+ \frac{b}{2\gamma_y}\widetilde{\theta}_y^2+ \frac{b}{2\gamma_r}\widetilde{\theta}_r^2. \]

Then

\[ \dot{V}=-a_m e^2+b\widetilde{\theta}_yye+ b\widetilde{\theta}_rre+ \frac{b}{\gamma_y}\widetilde{\theta}_y \dot{\hat{\theta} }_y+ \frac{b}{\gamma_r}\widetilde{\theta}_r \dot{\hat{\theta} }_r. \]

Substituting \( \dot{\hat{\theta} }_y=-\gamma_yye \) and \( \dot{\hat{\theta} }_r=-\gamma_rre \) gives equal and opposite cross terms. Thus

\[ \dot{V}=-a_m e^2\le0. \]

Problem 4 (Constant-Command Parameter Manifold): For a general reference-model steady-state ratio \( \rho=b_m/a_m \), derive the set of controller parameters that can maintain exact steady-state tracking of a nonzero constant command.

Solution: Exact steady-state tracking implies

\[ y=y_m=\rho r, \qquad \dot{y}=0. \]

Substitute these expressions into the plant:

\[ 0=-a\rho r+b(\hat{\theta}_y\rho r+ \hat{\theta}_r r). \]

For \( r\ne0 \),

\[ \rho\hat{\theta}_y+\hat{\theta}_r= \frac{a\rho}{b}. \]

This affine line contains infinitely many parameter pairs. Hence constant steady-state tracking alone cannot distinguish the two parameters.

Problem 5 (Effect of Adaptation Gain): At one instant, suppose \( y=0.8 \), \( r=1 \), and \( e=-0.1 \). Compute the parameter rates for \( \gamma_y=\gamma_r=4 \). What happens if both gains are doubled while the instantaneous signals are held fixed?

Solution:

\[ \dot{\hat{\theta} }_y=-4(0.8)(-0.1)=0.32, \qquad \dot{\hat{\theta} }_r=-4(1)(-0.1)=0.4. \]

Doubling both gains doubles the instantaneous rates to \( 0.64 \) and \( 0.8 \). This local calculation does not mean the full transient is exactly twice as fast, because the faster parameter change also alters \( y \), \( e \), and the future update direction.

Problem 6 (One Explicit Euler Step): At time \( t_k \), let \( y=0.8 \), \( y_m=0.9 \), \( r=1 \), \( \hat{\theta}_y=-0.2 \), \( \hat{\theta}_r=1.6 \), and \( h=0.01 \). Compute one explicit Euler update of all four principal states.

Solution: First,

\[ e=0.8-0.9=-0.1, \qquad u=(-0.2)(0.8)+(1.6)(1)=1.44. \]

The derivatives are

\[ \begin{aligned} \dot{y}&=-1.2(0.8)+0.8(1.44)=0.192,\\ \dot{y}_m&=-2(0.9)+2(1)=0.2,\\ \dot{\hat{\theta} }_y&=-4(0.8)(-0.1)=0.32,\\ \dot{\hat{\theta} }_r&=-4(1)(-0.1)=0.4. \end{aligned} \]

Therefore

\[ \begin{aligned} y_{k+1}&=0.8+0.01(0.192)=0.80192,\\ y_{m,k+1}&=0.9+0.01(0.2)=0.902,\\ \hat{\theta}_{y,k+1}&=-0.2+0.01(0.32)=-0.1968,\\ \hat{\theta}_{r,k+1}&=1.6+0.01(0.4)=1.604. \end{aligned} \]

Euler's method is useful for hand calculations, but the supplied programs use fourth-order Runge–Kutta for better accuracy at a practical step size.

A first-order plant, reference model, adaptive controller, and two parameter update laws were simulated as one coupled nonlinear system. The Lyapunov derivative \( \dot{V}=-a_m e^2 \) explains boundedness and asymptotic tracking. Parameter evolution is governed by the product of tracking error and the associated regressor: when the error becomes small, adaptation naturally slows or stops.

Constant steady-state tracking constrains only a combination of the two controller parameters, so accurate output tracking can occur before the individual estimates reach their ideal values. A dynamically changing command provides more information and, in this numerical example, moves the estimates much closer to the ideal pair. Finally, solver convergence, Lyapunov consistency, boundedness, and cumulative error should be checked whenever continuous-time adaptive laws are implemented numerically.