After completing this lesson, students should be able to:

• derive the matching conditions for a first-order direct MRAC;
• obtain the tracking-error model in terms of parameter errors;
• construct a composite Lyapunov function for the tracking and parameter errors;
• derive adaptive laws by cancelling cross terms in the Lyapunov derivative;
• prove boundedness and asymptotic tracking using Lyapunov theory and Barbalat’s lemma;
• implement the same continuous-time design using numerical integration and Simulink.

The lesson uses the Lyapunov-stability, combined state–parameter Lyapunov functions, and Barbalat arguments introduced in Chapter 3. It also uses the first-order plant, reference model, matching conditions, and error model developed in Chapter 5, Lessons 1–3.

Consider the uncertain scalar plant

\[ \dot{x}(t)=a\,x(t)+b\,u(t), \qquad b\neq 0, \]

where \( a \) and the magnitude of \( b \) are unknown. The sign \( s_b=\operatorname{sgn}(b) \) is assumed known. The desired dynamics are specified by the stable reference model

\[ \dot{x}_m(t)=a_m x_m(t)+b_m r(t), \qquad a_m < 0, \]

where \( r(t) \) is a bounded command. The direct adaptive controller is chosen as

\[ u(t)=\hat{k}_x(t)x(t)+\hat{k}_r(t)r(t). \]

The control objective is to make the tracking error \( e(t)=x(t)-x_m(t) \) converge to zero while keeping all internal signals bounded.

The known-sign assumption is essential for this elementary design. The unknown-control-direction problem is postponed to Chapter 16.

First imagine that the plant parameters were known and that constant gains \( k_x^\ast \) and \( k_r^\ast \) could be selected. Substitution of \( u=k_x^\ast x+k_r^\ast r \) into the plant gives

\[ \dot{x} = \left(a+b k_x^\ast\right)x + b k_r^\ast r. \]

Exact model matching requires the plant closed-loop coefficients to equal those of the reference model:

\[ a+b k_x^\ast=a_m, \qquad b k_r^\ast=b_m. \]

Therefore, the ideal gains are

\[ k_x^\ast=\frac{a_m-a}{b}, \qquad k_r^\ast=\frac{b_m}{b}. \]

These constants are used only in the analysis. The controller does not know them because \( a \) and \( b \) are unknown.

Define the parameter-estimation errors

\[ \tilde{k}_x=\hat{k}_x-k_x^\ast, \qquad \tilde{k}_r=\hat{k}_r-k_r^\ast. \]

Since the ideal parameters are constant, their derivatives vanish:

\[ \dot{\tilde{k}}_x=\dot{\hat{k}}_x, \qquad \dot{\tilde{k}}_r=\dot{\hat{k}}_r. \]

Starting from \( e=x-x_m \), subtract the reference-model equation from the plant equation:

\[ \begin{aligned} \dot{e} &= a x+b\left(\hat{k}_x x+\hat{k}_r r\right) -a_m x_m-b_m r \\\\ &= a x+b\left(k_x^\ast+\tilde{k}_x\right)x +b\left(k_r^\ast+\tilde{k}_r\right)r -a_m x_m-b_m r. \end{aligned} \]

Apply the matching equalities:

\[ \begin{aligned} \dot{e} &= a_m x+b\tilde{k}_x x+b\tilde{k}_r r-a_m x_m \\\\ &= a_m e+b\tilde{k}_x x+b\tilde{k}_r r. \end{aligned} \]

The terms \( b\tilde{k}_x x \) and \( b\tilde{k}_r r \) have indefinite signs. A Lyapunov design must remove their contribution from \( \dot{V} \).

The state-only candidate \( V_e=e^2/2 \) gives

\[ \dot{V}_e = a_m e^2 + b e\tilde{k}_x x + b e\tilde{k}_r r. \]

Although \( a_m e^2 \) is negative for nonzero error, the two parameter-error terms can be positive or negative. Therefore, \( V_e \) alone cannot establish stability.

Introduce positive adaptation gains \( \gamma_x > 0 \) and \( \gamma_r > 0 \), and choose the composite candidate

\[ V = \frac{1}{2}e^2 + \frac{|b|}{2\gamma_x}\tilde{k}_x^2 + \frac{|b|}{2\gamma_r}\tilde{k}_r^2. \]

For every nonzero \( b \) and positive adaptation gains, \( V \) is positive definite in \( (e,\tilde{k}_x,\tilde{k}_r) \). The unknown magnitude \( |b| \) appears only in the proof; it is not needed by the implemented update laws.

Differentiate the composite Lyapunov function:

\[ \dot{V} = e\dot{e} + \frac{|b|}{\gamma_x}\tilde{k}_x\dot{\hat{k}}_x + \frac{|b|}{\gamma_r}\tilde{k}_r\dot{\hat{k}}_r. \]

Substitute the error dynamics:

\[ \begin{aligned} \dot{V} &= a_m e^2 + b e\tilde{k}_x x + b e\tilde{k}_r r \\\\ &\quad+ \frac{|b|}{\gamma_x}\tilde{k}_x\dot{\hat{k}}_x + \frac{|b|}{\gamma_r}\tilde{k}_r\dot{\hat{k}}_r. \end{aligned} \]

Since \( b=s_b|b| \), collect the terms multiplying each parameter error:

\[ \dot{V} = a_m e^2 + |b|\tilde{k}_x \left( s_b e x+\frac{1}{\gamma_x}\dot{\hat{k}}_x \right) + |b|\tilde{k}_r \left( s_b e r+\frac{1}{\gamma_r}\dot{\hat{k}}_r \right). \]

Select the update laws so that both bracketed quantities are identically zero:

\[ \boxed{ \dot{\hat{k}}_x=-\gamma_x s_b e x }, \qquad \boxed{ \dot{\hat{k}}_r=-\gamma_r s_b e r }. \]

The Lyapunov derivative then reduces exactly to

\[ \boxed{ \dot{V}=a_m e^2=-|a_m|e^2 } \le 0. \]

For \( b > 0 \), use \( s_b=1 \). For \( b < 0 \), use \( s_b=-1 \), which reverses both update directions.

flowchart TD
  A["Plant and reference model"] --> B["Choose u = kx_hat*x + kr_hat*r"]
  B --> C["Derive e_dot = am*e + b*kx_tilde*x + b*kr_tilde*r"]
  C --> D["Choose V = error energy + parameter-error energy"]
  D --> E["Differentiate V"]
  E --> F["Select kx_hat_dot and kr_hat_dot to cancel cross terms"]
  F --> G["Obtain V_dot = am*e^2"]
  G --> H["Conclude boundedness and asymptotic tracking"]
        

Theorem. Suppose that:

1. the plant satisfies \( \dot{x}=a x+b u \) with constant unknown \( a \) and nonzero constant \( b \);
2. the sign \( s_b=\operatorname{sgn}(b) \) is known;
3. the reference model has \( a_m < 0 \);
4. the command \( r(t) \) is bounded and piecewise continuous;
5. the controller and adaptive laws are \( u=\hat{k}_x x+\hat{k}_r r \), \( \dot{\hat{k}}_x=-\gamma_x s_b e x \), and \( \dot{\hat{k}}_r=-\gamma_r s_b e r \), with positive adaptation gains.

Then all closed-loop signals are bounded, and \( \lim_{t→\infty}e(t)=0 \).

Proof, Part 1 — Non-increasing Lyapunov function.

From the previous section,

\[ \dot{V}=-|a_m|e^2\le 0. \]

Consequently,

\[ 0\le V(t)\le V(0), \qquad t\ge 0. \]

Because \( V \) is positive definite in \( e,\tilde{k}_x,\tilde{k}_r \), the tracking error and both parameter errors are bounded. Since the ideal parameters are finite constants, \( \hat{k}_x \) and \( \hat{k}_r \) are bounded.

Proof, Part 2 — Bounded plant and controller signals.

The stable reference model driven by bounded \( r(t) \) has bounded \( x_m(t) \). Since \( x=e+x_m \), the plant state is bounded. The control signal \( u=\hat{k}_x x+\hat{k}_r r \) is therefore bounded. The adaptive-law right-hand sides are products of bounded signals, so \( \dot{\hat{k}}_x \) and \( \dot{\hat{k}}_r \) are bounded.

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

Integrating the Lyapunov derivative gives

\[ V(t)-V(0) = -|a_m|\int_0^t e^2(\tau)\,d\tau. \]

Since \( V(t)\ge 0 \),

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

Thus \( e\in L_2 \).

Proof, Part 4 — Uniform continuity and Barbalat’s lemma.

The error derivative

\[ \dot{e} = a_m e+b\tilde{k}_x x+b\tilde{k}_r r \]

is bounded because every factor on its right-hand side is bounded. Therefore, \( e(t) \) is uniformly continuous. Since \( e\in L_2 \) and \( e \) is uniformly continuous, Barbalat’s lemma yields

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

This completes the proof.

flowchart TD
  R["Command r"] --> RM["Reference model: \nxm_dot = am*xm + bm*r"]
  R --> C["Controller: u = kx_hat*x + kr_hat*r"]
  C --> P["Unknown plant: x_dot = a*x + b*u"]
  P --> X["Plant state x"]
  RM --> XM["Model state xm"]
  X --> E["Tracking error e = x - xm"]
  XM --> E
  E --> AX["Update kx_hat using \n-gamma_x*sign_b*e*x"]
  X --> AX
  E --> AR["Update kr_hat using \n-gamma_r*sign_b*e*r"]
  R --> AR
  AX --> C
  AR --> C
        

The architecture is direct MRAC because the adjustable quantities are the controller gains themselves. No separate online estimate of \( a \) or \( b \) is formed.

9.1 The update law is an error–regressor product

The controller can be written in vector form as

\[ u=\hat{\boldsymbol{k}}^\mathsf{T}\boldsymbol{\phi}, \qquad \hat{\boldsymbol{k}} = \begin{bmatrix}\hat{k}_x\\\\\hat{k}_r\end{bmatrix}, \qquad \boldsymbol{\phi} = \begin{bmatrix}x\\\\r\end{bmatrix}. \]

The two scalar update laws become

\[ \dot{\hat{\boldsymbol{k}}} = -s_b\boldsymbol{\Gamma}\boldsymbol{\phi}e, \qquad \boldsymbol{\Gamma} = \begin{bmatrix} \gamma_x&0\\\\ 0&\gamma_r \end{bmatrix} \succ 0. \]

Thus the tracking error determines the direction of adaptation, while \( x \) and \( r \) distribute that correction among the feedback and feedforward gains.

9.2 Difference from the MIT rule

The MIT rule differentiates a cost such as \( J=e^2/2 \) with respect to adjustable parameters and requires sensitivity derivatives. Its stability is not automatic. The present design instead chooses a Lyapunov function first and then defines the adaptive laws to make \( \dot{V}\le 0 \). For this first-order matched plant, the resulting laws are simple error–regressor products, but their justification is a closed-loop stability proof rather than instantaneous gradient descent.

9.3 What the proof does and does not guarantee

• Guaranteed: bounded tracking error, bounded adaptive gains, bounded plant state, bounded control input, and asymptotic tracking.
• Not guaranteed here: convergence of each adaptive gain to its ideal value.

Tracking can become exact even when several parameter combinations produce the required behaviour along the particular trajectory. Stronger signal conditions for parameter convergence are studied later in the course.

1. Matching must be possible. For the scalar plant with nonzero \( b \), constant ideal gains exist. More complicated systems may violate the matching conditions.
2. The control direction must be known. Using the wrong value of \( s_b \) changes cancellation into reinforcement and invalidates the proof.
3. The reference model must be stable. The condition \( a_m < 0 \) is responsible for the negative term in \( \dot{V} \).
4. The elementary model excludes disturbances and unmodelled dynamics. Robust modifications are introduced in Chapter 8.
5. Actuator saturation is ignored. Saturation breaks the assumed plant equation because the commanded and applied control inputs differ.
6. Large adaptation gains are not universally better. They can accelerate parameter motion but increase sensitivity to noise, discretisation, and unmodelled high-frequency dynamics.

The simulations use the augmented state

\[ \boldsymbol{z} = \begin{bmatrix} x&x_m&\hat{k}_x&\hat{k}_r \end{bmatrix}^\mathsf{T}, \]

with differential equation

\[ \dot{\boldsymbol{z}} = \begin{bmatrix} a x+b(\hat{k}_x x+\hat{k}_r r)\\\\ a_m x_m+b_m r\\\\ -\gamma_x s_b(x-x_m)x\\\\ -\gamma_r s_b(x-x_m)r \end{bmatrix}. \]

A fourth-order Runge–Kutta method integrates all four states together. This is important: updating the plant and adaptive gains with inconsistent time levels can introduce artificial energy growth. The demonstration uses

\[ a=-0.4,\quad b=1.5,\quad a_m=-1.2,\quad b_m=1.2,\quad \gamma_x=\gamma_r=4. \]

The corresponding ideal gains are

\[ k_x^\ast=-0.533\overline{3}, \qquad k_r^\ast=0.8. \]

Python uses NumPy and Matplotlib. The C++ and Java versions use only their standard libraries and write CSV files. MATLAB uses base numerical and plotting functions. The Simulink script constructs the block model programmatically. Mathematica uses core Wolfram Language list operations and plotting functions.

Chapter5_Lesson4.py

"""
Chapter5_Lesson4.py
Lyapunov-based MRAC for a first-order plant.

Plant:
    x_dot = a*x + b*u
Reference model:
    xm_dot = am*xm + bm*r
Controller:
    u = kx_hat*x + kr_hat*r
Adaptive laws:
    kx_hat_dot = -gamma_x*sign(b)*e*x
    kr_hat_dot = -gamma_r*sign(b)*e*r
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    a: float = -0.4
    b: float = 1.5
    am: float = -1.2
    bm: float = 1.2
    gamma_x: float = 4.0
    gamma_r: float = 4.0
    dt: float = 0.002
    t_final: float = 30.0

    def validate(self) -> None:
        if self.b == 0.0:
            raise ValueError("The plant input gain b must be nonzero.")
        if self.am >= 0.0:
            raise ValueError("The reference-model pole am must be negative.")
        if self.gamma_x <= 0.0 or self.gamma_r <= 0.0:
            raise ValueError("Adaptation gains must be positive.")
        if self.dt <= 0.0 or self.t_final <= 0.0:
            raise ValueError("dt and t_final must be positive.")


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


def rhs(t: float, z: np.ndarray, p: Parameters) -> tuple[np.ndarray, float, float]:
    """Return z_dot, control input, and tracking error."""
    x, xm, kx_hat, kr_hat = z
    r = reference(t)
    e = x - xm
    u = kx_hat * x + kr_hat * r
    sign_b = math.copysign(1.0, p.b)

    dz = np.array(
        [
            p.a * x + p.b * u,
            p.am * xm + p.bm * r,
            -p.gamma_x * sign_b * e * x,
            -p.gamma_r * sign_b * e * r,
        ],
        dtype=float,
    )
    return dz, u, e


def rk4_step(t: float, z: np.ndarray, h: float, p: Parameters) -> np.ndarray:
    """One classical fourth-order Runge-Kutta step."""
    k1 = rhs(t, z, p)[0]
    k2 = rhs(t + 0.5 * h, z + 0.5 * h * k1, p)[0]
    k3 = rhs(t + 0.5 * h, z + 0.5 * h * k2, p)[0]
    k4 = rhs(t + h, z + h * k3, p)[0]
    return z + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(p: Parameters) -> dict[str, np.ndarray]:
    p.validate()
    n = int(round(p.t_final / p.dt)) + 1
    time = np.linspace(0.0, p.t_final, n)
    state = np.zeros((n, 4), dtype=float)
    control = np.zeros(n, dtype=float)
    command = np.zeros(n, dtype=float)
    error = np.zeros(n, dtype=float)
    lyapunov = np.zeros(n, dtype=float)

    # Nonzero initial mismatch makes the transient visible.
    state[0] = np.array([0.5, 0.0, 0.0, 0.0])

    kx_star = (p.am - p.a) / p.b
    kr_star = p.bm / p.b

    for i, t in enumerate(time):
        _, control[i], error[i] = rhs(t, state[i], p)
        command[i] = reference(t)
        kx_tilde = state[i, 2] - kx_star
        kr_tilde = state[i, 3] - kr_star
        lyapunov[i] = (
            0.5 * error[i] ** 2
            + abs(p.b) * kx_tilde**2 / (2.0 * p.gamma_x)
            + abs(p.b) * kr_tilde**2 / (2.0 * p.gamma_r)
        )
        if i + 1 < n:
            state[i + 1] = rk4_step(t, state[i], p.dt, p)

    return {
        "time": time,
        "x": state[:, 0],
        "xm": state[:, 1],
        "kx_hat": state[:, 2],
        "kr_hat": state[:, 3],
        "u": control,
        "r": command,
        "e": error,
        "V": lyapunov,
        "kx_star": np.array([kx_star]),
        "kr_star": np.array([kr_star]),
    }


def main() -> None:
    p = Parameters()
    data = simulate(p)

    print(f"Ideal kx*: {data['kx_star'][0]: .6f}")
    print(f"Ideal kr*: {data['kr_star'][0]: .6f}")
    print(f"Final tracking error: {data['e'][-1]: .6e}")
    print(f"Maximum V increase from numerical integration: {np.max(np.diff(data['V'])): .6e}")

    fig, axes = plt.subplots(4, 1, figsize=(10, 11), sharex=True)
    axes[0].plot(data["time"], data["x"], label="plant x")
    axes[0].plot(data["time"], data["xm"], "--", label="model xm")
    axes[0].plot(data["time"], data["r"], ":", label="command r")
    axes[0].set_ylabel("signals")
    axes[0].legend()
    axes[0].grid(True)

    axes[1].plot(data["time"], data["e"])
    axes[1].set_ylabel("error e")
    axes[1].grid(True)

    axes[2].plot(data["time"], data["kx_hat"], label="kx_hat")
    axes[2].plot(data["time"], data["kr_hat"], label="kr_hat")
    axes[2].axhline(data["kx_star"][0], linestyle="--", label="kx*")
    axes[2].axhline(data["kr_star"][0], linestyle=":", label="kr*")
    axes[2].set_ylabel("parameters")
    axes[2].legend()
    axes[2].grid(True)

    axes[3].plot(data["time"], data["V"], label="V")
    axes[3].plot(data["time"], data["u"], label="u")
    axes[3].set_xlabel("time (s)")
    axes[3].set_ylabel("V and u")
    axes[3].legend()
    axes[3].grid(True)

    fig.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Chapter5_Lesson4.cpp

/*
Chapter5_Lesson4.cpp
Lyapunov-based MRAC for a first-order plant.
Build:
    g++ -std=c++17 -O2 Chapter5_Lesson4.cpp -o Chapter5_Lesson4
*/

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

struct Parameters {
    double a = -0.4;
    double b = 1.5;
    double am = -1.2;
    double bm = 1.2;
    double gamma_x = 4.0;
    double gamma_r = 4.0;
    double dt = 0.002;
    double t_final = 30.0;
};

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

double reference(double t) {
    if (t < 10.0) {
        return 1.0;
    }
    if (t < 20.0) {
        return -0.5;
    }
    return 0.8;
}

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

State rhs(double t, const State& z, const Parameters& p) {
    const double x = z[0];
    const double xm = z[1];
    const double kx_hat = z[2];
    const double kr_hat = z[3];
    const double r = reference(t);
    const double e = x - xm;
    const double u = kx_hat * x + kr_hat * r;
    const double sign_b = (p.b > 0.0) ? 1.0 : -1.0;

    return State{
        p.a * x + p.b * u,
        p.am * xm + p.bm * r,
        -p.gamma_x * sign_b * e * x,
        -p.gamma_r * sign_b * e * r
    };
}

State rk4_step(double t, const State& z, double h, const Parameters& p) {
    const State k1 = rhs(t, z, p);
    const State k2 = rhs(t + 0.5 * h, add_scaled(z, k1, 0.5 * h), p);
    const State k3 = rhs(t + 0.5 * h, add_scaled(z, k2, 0.5 * h), p);
    const State k4 = rhs(t + h, add_scaled(z, k3, h), p);

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

void validate(const Parameters& p) {
    if (p.b == 0.0) {
        throw std::invalid_argument("The plant input gain b must be nonzero.");
    }
    if (p.am >= 0.0) {
        throw std::invalid_argument("The reference-model pole am must be negative.");
    }
    if (p.gamma_x <= 0.0 || p.gamma_r <= 0.0 || p.dt <= 0.0 || p.t_final <= 0.0) {
        throw std::invalid_argument("Adaptation gains, dt, and t_final must be positive.");
    }
}

int main() {
    try {
        const Parameters p;
        validate(p);

        const double kx_star = (p.am - p.a) / p.b;
        const double kr_star = p.bm / p.b;
        const long long steps = static_cast<long long>(std::llround(p.t_final / p.dt));

        State z{0.5, 0.0, 0.0, 0.0};
        std::ofstream csv("Chapter5_Lesson4_cpp.csv");
        if (!csv) {
            throw std::runtime_error("Could not create output CSV file.");
        }
        csv << "t,r,x,xm,e,u,kx_hat,kr_hat,V\n";
        csv << std::setprecision(12);

        double max_v_increase = -1.0e300;
        double previous_v = 0.0;

        for (long long i = 0; i <= steps; ++i) {
            const double t = i * p.dt;
            const double r = reference(t);
            const double x = z[0];
            const double xm = z[1];
            const double e = x - xm;
            const double u = z[2] * x + z[3] * r;
            const double kx_tilde = z[2] - kx_star;
            const double kr_tilde = z[3] - kr_star;
            const double V =
                0.5 * e * e
                + std::abs(p.b) * kx_tilde * kx_tilde / (2.0 * p.gamma_x)
                + std::abs(p.b) * kr_tilde * kr_tilde / (2.0 * p.gamma_r);

            if (i > 0) {
                max_v_increase = std::max(max_v_increase, V - previous_v);
            }
            previous_v = V;

            csv << t << ',' << r << ',' << x << ',' << xm << ',' << e << ','
                << u << ',' << z[2] << ',' << z[3] << ',' << V << '\n';

            if (i < steps) {
                z = rk4_step(t, z, p.dt, p);
            }
        }

        std::cout << std::setprecision(8)
                  << "Ideal kx*: " << kx_star << '\n'
                  << "Ideal kr*: " << kr_star << '\n'
                  << "Final tracking error: " << (z[0] - z[1]) << '\n'
                  << "Maximum V increase: " << max_v_increase << '\n'
                  << "Results written to Chapter5_Lesson4_cpp.csv\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter5_Lesson4.java

/*
Chapter5_Lesson4.java
Lyapunov-based MRAC for a first-order plant.

Build and run:
    javac Chapter5_Lesson4.java
    java Chapter5_Lesson4
*/

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

public final class Chapter5_Lesson4 {
    private static final class Parameters {
        final double a = -0.4;
        final double b = 1.5;
        final double am = -1.2;
        final double bm = 1.2;
        final double gammaX = 4.0;
        final double gammaR = 4.0;
        final double dt = 0.002;
        final double tFinal = 30.0;

        void validate() {
            if (b == 0.0) {
                throw new IllegalArgumentException("The plant input gain b must be nonzero.");
            }
            if (am >= 0.0) {
                throw new IllegalArgumentException("The reference-model pole am must be negative.");
            }
            if (gammaX <= 0.0 || gammaR <= 0.0 || dt <= 0.0 || tFinal <= 0.0) {
                throw new IllegalArgumentException(
                    "Adaptation gains, dt, and tFinal must be positive."
                );
            }
        }
    }

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

    private static double[] rhs(double t, double[] z, Parameters p) {
        double x = z[0];
        double xm = z[1];
        double kxHat = z[2];
        double krHat = z[3];
        double r = reference(t);
        double e = x - xm;
        double u = kxHat * x + krHat * r;
        double signB = Math.copySign(1.0, p.b);

        return new double[] {
            p.a * x + p.b * u,
            p.am * xm + p.bm * r,
            -p.gammaX * signB * e * x,
            -p.gammaR * signB * e * r
        };
    }

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

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

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

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Parameters p = new Parameters();
        p.validate();

        double kxStar = (p.am - p.a) / p.b;
        double krStar = p.bm / p.b;
        long steps = Math.round(p.tFinal / p.dt);
        double[] z = {0.5, 0.0, 0.0, 0.0};

        double previousV = 0.0;
        double maxVIncrease = -Double.MAX_VALUE;

        try (BufferedWriter writer =
                 new BufferedWriter(new FileWriter("Chapter5_Lesson4_java.csv"))) {
            writer.write("t,r,x,xm,e,u,kx_hat,kr_hat,V");
            writer.newLine();

            for (long i = 0; i <= steps; i++) {
                double t = i * p.dt;
                double r = reference(t);
                double x = z[0];
                double xm = z[1];
                double e = x - xm;
                double u = z[2] * x + z[3] * r;
                double kxTilde = z[2] - kxStar;
                double krTilde = z[3] - krStar;
                double V = 0.5 * e * e
                    + Math.abs(p.b) * kxTilde * kxTilde / (2.0 * p.gammaX)
                    + Math.abs(p.b) * krTilde * krTilde / (2.0 * p.gammaR);

                if (i > 0) {
                    maxVIncrease = Math.max(maxVIncrease, V - previousV);
                }
                previousV = V;

                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                    t, r, x, xm, e, u, z[2], z[3], V
                ));
                writer.newLine();

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

        System.out.printf("Ideal kx*: %.8f%n", kxStar);
        System.out.printf("Ideal kr*: %.8f%n", krStar);
        System.out.printf("Final tracking error: %.8e%n", z[0] - z[1]);
        System.out.printf("Maximum V increase: %.8e%n", maxVIncrease);
        System.out.println("Results written to Chapter5_Lesson4_java.csv");
    }
}

Chapter5_Lesson4.m

% Chapter5_Lesson4.m
% Lyapunov-based MRAC for a first-order plant using a fixed-step RK4 solver.
%
% Plant:          x_dot  = a*x + b*u
% Reference:      xm_dot = am*xm + bm*r
% Controller:     u      = kx_hat*x + kr_hat*r
% Adaptive laws:  kx_hat_dot = -gamma_x*sign(b)*e*x
%                 kr_hat_dot = -gamma_r*sign(b)*e*r

clear; clc; close all;

p.a = -0.4;
p.b = 1.5;
p.am = -1.2;
p.bm = 1.2;
p.gamma_x = 4.0;
p.gamma_r = 4.0;
p.dt = 0.002;
p.t_final = 30.0;

validate_parameters(p);

t = (0:p.dt:p.t_final).';
n = numel(t);
z = zeros(n, 4);              % [x, xm, kx_hat, kr_hat]
z(1, :) = [0.5, 0.0, 0.0, 0.0];
r = zeros(n, 1);
u = zeros(n, 1);
e = zeros(n, 1);
V = zeros(n, 1);

kx_star = (p.am - p.a) / p.b;
kr_star = p.bm / p.b;

for k = 1:n
    [~, u(k), e(k), r(k)] = mrac_rhs(t(k), z(k, :).', p);
    kx_tilde = z(k, 3) - kx_star;
    kr_tilde = z(k, 4) - kr_star;
    V(k) = 0.5 * e(k)^2 ...
        + abs(p.b) * kx_tilde^2 / (2 * p.gamma_x) ...
        + abs(p.b) * kr_tilde^2 / (2 * p.gamma_r);

    if k < n
        z(k + 1, :) = rk4_step(t(k), z(k, :).', p.dt, p).';
    end
end

fprintf('Ideal kx*: %.8f\n', kx_star);
fprintf('Ideal kr*: %.8f\n', kr_star);
fprintf('Final tracking error: %.8e\n', e(end));
fprintf('Maximum V increase: %.8e\n', max(diff(V)));

results = table(t, r, z(:,1), z(:,2), e, u, z(:,3), z(:,4), V, ...
    'VariableNames', {'t','r','x','xm','e','u','kx_hat','kr_hat','V'});
writetable(results, 'Chapter5_Lesson4_matlab.csv');

figure('Name', 'Chapter 5 Lesson 4 - Lyapunov MRAC');
tiledlayout(4, 1);

nexttile;
plot(t, z(:,1), 'LineWidth', 1.2); hold on;
plot(t, z(:,2), '--', 'LineWidth', 1.2);
plot(t, r, ':', 'LineWidth', 1.0);
grid on; ylabel('signals');
legend('plant x', 'model x_m', 'command r', 'Location', 'best');

nexttile;
plot(t, e, 'LineWidth', 1.2);
grid on; ylabel('error e');

nexttile;
plot(t, z(:,3), 'LineWidth', 1.2); hold on;
plot(t, z(:,4), 'LineWidth', 1.2);
yline(kx_star, '--');
yline(kr_star, ':');
grid on; ylabel('parameters');
legend('k_x hat', 'k_r hat', 'k_x^*', 'k_r^*', 'Location', 'best');

nexttile;
plot(t, V, 'LineWidth', 1.2); hold on;
plot(t, u, 'LineWidth', 1.0);
grid on; xlabel('time (s)'); ylabel('V and u');
legend('V', 'u', 'Location', 'best');

function validate_parameters(p)
    assert(p.b ~= 0, 'The plant input gain b must be nonzero.');
    assert(p.am < 0, 'The reference-model pole am must be negative.');
    assert(p.gamma_x > 0 && p.gamma_r > 0, ...
        'Adaptation gains must be positive.');
    assert(p.dt > 0 && p.t_final > 0, ...
        'dt and t_final must be positive.');
end

function r = command_signal(t)
    if t < 10
        r = 1.0;
    elseif t < 20
        r = -0.5;
    else
        r = 0.8;
    end
end

function [dz, u, e, r] = mrac_rhs(t, z, p)
    x = z(1);
    xm = z(2);
    kx_hat = z(3);
    kr_hat = z(4);

    r = command_signal(t);
    e = x - xm;
    u = kx_hat * x + kr_hat * r;
    sign_b = sign(p.b);

    dz = [
        p.a * x + p.b * u;
        p.am * xm + p.bm * r;
        -p.gamma_x * sign_b * e * x;
        -p.gamma_r * sign_b * e * r
    ];
end

function z_next = rk4_step(t, z, h, p)
    k1 = mrac_rhs(t, z, p);
    k2 = mrac_rhs(t + 0.5*h, z + 0.5*h*k1, p);
    k3 = mrac_rhs(t + 0.5*h, z + 0.5*h*k2, p);
    k4 = mrac_rhs(t + h, z + h*k3, p);
    z_next = z + (h/6) * (k1 + 2*k2 + 2*k3 + k4);
end

The following MATLAB script creates a Simulink model containing a vector integrator for \( [x,x_m,\hat{k}_x,\hat{k}_r]^\mathsf{T} \), a MATLAB Function block for the closed-loop differential equations, a clock, and workspace logging blocks. It then simulates and plots the result.

Chapter5_Lesson4_Simulink.m

% Chapter5_Lesson4_Simulink.m
% Programmatically builds and simulates a Simulink realization of the
% first-order Lyapunov MRAC. Requires Simulink.

clear; clc; close all;

model = 'Chapter5_Lesson4_Simulink_Model';

if ~license('test', 'Simulink')
    error('Simulink is required to run this script.');
end

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

new_system(model);
open_system(model);

% Blocks
add_block('simulink/Sources/Clock', [model '/Clock'], ...
    'Position', [40 75 70 105]);

add_block('simulink/Continuous/Integrator', [model '/State Integrator'], ...
    'InitialCondition', '[0.5;0;0;0]', ...
    'Position', [390 55 430 125]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/MRAC Dynamics'], ...
    'Position', [170 40 310 145]);

add_block('simulink/Sinks/To Workspace', [model '/State Log'], ...
    'VariableName', 'z_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 45 590 75]);

add_block('simulink/Sinks/To Workspace', [model '/Control Log'], ...
    'VariableName', 'u_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 90 590 120]);

add_block('simulink/Sinks/To Workspace', [model '/Error Log'], ...
    'VariableName', 'e_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 135 590 165]);

add_block('simulink/Sinks/To Workspace', [model '/Command Log'], ...
    'VariableName', 'r_log', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 180 590 210]);

% Configure the MATLAB Function block through the Stateflow API.
rt = sfroot;
chart = find(rt, '-isa', 'Stateflow.EMChart', ...
    'Path', [model '/MRAC Dynamics']);

chart.Script = sprintf([ ...
    'function [dz,u,e,r] = fcn(t,z)\n' ...
    '%%#codegen\n' ...
    'a = -0.4; b = 1.5; am = -1.2; bm = 1.2;\n' ...
    'gamma_x = 4.0; gamma_r = 4.0;\n' ...
    'if t < 10.0\n' ...
    '    r = 1.0;\n' ...
    'elseif t < 20.0\n' ...
    '    r = -0.5;\n' ...
    'else\n' ...
    '    r = 0.8;\n' ...
    'end\n' ...
    'x = z(1); xm = z(2); kx_hat = z(3); kr_hat = z(4);\n' ...
    'e = x - xm;\n' ...
    'u = kx_hat*x + kr_hat*r;\n' ...
    'sign_b = sign(b);\n' ...
    'dx = a*x + b*u;\n' ...
    'dxm = am*xm + bm*r;\n' ...
    'dkx = -gamma_x*sign_b*e*x;\n' ...
    'dkr = -gamma_r*sign_b*e*r;\n' ...
    'dz = [dx; dxm; dkx; dkr];\n' ...
    'end\n']);

% Connections
add_line(model, 'Clock/1', 'MRAC Dynamics/1', 'autorouting', 'on');
add_line(model, 'State Integrator/1', 'MRAC Dynamics/2', 'autorouting', 'on');
add_line(model, 'MRAC Dynamics/1', 'State Integrator/1', 'autorouting', 'on');
add_line(model, 'State Integrator/1', 'State Log/1', 'autorouting', 'on');
add_line(model, 'MRAC Dynamics/2', 'Control Log/1', 'autorouting', 'on');
add_line(model, 'MRAC Dynamics/3', 'Error Log/1', 'autorouting', 'on');
add_line(model, 'MRAC Dynamics/4', 'Command Log/1', 'autorouting', 'on');

set_param(model, ...
    'StopTime', '30', ...
    'SolverType', 'Variable-step', ...
    'Solver', 'ode45', ...
    'MaxStep', '0.01');

save_system(model);
sim(model);

t = z_log.time;
z = z_log.signals.values;
u = interp1(u_log.time, u_log.signals.values, t, 'linear', 'extrap');
e = interp1(e_log.time, e_log.signals.values, t, 'linear', 'extrap');
r = interp1(r_log.time, r_log.signals.values, t, 'previous', 'extrap');

figure('Name', 'Simulink Lyapunov MRAC');
tiledlayout(3,1);

nexttile;
plot(t, z(:,1), t, z(:,2), '--', t, r, ':', 'LineWidth', 1.1);
grid on; ylabel('signals');
legend('x', 'x_m', 'r', 'Location', 'best');

nexttile;
plot(t, e, 'LineWidth', 1.1);
grid on; ylabel('e');

nexttile;
plot(t, z(:,3), t, z(:,4), t, u, 'LineWidth', 1.1);
grid on; xlabel('time (s)'); ylabel('parameters / u');
legend('k_x hat', 'k_r hat', 'u', 'Location', 'best');

disp(['Saved Simulink model: ' model '.slx']);

Chapter5_Lesson4.nb

(* Chapter5_Lesson4.nb
   Lyapunov-based MRAC for a first-order plant.
   The notebook uses a fixed-step RK4 implementation from scratch. *)

ClearAll["Global`*"];

a = -0.4;
b = 1.5;
am = -1.2;
bm = 1.2;
gammaX = 4.0;
gammaR = 4.0;
dt = 0.002;
tFinal = 30.0;

If[b == 0, Print["b must be nonzero."]; Abort[]];
If[am >= 0, Print["am must be negative."]; Abort[]];

reference[t_?NumericQ] := Piecewise[{
    {1.0, t < 10.0},
    {-0.5, t < 20.0}
    }, 0.8];

rhs[t_?NumericQ, z_List] := Module[
    {x, xm, kxHat, krHat, r, e, u, signB},
    {x, xm, kxHat, krHat} = z;
    r = reference[t];
    e = x - xm;
    u = kxHat*x + krHat*r;
    signB = Sign[b];
    {
        a*x + b*u,
        am*xm + bm*r,
        -gammaX*signB*e*x,
        -gammaR*signB*e*r
    }
];

rk4Step[t_?NumericQ, z_List, h_?NumericQ] := Module[
    {k1, k2, k3, k4},
    k1 = rhs[t, z];
    k2 = rhs[t + h/2, z + (h/2)*k1];
    k3 = rhs[t + h/2, z + (h/2)*k2];
    k4 = rhs[t + h, z + h*k3];
    z + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
];

steps = Round[tFinal/dt];
time = N[Range[0, steps]*dt];
states = ConstantArray[0.0, {steps + 1, 4}];
states[[1]] = {0.5, 0.0, 0.0, 0.0};

Do[
    states[[k + 1]] = rk4Step[time[[k]], states[[k]], dt],
    {k, 1, steps}
];

kxStar = (am - a)/b;
krStar = bm/b;
command = reference /@ time;
error = states[[All, 1]] - states[[All, 2]];
control = MapThread[#1*#2 + #3*#4 &,
    {states[[All, 3]], states[[All, 1]],
     states[[All, 4]], command}];

lyapunov = MapThread[
    0.5*#1^2 + Abs[b]*(#2 - kxStar)^2/(2*gammaX)
      + Abs[b]*(#3 - krStar)^2/(2*gammaR) &,
    {error, states[[All, 3]], states[[All, 4]]}
];

Print["Ideal kx*: ", N[kxStar, 8]];
Print["Ideal kr*: ", N[krStar, 8]];
Print["Final tracking error: ", ScientificForm[Last[error], 6]];
Print["Maximum V increase: ", ScientificForm[Max[Differences[lyapunov]], 6]];

trackingPlot = ListLinePlot[
    {
        Transpose[{time, states[[All, 1]]}],
        Transpose[{time, states[[All, 2]]}],
        Transpose[{time, command}]
    },
    PlotLegends -> {"plant x", "model xm", "command r"},
    Frame -> True,
    FrameLabel -> {"time (s)", "signals"},
    PlotRange -> All,
    ImageSize -> Large
];

errorPlot = ListLinePlot[
    Transpose[{time, error}],
    Frame -> True,
    FrameLabel -> {"time (s)", "tracking error e"},
    PlotRange -> All,
    ImageSize -> Large
];

parameterPlot = ListLinePlot[
    {
        Transpose[{time, states[[All, 3]]}],
        Transpose[{time, states[[All, 4]]}],
        Transpose[{time, ConstantArray[kxStar, Length[time]]}],
        Transpose[{time, ConstantArray[krStar, Length[time]]}]
    },
    PlotLegends -> {"kxHat", "krHat", "kxStar", "krStar"},
    Frame -> True,
    FrameLabel -> {"time (s)", "parameters"},
    PlotRange -> All,
    ImageSize -> Large
];

lyapunovPlot = ListLinePlot[
    {
        Transpose[{time, lyapunov}],
        Transpose[{time, control}]
    },
    PlotLegends -> {"V", "u"},
    Frame -> True,
    FrameLabel -> {"time (s)", "V and u"},
    PlotRange -> All,
    ImageSize -> Large
];

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

1. The plant state initially differs from the model state, producing a nonzero tracking error.
2. The error–regressor products move the feedback and feedforward estimates.
3. The plant state approaches the reference-model state after each command change.
4. The computed Lyapunov function should be non-increasing in the ideal continuous-time system. A numerical solver may produce extremely small positive increments because of finite step size and round-off.
5. The final gain estimates need not equal the ideal values exactly, even when the tracking error is very small.

A useful numerical diagnostic is \( \max_k[V(t_{k+1})-V(t_k)] \). Reducing the integration step should reduce spurious positive increments caused by numerical error.

Problem 1 — Derive the adaptive laws for a positive input gain. Assume \( b > 0 \). Starting from \( \dot{e}=a_m e+b\tilde{k}_x x+b\tilde{k}_r r \) and

\[ V=\frac{1}{2}e^2+ \frac{b}{2\gamma_x}\tilde{k}_x^2+ \frac{b}{2\gamma_r}\tilde{k}_r^2, \]

derive update laws that give \( \dot{V}=a_m e^2 \).

Solution. Differentiate:

\[ \dot{V} = a_m e^2 + b\tilde{k}_x \left(e x+\frac{1}{\gamma_x}\dot{\hat{k}}_x\right) + b\tilde{k}_r \left(e r+\frac{1}{\gamma_r}\dot{\hat{k}}_r\right). \]

Set each bracket to zero:

\[ \dot{\hat{k}}_x=-\gamma_x e x, \qquad \dot{\hat{k}}_r=-\gamma_r e r. \]

Substitution cancels both cross terms and gives \( \dot{V}=a_m e^2\le 0 \).

Problem 2 — Negative plant input gain. Let \( b=-2 \), \( \gamma_x=3 \), and \( \gamma_r=5 \). Write the adaptive laws explicitly.

Solution. Since \( s_b=-1 \),

\[ \dot{\hat{k}}_x = -3(-1)e x = 3 e x, \qquad \dot{\hat{k}}_r = -5(-1)e r = 5 e r. \]

The signs are opposite to the positive-\( b \) case. The proof still uses the positive weights \( |b|/(2\gamma_x) \) and \( |b|/(2\gamma_r) \).

Problem 3 — Compute the ideal gains. For \( a=0.5 \), \( b=2 \), \( a_m=-1.5 \), and \( b_m=3 \), determine the ideal controller gains and verify the matched closed-loop model.

Solution.

\[ k_x^\ast = \frac{a_m-a}{b} = \frac{-1.5-0.5}{2} = -1, \qquad k_r^\ast = \frac{b_m}{b} = \frac{3}{2} = 1.5. \]

With these values,

\[ a+b k_x^\ast = 0.5+2(-1) = -1.5 = a_m, \qquad b k_r^\ast = 2(1.5) = 3 = b_m. \]

Hence the ideal plant closed loop is \( \dot{x}=-1.5x+3r \), identical to the reference model.

Problem 4 — Establish an integral error bound. Suppose \( a_m=-2 \) and \( V(0)=6 \). Find an upper bound on \( \int_0^\infty e^2(t)\,dt \).

Solution. Since

\[ \dot{V}=-|a_m|e^2=-2e^2, \]

integration gives

\[ 2\int_0^\infty e^2(t)\,dt \le V(0)=6. \]

Therefore,

\[ \int_0^\infty e^2(t)\,dt\le 3. \]

Problem 5 — Why tracking does not imply unique parameter convergence. Let \( r(t)=0 \), \( x(0)=x_m(0)=0 \), and let the initial adaptive gains be arbitrary finite constants. Determine the subsequent tracking error and parameter updates.

Solution. Initially,

\[ e(0)=0,\qquad x(0)=0,\qquad r(0)=0. \]

The control is zero, the plant and model remain at the origin, and hence \( e(t)=0 \) for all time. The update laws satisfy

\[ \dot{\hat{k}}_x=-\gamma_x s_b e x=0, \qquad \dot{\hat{k}}_r=-\gamma_r s_b e r=0. \]

Thus the gains remain at their arbitrary initial values. Perfect tracking has occurred, but the ideal gains have not been identified. This example separates the tracking objective from the parameter-identification objective.

Problem 6 — Discrete implementation check. A forward-Euler implementation uses

\[ \hat{k}_x[k+1] = \hat{k}_x[k] - h\gamma_x s_b e[k]x[k]. \]

Explain why the continuous-time proof does not automatically imply that the sampled implementation has a non-increasing Lyapunov sequence for every step size \( h \).

Solution. The continuous proof uses the exact derivative

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

A finite Euler step introduces truncation terms of order \( h^2 \) in the one-step change of \( V \). For a large step, these omitted terms can dominate the desired negative first-order term and cause \( V[k+1]-V[k] > 0 \). Therefore, the step must be selected sufficiently small, and a higher-order method or a separately derived discrete-time adaptive law may be preferable.

For the first-order uncertain plant, the direct MRAC controller \( u=\hat{k}_x x+\hat{k}_r r \) produces an error model containing products of tracking error, parameter error, and regressors. A composite Lyapunov function adds weighted parameter-error energy to the tracking-error energy. Choosing \( \dot{\hat{k}}_x=-\gamma_x s_b e x \) and \( \dot{\hat{k}}_r=-\gamma_r s_b e r \) cancels the indefinite cross terms exactly and gives \( \dot{V}=-|a_m|e^2 \). This proves boundedness of all closed-loop signals and, using Barbalat’s lemma, asymptotic convergence of the tracking error. The result depends on exact matching, a stable reference model, and knowledge of the plant input-gain sign.