After completing this lesson, students should be able to:

• formulate nonlinear tracking as a combined state-and-parameter error problem;
• construct a filtered error that reduces a second-order tracking problem to a scalar stabilization problem;
• derive a certainty-equivalent controller for a linearly parameterized nonlinear plant;
• derive the adaptive law by cancelling the parameter-error cross term in the Lyapunov derivative;
• prove boundedness of all closed-loop signals and asymptotic convergence of the tracking error;
• explain why successful tracking does not necessarily imply convergence of parameter estimates;
• translate the continuous-time design into numerically executable software.

The lesson uses Lyapunov stability, Barbalat's lemma, linear-in-the-parameters models, projection/normalization concepts, and persistent excitation, all of which were introduced in earlier chapters. Adaptive backstepping is not required for the principal derivation; the design is a one-step nonlinear adaptive cancellation argument.

Consider the second-order nonlinear plant

\[ \dot{x}_1=x_2, \qquad \dot{x}_2=\boldsymbol{\theta}^{T}\boldsymbol{\phi}(x_1,x_2)+b u, \]

where:

• \(x_1\) is the controlled output;
• \(x_2\) is its velocity-like state;
• \(u\) is the control input;
• \(b\) is known and satisfies \(b>0\);
• \(\boldsymbol{\theta}\in\mathbb{R}^{p}\) is an unknown constant parameter vector;
• \(\boldsymbol{\phi}:\mathbb{R}^{2}→\mathbb{R}^{p}\) is a known nonlinear regressor.

The plant is nonlinear because \(\boldsymbol{\phi}(x_1,x_2)\) can contain nonlinear functions such as sine, products, powers, or smooth friction models. It is nevertheless linear in the unknown parameters because the uncertainty appears as \(\boldsymbol{\theta}^{T}\boldsymbol{\phi}\).

Let the desired trajectory \(y_d(t)\) have bounded \(y_d\), \(\dot y_d\), and \(\ddot y_d\). Define

\[ e_1=x_1-y_d, \qquad e_2=x_2-\dot y_d. \]

The control objective is

\[ \lim_{t→\infty}e_1(t)=0, \qquad \lim_{t→\infty}e_2(t)=0, \]

while keeping every internal signal bounded.

flowchart TD
  R["Desired trajectory: yd, yd_dot, yd_ddot"] --> E["Tracking errors: e1, e2"]
  X["Measured states: x1, x2"] --> E
  E --> S["Filtered error: s = e2 + lambda e1"]
  X --> P["Known nonlinear regressor: phi(x)"]
  S --> C["Adaptive controller"]
  P --> C
  H["Parameter estimates: hat_theta"] --> C
  C --> U["Control input u"]
  U --> N["Uncertain nonlinear plant"]
  N --> X
  S --> A["Update law"]
  P --> A
  A --> H
        

Choose a design constant \(\lambda>0\) and define

\[ s=e_2+\lambda e_1. \]

Because \(e_2=\dot e_1\), this relation can be rearranged as

\[ \dot e_1=-\lambda e_1+s. \]

Thus, if \(s\) is bounded and converges to zero, then \(e_1\) is the output of the exponentially stable filter \(1/(p+\lambda)\) driven by a vanishing input. Consequently, \(e_1→0\). Since \(e_2=s-\lambda e_1\), it follows that \(e_2→0\) as well.

Differentiating the filtered error gives

\[ \begin{aligned} \dot s &=\dot e_2+\lambda\dot e_1\\ &=\boldsymbol{\theta}^{T}\boldsymbol{\phi}(x_1,x_2)+bu -\ddot y_d+\lambda e_2. \end{aligned} \]

The adaptive design will force this equation to behave like a stable first-order system despite the unknown parameter vector.

Let \(\hat{\boldsymbol{\theta} }\) denote the online estimate of \(\boldsymbol{\theta}\). Choose \(k>0\) and define the control law

\[ u=\frac{1}{b}\left[ \ddot y_d-\lambda e_2-k s -\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}(x_1,x_2) \right]. \]

The term \(-\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\) attempts to cancel the unknown nonlinear drift, while \(-ks\) provides stabilizing feedback. Define the parameter error using the convention

\[ \tilde{\boldsymbol{\theta} } =\boldsymbol{\theta}-\hat{\boldsymbol{\theta} }. \]

Substitution into the filtered-error dynamics yields

\[ \dot s=-ks+ \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}(x_1,x_2). \]

This equation exposes the only obstacle to a strictly stable error model: the cross term involving the unknown parameter error. The update law is selected specifically to cancel that term in the Lyapunov derivative.

Let \(\boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0\) be the adaptation-gain matrix. Consider

\[ V(s,\tilde{\boldsymbol{\theta} }) =\frac{1}{2}s^{2} +\frac{1}{2}\tilde{\boldsymbol{\theta} }^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta} }. \]

Because the true parameters are constant, \(\dot{\tilde{\boldsymbol{\theta} } } =-\dot{\hat{\boldsymbol{\theta} } }\). Therefore,

\[ \begin{aligned} \dot V &=s\dot s+ \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\tilde{\boldsymbol{\theta} } }\\ &=-ks^{2} +s\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} -\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \end{aligned} \]

Choose the gradient update law

\[ \boxed{ \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}(x_1,x_2)s } \]

Then

\[ \begin{aligned} \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } } &=\tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}s,\\ \dot V&=-ks^{2}\le0. \end{aligned} \]

The adaptive law is therefore not guessed independently. It is synthesized from the Lyapunov derivative so that the uncertain cross term disappears.

flowchart TD
  M["Write filtered-error model: s_dot = -k s + tilde_theta^T phi"] --> V["Choose V = 0.5 s^2 + 0.5 tilde_theta^T Gamma_inv tilde_theta"]
  V --> D["Differentiate V"]
  D --> X["Identify cross term: s tilde_theta^T phi"]
  X --> U["Choose hat_theta_dot = Gamma phi s"]
  U --> C["Cross term cancels"]
  C --> N["Obtain V_dot = -k s^2"]
  N --> B["Boundedness and s in L2"]
  B --> A["Barbalat argument: s tends to zero"]
  A --> T["Stable error filter gives e1 and e2 tend to zero"]
        

Theorem. Suppose that:

1. \(b>0\) is known;
2. \(\boldsymbol{\theta}\) is constant and finite;
3. \(\boldsymbol{\phi}(x_1,x_2)\) is locally Lipschitz;
4. \(y_d,\dot y_d,\ddot y_d\) are bounded;
5. \(k>0\), \(\lambda>0\), and \(\boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0\).

Under the control and adaptation laws derived above, all closed-loop signals are bounded and

\[ e_1(t)→0, \qquad e_2(t)→0, \qquad s(t)→0 \quad \text{as } t→\infty. \]

Proof, Step 1: Positive definiteness and monotonicity.

The Lyapunov function is positive definite in \((s,\tilde{\boldsymbol{\theta} })\), and \(\dot V=-ks^2\le0\). Hence \(V(t)\le V(0)\) for every \(t\ge0\). Therefore, \(s\in L_{\infty}\) and \(\tilde{\boldsymbol{\theta} }\in L_{\infty}\).

Proof, Step 2: Square integrability.

Integrating the Lyapunov derivative gives

\[ V(t)-V(0)=-k\int_{0}^{t}s^{2}(\tau)\,d\tau. \]

Since \(V(t)\ge0\),

\[ \int_{0}^{\infty}s^{2}(\tau)\,d\tau \le\frac{V(0)}{k}<\infty, \]

so \(s\in L_2\).

Proof, Step 3: Bounded tracking errors.

The stable filter \(\dot e_1=-\lambda e_1+s\) has bounded input \(s\). Therefore \(e_1\in L_{\infty}\). Because \(e_2=s-\lambda e_1\), also \(e_2\in L_{\infty}\). Bounded desired signals then imply bounded \(x_1\) and \(x_2\).

Proof, Step 4: Bounded derivative of the filtered error.

Local Lipschitz continuity of the regressor and bounded states imply bounded \(\boldsymbol{\phi}(x_1,x_2)\). Since \(\tilde{\boldsymbol{\theta} }\) is bounded,

\[ \dot s=-ks+ \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \]

is bounded. Thus \(s\) is uniformly continuous.

Proof, Step 5: Asymptotic convergence.

Since \(s\in L_2\cap L_{\infty}\) and \(\dot s\in L_{\infty}\), Barbalat's lemma gives \(s(t)→0\). The stable-filter relation then gives \(e_1(t)→0\), and \(e_2=s-\lambda e_1→0\). This proves asymptotic tracking. \(\square\)

The Lyapunov argument proves boundedness of \(\tilde{\boldsymbol{\theta} }\), not necessarily \(\tilde{\boldsymbol{\theta} }→\mathbf{0}\). When the tracking error becomes small, the update law

\[ \dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}s \]

also becomes small. The parameter estimate can therefore stop at any value that is compatible with the observed trajectory. Exact identification generally requires 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 \ge\alpha\mathbf{I}, \qquad \forall t\ge0. \]

Without this richness condition, asymptotic tracking may be excellent even though one or more estimates remain different from their true values. This distinction is central when interpreting simulation plots.

Consider the relative-degree-\(n\) model

\[ y^{(n)}= \boldsymbol{\theta}^{T}\boldsymbol{\phi}(\mathbf{x})+bu. \]

Let \(e=y-y_d\) and define a filtered error using a Hurwitz polynomial:

\[ s=e^{(n-1)}+a_{n-2}e^{(n-2)}+\cdots+a_1\dot e+a_0 e, \]

where \(p^{n-1}+a_{n-2}p^{n-2}+\cdots+a_0\) is Hurwitz. Differentiation gives

\[ \dot s= \boldsymbol{\theta}^{T}\boldsymbol{\phi}+bu-y_d^{(n)} +a_{n-2}e^{(n-1)}+\cdots+a_0\dot e. \]

The direct extension is

\[ \begin{aligned} u&=\frac{1}{b}\left[ y_d^{(n)}-a_{n-2}e^{(n-1)}-\cdots-a_0\dot e -ks-\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \right],\\ \dot{\hat{\boldsymbol{\theta} } } &=\boldsymbol{\Gamma}\boldsymbol{\phi}s. \end{aligned} \]

The same Lyapunov function gives \(\dot V=-ks^2\), while the Hurwitz error filter converts \(s→0\) into convergence of the complete tracking-error state. This extension requires the uncertainty to enter through the same channel as the control input; unmatched uncertainty needs a different design.

All implementations use

\[ \begin{aligned} \dot x_1&=x_2,\\ \dot x_2&=\theta_1\sin(x_1) +\theta_2 x_2|x_2|+bu, \end{aligned} \]

with

\[ \boldsymbol{\theta}=\begin{bmatrix}1.2&-0.15\end{bmatrix}^{T}, \qquad b=1.2, \qquad \boldsymbol{\phi}=\begin{bmatrix}\sin x_1&x_2|x_2|\end{bmatrix}^{T}. \]

The desired trajectory is

\[ y_d(t)=0.8\sin(0.6t)+0.3\sin(1.3t), \]

\[ \dot y_d(t)=0.48\cos(0.6t)+0.39\cos(1.3t), \]

\[ \ddot y_d(t)=-0.288\sin(0.6t)-0.507\sin(1.3t). \]

The gains and initial conditions are

\[ \lambda=2, \quad k=4, \quad \boldsymbol{\Gamma}=\operatorname{diag}(3,1), \quad \mathbf{x}(0)=\begin{bmatrix}1.4&-0.5\end{bmatrix}^{T}, \quad \hat{\boldsymbol{\theta} }(0)=\mathbf{0}. \]

With an RK4 step of \(0.002\,\mathrm{s}\), the verified Python, C++, and Java implementations produce a tail RMS tracking error of approximately \(8.03\times10^{-4}\). The second parameter estimate remains offset from its true value, illustrating that tracking convergence does not automatically imply parameter convergence.

Chapter18_Lesson3.py

This implementation uses NumPy for vector operations, a from-scratch fourth-order Runge–Kutta integrator, Matplotlib for plots, and the Python standard library for CSV export.

"""Chapter18_Lesson3.py
Lyapunov-based adaptive tracking for a linearly parameterized nonlinear plant.

Plant:
    x1_dot = x2
    x2_dot = theta1*sin(x1) + theta2*x2*abs(x2) + b*u

Controller:
    s = (x2 - yd_dot) + lambda*(x1 - yd)
    u = (yd_ddot - lambda*(x2 - yd_dot) - k*s - hat_theta.T*phi)/b
    hat_theta_dot = Gamma*phi*s
"""

from __future__ import annotations

import csv
from dataclasses import dataclass, field
from pathlib import Path
from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Settings:
    theta_true: np.ndarray
    b: float = 1.2
    lam: float = 2.0
    k: float = 4.0
    gamma: np.ndarray = field(default_factory=lambda: np.array([3.0, 1.0], dtype=float))
    dt: float = 0.002
    t_final: float = 30.0


def reference(t: float) -> Tuple[float, float, float]:
    yd = 0.8 * np.sin(0.6 * t) + 0.3 * np.sin(1.3 * t)
    yd_dot = 0.48 * np.cos(0.6 * t) + 0.39 * np.cos(1.3 * t)
    yd_ddot = -0.288 * np.sin(0.6 * t) - 0.507 * np.sin(1.3 * t)
    return float(yd), float(yd_dot), float(yd_ddot)


def closed_loop_rhs(t: float, z: np.ndarray, cfg: Settings) -> Tuple[np.ndarray, float, float, float]:
    x1, x2, hat1, hat2 = z
    yd, yd_dot, yd_ddot = reference(t)

    e1 = x1 - yd
    e2 = x2 - yd_dot
    s = e2 + cfg.lam * e1
    phi = np.array([np.sin(x1), x2 * abs(x2)], dtype=float)
    hat_theta = np.array([hat1, hat2], dtype=float)

    u = (yd_ddot - cfg.lam * e2 - cfg.k * s - float(hat_theta @ phi)) / cfg.b
    x1_dot = x2
    x2_dot = float(cfg.theta_true @ phi) + cfg.b * u
    hat_dot = cfg.gamma * phi * s

    dz = np.array([x1_dot, x2_dot, hat_dot[0], hat_dot[1]], dtype=float)
    return dz, u, yd, s


def rk4_step(t: float, z: np.ndarray, h: float, cfg: Settings) -> np.ndarray:
    k1 = closed_loop_rhs(t, z, cfg)[0]
    k2 = closed_loop_rhs(t + 0.5 * h, z + 0.5 * h * k1, cfg)[0]
    k3 = closed_loop_rhs(t + 0.5 * h, z + 0.5 * h * k2, cfg)[0]
    k4 = closed_loop_rhs(t + h, z + h * k3, cfg)[0]
    return z + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(cfg: Settings) -> dict[str, np.ndarray]:
    steps = int(round(cfg.t_final / cfg.dt)) + 1
    time = np.linspace(0.0, cfg.t_final, steps)
    state = np.zeros((steps, 4), dtype=float)
    state[0] = np.array([1.4, -0.5, 0.0, 0.0], dtype=float)

    control = np.zeros(steps)
    desired = np.zeros(steps)
    filtered_error = np.zeros(steps)

    for i, t in enumerate(time):
        _, control[i], desired[i], filtered_error[i] = closed_loop_rhs(t, state[i], cfg)
        if i + 1 < steps:
            state[i + 1] = rk4_step(t, state[i], cfg.dt, cfg)

    tracking_error = state[:, 0] - desired
    return {
        "time": time,
        "state": state,
        "control": control,
        "desired": desired,
        "tracking_error": tracking_error,
        "filtered_error": filtered_error,
    }


def save_results(data: dict[str, np.ndarray], output_dir: Path) -> None:
    output_dir.mkdir(parents=True, exist_ok=True)
    csv_path = output_dir / "Chapter18_Lesson3_results.csv"
    with csv_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(["t", "x1", "x2", "yd", "e1", "s", "hat_theta1", "hat_theta2", "u"])
        state = data["state"]
        for i, t in enumerate(data["time"]):
            writer.writerow([
                t,
                state[i, 0],
                state[i, 1],
                data["desired"][i],
                data["tracking_error"][i],
                data["filtered_error"][i],
                state[i, 2],
                state[i, 3],
                data["control"][i],
            ])

    plt.figure(figsize=(9, 5))
    plt.plot(data["time"], data["state"][:, 0], label="x1")
    plt.plot(data["time"], data["desired"], "--", label="yd")
    plt.xlabel("Time (s)")
    plt.ylabel("Position")
    plt.title("Adaptive nonlinear tracking")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter18_Lesson3_tracking.png", dpi=180)
    plt.close()

    plt.figure(figsize=(9, 5))
    plt.plot(data["time"], data["state"][:, 2], label="hat theta1")
    plt.plot(data["time"], data["state"][:, 3], label="hat theta2")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.title("Adaptive parameter estimates")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter18_Lesson3_parameters.png", dpi=180)
    plt.close()


def main() -> None:
    cfg = Settings(theta_true=np.array([1.2, -0.15], dtype=float))
    data = simulate(cfg)
    save_results(data, Path("."))

    e = data["tracking_error"]
    tail = e[int(0.8 * len(e)) :]
    print(f"Maximum absolute tracking error: {np.max(np.abs(e)):.6f}")
    print(f"Tail RMS tracking error: {np.sqrt(np.mean(tail**2)):.6e}")
    print("Final parameter estimates:", data["state"][-1, 2:])
    print("True parameters:", cfg.theta_true)
    print("Note: tracking can converge without parameter convergence unless the regressor is persistently exciting.")


if __name__ == "__main__":
    main()

Chapter18_Lesson3.cpp

This C++17 implementation uses only the standard library, performs RK4 integration, and writes the complete trajectory to a CSV file.

// Chapter18_Lesson3.cpp
// Lyapunov-based adaptive tracking for a linearly parameterized nonlinear plant.

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

struct Settings {
    std::array<double, 2> thetaTrue{1.2, -0.15};
    std::array<double, 2> gamma{3.0, 1.0};
    double b{1.2};
    double lambda{2.0};
    double k{4.0};
    double dt{0.002};
    double tFinal{30.0};
};

struct RefSignal {
    double y;
    double yd;
    double ydd;
};

struct Eval {
    std::array<double, 4> dz{};
    double u{};
    double desired{};
    double s{};
};

RefSignal reference(double t) {
    return {
        0.8 * std::sin(0.6 * t) + 0.3 * std::sin(1.3 * t),
        0.48 * std::cos(0.6 * t) + 0.39 * std::cos(1.3 * t),
        -0.288 * std::sin(0.6 * t) - 0.507 * std::sin(1.3 * t)
    };
}

Eval evaluate(double t, const std::array<double, 4>& z, const Settings& cfg) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double hat1 = z[2];
    const double hat2 = z[3];
    const RefSignal ref = reference(t);

    const double e1 = x1 - ref.y;
    const double e2 = x2 - ref.yd;
    const double s = e2 + cfg.lambda * e1;
    const double phi1 = std::sin(x1);
    const double phi2 = x2 * std::abs(x2);

    const double u = (ref.ydd - cfg.lambda * e2 - cfg.k * s - hat1 * phi1 - hat2 * phi2) / cfg.b;
    const double x1Dot = x2;
    const double x2Dot = cfg.thetaTrue[0] * phi1 + cfg.thetaTrue[1] * phi2 + cfg.b * u;
    const double hat1Dot = cfg.gamma[0] * phi1 * s;
    const double hat2Dot = cfg.gamma[1] * phi2 * s;

    return { { {x1Dot, x2Dot, hat1Dot, hat2Dot} }, u, ref.y, s};
}

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

std::array<double, 4> rk4Step(double t, const std::array<double, 4>& z, const Settings& cfg) {
    const auto k1 = evaluate(t, z, cfg).dz;
    const auto k2 = evaluate(t + 0.5 * cfg.dt, addScaled(z, k1, 0.5 * cfg.dt), cfg).dz;
    const auto k3 = evaluate(t + 0.5 * cfg.dt, addScaled(z, k2, 0.5 * cfg.dt), cfg).dz;
    const auto k4 = evaluate(t + cfg.dt, addScaled(z, k3, cfg.dt), cfg).dz;

    std::array<double, 4> next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (cfg.dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
    }
    return next;
}

int main() {
    try {
        const Settings cfg;
        std::array<double, 4> z{1.4, -0.5, 0.0, 0.0};
        const int steps = static_cast<int>(std::llround(cfg.tFinal / cfg.dt)) + 1;

        std::ofstream csv("Chapter18_Lesson3_cpp_results.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open output CSV file.");
        }
        csv << std::setprecision(12);
        csv << "t,x1,x2,yd,e1,s,hat_theta1,hat_theta2,u\n";

        double tailSquaredError = 0.0;
        int tailCount = 0;
        double maxAbsError = 0.0;
        const int tailStart = static_cast<int>(0.8 * steps);

        for (int i = 0; i < steps; ++i) {
            const double t = i * cfg.dt;
            const Eval ev = evaluate(t, z, cfg);
            const double e1 = z[0] - ev.desired;
            maxAbsError = std::max(maxAbsError, std::abs(e1));
            if (i >= tailStart) {
                tailSquaredError += e1 * e1;
                ++tailCount;
            }

            csv << t << ',' << z[0] << ',' << z[1] << ',' << ev.desired << ','
                << e1 << ',' << ev.s << ',' << z[2] << ',' << z[3] << ',' << ev.u << '\n';

            if (i + 1 < steps) {
                z = rk4Step(t, z, cfg);
            }
        }

        const double tailRms = std::sqrt(tailSquaredError / static_cast<double>(tailCount));
        std::cout << "Maximum absolute tracking error: " << maxAbsError << '\n';
        std::cout << "Tail RMS tracking error: " << tailRms << '\n';
        std::cout << "Final parameter estimates: [" << z[2] << ", " << z[3] << "]\n";
        std::cout << "True parameters: [" << cfg.thetaTrue[0] << ", " << cfg.thetaTrue[1] << "]\n";
        std::cout << "Tracking convergence does not by itself imply parameter convergence without persistent excitation.\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Chapter18_Lesson3.java

This Java implementation uses records for compact data structures, standard-library file I/O, and the same RK4 integration and adaptive law as the other versions.

// Chapter18_Lesson3.java
// Lyapunov-based adaptive tracking for a linearly parameterized nonlinear plant.

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 Chapter18_Lesson3 {
    private static final double[] THETA_TRUE = {1.2, -0.15};
    private static final double[] GAMMA = {3.0, 1.0};
    private static final double B = 1.2;
    private static final double LAMBDA = 2.0;
    private static final double K = 4.0;
    private static final double DT = 0.002;
    private static final double T_FINAL = 30.0;

    private Chapter18_Lesson3() {
    }

    private record Reference(double y, double yd, double ydd) {
    }

    private record Evaluation(double[] dz, double u, double desired, double s) {
    }

    private static Reference reference(double t) {
        double y = 0.8 * Math.sin(0.6 * t) + 0.3 * Math.sin(1.3 * t);
        double yd = 0.48 * Math.cos(0.6 * t) + 0.39 * Math.cos(1.3 * t);
        double ydd = -0.288 * Math.sin(0.6 * t) - 0.507 * Math.sin(1.3 * t);
        return new Reference(y, yd, ydd);
    }

    private static Evaluation evaluate(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double hat1 = z[2];
        double hat2 = z[3];
        Reference ref = reference(t);

        double e1 = x1 - ref.y();
        double e2 = x2 - ref.yd();
        double s = e2 + LAMBDA * e1;
        double phi1 = Math.sin(x1);
        double phi2 = x2 * Math.abs(x2);

        double u = (ref.ydd() - LAMBDA * e2 - K * s - hat1 * phi1 - hat2 * phi2) / B;
        double x1Dot = x2;
        double x2Dot = THETA_TRUE[0] * phi1 + THETA_TRUE[1] * phi2 + B * u;
        double hat1Dot = GAMMA[0] * phi1 * s;
        double hat2Dot = GAMMA[1] * phi2 * s;

        return new Evaluation(new double[]{x1Dot, x2Dot, hat1Dot, hat2Dot}, u, ref.y(), s);
    }

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

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

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 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);
        double[] z = {1.4, -0.5, 0.0, 0.0};
        int steps = (int) Math.round(T_FINAL / DT) + 1;
        int tailStart = (int) (0.8 * steps);
        double tailSquaredError = 0.0;
        int tailCount = 0;
        double maxAbsError = 0.0;

        Path output = Path.of("Chapter18_Lesson3_java_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("t,x1,x2,yd,e1,s,hat_theta1,hat_theta2,u\n");

            for (int i = 0; i < steps; i++) {
                double t = i * DT;
                Evaluation ev = evaluate(t, z);
                double e1 = z[0] - ev.desired();
                maxAbsError = Math.max(maxAbsError, Math.abs(e1));
                if (i >= tailStart) {
                    tailSquaredError += e1 * e1;
                    tailCount++;
                }

                writer.write(String.format(Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    t, z[0], z[1], ev.desired(), e1, ev.s(), z[2], z[3], ev.u()));

                if (i + 1 < steps) {
                    z = rk4Step(t, z);
                }
            }
        } catch (IOException ex) {
            System.err.println("Unable to write simulation results: " + ex.getMessage());
            System.exit(1);
        }

        double tailRms = Math.sqrt(tailSquaredError / tailCount);
        System.out.printf(Locale.US, "Maximum absolute tracking error: %.6f%n", maxAbsError);
        System.out.printf(Locale.US, "Tail RMS tracking error: %.6e%n", tailRms);
        System.out.printf(Locale.US, "Final parameter estimates: [%.6f, %.6f]%n", z[2], z[3]);
        System.out.printf(Locale.US, "True parameters: [%.6f, %.6f]%n", THETA_TRUE[0], THETA_TRUE[1]);
        System.out.println("Tracking convergence does not by itself imply parameter convergence without persistent excitation.");
    }
}

Chapter18_Lesson3.m

This MATLAB implementation uses ode45, local functions, tables, CSV export, and plots for the tracking and parameter trajectories.

% Chapter18_Lesson3.m
% Lyapunov-based adaptive tracking for a linearly parameterized nonlinear plant.

clear; clc; close all;

p.thetaTrue = [1.2; -0.15];
p.b = 1.2;
p.lambda = 2.0;
p.k = 4.0;
p.Gamma = diag([3.0, 1.0]);
p.tFinal = 30.0;

z0 = [1.4; -0.5; 0.0; 0.0];
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.01);
[t, z] = ode45(@(t, z) closedLoopDynamics(t, z, p), [0 p.tFinal], z0, options);

n = numel(t);
yd = zeros(n, 1);
s = zeros(n, 1);
u = zeros(n, 1);
for i = 1:n
    [~, u(i), yd(i), s(i)] = closedLoopDynamics(t(i), z(i, :).', p);
end

e1 = z(:, 1) - yd;
tail = e1(max(1, floor(0.8*n)):end);

fprintf('Maximum absolute tracking error: %.6f\n', max(abs(e1)));
fprintf('Tail RMS tracking error: %.6e\n', sqrt(mean(tail.^2)));
fprintf('Final parameter estimates: [%.6f, %.6f]\n', z(end, 3), z(end, 4));
fprintf('True parameters: [%.6f, %.6f]\n', p.thetaTrue(1), p.thetaTrue(2));
fprintf(['Tracking convergence does not by itself imply parameter convergence ' ...
         'without persistent excitation.\n']);

results = table(t, z(:,1), z(:,2), yd, e1, s, z(:,3), z(:,4), u, ...
    'VariableNames', {'t','x1','x2','yd','e1','s','hatTheta1','hatTheta2','u'});
writetable(results, 'Chapter18_Lesson3_matlab_results.csv');

figure;
plot(t, z(:,1), 'LineWidth', 1.3); hold on;
plot(t, yd, '--', 'LineWidth', 1.3);
grid on; xlabel('Time (s)'); ylabel('Position');
legend('x_1', 'y_d', 'Location', 'best');
title('Lyapunov-based adaptive nonlinear tracking');

figure;
plot(t, z(:,3), 'LineWidth', 1.3); hold on;
plot(t, z(:,4), 'LineWidth', 1.3);
yline(p.thetaTrue(1), '--');
yline(p.thetaTrue(2), '--');
grid on; xlabel('Time (s)'); ylabel('Parameter estimate');
legend('\hat{\theta}_1', '\hat{\theta}_2', '\theta_1', '\theta_2', 'Location', 'best');
title('Adaptive parameter estimates');

function [dz, u, yd, s] = closedLoopDynamics(t, z, p)
    x1 = z(1);
    x2 = z(2);
    hatTheta = z(3:4);

    [yd, ydDot, ydDDot] = referenceSignal(t);
    e1 = x1 - yd;
    e2 = x2 - ydDot;
    s = e2 + p.lambda*e1;
    phi = [sin(x1); x2*abs(x2)];

    u = (ydDDot - p.lambda*e2 - p.k*s - hatTheta.'*phi)/p.b;
    x1Dot = x2;
    x2Dot = p.thetaTrue.'*phi + p.b*u;
    hatThetaDot = p.Gamma*phi*s;

    dz = [x1Dot; x2Dot; hatThetaDot];
end

function [yd, ydDot, ydDDot] = referenceSignal(t)
    yd = 0.8*sin(0.6*t) + 0.3*sin(1.3*t);
    ydDot = 0.48*cos(0.6*t) + 0.39*cos(1.3*t);
    ydDDot = -0.288*sin(0.6*t) - 0.507*sin(1.3*t);
end

Chapter18_Lesson3_Simulink.m

This script programmatically creates a Simulink model containing continuous integrators and a MATLAB Function block, runs the model, and plots the tracking result. It requires Simulink and Stateflow.

% Chapter18_Lesson3_Simulink.m
% Programmatically builds a Simulink model for the adaptive nonlinear tracker.
% Requires Simulink and Stateflow (for the MATLAB Function block API).

clear; clc;
model = 'Chapter18_Lesson3_Model';

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

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

add_block('simulink/Sources/Clock', [model '/Clock'], ...
    'Position', [30 45 60 75]);

stateNames = {'x1','x2','hatTheta1','hatTheta2'};
initialConditions = {'1.4','-0.5','0','0'};
for i = 1:4
    y = 100 + 75*(i-1);
    add_block('simulink/Continuous/Integrator', [model '/' stateNames{i}], ...
        'InitialCondition', initialConditions{i}, ...
        'Position', [430 y 460 y+30]);
end

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/AdaptiveDynamics'], 'Position', [150 70 350 390]);
chart = find(sfroot, '-isa', 'Stateflow.EMChart', ...
    'Path', [model '/AdaptiveDynamics']);
chart.Script = sprintf([ ...
'function [dx1,dx2,dh1,dh2,u,yd,s] = fcn(t,x1,x2,h1,h2)\n' ...
'%%#codegen\n' ...
'b = 1.2; lambda = 2.0; k = 4.0;\n' ...
'theta1 = 1.2; theta2 = -0.15; gamma1 = 3.0; gamma2 = 1.0;\n' ...
'yd = 0.8*sin(0.6*t) + 0.3*sin(1.3*t);\n' ...
'ydDot = 0.48*cos(0.6*t) + 0.39*cos(1.3*t);\n' ...
'ydDDot = -0.288*sin(0.6*t) - 0.507*sin(1.3*t);\n' ...
'e1 = x1 - yd; e2 = x2 - ydDot; s = e2 + lambda*e1;\n' ...
'phi1 = sin(x1); phi2 = x2*abs(x2);\n' ...
'u = (ydDDot - lambda*e2 - k*s - h1*phi1 - h2*phi2)/b;\n' ...
'dx1 = x2;\n' ...
'dx2 = theta1*phi1 + theta2*phi2 + b*u;\n' ...
'dh1 = gamma1*phi1*s;\n' ...
'dh2 = gamma2*phi2*s;\n' ...
'end\n']);

add_block('simulink/Signal Routing/Mux', [model '/LogMux'], ...
    'Inputs', '7', 'Position', [540 105 545 315]);
add_block('simulink/Sinks/To Workspace', [model '/ToWorkspace'], ...
    'VariableName', 'adaptiveLog', 'SaveFormat', 'StructureWithTime', ...
    'Position', [600 185 700 215]);

add_line(model, 'Clock/1', 'AdaptiveDynamics/1');
for i = 1:4
    add_line(model, [stateNames{i} '/1'], ['AdaptiveDynamics/' num2str(i+1)]);
    add_line(model, ['AdaptiveDynamics/' num2str(i)], [stateNames{i} '/1']);
end

% Log x1, x2, hatTheta1, hatTheta2, u, yd, and s.
for i = 1:4
    add_line(model, [stateNames{i} '/1'], ['LogMux/' num2str(i)]);
end
add_line(model, 'AdaptiveDynamics/5', 'LogMux/5');
add_line(model, 'AdaptiveDynamics/6', 'LogMux/6');
add_line(model, 'AdaptiveDynamics/7', 'LogMux/7');
add_line(model, 'LogMux/1', 'ToWorkspace/1');

autoArrange = get_param(model, 'Location'); %#ok<NASGU>
Simulink.BlockDiagram.arrangeSystem(model);
save_system(model);

simOut = sim(model, 'ReturnWorkspaceOutputs', 'on');
logData = simOut.get('adaptiveLog');
t = logData.time;
y = logData.signals.values;

figure;
plot(t, y(:,1), 'LineWidth', 1.3); hold on;
plot(t, y(:,6), '--', 'LineWidth', 1.3);
grid on; xlabel('Time (s)'); ylabel('Position');
legend('x_1', 'y_d', 'Location', 'best');
title('Simulink adaptive nonlinear tracking');

fprintf('Created and simulated %s.slx\n', model);

The generated model contains four continuous integrators for \(x_1,x_2,\hat\theta_1,\hat\theta_2\). A MATLAB Function block evaluates the reference trajectory, regressor, controller, plant derivatives, and adaptive law. This arrangement mirrors the mathematical state equations directly and avoids algebraic loops.

Chapter18_Lesson3.nb

This Wolfram Notebook expression defines the closed-loop differential equations, solves them with NDSolveValue, and plots the tracking and parameter trajectories.

Notebook[{
  Cell[
    "Chapter 18, Lesson 3: Lyapunov-Based Adaptive Design for Nonlinear Tracking",
    "Title"
  ],
  Cell[
    "This notebook simulates the same nonlinear adaptive tracker used in the Python, C++, Java, and MATLAB implementations.",
    "Text"
  ],
  Cell[BoxData["ClearAll[\"Global`*\"];

thetaTrue = {1.2, -0.15};
b = 1.2;
lambda = 2.0;
k = 4.0;
gamma = {3.0, 1.0};

yd[t_] := 0.8 Sin[0.6 t] + 0.3 Sin[1.3 t];
ydDot[t_] := 0.48 Cos[0.6 t] + 0.39 Cos[1.3 t];
ydDDot[t_] := -0.288 Sin[0.6 t] - 0.507 Sin[1.3 t];

phi[q_, v_] := {Sin[q], v Abs[v]};
sErr[t_] := x2[t] - ydDot[t] + lambda (x1[t] - yd[t]);
u[t_] := (ydDDot[t] - lambda (x2[t] - ydDot[t]) - k sErr[t]
  - {h1[t], h2[t]}.phi[x1[t], x2[t]])/b;

eqns = {
  x1'[t] == x2[t],
  x2'[t] == thetaTrue.phi[x1[t], x2[t]] + b u[t],
  h1'[t] == gamma[[1]] phi[x1[t], x2[t]][[1]] sErr[t],
  h2'[t] == gamma[[2]] phi[x1[t], x2[t]][[2]] sErr[t],
  x1[0] == 1.4, x2[0] == -0.5, h1[0] == 0, h2[0] == 0
};

sol = NDSolveValue[
  eqns, {x1, x2, h1, h2}, {t, 0, 30},
  MaxStepFraction -> 1/3000
];

trackingPlot = Plot[
  Evaluate[{sol[[1]][t], yd[t]}], {t, 0, 30},
  PlotLegends -> {\"x1\", \"yd\"}, GridLines -> Automatic,
  PlotLabel -> \"Adaptive nonlinear tracking\"
];

parameterPlot = Plot[
  Evaluate[{sol[[3]][t], sol[[4]][t]}], {t, 0, 30},
  PlotLegends -> {\"hat theta1\", \"hat theta2\"}, GridLines -> Automatic,
  PlotLabel -> \"Adaptive parameter estimates\"
];

Column[{trackingPlot, parameterPlot}]"], "Input"],
  Cell[
    "Tracking may converge even when the parameter estimates do not equal the true parameters; persistent excitation is the additional condition needed for parameter convergence.",
    "Text"
  ]
},
  WindowTitle -> "Chapter18_Lesson3.nb",
  StyleDefinitions -> "Default.nb"
]        

16.1 Adaptation gain

Increasing \(\boldsymbol{\Gamma}\) accelerates parameter motion but also makes estimates more sensitive to discretization error, measurement noise, and high-frequency regressor components. The adaptation gain does not directly increase the negative term in \(\dot V=-ks^2\); that role belongs to \(k\).

16.2 Feedback gain

Larger \(k\) increases error dissipation, but excessive values can demand large control action. The design in this lesson assumes no actuator saturation; saturation will be treated explicitly in Chapter 19.

16.3 Error-filter bandwidth

The constant \(\lambda\) sets the desired decay rate from \(s\) to the position error. A large value can amplify velocity measurement noise and increase the control term \(-\lambda e_2\).

16.4 Integration method

The adaptive law is part of the continuous closed-loop dynamics and must be integrated at the same numerical stages as the plant. Updating parameters only once after a plant step is a different discretization and can degrade stability when the step size is not sufficiently small.

16.5 Robust modifications

The ideal proof assumes exact linearly parameterized dynamics. With noise or unmodeled dynamics, projection, normalization, dead zones, or leakage terms from Chapters 8 and 9 can be added. Those modifications generally change asymptotic convergence into a boundedness or ultimate-boundedness result.

Problem 1 — Sign consistency in the update law. Suppose the parameter error is defined instead as \(\bar{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}\). Derive the filtered-error model and the corresponding adaptive law.

Solution. Since

\[ \boldsymbol{\theta}-\hat{\boldsymbol{\theta} } =-\bar{\boldsymbol{\theta} }, \]

the error dynamics become

\[ \dot s=-ks- \bar{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]

Choose

\[ V=\frac12s^2+\frac12 \bar{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \bar{\boldsymbol{\theta} }. \]

Because \(\dot{\bar{\boldsymbol{\theta} } } =\dot{\hat{\boldsymbol{\theta} } }\),

\[ \dot V=-ks^2-s\bar{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} +\bar{\boldsymbol{\theta} }^{T}\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta} } }. \]

The same physical update law \(\dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}s\) cancels the cross term. The apparent sign depends on the parameter-error convention, but the actual estimator equation remains unchanged.

Problem 2 — Scalar nonlinear plant. Consider \(\dot x=\theta\sin x+u\) with desired trajectory \(x_d(t)\). Derive a first-order adaptive tracking controller.

Solution. Let \(e=x-x_d\). Then

\[ \dot e=\theta\sin x+u-\dot x_d. \]

Choose

\[ u=\dot x_d-ke-\hat\theta\sin x, \qquad \dot{\hat\theta}=\gamma e\sin x, \]

with \(k>0\) and \(\gamma>0\). Then

\[ \dot e=-ke+\tilde\theta\sin x. \]

For

\[ V=\frac12e^2+\frac{1}{2\gamma}\tilde\theta^2, \]

the chosen update law gives \(\dot V=-ke^2\). Under bounded reference signals and a bounded regressor, Barbalat's lemma yields \(e→0\).

Problem 3 — Why negative semidefiniteness is sufficient here. The derivative \(\dot V=-ks^2\) is not negative definite in \((s,\tilde{\boldsymbol{\theta} })\). Explain why tracking still converges and why parameter convergence does not follow.

Solution.

Negative semidefiniteness gives boundedness of both \(s\) and the parameter error and also gives \(s\in L_2\). The closed-loop equations then establish bounded \(\dot s\), so Barbalat's lemma yields \(s→0\). The stable error filter gives tracking convergence. However, the set \(s=0\) contains many possible parameter errors. Without a richness condition on the regressor, the trajectory cannot distinguish all parameter vectors; therefore the invariant set need not require \(\tilde{\boldsymbol{\theta} }=0\).

Problem 4 — Higher-order filtered error. For \(y^{(3)}=\boldsymbol{\theta}^{T}\boldsymbol{\phi}(\mathbf{x})+bu\), choose a filtered error and derive the control and update laws.

Solution. Let \(e=y-y_d\) and choose

\[ s=\ddot e+a_1\dot e+a_0e, \]

where \(p^2+a_1p+a_0\) is Hurwitz; equivalently, \(a_1>0\) and \(a_0>0\). Then

\[ \dot s= \boldsymbol{\theta}^{T}\boldsymbol{\phi}+bu-y_d^{(3)} +a_1\ddot e+a_0\dot e. \]

Choose

\[ \begin{aligned} u&=\frac1b\left[y_d^{(3)}-a_1\ddot e-a_0\dot e -ks-\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\right],\\ \dot{\hat{\boldsymbol{\theta} } } &=\boldsymbol{\Gamma}\boldsymbol{\phi}s. \end{aligned} \]

This gives \(\dot s=-ks+ \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\) and the same Lyapunov proof. Since the polynomial is Hurwitz, convergence of \(s\) implies convergence of \(e,\dot e,\ddot e\).

Problem 5 — Known negative input gain. Suppose the plant has the same form but \(b<0\) is known. Must the adaptive law change?

Solution.

No. The controller still divides by the known signed gain:

\[ u=\frac{1}{b}\left[ \ddot y_d-\lambda e_2-ks -\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\right]. \]

Substitution produces the same filtered-error equation because \(bu\) exactly reconstructs the bracketed term. Therefore the update law remains \(\dot{\hat{\boldsymbol{\theta} } } =\boldsymbol{\Gamma}\boldsymbol{\phi}s\). The difficult case is an unknown sign of \(b\), which requires the special methods studied in Chapter 16.

Problem 6 — Lyapunov performance bound. Derive an upper bound on the accumulated squared filtered error.

Solution. From

\[ V(t)-V(0)=-k\int_0^t s^2(\tau)\,d\tau \]

and \(V(t)\ge0\),

\[ \int_0^t s^2(\tau)\,d\tau \le\frac{V(0)}{k}. \]

Taking \(t→\infty\) gives

\[ \|s\|_{2}^{2}\le\frac{V(0)}{k}. \]

This inequality shows how a larger feedback gain reduces the theoretical upper bound on the total filtered-error energy, although actuator limits and noise prevent arbitrarily large gains in practice.

A nonlinear plant can be handled by classical Lyapunov adaptive control when its uncertainty is linear in a constant parameter vector and enters through the control channel. The filtered error converts tracking into stabilization; certainty-equivalent cancellation exposes a parameter-error cross term; and the parameter update law is selected to cancel that term in the Lyapunov derivative. The resulting identity \(\dot V=-ks^2\) establishes boundedness, square integrability, and—through Barbalat's lemma— asymptotic tracking. Parameter convergence remains a separate question that generally requires persistent excitation. The cross-language implementations demonstrate how the plant and estimator must be integrated as one coupled dynamical system.