After completing this lesson, students should be able to:

• derive the exact tracking-error model for a direct adaptive controller;
• construct a composite state–parameter Lyapunov function;
• derive an adaptive law by cancellation of indefinite cross terms;
• prove boundedness and asymptotic tracking using Barbalat's lemma;
• distinguish tracking convergence from parameter convergence; and
• verify the theory numerically in five programming environments.

The objective is not merely to observe stable-looking plots. We seek a proof valid for arbitrary initial plant states, model states, and controller gains under the assumptions stated below.

flowchart TD
  R["Bounded command r"] --> RM["Stable reference model"]
  RM --> XM["Desired state xm"]
  R --> REG["Build regressor"]
  X["Measured plant state x"] --> REG
  REG --> C["Controller with adaptive gains"]
  C --> U["Control input u"]
  U --> P["Uncertain plant"]
  P --> X
  X --> E["Tracking error e = x - xm"]
  XM --> E
  E --> L["Lyapunov gain update"]
  REG --> L
  L --> C
        

Continue with the second-order canonical plant from Lesson 2:

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=-a_0x_1-a_1x_2+b u,\qquad b\neq0. \]

Equivalently,

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

The coefficients \(a_0\), \(a_1\), and the magnitude of \(b\) are unknown. Only \(\operatorname{sgn}(b)\) is assumed known. The unknown-control-direction problem is postponed to Chapter 16.

The desired closed-loop dynamics are specified by

\[ \dot{\mathbf{x} }_m=A_m\mathbf{x}_m+B_mr, \]

\[ A_m= \begin{bmatrix} 0&1\\ -\omega_n^2&-2\zeta\omega_n \end{bmatrix},\qquad B_m= \begin{bmatrix} 0\\\omega_n^2 \end{bmatrix}, \]

with \(\omega_n>0\) and \(\zeta>0\), so \(A_m\) is Hurwitz.

The direct controller is

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

The matching assumption requires the existence of a constant ideal gain vector \(\boldsymbol{\theta}^*\) satisfying

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

Coefficient matching gives

\[ \theta_1^*=\frac{a_0-\omega_n^2}{b},\qquad \theta_2^*=\frac{a_1-2\zeta\omega_n}{b},\qquad \theta_r^*=\frac{\omega_n^2}{b}. \]

These values are used only in the proof and offline validation. The direct regulator neither computes them nor estimates the plant coefficients.

Define

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

Add and subtract the ideal control:

\[ u= \boldsymbol{\theta}^{*T}\boldsymbol{\phi} + \widetilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Then

\[ \begin{aligned} \dot{\mathbf{e} } &= A_p\mathbf{x}+Bb u-A_m\mathbf{x}_m-B_mr\\ &= A_p\mathbf{x} +Bb\boldsymbol{\theta}^{*T}\boldsymbol{\phi} -A_m\mathbf{x}_m-B_mr\\ &\quad +Bb\widetilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \end{aligned} \]

The matching identities reduce this to

\[ \boxed{ \dot{\mathbf{e} } =A_m\mathbf{e} +Bb\widetilde{\boldsymbol{\theta} }^T\boldsymbol{\phi} }. \]

The first term is exponentially stable. The second is an indefinite bilinear coupling between tracking error and gain error. The adaptive law will be chosen to cancel this coupling in the Lyapunov derivative.

Since \(A_m\) is Hurwitz, for every symmetric \(Q\succ0\), there exists a unique symmetric \(P\succ0\) satisfying

\[ A_m^TP+PA_m=-Q. \]

Let \(\Gamma=\Gamma^T\succ0\) be the adaptation-gain matrix and choose

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

The factor \(|b|\) is essential because it permits cancellation using the sign, rather than the unknown magnitude, of the input gain. Positive definiteness gives

\[ \lambda_{\min}(P)\|\mathbf{e}\|^2 + |b|\lambda_{\min}(\Gamma^{-1}) \|\widetilde{\boldsymbol{\theta} }\|^2 \le V, \]

\[ V\le \lambda_{\max}(P)\|\mathbf{e}\|^2 + |b|\lambda_{\max}(\Gamma^{-1}) \|\widetilde{\boldsymbol{\theta} }\|^2. \]

Because \(\boldsymbol{\theta}^*\) is constant, \(\dot{\widetilde{\boldsymbol{\theta} } } =\dot{\hat{\boldsymbol{\theta} } }\). Differentiating \(V\) gives

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

Using the Lyapunov equation,

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

Choose

\[ \boxed{ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma\boldsymbol{\phi} \mathbf{e}^TPB\operatorname{sgn}(b) }. \]

Then

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

Hence the cross terms cancel exactly:

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

flowchart TD
  A["Choose Hurwitz reference matrix Am"] --> B["Select positive definite Q"]
  B --> C["Solve Am' P + P Am = -Q"]
  C --> D["Construct state plus gain-error energy V"]
  D --> E["Differentiate V"]
  E --> F["Locate the mixed state-gain term"]
  F --> G["Choose gain update for exact cancellation"]
  G --> H["Obtain Vdot = -e' Q e"]
  H --> I["Prove boundedness"]
  I --> J["Apply Barbalat for asymptotic tracking"]
        

Theorem.

Assume that \(A_m\) is Hurwitz, the matching gains exist and are constant, \(\operatorname{sgn}(b)\) is known, the command is bounded, and the closed-loop differential equations have a unique solution. Then every closed-loop signal is bounded and \(\mathbf{e}(t)\) converges to zero.

Step 1: Lyapunov boundedness.

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

Therefore \(\mathbf{e}\in L_\infty\) and \(\widetilde{\boldsymbol{\theta} }\in L_\infty\). Since the ideal gains are constant, \(\hat{\boldsymbol{\theta} }\in L_\infty\).

Step 2: square integrability.

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

Consequently,

\[ \lambda_{\min}(Q) \int_0^t\|\mathbf{e}(\tau)\|^2d\tau \le V(0)-V(t)\le V(0), \]

so \(\mathbf{e}\in L_2\).

Step 3: bounded derivative.

A stable reference model with bounded input has bounded \(\mathbf{x}_m\). Since \(\mathbf{x}=\mathbf{e}+\mathbf{x}_m\), the plant state is bounded. Thus the regressor is bounded. The adaptive law then gives bounded \(\dot{\hat{\boldsymbol{\theta} } }\), and the error equation gives bounded \(\dot{\mathbf{e} }\).

Step 4: Barbalat's lemma.

Bounded \(\dot{\mathbf{e} }\) makes \(\mathbf{e}\) uniformly continuous. Together with \(\mathbf{e}\in L_2\cap L_\infty\), Barbalat's lemma gives

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

This proves global boundedness and asymptotic model following.

The derivative \(\dot{V}=-\mathbf{e}^TQ\mathbf{e}\) is negative definite in the tracking error but contains no negative-definite term in the gain error. Therefore the basic proof establishes

\[ \mathbf{e}(t)\rightarrow\mathbf{0},\qquad \widetilde{\boldsymbol{\theta} }(t)\text{ is bounded}, \]

but not necessarily

\[ \widetilde{\boldsymbol{\theta} }(t)\rightarrow\mathbf{0}. \]

A standard sufficient condition for parameter convergence is persistent excitation: there must exist \(T>0\) and \(\alpha>0\) such that

\[ \int_t^{t+T} \boldsymbol{\phi}(\tau)\boldsymbol{\phi}^T(\tau)d\tau \succeq\alpha I,\qquad t\ge0. \]

A step command can produce excellent tracking without independently exciting every controller direction. Therefore final adaptive gains may differ from the ideal vector even when the tracking error is nearly zero.

Choice of \(Q\).

The matrix \(Q\succ0\) determines \(P\). It shapes the Lyapunov metric and the scalar adaptation signal \(\mathbf{e}^TPB\).

Choice of \(\Gamma\).

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

Larger adaptation gains produce faster gain motion, but they also magnify noise, discretization error, and transient control effort.

Optional normalization.

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma \frac{\boldsymbol{\phi} } {1+\boldsymbol{\phi}^T\boldsymbol{\phi} } \mathbf{e}^TPB\operatorname{sgn}(b). \]

Any normalization, projection, leakage, or dead-zone modification changes the derivative calculation and must be analyzed rather than inserted without proof.

Continuous versus numerical time.

The theorem is continuous-time. A finite-step solver may produce tiny sampled increases in \(V\). Such increases should decrease as the integration step is reduced. A formal sampled-data result requires the discrete-time analysis of Chapter 14.

Use

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=-1.6x_1-1.1x_2+1.4u, \]

with

\[ \omega_n=2,\qquad\zeta=0.8,\qquad Q=I_2. \]

Then

\[ A_m= \begin{bmatrix} 0&1\\-4&-3.2 \end{bmatrix}. \]

Write \(P=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix}\). Expanding \(A_m^TP+PA_m=-I\) gives

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

Therefore,

\[ P= \begin{bmatrix} 1.18125&0.125\\ 0.125&0.1953125 \end{bmatrix}\succ0. \]

The validation-only ideal gains are

\[ \boldsymbol{\theta}^* = \begin{bmatrix} -1.7142857\\ -1.5\\ 2.8571429 \end{bmatrix}. \]

Choose

\[ \Gamma=\operatorname{diag}(10,10,5), \]

and use the bounded multisine command

\[ r(t)=0.8\sin(0.5t)+0.5\sin(1.3t)+0.3\sin(2.1t). \]

Every implementation below integrates the same seven coupled states and exports numerical results for independent verification.

Libraries: NumPy and Matplotlib. RK4 is implemented from scratch.

Chapter13_Lesson3.py

# Chapter13_Lesson3.py
"""
Direct self-tuning regulator stability demonstration.

Plant:
    x1_dot = x2
    x2_dot = -a0*x1 - a1*x2 + b*u

Reference model:
    xm1_dot = xm2
    xm2_dot = -wn^2*xm1 - 2*zeta*wn*xm2 + wn^2*r

Direct adaptive controller:
    u = theta_hat.T @ [x1, x2, r]

Lyapunov update:
    theta_hat_dot = -Gamma @ phi * (e.T @ P @ B) * sign(b)

The true ideal controller gains are used only to evaluate the Lyapunov
function and are never supplied to the controller.
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    a0: float = 1.6
    a1: float = 1.1
    b: float = 1.4
    wn: float = 2.0
    zeta: float = 0.8
    dt: float = 0.001
    tf: float = 40.0


P = np.array([[1.18125, 0.125], [0.125, 0.1953125]], dtype=float)
B = np.array([0.0, 1.0], dtype=float)
Q = np.eye(2)
GAMMA = np.diag([10.0, 10.0, 5.0])


def reference(t: float) -> float:
    """Bounded multisine command; it provides richer excitation than a step."""
    return (
        0.8 * np.sin(0.5 * t)
        + 0.5 * np.sin(1.3 * t)
        + 0.3 * np.sin(2.1 * t)
    )


def ideal_gains(p: Parameters) -> np.ndarray:
    return np.array(
        [
            (p.a0 - p.wn**2) / p.b,
            (p.a1 - 2.0 * p.zeta * p.wn) / p.b,
            p.wn**2 / p.b,
        ],
        dtype=float,
    )


def dynamics(t: float, z: np.ndarray, p: Parameters) -> np.ndarray:
    x = z[0:2]
    xm = z[2:4]
    theta_hat = z[4:7]

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

    x_dot = np.array(
        [x[1], -p.a0 * x[0] - p.a1 * x[1] + p.b * u], dtype=float
    )
    xm_dot = np.array(
        [xm[1], -(p.wn**2) * xm[0] - 2.0 * p.zeta * p.wn * xm[1]
         + (p.wn**2) * r],
        dtype=float,
    )

    e = x - xm
    scalar_error = float(e @ P @ B)
    theta_hat_dot = -GAMMA @ phi * scalar_error * np.sign(p.b)

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


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


def main() -> None:
    p = Parameters()
    n = int(round(p.tf / p.dt)) + 1
    t = np.linspace(0.0, p.tf, n)

    z = np.zeros((n, 7), dtype=float)
    z[0, 0:2] = np.array([0.4, -0.2])

    theta_star = ideal_gains(p)
    lyapunov = np.zeros(n)
    vdot_theory = np.zeros(n)
    control = np.zeros(n)
    command = np.zeros(n)

    gamma_inv = np.linalg.inv(GAMMA)

    for k in range(n):
        x = z[k, 0:2]
        xm = z[k, 2:4]
        theta_hat = z[k, 4:7]
        e = x - xm
        theta_tilde = theta_hat - theta_star
        command[k] = reference(t[k])
        phi = np.array([x[0], x[1], command[k]])
        control[k] = float(theta_hat @ phi)

        lyapunov[k] = float(
            e @ P @ e + abs(p.b) * theta_tilde @ gamma_inv @ theta_tilde
        )
        vdot_theory[k] = -float(e @ Q @ e)

        if k + 1 < n:
            z[k + 1] = rk4_step(dynamics, t[k], z[k], p.dt, p)

    e = z[:, 0:2] - z[:, 2:4]
    rms_last = float(np.sqrt(np.mean(e[int(0.75 * n):, 0] ** 2)))
    max_v_increase = float(np.max(np.diff(lyapunov)))

    print("Ideal gains (validation only):", theta_star)
    print("Final adaptive gains:", z[-1, 4:7])
    print("RMS position tracking error over final quarter:", rms_last)
    print("Maximum sampled increase in V:", max_v_increase)

    with open("Chapter13_Lesson3_results.csv", "w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow(
            [
                "t", "r", "x1", "x2", "xm1", "xm2", "e1", "e2",
                "theta1_hat", "theta2_hat", "theta_r_hat", "u", "V",
                "Vdot_theory",
            ]
        )
        for k in range(n):
            writer.writerow(
                [
                    t[k], command[k], z[k, 0], z[k, 1], z[k, 2], z[k, 3],
                    e[k, 0], e[k, 1], z[k, 4], z[k, 5], z[k, 6],
                    control[k], lyapunov[k], vdot_theory[k],
                ]
            )

    plt.figure()
    plt.plot(t, z[:, 0], label="plant x1")
    plt.plot(t, z[:, 2], "--", label="model xm1")
    plt.xlabel("Time (s)")
    plt.ylabel("Position-like state")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, e[:, 0], label="e1")
    plt.plot(t, e[:, 1], label="e2")
    plt.xlabel("Time (s)")
    plt.ylabel("Tracking error")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, z[:, 4], label="theta1_hat")
    plt.plot(t, z[:, 5], label="theta2_hat")
    plt.plot(t, z[:, 6], label="theta_r_hat")
    plt.xlabel("Time (s)")
    plt.ylabel("Adaptive gain")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, lyapunov, label="V")
    plt.xlabel("Time (s)")
    plt.ylabel("Lyapunov function")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.show()


if __name__ == "__main__":
    main()

The C++17 implementation uses only the standard library and exports CSV data.

Chapter13_Lesson3.cpp

// Chapter13_Lesson3.cpp
// C++17 implementation of the Lyapunov-stable direct STR.
// Compile: g++ -std=c++17 -O2 Chapter13_Lesson3.cpp -o Chapter13_Lesson3

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

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

constexpr double A0 = 1.6;
constexpr double A1 = 1.1;
constexpr double BGAIN = 1.4;
constexpr double WN = 2.0;
constexpr double ZETA = 0.8;
constexpr double DT = 0.001;
constexpr double TF = 40.0;

constexpr double P11 = 1.18125;
constexpr double P12 = 0.125;
constexpr double P22 = 0.1953125;
constexpr std::array<double, 3> GAMMA{10.0, 10.0, 5.0};

double reference(double t) {
    return 0.8 * std::sin(0.5 * t)
         + 0.5 * std::sin(1.3 * t)
         + 0.3 * std::sin(2.1 * t);
}

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

State dynamics(double t, const State& z) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double r = reference(t);

    const std::array<double, 3> phi{x1, x2, r};
    const double u = z[4] * phi[0] + z[5] * phi[1] + z[6] * phi[2];

    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;
    const double scalarError = P12 * e1 + P22 * e2;
    const double signB = (BGAIN >= 0.0) ? 1.0 : -1.0;

    State dz{};
    dz[0] = x2;
    dz[1] = -A0 * x1 - A1 * x2 + BGAIN * u;
    dz[2] = xm2;
    dz[3] = -WN * WN * xm1 - 2.0 * ZETA * WN * xm2 + WN * WN * r;
    dz[4] = -GAMMA[0] * phi[0] * scalarError * signB;
    dz[5] = -GAMMA[1] * phi[1] * scalarError * signB;
    dz[6] = -GAMMA[2] * phi[2] * scalarError * signB;
    return dz;
}

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

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

int main() {
    const std::array<double, 3> thetaStar{
        (A0 - WN * WN) / BGAIN,
        (A1 - 2.0 * ZETA * WN) / BGAIN,
        WN * WN / BGAIN
    };

    State z{0.4, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0};
    std::ofstream csv("Chapter13_Lesson3_results_cpp.csv");
    if (!csv) {
        throw std::runtime_error("Could not create CSV output.");
    }
    csv << std::setprecision(15);
    csv << "t,r,x1,x2,xm1,xm2,e1,e2,theta1_hat,theta2_hat,theta_r_hat,u,V\n";

    const long steps = static_cast<long>(std::llround(TF / DT));
    double squaredErrorTail = 0.0;
    long tailCount = 0;

    for (long k = 0; k <= steps; ++k) {
        const double t = k * DT;
        const double r = reference(t);
        const double e1 = z[0] - z[2];
        const double e2 = z[1] - z[3];
        const double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;

        const double d1 = z[4] - thetaStar[0];
        const double d2 = z[5] - thetaStar[1];
        const double d3 = z[6] - thetaStar[2];
        const double stateTerm = P11 * e1 * e1 + 2.0 * P12 * e1 * e2 + P22 * e2 * e2;
        const double parameterTerm = std::abs(BGAIN) *
            (d1 * d1 / GAMMA[0] + d2 * d2 / GAMMA[1] + d3 * d3 / GAMMA[2]);
        const double V = stateTerm + parameterTerm;

        csv << t << ',' << r << ',' << z[0] << ',' << z[1] << ','
            << z[2] << ',' << z[3] << ',' << e1 << ',' << e2 << ','
            << z[4] << ',' << z[5] << ',' << z[6] << ',' << u << ',' << V << '\n';

        if (k >= (3 * steps) / 4) {
            squaredErrorTail += e1 * e1;
            ++tailCount;
        }
        if (k < steps) {
            z = rk4Step(t, z);
        }
    }

    std::cout << std::setprecision(8);
    std::cout << "Ideal gains (validation only): "
              << thetaStar[0] << ", " << thetaStar[1] << ", " << thetaStar[2] << '\n';
    std::cout << "Final adaptive gains: "
              << z[4] << ", " << z[5] << ", " << z[6] << '\n';
    std::cout << "Tail RMS position error: "
              << std::sqrt(squaredErrorTail / static_cast<double>(tailCount)) << '\n';
    return 0;
}

The Java 17 implementation uses primitive arrays, a from-scratch RK4 integrator, and standard file I/O.

Chapter13_Lesson3.java

// Chapter13_Lesson3.java
// Java 17 implementation of the Lyapunov-stable direct STR.
// Compile: javac Chapter13_Lesson3.java
// Run:     java Chapter13_Lesson3

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

public final class Chapter13_Lesson3 {
    private static final double A0 = 1.6;
    private static final double A1 = 1.1;
    private static final double B_GAIN = 1.4;
    private static final double WN = 2.0;
    private static final double ZETA = 0.8;
    private static final double DT = 0.001;
    private static final double TF = 40.0;

    private static final double P11 = 1.18125;
    private static final double P12 = 0.125;
    private static final double P22 = 0.1953125;
    private static final double[] GAMMA = {10.0, 10.0, 5.0};

    private Chapter13_Lesson3() {}

    private static double reference(double t) {
        return 0.8 * Math.sin(0.5 * t)
             + 0.5 * Math.sin(1.3 * t)
             + 0.3 * Math.sin(2.1 * t);
    }

    private static double[] dynamics(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double r = reference(t);

        double[] phi = {x1, x2, r};
        double u = z[4] * phi[0] + z[5] * phi[1] + z[6] * phi[2];

        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double scalarError = P12 * e1 + P22 * e2;
        double signB = Math.copySign(1.0, B_GAIN);

        double[] dz = new double[7];
        dz[0] = x2;
        dz[1] = -A0 * x1 - A1 * x2 + B_GAIN * u;
        dz[2] = xm2;
        dz[3] = -WN * WN * xm1 - 2.0 * ZETA * WN * xm2 + WN * WN * r;
        dz[4] = -GAMMA[0] * phi[0] * scalarError * signB;
        dz[5] = -GAMMA[1] * phi[1] * scalarError * signB;
        dz[6] = -GAMMA[2] * phi[2] * scalarError * signB;
        return dz;
    }

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

    private static double[] rk4Step(double t, double[] z) {
        double[] k1 = dynamics(t, z);
        double[] k2 = dynamics(t + 0.5 * DT, addScaled(z, k1, 0.5 * DT));
        double[] k3 = dynamics(t + 0.5 * DT, addScaled(z, k2, 0.5 * DT));
        double[] k4 = dynamics(t + DT, addScaled(z, k3, DT));

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

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

        double[] thetaStar = {
            (A0 - WN * WN) / B_GAIN,
            (A1 - 2.0 * ZETA * WN) / B_GAIN,
            WN * WN / B_GAIN
        };
        double[] z = {0.4, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0};

        long steps = Math.round(TF / DT);
        double squaredErrorTail = 0.0;
        long tailCount = 0;

        try (BufferedWriter writer = Files.newBufferedWriter(
                Path.of("Chapter13_Lesson3_results_java.csv"))) {
            writer.write("t,r,x1,x2,xm1,xm2,e1,e2,theta1_hat,theta2_hat,theta_r_hat,u,V");
            writer.newLine();

            for (long k = 0; k <= steps; k++) {
                double t = k * DT;
                double r = reference(t);
                double e1 = z[0] - z[2];
                double e2 = z[1] - z[3];
                double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;

                double d1 = z[4] - thetaStar[0];
                double d2 = z[5] - thetaStar[1];
                double d3 = z[6] - thetaStar[2];
                double stateTerm = P11 * e1 * e1 + 2.0 * P12 * e1 * e2 + P22 * e2 * e2;
                double parameterTerm = Math.abs(B_GAIN) *
                    (d1 * d1 / GAMMA[0] + d2 * d2 / GAMMA[1] + d3 * d3 / GAMMA[2]);
                double v = stateTerm + parameterTerm;

                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f,%.12f,%.12f",
                    t, r, z[0], z[1], z[2], z[3], e1, e2,
                    z[4], z[5], z[6], u, v
                ));
                writer.newLine();

                if (k >= (3 * steps) / 4) {
                    squaredErrorTail += e1 * e1;
                    tailCount++;
                }
                if (k < steps) {
                    z = rk4Step(t, z);
                }
            }
        }

        System.out.printf(
            Locale.US,
            "Ideal gains (validation only): %.8f, %.8f, %.8f%n",
            thetaStar[0], thetaStar[1], thetaStar[2]
        );
        System.out.printf(
            Locale.US,
            "Final adaptive gains: %.8f, %.8f, %.8f%n",
            z[4], z[5], z[6]
        );
        System.out.printf(
            Locale.US,
            "Tail RMS position error: %.8e%n",
            Math.sqrt(squaredErrorTail / tailCount)
        );
    }
}

The MATLAB file uses ode45. The Simulink builder creates a seven-state continuous model with a vector Integrator and MATLAB Function block.

Chapter13_Lesson3.m

% Chapter13_Lesson3.m
% Lyapunov-stable direct self-tuning regulator.
% Requires base MATLAB only. The simulation uses ode45.

clear; clc; close all;

p.a0 = 1.6;
p.a1 = 1.1;
p.b = 1.4;
p.wn = 2.0;
p.zeta = 0.8;
p.P = [1.18125, 0.125; 0.125, 0.1953125];
p.B = [0; 1];
p.Q = eye(2);
p.Gamma = diag([10, 10, 5]);

thetaStar = [(p.a0 - p.wn^2)/p.b;
             (p.a1 - 2*p.zeta*p.wn)/p.b;
             p.wn^2/p.b];

z0 = [0.4; -0.2; 0; 0; 0; 0; 0];
tspan = linspace(0, 40, 40001);
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);

[t, z] = ode45(@(t,z) closedLoop(t, z, p), tspan, z0, opts);

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

V = zeros(size(t));
VdotTheory = zeros(size(t));
GammaInv = inv(p.Gamma);
for k = 1:numel(t)
    thetaTilde = thetaHat(k,:).' - thetaStar;
    V(k) = e(k,:)*p.P*e(k,:).' ...
         + abs(p.b)*thetaTilde.'*GammaInv*thetaTilde;
    VdotTheory(k) = -e(k,:)*p.Q*e(k,:).';
end

tail = round(0.75*numel(t)):numel(t);
fprintf('Ideal gains (validation only): %.8f %.8f %.8f\n', thetaStar);
fprintf('Final adaptive gains: %.8f %.8f %.8f\n', thetaHat(end,:));
fprintf('Tail RMS position error: %.8e\n', rms(e(tail,1)));
fprintf('Maximum sampled increase in V: %.8e\n', max(diff(V)));

results = table(t, r, x(:,1), x(:,2), xm(:,1), xm(:,2), ...
    e(:,1), e(:,2), thetaHat(:,1), thetaHat(:,2), thetaHat(:,3), ...
    u, V, VdotTheory, ...
    'VariableNames', {'t','r','x1','x2','xm1','xm2','e1','e2', ...
    'theta1_hat','theta2_hat','theta_r_hat','u','V','Vdot_theory'});
writetable(results, 'Chapter13_Lesson3_results_matlab.csv');

figure;
plot(t, x(:,1), t, xm(:,1), '--', 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('Position-like state');
legend('Plant x_1','Model x_{m1}','Location','best');
title('Direct STR tracking');

figure;
plot(t, e, 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('Tracking error');
legend('e_1','e_2','Location','best');

figure;
plot(t, thetaHat, 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('Adaptive gain');
legend('\theta_1','\theta_2','\theta_r','Location','best');

figure;
plot(t, V, 'LineWidth', 1.1);
grid on; xlabel('Time (s)'); ylabel('V');
title('Composite Lyapunov function');

function dz = closedLoop(t, z, p)
    x = z(1:2);
    xm = z(3:4);
    thetaHat = z(5:7);

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

    xdot = [x(2);
            -p.a0*x(1) - p.a1*x(2) + p.b*u];
    xmdot = [xm(2);
             -p.wn^2*xm(1) - 2*p.zeta*p.wn*xm(2) + p.wn^2*r];

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

    dz = [xdot; xmdot; thetaHatDot];
end

function r = referenceSignal(t)
    r = 0.8*sin(0.5*t) + 0.5*sin(1.3*t) + 0.3*sin(2.1*t);
end

Chapter13_Lesson3_Simulink.m

% Chapter13_Lesson3_Simulink.m
% Programmatically builds a continuous-time Simulink model for the
% seven-state adaptive system and runs a 40-second simulation.
%
% Requirements: Simulink.
%
% State vector:
% z = [x1; x2; xm1; xm2; theta1_hat; theta2_hat; theta_r_hat]

clear; clc;

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

new_system(model);
open_system(model);

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

add_block('simulink/Continuous/Integrator', [model '/StateIntegrator'], ...
    'InitialCondition', '[0.4; -0.2; 0; 0; 0; 0; 0]', ...
    'Position', [380 75 425 115]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/AdaptiveDerivatives'], ...
    'Position', [155 55 315 135]);

add_block('simulink/Sinks/To Workspace', [model '/StateLog'], ...
    'VariableName', 'zLog', 'SaveFormat', 'Structure With Time', ...
    'Position', [500 65 595 95]);

add_block('simulink/Sinks/Scope', [model '/StateScope'], ...
    'Position', [500 120 535 150]);

add_line(model, 'Clock/1', 'AdaptiveDerivatives/1', 'autorouting', 'on');
add_line(model, 'StateIntegrator/1', 'AdaptiveDerivatives/2', 'autorouting', 'on');
add_line(model, 'AdaptiveDerivatives/1', 'StateIntegrator/1', 'autorouting', 'on');
add_line(model, 'StateIntegrator/1', 'StateLog/1', 'autorouting', 'on');
add_line(model, 'StateIntegrator/1', 'StateScope/1', 'autorouting', 'on');

rt = sfroot;
chart = rt.find('-isa', 'Stateflow.EMChart', ...
                'Path', [model '/AdaptiveDerivatives']);

chart.Script = sprintf([ ...
'function dz = fcn(t,z)\n' ...
'%%#codegen\n' ...
'a0 = 1.6; a1 = 1.1; b = 1.4;\n' ...
'wn = 2.0; zeta = 0.8;\n' ...
'P = [1.18125 0.125; 0.125 0.1953125];\n' ...
'B = [0;1]; Gamma = diag([10 10 5]);\n' ...
'x = z(1:2); xm = z(3:4); thetaHat = z(5:7);\n' ...
'r = 0.8*sin(0.5*t) + 0.5*sin(1.3*t) + 0.3*sin(2.1*t);\n' ...
'phi = [x(1);x(2);r];\n' ...
'u = thetaHat.''*phi;\n' ...
'xdot = [x(2); -a0*x(1)-a1*x(2)+b*u];\n' ...
'xmdot = [xm(2); -wn^2*xm(1)-2*zeta*wn*xm(2)+wn^2*r];\n' ...
'e = x-xm;\n' ...
's = e.''*P*B;\n' ...
'thetaHatDot = -Gamma*phi*s*sign(b);\n' ...
'dz = [xdot;xmdot;thetaHatDot];\n' ...
'end\n']);

set_param(model, ...
    'StopTime', '40', ...
    'SolverType', 'Fixed-step', ...
    'Solver', 'ode4', ...
    'FixedStep', '0.001', ...
    'SaveTime', 'on', ...
    'TimeSaveName', 'tout');

save_system(model);
simOut = sim(model);

z = simOut.get('zLog');
tout = z.time;
states = z.signals.values;

figure;
plot(tout, states(:,1), tout, states(:,3), '--', 'LineWidth', 1.1);
grid on;
xlabel('Time (s)');
ylabel('Position-like state');
legend('Plant x_1','Model x_{m1}','Location','best');
title('Simulink direct STR tracking');

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

The textual notebook uses NDSolveValue and reconstructs the composite Lyapunov function.

Chapter13_Lesson3.nb

Notebook[{
Cell["Chapter 13, Lesson 3: Stability Analysis via Lyapunov Methods", "Title"],
Cell["Chapter13_Lesson3.nb", "Subtitle"],
Cell["The controller adapts gains directly. The ideal gains appear only in the validation calculation of V.", "Text"],
Cell[BoxData[ToBoxes[
ClearAll["Global`*"];
a0 = 1.6; a1 = 1.1; b = 1.4; wn = 2.0; zeta = 0.8;
P = { {1.18125, 0.125}, {0.125, 0.1953125} };
B0 = {0, 1};
Gamma = DiagonalMatrix[{10, 10, 5}];
thetaStar = {(a0 - wn^2)/b, (a1 - 2 zeta wn)/b, wn^2/b};
r[t_] := 0.8 Sin[0.5 t] + 0.5 Sin[1.3 t] + 0.3 Sin[2.1 t];
]], "Input"],
Cell[BoxData[ToBoxes[
sol = NDSolveValue[{
 x1'[t] == x2[t],
 x2'[t] == -a0 x1[t] - a1 x2[t] +
   b (th1[t] x1[t] + th2[t] x2[t] + thr[t] r[t]),
 xm1'[t] == xm2[t],
 xm2'[t] == -wn^2 xm1[t] - 2 zeta wn xm2[t] + wn^2 r[t],
 th1'[t] == -10 x1[t] ({x1[t] - xm1[t], x2[t] - xm2[t]}.P.B0) Sign[b],
 th2'[t] == -10 x2[t] ({x1[t] - xm1[t], x2[t] - xm2[t]}.P.B0) Sign[b],
 thr'[t] == -5 r[t] ({x1[t] - xm1[t], x2[t] - xm2[t]}.P.B0) Sign[b],
 x1[0] == 0.4, x2[0] == -0.2, xm1[0] == 0, xm2[0] == 0,
 th1[0] == 0, th2[0] == 0, thr[0] == 0
 },
 {x1, x2, xm1, xm2, th1, th2, thr}, {t, 0, 40},
 Method -> {"EquationSimplification" -> "Residual"}];
]], "Input"],
Cell[BoxData[ToBoxes[
Plot[Evaluate[{sol[[1]][t], sol[[3]][t]}], {t, 0, 40},
 PlotLegends -> {"Plant x1", "Model xm1"}, AxesLabel -> {"t", "state"}]
]], "Input"],
Cell[BoxData[ToBoxes[
Plot[Evaluate[{sol[[1]][t] - sol[[3]][t],
 sol[[2]][t] - sol[[4]][t]}], {t, 0, 40},
 PlotLegends -> {"e1", "e2"}, AxesLabel -> {"t", "error"}]
]], "Input"],
Cell[BoxData[ToBoxes[
Plot[Evaluate[{sol[[5]][t], sol[[6]][t], sol[[7]][t]}], {t, 0, 40},
 PlotLegends -> {"theta1Hat", "theta2Hat", "thetaRHat"},
 AxesLabel -> {"t", "adaptive gain"}]
]], "Input"],
Cell[BoxData[ToBoxes[
v[t_] := With[
 {e = {sol[[1]][t] - sol[[3]][t], sol[[2]][t] - sol[[4]][t]},
  tt = {sol[[5]][t], sol[[6]][t], sol[[7]][t]} - thetaStar},
 e.P.e + Abs[b] tt.Inverse[Gamma].tt];
Plot[v[t], {t, 0, 40}, AxesLabel -> {"t", "V"}]
]], "Input"],
Cell[BoxData[ToBoxes[
{"Ideal gains", thetaStar,
 "Final adaptive gains", {sol[[5]][40], sol[[6]][40], sol[[7]][40]} }
]], "Input"]
},
WindowSize -> {1200, 800},
StyleDefinitions -> "Default.nb"
]

Problem 1 (Lyapunov equation). For \(A_m=\begin{bmatrix}0&1\\-4&-3.2\end{bmatrix}\) and \(Q=I\), derive \(P\) and verify positive definiteness.

Solution.

\[ A_m^TP+PA_m= \begin{bmatrix} -8p_{12}&p_{11}-3.2p_{12}-4p_{22}\\ p_{11}-3.2p_{12}-4p_{22}&2p_{12}-6.4p_{22} \end{bmatrix}. \]

Equating this to \(-I\) gives \(p_{12}=0.125\), \(p_{22}=0.1953125\), and \(p_{11}=1.18125\). Moreover,

\[ p_{11}>0,\qquad \det(P)=0.215087890625>0, \]

so \(P\succ0\) by Sylvester's criterion.

Problem 2 (cross-term cancellation). Derive an adaptive law that makes \(\dot{V}=-\mathbf{e}^TQ\mathbf{e}\).

Solution.

The indefinite term is

\[ 2b\mathbf{e}^TPB \widetilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Choose

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma\boldsymbol{\phi} \mathbf{e}^TPB\operatorname{sgn}(b). \]

Because \(|b|\operatorname{sgn}(b)=b\), the parameter-energy derivative is the exact negative of the mixed term.

Problem 3 (gain convergence). Why does asymptotic tracking not automatically imply convergence to the ideal gain vector?

Solution.

The derivative contains no negative-definite gain-error term. When tracking error becomes small, adaptation stops even if the gain vector is nonideal. Multiple gain vectors can be indistinguishable along a nonexciting trajectory. Persistent excitation supplies the missing identifiability.

Problem 4 (unknown input-gain sign). Explain why the proof fails if \(\operatorname{sgn}(b)\) is unknown.

Solution.

The cancellation requires the correct sign. With the wrong sign, the parameter-energy term reinforces the mixed term, so \(\dot{V}\) need not be nonpositive. This is a structural uncertainty, not a tuning problem.

Problem 5 (Barbalat prerequisites). Show why \(\dot{\mathbf{e} }\) is bounded.

Solution.

Nonincrease of \(V\) bounds \(\mathbf{e}\) and \(\hat{\boldsymbol{\theta} }\). Stability of the reference model and bounded \(r\) bound \(\mathbf{x}_m\). Therefore \(\mathbf{x}=\mathbf{e}+\mathbf{x}_m\) and \(\boldsymbol{\phi}\) are bounded. Every term in

\[ \dot{\mathbf{e} } =A_m\mathbf{e} +Bb\widetilde{\boldsymbol{\theta} }^T\boldsymbol{\phi} \]

is bounded, hence \(\dot{\mathbf{e} }\in L_\infty\). This establishes uniform continuity for Barbalat's lemma.

Matching reduces the direct STR tracking dynamics to a stable reference-model term plus a bilinear gain-error term. The composite Lyapunov function

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

leads to

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma\boldsymbol{\phi} \mathbf{e}^TPB\operatorname{sgn}(b), \qquad \dot{V}=-\mathbf{e}^TQ\mathbf{e}\le0. \]

The result proves boundedness and asymptotic tracking. Parameter convergence is distinct and generally requires persistent excitation. Lesson 4 compares this direct architecture with indirect STRs.