In ordinary backstepping, the next state is temporarily treated as a control input for the preceding subsystem. The stabilizing value assigned to that state is called a virtual control. In adaptive backstepping, the subsystem may contain an unknown constant parameter, so the virtual control must normally use an estimate of that parameter.

Consider the first subsystem of a parametric strict-feedback plant,

\[ \dot{x}_1 = x_2 + \theta\,\phi_1(x_1), \qquad \theta \in \mathbb{R}\ \text{unknown and constant}. \]

If \( x_2 \) could be commanded directly, a natural virtual control would cancel the uncertain term by replacing \( \theta \) with an estimate \( \hat{\theta} \). Consequently, the virtual control is no longer a function of measured states and references alone; it is a function of the adaptive state as well:

\[ \alpha_1 = \alpha_1(x_1,r,\dot r,\hat{\theta}). \]

This dependence matters because the next recursive step requires \( \dot{\alpha}_1 \). The chain rule then introduces \( \dot{\hat{\theta}} \) and may reintroduce the true unknown parameter through \( \dot{x}_1 \). This is the central structural issue of the lesson.

flowchart TD
  R["Reference r, r_dot, r_ddot"] --> E1["z1 = x1 - r"]
  X1["Measured x1"] --> E1
  H1["Estimate theta_hat_1"] --> A1["Virtual control alpha1"]
  E1 --> A1
  R --> A1
  A1 --> E2["z2 = x2 - alpha1"]
  X2["Measured x2"] --> E2
  E2 --> U["Actual control u"]
  H2["Derivative estimate theta_hat_2"] --> U
  U --> P["Strict-feedback plant"]
  P --> X1
  P --> X2
        

To isolate parameter dependence without introducing unnecessary nonlinearities, consider the second-order strict-feedback system

\[ \begin{aligned} \dot{x}_1 &= x_2 + \theta x_1,\\ \dot{x}_2 &= u,\\ y &= x_1, \end{aligned} \qquad \theta \in \mathbb{R}\ \text{unknown}. \]

The desired trajectory \( r(t) \) is assumed bounded with known bounded derivatives \( \dot r(t) \) and \( \ddot r(t) \). The objective is to make \( x_1(t) \) track \( r(t) \) while keeping all closed-loop signals bounded. Define

\[ z_1 = x_1-r. \]

The true parameter is represented by two estimation errors later in the design:

\[ \tilde{\theta}_1 = \theta-\hat{\theta}_1, \qquad \tilde{\theta}_2 = \theta-\hat{\theta}_2. \]

Using two estimates for one physical parameter is called overparameterization. It is introduced here deliberately so the derivative of the parameter-dependent virtual control remains implementable and the Lyapunov cancellation is transparent.

Differentiating the first tracking error gives

\[ \dot z_1 = x_2 + \theta x_1-\dot r. \]

Treat \( x_2 \) as a virtual input and select

\[ \boxed{\alpha_1 = \dot r-k_1z_1-\hat{\theta}_1x_1}, \qquad k_1 > 0. \]

The second backstepping coordinate measures how far the physical state \( x_2 \) is from the desired virtual value:

\[ z_2=x_2-\alpha_1. \]

Substituting \( x_2=z_2+\alpha_1 \) into \( \dot z_1 \) yields

\[ \begin{aligned} \dot z_1 &=z_2+\dot r-k_1z_1-\hat{\theta}_1x_1+\theta x_1-\dot r\\ &=-k_1z_1+z_2+\tilde{\theta}_1x_1. \end{aligned} \]

For the first-stage Lyapunov function

\[ V_1=\frac{1}{2}z_1^2+\frac{1}{2\gamma_1}\tilde{\theta}_1^2, \qquad \gamma_1 > 0, \]

its derivative is

\[ \dot V_1=-k_1z_1^2+z_1z_2+\tilde{\theta}_1x_1z_1 -\frac{1}{\gamma_1}\tilde{\theta}_1\dot{\hat{\theta}}_1. \]

The update law

\[ \boxed{\dot{\hat{\theta}}_1=\gamma_1x_1z_1} \]

cancels the first parameter-error product and leaves

\[ \dot V_1=-k_1z_1^2+z_1z_2. \]

The remaining cross term \( z_1z_2 \) is expected. It signals that the virtual control cannot be imposed exactly and must be completed by the actual input design.

Rewrite the virtual control in a form convenient for differentiation:

\[ \alpha_1=\dot r+k_1r-(k_1+\hat{\theta}_1)x_1. \]

The exact derivative is

\[ \dot\alpha_1=\ddot r+k_1\dot r -(k_1+\hat{\theta}_1)\dot x_1 -\dot{\hat{\theta}}_1x_1. \]

Substitution of \( \dot x_1=x_2+\theta x_1 \) gives

\[ \dot\alpha_1=\ddot r+k_1\dot r -(k_1+\hat{\theta}_1)x_2 -(k_1+\hat{\theta}_1)\theta x_1 -\dot{\hat{\theta}}_1x_1. \]

This expression is not implementable because it contains the true unknown \( \theta \). Merely replacing it by \( \hat{\theta}_1 \) would create a residual product involving the same estimate that already appears inside the virtual control. The classical early adaptive-backstepping construction instead introduces a second estimate:

\[ \boxed{\widehat{\dot\alpha}_1=\ddot r+k_1\dot r -(k_1+\hat{\theta}_1)x_2 -(k_1+\hat{\theta}_1)\hat{\theta}_2x_1 -\dot{\hat{\theta}}_1x_1}. \]

The derivative mismatch is exactly

\[ \begin{aligned} \dot\alpha_1-\widehat{\dot\alpha}_1 &=-(k_1+\hat{\theta}_1)(\theta-\hat{\theta}_2)x_1\\ &=-(k_1+\hat{\theta}_1)\tilde{\theta}_2x_1. \end{aligned} \]

Thus, parameter dependence has been converted into a linearly appearing second parameter error. This structure can be canceled by another adaptation law.

Since \( z_2=x_2-\alpha_1 \), its dynamics satisfy

\[ \dot z_2=u-\dot\alpha_1. \]

Select the actual input

\[ \boxed{u=\widehat{\dot\alpha}_1-z_1-k_2z_2}, \qquad k_2 > 0. \]

Using the derivative mismatch from the previous section gives

\[ \dot z_2=-z_1-k_2z_2 +(k_1+\hat{\theta}_1)\tilde{\theta}_2x_1. \]

Extend the Lyapunov function by including the second error coordinate and second estimation error:

\[ V_2=\frac{1}{2}z_1^2+\frac{1}{2}z_2^2 +\frac{1}{2\gamma_1}\tilde{\theta}_1^2 +\frac{1}{2\gamma_2}\tilde{\theta}_2^2, \qquad \gamma_2 > 0. \]

The regressor multiplying \( \tilde{\theta}_2 \) in \( \dot V_2 \) is \( x_1(k_1+\hat{\theta}_1)z_2 \). Therefore choose

\[ \boxed{\dot{\hat{\theta}}_2= \gamma_2x_1(k_1+\hat{\theta}_1)z_2}. \]

The update law depends on the preceding estimate. This is not an accidental implementation detail; it is the algebraic footprint of the parameter dependence inherited from \( \alpha_1 \).

Differentiate \( V_2 \). Because \( \dot{\tilde{\theta}}_i=-\dot{\hat{\theta}}_i \),

\[ \begin{aligned} \dot V_2 &=z_1\dot z_1+z_2\dot z_2 -\frac{1}{\gamma_1}\tilde{\theta}_1\dot{\hat{\theta}}_1 -\frac{1}{\gamma_2}\tilde{\theta}_2\dot{\hat{\theta}}_2\\ &=z_1\left(-k_1z_1+z_2+\tilde{\theta}_1x_1\right)\\ &\quad+z_2\left[-z_1-k_2z_2 +(k_1+\hat{\theta}_1)\tilde{\theta}_2x_1\right]\\ &\quad-\frac{1}{\gamma_1}\tilde{\theta}_1\dot{\hat{\theta}}_1 -\frac{1}{\gamma_2}\tilde{\theta}_2\dot{\hat{\theta}}_2. \end{aligned} \]

First, the recursive cross terms cancel: \( z_1z_2-z_2z_1=0 \). Next, substitute the two update laws:

\[ \begin{aligned} \tilde{\theta}_1x_1z_1 -\frac{1}{\gamma_1}\tilde{\theta}_1 (\gamma_1x_1z_1)&=0,\\ \tilde{\theta}_2x_1(k_1+\hat{\theta}_1)z_2 -\frac{1}{\gamma_2}\tilde{\theta}_2 \left[\gamma_2x_1(k_1+\hat{\theta}_1)z_2\right]&=0. \end{aligned} \]

Therefore

\[ \boxed{\dot V_2=-k_1z_1^2-k_2z_2^2\leq 0}. \]

Hence \( V_2(t)\leq V_2(0) \), which proves boundedness of \( z_1,z_2,\tilde{\theta}_1,\tilde{\theta}_2 \). Bounded reference signals imply bounded states and controller signals. Moreover,

\[ \int_0^{\infty}\left(k_1z_1^2+k_2z_2^2\right)dt \leq V_2(0), \]

so \( z_1,z_2\in L_2 \). Under the stated smoothness and boundedness assumptions, their derivatives are bounded. Barbalat's lemma then gives

\[ z_1(t)→0,\qquad z_2(t)→0 \qquad \text{as}\quad t→\infty. \]

Tracking follows because \( z_1=x_1-r \). Parameter convergence is not implied by this proof; it requires sufficient excitation of the associated regressors, as studied earlier in the course.

The example reveals four general facts that persist in higher-order adaptive backstepping.

1. Adaptive states enter coordinate transformations. Because \( z_2=x_2-\alpha_1(\hat\theta_1) \), the error coordinates themselves depend on online estimates.
2. Differentiation propagates adaptation. Every derivative of a virtual control contains terms of the form \( \partial\alpha_i/\partial\hat\theta\cdot \dot{\hat\theta} \).
3. Unknown plant dynamics can reappear. Through \( \partial\alpha_i/\partial x_j\cdot\dot x_j \), the true parameter may enter the next step even though it was apparently canceled in the preceding step.
4. The Lyapunov function dictates the update regressor. The adaptation law is selected from the exact coefficient multiplying each parameter error in \( \dot V \), not from an independent estimation objective.

For an \( n \)-stage system, repeating the simple overparameterized construction can create many estimates of the same physical parameter. Later adaptive-backstepping methods introduce tuning functions so that a single estimate can be shared recursively. That recursive mechanism is developed in the next lesson.

flowchart TD
  S1["Choose alpha1 with theta_hat_1"] --> D1["Differentiate alpha1"]
  D1 --> C1["Chain-rule term: theta_hat_1_dot"]
  D1 --> U1["Unknown theta reappears through x1_dot"]
  U1 --> H2["Introduce theta_hat_2 in alpha1_dot_hat"]
  H2 --> Z2["Form z2 dynamics"]
  Z2 --> L2["Add z2 and parameter error to Lyapunov function"]
  L2 --> A2["Select theta_hat_2_dot from cancellation condition"]
  A2 --> N["Obtain V_dot = -k1 z1^2 - k2 z2^2"]
        

At each integration instant, evaluate the following equations in order:

\[ \begin{aligned} z_1 &= x_1-r,\\ \dot{\hat\theta}_1 &= \gamma_1x_1z_1,\\ \alpha_1 &= \dot r-k_1z_1-\hat\theta_1x_1,\\ z_2 &= x_2-\alpha_1,\\ \dot{\hat\theta}_2 &= \gamma_2x_1(k_1+\hat\theta_1)z_2,\\ \widehat{\dot\alpha}_1 &=\ddot r+k_1\dot r-(k_1+\hat\theta_1)x_2\\ &\quad-(k_1+\hat\theta_1)\hat\theta_2x_1 -\dot{\hat\theta}_1x_1,\\ u &= \widehat{\dot\alpha}_1-z_1-k_2z_2. \end{aligned} \]

The order is important: \( \dot{\hat\theta}_1 \) is needed inside \( \widehat{\dot\alpha}_1 \). In software, compute the controller algebraically from the current state and then pass the resulting derivatives to the numerical integrator.

The Python program uses NumPy for state vectors, Matplotlib for visualization, and a from-scratch RK4 solver so every closed-loop equation remains visible.

Chapter17_Lesson2.py

"""Chapter17_Lesson2.py

Adaptive backstepping demonstration for a two-state strict-feedback plant.
The example highlights how parameter dependence in the virtual control creates
an additional estimated derivative term and a second parameter estimate.
"""

from __future__ import annotations

from dataclasses import dataclass
import csv
import math
from pathlib import Path
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Settings:
    theta: float = 1.2
    k1: float = 2.0
    k2: float = 3.0
    gamma1: float = 4.0
    gamma2: float = 2.0
    t_final: float = 20.0
    dt: float = 0.001


def reference(t: float) -> tuple[float, float, float]:
    """Return r(t), r_dot(t), and r_ddot(t)."""
    omega = 0.7
    amplitude = 0.5
    r = amplitude * math.sin(omega * t)
    r_dot = amplitude * omega * math.cos(omega * t)
    r_ddot = -amplitude * omega**2 * math.sin(omega * t)
    return r, r_dot, r_ddot


def closed_loop_rhs(
    t: float,
    state: np.ndarray,
    cfg: Settings,
) -> tuple[np.ndarray, dict[str, float]]:
    """Evaluate the plant, controller, and adaptation dynamics."""
    x1, x2, theta_hat_1, theta_hat_2 = state
    r, r_dot, r_ddot = reference(t)

    z1 = x1 - r
    theta_hat_1_dot = cfg.gamma1 * x1 * z1

    alpha1 = r_dot - cfg.k1 * z1 - theta_hat_1 * x1
    z2 = x2 - alpha1

    theta_hat_2_dot = (
        cfg.gamma2 * x1 * (cfg.k1 + theta_hat_1) * z2
    )

    alpha1_dot_hat = (
        r_ddot
        + cfg.k1 * r_dot
        - (cfg.k1 + theta_hat_1) * x2
        - (cfg.k1 + theta_hat_1) * theta_hat_2 * x1
        - theta_hat_1_dot * x1
    )

    u = alpha1_dot_hat - z1 - cfg.k2 * z2

    x1_dot = x2 + cfg.theta * x1
    x2_dot = u

    derivative = np.array(
        [x1_dot, x2_dot, theta_hat_1_dot, theta_hat_2_dot],
        dtype=float,
    )
    outputs = {
        "r": r,
        "z1": z1,
        "z2": z2,
        "alpha1": alpha1,
        "alpha1_dot_hat": alpha1_dot_hat,
        "u": u,
    }
    return derivative, outputs


def rk4_step(
    rhs: Callable[[float, np.ndarray, Settings], tuple[np.ndarray, dict[str, float]]],
    t: float,
    state: np.ndarray,
    dt: float,
    cfg: Settings,
) -> np.ndarray:
    """Advance one classical fourth-order Runge-Kutta step."""
    k1, _ = rhs(t, state, cfg)
    k2, _ = rhs(t + 0.5 * dt, state + 0.5 * dt * k1, cfg)
    k3, _ = rhs(t + 0.5 * dt, state + 0.5 * dt * k2, cfg)
    k4, _ = rhs(t + dt, state + dt * k3, cfg)
    return state + (dt / 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.array([1.2, -0.2, 0.0, 0.0], dtype=float)

    data = {
        name: np.zeros(steps, dtype=float)
        for name in (
            "x1",
            "x2",
            "theta_hat_1",
            "theta_hat_2",
            "r",
            "z1",
            "z2",
            "alpha1",
            "u",
            "V",
        )
    }

    for index, t in enumerate(time):
        _, outputs = closed_loop_rhs(t, state, cfg)
        x1, x2, theta_hat_1, theta_hat_2 = state
        tilde_1 = cfg.theta - theta_hat_1
        tilde_2 = cfg.theta - theta_hat_2
        lyapunov = (
            0.5 * outputs["z1"] ** 2
            + 0.5 * outputs["z2"] ** 2
            + 0.5 * tilde_1**2 / cfg.gamma1
            + 0.5 * tilde_2**2 / cfg.gamma2
        )

        data["x1"][index] = x1
        data["x2"][index] = x2
        data["theta_hat_1"][index] = theta_hat_1
        data["theta_hat_2"][index] = theta_hat_2
        data["r"][index] = outputs["r"]
        data["z1"][index] = outputs["z1"]
        data["z2"][index] = outputs["z2"]
        data["alpha1"][index] = outputs["alpha1"]
        data["u"][index] = outputs["u"]
        data["V"][index] = lyapunov

        if index + 1 < steps:
            state = rk4_step(closed_loop_rhs, t, state, cfg.dt, cfg)

    data["time"] = time
    return data


def write_csv(data: dict[str, np.ndarray], path: Path) -> None:
    fields = [
        "time",
        "x1",
        "r",
        "x2",
        "z1",
        "z2",
        "theta_hat_1",
        "theta_hat_2",
        "alpha1",
        "u",
        "V",
    ]
    with path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(fields)
        for row in zip(*(data[field] for field in fields), strict=True):
            writer.writerow(row)


def plot_results(data: dict[str, np.ndarray]) -> None:
    time = data["time"]

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

    plt.figure()
    plt.plot(time, data["z1"], label="z1")
    plt.plot(time, data["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(time, data["theta_hat_1"], label="theta_hat_1")
    plt.plot(time, data["theta_hat_2"], label="theta_hat_2")
    plt.xlabel("Time [s]")
    plt.ylabel("Estimate")
    plt.title("Overparameterized estimates")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(time, data["u"])
    plt.xlabel("Time [s]")
    plt.ylabel("u")
    plt.title("Control input")
    plt.grid(True)

    plt.show()


def main() -> None:
    cfg = Settings()
    data = simulate(cfg)
    csv_path = Path("Chapter17_Lesson2_python_results.csv")
    write_csv(data, csv_path)

    print(f"Final z1 = {data['z1'][-1]: .6e}")
    print(f"Final z2 = {data['z2'][-1]: .6e}")
    print(f"Final theta_hat_1 = {data['theta_hat_1'][-1]: .6f}")
    print(f"Final theta_hat_2 = {data['theta_hat_2'][-1]: .6f}")
    print(f"Results written to {csv_path.resolve()}")

    plot_results(data)


if __name__ == "__main__":
    main()

The C++17 version uses only the standard library. It integrates the same four-state closed loop by RK4 and writes a CSV file suitable for plotting with Python, MATLAB, GNUplot, or spreadsheet software.

Chapter17_Lesson2.cpp

// Chapter17_Lesson2.cpp
// Adaptive backstepping with parameter-dependent virtual control.

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

namespace {

constexpr double kTheta = 1.2;
constexpr double kGain1 = 2.0;
constexpr double kGain2 = 3.0;
constexpr double kGamma1 = 4.0;
constexpr double kGamma2 = 2.0;
constexpr double kTimeStep = 0.001;
constexpr double kFinalTime = 20.0;

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

struct Outputs {
    double r{};
    double z1{};
    double z2{};
    double alpha1{};
    double u{};
};

std::array<double, 3> reference(const double t) {
    constexpr double amplitude = 0.5;
    constexpr double omega = 0.7;
    return {
        amplitude * std::sin(omega * t),
        amplitude * omega * std::cos(omega * t),
        -amplitude * omega * omega * std::sin(omega * t)
    };
}

State rhs(const double t, const State& state, Outputs* outputs = nullptr) {
    const auto ref = reference(t);
    const double r = ref[0];
    const double rDot = ref[1];
    const double rDDot = ref[2];

    const double x1 = state[0];
    const double x2 = state[1];
    const double thetaHat1 = state[2];
    const double thetaHat2 = state[3];

    const double z1 = x1 - r;
    const double thetaHat1Dot = kGamma1 * x1 * z1;
    const double alpha1 = rDot - kGain1 * z1 - thetaHat1 * x1;
    const double z2 = x2 - alpha1;
    const double thetaHat2Dot =
        kGamma2 * x1 * (kGain1 + thetaHat1) * z2;

    const double alpha1DotHat =
        rDDot + kGain1 * rDot
        - (kGain1 + thetaHat1) * x2
        - (kGain1 + thetaHat1) * thetaHat2 * x1
        - thetaHat1Dot * x1;

    const double u = alpha1DotHat - z1 - kGain2 * z2;

    if (outputs != nullptr) {
        outputs->r = r;
        outputs->z1 = z1;
        outputs->z2 = z2;
        outputs->alpha1 = alpha1;
        outputs->u = u;
    }

    return {
        x2 + kTheta * x1,
        u,
        thetaHat1Dot,
        thetaHat2Dot
    };
}

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

State rk4Step(const double t, const State& state, const double dt) {
    const State k1 = rhs(t, state);
    const State k2 = rhs(t + 0.5 * dt, addScaled(state, k1, 0.5 * dt));
    const State k3 = rhs(t + 0.5 * dt, addScaled(state, k2, 0.5 * dt));
    const State k4 = rhs(t + dt, addScaled(state, k3, dt));

    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;
}

}  // namespace

int main() {
    std::ofstream csv("Chapter17_Lesson2_cpp_results.csv");
    if (!csv) {
        throw std::runtime_error("Unable to create output CSV file.");
    }

    csv << "time,x1,r,x2,z1,z2,theta_hat_1,theta_hat_2,alpha1,u,V\n";
    csv << std::setprecision(12);

    State state{1.2, -0.2, 0.0, 0.0};
    const int steps = static_cast<int>(std::llround(kFinalTime / kTimeStep));

    for (int index = 0; index <= steps; ++index) {
        const double t = index * kTimeStep;
        Outputs outputs{};
        rhs(t, state, &outputs);

        const double tilde1 = kTheta - state[2];
        const double tilde2 = kTheta - state[3];
        const double lyapunov =
            0.5 * outputs.z1 * outputs.z1
            + 0.5 * outputs.z2 * outputs.z2
            + 0.5 * tilde1 * tilde1 / kGamma1
            + 0.5 * tilde2 * tilde2 / kGamma2;

        csv << t << ',' << state[0] << ',' << outputs.r << ',' << state[1]
            << ',' << outputs.z1 << ',' << outputs.z2 << ',' << state[2]
            << ',' << state[3] << ',' << outputs.alpha1 << ',' << outputs.u
            << ',' << lyapunov << '\n';

        if (index < steps) {
            state = rk4Step(t, state, kTimeStep);
        }
    }

    Outputs finalOutputs{};
    rhs(kFinalTime, state, &finalOutputs);
    std::cout << std::scientific << std::setprecision(6)
              << "Final z1 = " << finalOutputs.z1 << '\n'
              << "Final z2 = " << finalOutputs.z2 << '\n'
              << "Final theta_hat_1 = " << state[2] << '\n'
              << "Final theta_hat_2 = " << state[3] << '\n';
    return 0;
}

The Java implementation requires Java 17 or newer because it uses a record for immutable controller outputs. No external control library is required; the RK4 integrator is implemented directly.

Chapter17_Lesson2.java

// Chapter17_Lesson2.java
// Adaptive backstepping with parameter-dependent virtual control.

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 Chapter17_Lesson2 {
    private static final double THETA = 1.2;
    private static final double K1 = 2.0;
    private static final double K2 = 3.0;
    private static final double GAMMA1 = 4.0;
    private static final double GAMMA2 = 2.0;
    private static final double DT = 0.001;
    private static final double T_FINAL = 20.0;

    private Chapter17_Lesson2() {
    }

    private record Evaluation(double[] derivative, double r, double z1,
                              double z2, double alpha1, double u) {
    }

    private static double[] reference(double t) {
        double amplitude = 0.5;
        double omega = 0.7;
        return new double[] {
            amplitude * Math.sin(omega * t),
            amplitude * omega * Math.cos(omega * t),
            -amplitude * omega * omega * Math.sin(omega * t)
        };
    }

    private static Evaluation evaluate(double t, double[] state) {
        double[] ref = reference(t);
        double r = ref[0];
        double rDot = ref[1];
        double rDDot = ref[2];

        double x1 = state[0];
        double x2 = state[1];
        double thetaHat1 = state[2];
        double thetaHat2 = state[3];

        double z1 = x1 - r;
        double thetaHat1Dot = GAMMA1 * x1 * z1;
        double alpha1 = rDot - K1 * z1 - thetaHat1 * x1;
        double z2 = x2 - alpha1;
        double thetaHat2Dot = GAMMA2 * x1 * (K1 + thetaHat1) * z2;

        double alpha1DotHat = rDDot + K1 * rDot
            - (K1 + thetaHat1) * x2
            - (K1 + thetaHat1) * thetaHat2 * x1
            - thetaHat1Dot * x1;

        double u = alpha1DotHat - z1 - K2 * z2;

        double[] derivative = {
            x2 + THETA * x1,
            u,
            thetaHat1Dot,
            thetaHat2Dot
        };
        return new Evaluation(derivative, r, z1, z2, alpha1, u);
    }

    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 = evaluate(t, state).derivative();
        double[] k2 = evaluate(t + 0.5 * dt,
            addScaled(state, k1, 0.5 * dt)).derivative();
        double[] k3 = evaluate(t + 0.5 * dt,
            addScaled(state, k2, 0.5 * dt)).derivative();
        double[] k4 = evaluate(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) throws IOException {
        Locale.setDefault(Locale.US);
        Path output = Path.of("Chapter17_Lesson2_java_results.csv");
        double[] state = {1.2, -0.2, 0.0, 0.0};
        int steps = (int) Math.round(T_FINAL / DT);

        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write("time,x1,r,x2,z1,z2,theta_hat_1,theta_hat_2,alpha1,u,V");
            writer.newLine();

            for (int index = 0; index <= steps; index++) {
                double t = index * DT;
                Evaluation evaluation = evaluate(t, state);
                double tilde1 = THETA - state[2];
                double tilde2 = THETA - state[3];
                double lyapunov = 0.5 * evaluation.z1() * evaluation.z1()
                    + 0.5 * evaluation.z2() * evaluation.z2()
                    + 0.5 * tilde1 * tilde1 / GAMMA1
                    + 0.5 * tilde2 * tilde2 / GAMMA2;

                writer.write(String.format(Locale.US,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                    t, state[0], evaluation.r(), state[1], evaluation.z1(),
                    evaluation.z2(), state[2], state[3], evaluation.alpha1(),
                    evaluation.u(), lyapunov));
                writer.newLine();

                if (index < steps) {
                    state = rk4Step(t, state, DT);
                }
            }
        }

        Evaluation finalEvaluation = evaluate(T_FINAL, state);
        System.out.printf(Locale.US, "Final z1 = %.6e%n", finalEvaluation.z1());
        System.out.printf(Locale.US, "Final z2 = %.6e%n", finalEvaluation.z2());
        System.out.printf(Locale.US, "Final theta_hat_1 = %.6f%n", state[2]);
        System.out.printf(Locale.US, "Final theta_hat_2 = %.6f%n", state[3]);
        System.out.println("Results written to " + output.toAbsolutePath());
    }
}

MATLAB uses ode45 with tight tolerances and reconstructs all controller signals for plotting. In Simulink, use four Integrator blocks for \( x_1,x_2,\hat\theta_1,\hat\theta_2 \), a Clock block for \( t \), and a MATLAB Function block containing the equations in controllerSignals. The comments at the end of the file give the block-level map.

Chapter17_Lesson2.m

% Chapter17_Lesson2.m
% Adaptive backstepping with a parameter-dependent virtual control.
% The same equations can be implemented in Simulink with four Integrator
% blocks and one MATLAB Function block containing controllerSignals().

clear; close all; clc;

p.theta = 1.2;
p.k1 = 2.0;
p.k2 = 3.0;
p.gamma1 = 4.0;
p.gamma2 = 2.0;

x0 = [1.2; -0.2; 0.0; 0.0];
tspan = [0.0, 20.0];
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-10);

[t, state] = ode45(@(t, x) closedLoopDynamics(t, x, p), ...
                   tspan, x0, options);

n = numel(t);
r = zeros(n, 1);
z1 = zeros(n, 1);
z2 = zeros(n, 1);
alpha1 = zeros(n, 1);
u = zeros(n, 1);
V = zeros(n, 1);

for i = 1:n
    signals = controllerSignals(t(i), state(i, :).', p);
    r(i) = signals.r;
    z1(i) = signals.z1;
    z2(i) = signals.z2;
    alpha1(i) = signals.alpha1;
    u(i) = signals.u;

    tilde1 = p.theta - state(i, 3);
    tilde2 = p.theta - state(i, 4);
    V(i) = 0.5 * z1(i)^2 + 0.5 * z2(i)^2 ...
         + 0.5 * tilde1^2 / p.gamma1 ...
         + 0.5 * tilde2^2 / p.gamma2;
end

results = table(t, state(:,1), r, state(:,2), z1, z2, state(:,3), ...
    state(:,4), alpha1, u, V, 'VariableNames', ...
    {'time','x1','reference','x2','z1','z2','theta_hat_1', ...
     'theta_hat_2','alpha1','u','V'});
writetable(results, 'Chapter17_Lesson2_matlab_results.csv');

fprintf('Final z1 = %.6e\n', z1(end));
fprintf('Final z2 = %.6e\n', z2(end));
fprintf('Final theta_hat_1 = %.6f\n', state(end,3));
fprintf('Final theta_hat_2 = %.6f\n', state(end,4));

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

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

figure;
plot(t, state(:,3), 'LineWidth', 1.4); hold on;
plot(t, state(:,4), 'LineWidth', 1.4);
grid on; xlabel('Time [s]'); ylabel('Estimate');
legend('theta hat 1', 'theta hat 2', 'Location', 'best');
title('Overparameterized estimates');

figure;
plot(t, u, 'LineWidth', 1.4);
grid on; xlabel('Time [s]'); ylabel('u');
title('Control input');

% Simulink implementation map:
% 1. Integrate x1_dot, x2_dot, theta_hat_1_dot, theta_hat_2_dot.
% 2. Feed the four states and Clock time into a MATLAB Function block.
% 3. Place the equations in controllerSignals() inside that block.
% 4. Return u, theta_hat_1_dot, theta_hat_2_dot, z1, and z2.
% 5. Form x1_dot = x2 + theta*x1 and x2_dot = u with Sum/Gain blocks.

function dx = closedLoopDynamics(t, x, p)
    s = controllerSignals(t, x, p);
    dx = [x(2) + p.theta * x(1); ...
          s.u; ...
          s.theta_hat_1_dot; ...
          s.theta_hat_2_dot];
end

function s = controllerSignals(t, x, p)
    x1 = x(1);
    x2 = x(2);
    thetaHat1 = x(3);
    thetaHat2 = x(4);

    amplitude = 0.5;
    omega = 0.7;
    r = amplitude * sin(omega * t);
    rDot = amplitude * omega * cos(omega * t);
    rDDot = -amplitude * omega^2 * sin(omega * t);

    z1 = x1 - r;
    thetaHat1Dot = p.gamma1 * x1 * z1;
    alpha1 = rDot - p.k1 * z1 - thetaHat1 * x1;
    z2 = x2 - alpha1;
    thetaHat2Dot = p.gamma2 * x1 * (p.k1 + thetaHat1) * z2;

    alpha1DotHat = rDDot + p.k1 * rDot ...
        - (p.k1 + thetaHat1) * x2 ...
        - (p.k1 + thetaHat1) * thetaHat2 * x1 ...
        - thetaHat1Dot * x1;

    u = alpha1DotHat - z1 - p.k2 * z2;

    s = struct('r', r, 'z1', z1, 'z2', z2, 'alpha1', alpha1, ...
        'u', u, 'theta_hat_1_dot', thetaHat1Dot, ...
        'theta_hat_2_dot', thetaHat2Dot);
end

The notebook source uses NDSolveValue for the nonlinear adaptive differential equations, constructs tracking and parameter plots, and exports sampled closed-loop data to CSV.

Chapter17_Lesson2.nb

Notebook[{
 Cell["Chapter 17, Lesson 2: Virtual Controls and Parameter Dependence", "Title"],
 Cell["Adaptive backstepping simulation with an overparameterized pair of estimates.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];

theta = 1.2;
k1 = 2.0;
k2 = 3.0;
gamma1 = 4.0;
gamma2 = 2.0;
tFinal = 20.0;

r[t_] := 0.5 Sin[0.7 t];
rDot[t_] := 0.35 Cos[0.7 t];
rDDot[t_] := -0.245 Sin[0.7 t];

z1[t_] := x1[t] - r[t];
alpha1[t_] := rDot[t] - k1 z1[t] - thetaHat1[t] x1[t];
z2[t_] := x2[t] - alpha1[t];

thetaHat1Dot[t_] := gamma1 x1[t] z1[t];
thetaHat2Dot[t_] := gamma2 x1[t] (k1 + thetaHat1[t]) z2[t];

alpha1DotHat[t_] := rDDot[t] + k1 rDot[t]
  - (k1 + thetaHat1[t]) x2[t]
  - (k1 + thetaHat1[t]) thetaHat2[t] x1[t]
  - thetaHat1Dot[t] x1[t];

u[t_] := alpha1DotHat[t] - z1[t] - k2 z2[t];

eqns = {
  x1'[t] == x2[t] + theta x1[t],
  x2'[t] == u[t],
  thetaHat1'[t] == thetaHat1Dot[t],
  thetaHat2'[t] == thetaHat2Dot[t],
  x1[0] == 1.2,
  x2[0] == -0.2,
  thetaHat1[0] == 0.0,
  thetaHat2[0] == 0.0
};

solution = NDSolveValue[
  eqns,
  {x1, x2, thetaHat1, thetaHat2},
  {t, 0, tFinal},
  Method -> {\"TimeIntegration\" -> {\"ExplicitRungeKutta\", \"DifferenceOrder\" -> 4}},
  AccuracyGoal -> 10,
  PrecisionGoal -> 10
];

{x1s, x2s, thetaHat1s, thetaHat2s} = solution;

trackingPlot = Plot[
  Evaluate[{x1s[t], r[t]}],
  {t, 0, tFinal},
  PlotLegends -> {\"x1\", \"reference\"},
  PlotLabel -> \"Adaptive backstepping tracking\",
  AxesLabel -> {\"t\", \"output\"},
  GridLines -> Automatic
];

errorPlot = Plot[
  Evaluate[{
    x1s[t] - r[t],
    x2s[t] - (rDot[t] - k1 (x1s[t] - r[t]) - thetaHat1s[t] x1s[t])
  }],
  {t, 0, tFinal},
  PlotLegends -> {\"z1\", \"z2\"},
  PlotLabel -> \"Recursive error coordinates\",
  AxesLabel -> {\"t\", \"error\"},
  GridLines -> Automatic
];

estimatePlot = Plot[
  Evaluate[{thetaHat1s[t], thetaHat2s[t]}],
  {t, 0, tFinal},
  PlotLegends -> {\"thetaHat1\", \"thetaHat2\"},
  PlotLabel -> \"Overparameterized estimates\",
  AxesLabel -> {\"t\", \"estimate\"},
  GridLines -> Automatic
];

sampleTimes = Subdivide[0, tFinal, 4000];
data = Table[
  Module[{x1v, x2v, h1v, h2v, rv, z1v, a1v, z2v, h1dv, a1dh, uv, vv},
    x1v = x1s[tau]; x2v = x2s[tau];
    h1v = thetaHat1s[tau]; h2v = thetaHat2s[tau];
    rv = r[tau]; z1v = x1v - rv;
    h1dv = gamma1 x1v z1v;
    a1v = rDot[tau] - k1 z1v - h1v x1v;
    z2v = x2v - a1v;
    a1dh = rDDot[tau] + k1 rDot[tau]
      - (k1 + h1v) x2v - (k1 + h1v) h2v x1v - h1dv x1v;
    uv = a1dh - z1v - k2 z2v;
    vv = 0.5 z1v^2 + 0.5 z2v^2
      + 0.5 (theta - h1v)^2/gamma1
      + 0.5 (theta - h2v)^2/gamma2;
    {tau, x1v, rv, x2v, z1v, z2v, h1v, h2v, a1v, uv, vv}
  ],
  {tau, sampleTimes}
];

Export[
  \"Chapter17_Lesson2_mathematica_results.csv\",
  Prepend[data, {\"time\", \"x1\", \"reference\", \"x2\", \"z1\", \"z2\",
    \"theta_hat_1\", \"theta_hat_2\", \"alpha1\", \"u\", \"V\"}]
];

Print[\"Final z1 = \", x1s[tFinal] - r[tFinal]];
Print[\"Final thetaHat1 = \", thetaHat1s[tFinal]];
Print[\"Final thetaHat2 = \", thetaHat2s[tFinal]];
Column[{trackingPlot, errorPlot, estimatePlot}]"], "Input"]
}, WindowTitle -> "Chapter17_Lesson2"]

With \( \theta=1.2 \), \( k_1=2 \), \( k_2=3 \), \( \gamma_1=4 \), and \( \gamma_2=2 \), the simulations should exhibit the following behavior:

• \( z_1 \) and \( z_2 \) decay toward zero, so \( x_1 \) tracks the sinusoidal reference.
• The two estimates generally follow different transients because they cancel different terms in the recursive Lyapunov derivative.
• Both estimates may approach the true value in this particular simulation, but the stability proof requires only boundedness, not exact identification.
• Larger adaptation gains accelerate estimate motion but can increase high-frequency control activity and numerical stiffness.

The computed Lyapunov function should be nonincreasing up to numerical integration error. A coarse integration step can violate this diagnostic even when the continuous-time proof is correct; reduce the step size or tighten solver tolerances before changing the controller.

Problem 1 — Derive the first error dynamics. For \( \dot x_1=x_2+\theta x_1 \), \( z_1=x_1-r \), and \( \alpha_1=\dot r-k_1z_1-\hat\theta_1x_1 \), derive \( \dot z_1 \) in terms of \( z_1,z_2,\tilde\theta_1 \).

Solution:

\[ \begin{aligned} \dot z_1&=\dot x_1-\dot r=x_2+\theta x_1-\dot r\\ &=z_2+\alpha_1+\theta x_1-\dot r\\ &=z_2-k_1z_1+(\theta-\hat\theta_1)x_1\\ &=-k_1z_1+z_2+\tilde\theta_1x_1. \end{aligned} \]

Problem 2 — Identify every chain-rule contribution. Starting from \( \alpha_1=\dot r+k_1r-(k_1+\hat\theta_1)x_1 \), find \( \dot\alpha_1 \) before substituting the plant.

Solution:

\[ \dot\alpha_1=\ddot r+k_1\dot r -(k_1+\hat\theta_1)\dot x_1 -\dot{\hat\theta}_1x_1. \]

The last term appears because the estimate is a dynamic state. Omitting it breaks the exact derivative and invalidates the subsequent Lyapunov cancellation.

Problem 3 — Prove the derivative mismatch. Replace \( \theta \) by \( \hat\theta_2 \) in the unimplementable part of \( \dot\alpha_1 \). Show that the exact-minus-estimated derivative equals \( -(k_1+\hat\theta_1)\tilde\theta_2x_1 \).

Solution:

\[ \begin{aligned} \dot\alpha_1-\widehat{\dot\alpha}_1 &=-(k_1+\hat\theta_1)\theta x_1 +(k_1+\hat\theta_1)\hat\theta_2x_1\\ &=-(k_1+\hat\theta_1)(\theta-\hat\theta_2)x_1\\ &=-(k_1+\hat\theta_1)\tilde\theta_2x_1. \end{aligned} \]

Problem 4 — Recover the second update law. Suppose the derivative of the augmented Lyapunov function contains

\[ \tilde\theta_2x_1(k_1+\hat\theta_1)z_2 -\frac{1}{\gamma_2}\tilde\theta_2\dot{\hat\theta}_2. \]

Find an update law that cancels this expression.

Solution:

\[ \dot{\hat\theta}_2= \gamma_2x_1(k_1+\hat\theta_1)z_2. \]

Substitution makes the two terms identical with opposite signs for every value of \( \tilde\theta_2 \).

Problem 5 — Why does tracking not prove parameter identification? Explain why \( z_1,z_2→0 \) does not necessarily imply \( \hat\theta_1,\hat\theta_2→\theta \).

Solution:

The adaptation laws stop changing when their regressors vanish. If the reference and state signals do not excite the uncertain dynamics over a sufficiently rich interval, multiple parameter values can produce the same asymptotic tracking behavior. The Lyapunov function establishes bounded estimation errors, while exact convergence requires an additional excitation condition.

Problem 6 — Constant-reference specialization. Let \( r(t)=r_0 \). Simplify the virtual control, its estimated derivative, and the control law.

Solution: Since \( \dot r=\ddot r=0 \),

\[ \begin{aligned} \alpha_1&=-k_1(x_1-r_0)-\hat\theta_1x_1,\\ \widehat{\dot\alpha}_1 &=-(k_1+\hat\theta_1)x_2 -(k_1+\hat\theta_1)\hat\theta_2x_1 -\dot{\hat\theta}_1x_1,\\ u&=\widehat{\dot\alpha}_1-z_1-k_2z_2. \end{aligned} \]

The Lyapunov proof is unchanged because the reference derivatives cancel algebraically and do not appear in \( \dot V_2 \).

A virtual control in adaptive backstepping generally depends on online parameter estimates. Differentiating it introduces estimate derivatives and can reintroduce unknown plant parameters through the state dynamics. In the worked two-state design, a second estimate converts the unknown derivative contribution into a cancelable parameter-error term. The resulting control and two update laws produce \( \dot V_2=-k_1z_1^2-k_2z_2^2 \), proving boundedness and asymptotic tracking. The price is overparameterization, which motivates the recursive tuning-function construction developed next.