Chapter 25: Applications in Robotics and Mechatronics

Lesson 5: Lessons Learned from Lab/Industrial Robotic Applications

This lesson consolidates the mathematical and implementation lessons that determine whether an adaptive robot controller remains useful outside a clean simulation. We connect Euler–Lagrange stability analysis to sampled execution, noisy sensing, actuator saturation, parameter projection, adaptation gating, safety supervision, commissioning tests, and quantitative acceptance criteria. A common industrial-style benchmark is implemented in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

1. Why a Stable Adaptive Law Can Still Fail on a Robot

A continuous-time Lyapunov proof establishes properties of an idealized closed loop. A laboratory or industrial robot adds sampled sensing, quantization, velocity reconstruction, communication delay, torque and current limits, drivetrain friction, flexible modes, thermal variation, collision logic, and software faults. These effects do not invalidate adaptive control as a field, but they change the engineering claim from “asymptotic convergence under an exact model” to “bounded, measurable, recoverable performance under a documented operating envelope.”

For an \(n\)-joint rigid manipulator, the standard model is

\[ \mathbf{ {M} }(\mathbf{ {q} })\ddot{ {\mathbf{ {q} } } } + \mathbf{ {C} }(\mathbf{ {q} },\dot{ {\mathbf{ {q} } } })\dot{ {\mathbf{ {q} } } } + \mathbf{ {g} }(\mathbf{ {q} }) + \mathbf{ {f} }(\dot{ {\mathbf{ {q} } } }) = \boldsymbol{ {\tau} } + \mathbf{ {d} }. \]

Practical deployment requires explicit bounds for the disturbance \(\mathbf{ {d} }\), unmodeled friction and elasticity, measurement errors, command delay, and the difference between requested and applied torque. The central lesson is that every omitted mechanism eventually reappears as an effective disturbance in the filtered-error dynamics.

2. Baseline Adaptive Robot Controller

Let the position error and reference velocity be \(\mathbf{ {e} }=\mathbf{ {q} }-\mathbf{ {q} }_d\) and

\[ \dot{ {\mathbf{ {q} } } }_r = \dot{ {\mathbf{ {q} } } }_d - \boldsymbol{ {\Lambda} }\mathbf{ {e} }, \qquad \mathbf{ {s} } = \dot{ {\mathbf{ {q} } } }-\dot{ {\mathbf{ {q} } } }_r = \dot{ {\mathbf{ {e} } } }+\boldsymbol{ {\Lambda} }\mathbf{ {e} }, \]

where \(\boldsymbol{ {\Lambda} }=\boldsymbol{ {\Lambda} }^T\succ0\). Robot dynamics are linear in a constant physical parameter vector:

\[ \mathbf{ {M} }(\mathbf{ {q} })\ddot{ {\mathbf{ {q} } } }_r + \mathbf{ {C} }(\mathbf{ {q} },\dot{ {\mathbf{ {q} } } })\dot{ {\mathbf{ {q} } } }_r + \mathbf{ {g} }(\mathbf{ {q} }) = \mathbf{ {Y} }(\mathbf{ {q} },\dot{ {\mathbf{ {q} } } }, \dot{ {\mathbf{ {q} } } }_r,\ddot{ {\mathbf{ {q} } } }_r)\boldsymbol{ {\theta} }. \]

The nominal Slotine–Li-type controller is

\[ \boldsymbol{ {\tau} }_c = \mathbf{ {Y} }\hat{ {\boldsymbol{ {\theta} } } }- \mathbf{ {K} }_s\mathbf{ {s} }, \qquad \dot{ {\hat{ {\boldsymbol{ {\theta} } } } } } = -\boldsymbol{ {\Gamma} }\mathbf{ {Y} }^T\mathbf{ {s} }. \]

With exact torque realization and the usual skew-symmetry property of \(\dot{ {\mathbf{ {M} } } }-2\mathbf{ {C} }\), the Lyapunov function

\[ V = \frac{ {1} }{ {2} }\mathbf{ {s} }^T\mathbf{ {M} }\mathbf{ {s} }+ \frac{ {1} }{ {2} }\tilde{ {\boldsymbol{ {\theta} } } }^T \boldsymbol{ {\Gamma} }^{ {-1} }\tilde{ {\boldsymbol{ {\theta} } } }, \qquad \tilde{ {\boldsymbol{ {\theta} } } }=\hat{ {\boldsymbol{ {\theta} } } }- \boldsymbol{ {\theta} }, \]

gives \(\dot V=-\mathbf{ {s} }^T\mathbf{ {K} }_s\mathbf{ {s} }\).

3. Robust Practical Modification and Ultimate-Bound Derivation

In a real actuator, \(\boldsymbol{ {\tau} }\neq\boldsymbol{ {\tau} }_c\). Define the implementation error \(\Delta\boldsymbol{ {\tau} }=\boldsymbol{ {\tau} }- \boldsymbol{ {\tau} }_c\) and collect disturbances, estimation errors, and unmodeled effects into \(\mathbf{ {w} }\). A deployment-oriented adaptive law uses projection and leakage:

\[ \dot{ {\hat{ {\boldsymbol{ {\theta} } } } } } = \operatorname{ {Proj} }_{ {\Omega} }\!\left[ -\boldsymbol{ {\Gamma} }\mathbf{ {Y} }^T\mathbf{ {s} }- \sigma\boldsymbol{ {\Gamma} }(\hat{ {\boldsymbol{ {\theta} } } }- \boldsymbol{ {\theta} }_0)\right], \qquad \sigma>0. \]

Projection keeps estimates in a physically defensible compact set \(\Omega\). Under the projection inequality, bounded \(\|\mathbf{ {w} }+\Delta\boldsymbol{ {\tau} }\|\le\bar w\), and constant true parameters, the derivative satisfies

\[ \dot V \le -\lambda_{ {\min} }(\mathbf{ {K} }_s)\|\mathbf{ {s} }\|^2 +\|\mathbf{ {s} }\|\bar w -\sigma\|\tilde{ {\boldsymbol{ {\theta} } } }\|^2 -\sigma\tilde{ {\boldsymbol{ {\theta} } } }^T (\boldsymbol{ {\theta} }-\boldsymbol{ {\theta} }_0). \]

Apply Young's inequality with \(\eta>0\):

\[ \|\mathbf{ {s} }\|\bar w \le \frac{ {\eta} }{ {2} }\|\mathbf{ {s} }\|^2+ \frac{ {\bar w^2} }{ {2\eta} }, \qquad -\sigma\tilde{ {\boldsymbol{ {\theta} } } }^T (\boldsymbol{ {\theta} }-\boldsymbol{ {\theta} }_0) \le \frac{ {\sigma} }{ {2} }\|\tilde{ {\boldsymbol{ {\theta} } } }\|^2+ \frac{ {\sigma} }{ {2} }\|\boldsymbol{ {\theta} }- \boldsymbol{ {\theta} }_0\|^2. \]

Therefore, if \(\eta<2\lambda_{ {\min} }(\mathbf{ {K} }_s)\),

\[ \dot V \le -\alpha_s\|\mathbf{ {s} }\|^2- \frac{ {\sigma} }{ {2} }\|\tilde{ {\boldsymbol{ {\theta} } } }\|^2+c, \quad \alpha_s=\lambda_{ {\min} }(\mathbf{ {K} }_s)-\frac{ {\eta} }{ {2} }, \quad c=\frac{ {\bar w^2} }{ {2\eta} }+ \frac{ {\sigma} }{ {2} }\|\boldsymbol{ {\theta} }- \boldsymbol{ {\theta} }_0\|^2. \]

The closed loop is uniformly ultimately bounded. This is the correct practical interpretation: tighter tracking requires lower implementation error, greater damping, and carefully selected leakage—not unlimited adaptation gain.

4. Sampled Implementation: What Changes

With sample period \(T_s\), a direct Euler update is

\[ \hat{ {\boldsymbol{ {\theta} } } }_{ {k+1} }= \Pi_{ {\Omega} }\!\left(\hat{ {\boldsymbol{ {\theta} } } }_k + T_s \mathbf{ {g} }_k\right), \quad \mathbf{ {g} }_k=-\boldsymbol{ {\Gamma} }\mathbf{ {Y} }_k^T\mathbf{ {s} }_k- \sigma\boldsymbol{ {\Gamma} }(\hat{ {\boldsymbol{ {\theta} } } }_k- \boldsymbol{ {\theta} }_0). \]

The effective parameter step is proportional to \(T_s\|\boldsymbol{ {\Gamma} }\|\|\mathbf{ {Y} }_k\| \|\mathbf{ {s} }_k\|\). Large regressors, noise spikes, or missed deadlines can therefore produce abrupt estimate changes even when the continuous-time law appears benign. A practical rate bound is

\[ \|\hat{ {\boldsymbol{ {\theta} } } }_{ {k+1} }- \hat{ {\boldsymbol{ {\theta} } } }_k\| \le \Delta_{ {\theta,\max} }, \]

enforced by projection, explicit derivative clipping, or both. The applied torque should also satisfy amplitude and slew constraints:

\[ \boldsymbol{ {\tau} }_k= \operatorname{ {RateLimit} }\!\left( \operatorname{ {Sat} }_{ {\tau_{ {\max} } } }(\boldsymbol{ {\tau} }_{ {c,k} }), \boldsymbol{ {\tau} }_{ {k-1} },\dot{ {\boldsymbol{ {\tau} } } }_{ {\max} },T_s \right). \]

Adaptation should be frozen or modified while the commanded torque is not physically realized; otherwise the estimator attempts to explain actuator saturation as a change in robot parameters.

5. Commissioning Workflow from Simulation to Production

flowchart TD
  A["Offline model and parameter bounds"] --> B["Software-in-the-loop tests"]
  B --> C["Hardware-in-the-loop with real sample timing"]
  C --> D["Low-energy single-axis commissioning"]
  D --> E["Adaptation disabled: verify baseline controller"]
  E --> F["Enable projection, leakage, dead zone, logging"]
  F --> G["Inject bounded payload and disturbance tests"]
  G --> H["Review metrics, saturation, timing, safety trips"]
  H --> I{"Acceptance criteria \nsatisfied?"}
  I -->|No| J["Retune, revise bounds, or \nreduce operating envelope"]
  J --> B
  I -->|Yes| K["Multi-axis and task-level validation"]
  K --> L["Release with configuration control and rollback"]
        

The baseline fixed-gain controller must be demonstrated before adaptation is enabled. Adaptation is an augmentation, not a substitute for correct sign conventions, encoder calibration, torque scaling, gravity direction, emergency-stop logic, or a stable nominal loop.

6. Failure Modes and Their Mathematical Signatures

Noise-driven drift. If velocity is reconstructed from position,

\[ \hat{ {\dot q} }_k = \alpha\hat{ {\dot q} }_{ {k-1} }+ (1-\alpha)\frac{ {q_k-q_{ {k-1} } } }{ {T_s} }, \quad 0<\alpha<1, \]

then differentiation amplifies encoder noise. Since the update contains \(\mathbf{ {Y} }^T\mathbf{ {s} }\), noisy \(\mathbf{ {s} }\) can sustain parameter motion near zero tracking error. Dead-zone adaptation sets the gradient term to zero when \(\|\mathbf{ {s} }\|\le s_0\).

Insufficient excitation. Parameter convergence requires informative motion:

\[ \exists\,T,\mu>0:\quad \int_t^{ {t+T} }\mathbf{ {Y} }^T(\tau)\mathbf{ {Y} }(\tau)d\tau \succeq \mu\mathbf{ {I} }. \]

A production task may track well without satisfying this condition. In that case, parameter estimates are internal controller states, not certified physical measurements.

Saturation-induced false learning. When \(\Delta\boldsymbol{ {\tau} }\) is large and correlated with \(\mathbf{ {s} }\), the adaptive law attributes the missing actuator effort to parameter error. Gating, command governors, and anti-windup-like compensation are required.

Unmodeled flexible dynamics. Increasing adaptation bandwidth can inject energy near structural resonances. The controller bandwidth, velocity filter, and adaptation rate must remain below validated flexible modes unless those modes are explicitly modeled.

Estimator reset transients. Resetting estimates changes feedforward torque immediately. A reset must be bumpless, rate-limited, and performed only in a safe state.

7. Safety Architecture and Adaptation Authority

flowchart TD
  R["Trajectory reference"] --> C["Nominal plus \nadaptive controller"]
  S["Sensors and state estimates"] --> C
  C --> L["Torque amplitude \nand rate limits"]
  L --> A["Actuator command"]
  A --> P["Robot plant"]
  P --> S
  S --> M["Independent safety monitor"]
  C --> M
  M --> G{"Safe state?"}
  G -->|Yes| E["Enable bounded adaptation"]
  E --> C
  G -->|No| F["Freeze estimates and \ncommand safe torque"]
  F --> L
  M --> D["Event log and fault code"]
        

The safety monitor must not depend solely on the adaptive model. It should independently check joint position, velocity, following error, torque, temperature, communication health, numerical finiteness, parameter bounds, and execution timing. Adaptation authority can be represented by a gate \(\rho_k\in\{0,1\}\):

\[ \hat{ {\boldsymbol{ {\theta} } } }_{ {k+1} }= \Pi_{ {\Omega} }\!\left(\hat{ {\boldsymbol{ {\theta} } } }_k+ \rho_k T_s\mathbf{ {g} }_k\right), \]

where \(\rho_k=0\) during saturation, safety trips, invalid sensing, startup, shutdown, manual jogging, or any mode not covered by the validation evidence.

8. Quantitative Acceptance Metrics

Tracking performance should be reported with reproducible metrics:

\[ e_{ {\mathrm{ {RMS} } } }= \sqrt{ {\frac{ {1} }{ {T_2-T_1} }\int_{ {T_1} }^{ {T_2} }\|\mathbf{ {e} }(t)\|^2dt} }, \qquad e_{ {\max} }=\sup_{ {t\in[T_1,T_2]} }\|\mathbf{ {e} }(t)\|. \]

Control and adaptation burden can be measured by

\[ J_\tau=\int_{ {T_1} }^{ {T_2} }\|\boldsymbol{ {\tau} }(t)\|^2dt, \qquad J_{ {\dot\theta} }=\int_{ {T_1} }^{ {T_2} } \|\dot{ {\hat{ {\boldsymbol{ {\theta} } } } } }(t)\|^2dt. \]

Additional release criteria should include saturation fraction, maximum torque slew, missed-deadline count, safety-trip count, parameter-bound contact time, recovery time after payload change, and performance with adaptation disabled. A useful comparison is the normalized improvement

\[ I_e = \frac{ {e_{ {\mathrm{ {RMS} },\,fixed} }- e_{ {\mathrm{ {RMS} },\,adaptive} } } }{ {e_{ {\mathrm{ {RMS} },\,fixed} } } }, \]

reported together with any increase in control effort or safety-margin consumption. An adaptive controller is not better merely because its parameter estimates move.

9. Common Benchmark Used in All Implementations

To keep five language implementations directly comparable, we use one revolute joint:

\[ J\ddot q+b\dot q+g\sin q=\tau+d(t), \qquad Y=[\ddot q_r,\dot q_r,\sin q], \qquad \theta=[J,b,g]^T. \]

At \(t=6\,\mathrm{ {s} }\), the inertia, viscous coefficient, and gravity coefficient change to emulate tool pickup. The implementation includes noisy position sensing, filtered velocity reconstruction, torque amplitude and slew limits, dead-zone adaptation, leakage, elementwise projection, saturation gating, a disturbance pulse, and a safety trip.

10. Python Implementation

Chapter25_Lesson5.py

"""Chapter25_Lesson5.py
Industrial-style adaptive control benchmark for one robot joint.
Dependencies: numpy, matplotlib.
"""
from __future__ import annotations

import math
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np


@dataclass
class Config:
    dt: float = 0.001
    duration: float = 12.0
    lam: float = 5.0
    k_s: float = 10.0
    gamma: tuple[float, float, float] = (4.0, 1.5, 2.0)
    sigma: float = 0.04
    dead_zone: float = 0.02
    torque_limit: float = 18.0
    torque_rate_limit: float = 1500.0
    velocity_cutoff_hz: float = 15.0
    q_limit: float = 2.5
    qdot_limit: float = 6.0


def desired_motion(t: float) -> tuple[float, float, float]:
    q_d = 0.6 * math.sin(0.8 * t) + 0.2 * math.sin(1.7 * t)
    qd_d = 0.48 * math.cos(0.8 * t) + 0.34 * math.cos(1.7 * t)
    qdd_d = -0.384 * math.sin(0.8 * t) - 0.578 * math.sin(1.7 * t)
    return q_d, qd_d, qdd_d


def project_rate(theta: np.ndarray, rate: np.ndarray,
                 lower: np.ndarray, upper: np.ndarray) -> np.ndarray:
    projected = rate.copy()
    for i in range(theta.size):
        if theta[i] <= lower[i] and projected[i] < 0.0:
            projected[i] = 0.0
        if theta[i] >= upper[i] and projected[i] > 0.0:
            projected[i] = 0.0
    return projected


def clamp_rate(command: float, previous: float, limit: float, dt: float) -> float:
    step = limit * dt
    return float(np.clip(command, previous - step, previous + step))


def run_simulation(cfg: Config) -> dict[str, np.ndarray | float]:
    steps = int(cfg.duration / cfg.dt) + 1
    t_log = np.linspace(0.0, cfg.duration, steps)

    q = 0.0
    qdot = 0.0
    q_meas_prev = 0.0
    qdot_est = 0.0
    tau_prev = 0.0
    theta_hat = np.array([1.4, 0.20, 6.5], dtype=float)
    theta_nominal = theta_hat.copy()
    theta_lower = np.array([0.6, 0.02, 2.0])
    theta_upper = np.array([4.0, 1.20, 14.0])
    gamma = np.array(cfg.gamma)
    alpha = math.exp(-2.0 * math.pi * cfg.velocity_cutoff_hz * cfg.dt)
    rng = np.random.default_rng(7)

    q_log = np.zeros(steps)
    qd_log = np.zeros(steps)
    e_log = np.zeros(steps)
    tau_log = np.zeros(steps)
    theta_log = np.zeros((steps, 3))
    gate_log = np.zeros(steps)
    saturation_log = np.zeros(steps)
    safety_log = np.zeros(steps)

    safety_trip = False

    for k, t in enumerate(t_log):
        # Payload change at t = 6 s emulates a tool pickup.
        if t < 6.0:
            inertia, viscous, gravity = 1.8, 0.30, 8.0
        else:
            inertia, viscous, gravity = 2.4, 0.38, 10.0

        disturbance = 0.30 * math.sin(4.0 * t)
        if 8.0 <= t <= 8.2:
            disturbance += 0.8

        q_d, qdot_d, qddot_d = desired_motion(t)
        q_meas = q + rng.normal(0.0, 0.0003)
        raw_derivative = (q_meas - q_meas_prev) / cfg.dt
        qdot_est = alpha * qdot_est + (1.0 - alpha) * raw_derivative
        q_meas_prev = q_meas

        e = q_meas - q_d
        edot = qdot_est - qdot_d
        qdot_r = qdot_d - cfg.lam * e
        qddot_r = qddot_d - cfg.lam * edot
        s = qdot_est - qdot_r
        regressor = np.array([qddot_r, qdot_r, math.sin(q_meas)])

        tau_raw = float(regressor @ theta_hat - cfg.k_s * s)
        tau_limited = float(np.clip(tau_raw, -cfg.torque_limit, cfg.torque_limit))
        tau = clamp_rate(tau_limited, tau_prev, cfg.torque_rate_limit, cfg.dt)
        saturated = abs(tau_raw - tau) > 1.0e-9

        finite = np.isfinite([q_meas, qdot_est, tau, *theta_hat]).all()
        if (not finite or abs(q_meas) > cfg.q_limit or
                abs(qdot_est) > cfg.qdot_limit):
            safety_trip = True

        adaptation_gate = (
            (not safety_trip)
            and (not saturated)
            and abs(s) > cfg.dead_zone
            and abs(tau_raw) < 0.98 * cfg.torque_limit
        )

        if adaptation_gate:
            theta_rate = -gamma * regressor * s - cfg.sigma * gamma * (theta_hat - theta_nominal)
            theta_rate = project_rate(theta_hat, theta_rate, theta_lower, theta_upper)
            theta_hat += cfg.dt * theta_rate
            theta_hat = np.clip(theta_hat, theta_lower, theta_upper)

        if safety_trip:
            tau = 0.0

        qddot = (tau + disturbance - viscous * qdot - gravity * math.sin(q)) / inertia
        qdot += cfg.dt * qddot
        q += cfg.dt * qdot
        tau_prev = tau

        q_log[k] = q
        qd_log[k] = q_d
        e_log[k] = q - q_d
        tau_log[k] = tau
        theta_log[k, :] = theta_hat
        gate_log[k] = float(adaptation_gate)
        saturation_log[k] = float(saturated)
        safety_log[k] = float(safety_trip)

    evaluation = t_log >= 2.0
    rmse = float(np.sqrt(np.mean(e_log[evaluation] ** 2)))
    saturation_fraction = float(np.mean(saturation_log))

    return {
        "t": t_log,
        "q": q_log,
        "qd": qd_log,
        "e": e_log,
        "tau": tau_log,
        "theta": theta_log,
        "gate": gate_log,
        "saturation": saturation_log,
        "safety": safety_log,
        "rmse": rmse,
        "saturation_fraction": saturation_fraction,
    }


def main() -> None:
    result = run_simulation(Config())
    print(f"Tracking RMSE after 2 s: {result['rmse']:.6f} rad")
    print(f"Torque saturation fraction: {100.0 * result['saturation_fraction']:.2f}%")
    print("Final parameter estimate [J, b, g]:", result["theta"][-1])
    print("Safety trip occurred:", bool(np.any(result["safety"])))

    plt.figure()
    plt.plot(result["t"], result["qd"], label="desired")
    plt.plot(result["t"], result["q"], label="actual")
    plt.xlabel("Time [s]")
    plt.ylabel("Position [rad]")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(result["t"], result["e"])
    plt.xlabel("Time [s]")
    plt.ylabel("Tracking error [rad]")
    plt.grid(True)

    plt.figure()
    plt.plot(result["t"], result["tau"])
    plt.xlabel("Time [s]")
    plt.ylabel("Torque [N m]")
    plt.grid(True)

    plt.figure()
    plt.plot(result["t"], result["theta"][:, 0], label="J estimate")
    plt.plot(result["t"], result["theta"][:, 1], label="b estimate")
    plt.plot(result["t"], result["theta"][:, 2], label="g estimate")
    plt.xlabel("Time [s]")
    plt.ylabel("Parameter estimate")
    plt.legend()
    plt.grid(True)
    plt.show()


if __name__ == "__main__":
    main()

11. C++ Implementation

Chapter25_Lesson5.cpp

// Chapter25_Lesson5.cpp
// Industrial-style adaptive control benchmark for one robot joint.
// Build: g++ -std=c++17 -O2 Chapter25_Lesson5.cpp -o Chapter25_Lesson5
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>

struct DesiredMotion {
    double q;
    double qdot;
    double qddot;
};

DesiredMotion desiredMotion(double t) {
    return {
        0.6 * std::sin(0.8 * t) + 0.2 * std::sin(1.7 * t),
        0.48 * std::cos(0.8 * t) + 0.34 * std::cos(1.7 * t),
        -0.384 * std::sin(0.8 * t) - 0.578 * std::sin(1.7 * t)
    };
}

double clampValue(double x, double lo, double hi) {
    return std::max(lo, std::min(x, hi));
}

double clampRate(double command, double previous, double limit, double dt) {
    const double step = limit * dt;
    return clampValue(command, previous - step, previous + step);
}

int main() {
    constexpr double dt = 0.001;
    constexpr double duration = 12.0;
    constexpr double lambda = 5.0;
    constexpr double ks = 10.0;
    constexpr double sigma = 0.04;
    constexpr double deadZone = 0.02;
    constexpr double torqueLimit = 18.0;
    constexpr double torqueRateLimit = 1500.0;
    constexpr double velocityCutoffHz = 15.0;
    constexpr double qLimit = 2.5;
    constexpr double qdotLimit = 6.0;
    constexpr int steps = static_cast<int>(duration / dt) + 1;

    std::array<double, 3> theta{1.4, 0.20, 6.5};
    const std::array<double, 3> thetaNominal = theta;
    const std::array<double, 3> lower{0.6, 0.02, 2.0};
    const std::array<double, 3> upper{4.0, 1.20, 14.0};
    const std::array<double, 3> gamma{4.0, 1.5, 2.0};

    double q = 0.0;
    double qdot = 0.0;
    double qMeasuredPrevious = 0.0;
    double qdotEstimate = 0.0;
    double previousTorque = 0.0;
    bool safetyTrip = false;
    const double alpha = std::exp(-2.0 * 3.14159265358979323846 * velocityCutoffHz * dt);

    std::mt19937_64 generator(7);
    std::normal_distribution<double> noise(0.0, 0.0003);
    std::ofstream csv("Chapter25_Lesson5_cpp_results.csv");
    if (!csv) {
        std::cerr << "Cannot open output CSV file.\n";
        return 1;
    }
    csv << "t,q,qd,error,tau,Jhat,bhat,ghat,gate,saturated,safety\n";

    double squaredErrorSum = 0.0;
    int evaluationCount = 0;
    int saturationCount = 0;

    for (int k = 0; k < steps; ++k) {
        const double t = k * dt;
        const double inertia = (t < 6.0) ? 1.8 : 2.4;
        const double viscous = (t < 6.0) ? 0.30 : 0.38;
        const double gravity = (t < 6.0) ? 8.0 : 10.0;
        double disturbance = 0.30 * std::sin(4.0 * t);
        if (t >= 8.0 && t <= 8.2) {
            disturbance += 0.8;
        }

        const DesiredMotion desired = desiredMotion(t);
        const double qMeasured = q + noise(generator);
        const double rawDerivative = (qMeasured - qMeasuredPrevious) / dt;
        qdotEstimate = alpha * qdotEstimate + (1.0 - alpha) * rawDerivative;
        qMeasuredPrevious = qMeasured;

        const double error = qMeasured - desired.q;
        const double errorDot = qdotEstimate - desired.qdot;
        const double qdotReference = desired.qdot - lambda * error;
        const double qddotReference = desired.qddot - lambda * errorDot;
        const double s = qdotEstimate - qdotReference;
        const std::array<double, 3> y{qddotReference, qdotReference, std::sin(qMeasured)};

        const double torqueRaw = y[0] * theta[0] + y[1] * theta[1] + y[2] * theta[2] - ks * s;
        const double torqueAmplitudeLimited = clampValue(torqueRaw, -torqueLimit, torqueLimit);
        double torque = clampRate(torqueAmplitudeLimited, previousTorque, torqueRateLimit, dt);
        const bool saturated = std::abs(torqueRaw - torque) > 1.0e-9;
        saturationCount += saturated ? 1 : 0;

        const bool finite = std::isfinite(qMeasured) && std::isfinite(qdotEstimate) &&
                            std::isfinite(torque) && std::isfinite(theta[0]) &&
                            std::isfinite(theta[1]) && std::isfinite(theta[2]);
        if (!finite || std::abs(qMeasured) > qLimit || std::abs(qdotEstimate) > qdotLimit) {
            safetyTrip = true;
        }

        const bool gate = !safetyTrip && !saturated && std::abs(s) > deadZone &&
                          std::abs(torqueRaw) < 0.98 * torqueLimit;
        if (gate) {
            for (std::size_t i = 0; i < theta.size(); ++i) {
                double rate = -gamma[i] * y[i] * s - sigma * gamma[i] * (theta[i] - thetaNominal[i]);
                if ((theta[i] <= lower[i] && rate < 0.0) ||
                    (theta[i] >= upper[i] && rate > 0.0)) {
                    rate = 0.0;
                }
                theta[i] = clampValue(theta[i] + dt * rate, lower[i], upper[i]);
            }
        }

        if (safetyTrip) {
            torque = 0.0;
        }

        const double qddot = (torque + disturbance - viscous * qdot - gravity * std::sin(q)) / inertia;
        qdot += dt * qddot;
        q += dt * qdot;
        previousTorque = torque;

        const double trueError = q - desired.q;
        if (t >= 2.0) {
            squaredErrorSum += trueError * trueError;
            ++evaluationCount;
        }

        csv << std::setprecision(12) << t << ',' << q << ',' << desired.q << ','
            << trueError << ',' << torque << ',' << theta[0] << ',' << theta[1] << ','
            << theta[2] << ',' << (gate ? 1 : 0) << ',' << (saturated ? 1 : 0) << ','
            << (safetyTrip ? 1 : 0) << '\n';
    }

    const double rmse = std::sqrt(squaredErrorSum / std::max(1, evaluationCount));
    std::cout << std::fixed << std::setprecision(6)
              << "Tracking RMSE after 2 s: " << rmse << " rad\n"
              << "Torque saturation fraction: "
              << 100.0 * static_cast<double>(saturationCount) / steps << "%\n"
              << "Final parameter estimate [J, b, g]: [" << theta[0] << ", "
              << theta[1] << ", " << theta[2] << "]\n"
              << "Safety trip occurred: " << (safetyTrip ? "true" : "false") << '\n'
              << "Results written to Chapter25_Lesson5_cpp_results.csv\n";
    return 0;
}

12. Java Implementation

Chapter25_Lesson5.java

// Chapter25_Lesson5.java
// Industrial-style adaptive control benchmark for one robot joint.
// Build and run: javac Chapter25_Lesson5.java && java Chapter25_Lesson5
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;
import java.util.Random;

public final class Chapter25_Lesson5 {
    private static final double DT = 0.001;
    private static final double DURATION = 12.0;
    private static final double LAMBDA = 5.0;
    private static final double KS = 10.0;
    private static final double SIGMA = 0.04;
    private static final double DEAD_ZONE = 0.02;
    private static final double TORQUE_LIMIT = 18.0;
    private static final double TORQUE_RATE_LIMIT = 1500.0;
    private static final double VELOCITY_CUTOFF_HZ = 15.0;
    private static final double Q_LIMIT = 2.5;
    private static final double QDOT_LIMIT = 6.0;

    private Chapter25_Lesson5() {}

    private static double clamp(double x, double lower, double upper) {
        return Math.max(lower, Math.min(x, upper));
    }

    private static double clampRate(double command, double previous, double limit, double dt) {
        double step = limit * dt;
        return clamp(command, previous - step, previous + step);
    }

    private static double[] desiredMotion(double t) {
        return new double[] {
            0.6 * Math.sin(0.8 * t) + 0.2 * Math.sin(1.7 * t),
            0.48 * Math.cos(0.8 * t) + 0.34 * Math.cos(1.7 * t),
            -0.384 * Math.sin(0.8 * t) - 0.578 * Math.sin(1.7 * t)
        };
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        int steps = (int) (DURATION / DT) + 1;
        double[] theta = {1.4, 0.20, 6.5};
        double[] thetaNominal = theta.clone();
        double[] lower = {0.6, 0.02, 2.0};
        double[] upper = {4.0, 1.20, 14.0};
        double[] gamma = {4.0, 1.5, 2.0};

        double q = 0.0;
        double qdot = 0.0;
        double qMeasuredPrevious = 0.0;
        double qdotEstimate = 0.0;
        double previousTorque = 0.0;
        boolean safetyTrip = false;
        double alpha = Math.exp(-2.0 * Math.PI * VELOCITY_CUTOFF_HZ * DT);
        Random random = new Random(7L);

        double squaredErrorSum = 0.0;
        int evaluationCount = 0;
        int saturationCount = 0;

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter25_Lesson5_java_results.csv"))) {
            writer.write("t,q,qd,error,tau,Jhat,bhat,ghat,gate,saturated,safety\n");

            for (int k = 0; k < steps; ++k) {
                double t = k * DT;
                double inertia = t < 6.0 ? 1.8 : 2.4;
                double viscous = t < 6.0 ? 0.30 : 0.38;
                double gravity = t < 6.0 ? 8.0 : 10.0;
                double disturbance = 0.30 * Math.sin(4.0 * t);
                if (t >= 8.0 && t <= 8.2) {
                    disturbance += 0.8;
                }

                double[] desired = desiredMotion(t);
                double qMeasured = q + 0.0003 * random.nextGaussian();
                double rawDerivative = (qMeasured - qMeasuredPrevious) / DT;
                qdotEstimate = alpha * qdotEstimate + (1.0 - alpha) * rawDerivative;
                qMeasuredPrevious = qMeasured;

                double error = qMeasured - desired[0];
                double errorDot = qdotEstimate - desired[1];
                double qdotReference = desired[1] - LAMBDA * error;
                double qddotReference = desired[2] - LAMBDA * errorDot;
                double s = qdotEstimate - qdotReference;
                double[] y = {qddotReference, qdotReference, Math.sin(qMeasured)};

                double torqueRaw = y[0] * theta[0] + y[1] * theta[1]
                        + y[2] * theta[2] - KS * s;
                double torqueAmplitudeLimited = clamp(torqueRaw, -TORQUE_LIMIT, TORQUE_LIMIT);
                double torque = clampRate(torqueAmplitudeLimited, previousTorque,
                        TORQUE_RATE_LIMIT, DT);
                boolean saturated = Math.abs(torqueRaw - torque) > 1.0e-9;
                if (saturated) {
                    ++saturationCount;
                }

                boolean finite = Double.isFinite(qMeasured) && Double.isFinite(qdotEstimate)
                        && Double.isFinite(torque) && Double.isFinite(theta[0])
                        && Double.isFinite(theta[1]) && Double.isFinite(theta[2]);
                if (!finite || Math.abs(qMeasured) > Q_LIMIT ||
                        Math.abs(qdotEstimate) > QDOT_LIMIT) {
                    safetyTrip = true;
                }

                boolean gate = !safetyTrip && !saturated && Math.abs(s) > DEAD_ZONE
                        && Math.abs(torqueRaw) < 0.98 * TORQUE_LIMIT;
                if (gate) {
                    for (int i = 0; i < theta.length; ++i) {
                        double rate = -gamma[i] * y[i] * s
                                - SIGMA * gamma[i] * (theta[i] - thetaNominal[i]);
                        if ((theta[i] <= lower[i] && rate < 0.0)
                                || (theta[i] >= upper[i] && rate > 0.0)) {
                            rate = 0.0;
                        }
                        theta[i] = clamp(theta[i] + DT * rate, lower[i], upper[i]);
                    }
                }

                if (safetyTrip) {
                    torque = 0.0;
                }

                double qddot = (torque + disturbance - viscous * qdot
                        - gravity * Math.sin(q)) / inertia;
                qdot += DT * qddot;
                q += DT * qdot;
                previousTorque = torque;

                double trueError = q - desired[0];
                if (t >= 2.0) {
                    squaredErrorSum += trueError * trueError;
                    ++evaluationCount;
                }

                writer.write(String.format(Locale.US,
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%d,%d,%d%n",
                        t, q, desired[0], trueError, torque, theta[0], theta[1], theta[2],
                        gate ? 1 : 0, saturated ? 1 : 0, safetyTrip ? 1 : 0));
            }
        }

        double rmse = Math.sqrt(squaredErrorSum / Math.max(1, evaluationCount));
        System.out.printf("Tracking RMSE after 2 s: %.6f rad%n", rmse);
        System.out.printf("Torque saturation fraction: %.2f%%%n",
                100.0 * saturationCount / steps);
        System.out.printf("Final parameter estimate [J, b, g]: [%.6f, %.6f, %.6f]%n",
                theta[0], theta[1], theta[2]);
        System.out.println("Safety trip occurred: " + safetyTrip);
        System.out.println("Results written to Chapter25_Lesson5_java_results.csv");
    }
}

13. MATLAB and Simulink Implementation

The first script runs the numerical benchmark using base MATLAB. The second script creates a Simulink safety shell showing the required ordering of saturation, torque-rate limiting, and safety override blocks.

Chapter25_Lesson5.m

% Chapter25_Lesson5.m
% Industrial-style adaptive control benchmark for one robot joint.
% This script uses base MATLAB only. The same signal flow can be placed in
% Simulink using Discrete-Time Integrator, Saturation, Rate Limiter,
% MATLAB Function, and Assertion blocks.
clear; clc; close all;
rng(7);

dt = 0.001;
T = 12.0;
t = (0:dt:T)';
N = numel(t);
lambda = 5.0;
ks = 10.0;
Gamma = [4.0; 1.5; 2.0];
sigma = 0.04;
deadZone = 0.02;
torqueLimit = 18.0;
torqueRateLimit = 1500.0;
velocityCutoffHz = 15.0;
qLimit = 2.5;
qdotLimit = 6.0;
alpha = exp(-2*pi*velocityCutoffHz*dt);

thetaHat = [1.4; 0.20; 6.5];
thetaNominal = thetaHat;
thetaLower = [0.6; 0.02; 2.0];
thetaUpper = [4.0; 1.20; 14.0];

q = 0.0;
qdot = 0.0;
qMeasuredPrevious = 0.0;
qdotEstimate = 0.0;
previousTorque = 0.0;
safetyTrip = false;

qLog = zeros(N,1);
qdLog = zeros(N,1);
errorLog = zeros(N,1);
torqueLog = zeros(N,1);
thetaLog = zeros(N,3);
gateLog = false(N,1);
saturationLog = false(N,1);
safetyLog = false(N,1);

for k = 1:N
    tk = t(k);
    if tk < 6.0
        inertia = 1.8; viscous = 0.30; gravity = 8.0;
    else
        inertia = 2.4; viscous = 0.38; gravity = 10.0;
    end

    disturbance = 0.30*sin(4.0*tk);
    if tk >= 8.0 && tk <= 8.2
        disturbance = disturbance + 0.8;
    end

    qDesired = 0.6*sin(0.8*tk) + 0.2*sin(1.7*tk);
    qdotDesired = 0.48*cos(0.8*tk) + 0.34*cos(1.7*tk);
    qddotDesired = -0.384*sin(0.8*tk) - 0.578*sin(1.7*tk);

    qMeasured = q + 0.0003*randn;
    rawDerivative = (qMeasured - qMeasuredPrevious)/dt;
    qdotEstimate = alpha*qdotEstimate + (1-alpha)*rawDerivative;
    qMeasuredPrevious = qMeasured;

    error = qMeasured - qDesired;
    errorDot = qdotEstimate - qdotDesired;
    qdotReference = qdotDesired - lambda*error;
    qddotReference = qddotDesired - lambda*errorDot;
    s = qdotEstimate - qdotReference;
    Y = [qddotReference, qdotReference, sin(qMeasured)];

    torqueRaw = Y*thetaHat - ks*s;
    torqueAmplitudeLimited = min(max(torqueRaw, -torqueLimit), torqueLimit);
    torqueStep = torqueRateLimit*dt;
    torque = min(max(torqueAmplitudeLimited, previousTorque-torqueStep), ...
        previousTorque+torqueStep);
    saturated = abs(torqueRaw - torque) > 1e-9;

    finiteSignals = all(isfinite([qMeasured, qdotEstimate, torque, thetaHat']));
    if ~finiteSignals || abs(qMeasured) > qLimit || abs(qdotEstimate) > qdotLimit
        safetyTrip = true;
    end

    gate = ~safetyTrip && ~saturated && abs(s) > deadZone && abs(torqueRaw) < 0.98*torqueLimit;
    if gate
        thetaRate = -Gamma.*Y'.*s - sigma*Gamma.*(thetaHat-thetaNominal);
        for i = 1:3
            if (thetaHat(i) <= thetaLower(i) && thetaRate(i) < 0) || ...
                    (thetaHat(i) >= thetaUpper(i) && thetaRate(i) > 0)
                thetaRate(i) = 0;
            end
        end
        thetaHat = thetaHat + dt*thetaRate;
        thetaHat = min(max(thetaHat, thetaLower), thetaUpper);
    end

    if safetyTrip
        torque = 0.0;
    end

    qddot = (torque + disturbance - viscous*qdot - gravity*sin(q))/inertia;
    qdot = qdot + dt*qddot;
    q = q + dt*qdot;
    previousTorque = torque;

    qLog(k) = q;
    qdLog(k) = qDesired;
    errorLog(k) = q - qDesired;
    torqueLog(k) = torque;
    thetaLog(k,:) = thetaHat';
    gateLog(k) = gate;
    saturationLog(k) = saturated;
    safetyLog(k) = safetyTrip;
end

evaluation = t >= 2.0;
rmse = sqrt(mean(errorLog(evaluation).^2));
fprintf('Tracking RMSE after 2 s: %.6f rad\n', rmse);
fprintf('Torque saturation fraction: %.2f%%\n', 100*mean(saturationLog));
fprintf('Final parameter estimate [J, b, g]: [%.6f %.6f %.6f]\n', thetaHat);
fprintf('Safety trip occurred: %d\n', any(safetyLog));

figure; plot(t, qdLog, '--', t, qLog, '-'); grid on;
xlabel('Time [s]'); ylabel('Position [rad]'); legend('desired','actual');
figure; plot(t, errorLog); grid on;
xlabel('Time [s]'); ylabel('Tracking error [rad]');
figure; plot(t, torqueLog); grid on;
xlabel('Time [s]'); ylabel('Torque [N m]');
figure; plot(t, thetaLog); grid on;
xlabel('Time [s]'); ylabel('Parameter estimate'); legend('J','b','g');

results = table(t, qLog, qdLog, errorLog, torqueLog, thetaLog(:,1), ...
    thetaLog(:,2), thetaLog(:,3), gateLog, saturationLog, safetyLog, ...
    'VariableNames', {'t','q','qd','error','tau','Jhat','bhat','ghat', ...
    'gate','saturated','safety'});
writetable(results, 'Chapter25_Lesson5_matlab_results.csv');

% Simulink implementation map:
% 1. Put desired trajectory and controller equations in a MATLAB Function block.
% 2. Use a Unit Delay for thetaHat and previousTorque.
% 3. Use Saturation and Rate Limiter blocks before the plant.
% 4. Use Discrete-Time Integrator blocks for qdot and q.
% 5. Use Assertion blocks for q, qdot, finite values, and parameter bounds.
% 6. Gate the update law with Switch blocks during saturation or safety trips.

Chapter25_Lesson5_Simulink.m

% Chapter25_Lesson5_Simulink.m
% Programmatically creates a compact Simulink safety shell for the lesson.
% It illustrates the industrial signal ordering around a controller subsystem.
% The detailed adaptive equations remain in Chapter25_Lesson5.m.
model = 'Chapter25_Lesson5_SafetyShell';
if bdIsLoaded(model)
    close_system(model, 0);
end
new_system(model);
open_system(model);

add_block('simulink/Sources/In1', [model '/RawTorque'], 'Position', [30 55 60 75]);
add_block('simulink/Discontinuities/Saturation', [model '/TorqueSaturation'], ...
    'UpperLimit', '18', 'LowerLimit', '-18', 'Position', [110 45 180 85]);
add_block('simulink/Discontinuities/Rate Limiter', [model '/TorqueRateLimiter'], ...
    'RisingSlewLimit', '1500', 'FallingSlewLimit', '-1500', ...
    'Position', [230 45 310 85]);
add_block('simulink/Signal Routing/Switch', [model '/SafetySwitch'], ...
    'Threshold', '0.5', 'Position', [380 40 430 90]);
add_block('simulink/Sources/Constant', [model '/ZeroTorque'], ...
    'Value', '0', 'Position', [315 115 345 135]);
add_block('simulink/Sources/In1', [model '/SafetyTrip'], ...
    'Port', '2', 'Position', [315 75 345 95]);
add_block('simulink/Sinks/Out1', [model '/AppliedTorque'], ...
    'Position', [500 55 530 75]);

add_line(model, 'RawTorque/1', 'TorqueSaturation/1');
add_line(model, 'TorqueSaturation/1', 'TorqueRateLimiter/1');
add_line(model, 'TorqueRateLimiter/1', 'SafetySwitch/1');
add_line(model, 'SafetyTrip/1', 'SafetySwitch/2');
add_line(model, 'ZeroTorque/1', 'SafetySwitch/3');
add_line(model, 'SafetySwitch/1', 'AppliedTorque/1');

set_param(model, 'StopTime', '12');
save_system(model, [model '.slx']);
fprintf('Created %s.slx\n', model);

14. Wolfram Mathematica Implementation

Chapter25_Lesson5.nb

Notebook[{
 Cell["Chapter 25, Lesson 5: Industrial Adaptive Robot Benchmark", "Title"],
 Cell["Chapter25_Lesson5.nb", "Subtitle"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"dt", "=", "0.001", ";", "T", "=", "12.0", ";", 
   "lambda", "=", "5.0", ";", "ks", "=", "10.0", ";", 
   "sigma", "=", "0.04", ";", "deadZone", "=", "0.02", ";", 
   "torqueLimit", "=", "18.0", ";", "torqueRateLimit", "=", "1500.0", ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{"desired", "[", "t_", "]", ":=", 
   RowBox[{"{", 
    RowBox[{"0.6 Sin[0.8 t]+0.2 Sin[1.7 t]", ",", 
     "0.48 Cos[0.8 t]+0.34 Cos[1.7 t]", ",", 
     "-0.384 Sin[0.8 t]-0.578 Sin[1.7 t]"}], "}"}]}]], "Input"],
 Cell[BoxData[
  RowBox[{"clip", "[", "x_,lo_,hi_", "]", ":=", "Min[Max[x,lo],hi]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"SeedRandom", "[", "7", "]", ";", 
   "theta", "=", "{1.4,0.20,6.5}", ";", 
   "theta0", "=", "theta", ";", 
   "lower", "=", "{0.6,0.02,2.0}", ";", 
   "upper", "=", "{4.0,1.20,14.0}", ";", 
   "gamma", "=", "{4.0,1.5,2.0}", ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{"q", "=", "0.0", ";", "qdot", "=", "0.0", ";", 
   "qPrev", "=", "0.0", ";", "qdotEst", "=", "0.0", ";", 
   "tauPrev", "=", "0.0", ";", "trip", "=", "False", ";", 
   "alpha", "=", "Exp[-2 Pi 15.0 dt]", ";", "log", "=", "{}", ";"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Do", "[", 
   RowBox[{"inertia=If[t<6.0,1.8,2.4];", 
    "viscous=If[t<6.0,0.30,0.38];", 
    "gravity=If[t<6.0,8.0,10.0];", 
    "dist=0.30 Sin[4.0 t]+If[8.0<=t<=8.2,0.8,0.0];", 
    "{qd,qdDot,qdDDot}=desired[t];", 
    "qMeas=q+RandomVariate[NormalDistribution[0,0.0003]];", 
    "raw=(qMeas-qPrev)/dt; qdotEst=alpha qdotEst+(1-alpha) raw; qPrev=qMeas;", 
    "e=qMeas-qd; eDot=qdotEst-qdDot;", 
    "qdotR=qdDot-lambda e; qddotR=qdDDot-lambda eDot;", 
    "s=qdotEst-qdotR; y={qddotR,qdotR,Sin[qMeas]};", 
    "tauRaw=y.theta-ks s;", 
    "tauAmp=clip[tauRaw,-torqueLimit,torqueLimit];", 
    "tau=clip[tauAmp,tauPrev-torqueRateLimit dt,tauPrev+torqueRateLimit dt];", 
    "sat=Abs[tauRaw-tau]>10^-9;", 
    "If[!And@@(NumericQ/@Join[{qMeas,qdotEst,tau},theta])||Abs[qMeas]>2.5||Abs[qdotEst]>6.0,trip=True];", 
    "gate=!trip&&!sat&&Abs[s]>deadZone&&Abs[tauRaw]<0.98 torqueLimit;", 
    "If[gate,rate=-gamma y s-sigma gamma (theta-theta0);", 
    "Do[If[(theta[[i]]<=lower[[i]]&&rate[[i]]<0)||(theta[[i]]>=upper[[i]]&&rate[[i]]>0),rate[[i]]=0],{i,3}];", 
    "theta=MapThread[clip,{theta+dt rate,lower,upper}]];", 
    "If[trip,tau=0.0];", 
    "qddot=(tau+dist-viscous qdot-gravity Sin[q])/inertia;", 
    "qdot=qdot+dt qddot; q=q+dt qdot; tauPrev=tau;", 
    "AppendTo[log,{t,q,qd,q-qd,tau,theta[[1]],theta[[2]],theta[[3]],Boole[gate],Boole[sat],Boole[trip]}]", 
    ",", "{t,0.0,T,dt}"}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"rmse", "=", "Sqrt[Mean[Select[log,#[[1]]>=2.0&][[All,4]]^2]]", ";", 
   "Print[\"Tracking RMSE after 2 s: \",N[rmse],\" rad\"]", ";", 
   "Print[\"Final parameter estimate: \",theta]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"ListLinePlot", "[", 
   RowBox[{"{log[[All,{1,2}]],log[[All,{1,3}]]}", ",", 
    "PlotLegends->{\"actual\",\"desired\"}", ",", "Frame->True"}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"Export", "[", 
   RowBox[{"\"Chapter25_Lesson5_mathematica_results.csv\"", ",", 
    "Prepend[log,{\"t\",\"q\",\"qd\",\"error\",\"tau\",\"Jhat\",\"bhat\",\"ghat\",\"gate\",\"saturated\",\"safety\"}]"}], "]"}]], "Input"]
}, WindowSize->{1100, 800}, StyleDefinitions->"Default.nb"]

15. Problems and Solutions

Problem 1 (Ultimate Bound): Suppose

\[ \dot V\le-a\|\mathbf{ {s} }\|^2-b\|\tilde{ {\boldsymbol{ {\theta} } } }\|^2+c, \qquad a,b,c>0. \]

Show that trajectories enter a compact residual set.

Solution: Outside the set

\[ \mathcal{ {B} }=\left\{ {(\mathbf{ {s} },\tilde{ {\boldsymbol{ {\theta} } } }): a\|\mathbf{ {s} }\|^2+b\|\tilde{ {\boldsymbol{ {\theta} } } }\|^2\le c}\right\}, \]

one has \(\dot V<0\). Since \(V\) is positive definite and radially unbounded on the projected parameter set, solutions approach and remain in \(\mathcal{ {B} }\). A conservative componentwise bound is \(\|\mathbf{ {s} }\|\le\sqrt{ {c/a} }\) and \(\|\tilde{ {\boldsymbol{ {\theta} } } }\|\le\sqrt{ {c/b} }\).

Problem 2 (Saturation Bias): For scalar filtered-error dynamics

\[ J\dot s=-k_s s+Y\tilde\theta+\Delta\tau+d, \]

explain why continued gradient adaptation during persistent negative \(\Delta\tau\) can produce false parameter learning.

Solution: With the ideal law, the cross term \(sY\tilde\theta\) is canceled by the parameter update, but \(s\Delta\tau\) remains. During positive tracking demand and upper saturation, \(\Delta\tau<0\). The resulting error is correlated with the regressor, so the gradient law changes \(\hat\theta\) as though plant parameters were wrong. Freezing adaptation while saturation is active prevents this attribution error.

Problem 3 (Dead-Zone Selection): Position noise has standard deviation \(\sigma_q\). Approximate the standard deviation of the unfiltered backward-difference velocity noise and state a dead-zone rule.

Solution: Independent consecutive samples give

\[ \operatorname{ {Std} }\!\left(\frac{ {n_k-n_{ {k-1} } } }{ {T_s} }\right) =\frac{ {\sqrt{ {2} }\sigma_q} }{ {T_s} }. \]

After filtering, measure the actual standard deviation \(\sigma_s\) of the filtered error while the robot is safely stationary. A common engineering rule is \(s_0\approx3\sigma_s\), followed by validation that real payload changes still exceed the threshold.

Problem 4 (Sampling Constraint): The maximum observed values are \(\|Y_k\|\le20\), \(|s_k|\le0.4\), and \(\|\Gamma\|=5\). Ignoring leakage, choose \(T_s\) so that one update changes the parameter vector by at most \(0.01\).

Solution:

\[ \|\Delta\hat\theta\|\le T_s\|\Gamma\|\|Y_k\||s_k| \le T_s(5)(20)(0.4)=40T_s. \]

Therefore \(40T_s\le0.01\), giving \(T_s\le2.5\times10^{ {-4} }\,\mathrm{ {s} }\). If the actual controller runs more slowly, reduce \(\Gamma\), normalize the update, or impose an explicit parameter-rate limit.

Problem 5 (Acceptance Decision): A fixed controller has \(e_{ {\mathrm{ {RMS} } } }=0.030\,\mathrm{ {rad} }\) and an adaptive controller has \(0.018\,\mathrm{ {rad} }\). The adaptive version increases \(J_\tau\) by 35% and spends 12% of the test in torque saturation. Compute \(I_e\) and interpret the result.

Solution:

\[ I_e=\frac{ {0.030-0.018} }{ {0.030} }=0.40. \]

Tracking RMS improves by 40%, but release is not automatic. The high saturation fraction means the experiment may be outside the validated actuator envelope, while the control-energy increase can affect heating and gearbox life. The task should be repeated with a command governor or lower bandwidth, and all safety and thermal limits must pass.

16. Summary

Successful robotic adaptive control depends on the complete closed-loop implementation, not only the continuous-time update law. The most reusable lessons are to prove a bounded-error result that includes implementation error; establish a stable nominal controller first; constrain estimates to physical ranges; filter measured signals; use leakage and dead zones; suspend adaptation during saturation and unsafe modes; enforce torque and parameter-rate limits; separate the safety monitor from the adaptive model; commission from low energy to full tasks; and release only against logged, quantitative acceptance criteria. Parameter convergence is desirable only when excitation supports it and should never be confused with the primary objective of safe tracking.

17. References

  1. Slotine, J.-J.E., & Li, W. (1987). On the adaptive control of robot manipulators. The International Journal of Robotics Research, 6(3), 49–59.
  2. Ortega, R., & Spong, M.W. (1989). Adaptive motion control of rigid robots: A tutorial. Automatica, 25(6), 877–888.
  3. Tomei, P. (1991). Adaptive PD controller for robot manipulators. IEEE Transactions on Robotics and Automation, 7(4), 565–570.
  4. Tao, G. (1992). On robust adaptive control of robot manipulators. Automatica, 28(4), 803–807.
  5. Spong, M.W. (1989). Adaptive control of flexible joint manipulators. Systems & Control Letters, 13(1), 15–21.
  6. Canudas de Wit, C., & Slotine, J.-J.E. (1991). Sliding observers for robot manipulators. Automatica, 27(5), 859–864.
  7. Reed, J.S., & Ioannou, P.A. (1989). Instability analysis and robust adaptive control of robotic manipulators. IEEE Transactions on Robotics and Automation, 5(3), 381–386.
  8. Panteley, E., Ortega, R., & Gäfvert, M. (1998). An adaptive friction compensator for global tracking in robot manipulators. Systems & Control Letters, 33(5), 307–313.
Support CaaT Academy

Help keep these engineering tutorials free and growing

If these lessons, examples, and project pages help you, a small donation supports the continued creation and improvement of free control, robotics, software, and engineering education resources.

Created and maintained by Abolfazl Mohammadijoo.