After completing this lesson, students should be able to:

• construct recursive state-error coordinates for a strict-feedback plant;
• build a positive-definite Lyapunov function containing state and parameter errors;
• identify undesirable cross terms created at each backstepping step;
• choose virtual controls, the actual control, and an adaptive law that cancel those terms;
• distinguish tracking convergence from parameter convergence; and
• implement and numerically verify an adaptive backstepping controller.

Lessons 1–3 introduced strict-feedback structure, virtual controls, and recursive adaptive laws. The present lesson focuses on the proof mechanism: every design choice is made so that the derivative of one composite Lyapunov function has a controlled sign.

Consider the linearly parameterized strict-feedback system

\[ \begin{aligned} \dot{x}_1 &= x_2 + f_1(x_1)+\boldsymbol{\phi}_1^{\mathsf T}(x_1)\boldsymbol{\theta},\\ \dot{x}_2 &= x_3 + f_2(\bar{\mathbf{x}}_2)+\boldsymbol{\phi}_2^{\mathsf T}(\bar{\mathbf{x}}_2)\boldsymbol{\theta},\\ &\hspace{1.5cm}\vdots\\ \dot{x}_{n-1} &= x_n+f_{n-1}(\bar{\mathbf{x}}_{n-1})+ \boldsymbol{\phi}_{n-1}^{\mathsf T}(\bar{\mathbf{x}}_{n-1})\boldsymbol{\theta},\\ \dot{x}_n &= b(\mathbf{x})u+f_n(\mathbf{x})+ \boldsymbol{\phi}_n^{\mathsf T}(\mathbf{x})\boldsymbol{\theta}, \end{aligned} \]

where \( \bar{\mathbf{x}}_i=[x_1,\ldots,x_i]^{\mathsf T} \), the constant vector \( \boldsymbol{\theta}\in\mathbb{R}^p \) is unknown, the functions \(f_i\) and regressors \(\boldsymbol{\phi}_i\) are known, and the input gain satisfies

\[ 0<b_0\le b(\mathbf{x})\le b_1. \]

The control objective is to make \(x_1(t)\) track a bounded reference \(y_d(t)\) whose required derivatives are bounded. The parameter estimate and estimation error are

\[ \hat{\boldsymbol{\theta}}(t)\in\mathbb{R}^p, \qquad \tilde{\boldsymbol{\theta}} =\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}. \]

Define the first tracking error and subsequent backstepping errors by

\[ z_1=x_1-y_d,\qquad z_i=x_i-\alpha_{i-1},\quad i=2,\ldots,n, \]

where \(\alpha_i\) is the virtual control selected at step \(i\). The recursive target is to preserve the derivative form

\[ \dot V_i = -\sum_{j=1}^{i}c_j z_j^2 +z_i z_{i+1} +\tilde{\boldsymbol{\theta}}^{\mathsf T} \left(\boldsymbol{\tau}_i- \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}}\right), \]

for \(i<n\), with \(c_j>0\) and \(\boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\mathsf T}>0\). The term \(z_i z_{i+1}\) is intentionally left for the next recursive step. The tuning function \(\boldsymbol{\tau}_i\) accumulates the effective regressors generated by the chain rule.

flowchart TD
  A["Start: z1 and V1"] --> B["Differentiate Vi"]
  B --> C["Identify known terms, zi*z(i+1), and parameter-error terms"]
  C --> D["Choose virtual control alpha_i"]
  D --> E["Augment: V(i+1) = Vi + 0.5*z(i+1)^2"]
  E --> F["Update tuning function tau_i"]
  F --> G{"Final state reached?"}
  G -->|"no"| B
  G -->|"yes"| H["Choose actual u and theta-hat-dot"]
  H --> I["Obtain V-dot = negative quadratic form"]
        

A state-only Lyapunov function cannot directly absorb products involving the unknown parameter error. The standard composite candidate is

\[ V_n(\mathbf{z},\tilde{\boldsymbol{\theta}}) =\frac{1}{2}\sum_{i=1}^{n}z_i^2 +\frac{1}{2}\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}. \]

Because the true parameter is constant,

\[ \dot{\tilde{\boldsymbol{\theta}}} =-\dot{\hat{\boldsymbol{\theta}}}, \]

and therefore the parameter-energy derivative is

\[ \frac{d}{dt}\left( \frac{1}{2}\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}} \right) =-\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}}. \]

This expression supplies the exact algebraic channel needed to cancel terms such as \(\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\tau}_n\). The matrix \(\boldsymbol{\Gamma}\) is an adaptation-gain matrix; it changes adaptation speed but not the sign of the ideal Lyapunov derivative.

The parameter-error energy is included once in the composite function. Repeating it at every recursive step would incorrectly count the same unknown parameter multiple times and would create artificial update laws.

For the first subsystem, define

\[ V_{z,1}=\frac{1}{2}z_1^2. \]

Its derivative has the generic structure

\[ \dot V_{z,1}=z_1\left[x_2+f_1(x_1) +\boldsymbol{\phi}_1^{\mathsf T}(x_1)\boldsymbol{\theta} -\dot y_d\right]. \]

Add and subtract an estimate-based virtual control. A typical choice is

\[ \alpha_1=-c_1z_1-f_1(x_1) -\boldsymbol{\phi}_1^{\mathsf T}(x_1) \hat{\boldsymbol{\theta}}+\dot y_d. \]

Since \(x_2=z_2+\alpha_1\), substitution gives

\[ \dot V_{z,1}=-c_1z_1^2+z_1z_2 +z_1\boldsymbol{\phi}_1^{\mathsf T}(x_1) \tilde{\boldsymbol{\theta}}. \]

Three components are visible: desired dissipation, a recursive coupling to the next error, and a parameter-error cross term. The later steps are built to preserve this decomposition rather than trying to eliminate everything at once.

At step \(i\), augment the state portion as

\[ V_{z,i}=V_{z,i-1}+\frac{1}{2}z_i^2. \]

Differentiation gives

\[ \dot V_{z,i}=\dot V_{z,i-1}+z_i\dot z_i, \qquad \dot z_i=\dot x_i-\dot\alpha_{i-1}. \]

The chain rule in \(\dot\alpha_{i-1}\) introduces known state terms, reference derivatives, parameter-estimate derivatives, and an effective parameter regressor. Denote that effective regressor by \(\boldsymbol{\omega}_i\). Before the tuning-function corrections, its principal parametric part is

\[ \boldsymbol{\omega}_i =\boldsymbol{\phi}_i -\sum_{k=1}^{i-1} \frac{\partial\alpha_{i-1}}{\partial x_k} \boldsymbol{\phi}_k. \]

The cumulative tuning function is organized recursively as

\[ \boldsymbol{\tau}_0=\mathbf{0},\qquad \boldsymbol{\tau}_i=\boldsymbol{\tau}_{i-1} +z_i\boldsymbol{\omega}_i. \]

The virtual control \(\alpha_i\) is then selected to cancel known drift terms, cancel the predecessor coupling \(z_{i-1}z_i\), add \(-c_i z_i\), and compensate the terms produced by \(\partial\alpha_{i-1}/\partial\hat{\boldsymbol{\theta}}\). After those algebraic choices, the design invariant becomes

\[ \dot V_i=-\sum_{j=1}^{i}c_jz_j^2+z_i z_{i+1} +\tilde{\boldsymbol{\theta}}^{\mathsf T} \left(\boldsymbol{\tau}_i- \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}}\right). \]

At the final step, the actual input removes the last state coupling and the adaptive law is chosen as

\[ \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}\boldsymbol{\tau}_n, \]

which yields the ideal derivative

\[ \dot V_n=-\sum_{j=1}^{n}c_jz_j^2\le 0. \]

Consider the nonlinear strict-feedback plant

\[ \dot x_1=x_2,\qquad \dot x_2=u+\theta\sin(x_1), \]

where \(\theta\) is an unknown constant. Let \(r(t)\), \(\dot r(t)\), and \(\ddot r(t)\) be bounded.

7.1 Step 1: Virtual Control

\[ z_1=x_1-r,\qquad V_1=\frac{1}{2}z_1^2. \]

Since \(\dot z_1=x_2-\dot r\), select

\[ \alpha_1=\dot r-c_1z_1,\qquad z_2=x_2-\alpha_1. \]

Then

\[ \dot z_1=-c_1z_1+z_2,\qquad \dot V_1=-c_1z_1^2+z_1z_2. \]

7.2 Step 2: Actual Control and Parameter Adaptation

The virtual-control derivative is implementable:

\[ \dot\alpha_1=\ddot r-c_1(x_2-\dot r). \]

The second error satisfies

\[ \dot z_2=u+\theta\sin(x_1)-\dot\alpha_1. \]

Choose the actual control

\[ u=\dot\alpha_1-z_1-c_2z_2 -\hat\theta\sin(x_1). \]

Using \(\tilde\theta=\theta-\hat\theta\),

\[ \dot z_2=-z_1-c_2z_2+\tilde\theta\sin(x_1). \]

Now introduce the composite candidate

\[ V=\frac{1}{2}z_1^2+\frac{1}{2}z_2^2 +\frac{1}{2\gamma}\tilde\theta^2, \qquad \gamma>0. \]

Its derivative is

\[ \begin{aligned} \dot V &=z_1(-c_1z_1+z_2) +z_2[-z_1-c_2z_2+\tilde\theta\sin(x_1)] -\frac{1}{\gamma}\tilde\theta\dot{\hat\theta}\\ &=-c_1z_1^2-c_2z_2^2 +\tilde\theta\left[z_2\sin(x_1) -\frac{1}{\gamma}\dot{\hat\theta}\right]. \end{aligned} \]

Select the adaptive law

\[ \dot{\hat\theta}=\gamma z_2\sin(x_1). \]

The final derivative is therefore

\[ \boxed{\dot V=-c_1z_1^2-c_2z_2^2\le 0}. \]

Theorem. Suppose:

1. \(c_1>0\), \(c_2>0\), and \(\gamma>0\);
2. \(r\), \(\dot r\), and \(\ddot r\) are bounded; and
3. the plant and controller equations of Section 7 hold without actuator saturation.

Then all closed-loop signals are bounded and \(z_1(t)\to0\), \(z_2(t)\to0\) as \(t\to\infty\).

Proof. The derived identity gives

\[ V(t)+\int_0^t\left(c_1z_1^2(\sigma)+c_2z_2^2(\sigma)\right)d\sigma =V(0). \]

Consequently,

\[ V(t)\le V(0),\qquad z_1,z_2,\tilde\theta\in\mathcal{L}_{\infty},\qquad z_1,z_2\in\mathcal{L}_2. \]

Because \(r\) is bounded and \(x_1=z_1+r\), \(x_1\) is bounded. The virtual control \(\alpha_1=\dot r-c_1z_1\) is bounded, so \(x_2=z_2+\alpha_1\) is bounded. The estimate \(\hat\theta=\theta-\tilde\theta\) is bounded. Hence \(\sin(x_1)\), \(\dot{\hat\theta}\), \(u\), \(\dot z_1\), and \(\dot z_2\) are bounded.

Therefore \(z_1\) and \(z_2\) are uniformly continuous. Since they also belong to \(\mathcal{L}_2\), Barbalat's lemma gives

\[ \lim_{t\to\infty}z_1(t)=0,\qquad \lim_{t\to\infty}z_2(t)=0. \]

This proves asymptotic tracking and boundedness. \(\square\)

The derivative is negative definite in the tracking errors but only negative semidefinite in the complete state–parameter vector. It does not contain a term such as \(-k_\theta\tilde\theta^2\). Therefore the proof guarantees bounded parameter error, not parameter convergence.

\[ z_1(t)\to0,\quad z_2(t)\to0 \quad\not\Rightarrow\quad \tilde\theta(t)\to0. \]

Exact parameter convergence requires additional excitation or identifiability conditions. A controller can track perfectly while retaining an incorrect parameter estimate if the closed-loop trajectory does not reveal enough information about the unknown coefficient.

This distinction is central in adaptive control: the Lyapunov function is designed first to establish stability and tracking. Identification is a stronger objective and must be analyzed separately.

The implementable signal sequence for the worked example is:

flowchart TD
  R["r, r_dot, r_ddot"] --> E1["z1 = x1 - r"]
  X1["x1"] --> E1
  E1 --> A1["alpha1 = r_dot - c1*z1"]
  X2["x2"] --> E2["z2 = x2 - alpha1"]
  A1 --> E2
  E2 --> AD["theta_hat_dot = gamma*z2*sin(x1)"]
  AD --> INT["Parameter estimator integrator"]
  INT --> TH["theta_hat"]
  E1 --> U["u = alpha1_dot - z1 - c2*z2 - theta_hat*sin(x1)"]
  E2 --> U
  TH --> U
  U --> P["Strict-feedback plant"]
  P --> X1
  P --> X2
        

Recommended numerical verification steps are:

1. compute \(z_1\), \(z_2\), and \(V\) at every integration step;
2. confirm that the numerical trajectory of \(V\) is nonincreasing up to integration tolerance;
3. compare a numerical derivative of \(V\) with \(-c_1z_1^2-c_2z_2^2\);
4. repeat with several initial estimates and adaptation gains; and
5. avoid interpreting parameter convergence in one simulation as a general theorem.

For the supplied test case \(\theta=1.5\), \(c_1=2\), \(c_2=2.5\), \(\gamma=5\), and \(r(t)=0.8\sin(0.6t)\), the tested Python, C++, and Java implementations produce final tracking errors on the order of \(10^{-8}\).

This implementation uses scipy.integrate.solve_ivp for high-accuracy integration, NumPy for vector operations, and Matplotlib for state, error, parameter, and Lyapunov plots.

Chapter17_Lesson4.py

"""Chapter17_Lesson4.py
Adaptive backstepping simulation for a second-order strict-feedback system.

Plant:
    x1_dot = x2
    x2_dot = u + theta * sin(x1)

Controller:
    z1 = x1 - r
    alpha1 = r_dot - c1*z1
    z2 = x2 - alpha1
    u = alpha1_dot - z1 - c2*z2 - theta_hat*sin(x1)
    theta_hat_dot = gamma*z2*sin(x1)
"""

from __future__ import annotations

import math
from dataclasses import dataclass

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


@dataclass(frozen=True)
class Config:
    theta: float = 1.5
    c1: float = 2.0
    c2: float = 2.5
    gamma: float = 5.0
    t_final: float = 25.0


def reference(t: float) -> tuple[float, float, float]:
    """Return r(t), r_dot(t), and r_ddot(t)."""
    r = 0.8 * math.sin(0.6 * t)
    r_dot = 0.48 * math.cos(0.6 * t)
    r_ddot = -0.288 * math.sin(0.6 * t)
    return r, r_dot, r_ddot


def closed_loop(t: float, state: np.ndarray, cfg: Config) -> np.ndarray:
    x1, x2, theta_hat = state
    r, r_dot, r_ddot = reference(t)

    z1 = x1 - r
    alpha1 = r_dot - cfg.c1 * z1
    z2 = x2 - alpha1
    phi = math.sin(x1)

    alpha1_dot = r_ddot - cfg.c1 * (x2 - r_dot)
    u = alpha1_dot - z1 - cfg.c2 * z2 - theta_hat * phi
    theta_hat_dot = cfg.gamma * z2 * phi

    x1_dot = x2
    x2_dot = u + cfg.theta * phi
    return np.array([x1_dot, x2_dot, theta_hat_dot], dtype=float)


def main() -> None:
    cfg = Config()
    initial_state = np.array([1.2, -0.4, 0.0], dtype=float)
    sample_times = np.linspace(0.0, cfg.t_final, 2501)

    solution = solve_ivp(
        fun=lambda t, y: closed_loop(t, y, cfg),
        t_span=(0.0, cfg.t_final),
        y0=initial_state,
        t_eval=sample_times,
        rtol=1.0e-9,
        atol=1.0e-11,
    )
    if not solution.success:
        raise RuntimeError(solution.message)

    t = solution.t
    x1, x2, theta_hat = solution.y
    r = 0.8 * np.sin(0.6 * t)
    r_dot = 0.48 * np.cos(0.6 * t)
    z1 = x1 - r
    alpha1 = r_dot - cfg.c1 * z1
    z2 = x2 - alpha1
    theta_tilde = cfg.theta - theta_hat
    lyapunov = 0.5 * z1**2 + 0.5 * z2**2 + 0.5 * theta_tilde**2 / cfg.gamma

    print(f"maximum |z1|: {np.max(np.abs(z1)):.6f}")
    print(f"final |z1|:   {abs(z1[-1]):.6e}")
    print(f"final |z2|:   {abs(z2[-1]):.6e}")
    print(f"final theta_hat: {theta_hat[-1]:.6f}")
    print("Tracking convergence does not, by itself, guarantee theta_hat -> theta.")

    plt.figure()
    plt.plot(t, x1, label="x1")
    plt.plot(t, r, "--", label="reference")
    plt.xlabel("Time [s]")
    plt.ylabel("Position")
    plt.title("Adaptive backstepping tracking")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(t, z1, label="z1")
    plt.plot(t, z2, label="z2")
    plt.xlabel("Time [s]")
    plt.ylabel("Backstepping errors")
    plt.title("Recursive error coordinates")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(t, theta_hat, label="theta_hat")
    plt.axhline(cfg.theta, linestyle="--", label="theta")
    plt.xlabel("Time [s]")
    plt.ylabel("Parameter")
    plt.title("Adaptive parameter estimate")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(t, lyapunov)
    plt.xlabel("Time [s]")
    plt.ylabel("V")
    plt.title("Composite Lyapunov function")
    plt.grid(True)

    plt.show()


if __name__ == "__main__":
    main()

The C++ version uses only the standard library and a fixed-step fourth-order Runge–Kutta integrator. It writes all signals to a CSV file for plotting or independent verification.

Chapter17_Lesson4.cpp

// Chapter17_Lesson4.cpp
// Standalone RK4 simulation of adaptive backstepping.

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

struct Config {
    double theta = 1.5;
    double c1 = 2.0;
    double c2 = 2.5;
    double gamma = 5.0;
    double dt = 0.002;
    double tFinal = 25.0;
};

using State = std::array<double, 3>; // x1, x2, thetaHat

struct Reference {
    double r;
    double rDot;
    double rDDot;
};

Reference reference(double t) {
    return {
        0.8 * std::sin(0.6 * t),
        0.48 * std::cos(0.6 * t),
        -0.288 * std::sin(0.6 * t)
    };
}

State dynamics(double t, const State& y, const Config& cfg) {
    const double x1 = y[0];
    const double x2 = y[1];
    const double thetaHat = y[2];
    const Reference ref = reference(t);

    const double z1 = x1 - ref.r;
    const double alpha1 = ref.rDot - cfg.c1 * z1;
    const double z2 = x2 - alpha1;
    const double phi = std::sin(x1);
    const double alpha1Dot = ref.rDDot - cfg.c1 * (x2 - ref.rDot);
    const double u = alpha1Dot - z1 - cfg.c2 * z2 - thetaHat * phi;

    return {
        x2,
        u + cfg.theta * phi,
        cfg.gamma * z2 * phi
    };
}

State addScaled(const State& a, const State& b, double scale) {
    return {a[0] + scale * b[0], a[1] + scale * b[1], a[2] + scale * b[2]};
}

State rk4Step(double t, const State& y, const Config& cfg) {
    const State k1 = dynamics(t, y, cfg);
    const State k2 = dynamics(t + 0.5 * cfg.dt, addScaled(y, k1, 0.5 * cfg.dt), cfg);
    const State k3 = dynamics(t + 0.5 * cfg.dt, addScaled(y, k2, 0.5 * cfg.dt), cfg);
    const State k4 = dynamics(t + cfg.dt, addScaled(y, k3, cfg.dt), cfg);

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

int main() {
    const Config cfg;
    State y{1.2, -0.4, 0.0};
    std::ofstream csv("Chapter17_Lesson4_cpp.csv");
    if (!csv) {
        throw std::runtime_error("Cannot open output CSV file.");
    }
    csv << "t,x1,x2,r,z1,z2,theta_hat,V\n";
    csv << std::setprecision(12);

    const int steps = static_cast<int>(std::ceil(cfg.tFinal / cfg.dt));
    for (int k = 0; k <= steps; ++k) {
        const double t = k * cfg.dt;
        const Reference ref = reference(t);
        const double z1 = y[0] - ref.r;
        const double alpha1 = ref.rDot - cfg.c1 * z1;
        const double z2 = y[1] - alpha1;
        const double thetaTilde = cfg.theta - y[2];
        const double V = 0.5 * z1 * z1 + 0.5 * z2 * z2
                       + 0.5 * thetaTilde * thetaTilde / cfg.gamma;
        csv << t << ',' << y[0] << ',' << y[1] << ',' << ref.r << ','
            << z1 << ',' << z2 << ',' << y[2] << ',' << V << '\n';

        if (k < steps) {
            y = rk4Step(t, y, cfg);
        }
    }

    const Reference finalRef = reference(cfg.tFinal);
    const double finalZ1 = y[0] - finalRef.r;
    const double finalAlpha1 = finalRef.rDot - cfg.c1 * finalZ1;
    const double finalZ2 = y[1] - finalAlpha1;
    std::cout << "Final |z1| = " << std::abs(finalZ1) << '\n';
    std::cout << "Final |z2| = " << std::abs(finalZ2) << '\n';
    std::cout << "Final theta_hat = " << y[2] << '\n';
    std::cout << "Data written to Chapter17_Lesson4_cpp.csv\n";
    return 0;
}

The Java implementation also uses a standalone RK4 solver and writes a CSV record. No external control or numerical library is required.

Chapter17_Lesson4.java

// Chapter17_Lesson4.java
// Standalone RK4 simulation of adaptive backstepping.

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

public final class Chapter17_Lesson4 {
    private static final double THETA = 1.5;
    private static final double C1 = 2.0;
    private static final double C2 = 2.5;
    private static final double GAMMA = 5.0;
    private static final double DT = 0.002;
    private static final double T_FINAL = 25.0;

    private Chapter17_Lesson4() {
    }

    private static double[] reference(double t) {
        return new double[] {
            0.8 * Math.sin(0.6 * t),
            0.48 * Math.cos(0.6 * t),
            -0.288 * Math.sin(0.6 * t)
        };
    }

    private static double[] dynamics(double t, double[] y) {
        double x1 = y[0];
        double x2 = y[1];
        double thetaHat = y[2];
        double[] ref = reference(t);

        double z1 = x1 - ref[0];
        double alpha1 = ref[1] - C1 * z1;
        double z2 = x2 - alpha1;
        double phi = Math.sin(x1);
        double alpha1Dot = ref[2] - C1 * (x2 - ref[1]);
        double u = alpha1Dot - z1 - C2 * z2 - thetaHat * phi;

        return new double[] {
            x2,
            u + THETA * phi,
            GAMMA * z2 * phi
        };
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        return new double[] {
            a[0] + scale * b[0],
            a[1] + scale * b[1],
            a[2] + scale * b[2]
        };
    }

    private static double[] rk4Step(double t, double[] y) {
        double[] k1 = dynamics(t, y);
        double[] k2 = dynamics(t + 0.5 * DT, addScaled(y, k1, 0.5 * DT));
        double[] k3 = dynamics(t + 0.5 * DT, addScaled(y, k2, 0.5 * DT));
        double[] k4 = dynamics(t + DT, addScaled(y, k3, DT));
        return new double[] {
            y[0] + DT * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            y[1] + DT * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
            y[2] + DT * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0
        };
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        double[] y = {1.2, -0.4, 0.0};
        int steps = (int) Math.ceil(T_FINAL / DT);

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter17_Lesson4_java.csv"))) {
            writer.write("t,x1,x2,r,z1,z2,theta_hat,V\n");
            for (int k = 0; k <= steps; ++k) {
                double t = k * DT;
                double[] ref = reference(t);
                double z1 = y[0] - ref[0];
                double alpha1 = ref[1] - C1 * z1;
                double z2 = y[1] - alpha1;
                double thetaTilde = THETA - y[2];
                double v = 0.5 * z1 * z1 + 0.5 * z2 * z2
                         + 0.5 * thetaTilde * thetaTilde / GAMMA;
                writer.write(String.format(Locale.US,
                        "%.9f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        t, y[0], y[1], ref[0], z1, z2, y[2], v));

                if (k < steps) {
                    y = rk4Step(t, y);
                }
            }
        }

        double[] finalRef = reference(T_FINAL);
        double finalZ1 = y[0] - finalRef[0];
        double finalAlpha1 = finalRef[1] - C1 * finalZ1;
        double finalZ2 = y[1] - finalAlpha1;
        System.out.printf("Final |z1| = %.8e%n", Math.abs(finalZ1));
        System.out.printf("Final |z2| = %.8e%n", Math.abs(finalZ2));
        System.out.printf("Final theta_hat = %.8f%n", y[2]);
        System.out.println("Data written to Chapter17_Lesson4_java.csv");
    }
}

The MATLAB script uses ode45. For Simulink, use three Integrator blocks for \(x_1\), \(x_2\), and \(\hat\theta\), and place the equations in controllerSignals inside a MATLAB Function block. The script comments specify the required signal connections.

Chapter17_Lesson4.m

% Chapter17_Lesson4.m
% Adaptive backstepping for a second-order strict-feedback system.
% Requires base MATLAB. The same equations can be placed in a Simulink
% MATLAB Function block with three Integrator blocks for x1, x2, theta_hat.

clear; clc; close all;

theta = 1.5;
c1 = 2.0;
c2 = 2.5;
gamma = 5.0;
tspan = [0 25];
y0 = [1.2; -0.4; 0.0];

rhs = @(t,y) closedLoop(t, y, theta, c1, c2, gamma);
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, y] = ode45(rhs, tspan, y0, options);

x1 = y(:,1);
x2 = y(:,2);
thetaHat = y(:,3);
r = 0.8*sin(0.6*t);
rDot = 0.48*cos(0.6*t);
z1 = x1 - r;
alpha1 = rDot - c1*z1;
z2 = x2 - alpha1;
thetaTilde = theta - thetaHat;
V = 0.5*z1.^2 + 0.5*z2.^2 + 0.5*(thetaTilde.^2)/gamma;

fprintf('Final |z1| = %.8e\n', abs(z1(end)));
fprintf('Final |z2| = %.8e\n', abs(z2(end)));
fprintf('Final theta_hat = %.8f\n', thetaHat(end));

figure;
plot(t, x1, 'LineWidth', 1.4); hold on;
plot(t, r, '--', 'LineWidth', 1.4);
grid on; xlabel('Time [s]'); ylabel('Position');
title('Adaptive backstepping tracking');
legend('x_1', 'reference', 'Location', 'best');

figure;
plot(t, z1, 'LineWidth', 1.4); hold on;
plot(t, z2, 'LineWidth', 1.4);
grid on; xlabel('Time [s]'); ylabel('Error');
title('Recursive error coordinates');
legend('z_1', 'z_2', 'Location', 'best');

figure;
plot(t, thetaHat, 'LineWidth', 1.4); hold on;
yline(theta, '--', 'LineWidth', 1.4);
grid on; xlabel('Time [s]'); ylabel('Parameter');
title('Adaptive parameter estimate');
legend('\theta-hat', '\theta', 'Location', 'best');

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

% Simulink realization:
% 1. Use three Integrator blocks for x1, x2, and thetaHat.
% 2. Put the equations in controllerSignals below into a MATLAB Function block.
% 3. Feed x2 to x1_dot, u + theta*sin(x1) to x2_dot, and thetaHatDot
%    to the parameter-estimate integrator.

function dy = closedLoop(t, y, theta, c1, c2, gamma)
    x1 = y(1);
    x2 = y(2);
    thetaHat = y(3);
    [u, thetaHatDot] = controllerSignals(t, x1, x2, thetaHat, c1, c2, gamma);
    phi = sin(x1);
    dy = [x2; u + theta*phi; thetaHatDot];
end

function [u, thetaHatDot] = controllerSignals(t, x1, x2, thetaHat, c1, c2, gamma)
    r = 0.8*sin(0.6*t);
    rDot = 0.48*cos(0.6*t);
    rDDot = -0.288*sin(0.6*t);
    z1 = x1 - r;
    alpha1 = rDot - c1*z1;
    z2 = x2 - alpha1;
    phi = sin(x1);
    alpha1Dot = rDDot - c1*(x2 - rDot);
    u = alpha1Dot - z1 - c2*z2 - thetaHat*phi;
    thetaHatDot = gamma*z2*phi;
end

The Wolfram Language source uses NDSolveValue and directly defines the recursive errors, controller, adaptive law, and composite Lyapunov function.

Chapter17_Lesson4.nb

(* Chapter17_Lesson4.nb
   Wolfram Language source that can be pasted into a Mathematica notebook.
   Adaptive backstepping for a second-order strict-feedback system. *)

ClearAll["Global`*"];

theta = 1.5;
c1 = 2.0;
c2 = 2.5;
gamma = 5.0;
tFinal = 25.0;

r[t_] := 0.8 Sin[0.6 t];
rDot[t_] := 0.48 Cos[0.6 t];
rDDot[t_] := -0.288 Sin[0.6 t];
phi[x_] := Sin[x];

z1[t_] := x1[t] - r[t];
alpha1[t_] := rDot[t] - c1 z1[t];
z2[t_] := x2[t] - alpha1[t];
alpha1Dot[t_] := rDDot[t] - c1 (x2[t] - rDot[t]);
u[t_] := alpha1Dot[t] - z1[t] - c2 z2[t] - thetaHat[t] phi[x1[t]];

eqns = {
   x1'[t] == x2[t],
   x2'[t] == u[t] + theta phi[x1[t]],
   thetaHat'[t] == gamma z2[t] phi[x1[t]],
   x1[0] == 1.2,
   x2[0] == -0.4,
   thetaHat[0] == 0.0
};

sol = NDSolveValue[eqns, {x1, x2, thetaHat}, {t, 0, tFinal},
   Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5}},
   AccuracyGoal -> 10, PrecisionGoal -> 10];

x1s[t_] := sol[[1]][t];
x2s[t_] := sol[[2]][t];
thetaHats[t_] := sol[[3]][t];
z1s[t_] := x1s[t] - r[t];
alpha1s[t_] := rDot[t] - c1 z1s[t];
z2s[t_] := x2s[t] - alpha1s[t];
Vs[t_] := 1/2 z1s[t]^2 + 1/2 z2s[t]^2 + (theta - thetaHats[t])^2/(2 gamma);

Print["Final |z1| = ", N[Abs[z1s[tFinal]], 8]];
Print["Final |z2| = ", N[Abs[z2s[tFinal]], 8]];
Print["Final thetaHat = ", N[thetaHats[tFinal], 8]];

trackingPlot = Plot[{x1s[t], r[t]}, {t, 0, tFinal},
   PlotLegends -> {"x1", "reference"}, AxesLabel -> {"t", "position"},
   PlotLabel -> "Adaptive backstepping tracking", PlotRange -> All];

errorPlot = Plot[{z1s[t], z2s[t]}, {t, 0, tFinal},
   PlotLegends -> {"z1", "z2"}, AxesLabel -> {"t", "error"},
   PlotLabel -> "Recursive error coordinates", PlotRange -> All];

parameterPlot = Plot[{thetaHats[t], theta}, {t, 0, tFinal},
   PlotLegends -> {"thetaHat", "theta"}, AxesLabel -> {"t", "parameter"},
   PlotLabel -> "Adaptive parameter estimate", PlotRange -> All];

lyapunovPlot = Plot[Vs[t], {t, 0, tFinal}, AxesLabel -> {"t", "V"},
   PlotLabel -> "Composite Lyapunov function", PlotRange -> All];

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

1. Inconsistent parameter-error sign: choosing \(\tilde\theta=\hat\theta-\theta\) but using the update law derived for \(\tilde\theta=\theta-\hat\theta\) reverses the cancellation term.
2. Omitting the predecessor cancellation: the actual control must contain \(-z_1\) in the second-order example so that \(z_1z_2\) cancels.
3. Differentiating an unavailable signal: every derivative in \(\dot\alpha_i\) must be generated from known models, measured states, filters, or command-filter methods.
4. Adding saturation after the proof: clipping the control changes the error dynamics and invalidates the exact Lyapunov identity.
5. Claiming exponential convergence of all variables: \(\dot V\) is not negative definite in \(\tilde\theta\) under the basic law.
6. Using excessive adaptation gain: the ideal proof permits any positive gain, but measurement noise and discretization can make a large gain impractical.

Problem 1 (Scalar Cancellation Law). Suppose

\[ \dot V=-a z^2+\tilde\theta z\phi -\frac{1}{\gamma}\tilde\theta\dot{\hat\theta}, \qquad a>0. \]

Find an adaptive law that makes \(\dot V=-az^2\).

Solution. Factor the parameter-error terms:

\[ \dot V=-az^2+\tilde\theta \left(z\phi-\frac{1}{\gamma}\dot{\hat\theta}\right). \]

Set

\[ \dot{\hat\theta}=\gamma z\phi. \]

The parenthesized term vanishes identically.

Problem 2 (Reversed Error Convention). Redefine the parameter error as \(\bar\theta=\hat\theta-\theta\). Determine the corresponding update law for the worked example.

Solution. Now \(\dot{\bar\theta}=\dot{\hat\theta}\), and the second-error dynamics become

\[ \dot z_2=-z_1-c_2z_2-\bar\theta\sin(x_1). \]

The Lyapunov derivative contains

\[ \bar\theta\left[-z_2\sin(x_1) +\frac{1}{\gamma}\dot{\hat\theta}\right]. \]

Thus the update law remains numerically

\[ \dot{\hat\theta}=\gamma z_2\sin(x_1), \]

but the intermediate signs differ. This illustrates why the convention must be stated before differentiating the Lyapunov function.

Problem 3 (Vector Parameter). Replace the uncertain term by \(\boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}(x_1)\), where \(\boldsymbol{\theta}\in\mathbb{R}^p\). Derive the update law.

Solution. Choose

\[ V=\frac{1}{2}z_1^2+\frac{1}{2}z_2^2 +\frac{1}{2}\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}, \]

which produces the parameter term

\[ \tilde{\boldsymbol{\theta}}^{\mathsf T} \left[z_2\boldsymbol{\phi}(x_1) -\boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}\right]. \]

Therefore

\[ \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}z_2\boldsymbol{\phi}(x_1). \]

Problem 4 (Integral Bound). Show that the worked controller implies a finite weighted tracking-energy integral.

Solution. Integrate the Lyapunov identity:

\[ \int_0^T\left(c_1z_1^2+c_2z_2^2\right)dt =V(0)-V(T)\le V(0). \]

Letting \(T\to\infty\) gives

\[ \int_0^\infty\left(c_1z_1^2+c_2z_2^2\right)dt\le V(0). \]

Problem 5 (Why Tracking Does Not Identify the Parameter). Consider a constant reference \(r=0\) and initial condition \(x_1(0)=x_2(0)=0\). What happens for an arbitrary initial estimate?

Solution. At the origin,

\[ z_1=z_2=0,\qquad \sin(x_1)=0, \qquad \dot{\hat\theta}=0. \]

Tracking is exact, but the estimate remains at its initial value. The data contain no information about \(\theta\), so parameter convergence is impossible from this trajectory.

Problem 6 (Effect of Control Saturation). Suppose the implemented input is \(u_s=\operatorname{sat}(u)\). Write the modified Lyapunov derivative.

Solution. Define \(\Delta_u=u_s-u\). The second error dynamics gain the term \(\Delta_u\), so

\[ \dot V=-c_1z_1^2-c_2z_2^2+z_2\Delta_u. \]

The final term has no guaranteed sign. Therefore the unsaturated proof no longer establishes monotonic decrease of \(V\); saturation compensation or a constrained analysis is required.

Adaptive backstepping constructs stability recursively. Each new squared error term absorbs one additional strict-feedback state, each virtual control creates dissipation and cancels the previous coupling, and one parameter-error energy term creates the channel for an adaptive cancellation law. For the worked nonlinear plant, the resulting identity \(\dot V=-c_1z_1^2-c_2z_2^2\) establishes boundedness and asymptotic tracking through Lyapunov arguments and Barbalat's lemma. It does not, without excitation, establish convergence of the parameter estimate.