After this lesson, students should be able to:

• model a hard actuator limit as a saturation map and a saturation-deficiency signal;
• derive a linearly parameterized adaptive position controller;
• construct an auxiliary anti-windup state and a compensated tracking surface;
• apply a projection operator without destroying the Lyapunov inequality;
• distinguish convergence of compensated error from physically feasible tracking;
• implement and test the same design in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The derivation uses only concepts introduced earlier in the course: tracking surfaces, linear parameterization, gradient adaptation, projection, Barbalat's lemma, robust modification concepts, and actuator-saturation handling.

The controller first computes an unconstrained command \( u_c \). The actuator applies the saturated signal \( u \). Their difference is injected into a stable auxiliary filter. Adaptation uses the compensated surface rather than the raw surface.

flowchart TD
  R["Reference qd, qd_dot, qd_ddot"] --> E["Tracking errors e and e_dot"]
  X["Measured q and q_dot"] --> E
  E --> S["Surface s = e_dot + lambda e"]
  S --> C["Adaptive command uc"]
  P["Projected parameter estimates"] --> C
  C --> SAT["Actuator saturation u = sat(uc)"]
  SAT --> PLANT["Uncertain position plant"]
  PLANT --> X
  C --> D["Deficiency delta_u = u - uc"]
  SAT --> D
  D --> AUX["Stable auxiliary state chi"]
  S --> SC["Compensated surface sc = s - chi"]
  AUX --> SC
  SC --> ADAPT["Projection-based update law"]
  ADAPT --> P
        

Consider the second-order SISO plant

\[ \ddot q = \boldsymbol{\theta}^{T}\boldsymbol{\phi}(q,\dot q,t)+b u, \qquad b>0, \]

where \( q \) is position, \( b \) is a known positive input gain, and \( \boldsymbol{\theta}\in\mathbb{R}^{p} \) is an unknown constant parameter vector. The regressor is known and measured. In the numerical example,

\[ \boldsymbol{\phi}(q,\dot q)= \begin{bmatrix}\sin q & \dot q\end{bmatrix}^{T}, \qquad \boldsymbol{\theta}= \begin{bmatrix}\theta_1 & \theta_2\end{bmatrix}^{T}. \]

The actuator magnitude is limited by \( u_{\max} \):

\[ u=\operatorname{sat}_{u_{\max}}(u_c)= \begin{cases} u_{\max}, & u_c>u_{\max},\\ u_c, & |u_c|\le u_{\max},\\ -u_{\max}, & u_c<-u_{\max}. \end{cases} \]

Define the saturation deficiency

\[ \Delta u \triangleq u-u_c. \]

The deficiency is zero in the unsaturated region and opposes the command outside that region. In particular,

\[ \Delta u\,u_c\le 0, \qquad |u|\le u_{\max}. \]

Let the desired trajectory be twice differentiable and define

\[ e=q-q_d, \qquad \dot e=\dot q-\dot q_d, \qquad s=\dot e+\lambda e, \quad \lambda>0. \]

The nominal adaptive command is

\[ u_c=\frac{1}{b}\left[-\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +\ddot q_d-\lambda\dot e-k s\right], \qquad k>0. \]

With the parameter error \( \tilde{\boldsymbol{\theta}}=\hat{\boldsymbol{\theta}}-\boldsymbol{\theta} \), substitution into the plant gives

\[ \dot s=-k s-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}+b\Delta u. \]

If adaptation were driven directly by \( s \), the term \( b\Delta u \) would be interpreted as parameter mismatch. During saturation, the estimator would therefore try to compensate for control authority that the actuator physically cannot deliver. This is the adaptive analogue of integrator windup.

Introduce the stable auxiliary state

\[ \dot\chi=-k\chi+b\Delta u, \qquad \chi(0)=0, \]

and define the compensated surface

\[ s_c=s-\chi. \]

Subtracting the auxiliary dynamics from the raw surface dynamics yields

\[ \dot s_c=-k s_c-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}. \]

The saturation term cancels exactly from the adaptation error model. The state \( \chi \) records the component of the raw tracking surface attributable to unavailable actuator action. Consequently, adaptation reacts to model uncertainty but does not chase the saturation deficiency.

This cancellation does not create control authority. It prevents estimator corruption while saturation is active. Actual tracking still depends on whether the desired motion is feasible for the bounded actuator.

Assume the true parameter belongs to a known compact convex set \( \Theta \). Use

\[ \dot{\hat{\boldsymbol{\theta}}}= \operatorname{Proj}_{\Theta}\!\left( \hat{\boldsymbol{\theta}},\;\boldsymbol{\Gamma}\boldsymbol{\phi}s_c \right), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0. \]

The projection operator is selected to satisfy

\[ \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\Gamma}^{-1} \left[\operatorname{Proj}_{\Theta}(\hat{\boldsymbol{\theta}},\mathbf{y})- \mathbf{y}\right]\le 0, \qquad \boldsymbol{\theta}\in\Theta. \]

For a rectangular set, a simple implementation removes any update component that points outward when the corresponding estimate is on a boundary. Projection guarantees \( \hat{\boldsymbol{\theta}}(t)\in\Theta \) for all time and bounds the parameter-dependent part of the command.

Theorem 1 (compensated-error convergence). Suppose:

1. \( \boldsymbol{\theta}\in\Theta \) and the projection property above holds;
2. \( q_d,\dot q_d,\ddot q_d \) and the regressor are bounded along the solution;
3. the closed-loop solution exists for all \( t\ge 0 \).

Then \( s_c \) and \( \tilde{\boldsymbol{\theta}} \) are bounded, \( s_c\in L_2\cap L_\infty \), and \( s_c(t)\to 0 \). If saturation ceases after a finite time, then \( \chi(t)\to 0 \), \( s(t)\to 0 \), and \( e(t)\to 0 \).

Proof. Choose

\[ V=\frac{1}{2}s_c^2+\frac{1}{2} \tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}}. \]

Differentiating and using the compensated dynamics gives

\[ \begin{aligned} \dot V &=s_c\left(-k s_c-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}\right) +\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\Gamma}^{-1} \operatorname{Proj}_{\Theta}\!\left( \hat{\boldsymbol{\theta}},\boldsymbol{\Gamma}\boldsymbol{\phi}s_c\right)\\ &\le -k s_c^2-s_c\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}s_c\\ &=-k s_c^2\le 0. \end{aligned} \]

Therefore,

\[ V(t)+k\int_{0}^{t}s_c^2(\tau)\,d\tau\le V(0), \]

so \( s_c \) and the parameter error are bounded and \( s_c\in L_2 \). Under the bounded-signal assumptions, \( \dot s_c \) is bounded; hence Barbalat's lemma gives \( s_c(t)\to 0 \).

If \( \Delta u(t)=0 \) after some finite time \( T_s \), then

\[ \chi(t)=\chi(T_s)e^{-k(t-T_s)}\to 0. \]

Thus \( s=s_c+\chi\to 0 \). Finally, \( \dot e=-\lambda e+s \) is an exponentially stable linear filter driven by a vanishing input, so \( e(t)\to 0 \). \( \square \)

The theorem intentionally does not claim zero physical tracking error while the requested motion remains infeasible. If the deficiency approaches a nonzero constant \( \Delta u_0 \), then the auxiliary state approaches

\[ \chi_{\mathrm{ss}}=\frac{b\Delta u_0}{k}. \]

Since \( s_c\to 0 \), the raw surface approaches \( s\to\chi_{\mathrm{ss}} \). For a constant limiting surface, the stable error filter gives

\[ e_{\mathrm{ss}}=\frac{\chi_{\mathrm{ss}}}{\lambda} =\frac{b\Delta u_0}{k\lambda}. \]

This result is not a failure of the adaptive law. It quantifies the tracking offset caused by insufficient actuator authority. The proper remedies are to reduce reference acceleration, increase available actuator authority, or use a supervisory reference governor; increasing adaptation gain cannot remove a hard input limit.

The command remains unsaturated when

\[ \left| -\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi} +\ddot q_d-\lambda\dot e-k s \right|\le b u_{\max}. \]

A conservative online reference-acceleration margin can be computed from the projection set. Let \( \bar\theta=\max_{\boldsymbol{\vartheta}\in\Theta}\|\boldsymbol{\vartheta}\| \). A sufficient condition is

\[ |\ddot q_d|\le b u_{\max}-\bar\theta\|\boldsymbol{\phi}\| -\lambda|\dot e|-k|s|. \]

A reference governor can reduce \( \ddot q_d \) whenever the right-hand side becomes small or negative.

• \( \lambda \): sets the nominal error-filter bandwidth.
• \( k \): increases surface damping and shortens the auxiliary-state time constant \( 1/k \).
• \( \boldsymbol{\Gamma} \): sets adaptation speed; excessive values amplify noise and discretization error.
• Projection set: should contain all physically plausible parameters but should not be unnecessarily large.
• Reference limits: should be enforced before the adaptive loop, not left to saturation alone.
• Independent safety layer: must enforce position, velocity, current, temperature, and emergency-stop limits.

The implementations use

\[ \boldsymbol{\theta}= \begin{bmatrix}-2 & -0.8\end{bmatrix}^{T}, \quad b=1, \quad u_{\max}=2.5, \quad \lambda=2, \quad k=4, \]

\[ \boldsymbol{\Gamma}=\operatorname{diag}(4,2), \qquad q_d(t)=1.2\sin(0.5t). \]

The parameter box is

\[ -5\le\hat\theta_1\le 1, \qquad -3\le\hat\theta_2\le 1. \]

A fixed step of \( h=10^{-3}\,\mathrm{s} \) is used. The actuator reaches its limit during the initial transient, after which the command becomes feasible and the ordinary tracking error converges close to zero.

flowchart TD
  A["Read q, q_dot and reference"] --> B["Compute e, e_dot and surface s"]
  B --> C["Compute unconstrained command uc"]
  C --> D["Apply u = sat(uc)"]
  D --> E["Compute deficiency delta_u = u - uc"]
  E --> F["Update auxiliary state chi"]
  F --> G["Compute compensated surface sc = s - chi"]
  G --> H["Compute gradient parameter update"]
  H --> I["Project update onto parameter box"]
  I --> J["Integrate plant and repeat"]
        

This implementation uses NumPy for array operations and Matplotlib for plots. The projection and simulation loop are written explicitly so that every adaptive control calculation remains visible.

Chapter19_Lesson5.py

"""Chapter19_Lesson5.py
Saturation-aware adaptive position control with an auxiliary anti-windup state.

Dependencies:
    pip install numpy matplotlib
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Config:
    dt: float = 1.0e-3
    final_time: float = 20.0
    input_gain: float = 1.0
    u_max: float = 2.5
    lambda_error: float = 2.0
    k_surface: float = 4.0
    gamma_1: float = 4.0
    gamma_2: float = 2.0
    theta_lower_1: float = -5.0
    theta_lower_2: float = -3.0
    theta_upper_1: float = 1.0
    theta_upper_2: float = 1.0


def saturate(value: float, limit: float) -> float:
    return max(-limit, min(limit, value))


def projected_derivative(
    estimate: np.ndarray,
    raw_derivative: np.ndarray,
    lower: np.ndarray,
    upper: np.ndarray,
) -> np.ndarray:
    """Box projection: remove derivative components pointing outside the set."""
    result = raw_derivative.copy()
    for index in range(estimate.size):
        at_lower_and_outward = estimate[index] <= lower[index] and result[index] < 0.0
        at_upper_and_outward = estimate[index] >= upper[index] and result[index] > 0.0
        if at_lower_and_outward or at_upper_and_outward:
            result[index] = 0.0
    return result


def reference(t: float) -> tuple[float, float, float]:
    q_d = 1.2 * math.sin(0.5 * t)
    q_d_dot = 0.6 * math.cos(0.5 * t)
    q_d_ddot = -0.3 * math.sin(0.5 * t)
    return q_d, q_d_dot, q_d_ddot


def simulate(config: Config) -> dict[str, np.ndarray]:
    steps = int(round(config.final_time / config.dt)) + 1
    time = np.linspace(0.0, config.final_time, steps)

    q = np.zeros(steps)
    q_dot = np.zeros(steps)
    q_d = np.zeros(steps)
    u = np.zeros(steps)
    u_command = np.zeros(steps)
    chi = np.zeros(steps)
    s_compensated = np.zeros(steps)
    theta_hat = np.zeros((steps, 2))

    theta_true = np.array([-2.0, -0.8])
    theta_lower = np.array([config.theta_lower_1, config.theta_lower_2])
    theta_upper = np.array([config.theta_upper_1, config.theta_upper_2])
    gamma = np.array([config.gamma_1, config.gamma_2])

    for index in range(steps - 1):
        t = time[index]
        q_d[index], q_d_dot, q_d_ddot = reference(t)

        error = q[index] - q_d[index]
        error_dot = q_dot[index] - q_d_dot
        surface = error_dot + config.lambda_error * error
        regressor = np.array([math.sin(q[index]), q_dot[index]])

        u_command[index] = (
            -float(theta_hat[index] @ regressor)
            + q_d_ddot
            - config.lambda_error * error_dot
            - config.k_surface * surface
        ) / config.input_gain

        u[index] = saturate(u_command[index], config.u_max)
        saturation_deficiency = u[index] - u_command[index]

        chi_dot = -config.k_surface * chi[index] + config.input_gain * saturation_deficiency
        s_compensated[index] = surface - chi[index]

        raw_theta_dot = gamma * regressor * s_compensated[index]
        theta_dot = projected_derivative(
            theta_hat[index], raw_theta_dot, theta_lower, theta_upper
        )

        q_ddot = float(theta_true @ regressor) + config.input_gain * u[index]

        # Semi-implicit Euler integration is robust for this small fixed time step.
        q_dot[index + 1] = q_dot[index] + config.dt * q_ddot
        q[index + 1] = q[index] + config.dt * q_dot[index + 1]
        chi[index + 1] = chi[index] + config.dt * chi_dot
        theta_hat[index + 1] = np.clip(
            theta_hat[index] + config.dt * theta_dot,
            theta_lower,
            theta_upper,
        )

    q_d[-1], q_d_dot, q_d_ddot = reference(time[-1])
    error = q[-1] - q_d[-1]
    error_dot = q_dot[-1] - q_d_dot
    surface = error_dot + config.lambda_error * error
    regressor = np.array([math.sin(q[-1]), q_dot[-1]])
    u_command[-1] = (
        -float(theta_hat[-1] @ regressor)
        + q_d_ddot
        - config.lambda_error * error_dot
        - config.k_surface * surface
    ) / config.input_gain
    u[-1] = saturate(u_command[-1], config.u_max)
    s_compensated[-1] = surface - chi[-1]

    return {
        "time": time,
        "q": q,
        "q_dot": q_dot,
        "q_d": q_d,
        "u": u,
        "u_command": u_command,
        "chi": chi,
        "s_compensated": s_compensated,
        "theta_hat_1": theta_hat[:, 0],
        "theta_hat_2": theta_hat[:, 1],
    }


def save_csv(results: dict[str, np.ndarray], path: Path) -> None:
    keys = list(results.keys())
    with path.open("w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(keys)
        for row in zip(*(results[key] for key in keys)):
            writer.writerow(row)


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

    figure_1 = plt.figure()
    plt.plot(time, results["q_d"], "--", label="desired position")
    plt.plot(time, results["q"], label="actual position")
    plt.xlabel("time [s]")
    plt.ylabel("position")
    plt.grid(True)
    plt.legend()
    figure_1.tight_layout()

    figure_2 = plt.figure()
    plt.plot(time, results["u_command"], label="commanded input")
    plt.plot(time, results["u"], label="applied saturated input")
    plt.axhline(2.5, linestyle=":")
    plt.axhline(-2.5, linestyle=":")
    plt.xlabel("time [s]")
    plt.ylabel("control input")
    plt.grid(True)
    plt.legend()
    figure_2.tight_layout()

    figure_3 = plt.figure()
    plt.plot(time, results["theta_hat_1"], label="theta_hat_1")
    plt.plot(time, results["theta_hat_2"], label="theta_hat_2")
    plt.axhline(-2.0, linestyle=":", label="theta_1 true")
    plt.axhline(-0.8, linestyle=":", label="theta_2 true")
    plt.xlabel("time [s]")
    plt.ylabel("parameter estimate")
    plt.grid(True)
    plt.legend()
    figure_3.tight_layout()

    figure_4 = plt.figure()
    plt.plot(time, results["s_compensated"], label="compensated surface")
    plt.plot(time, results["chi"], label="saturation state chi")
    plt.xlabel("time [s]")
    plt.ylabel("error coordinates")
    plt.grid(True)
    plt.legend()
    figure_4.tight_layout()

    plt.show()


def main() -> None:
    config = Config()
    results = simulate(config)
    output_path = Path(__file__).with_name("Chapter19_Lesson5_python_results.csv")
    save_csv(results, output_path)
    print(f"Saved simulation data to: {output_path}")
    print(f"Final position error: {results['q'][-1] - results['q_d'][-1]:.6f}")
    print(f"Peak applied input: {np.max(np.abs(results['u'])):.6f}")
    plot_results(results)


if __name__ == "__main__":
    main()

The C++ version uses only the standard library and writes results to CSV. Eigen can be substituted for fixed-size arrays when the parameter vector is larger.

Chapter19_Lesson5.cpp

// Chapter19_Lesson5.cpp
// Saturation-aware adaptive position control using only the C++ standard library.

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

namespace {
constexpr double kDt = 1.0e-3;
constexpr double kFinalTime = 20.0;
constexpr double kInputGain = 1.0;
constexpr double kUMax = 2.5;
constexpr double kLambda = 2.0;
constexpr double kSurfaceGain = 4.0;
constexpr std::array<double, 2> kGamma{4.0, 2.0};
constexpr std::array<double, 2> kThetaTrue{-2.0, -0.8};
constexpr std::array<double, 2> kThetaLower{-5.0, -3.0};
constexpr std::array<double, 2> kThetaUpper{1.0, 1.0};

double saturate(const double value, const double limit) {
    return std::clamp(value, -limit, limit);
}

std::array<double, 3> reference(const double t) {
    return {
        1.2 * std::sin(0.5 * t),
        0.6 * std::cos(0.5 * t),
        -0.3 * std::sin(0.5 * t),
    };
}

double projectedDerivative(
    const double estimate,
    const double rawDerivative,
    const double lower,
    const double upper) {
    const bool lowerOutward = estimate <= lower && rawDerivative < 0.0;
    const bool upperOutward = estimate >= upper && rawDerivative > 0.0;
    return (lowerOutward || upperOutward) ? 0.0 : rawDerivative;
}
}  // namespace

int main() {
    std::ofstream csv("Chapter19_Lesson5_cpp_results.csv");
    if (!csv) {
        std::cerr << "Could not open the output CSV file.\n";
        return 1;
    }
    csv << std::setprecision(12);
    csv << "time,q,q_desired,u,u_command,chi,s_compensated,theta_hat_1,theta_hat_2\n";

    double q = 0.0;
    double qDot = 0.0;
    double chi = 0.0;
    std::array<double, 2> thetaHat{0.0, 0.0};

    const int steps = static_cast<int>(std::lround(kFinalTime / kDt));
    double finalError = 0.0;
    double peakAppliedInput = 0.0;

    for (int index = 0; index <= steps; ++index) {
        const double t = index * kDt;
        const auto desired = reference(t);
        const double qDesired = desired[0];
        const double qDesiredDot = desired[1];
        const double qDesiredDDot = desired[2];

        const double error = q - qDesired;
        const double errorDot = qDot - qDesiredDot;
        const double surface = errorDot + kLambda * error;
        const std::array<double, 2> phi{std::sin(q), qDot};

        const double estimatedDynamics = thetaHat[0] * phi[0] + thetaHat[1] * phi[1];
        const double uCommand =
            (-estimatedDynamics + qDesiredDDot - kLambda * errorDot - kSurfaceGain * surface)
            / kInputGain;
        const double u = saturate(uCommand, kUMax);
        const double saturationDeficiency = u - uCommand;
        const double compensatedSurface = surface - chi;

        csv << t << ',' << q << ',' << qDesired << ',' << u << ',' << uCommand << ','
            << chi << ',' << compensatedSurface << ',' << thetaHat[0] << ',' << thetaHat[1]
            << '\n';

        finalError = error;
        peakAppliedInput = std::max(peakAppliedInput, std::abs(u));
        if (index == steps) {
            break;
        }

        const double chiDot = -kSurfaceGain * chi + kInputGain * saturationDeficiency;
        std::array<double, 2> thetaDot{};
        for (std::size_t j = 0; j < thetaDot.size(); ++j) {
            const double rawDerivative = kGamma[j] * phi[j] * compensatedSurface;
            thetaDot[j] = projectedDerivative(
                thetaHat[j], rawDerivative, kThetaLower[j], kThetaUpper[j]);
        }

        const double qDDot = kThetaTrue[0] * phi[0] + kThetaTrue[1] * phi[1] + kInputGain * u;

        qDot += kDt * qDDot;
        q += kDt * qDot;
        chi += kDt * chiDot;
        for (std::size_t j = 0; j < thetaHat.size(); ++j) {
            thetaHat[j] = std::clamp(
                thetaHat[j] + kDt * thetaDot[j], kThetaLower[j], kThetaUpper[j]);
        }
    }

    std::cout << std::fixed << std::setprecision(6)
              << "Final position error: " << finalError << '\n'
              << "Peak applied input: " << peakAppliedInput << '\n'
              << "Saved: Chapter19_Lesson5_cpp_results.csv\n";
    return 0;
}

The Java version uses primitive arrays and the standard NIO file API. For larger MIMO designs, EJML is a suitable matrix library, but no external dependency is required for this two-parameter example.

Chapter19_Lesson5.java

// Chapter19_Lesson5.java
// Saturation-aware adaptive position control using the Java standard library.

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 Chapter19_Lesson5 {
    private static final double DT = 1.0e-3;
    private static final double FINAL_TIME = 20.0;
    private static final double INPUT_GAIN = 1.0;
    private static final double U_MAX = 2.5;
    private static final double LAMBDA = 2.0;
    private static final double SURFACE_GAIN = 4.0;

    private static final double[] GAMMA = {4.0, 2.0};
    private static final double[] THETA_TRUE = {-2.0, -0.8};
    private static final double[] THETA_LOWER = {-5.0, -3.0};
    private static final double[] THETA_UPPER = {1.0, 1.0};

    private Chapter19_Lesson5() {
    }

    private static double saturate(double value, double limit) {
        return Math.max(-limit, Math.min(limit, value));
    }

    private static double[] reference(double time) {
        return new double[] {
            1.2 * Math.sin(0.5 * time),
            0.6 * Math.cos(0.5 * time),
            -0.3 * Math.sin(0.5 * time)
        };
    }

    private static double projectedDerivative(
            double estimate,
            double rawDerivative,
            double lower,
            double upper) {
        boolean lowerOutward = estimate <= lower && rawDerivative < 0.0;
        boolean upperOutward = estimate >= upper && rawDerivative > 0.0;
        return (lowerOutward || upperOutward) ? 0.0 : rawDerivative;
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        Path output = Path.of("Chapter19_Lesson5_java_results.csv");

        double q = 0.0;
        double qDot = 0.0;
        double chi = 0.0;
        double[] thetaHat = {0.0, 0.0};

        int steps = (int) Math.round(FINAL_TIME / DT);
        double finalError = 0.0;
        double peakAppliedInput = 0.0;

        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("time,q,q_desired,u,u_command,chi,s_compensated,theta_hat_1,theta_hat_2");
            writer.newLine();

            for (int index = 0; index <= steps; index++) {
                double time = index * DT;
                double[] desired = reference(time);
                double qDesired = desired[0];
                double qDesiredDot = desired[1];
                double qDesiredDDot = desired[2];

                double error = q - qDesired;
                double errorDot = qDot - qDesiredDot;
                double surface = errorDot + LAMBDA * error;
                double[] phi = {Math.sin(q), qDot};

                double estimatedDynamics = thetaHat[0] * phi[0] + thetaHat[1] * phi[1];
                double uCommand = (
                        -estimatedDynamics
                                + qDesiredDDot
                                - LAMBDA * errorDot
                                - SURFACE_GAIN * surface)
                        / INPUT_GAIN;
                double u = saturate(uCommand, U_MAX);
                double saturationDeficiency = u - uCommand;
                double compensatedSurface = surface - chi;

                writer.write(String.format(
                        Locale.US,
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f",
                        time,
                        q,
                        qDesired,
                        u,
                        uCommand,
                        chi,
                        compensatedSurface,
                        thetaHat[0],
                        thetaHat[1]));
                writer.newLine();

                finalError = error;
                peakAppliedInput = Math.max(peakAppliedInput, Math.abs(u));
                if (index == steps) {
                    break;
                }

                double chiDot = -SURFACE_GAIN * chi + INPUT_GAIN * saturationDeficiency;
                double[] thetaDot = new double[2];
                for (int j = 0; j < thetaDot.length; j++) {
                    double rawDerivative = GAMMA[j] * phi[j] * compensatedSurface;
                    thetaDot[j] = projectedDerivative(
                            thetaHat[j], rawDerivative, THETA_LOWER[j], THETA_UPPER[j]);
                }

                double qDDot = THETA_TRUE[0] * phi[0]
                        + THETA_TRUE[1] * phi[1]
                        + INPUT_GAIN * u;

                qDot += DT * qDDot;
                q += DT * qDot;
                chi += DT * chiDot;
                for (int j = 0; j < thetaHat.length; j++) {
                    thetaHat[j] = Math.max(
                            THETA_LOWER[j],
                            Math.min(THETA_UPPER[j], thetaHat[j] + DT * thetaDot[j]));
                }
            }
        }

        System.out.printf(Locale.US, "Final position error: %.6f%n", finalError);
        System.out.printf(Locale.US, "Peak applied input: %.6f%n", peakAppliedInput);
        System.out.println("Saved: " + output.toAbsolutePath());
    }
}

The MATLAB script performs the numerical simulation, creates plots, and contains an optional programmatic Simulink builder. When the flag buildSimulink is set to true, the script creates a fixed-step model whose MATLAB Function block implements the same controller and plant update.

Chapter19_Lesson5.m

% Chapter19_Lesson5.m
% Saturation-aware adaptive position control and optional Simulink model builder.

clear; close all; clc;

cfg.dt = 1.0e-3;
cfg.finalTime = 20.0;
cfg.inputGain = 1.0;
cfg.uMax = 2.5;
cfg.lambda = 2.0;
cfg.k = 4.0;
cfg.gamma = [4.0; 2.0];
cfg.thetaTrue = [-2.0; -0.8];
cfg.thetaLower = [-5.0; -3.0];
cfg.thetaUpper = [1.0; 1.0];

N = round(cfg.finalTime / cfg.dt) + 1;
time = (0:N-1)' * cfg.dt;
q = zeros(N, 1);
qDot = zeros(N, 1);
qDesired = zeros(N, 1);
u = zeros(N, 1);
uCommand = zeros(N, 1);
chi = zeros(N, 1);
sCompensated = zeros(N, 1);
thetaHat = zeros(N, 2);

for index = 1:N-1
    t = time(index);
    qDesired(index) = 1.2 * sin(0.5 * t);
    qDesiredDot = 0.6 * cos(0.5 * t);
    qDesiredDDot = -0.3 * sin(0.5 * t);

    error = q(index) - qDesired(index);
    errorDot = qDot(index) - qDesiredDot;
    surface = errorDot + cfg.lambda * error;
    phi = [sin(q(index)); qDot(index)];

    uCommand(index) = (-thetaHat(index, :) * phi + qDesiredDDot ...
        - cfg.lambda * errorDot - cfg.k * surface) / cfg.inputGain;
    u(index) = min(cfg.uMax, max(-cfg.uMax, uCommand(index)));
    saturationDeficiency = u(index) - uCommand(index);

    chiDot = -cfg.k * chi(index) + cfg.inputGain * saturationDeficiency;
    sCompensated(index) = surface - chi(index);

    rawThetaDot = cfg.gamma .* phi * sCompensated(index);
    thetaDot = projectedDerivative(thetaHat(index, :)', rawThetaDot, ...
        cfg.thetaLower, cfg.thetaUpper);

    qDDot = cfg.thetaTrue' * phi + cfg.inputGain * u(index);

    qDot(index + 1) = qDot(index) + cfg.dt * qDDot;
    q(index + 1) = q(index) + cfg.dt * qDot(index + 1);
    chi(index + 1) = chi(index) + cfg.dt * chiDot;
    thetaHat(index + 1, :) = min(cfg.thetaUpper, max(cfg.thetaLower, ...
        thetaHat(index, :)' + cfg.dt * thetaDot))';
end

qDesired(end) = 1.2 * sin(0.5 * time(end));
qDesiredDot = 0.6 * cos(0.5 * time(end));
qDesiredDDot = -0.3 * sin(0.5 * time(end));
error = q(end) - qDesired(end);
errorDot = qDot(end) - qDesiredDot;
surface = errorDot + cfg.lambda * error;
phi = [sin(q(end)); qDot(end)];
uCommand(end) = (-thetaHat(end, :) * phi + qDesiredDDot ...
    - cfg.lambda * errorDot - cfg.k * surface) / cfg.inputGain;
u(end) = min(cfg.uMax, max(-cfg.uMax, uCommand(end)));
sCompensated(end) = surface - chi(end);

results = table(time, q, qDesired, u, uCommand, chi, sCompensated, ...
    thetaHat(:, 1), thetaHat(:, 2), ...
    'VariableNames', {'time', 'q', 'qDesired', 'u', 'uCommand', 'chi', ...
    'sCompensated', 'thetaHat1', 'thetaHat2'});
writetable(results, 'Chapter19_Lesson5_matlab_results.csv');

figure;
plot(time, qDesired, '--', time, q, 'LineWidth', 1.2);
grid on; xlabel('time [s]'); ylabel('position');
legend('desired position', 'actual position', 'Location', 'best');

figure;
plot(time, uCommand, time, u, 'LineWidth', 1.1); hold on;
yline(cfg.uMax, ':'); yline(-cfg.uMax, ':');
grid on; xlabel('time [s]'); ylabel('control input');
legend('commanded input', 'applied saturated input', 'Location', 'best');

figure;
plot(time, thetaHat(:, 1), time, thetaHat(:, 2), 'LineWidth', 1.1); hold on;
yline(cfg.thetaTrue(1), ':'); yline(cfg.thetaTrue(2), ':');
grid on; xlabel('time [s]'); ylabel('parameter estimate');
legend('thetaHat1', 'thetaHat2', 'theta1 true', 'theta2 true', 'Location', 'best');

figure;
plot(time, sCompensated, time, chi, 'LineWidth', 1.1);
grid on; xlabel('time [s]'); ylabel('error coordinates');
legend('compensated surface', 'saturation state chi', 'Location', 'best');

fprintf('Final position error: %.6f\n', q(end) - qDesired(end));
fprintf('Peak applied input: %.6f\n', max(abs(u)));

% Set this flag to true when Simulink is installed. The generated model uses
% one discrete MATLAB Function block so that the numerical algorithm remains
% identical to the script implementation.
buildSimulink = false;
if buildSimulink
    buildChapter19Lesson5Simulink(cfg);
end

function derivative = projectedDerivative(estimate, rawDerivative, lower, upper)
    derivative = rawDerivative;
    for index = 1:numel(estimate)
        lowerOutward = estimate(index) <= lower(index) && derivative(index) < 0;
        upperOutward = estimate(index) >= upper(index) && derivative(index) > 0;
        if lowerOutward || upperOutward
            derivative(index) = 0;
        end
    end
end

function buildChapter19Lesson5Simulink(cfg)
    model = 'Chapter19_Lesson5_Simulink';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    if isfile([model '.slx'])
        delete([model '.slx']);
    end

    new_system(model);
    open_system(model);
    set_param(model, 'SolverType', 'Fixed-step', ...
        'Solver', 'FixedStepDiscrete', ...
        'FixedStep', num2str(cfg.dt), ...
        'StopTime', num2str(cfg.finalTime));

    add_block('simulink/Sources/Clock', [model '/Clock'], ...
        'Position', [40 80 70 110]);
    add_block('simulink/User-Defined Functions/MATLAB Function', ...
        [model '/AdaptivePlantStep'], 'Position', [130 45 340 145]);
    add_block('simulink/Signal Routing/Demux', [model '/Demux'], ...
        'Outputs', '6', 'Position', [390 45 395 175]);

    outputNames = {'q', 'qDesired', 'u', 'uCommand', 'thetaHat1', 'thetaHat2'};
    for index = 1:numel(outputNames)
        add_block('simulink/Sinks/To Workspace', ...
            [model '/' outputNames{index}], ...
            'VariableName', outputNames{index}, ...
            'SaveFormat', 'Array', ...
            'Position', [470 20 + 45 * index 585 45 + 45 * index]);
        add_line(model, ['Demux/' num2str(index)], ...
            [outputNames{index} '/1'], 'autorouting', 'on');
    end
    add_line(model, 'Clock/1', 'AdaptivePlantStep/1', 'autorouting', 'on');
    add_line(model, 'AdaptivePlantStep/1', 'Demux/1', 'autorouting', 'on');

    root = sfroot;
    chart = root.find('-isa', 'Stateflow.EMChart', ...
        'Path', [model '/AdaptivePlantStep']);
    chart.Script = sprintf([ ...
        'function y = fcn(t)\n' ...
        '%%#codegen\n' ...
        'persistent q qDot chi thetaHat initialized\n' ...
        'dt = %.17g; b = %.17g; uMax = %.17g; lambda = %.17g; k = %.17g;\n' ...
        'gamma = [%.17g; %.17g]; thetaTrue = [%.17g; %.17g];\n' ...
        'lower = [%.17g; %.17g]; upper = [%.17g; %.17g];\n' ...
        'if isempty(initialized)\n' ...
        '  q = 0; qDot = 0; chi = 0; thetaHat = [0;0]; initialized = true;\n' ...
        'end\n' ...
        'qd = 1.2*sin(0.5*t); qdDot = 0.6*cos(0.5*t); qdDDot = -0.3*sin(0.5*t);\n' ...
        'e = q-qd; eDot = qDot-qdDot; s = eDot+lambda*e; phi = [sin(q); qDot];\n' ...
        'uc = (-thetaHat''*phi+qdDDot-lambda*eDot-k*s)/b;\n' ...
        'u = min(uMax,max(-uMax,uc)); delta = u-uc; sa = s-chi;\n' ...
        'raw = gamma.*phi*sa; dtheta = raw;\n' ...
        'for j=1:2\n' ...
        ' if (thetaHat(j)<=lower(j) && dtheta(j)<0) || (thetaHat(j)>=upper(j) && dtheta(j)>0)\n' ...
        '  dtheta(j)=0;\n' ...
        ' end\n' ...
        'end\n' ...
        'qDDot = thetaTrue''*phi+b*u; qDot=qDot+dt*qDDot; q=q+dt*qDot;\n' ...
        'chi=chi+dt*(-k*chi+b*delta); thetaHat=min(upper,max(lower,thetaHat+dt*dtheta));\n' ...
        'y=[q;qd;u;uc;thetaHat(1);thetaHat(2)];\n' ...
        'end\n'], ...
        cfg.dt, cfg.inputGain, cfg.uMax, cfg.lambda, cfg.k, ...
        cfg.gamma(1), cfg.gamma(2), cfg.thetaTrue(1), cfg.thetaTrue(2), ...
        cfg.thetaLower(1), cfg.thetaLower(2), cfg.thetaUpper(1), cfg.thetaUpper(2));

    save_system(model);
    fprintf('Created %s.slx\n', model);
end

The notebook uses built-in list operations, Clip, CSV export, and ListLinePlot. The downloadable archive also includes the readable Wolfram source separately from the notebook wrapper.

Chapter19_Lesson5.nb

(* Chapter19_Lesson5.nb
   Wolfram Language implementation of saturation-aware adaptive position control.
   The file is intentionally plain-text Wolfram Language and can be evaluated in
   Mathematica; it is also wrapped as a notebook by the generated Notebook expression.
*)

ClearAll["Global`*"];

dt = 0.001;
finalTime = 20.0;
inputGain = 1.0;
uMax = 2.5;
lambda = 2.0;
kSurface = 4.0;
gamma = {4.0, 2.0};
thetaTrue = {-2.0, -0.8};
thetaLower = {-5.0, -3.0};
thetaUpper = {1.0, 1.0};

saturate[value_] := Clip[value, {-uMax, uMax}];

reference[t_] := {
  1.2 Sin[0.5 t],
  0.6 Cos[0.5 t],
  -0.3 Sin[0.5 t]
};

projectedDerivative[estimate_, raw_] := MapThread[
  Function[{est, deriv, lower, upper},
    If[(est <= lower && deriv < 0) || (est >= upper && deriv > 0), 0.0, deriv]
  ],
  {estimate, raw, thetaLower, thetaUpper}
];

steps = Round[finalTime/dt];
state = <|
  "q" -> 0.0,
  "qDot" -> 0.0,
  "chi" -> 0.0,
  "thetaHat" -> {0.0, 0.0}
|>;

history = Reap[
  Do[
    t = index dt;
    {qDesired, qDesiredDot, qDesiredDDot} = reference[t];
    q = state["q"];
    qDot = state["qDot"];
    chi = state["chi"];
    thetaHat = state["thetaHat"];

    error = q - qDesired;
    errorDot = qDot - qDesiredDot;
    surface = errorDot + lambda error;
    phi = {Sin[q], qDot};

    uCommand = (-thetaHat . phi + qDesiredDDot - lambda errorDot -
        kSurface surface)/inputGain;
    u = saturate[uCommand];
    saturationDeficiency = u - uCommand;
    compensatedSurface = surface - chi;

    Sow[{t, q, qDesired, u, uCommand, chi, compensatedSurface,
      thetaHat[[1]], thetaHat[[2]]}];

    If[index < steps,
      chiDot = -kSurface chi + inputGain saturationDeficiency;
      rawThetaDot = gamma phi compensatedSurface;
      thetaDot = projectedDerivative[thetaHat, rawThetaDot];
      qDDot = thetaTrue . phi + inputGain u;

      qDotNext = qDot + dt qDDot;
      qNext = q + dt qDotNext;
      chiNext = chi + dt chiDot;
      thetaNext = Clip[thetaHat + dt thetaDot, {thetaLower, thetaUpper}];
      state = <|"q" -> qNext, "qDot" -> qDotNext, "chi" -> chiNext,
        "thetaHat" -> thetaNext|>;
    ],
    {index, 0, steps}
  ]
][[2, 1]];

header = {"time", "q", "qDesired", "u", "uCommand", "chi",
  "sCompensated", "thetaHat1", "thetaHat2"};
Export["Chapter19_Lesson5_mathematica_results.csv", Prepend[history, header]];

positionPlot = ListLinePlot[
  {history[[All, {1, 2}]], history[[All, {1, 3}]]},
  PlotLegends -> {"actual position", "desired position"},
  Frame -> True,
  FrameLabel -> {"time [s]", "position"},
  PlotRange -> All,
  ImageSize -> Large
];

inputPlot = ListLinePlot[
  {history[[All, {1, 4}]], history[[All, {1, 5}]]},
  PlotLegends -> {"applied saturated input", "commanded input"},
  Frame -> True,
  FrameLabel -> {"time [s]", "control input"},
  PlotRange -> All,
  ImageSize -> Large
];

parameterPlot = ListLinePlot[
  {history[[All, {1, 8}]], history[[All, {1, 9}]]},
  PlotLegends -> {"thetaHat1", "thetaHat2"},
  Frame -> True,
  FrameLabel -> {"time [s]", "parameter estimate"},
  PlotRange -> All,
  ImageSize -> Large
];

Print["Final position error: ", Last[history][[2]] - Last[history][[3]]];
Print["Peak applied input: ", Max[Abs[history[[All, 4]]]]];
Column[{positionPlot, inputPlot, parameterPlot}]

A correct implementation should satisfy all of the following:

1. Hard input bound: numerically verify \( |u(t)|\le 2.5 \) at every sample.
2. Projection invariance: verify each parameter estimate remains inside its specified interval.
3. Saturation isolation: during saturation, \( \chi \) becomes nonzero while \( s_c \) remains the adaptation-driving error.
4. Post-saturation recovery: after \( \Delta u \) returns to zero, \( \chi \) decays exponentially.
5. Tracking: the position follows the sinusoidal reference after the transient.
6. No parameter-convergence claim without excitation: tracking may converge even when estimates do not exactly equal the true parameters.

On the tested implementation, the peak applied input is exactly the specified limit and the final position error is approximately \( 5.5\times10^{-3} \). Exact values can vary slightly with integration method and floating-point implementation.

Problem 1 (derive the compensated dynamics). Starting from

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

and \( \dot\chi=-k\chi+b\Delta u \), derive the dynamics of \( s_c=s-\chi \).

Solution.

\[ \begin{aligned} \dot s_c&=\dot s-\dot\chi\\ &=\left(-ks-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}+b\Delta u\right) -\left(-k\chi+b\Delta u\right)\\ &=-k(s-\chi)-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}\\ &=-ks_c-\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}. \end{aligned} \]

The saturation deficiency cancels exactly.

Problem 2 (persistent saturation offset). Suppose \( \Delta u(t)\to\Delta u_0 \) and \( s_c(t)\to0 \). Find the limiting values of \( \chi \), \( s \), and \( e \).

Solution. The stable auxiliary filter has equilibrium

\[ \chi_{\mathrm{ss}}=\frac{b\Delta u_0}{k}. \]

Since \( s=s_c+\chi \),

\[ s_{\mathrm{ss}}=\frac{b\Delta u_0}{k}. \]

At equilibrium, \( 0=-\lambda e_{\mathrm{ss}}+s_{\mathrm{ss}} \), hence

\[ e_{\mathrm{ss}}=\frac{b\Delta u_0}{k\lambda}. \]

Problem 3 (projection at a boundary). Let \( \hat\theta_1=1 \), with admissible interval \( -5\le\hat\theta_1\le1 \). Determine the projected derivative for raw updates \( y_1=0.7 \) and \( y_1=-0.7 \).

Solution. At the upper boundary, a positive derivative points outside:

\[ \operatorname{Proj}(1,0.7)=0. \]

A negative derivative points back into the set and is retained:

\[ \operatorname{Proj}(1,-0.7)=-0.7. \]

Problem 4 (reference-acceleration feasibility). Derive a sufficient inequality on \( |\ddot q_d| \) that guarantees \( |u_c|\le u_{\max} \), using the bound \( |\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}| \le\bar\theta\|\boldsymbol{\phi}\| \).

Solution. From the command equation and the triangle inequality,

\[ b|u_c|\le \bar\theta\|\boldsymbol{\phi}\|+|\ddot q_d| +\lambda|\dot e|+k|s|. \]

Therefore a sufficient condition for \( |u_c|\le u_{\max} \) is

\[ |\ddot q_d|\le b u_{\max}-\bar\theta\|\boldsymbol{\phi}\| -\lambda|\dot e|-k|s|. \]

Problem 5 (discrete auxiliary-state stability). Forward Euler discretization gives

\[ \chi_{k+1}=(1-kh)\chi_k+h b\Delta u_k. \]

Find the step-size condition for stability of the homogeneous recursion.

Solution. Discrete-time stability requires

\[ |1-kh|<1. \]

Thus

\[ 0<h<\frac{2}{k}. \]

This condition concerns only the auxiliary filter. In practice, the plant and adaptation dynamics generally require a substantially smaller step.

Problem 6 (parameter convergence). Does \( s_c(t)\to0 \) imply \( \hat{\boldsymbol{\theta}}(t)\to\boldsymbol{\theta} \)?

Solution. No. The Lyapunov proof establishes bounded parameter error and compensated tracking convergence. Exact parameter convergence requires sufficient excitation of the regressor, typically a persistent-excitation condition. A single low-frequency reference may not excite all parameter directions even though tracking is accurate.

A hard actuator limit breaks the nominal adaptive error equation by introducing a saturation-deficiency term. The auxiliary state \( \dot\chi=-k\chi+b\Delta u \) removes that term from the adaptation-driving surface. Projection keeps parameter estimates physically bounded, and the Lyapunov derivative satisfies \( \dot V\le-k s_c^2 \). Therefore the compensated surface converges to zero under standard boundedness assumptions. Zero physical tracking error follows when the reference eventually becomes feasible and saturation ends. Persistent infeasibility produces a quantifiable offset and must be handled through reference management or greater actuator authority, not faster adaptation.