After completing this lesson, students should be able to:

• recognize a second-order parameter strict-feedback model;
• construct the first backstepping error and a stabilizing virtual control;
• augment the Lyapunov function with a parameter-estimation error term;
• derive the actual control law and adaptive update law by cancellation;
• prove global boundedness and asymptotic tracking;
• distinguish tracking convergence from parameter convergence; and
• implement the same continuous-time design with several numerical tools.

flowchart TD
  A["Specify strict-feedback plant and reference"] --> B["Define z1 = x1 - r"]
  B --> C["Choose virtual control alpha"]
  C --> D["Define z2 = x2 - alpha"]
  D --> E["Augment Lyapunov function with parameter error"]
  E --> F["Choose actual control u"]
  F --> G["Choose parameter update law"]
  G --> H["Obtain negative-semidefinite Vdot"]
  H --> I["Use boundedness and Barbalat to prove tracking"]
  I --> J["Simulate and inspect excitation"]
        

Consider the second-order nonlinear plant

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

with the known nonlinear regressor

\[ \phi(x_1,x_2)=x_1^3+\frac{1}{2}\sin(x_2), \]

where \(x_1\) and \(x_2\) are measured, \(u\) is the control input, and \(\theta\in\mathbb{R}\) is an unknown constant. The uncertainty is linearly parameterized even though the state dependence is nonlinear.

The desired output is \(y=x_1\), and the objective is

\[ \lim_{t→\infty}\bigl(x_1(t)-r(t)\bigr)=0, \]

for a reference \(r(t)\) whose first two derivatives are available and bounded. For simulation, we use

\[ r(t)=0.8\sin(0.6t),\quad \dot r(t)=0.48\cos(0.6t),\quad \ddot r(t)=-0.288\sin(0.6t). \]

Assumptions

1. \(\theta\) is constant but unknown.
2. \(x_1,x_2\) are available for feedback.
3. \(r,\dot r,\ddot r\) are bounded and known.
4. The actuator can realize the derived continuous control law.
5. The design gains satisfy \(k_1,k_2,\gamma > 0\).

Define the first tracking error

\[ z_1=x_1-r. \]

Its dynamics are

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

At this stage, \(x_2\) is treated as a virtual control. Choose the stabilizing function

\[ \alpha=\dot r-k_1z_1. \]

Introduce the second backstepping error

\[ z_2=x_2-\alpha. \]

Since \(x_2=z_2+\alpha\),

\[ \dot z_1=-k_1z_1+z_2. \]

The term \(z_2\) measures the failure of the physical state \(x_2\) to equal the desired virtual control \(\alpha\).

First-step Lyapunov calculation

For \(V_1=\frac12 z_1^2\),

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

The negative term is desirable; the cross term \(z_1z_2\) must be cancelled at the next step.

Differentiate the virtual control:

\[ \dot\alpha=\ddot r-k_1\dot z_1 =\ddot r-k_1(x_2-\dot r). \]

The second error satisfies

\[ \dot z_2=\dot x_2-\dot\alpha =u+\theta\phi(x_1,x_2)-\dot\alpha. \]

Let \(\hat\theta\) be the online estimate and define the parameter error using the convention

\[ \tilde\theta=\theta-\hat\theta. \]

Choose the actual controller

\[ \boxed{ u=\dot\alpha-z_1-k_2z_2-\hat\theta\phi(x_1,x_2) } \]

which yields

\[ \dot z_2=-z_1-k_2z_2+\tilde\theta\phi(x_1,x_2). \]

Composite Lyapunov function

Augment the state-error energy with parameter-error energy:

\[ V=\frac12z_1^2+\frac12z_2^2+\frac{1}{2\gamma}\tilde\theta^2. \]

Because \(\theta\) is constant,

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

Differentiation gives

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

Choose the adaptation law

\[ \boxed{ \dot{\hat\theta}=\gamma z_2\phi(x_1,x_2) } \]

to cancel the indefinite parameter-error term. Therefore,

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

flowchart TD
  R["Reference r, r_dot, r_ddot"] --> E1["z1 = x1 - r"]
  X1["Measured x1"] --> E1
  E1 --> A["alpha = r_dot - k1 z1"]
  X2["Measured x2"] --> E2["z2 = x2 - alpha"]
  A --> E2
  X1 --> P["phi = x1^3 + 0.5 sin(x2)"]
  X2 --> P
  E2 --> AD["thetaHat_dot = gamma z2 phi"]
  P --> AD
  AD --> EST["Integrator gives thetaHat"]
  E1 --> U["u = alpha_dot - z1 - k2 z2 - thetaHat phi"]
  E2 --> U
  P --> U
  EST --> U
  U --> PLANT["Nonlinear plant"]
  PLANT --> X1
  PLANT --> X2
        

The estimator is driven by the product of the final backstepping error and the nonlinear regressor. The estimator is not an independent identifier; it is selected to make the composite Lyapunov derivative nonpositive.

Theorem

Consider the plant in Section 2 under the stated assumptions. Apply the virtual control, actual control, and update law

\[ \begin{aligned} \alpha&=\dot r-k_1z_1,\\ u&=\dot\alpha-z_1-k_2z_2-\hat\theta\phi,\\ \dot{\hat\theta}&=\gamma z_2\phi, \end{aligned} \qquad k_1,k_2,\gamma >0. \]

Then all closed-loop signals are bounded and the tracking errors satisfy \(z_1(t)→ 0\) and \(z_2(t)→ 0\) as \(t→\infty\).

Proof

Step 1: Positive definiteness. The function

\[ V=\frac12z_1^2+\frac12z_2^2+\frac{1}{2\gamma}\tilde\theta^2 \]

is positive definite and radially unbounded in \((z_1,z_2,\tilde\theta)\) because \(\gamma >0\).

Step 2: Monotonicity. The controller and adaptive law give

\[ \dot V=-k_1z_1^2-k_2z_2^2\leq0. \]

Thus \(V(t)\leq V(0)\) for every \(t\geq0\).

Step 3: Boundedness. Since each term in \(V\) is nonnegative,

\[ |z_1(t)|\leq\sqrt{2V(0)},\quad |z_2(t)|\leq\sqrt{2V(0)},\quad |\tilde\theta(t)|\leq\sqrt{2\gamma V(0)}. \]

Hence \(z_1,z_2,\hat\theta\) are bounded. Because the reference and its derivatives are bounded, \(x_1=z_1+r\), \(\alpha=\dot r-k_1z_1\), and \(x_2=z_2+\alpha\) are bounded. Smoothness of \(\phi\) then implies bounded \(\phi,\dot{\hat\theta},u\).

Step 4: Square integrability. Integrating the Lyapunov derivative,

\[ \int_0^T\left(k_1z_1^2+k_2z_2^2\right)d\tau =V(0)-V(T)\leq V(0). \]

Letting \(T→\infty\) gives

\[ z_1,z_2\in L_2\cap L_\infty. \]

Step 5: Uniform continuity. The closed-loop error equations

\[ \dot z_1=-k_1z_1+z_2,\qquad \dot z_2=-z_1-k_2z_2+\tilde\theta\phi \]

have bounded right-hand sides, so \(\dot z_1\) and \(\dot z_2\) are bounded. Therefore \(z_1\) and \(z_2\) are uniformly continuous.

Step 6: Barbalat's lemma. Uniform continuity together with square integrability implies

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

Since \(z_1=x_1-r\), asymptotic output tracking follows. This completes the proof. \(\square\)

Guaranteed conclusions

• global boundedness of the error coordinates and parameter estimate;
• finite integral of the squared tracking errors;
• asymptotic convergence of \(z_1\) and \(z_2\);
• asymptotic tracking of \(x_1\) to \(r\).

Not automatically guaranteed

The proof does not by itself imply \(\hat\theta(t)→\theta\). Parameter convergence requires sufficient excitation. In this scalar example, a useful condition is that the regressor be persistently exciting: there must exist constants \(T_0>0\) and \(\mu>0\) such that

\[ \int_t^{t+T_0}\phi^2\bigl(x_1(\tau),x_2(\tau)\bigr)d\tau\geq\mu \quad\text{for all }t\geq0. \]

If the reference tends to a constant and the states settle near an equilibrium where \(\phi=0\), tracking may be excellent while the estimate converges to a value different from the true parameter.

The implementations use the following numerical values:

\[ \theta=1.4,\quad k_1=2,\quad k_2=3,\quad \gamma=4, \]

\[ x_1(0)=1.2,\quad x_2(0)=-0.4,\quad \hat\theta(0)=0. \]

The Python, C++, and Java programs use a fixed-step fourth-order Runge–Kutta method with \(h=0.002\) seconds. MATLAB uses ode45, and Mathematica uses an explicit Runge–Kutta time integrator. Each implementation evaluates the controller at every internal integration stage, which is essential for a nonlinear closed loop.

Signals to inspect

1. Tracking: compare \(x_1(t)\) with \(r(t)\).
2. Error coordinates: inspect \(z_1,z_2\).
3. Estimate: compare \(\hat\theta\) with the true parameter.
4. Control effort: check peak and steady-state values of \(u\).
5. Excitation: inspect whether \(\phi\) remains informative.

Code file: Chapter17_Lesson5.py

"""Chapter 17, Lesson 5: Adaptive backstepping for a second-order nonlinear system."""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Parameters:
    theta_true: float = 1.4
    k1: float = 2.0
    k2: float = 3.0
    gamma: float = 4.0
    amplitude: float = 0.8
    omega: float = 0.6


def reference(t: float, p: Parameters) -> tuple[float, float, float]:
    r = p.amplitude * np.sin(p.omega * t)
    r_dot = p.amplitude * p.omega * np.cos(p.omega * t)
    r_ddot = -p.amplitude * p.omega**2 * np.sin(p.omega * t)
    return float(r), float(r_dot), float(r_ddot)


def regressor(x1: float, x2: float) -> float:
    return x1**3 + 0.5 * np.sin(x2)


def closed_loop_rhs(t: float, state: np.ndarray, p: Parameters) -> np.ndarray:
    x1, x2, theta_hat = state
    r, r_dot, r_ddot = reference(t, p)

    z1 = x1 - r
    alpha = r_dot - p.k1 * z1
    z2 = x2 - alpha
    alpha_dot = r_ddot - p.k1 * (x2 - r_dot)
    phi = regressor(x1, x2)

    u = alpha_dot - z1 - p.k2 * z2 - theta_hat * phi
    theta_hat_dot = p.gamma * z2 * phi

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


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


def simulate(
    p: Parameters,
    t_final: float = 30.0,
    dt: float = 0.002,
    x0: np.ndarray | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    if dt <= 0.0 or t_final <= 0.0:
        raise ValueError("dt and t_final must be positive")

    state = np.array([1.2, -0.4, 0.0], dtype=float) if x0 is None else np.array(x0, dtype=float)
    if state.shape != (3,):
        raise ValueError("x0 must contain [x1, x2, theta_hat]")

    steps = int(round(t_final / dt))
    time = np.linspace(0.0, steps * dt, steps + 1)
    states = np.empty((steps + 1, 3), dtype=float)
    controls = np.empty(steps + 1, dtype=float)
    states[0] = state

    for i, t in enumerate(time):
        x1, x2, theta_hat = state
        r, r_dot, r_ddot = reference(float(t), p)
        z1 = x1 - r
        alpha = r_dot - p.k1 * z1
        z2 = x2 - alpha
        alpha_dot = r_ddot - p.k1 * (x2 - r_dot)
        controls[i] = alpha_dot - z1 - p.k2 * z2 - theta_hat * regressor(x1, x2)

        if i < steps:
            state = rk4_step(closed_loop_rhs, float(t), state, dt, p)
            if not np.all(np.isfinite(state)):
                raise FloatingPointError(f"non-finite state at t={t:.6f}")
            states[i + 1] = state

    return time, states, controls


def save_results(time: np.ndarray, states: np.ndarray, controls: np.ndarray, p: Parameters) -> None:
    output_dir = Path(__file__).resolve().parent
    csv_path = output_dir / "Chapter17_Lesson5_results.csv"

    with csv_path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.writer(handle)
        writer.writerow(["t", "x1", "x2", "theta_hat", "reference", "tracking_error", "u"])
        for t, state, u in zip(time, states, controls):
            r, _, _ = reference(float(t), p)
            writer.writerow([t, state[0], state[1], state[2], r, state[0] - r, u])

    reference_values = np.array([reference(float(t), p)[0] for t in time])

    plt.figure(figsize=(9, 5))
    plt.plot(time, states[:, 0], label="x1")
    plt.plot(time, reference_values, "--", label="r")
    plt.xlabel("Time (s)")
    plt.ylabel("Position / reference")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter17_Lesson5_tracking.png", dpi=180)

    plt.figure(figsize=(9, 5))
    plt.plot(time, states[:, 2], label="theta_hat")
    plt.axhline(p.theta_true, linestyle="--", label="theta true")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter17_Lesson5_parameter.png", dpi=180)

    plt.figure(figsize=(9, 5))
    plt.plot(time, controls, label="u")
    plt.xlabel("Time (s)")
    plt.ylabel("Control input")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter17_Lesson5_control.png", dpi=180)

    plt.show()


if __name__ == "__main__":
    parameters = Parameters()
    t_grid, state_history, control_history = simulate(parameters)
    save_results(t_grid, state_history, control_history, parameters)

    final_reference, _, _ = reference(float(t_grid[-1]), parameters)
    print(f"Final tracking error: {state_history[-1, 0] - final_reference:.6e}")
    print(f"Final parameter estimate: {state_history[-1, 2]:.6f}")

Required libraries are numpy and matplotlib. The script writes a CSV file and saves tracking, parameter, and control plots.

Code file: Chapter17_Lesson5.cpp

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

struct Parameters {
    double thetaTrue = 1.4;
    double k1 = 2.0;
    double k2 = 3.0;
    double gamma = 4.0;
    double amplitude = 0.8;
    double omega = 0.6;
};

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

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

Reference referenceSignal(double t, const Parameters& p) {
    return {
        p.amplitude * std::sin(p.omega * t),
        p.amplitude * p.omega * std::cos(p.omega * t),
        -p.amplitude * p.omega * p.omega * std::sin(p.omega * t)
    };
}

double regressor(double x1, double x2) {
    return x1 * x1 * x1 + 0.5 * std::sin(x2);
}

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 rhs(double t, const State& x, const Parameters& p) {
    const double x1 = x[0];
    const double x2 = x[1];
    const double thetaHat = x[2];
    const Reference ref = referenceSignal(t, p);

    const double z1 = x1 - ref.r;
    const double alpha = ref.rDot - p.k1 * z1;
    const double z2 = x2 - alpha;
    const double alphaDot = ref.rDDot - p.k1 * (x2 - ref.rDot);
    const double phi = regressor(x1, x2);

    const double u = alphaDot - z1 - p.k2 * z2 - thetaHat * phi;
    const double thetaHatDot = p.gamma * z2 * phi;

    return {x2, u + p.thetaTrue * phi, thetaHatDot};
}

State rk4Step(double t, const State& x, double h, const Parameters& p) {
    const State k1 = rhs(t, x, p);
    const State k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h), p);
    const State k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h), p);
    const State k4 = rhs(t + h, addScaled(x, k3, h), p);

    return {
        x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
        x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]),
        x[2] + (h / 6.0) * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2])
    };
}

double controlInput(double t, const State& x, const Parameters& p) {
    const Reference ref = referenceSignal(t, p);
    const double z1 = x[0] - ref.r;
    const double alpha = ref.rDot - p.k1 * z1;
    const double z2 = x[1] - alpha;
    const double alphaDot = ref.rDDot - p.k1 * (x[1] - ref.rDot);
    return alphaDot - z1 - p.k2 * z2 - x[2] * regressor(x[0], x[1]);
}

int main() {
    try {
        const Parameters p;
        const double dt = 0.002;
        const double tFinal = 30.0;
        const int steps = static_cast<int>(std::llround(tFinal / dt));
        State state{1.2, -0.4, 0.0};

        std::ofstream output("Chapter17_Lesson5_results_cpp.csv");
        if (!output) {
            throw std::runtime_error("Unable to open output CSV file");
        }
        output << "t,x1,x2,theta_hat,reference,tracking_error,u\n";
        output << std::setprecision(15);

        for (int i = 0; i <= steps; ++i) {
            const double t = i * dt;
            const Reference ref = referenceSignal(t, p);
            const double u = controlInput(t, state, p);
            output << t << ',' << state[0] << ',' << state[1] << ',' << state[2] << ','
                   << ref.r << ',' << state[0] - ref.r << ',' << u << '\n';

            if (i < steps) {
                state = rk4Step(t, state, dt, p);
                if (!std::isfinite(state[0]) || !std::isfinite(state[1]) || !std::isfinite(state[2])) {
                    throw std::runtime_error("Simulation produced a non-finite state");
                }
            }
        }

        const Reference finalRef = referenceSignal(tFinal, p);
        std::cout << "Final tracking error: " << state[0] - finalRef.r << '\n';
        std::cout << "Final parameter estimate: " << state[2] << '\n';
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

Compile with a C++17 compiler, for example:

g++ -std=c++17 -O2 Chapter17_Lesson5.cpp -o Chapter17_Lesson5
./Chapter17_Lesson5

Code file: Chapter17_Lesson5.java

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_Lesson5 {
    private static final class Parameters {
        final double thetaTrue = 1.4;
        final double k1 = 2.0;
        final double k2 = 3.0;
        final double gamma = 4.0;
        final double amplitude = 0.8;
        final double omega = 0.6;
    }

    private record Reference(double r, double rDot, double rDDot) {}

    private static Reference referenceSignal(double t, Parameters p) {
        return new Reference(
            p.amplitude * Math.sin(p.omega * t),
            p.amplitude * p.omega * Math.cos(p.omega * t),
            -p.amplitude * p.omega * p.omega * Math.sin(p.omega * t)
        );
    }

    private static double regressor(double x1, double x2) {
        return x1 * x1 * x1 + 0.5 * Math.sin(x2);
    }

    private static double[] rhs(double t, double[] state, Parameters p) {
        double x1 = state[0];
        double x2 = state[1];
        double thetaHat = state[2];
        Reference ref = referenceSignal(t, p);

        double z1 = x1 - ref.r();
        double alpha = ref.rDot() - p.k1 * z1;
        double z2 = x2 - alpha;
        double alphaDot = ref.rDDot() - p.k1 * (x2 - ref.rDot());
        double phi = regressor(x1, x2);

        double u = alphaDot - z1 - p.k2 * z2 - thetaHat * phi;
        double thetaHatDot = p.gamma * z2 * phi;

        return new double[] {x2, u + p.thetaTrue * phi, thetaHatDot};
    }

    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[] state, double h, Parameters p) {
        double[] k1 = rhs(t, state, p);
        double[] k2 = rhs(t + 0.5 * h, addScaled(state, k1, 0.5 * h), p);
        double[] k3 = rhs(t + 0.5 * h, addScaled(state, k2, 0.5 * h), p);
        double[] k4 = rhs(t + h, addScaled(state, k3, h), p);

        return new double[] {
            state[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
            state[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]),
            state[2] + (h / 6.0) * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2])
        };
    }

    private static double controlInput(double t, double[] state, Parameters p) {
        Reference ref = referenceSignal(t, p);
        double z1 = state[0] - ref.r();
        double alpha = ref.rDot() - p.k1 * z1;
        double z2 = state[1] - alpha;
        double alphaDot = ref.rDDot() - p.k1 * (state[1] - ref.rDot());
        return alphaDot - z1 - p.k2 * z2 - state[2] * regressor(state[0], state[1]);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Parameters p = new Parameters();
        double dt = 0.002;
        double tFinal = 30.0;
        int steps = (int) Math.round(tFinal / dt);
        double[] state = {1.2, -0.4, 0.0};
        Path outputPath = Path.of("Chapter17_Lesson5_results_java.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(outputPath, StandardCharsets.UTF_8)) {
            writer.write("t,x1,x2,theta_hat,reference,tracking_error,u\n");
            for (int i = 0; i <= steps; i++) {
                double t = i * dt;
                Reference ref = referenceSignal(t, p);
                double u = controlInput(t, state, p);
                writer.write(String.format(
                    Locale.US,
                    "%.15g,%.15g,%.15g,%.15g,%.15g,%.15g,%.15g%n",
                    t, state[0], state[1], state[2], ref.r(), state[0] - ref.r(), u
                ));

                if (i < steps) {
                    state = rk4Step(t, state, dt, p);
                    if (!Double.isFinite(state[0]) || !Double.isFinite(state[1]) || !Double.isFinite(state[2])) {
                        throw new ArithmeticException("Simulation produced a non-finite state");
                    }
                }
            }

            Reference finalRef = referenceSignal(tFinal, p);
            System.out.printf(Locale.US, "Final tracking error: %.6e%n", state[0] - finalRef.r());
            System.out.printf(Locale.US, "Final parameter estimate: %.6f%n", state[2]);
        } catch (IOException | ArithmeticException ex) {
            System.err.println("Error: " + ex.getMessage());
            System.exit(1);
        }
    }
}

Compile and execute with Java 17 or later:

javac Chapter17_Lesson5.java
java Chapter17_Lesson5

Code file: Chapter17_Lesson5.m

%% Chapter17_Lesson5.m
% Adaptive backstepping for a second-order nonlinear strict-feedback plant.
clear; clc; close all;

p.thetaTrue = 1.4;
p.k1 = 2.0;
p.k2 = 3.0;
p.gamma = 4.0;
p.amplitude = 0.8;
p.omega = 0.6;

x0 = [1.2; -0.4; 0.0]; % [x1; x2; thetaHat]
tSpan = [0, 30];
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10, 'MaxStep', 0.02);

[t, x] = ode45(@(t, x) closedLoopRhs(t, x, p), tSpan, x0, options);

r = p.amplitude * sin(p.omega * t);
rDot = p.amplitude * p.omega * cos(p.omega * t);
rDDot = -p.amplitude * p.omega^2 * sin(p.omega * t);
z1 = x(:,1) - r;
alpha = rDot - p.k1 .* z1;
z2 = x(:,2) - alpha;
phi = x(:,1).^3 + 0.5 .* sin(x(:,2));
alphaDot = rDDot - p.k1 .* (x(:,2) - rDot);
u = alphaDot - z1 - p.k2 .* z2 - x(:,3) .* phi;

results = table(t, x(:,1), x(:,2), x(:,3), r, z1, z2, u, ...
    'VariableNames', {'t','x1','x2','thetaHat','reference','z1','z2','u'});
writetable(results, 'Chapter17_Lesson5_results_matlab.csv');

figure('Name', 'Tracking');
plot(t, x(:,1), 'LineWidth', 1.4); hold on;
plot(t, r, '--', 'LineWidth', 1.4);
grid on; xlabel('Time (s)'); ylabel('x_1 and r');
legend('x_1', 'r', 'Location', 'best');

figure('Name', 'Parameter estimate');
plot(t, x(:,3), 'LineWidth', 1.4); hold on;
yline(p.thetaTrue, '--', 'LineWidth', 1.4);
grid on; xlabel('Time (s)'); ylabel('\theta estimate');
legend('\hat{\theta}', '\theta', 'Location', 'best');

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

fprintf('Final tracking error: %.6e\n', z1(end));
fprintf('Final parameter estimate: %.6f\n', x(end,3));

function dx = closedLoopRhs(t, x, p)
    x1 = x(1);
    x2 = x(2);
    thetaHat = x(3);

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

    z1 = x1 - r;
    alpha = rDot - p.k1 * z1;
    z2 = x2 - alpha;
    alphaDot = rDDot - p.k1 * (x2 - rDot);
    phi = x1^3 + 0.5 * sin(x2);

    u = alphaDot - z1 - p.k2 * z2 - thetaHat * phi;
    thetaHatDot = p.gamma * z2 * phi;

    dx = [x2;
          u + p.thetaTrue * phi;
          thetaHatDot];
end

Build the plant with two Integrator blocks. The second integrator output is \(x_2\), the first integrator output is \(x_1\), and the acceleration input is \(u+\theta\phi\). Use a third Integrator block for \(\hat\theta\). Place the following code in a MATLAB Function block.

Code file: Chapter17_Lesson5_SimulinkController.m

function [u, thetaHatDot, z1, z2] = Chapter17_Lesson5_SimulinkController( ...
    t, x1, x2, thetaHat, k1, k2, gamma, amplitude, omega)
%#codegen
% MATLAB Function block implementation for the adaptive backstepping law.

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

z1 = x1 - r;
alpha = rDot - k1 * z1;
z2 = x2 - alpha;
alphaDot = rDDot - k1 * (x2 - rDot);
phi = x1^3 + 0.5 * sin(x2);

u = alphaDot - z1 - k2 * z2 - thetaHat * phi;
thetaHatDot = gamma * z2 * phi;
end

Recommended block connections

1. Use a Clock block for \(t\).
2. Feed \(x_1,x_2,\hat\theta\) into the MATLAB Function block.
3. Integrate \(\dot{\hat\theta}\) to produce \(\hat\theta\).
4. Compute the plant regressor in a second MATLAB Function block or inside the plant subsystem.
5. Use a variable-step solver such as ode45 initially; then compare with a fixed-step solver for deployment studies.

Code file: Chapter17_Lesson5.nb

(* Chapter17_Lesson5.nb -- Wolfram Language source that may be opened as a notebook. *)
ClearAll["Global`*"];

params = {
  thetaTrue -> 1.4,
  k1 -> 2.0,
  k2 -> 3.0,
  gamma -> 4.0,
  amplitude -> 0.8,
  omega -> 0.6
};

r[t_] := amplitude Sin[omega t];
rd[t_] := amplitude omega Cos[omega t];
rdd[t_] := -amplitude omega^2 Sin[omega t];
phi[x1_, x2_] := x1^3 + 0.5 Sin[x2];

z1[t_] := x1[t] - r[t];
alpha[t_] := rd[t] - k1 z1[t];
z2[t_] := x2[t] - alpha[t];
alphaDot[t_] := rdd[t] - k1 (x2[t] - rd[t]);
control[t_] := alphaDot[t] - z1[t] - k2 z2[t] - thetaHat[t] phi[x1[t], x2[t]];

solution = NDSolveValue[
  {
    x1'[t] == x2[t],
    x2'[t] == control[t] + thetaTrue phi[x1[t], x2[t]],
    thetaHat'[t] == gamma z2[t] phi[x1[t], x2[t]],
    x1[0] == 1.2,
    x2[0] == -0.4,
    thetaHat[0] == 0.0
  } /. params,
  {x1, x2, thetaHat},
  {t, 0, 30},
  Method -> {"TimeIntegration" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}},
  MaxStepFraction -> 1/1500
];

{x1Sol, x2Sol, thetaHatSol} = solution;

trackingPlot = Plot[
  Evaluate[{x1Sol[t], r[t] /. params}],
  {t, 0, 30},
  PlotLegends -> {"x1", "r"},
  AxesLabel -> {"t", "tracking"},
  PlotRange -> All,
  GridLines -> Automatic
];

parameterPlot = Plot[
  Evaluate[{thetaHatSol[t], thetaTrue /. params}],
  {t, 0, 30},
  PlotLegends -> {"thetaHat", "theta true"},
  AxesLabel -> {"t", "parameter"},
  PlotRange -> All,
  GridLines -> Automatic
];

controlPlot = Plot[
  Evaluate[control[t] /. params /. {
      x1[t] -> x1Sol[t], x2[t] -> x2Sol[t], thetaHat[t] -> thetaHatSol[t]
    }],
  {t, 0, 30},
  PlotLegends -> {"u"},
  AxesLabel -> {"t", "control"},
  PlotRange -> All,
  GridLines -> Automatic
];

sampleTimes = Range[0, 30, 0.01];
data = Table[
  With[
    {
      x1v = x1Sol[tau],
      x2v = x2Sol[tau],
      thv = thetaHatSol[tau],
      rv = r[tau] /. params
    },
    {tau, x1v, x2v, thv, rv, x1v - rv}
  ],
  {tau, sampleTimes}
];

Export["Chapter17_Lesson5_results_mathematica.csv",
  Prepend[data, {"t", "x1", "x2", "thetaHat", "reference", "trackingError"}]
];

Print["Final tracking error: ", x1Sol[30] - (r[30] /. params)];
Print["Final parameter estimate: ", thetaHatSol[30]];
GraphicsGrid[{{trackingPlot}, {parameterPlot}, {controlPlot}}]

Effect of \(k_1\)

The first gain shapes the virtual first-order target dynamics \(\dot z_1=-k_1z_1+z_2\). A larger value accelerates the nominal decay of \(z_1\), but increases \(\alpha\) and its derivative.

Effect of \(k_2\)

The second gain damps the final backstepping error. Increasing it normally reduces \(z_2\), but may amplify control activity and measurement noise.

Effect of \(\gamma\)

The adaptation gain controls how aggressively the estimate responds to \(z_2\phi\). A large value can produce fast transient adjustment, but also magnifies numerical and measurement effects. The Lyapunov proof permits any positive value; practical tuning remains a bandwidth and noise trade-off.

Recommended diagnostic checks

• repeat the simulation with smaller integration steps;
• plot the numerical Lyapunov function and verify that it is nearly nonincreasing;
• check that controller evaluation occurs at every RK4 stage;
• avoid introducing saturation without also analyzing its effect on the adaptive loop;
• compare several references to distinguish tracking from identification behavior.

Problem 1 (Derive the virtual-control dynamics): Starting from \(z_1=x_1-r\) and \(\alpha=\dot r-k_1z_1\), prove that \(\dot z_1=-k_1z_1+z_2\) when \(z_2=x_2-\alpha\).

Solution:

\[ \dot z_1=x_2-\dot r=(z_2+\alpha)-\dot r =z_2+(\dot r-k_1z_1)-\dot r=-k_1z_1+z_2. \]

Problem 2 (Sign convention): Suppose the parameter error is instead defined as \(\bar\theta=\hat\theta-\theta\). Determine the adaptation law that produces the same Lyapunov derivative.

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

\[ \dot z_2=-z_1-k_2z_2-\bar\theta\phi. \]

With \(V=\frac12z_1^2+\frac12z_2^2+\frac{1}{2\gamma}\bar\theta^2\),

\[ \dot V=-k_1z_1^2-k_2z_2^2 +\bar\theta\left(-z_2\phi+\frac{1}{\gamma}\dot{\hat\theta}\right). \]

Thus the same update law is obtained:

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

Problem 3 (Constant reference): Let \(r(t)=r_0\). Write the controller explicitly.

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

\[ z_1=x_1-r_0,\quad \alpha=-k_1z_1,\quad z_2=x_2+k_1z_1, \]

\[ \dot\alpha=-k_1x_2, \]

\[ u=-k_1x_2-z_1-k_2z_2-\hat\theta\phi, \qquad \dot{\hat\theta}=\gamma z_2\phi. \]

Problem 4 (Vector parameter extension): Replace the scalar uncertainty by \(\boldsymbol\theta^T\boldsymbol\phi(x)\). Derive the corresponding controller and update law.

Solution: Let

\[ \dot x_2=u+\boldsymbol\theta^T\boldsymbol\phi(x),\qquad \tilde{\boldsymbol\theta}=\boldsymbol\theta-\hat{\boldsymbol\theta}. \]

Choose

\[ u=\dot\alpha-z_1-k_2z_2- \hat{\boldsymbol\theta}^{T}\boldsymbol\phi(x), \]

\[ \dot{\hat{\boldsymbol\theta}} =\boldsymbol\Gamma\boldsymbol\phi(x)z_2, \qquad \boldsymbol\Gamma=\boldsymbol\Gamma^T>0. \]

For

\[ V=\frac12z_1^2+\frac12z_2^2+ \frac12\tilde{\boldsymbol\theta}^{T} \boldsymbol\Gamma^{-1}\tilde{\boldsymbol\theta}, \]

the parameter cross term cancels and

\[ \dot V=-k_1z_1^2-k_2z_2^2. \]

Problem 5 (Lyapunov bound): Given \(V(0)=5\), \(\gamma=4\), and \(\theta=1.4\), find a guaranteed bound on \(|\hat\theta(t)|\).

Solution: Since

\[ \frac{1}{2\gamma}\tilde\theta^2\leq V(t)\leq V(0)=5, \]

\[ |\tilde\theta(t)|\leq\sqrt{2\gamma V(0)} =\sqrt{40}. \]

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

\[ |\hat\theta(t)|\leq|\theta|+|\tilde\theta(t)| \leq1.4+\sqrt{40}\approx7.7249. \]

Problem 6 (Programming exercise): Modify the Python program so that it computes

\[ V(t)=\frac12z_1^2+\frac12z_2^2+ \frac{1}{2\gamma}(\theta-\hat\theta)^2 \]

and compares the finite-difference derivative of \(V\) with \(-k_1z_1^2-k_2z_2^2\).

Solution: Add the following code after the simulation:

Code file: Chapter17_Lesson5_Ex1.py

r_values = np.array([reference(float(t), parameters)[0] for t in t_grid])
r_dot_values = np.array([reference(float(t), parameters)[1] for t in t_grid])
z1 = state_history[:, 0] - r_values
alpha = r_dot_values - parameters.k1 * z1
z2 = state_history[:, 1] - alpha
theta_tilde = parameters.theta_true - state_history[:, 2]
V = 0.5 * z1**2 + 0.5 * z2**2 + theta_tilde**2 / (2.0 * parameters.gamma)
V_dot_numeric = np.gradient(V, t_grid)
V_dot_theory = -parameters.k1 * z1**2 - parameters.k2 * z2**2
print("Maximum derivative mismatch:", np.max(np.abs(V_dot_numeric - V_dot_theory)))

A small mismatch away from the initial and final samples is expected because np.gradient approximates a derivative from discrete data.

The example demonstrated the complete adaptive-backstepping mechanism for a second-order nonlinear strict-feedback system. The first recursion generated a virtual control and a residual cross term. The second recursion selected the physical input, while the adaptive law cancelled the parameter-error product in the composite Lyapunov derivative. The resulting derivative is negative semidefinite, providing global boundedness and, through Barbalat's lemma, asymptotic tracking. Parameter convergence remains a separate excitation-dependent question.