In Lesson 1, strict-feedback systems with parametric uncertainty were introduced. In Lesson 2, virtual controls were allowed to depend on the current parameter estimate. The unresolved difficulty is that differentiating such a virtual control produces \( \dot{\hat{\boldsymbol{\theta}}} \) inside the next backstepping equation. Therefore, the adaptive law cannot be designed as an unrelated afterthought; it must be constructed recursively with the virtual controls.

After completing this lesson, students should be able to:

• derive the propagated parameter regressor at each backstepping stage;
• explain why independent stagewise parameter estimates cause overparameterization;
• construct recursive tuning functions \( \boldsymbol{\tau}_i \) and the terminal update law;
• prove cancellation of state–parameter cross terms for a second-order strict-feedback plant;
• implement the resulting closed loop with fixed-step Runge–Kutta integration.

Consider the parametric strict-feedback chain

\[ \begin{aligned} \dot{x}_1 &= x_2 + \boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_1(x_1),\\ \dot{x}_2 &= x_3 + \boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_2(x_1,x_2),\\ &\ \vdots\\ \dot{x}_{n-1} &= x_n + \boldsymbol{\theta}^{\mathsf T} \boldsymbol{\phi}_{n-1}(\mathbf{x}_{1:n-1}),\\ \dot{x}_n &= u + \boldsymbol{\theta}^{\mathsf T} \boldsymbol{\phi}_n(\mathbf{x}). \end{aligned} \]

The unknown vector \( \boldsymbol{\theta}\in\mathbb{R}^{p} \) is constant, all regressors \( \boldsymbol{\phi}_i \) are known and sufficiently smooth, and the full state is measured. The reference \( y_d(t) \) and the derivatives needed by the design are bounded and available. The parameter estimate and error are

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

A positive-definite adaptation-gain matrix \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\mathsf T}>0 \) is selected. The key structural requirement is that the same physical parameter vector appears throughout the chain. Consequently, only one estimate vector should be updated.

Ordinary backstepping differentiates each virtual control when proceeding to the next state equation. In adaptive backstepping, a virtual control generally contains \( \hat{\boldsymbol{\theta}} \). Hence its total derivative contains both plant uncertainty and the yet unknown update signal:

\[ \dot{\alpha}_i =\frac{\partial\alpha_i}{\partial t} +\sum_{j=1}^{i} \frac{\partial\alpha_i}{\partial x_j}\dot{x}_j +\frac{\partial\alpha_i}{\partial\hat{\boldsymbol{\theta}}} \dot{\hat{\boldsymbol{\theta}}}. \]

Three objects must therefore be propagated from one stage to the next: the state-error coordinate, the effective parameter regressor, and the dependence on the adaptive update. A tuning function stores the cumulative parameter-error cancellation requested by all completed stages.

flowchart TD
  A["Stage i: define error zi"] --> B["Differentiate zi"]
  B --> C["Collect known drift terms"]
  B --> D["Collect shared parameter \nregressor omega_i"]
  B --> E["Collect dependence \non theta_hat_dot"]
  C --> F["Choose virtual control alpha_i"]
  D --> G["Update tuning function tau_i"]
  E --> F
  G --> H["Proceed to stage i+1"]
  F --> H
  H --> I["At final stage set theta_hat_dot = tau_n"]
        

A naïve alternative would introduce a separate estimate \( \hat{\boldsymbol{\theta}}_i \) at every stage. For a plant with \( p \) unknown constants and \( n \) states, that choice can create \( np \) adaptive states even though only \( p \) physical constants are unknown. Tuning functions avoid this duplication.

The central derivation uses the two-state plant

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

The objective is asymptotic tracking of a smooth bounded reference \( y_d(t) \). Choose gains \( c_1>0 \) and \( c_2>0 \). Define the first error coordinate

\[ z_1=x_1-y_d. \]

Treating \( x_2 \) as the first virtual input, select

\[ \alpha_1=-c_1z_1- \hat{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1(x_1) +\dot{y}_d, \qquad z_2=x_2-\alpha_1. \]

Substitution gives the first transformed error equation:

\[ \dot{z}_1=-c_1z_1+z_2+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1. \]

Introduce the partial Lyapunov function

\[ V_1=\frac{1}{2}z_1^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}}} \). Therefore,

\[ \dot{V}_1=-c_1z_1^2+z_1z_2+ \tilde{\boldsymbol{\theta}}^{\mathsf T} \left( \boldsymbol{\phi}_1z_1- \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}} \right). \]

At this point it would be premature to set the actual adaptive law equal to \( \boldsymbol{\Gamma}\boldsymbol{\phi}_1z_1 \), because the second error equation will generate another state–parameter product. Instead, define the first tuning function

\[ \boldsymbol{\tau}_1 =\boldsymbol{\Gamma}\boldsymbol{\phi}_1z_1. \]

This is a bookkeeping signal, not yet the terminal update law. It records the adaptation contribution requested by stage 1.

Let

\[ a_1(x_1,\hat{\boldsymbol{\theta}}) =\frac{\partial\alpha_1}{\partial x_1} =-c_1- \hat{\boldsymbol{\theta}}^{\mathsf T} \frac{\partial\boldsymbol{\phi}_1}{\partial x_1}. \]

Since \( \partial\alpha_1/ \partial\hat{\boldsymbol{\theta}} =-\boldsymbol{\phi}_1^{\mathsf T} \), the total derivative is

\[ \dot{\alpha}_1 =a_1\left(x_2+ \boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_1\right) +c_1\dot{y}_d+\ddot{y}_d -\boldsymbol{\phi}_1^{\mathsf T} \dot{\hat{\boldsymbol{\theta}}}. \]

Differentiating \( z_2=x_2-\alpha_1 \) gives

\[ \begin{aligned} \dot{z}_2={}&u-a_1x_2-c_1\dot{y}_d-\ddot{y}_d\\ &+\boldsymbol{\theta}^{\mathsf T} \left(\boldsymbol{\phi}_2-a_1\boldsymbol{\phi}_1\right) +\boldsymbol{\phi}_1^{\mathsf T} \dot{\hat{\boldsymbol{\theta}}}. \end{aligned} \]

The effective stage-2 regressor is therefore

\[ \boldsymbol{\omega}_2 =\boldsymbol{\phi}_2-a_1\boldsymbol{\phi}_1. \]

This formula is the simplest example of regressor propagation. The second plant regressor is modified by the sensitivity of the previous virtual control to the preceding state.

Add the stage-2 request to the existing tuning function:

\[ \boldsymbol{\tau}_2 =\boldsymbol{\tau}_1+ \boldsymbol{\Gamma}\boldsymbol{\omega}_2z_2 =\boldsymbol{\Gamma} \left(\boldsymbol{\phi}_1z_1+ \boldsymbol{\omega}_2z_2\right). \]

Because this is the final backstepping stage, set the actual adaptive law equal to the final tuning function:

\[ \dot{\hat{\boldsymbol{\theta}}}=\boldsymbol{\tau}_2. \]

Select the physical control input as

\[ \begin{aligned} u={}&a_1x_2+c_1\dot{y}_d+\ddot{y}_d-c_2z_2-z_1\\ &-\hat{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_2 -\boldsymbol{\phi}_1^{\mathsf T}\boldsymbol{\tau}_2. \end{aligned} \]

The last term is not optional. It compensates the explicit appearance of \( \dot{\hat{\boldsymbol{\theta}}} \) in \( \dot{\alpha}_1 \). Substituting the controller and adaptive law produces

\[ \dot{z}_2=-c_2z_2-z_1+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_2. \]

Use the composite candidate

\[ V_2=\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}}. \]

Its derivative along the transformed dynamics is

\[ \begin{aligned} \dot{V}_2={}&z_1\left(-c_1z_1+z_2+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1\right)\\ &+z_2\left(-c_2z_2-z_1+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_2\right)\\ &-\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\dot{\hat{\boldsymbol{\theta}}}. \end{aligned} \]

The terms \( z_1z_2 \) and \( -z_1z_2 \) cancel. Grouping all terms that multiply the parameter error gives

\[ \dot{V}_2=-c_1z_1^2-c_2z_2^2+ \tilde{\boldsymbol{\theta}}^{\mathsf T} \left( \boldsymbol{\phi}_1z_1+ \boldsymbol{\omega}_2z_2- \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}} \right). \]

Using \( \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma} (\boldsymbol{\phi}_1z_1+\boldsymbol{\omega}_2z_2) \) annihilates the entire bracket. Thus

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

It follows that \( z_1 \), \( z_2 \), and \( \tilde{\boldsymbol{\theta}} \) are bounded. Under the smoothness and bounded-reference assumptions, the closed-loop derivatives are bounded; the standard Barbalat argument introduced earlier in the course then yields convergence of the tracking coordinates to zero. Parameter convergence is a separate issue and requires sufficient excitation of the combined regressor.

For the general strict-feedback chain, define

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

At each stage, differentiate \( z_i \), substitute all previous virtual controls, and collect the coefficient of the common parameter error. Denote the resulting propagated regressor by \( \boldsymbol{\omega}_i \), with \( \boldsymbol{\omega}_1=\boldsymbol{\phi}_1 \). The tuning functions are accumulated as

\[ \boldsymbol{\tau}_0=\mathbf{0},\qquad \boldsymbol{\tau}_i =\boldsymbol{\tau}_{i-1} +\boldsymbol{\Gamma}\boldsymbol{\omega}_iz_i, \quad i=1,\ldots,n. \]

The virtual-control recursion is chosen so that, after compensating all derivatives of earlier virtual controls, the transformed coordinates have the target structure

\[ \begin{aligned} \dot{z}_1&=-c_1z_1+z_2+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_1,\\ \dot{z}_i&=-c_iz_i-z_{i-1}+z_{i+1}+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_i, \quad i=2,\ldots,n-1,\\ \dot{z}_n&=-c_nz_n-z_{n-1}+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\omega}_n. \end{aligned} \]

The detailed algebra that creates the compensating terms depends on the nonlinear regressors, but the invariant above is the design target. With

\[ V_n=\frac{1}{2}\sum_{i=1}^{n}z_i^2+ \frac{1}{2}\tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1}\tilde{\boldsymbol{\theta}}, \qquad \dot{\hat{\boldsymbol{\theta}}}=\boldsymbol{\tau}_n, \]

the chain cross terms telescope:

\[ \sum_{i=1}^{n-1}z_iz_{i+1} -\sum_{i=2}^{n}z_iz_{i-1}=0. \]

The remaining parameter terms also cancel because

\[ \boldsymbol{\Gamma}^{-1}\boldsymbol{\tau}_n =\sum_{i=1}^{n}\boldsymbol{\omega}_iz_i. \]

Consequently, the terminal derivative has the recursive invariant

\[ \boxed{ \dot{V}_n=-\sum_{i=1}^{n}c_iz_i^2\leq 0.} \]

1. Write every uncertain term with the same parameter vector \( \boldsymbol{\theta} \).
2. Define \( z_1=x_1-y_d \) and select \( \alpha_1 \) using certainty-equivalent cancellation.
3. Differentiate the current virtual control using the total derivative, including its dependence on \( \hat{\boldsymbol{\theta}} \).
4. Collect all unknown-parameter coefficients into one propagated regressor \( \boldsymbol{\omega}_i \).
5. Extend the tuning function by \( \boldsymbol{\Gamma}\boldsymbol{\omega}_iz_i \).
6. Insert compensation for every occurrence of a provisional adaptive update in the differentiated virtual controls.
7. Continue until the physical input stage, then impose \( \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\tau}_n \).
8. Verify the final Lyapunov derivative by explicitly checking both telescoping state-error products and parameter-error cancellation.

flowchart TD
  P["Plant regressors phi_1 ... phi_n"] --> R["Recursive regressor propagation"]
  R --> W["omega_1 ... omega_n"]
  W --> T["tau_i = tau_(i-1) + Gamma omega_i zi"]
  T --> U["Final adaptive law theta_hat_dot = tau_n"]
  V["Virtual controls alpha_1 ... alpha_(n-1)"] --> R
  U --> C["Compensation terms in final control"]
  C --> L["V_dot = -sum ci zi^2"]
        

The code laboratories use

\[ \boldsymbol{\theta}= \begin{bmatrix}0.8\\-0.6\end{bmatrix},\qquad \boldsymbol{\phi}_1= \begin{bmatrix}x_1\\\sin x_1\end{bmatrix},\qquad \boldsymbol{\phi}_2= \begin{bmatrix}x_1x_2\\x_2\end{bmatrix}. \]

The reference and derivatives are

\[ y_d=0.8\sin(0.5t),\qquad \dot{y}_d=0.4\cos(0.5t),\qquad \ddot{y}_d=-0.2\sin(0.5t). \]

The selected design parameters are

\[ c_1=2,\qquad c_2=3,\qquad \boldsymbol{\Gamma}=\operatorname{diag}(1.2,0.8). \]

All programs integrate the plant and adaptive states with a fourth-order Runge–Kutta method, save the time history to Chapter17_Lesson3_output.csv, and report the root-mean-square tracking error. The parameter estimates are not expected to equal the true parameters unless the combined regressor is persistently exciting. Good tracking with nonconvergent estimates is therefore not a contradiction.

This implementation uses NumPy for vector operations and an optional Matplotlib visualization. The controller and plant derivatives are evaluated together so that the Runge–Kutta stages remain dynamically consistent.

Chapter17_Lesson3.py

"""Chapter17_Lesson3.py

Adaptive backstepping with recursive tuning functions for a second-order
parametric strict-feedback system.

Dependencies:
    numpy
Optional:
    matplotlib
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from typing import Tuple

import numpy as np


@dataclass(frozen=True)
class Design:
    theta_true: np.ndarray
    gamma: np.ndarray
    c1: float = 2.0
    c2: float = 3.0


def reference(t: float) -> Tuple[float, float, float]:
    """Return y_d, first derivative, and second derivative."""
    yd = 0.8 * np.sin(0.5 * t)
    yd_dot = 0.4 * np.cos(0.5 * t)
    yd_ddot = -0.2 * np.sin(0.5 * t)
    return float(yd), float(yd_dot), float(yd_ddot)


def closed_loop_rhs(
    t: float, state: np.ndarray, design: Design
) -> Tuple[np.ndarray, Tuple[float, float, float]]:
    """Evaluate plant and adaptation dynamics.

    State ordering is [x1, x2, theta_hat_1, theta_hat_2].
    """
    x1, x2, theta_hat_1, theta_hat_2 = state
    theta_hat = np.array([theta_hat_1, theta_hat_2], dtype=float)

    yd, yd_dot, yd_ddot = reference(t)

    phi1 = np.array([x1, np.sin(x1)], dtype=float)
    dphi1_dx1 = np.array([1.0, np.cos(x1)], dtype=float)
    phi2 = np.array([x1 * x2, x2], dtype=float)

    z1 = x1 - yd
    alpha1 = -design.c1 * z1 - theta_hat @ phi1 + yd_dot
    z2 = x2 - alpha1

    # Partial derivative of alpha1 with respect to x1.
    a1 = -design.c1 - theta_hat @ dphi1_dx1

    # Regressor propagated from step 1 into step 2.
    omega2 = phi2 - a1 * phi1

    # Recursive tuning function and terminal adaptive law.
    tau2 = design.gamma @ (phi1 * z1 + omega2 * z2)
    theta_hat_dot = tau2

    # The term -phi1^T tau2 compensates the dependence of alpha1_dot on
    # theta_hat_dot. This is essential for the Lyapunov cancellation.
    u = (
        a1 * x2
        + design.c1 * yd_dot
        + yd_ddot
        - design.c2 * z2
        - z1
        - theta_hat @ omega2
        - phi1 @ tau2
    )

    x1_dot = x2 + design.theta_true @ phi1
    x2_dot = u + design.theta_true @ phi2

    derivative = np.array(
        [x1_dot, x2_dot, theta_hat_dot[0], theta_hat_dot[1]], dtype=float
    )
    return derivative, (float(u), float(z1), float(z2))


def rk4_step(t: float, state: np.ndarray, dt: float, design: Design) -> np.ndarray:
    k1, _ = closed_loop_rhs(t, state, design)
    k2, _ = closed_loop_rhs(t + 0.5 * dt, state + 0.5 * dt * k1, design)
    k3, _ = closed_loop_rhs(t + 0.5 * dt, state + 0.5 * dt * k2, design)
    k4, _ = closed_loop_rhs(t + dt, state + dt * k3, design)
    return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def main() -> None:
    design = Design(
        theta_true=np.array([0.8, -0.6], dtype=float),
        gamma=np.diag([1.2, 0.8]),
    )

    dt = 0.002
    final_time = 20.0
    steps = int(round(final_time / dt))
    state = np.array([1.0, -0.5, 0.0, 0.0], dtype=float)

    rows = []
    squared_tracking_error = 0.0

    for k in range(steps + 1):
        t = k * dt
        _, (u, z1, z2) = closed_loop_rhs(t, state, design)
        yd, _, _ = reference(t)
        rows.append([t, *state.tolist(), yd, u, z1, z2])
        squared_tracking_error += z1 * z1

        if k < steps:
            state = rk4_step(t, state, dt, design)

    rms_error = np.sqrt(squared_tracking_error / (steps + 1))
    print(f"RMS tracking error: {rms_error:.6f}")
    print(f"Final parameter estimate: [{state[2]:.6f}, {state[3]:.6f}]")
    print("True parameters need not be recovered without persistent excitation.")

    output_name = "Chapter17_Lesson3_output.csv"
    with open(output_name, "w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(
            ["t", "x1", "x2", "theta_hat_1", "theta_hat_2", "yd", "u", "z1", "z2"]
        )
        writer.writerows(rows)
    print(f"Saved simulation data to {output_name}")

    try:
        import matplotlib.pyplot as plt

        data = np.asarray(rows)
        plt.figure()
        plt.plot(data[:, 0], data[:, 1], label="x1")
        plt.plot(data[:, 0], data[:, 5], "--", label="yd")
        plt.xlabel("Time (s)")
        plt.ylabel("Output")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        plt.figure()
        plt.plot(data[:, 0], data[:, 3], label="theta_hat_1")
        plt.plot(data[:, 0], data[:, 4], label="theta_hat_2")
        plt.xlabel("Time (s)")
        plt.ylabel("Parameter estimates")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()
    except ImportError:
        print("matplotlib is not installed; CSV output was still generated.")


if __name__ == "__main__":
    main()

The C++17 program uses fixed-size std::array objects and the standard library only. This makes the numerical structure transparent and suitable for later migration to embedded or real-time frameworks.

Chapter17_Lesson3.cpp

// Chapter17_Lesson3.cpp
// Adaptive backstepping with recursive tuning functions.
// Build: g++ -std=c++17 -O2 Chapter17_Lesson3.cpp -o Chapter17_Lesson3

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

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

struct Design {
    Vec2 theta_true{0.8, -0.6};
    Vec2 gamma_diagonal{1.2, 0.8};
    double c1{2.0};
    double c2{3.0};
};

struct Auxiliary {
    double u{};
    double z1{};
    double z2{};
    double yd{};
};

double dot(const Vec2& a, const Vec2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

std::array<double, 3> reference(double t) {
    return {
        0.8 * std::sin(0.5 * t),
        0.4 * std::cos(0.5 * t),
        -0.2 * std::sin(0.5 * t)
    };
}

std::pair<State, Auxiliary> closed_loop_rhs(
    double t, const State& state, const Design& design
) {
    const double x1 = state[0];
    const double x2 = state[1];
    const Vec2 theta_hat{state[2], state[3]};

    const auto ref = reference(t);
    const double yd = ref[0];
    const double yd_dot = ref[1];
    const double yd_ddot = ref[2];

    const Vec2 phi1{x1, std::sin(x1)};
    const Vec2 dphi1_dx1{1.0, std::cos(x1)};
    const Vec2 phi2{x1 * x2, x2};

    const double z1 = x1 - yd;
    const double alpha1 = -design.c1 * z1 - dot(theta_hat, phi1) + yd_dot;
    const double z2 = x2 - alpha1;

    const double a1 = -design.c1 - dot(theta_hat, dphi1_dx1);
    const Vec2 omega2{
        phi2[0] - a1 * phi1[0],
        phi2[1] - a1 * phi1[1]
    };

    const Vec2 tau2{
        design.gamma_diagonal[0] * (phi1[0] * z1 + omega2[0] * z2),
        design.gamma_diagonal[1] * (phi1[1] * z1 + omega2[1] * z2)
    };

    const double u =
        a1 * x2
        + design.c1 * yd_dot
        + yd_ddot
        - design.c2 * z2
        - z1
        - dot(theta_hat, omega2)
        - dot(phi1, tau2);

    const double x1_dot = x2 + dot(design.theta_true, phi1);
    const double x2_dot = u + dot(design.theta_true, phi2);

    return {
        State{x1_dot, x2_dot, tau2[0], tau2[1]},
        Auxiliary{u, z1, z2, yd}
    };
}

State add_scaled(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 rk4_step(double t, const State& state, double dt, const Design& design) {
    const State k1 = closed_loop_rhs(t, state, design).first;
    const State k2 = closed_loop_rhs(
        t + 0.5 * dt, add_scaled(state, k1, 0.5 * dt), design
    ).first;
    const State k3 = closed_loop_rhs(
        t + 0.5 * dt, add_scaled(state, k2, 0.5 * dt), design
    ).first;
    const State k4 = closed_loop_rhs(
        t + dt, add_scaled(state, k3, dt), design
    ).first;

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

int main() {
    try {
        const Design design{};
        const double dt = 0.002;
        const double final_time = 20.0;
        const int steps = static_cast<int>(std::lround(final_time / dt));
        State state{1.0, -0.5, 0.0, 0.0};

        std::ofstream output("Chapter17_Lesson3_output.csv");
        if (!output) {
            throw std::runtime_error("Could not create output CSV file.");
        }
        output << "t,x1,x2,theta_hat_1,theta_hat_2,yd,u,z1,z2\n";
        output << std::setprecision(12);

        double squared_tracking_error = 0.0;
        for (int k = 0; k <= steps; ++k) {
            const double t = k * dt;
            const auto evaluation = closed_loop_rhs(t, state, design);
            const Auxiliary& aux = evaluation.second;

            output << t << ',' << state[0] << ',' << state[1] << ','
                   << state[2] << ',' << state[3] << ',' << aux.yd << ','
                   << aux.u << ',' << aux.z1 << ',' << aux.z2 << '\n';

            squared_tracking_error += aux.z1 * aux.z1;
            if (k < steps) {
                state = rk4_step(t, state, dt, design);
            }
        }

        const double rms_error = std::sqrt(
            squared_tracking_error / static_cast<double>(steps + 1)
        );
        std::cout << std::fixed << std::setprecision(6)
                  << "RMS tracking error: " << rms_error << '\n'
                  << "Final parameter estimate: [" << state[2] << ", "
                  << state[3] << "]\n"
                  << "True parameters need not be recovered without persistent excitation.\n"
                  << "Saved simulation data to Chapter17_Lesson3_output.csv\n";
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
    return 0;
}

The Java implementation uses primitive arrays, deterministic fixed-step integration, and buffered CSV output. No external numerical package is required for this low-dimensional example.

Chapter17_Lesson3.java

// Chapter17_Lesson3.java
// Adaptive backstepping with recursive tuning functions.
// Build and run: javac Chapter17_Lesson3.java && java Chapter17_Lesson3

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

public final class Chapter17_Lesson3 {
    private static final double[] THETA_TRUE = {0.8, -0.6};
    private static final double[] GAMMA_DIAGONAL = {1.2, 0.8};
    private static final double C1 = 2.0;
    private static final double C2 = 3.0;

    private Chapter17_Lesson3() {
    }

    private static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1];
    }

    private static double[] reference(double t) {
        return new double[] {
            0.8 * Math.sin(0.5 * t),
            0.4 * Math.cos(0.5 * t),
            -0.2 * Math.sin(0.5 * t)
        };
    }

    private static Evaluation closedLoopRhs(double t, double[] state) {
        double x1 = state[0];
        double x2 = state[1];
        double[] thetaHat = {state[2], state[3]};

        double[] ref = reference(t);
        double yd = ref[0];
        double ydDot = ref[1];
        double ydDdot = ref[2];

        double[] phi1 = {x1, Math.sin(x1)};
        double[] dphi1Dx1 = {1.0, Math.cos(x1)};
        double[] phi2 = {x1 * x2, x2};

        double z1 = x1 - yd;
        double alpha1 = -C1 * z1 - dot(thetaHat, phi1) + ydDot;
        double z2 = x2 - alpha1;

        double a1 = -C1 - dot(thetaHat, dphi1Dx1);
        double[] omega2 = {
            phi2[0] - a1 * phi1[0],
            phi2[1] - a1 * phi1[1]
        };

        double[] tau2 = {
            GAMMA_DIAGONAL[0] * (phi1[0] * z1 + omega2[0] * z2),
            GAMMA_DIAGONAL[1] * (phi1[1] * z1 + omega2[1] * z2)
        };

        double u = a1 * x2
            + C1 * ydDot
            + ydDdot
            - C2 * z2
            - z1
            - dot(thetaHat, omega2)
            - dot(phi1, tau2);

        double x1Dot = x2 + dot(THETA_TRUE, phi1);
        double x2Dot = u + dot(THETA_TRUE, phi2);

        return new Evaluation(
            new double[] {x1Dot, x2Dot, tau2[0], tau2[1]},
            u,
            z1,
            z2,
            yd
        );
    }

    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[] state, double dt) {
        double[] k1 = closedLoopRhs(t, state).derivative;
        double[] k2 = closedLoopRhs(t + 0.5 * dt, addScaled(state, k1, 0.5 * dt)).derivative;
        double[] k3 = closedLoopRhs(t + 0.5 * dt, addScaled(state, k2, 0.5 * dt)).derivative;
        double[] k4 = closedLoopRhs(t + dt, addScaled(state, k3, dt)).derivative;

        double[] next = new double[state.length];
        for (int i = 0; i < state.length; i++) {
            next[i] = state[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 dt = 0.002;
        double finalTime = 20.0;
        int steps = (int) Math.round(finalTime / dt);
        double[] state = {1.0, -0.5, 0.0, 0.0};
        double squaredTrackingError = 0.0;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter17_Lesson3_output.csv"))) {
            writer.write("t,x1,x2,theta_hat_1,theta_hat_2,yd,u,z1,z2\n");

            for (int k = 0; k <= steps; k++) {
                double t = k * dt;
                Evaluation evaluation = closedLoopRhs(t, state);
                writer.write(String.format(
                    Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    t,
                    state[0],
                    state[1],
                    state[2],
                    state[3],
                    evaluation.yd,
                    evaluation.u,
                    evaluation.z1,
                    evaluation.z2
                ));

                squaredTrackingError += evaluation.z1 * evaluation.z1;
                if (k < steps) {
                    state = rk4Step(t, state, dt);
                }
            }
        } catch (IOException error) {
            System.err.println("Could not write output CSV: " + error.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(squaredTrackingError / (steps + 1.0));
        System.out.printf(Locale.US, "RMS tracking error: %.6f%n", rmsError);
        System.out.printf(
            Locale.US,
            "Final parameter estimate: [%.6f, %.6f]%n",
            state[2],
            state[3]
        );
        System.out.println("True parameters need not be recovered without persistent excitation.");
        System.out.println("Saved simulation data to Chapter17_Lesson3_output.csv");
    }

    private static final class Evaluation {
        private final double[] derivative;
        private final double u;
        private final double z1;
        private final double z2;
        private final double yd;

        private Evaluation(double[] derivative, double u, double z1, double z2, double yd) {
            this.derivative = derivative;
            this.u = u;
            this.z1 = z1;
            this.z2 = z2;
            this.yd = yd;
        }
    }
}

The MATLAB script contains local functions for the reference, closed-loop derivative, and RK4 step. In Simulink, the same equations can be placed in a MATLAB Function block with integrators for \( x_1 \), \( x_2 \), and the two parameter estimates. The tuning function is the derivative input to the parameter-estimate integrators.

Chapter17_Lesson3.m

% Chapter17_Lesson3.m
% Adaptive backstepping with recursive tuning functions.
% Requires only base MATLAB. Local functions require MATLAB R2016b or newer.

clear; clc; close all;

parameters.thetaTrue = [0.8; -0.6];
parameters.Gamma = diag([1.2, 0.8]);
parameters.c1 = 2.0;
parameters.c2 = 3.0;

dt = 0.002;
finalTime = 20.0;
steps = round(finalTime / dt);
state = [1.0; -0.5; 0.0; 0.0];

data = zeros(steps + 1, 9);
squaredTrackingError = 0.0;

for k = 1:(steps + 1)
    t = (k - 1) * dt;
    [~, auxiliary] = closedLoopRhs(t, state, parameters);
    data(k, :) = [t, state.', auxiliary.yd, auxiliary.u, auxiliary.z1, auxiliary.z2];
    squaredTrackingError = squaredTrackingError + auxiliary.z1^2;

    if k <= steps
        state = rk4Step(t, state, dt, parameters);
    end
end

rmsError = sqrt(squaredTrackingError / (steps + 1));
fprintf('RMS tracking error: %.6f\n', rmsError);
fprintf('Final parameter estimate: [%.6f, %.6f]\n', state(3), state(4));
fprintf('True parameters need not be recovered without persistent excitation.\n');

variableNames = {'t', 'x1', 'x2', 'theta_hat_1', 'theta_hat_2', ...
                 'yd', 'u', 'z1', 'z2'};
resultTable = array2table(data, 'VariableNames', variableNames);
writetable(resultTable, 'Chapter17_Lesson3_output.csv');
fprintf('Saved simulation data to Chapter17_Lesson3_output.csv\n');

figure;
plot(data(:, 1), data(:, 2), 'LineWidth', 1.2);
hold on;
plot(data(:, 1), data(:, 6), '--', 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Output');
legend('x_1', 'y_d', 'Location', 'best');
title('Adaptive Backstepping Tracking');

figure;
plot(data(:, 1), data(:, 4), 'LineWidth', 1.2);
hold on;
plot(data(:, 1), data(:, 5), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Parameter estimates');
legend('\theta hat_1', '\theta hat_2', 'Location', 'best');
title('Adaptive Parameter Estimates');

function [derivative, auxiliary] = closedLoopRhs(t, state, parameters)
    x1 = state(1);
    x2 = state(2);
    thetaHat = state(3:4);

    [yd, ydDot, ydDdot] = referenceSignal(t);

    phi1 = [x1; sin(x1)];
    dphi1Dx1 = [1.0; cos(x1)];
    phi2 = [x1 * x2; x2];

    z1 = x1 - yd;
    alpha1 = -parameters.c1 * z1 - thetaHat.' * phi1 + ydDot;
    z2 = x2 - alpha1;

    a1 = -parameters.c1 - thetaHat.' * dphi1Dx1;
    omega2 = phi2 - a1 * phi1;

    tau2 = parameters.Gamma * (phi1 * z1 + omega2 * z2);
    thetaHatDot = tau2;

    u = a1 * x2 ...
        + parameters.c1 * ydDot ...
        + ydDdot ...
        - parameters.c2 * z2 ...
        - z1 ...
        - thetaHat.' * omega2 ...
        - phi1.' * tau2;

    x1Dot = x2 + parameters.thetaTrue.' * phi1;
    x2Dot = u + parameters.thetaTrue.' * phi2;

    derivative = [x1Dot; x2Dot; thetaHatDot];
    auxiliary = struct('u', u, 'z1', z1, 'z2', z2, 'yd', yd);
end

function nextState = rk4Step(t, state, dt, parameters)
    k1 = closedLoopRhs(t, state, parameters);
    k2 = closedLoopRhs(t + 0.5 * dt, state + 0.5 * dt * k1, parameters);
    k3 = closedLoopRhs(t + 0.5 * dt, state + 0.5 * dt * k2, parameters);
    k4 = closedLoopRhs(t + dt, state + dt * k3, parameters);
    nextState = state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
end

function [yd, ydDot, ydDdot] = referenceSignal(t)
    yd = 0.8 * sin(0.5 * t);
    ydDot = 0.4 * cos(0.5 * t);
    ydDdot = -0.2 * sin(0.5 * t);
end

The notebook file loads a companion Wolfram Language source file from the same directory. The source implements the same fixed-step recursion, exports the same CSV columns, and uses ListLinePlot for the tracking and parameter histories.

Chapter17_Lesson3.nb

Notebook[{
  Cell["Chapter 17, Lesson 3: Recursive Design of Adaptive Laws in Backstepping", "Title"],
  Cell["Evaluate the input cell below. It loads the companion Wolfram Language implementation from the same directory.", "Text"],
  Cell[
    BoxData[
      RowBox[{
        "Get", "[",
        RowBox[{
          "FileNameJoin", "[",
          RowBox[{
            "{",
            RowBox[{
              RowBox[{"NotebookDirectory", "[", "]"}],
              ",",
              "\"Chapter17_Lesson3.wl\""
            }],
            "}"
          }],
          "]"
        }],
        "]"
      }]
    ],
    "Input"
  ]
},
WindowTitle -> "Chapter17_Lesson3",
StyleDefinitions -> "Default.nb"
]

Chapter17_Lesson3.wl

(* Chapter17_Lesson3.wl
   Adaptive backstepping with recursive tuning functions.
   Wolfram Language implementation loaded by Chapter17_Lesson3.nb.
*)

ClearAll["Global`*"];

thetaTrue = {0.8, -0.6};
gamma = DiagonalMatrix[{1.2, 0.8}];
c1 = 2.0;
c2 = 3.0;

reference[t_?NumericQ] := {
  0.8 Sin[0.5 t],
  0.4 Cos[0.5 t],
  -0.2 Sin[0.5 t]
};

closedLoopRhs[t_?NumericQ, state_List] := Module[
  {
    x1, x2, thetaHat, yd, ydDot, ydDdot, phi1, dphi1Dx1,
    phi2, z1, alpha1, z2, a1, omega2, tau2, u, x1Dot, x2Dot
  },
  x1 = state[[1]];
  x2 = state[[2]];
  thetaHat = state[[3 ;; 4]];
  {yd, ydDot, ydDdot} = reference[t];

  phi1 = {x1, Sin[x1]};
  dphi1Dx1 = {1.0, Cos[x1]};
  phi2 = {x1 x2, x2};

  z1 = x1 - yd;
  alpha1 = -c1 z1 - thetaHat . phi1 + ydDot;
  z2 = x2 - alpha1;

  a1 = -c1 - thetaHat . dphi1Dx1;
  omega2 = phi2 - a1 phi1;

  tau2 = gamma . (phi1 z1 + omega2 z2);

  u = a1 x2 + c1 ydDot + ydDdot - c2 z2 - z1
      - thetaHat . omega2 - phi1 . tau2;

  x1Dot = x2 + thetaTrue . phi1;
  x2Dot = u + thetaTrue . phi2;

  {{x1Dot, x2Dot, tau2[[1]], tau2[[2]]}, {u, z1, z2, yd}}
];

rk4Step[t_?NumericQ, state_List, dt_?NumericQ] := Module[
  {k1, k2, k3, k4},
  k1 = First[closedLoopRhs[t, state]];
  k2 = First[closedLoopRhs[t + dt/2.0, state + dt k1/2.0]];
  k3 = First[closedLoopRhs[t + dt/2.0, state + dt k2/2.0]];
  k4 = First[closedLoopRhs[t + dt, state + dt k3]];
  state + dt (k1 + 2.0 k2 + 2.0 k3 + k4)/6.0
];

dt = 0.002;
finalTime = 20.0;
steps = Round[finalTime/dt];
state = {1.0, -0.5, 0.0, 0.0};

simulationData = Reap[
  Do[
    t = k dt;
    auxiliary = Last[closedLoopRhs[t, state]];
    Sow[Join[{t}, state, {auxiliary[[4]], auxiliary[[1]], auxiliary[[2]], auxiliary[[3]]}]];
    If[k < steps, state = rk4Step[t, state, dt]],
    {k, 0, steps}
  ]
][[2, 1]];

headers = {
  "t", "x1", "x2", "theta_hat_1", "theta_hat_2",
  "yd", "u", "z1", "z2"
};
Export["Chapter17_Lesson3_output.csv", Prepend[simulationData, headers]];

rmsError = Sqrt[Mean[simulationData[[All, 8]]^2]];
Print["RMS tracking error: ", NumberForm[rmsError, {8, 6}]];
Print["Final parameter estimate: ", NumberForm[state[[3 ;; 4]], {8, 6}]];
Print["True parameters need not be recovered without persistent excitation."];

trackingPlot = ListLinePlot[
  {
    simulationData[[All, {1, 2}]],
    simulationData[[All, {1, 6}]]
  },
  PlotLegends -> {"x1", "yd"},
  Frame -> True,
  FrameLabel -> {"Time (s)", "Output"},
  PlotRange -> All,
  ImageSize -> Large
];

parameterPlot = ListLinePlot[
  {
    simulationData[[All, {1, 4}]],
    simulationData[[All, {1, 5}]]
  },
  PlotLegends -> {"thetaHat1", "thetaHat2"},
  Frame -> True,
  FrameLabel -> {"Time (s)", "Parameter estimates"},
  PlotRange -> All,
  ImageSize -> Large
];

Column[{trackingPlot, parameterPlot}]

17.1 Use the Total Derivative

Omitting \( (\partial\alpha_i/ \partial\hat{\boldsymbol{\theta}}) \dot{\hat{\boldsymbol{\theta}}} \) changes the actual closed-loop system and invalidates the Lyapunov cancellation.

17.2 Do Not Update Intermediate Tuning Functions as Separate States

Only \( \hat{\boldsymbol{\theta}} \) is integrated. The signals \( \boldsymbol{\tau}_1,\ldots, \boldsymbol{\tau}_n \) are algebraic design quantities.

17.3 Keep Dimensions Explicit

Each \( \boldsymbol{\omega}_i \) and \( \boldsymbol{\tau}_i \) lies in \( \mathbb{R}^{p} \), while each \( z_i \) is scalar in the SISO strict-feedback chain. Dimension checks catch many transcription errors.

17.4 Distinguish Tracking from Identification

The negative-semidefinite Lyapunov derivative establishes tracking-error convergence under the stated assumptions, but it does not by itself imply \( \tilde{\boldsymbol{\theta}}=\mathbf{0} \). Parameter convergence requires the excitation conditions studied earlier.

17.5 Numerical Gain Selection

Larger \( c_i \) values increase nominal error damping, while larger entries of \( \boldsymbol{\Gamma} \) speed adaptation but amplify sensitivity to discretization and measurement noise. The ideal continuous-time proof should therefore be accompanied by a sampling-rate and numerical-stability check.

Problem 1 — First-stage transformed dynamics. For \( \dot{x}_1=x_2+ \boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_1(x_1) \), define \( z_1=x_1-y_d \) and \( \alpha_1=-c_1z_1- \hat{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1+ \dot{y}_d \). Derive \( \dot{z}_1 \) in terms of \( z_1,z_2 \), and the parameter error.

Solution: Since \( x_2=z_2+\alpha_1 \),

\[ \begin{aligned} \dot{z}_1 &=x_2+\boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_1- \dot{y}_d\\ &=z_2-c_1z_1- \hat{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1+ \boldsymbol{\theta}^{\mathsf T}\boldsymbol{\phi}_1\\ &=-c_1z_1+z_2+ \tilde{\boldsymbol{\theta}}^{\mathsf T}\boldsymbol{\phi}_1. \end{aligned} \]

Problem 2 — Propagated regressor. Let \( \boldsymbol{\phi}_1=[x_1,\sin x_1]^{\mathsf T} \) and \( \boldsymbol{\phi}_2=[x_1x_2,x_2]^{\mathsf T} \). Compute \( a_1 \) and \( \boldsymbol{\omega}_2 \) explicitly.

Solution:

\[ \frac{\partial\boldsymbol{\phi}_1}{\partial x_1} =\begin{bmatrix}1\\\cos x_1\end{bmatrix}, \qquad a_1=-c_1-\hat{\theta}_1-\hat{\theta}_2\cos x_1. \]

\[ \boldsymbol{\omega}_2= \begin{bmatrix} x_1x_2-a_1x_1\\ x_2-a_1\sin x_1 \end{bmatrix}. \]

Problem 3 — Necessity of the adaptive-derivative compensation. Suppose the controller omits \( -\boldsymbol{\phi}_1^{\mathsf T} \boldsymbol{\tau}_2 \). Determine the residual term in \( \dot{z}_2 \) after imposing \( \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\tau}_2 \).

Solution: Before control substitution,

\[ \dot{z}_2=\cdots+ \boldsymbol{\phi}_1^{\mathsf T} \dot{\hat{\boldsymbol{\theta}}}. \]

Without the compensating control term, substitution of the terminal update leaves \( \boldsymbol{\phi}_1^{\mathsf T} \boldsymbol{\tau}_2 \) in the state dynamics. This term has no definite sign in the Lyapunov derivative and destroys the stated cancellation proof.

Problem 4 — Full derivative cancellation. Starting from the two transformed error equations, prove that the proposed tuning law yields \( \dot{V}_2=-c_1z_1^2-c_2z_2^2 \).

Solution: Substitute the error equations into

\[ \dot{V}_2=z_1\dot{z}_1+z_2\dot{z}_2- \tilde{\boldsymbol{\theta}}^{\mathsf T} \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}. \]

The state cross products cancel pairwise. The remaining uncertain part is

\[ \tilde{\boldsymbol{\theta}}^{\mathsf T} \left(\boldsymbol{\phi}_1z_1+ \boldsymbol{\omega}_2z_2- \boldsymbol{\Gamma}^{-1} \dot{\hat{\boldsymbol{\theta}}}\right). \]

The terminal tuning law makes the bracket exactly zero, leaving only the two negative quadratic damping terms.

Problem 5 — Quantifying overparameterization. A fourth-order strict-feedback plant contains a common unknown vector \( \boldsymbol{\theta}\in\mathbb{R}^{3} \). Compare the number of adaptive states in a design with one independent estimate at each stage and a tuning-function design.

Solution: Independent stagewise estimates require \( 4\times 3=12 \) adaptive states. A tuning-function design updates one three-dimensional estimate and therefore requires only \( 3 \) adaptive states. The reduction is nine states, and all stages remain tied to the same physical parameter vector.

Problem 6 — Tracking without parameter convergence. Explain how \( z_1 \) and \( z_2 \) can converge to zero while \( \tilde{\boldsymbol{\theta}} \) approaches a nonzero constant.

Solution: The Lyapunov derivative penalizes the transformed tracking errors but contains no negative-definite term in \( \tilde{\boldsymbol{\theta}} \). Once the tracking errors become small, the adaptive update \( \boldsymbol{\Gamma} (\boldsymbol{\phi}_1z_1+\boldsymbol{\omega}_2z_2) \) also becomes small. Without persistent excitation, several parameter vectors can be compatible with the observed trajectory, so the estimate can stop at a bounded nontrue value while the controller still tracks successfully.

Adaptive backstepping is recursive in both the stabilizing controls and the parameter update. Differentiating a parameter-dependent virtual control propagates a new effective regressor and introduces the adaptive derivative into the next equation. Tuning functions accumulate each stage's requested cancellation while retaining one shared parameter estimate. At the terminal stage, the update law is set equal to the final tuning function, state-error cross terms telescope, and all parameter-error products cancel. The resulting design guarantees bounded closed-loop signals and asymptotic tracking under the ideal assumptions, while parameter convergence still depends on excitation.