A conventional position loop may remain stable while failing its precision requirement. In a micrometer- or nanometer-class system, small changes in payload, cable force, lubricant state, thermal expansion, cogging force, current-loop gain, or structural resonance can dominate the allowed error budget. It is therefore useful to separate the error into deterministic, parameter-dependent, and unstructured components:

\[ e(t)=e_{\mathrm{par} }(t)+e_{\mathrm{per} }(t)+e_{\mathrm{unm} }(t) +e_{\mathrm{sens} }(t), \]

where \(e_{\mathrm{par} }\) is caused by uncertain physical parameters, \(e_{\mathrm{per} }\) by repeatable periodic disturbances such as force ripple, \(e_{\mathrm{unm} }\) by unmodeled dynamics and nonrepeatable disturbances, and \(e_{\mathrm{sens} }\) by sensing, quantization, and estimation. A useful scalar performance measure is the post-transient RMS error

\[ e_{\mathrm{RMS} }= \sqrt{\frac{1}{T_2-T_1}\int_{T_1}^{T_2}e^2(t)\,dt}, \]

but a complete evaluation also records peak error, settling time, control force, saturation duration, spectral error near ripple harmonics, and parameter-estimate boundedness. Adaptive control is most useful here when the dominant uncertainty can be expressed with known basis functions and unknown coefficients. It does not eliminate sensor noise, delays, or neglected high-frequency modes; those effects require filtering, bandwidth separation, robust modification, and safety logic.

Consider a single translational axis. The same construction can be applied joint-by-joint or embedded in the multivariable structures introduced in Chapter 23. A useful low-frequency model is

\[ m\ddot q+b\dot q+kq+f_c\tanh\!\left(\frac{\dot q}{v_s}\right) +a_s\sin\!\left(\frac{2\pi q}{p}\right) +a_c\cos\!\left(\frac{2\pi q}{p}\right)=u+d(t), \]

where \(q\) is position, \(u\) is actuator force, \(m\) is effective moving mass, \(b\) is viscous damping, \(k\) is the effective stiffness seen in the controlled bandwidth, and \(f_c\) is a smooth approximation of Coulomb friction. The pair \((a_s,a_c)\) represents the first spatial harmonic of cogging or force ripple with pitch \(p\). The remainder \(d(t)\) contains bounded unmodeled effects.

The unknown coefficients form the parameter vector

\[ \boldsymbol\theta= \begin{bmatrix}m&b&k&f_c&a_s&a_c\end{bmatrix}^{\!T}. \]

Although friction and ripple are nonlinear functions of the measured signals, the model is linear in the unknown parameters. This is the structural property required by the Lyapunov adaptive design. Physical information should be encoded as a compact admissible set

\[ \Omega=\left\{\boldsymbol\theta: \underline{\boldsymbol\theta}\preceq\boldsymbol\theta \preceq\overline{\boldsymbol\theta},\quad m\geq m_{\min}>0\right\}, \]

so that projection prevents nonphysical estimates such as negative mass or excessive friction compensation.

flowchart TD
  R["Reference qd, qd_dot, qd_ddot"] --> E["Filtered error: s = e_dot + lambda e"]
  Q["Encoder position"] --> V["Quantization-aware velocity filter"]
  V --> E
  R --> Y["Regressor construction"]
  Q --> Y
  V --> Y
  H["Parameter estimates"] --> FF["Adaptive model compensation"]
  Y --> FF
  E --> FB["Damping feedback: -ks s"]
  FF --> SUM["Control summation"]
  FB --> SUM
  SUM --> SAT["Force limit and adaptation gate"]
  SAT --> P["Precision stage"]
  P --> Q
  E --> A["Normalized projection update"]
  Y --> A
  A --> H
        

Let the position-tracking error be

\[ e=q-q_d. \]

Introduce a stable first-order error filter with \(\lambda>0\):

\[ \dot q_r=\dot q_d-\lambda e,\qquad s=\dot q-\dot q_r=\dot e+\lambda e. \]

The corresponding reference acceleration is

\[ \ddot q_r=\ddot q_d-\lambda\dot e. \]

The filtered error is chosen because convergence of \(s\) directly implies convergence of the original tracking error. Indeed,

\[ \dot e+\lambda e=s \quad\Longrightarrow\quad e(t)=e^{-\lambda t}e(0)+ \int_0^t e^{-\lambda(t-\tau)}s(\tau)\,d\tau. \]

Therefore, if \(s(t)\to0\) and is bounded, the stable filter guarantees \(e(t)\to0\). In implementation, \(\dot q\) is normally estimated from the encoder through a filtered differentiator rather than obtained by raw finite differences.

Define the precision-motion regressor

\[ \mathbf{Y}(q,\dot q,q_d)= \begin{bmatrix} \ddot q_r\\ \dot q_r\\ q\\ \tanh(\dot q/v_s)\\ \sin(2\pi q/p)\\ \cos(2\pi q/p) \end{bmatrix}. \]

The control law is

\[ u=\hat{\boldsymbol\theta}^{T}\mathbf{Y}-k_s s, \qquad k_s>0. \]

The first term is adaptive model compensation; the second injects damping into the filtered-error dynamics. Since \(\dot q=s+\dot q_r\) and \(\ddot q=\dot s+\ddot q_r\), substitution into the plant gives

\[ m\dot s=-(b+k_s)s+ \tilde{\boldsymbol\theta}^{T}\mathbf{Y}+d(t), \qquad \tilde{\boldsymbol\theta}=\hat{\boldsymbol\theta}-\boldsymbol\theta. \]

The ideal continuous-time update used for the nominal proof is

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

For a sampled precision system, a more practical update combines normalization, projection, leakage, and an anti-saturation gate:

\[ \dot{\hat{\boldsymbol\theta} }= \operatorname{Proj}_{\Omega}\!\left[ -\frac{\mathbf{\Gamma}\mathbf{Y}s} {1+\mathbf{Y}^{T}\mathbf{Y} } -\sigma\left(\hat{\boldsymbol\theta}-\boldsymbol\theta_0\right) \right], \qquad \sigma\geq0. \]

Here \(\boldsymbol\theta_0\) is a physically plausible nominal vector. The code freezes the update whenever the requested force exceeds the actuator limit, because saturated data no longer correspond to the control law assumed in the nominal error model.

First set \(d(t)=0\), omit leakage and normalization, and assume that the true parameter vector lies in the interior of the projection set. Choose

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

Because the true parameters are constant, \(\dot{\tilde{\boldsymbol\theta} }= \dot{\hat{\boldsymbol\theta} }\). Along the closed-loop system,

\[ \begin{aligned} \dot V &=ms\dot s+ \tilde{\boldsymbol\theta}^{T}\mathbf{\Gamma}^{-1} \dot{\hat{\boldsymbol\theta} }\\ &=s\left[-(b+k_s)s+ \tilde{\boldsymbol\theta}^{T}\mathbf{Y}\right] +\tilde{\boldsymbol\theta}^{T}\mathbf{\Gamma}^{-1} \left(-\mathbf{\Gamma}\mathbf{Y}s\right)\\ &=-(b+k_s)s^2\leq0. \end{aligned} \]

Consequently, \(s\) and \(\tilde{\boldsymbol\theta}\) are bounded, and \(s\in L_2\). If the reference and regressor are bounded, then \(\dot s\) is bounded. Barbalat's lemma yields \(s(t)\to0\), and the stable filtered-error equation gives \(e(t)\to0\).

This proof establishes tracking convergence, not automatic convergence of every physical parameter. Parameter convergence requires sufficient excitation of the corresponding regressor directions. A precision stage can track a narrow trajectory accurately while several estimates remain nonunique or close to their initial values.

Suppose \(|d(t)|\leq\bar d\). Ignoring the bounded projection correction for the moment, the disturbance contribution is \(s d\). For any \(\varepsilon>0\), Young's inequality gives

\[ |s d|\leq\frac{\varepsilon}{2}s^2+ \frac{\bar d^2}{2\varepsilon}. \]

Hence, without leakage,

\[ \dot V\leq- \left(b+k_s-\frac{\varepsilon}{2}\right)s^2+ \frac{\bar d^2}{2\varepsilon}. \]

Selecting \(0<\varepsilon<2(b+k_s)\) proves that the filtered error enters an ultimate neighborhood whose size decreases with stronger damping and increases with disturbance magnitude. Leakage adds a restoring term toward the nominal estimate. Using

\[ -\sigma\tilde{\boldsymbol\theta}^{T} \mathbf{\Gamma}^{-1} (\hat{\boldsymbol\theta}-\boldsymbol\theta_0) \leq-c_\theta\|\tilde{\boldsymbol\theta}\|^2+C_\theta, \]

for suitable positive constants \(c_\theta\) and \(C_\theta\), one obtains uniform ultimate boundedness of both tracking and parameter errors. Projection preserves the essential Lyapunov inequality because its correction does not point outward relative to an admissible true parameter vector. The practical design objective is therefore not unlimited adaptation gain; it is the smallest certified ultimate error compatible with noise, sampling, resonance margins, and actuator limits.

Direct-drive motors and screw-driven stages often exhibit disturbances that repeat with electrical angle, magnet pitch, screw pitch, or bearing rotation. A spatial Fourier basis extends the regressor without changing the linear-in-parameters structure:

\[ d_p(q)=\sum_{j=1}^{N_h} \left[a_{s,j}\sin\!\left(\frac{2\pi j q}{p}\right)+ a_{c,j}\cos\!\left(\frac{2\pi j q}{p}\right)\right] =\mathbf{w}^{*T}\boldsymbol\psi_p(q). \]

The adaptive feedforward term is \(\hat{\mathbf w}^{T}\boldsymbol\psi_p(q)\), with update

\[ \dot{\hat{\mathbf w} }=- \frac{\mathbf{\Gamma}_w\boldsymbol\psi_p(q)s} {1+\boldsymbol\psi_p^T(q)\boldsymbol\psi_p(q)}. \]

The basis must match the physical repeatability variable. A disturbance fixed to motor position should use position harmonics; one fixed to time or spindle speed may require phase generated from an encoder or clock. Adding many harmonics increases approximation capability but also increases excitation requirements and susceptibility to noise. The benchmark uses only the first sine/cosine pair so the contribution of adaptive ripple cancellation remains transparent.

A continuous-time proof must be translated into a deterministic sampled implementation. The following practices are especially important.

Encoder processing. Quantize the simulated measurement at the actual encoder resolution. Estimate velocity with a causal low-pass differentiator; raw differencing converts position quantization directly into high-frequency velocity noise.

Bandwidth separation. Keep the adaptation dynamics slower than neglected flexible modes and the inner current loop. The adaptive algorithm should modify low-frequency compensation, not attempt to chase every sample of sensor noise.

Normalization and projection. Normalize by regressor energy to reduce trajectory-dependent update magnitude, and enforce physical parameter bounds at every sample.

Saturation management. Apply the force limit before the plant and suspend or modify adaptation when saturation occurs. Otherwise, the estimator interprets actuator deficiency as a parameter error and can drift toward its bounds.

Monitoring. Log reference, measured position, filtered velocity, control force, saturation status, filtered error, parameter estimates, RMS/peak error, and spectral lines at known ripple frequencies.

flowchart TD
  A["Sample encoder"] --> B["Apply resolution and validity checks"]
  B --> C["Filtered velocity estimate"]
  C --> D["Compute e, e_dot, s"]
  D --> E["Build physical and periodic regressors"]
  E --> F["Compute unsaturated force"]
  F --> G["Apply force and rate limits"]
  G --> H["Send command to current loop"]
  G --> I["Saturated?"]
  I -->|"no"| J["Normalized projection update"]
  I -->|"yes"| K["Freeze or safely \nmodify adaptation"]
  J --> L["Log errors, estimates, \nand safety flags"]
  K --> L
  L --> A
        

The supplied programs simulate a two-second trajectory containing 2 Hz and 5 Hz position components. The true stage parameters and deliberately inaccurate nominal estimates are:

The simulation also includes 10 nm encoder quantization, deterministic nanometer-scale sensor noise, a 73 Hz bounded disturbance, a filtered velocity estimate, and a ±25 N force limit. With the common settings in the five implementations, the tested Python, C++, and Java programs produce:

Fixed:    RMS error = 21.877 um, peak error = 37.940 um, max |u| = 2.123 N
Adaptive: RMS error =  6.550 um, peak error = 22.414 um, max |u| = 2.124 N

This particular benchmark therefore reduces post-transient RMS error by approximately 70%. The result is not a universal performance guarantee; it is a reproducible illustration of how adaptation identifies the dominant friction/ripple compensation needed by the executed trajectory. Notice that all physical estimates do not converge to their true values. Accurate tracking can occur without full parameter identification because the trajectory may not persistently excite every regressor direction.

Chapter25_Lesson4.py

Requires NumPy. Matplotlib is optional; when available, the script saves a tracking-error plot in addition to the CSV file.

"""Chapter 25, Lesson 4: Adaptive Control in High-Precision Motion Systems.

Simulates a precision stage with uncertain mass, damping, stiffness, Coulomb
friction, position-periodic force ripple, encoder quantization, and measurement
noise. It compares a fixed nominal controller with a projection-based adaptive
controller using normalized adaptation and sigma modification.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path

import numpy as np


@dataclass
class Result:
    time: np.ndarray
    reference: np.ndarray
    position: np.ndarray
    error: np.ndarray
    control: np.ndarray
    parameters: np.ndarray


def reference_trajectory(t: float) -> tuple[float, float, float]:
    """Return desired position, velocity, and acceleration."""
    a1, f1 = 4.0e-4, 2.0
    a2, f2 = 1.5e-4, 5.0
    w1, w2 = 2.0 * np.pi * f1, 2.0 * np.pi * f2
    xd = a1 * np.sin(w1 * t) + a2 * np.sin(w2 * t)
    vd = a1 * w1 * np.cos(w1 * t) + a2 * w2 * np.cos(w2 * t)
    ad = -a1 * w1**2 * np.sin(w1 * t) - a2 * w2**2 * np.sin(w2 * t)
    return float(xd), float(vd), float(ad)


def simulate(adaptive: bool, dt: float = 1.0e-4, duration: float = 2.0) -> Result:
    """Run the closed-loop simulation."""
    steps = int(round(duration / dt))
    time = np.arange(steps, dtype=float) * dt

    # theta = [mass, viscous damping, stiffness, Coulomb friction,
    #          sine ripple coefficient, cosine ripple coefficient]
    theta_true = np.array([0.65, 18.0, 250.0, 0.55, 0.32, -0.20])
    theta_hat = np.array([0.45, 10.0, 180.0, 0.20, 0.0, 0.0])
    theta_nominal = theta_hat.copy()
    lower = np.array([0.20, 0.0, 50.0, 0.0, -1.0, -1.0])
    upper = np.array([1.20, 50.0, 500.0, 1.5, 1.0, 1.0])

    gamma = np.array([2.0e2, 5.0e3, 5.0e2, 2.0e3, 4.0e2, 4.0e2])
    sigma = 0.05
    lambda_e = 300.0
    k_s = 55.0
    velocity_scale = 2.0e-4
    ripple_pitch = 5.0e-3
    control_limit = 25.0
    encoder_resolution = 1.0e-8

    # First-order low-pass filter on the differentiated encoder signal.
    alpha_v = float(np.exp(-2.0 * np.pi * 400.0 * dt))

    position = np.zeros(steps)
    reference = np.zeros(steps)
    error = np.zeros(steps)
    control = np.zeros(steps)
    parameters = np.zeros((steps, theta_hat.size))

    x = 0.0
    velocity = 0.0
    velocity_estimate = 0.0
    previous_measurement = 0.0

    for k, t in enumerate(time):
        xd, vd, ad = reference_trajectory(float(t))

        # Deterministic nanometer-scale measurement noise plus encoder quantization.
        sensor_noise = (
            4.0e-9 * np.sin(2.0 * np.pi * 997.0 * t)
            + 2.0e-9 * np.cos(2.0 * np.pi * 619.0 * t)
        )
        measured_position = (
            np.round((x + sensor_noise) / encoder_resolution) * encoder_resolution
        )
        raw_velocity = (measured_position - previous_measurement) / dt
        velocity_estimate = (
            alpha_v * velocity_estimate + (1.0 - alpha_v) * raw_velocity
        )
        previous_measurement = measured_position

        e = measured_position - xd
        e_dot = velocity_estimate - vd
        sliding_error = e_dot + lambda_e * e
        reference_velocity = vd - lambda_e * e
        reference_acceleration = ad - lambda_e * e_dot

        regressor = np.array(
            [
                reference_acceleration,
                reference_velocity,
                measured_position,
                np.tanh(velocity_estimate / velocity_scale),
                np.sin(2.0 * np.pi * measured_position / ripple_pitch),
                np.cos(2.0 * np.pi * measured_position / ripple_pitch),
            ]
        )

        unsaturated_control = float(theta_hat @ regressor - k_s * sliding_error)
        applied_control = float(
            np.clip(unsaturated_control, -control_limit, control_limit)
        )

        # Freeze adaptation while saturated; otherwise use normalized adaptation,
        # sigma modification, and componentwise projection.
        if adaptive and abs(unsaturated_control) <= control_limit:
            normalization = 1.0 + float(regressor @ regressor)
            theta_dot = (
                -gamma * regressor * sliding_error / normalization
                - sigma * (theta_hat - theta_nominal)
            )
            theta_hat = np.clip(theta_hat + dt * theta_dot, lower, upper)

        periodic_friction = (
            theta_true[3] * np.tanh(velocity / velocity_scale)
            + theta_true[4] * np.sin(2.0 * np.pi * x / ripple_pitch)
            + theta_true[5] * np.cos(2.0 * np.pi * x / ripple_pitch)
        )
        external_disturbance = 0.04 * np.sin(2.0 * np.pi * 73.0 * t)
        acceleration = (
            applied_control
            - theta_true[1] * velocity
            - theta_true[2] * x
            - periodic_friction
            + external_disturbance
        ) / theta_true[0]

        velocity += dt * acceleration
        x += dt * velocity

        reference[k] = xd
        position[k] = x
        error[k] = e
        control[k] = applied_control
        parameters[k, :] = theta_hat

    return Result(time, reference, position, error, control, parameters)


def report(name: str, result: Result, settling_time: float = 0.2) -> None:
    mask = result.time >= settling_time
    rms_um = 1.0e6 * np.sqrt(np.mean(result.error[mask] ** 2))
    peak_um = 1.0e6 * np.max(np.abs(result.error[mask]))
    max_force = np.max(np.abs(result.control))
    print(
        f"{name:>8s}: RMS error = {rms_um:8.3f} um, "
        f"peak error = {peak_um:8.3f} um, max |u| = {max_force:6.3f} N"
    )


def save_csv(path: Path, fixed: Result, adaptive: Result) -> None:
    data = np.column_stack(
        [
            fixed.time,
            fixed.reference,
            fixed.position,
            adaptive.position,
            fixed.error,
            adaptive.error,
            fixed.control,
            adaptive.control,
        ]
    )
    header = (
        "time_s,reference_m,fixed_position_m,adaptive_position_m,"
        "fixed_error_m,adaptive_error_m,fixed_control_N,adaptive_control_N"
    )
    np.savetxt(path, data, delimiter=",", header=header, comments="")


def main() -> None:
    fixed = simulate(adaptive=False)
    adaptive = simulate(adaptive=True)
    report("Fixed", fixed)
    report("Adaptive", adaptive)

    output_path = Path("Chapter25_Lesson4_results.csv")
    save_csv(output_path, fixed, adaptive)
    print(f"Saved {output_path.resolve()}")
    print("Final adaptive parameter estimate:", adaptive.parameters[-1])

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(fixed.time, 1.0e6 * fixed.error, label="Fixed nominal")
        plt.plot(adaptive.time, 1.0e6 * adaptive.error, label="Adaptive")
        plt.xlabel("Time (s)")
        plt.ylabel("Tracking error (um)")
        plt.title("High-Precision Stage Tracking Error")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        figure_path = Path("Chapter25_Lesson4_tracking_error.png")
        plt.savefig(figure_path, dpi=160)
        print(f"Saved {figure_path.resolve()}")
        if "agg" not in plt.get_backend().lower():
            plt.show()
    except ImportError:
        print("Matplotlib is not installed; CSV output is still available.")


if __name__ == "__main__":
    main()

Chapter25_Lesson4.cpp

This implementation uses only the C++17 standard library and writes a CSV file suitable for plotting in any numerical package.

// Chapter 25, Lesson 4: Adaptive Control in High-Precision Motion Systems.
// Standard-library-only simulation; writes Chapter25_Lesson4_results_cpp.csv.

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

namespace {
constexpr double kPi = 3.14159265358979323846;
constexpr std::size_t kParameterCount = 6;
using Vector6 = std::array<double, kParameterCount>;

struct Result {
    std::vector<double> time;
    std::vector<double> reference;
    std::vector<double> position;
    std::vector<double> error;
    std::vector<double> control;
    Vector6 final_parameters{};
};

std::array<double, 3> referenceTrajectory(double t) {
    const double a1 = 4.0e-4;
    const double f1 = 2.0;
    const double a2 = 1.5e-4;
    const double f2 = 5.0;
    const double w1 = 2.0 * kPi * f1;
    const double w2 = 2.0 * kPi * f2;
    const double xd = a1 * std::sin(w1 * t) + a2 * std::sin(w2 * t);
    const double vd = a1 * w1 * std::cos(w1 * t) + a2 * w2 * std::cos(w2 * t);
    const double ad = -a1 * w1 * w1 * std::sin(w1 * t)
                    - a2 * w2 * w2 * std::sin(w2 * t);
    return {xd, vd, ad};
}

double dot(const Vector6& a, const Vector6& b) {
    double value = 0.0;
    for (std::size_t i = 0; i < kParameterCount; ++i) {
        value += a[i] * b[i];
    }
    return value;
}

Result simulate(bool adaptive, double dt = 1.0e-4, double duration = 2.0) {
    const int steps = static_cast<int>(std::llround(duration / dt));
    const Vector6 theta_true{0.65, 18.0, 250.0, 0.55, 0.32, -0.20};
    Vector6 theta_hat{0.45, 10.0, 180.0, 0.20, 0.0, 0.0};
    const Vector6 theta_nominal = theta_hat;
    const Vector6 lower{0.20, 0.0, 50.0, 0.0, -1.0, -1.0};
    const Vector6 upper{1.20, 50.0, 500.0, 1.5, 1.0, 1.0};
    const Vector6 gamma{2.0e2, 5.0e3, 5.0e2, 2.0e3, 4.0e2, 4.0e2};

    const double sigma = 0.05;
    const double lambda_e = 300.0;
    const double k_s = 55.0;
    const double velocity_scale = 2.0e-4;
    const double ripple_pitch = 5.0e-3;
    const double control_limit = 25.0;
    const double encoder_resolution = 1.0e-8;
    const double alpha_v = std::exp(-2.0 * kPi * 400.0 * dt);

    Result result;
    result.time.resize(steps);
    result.reference.resize(steps);
    result.position.resize(steps);
    result.error.resize(steps);
    result.control.resize(steps);

    double x = 0.0;
    double velocity = 0.0;
    double velocity_estimate = 0.0;
    double previous_measurement = 0.0;

    for (int k = 0; k < steps; ++k) {
        const double t = static_cast<double>(k) * dt;
        const auto ref = referenceTrajectory(t);
        const double xd = ref[0];
        const double vd = ref[1];
        const double ad = ref[2];

        const double sensor_noise = 4.0e-9 * std::sin(2.0 * kPi * 997.0 * t)
                                  + 2.0e-9 * std::cos(2.0 * kPi * 619.0 * t);
        const double measured_position =
            std::round((x + sensor_noise) / encoder_resolution) * encoder_resolution;
        const double raw_velocity = (measured_position - previous_measurement) / dt;
        velocity_estimate = alpha_v * velocity_estimate
                          + (1.0 - alpha_v) * raw_velocity;
        previous_measurement = measured_position;

        const double e = measured_position - xd;
        const double e_dot = velocity_estimate - vd;
        const double sliding_error = e_dot + lambda_e * e;
        const double reference_velocity = vd - lambda_e * e;
        const double reference_acceleration = ad - lambda_e * e_dot;

        const Vector6 regressor{
            reference_acceleration,
            reference_velocity,
            measured_position,
            std::tanh(velocity_estimate / velocity_scale),
            std::sin(2.0 * kPi * measured_position / ripple_pitch),
            std::cos(2.0 * kPi * measured_position / ripple_pitch)};

        const double unsaturated_control = dot(theta_hat, regressor) - k_s * sliding_error;
        const double applied_control =
            std::clamp(unsaturated_control, -control_limit, control_limit);

        if (adaptive && std::abs(unsaturated_control) <= control_limit) {
            const double normalization = 1.0 + dot(regressor, regressor);
            for (std::size_t i = 0; i < kParameterCount; ++i) {
                const double theta_dot =
                    -gamma[i] * regressor[i] * sliding_error / normalization
                    - sigma * (theta_hat[i] - theta_nominal[i]);
                theta_hat[i] = std::clamp(theta_hat[i] + dt * theta_dot,
                                          lower[i], upper[i]);
            }
        }

        const double periodic_friction =
            theta_true[3] * std::tanh(velocity / velocity_scale)
            + theta_true[4] * std::sin(2.0 * kPi * x / ripple_pitch)
            + theta_true[5] * std::cos(2.0 * kPi * x / ripple_pitch);
        const double external_disturbance = 0.04 * std::sin(2.0 * kPi * 73.0 * t);
        const double acceleration =
            (applied_control - theta_true[1] * velocity - theta_true[2] * x
             - periodic_friction + external_disturbance) / theta_true[0];

        velocity += dt * acceleration;
        x += dt * velocity;

        result.time[k] = t;
        result.reference[k] = xd;
        result.position[k] = x;
        result.error[k] = e;
        result.control[k] = applied_control;
    }

    result.final_parameters = theta_hat;
    return result;
}

void report(const std::string& name, const Result& result, double settling_time = 0.2) {
    double squared_sum = 0.0;
    double peak = 0.0;
    double max_force = 0.0;
    std::size_t count = 0;
    for (std::size_t i = 0; i < result.time.size(); ++i) {
        max_force = std::max(max_force, std::abs(result.control[i]));
        if (result.time[i] >= settling_time) {
            squared_sum += result.error[i] * result.error[i];
            peak = std::max(peak, std::abs(result.error[i]));
            ++count;
        }
    }
    if (count == 0) {
        throw std::runtime_error("No samples remain after the settling interval.");
    }
    const double rms_um = 1.0e6 * std::sqrt(squared_sum / static_cast<double>(count));
    const double peak_um = 1.0e6 * peak;
    std::cout << std::setw(8) << name << ": RMS error = " << std::fixed
              << std::setprecision(3) << std::setw(8) << rms_um
              << " um, peak error = " << std::setw(8) << peak_um
              << " um, max |u| = " << std::setw(6) << max_force << " N\n";
}

void saveCsv(const std::string& filename, const Result& fixed, const Result& adaptive) {
    std::ofstream file(filename);
    if (!file) {
        throw std::runtime_error("Could not open output CSV file.");
    }
    file << "time_s,reference_m,fixed_position_m,adaptive_position_m,"
            "fixed_error_m,adaptive_error_m,fixed_control_N,adaptive_control_N\n";
    file << std::setprecision(12);
    for (std::size_t i = 0; i < fixed.time.size(); ++i) {
        file << fixed.time[i] << ',' << fixed.reference[i] << ','
             << fixed.position[i] << ',' << adaptive.position[i] << ','
             << fixed.error[i] << ',' << adaptive.error[i] << ','
             << fixed.control[i] << ',' << adaptive.control[i] << '\n';
    }
}
}  // namespace

int main() {
    try {
        const Result fixed = simulate(false);
        const Result adaptive = simulate(true);
        report("Fixed", fixed);
        report("Adaptive", adaptive);
        saveCsv("Chapter25_Lesson4_results_cpp.csv", fixed, adaptive);

        std::cout << "Final adaptive parameter estimate:";
        for (double value : adaptive.final_parameters) {
            std::cout << ' ' << std::setprecision(8) << value;
        }
        std::cout << "\nSaved Chapter25_Lesson4_results_cpp.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Chapter25_Lesson4.java

The Java version uses only the standard library and reproduces the same controller, projection bounds, quantization, and benchmark metrics.

// Chapter 25, Lesson 4: Adaptive Control in High-Precision Motion Systems.
// Standard-library-only simulation; writes Chapter25_Lesson4_results_java.csv.

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter25_Lesson4 {
    private static final double PI = Math.PI;
    private static final int PARAMETER_COUNT = 6;

    private Chapter25_Lesson4() {}

    private static final class Result {
        final double[] time;
        final double[] reference;
        final double[] position;
        final double[] error;
        final double[] control;
        final double[] finalParameters;

        Result(int steps, double[] finalParameters) {
            this.time = new double[steps];
            this.reference = new double[steps];
            this.position = new double[steps];
            this.error = new double[steps];
            this.control = new double[steps];
            this.finalParameters = finalParameters;
        }
    }

    private static double[] referenceTrajectory(double t) {
        double a1 = 4.0e-4;
        double f1 = 2.0;
        double a2 = 1.5e-4;
        double f2 = 5.0;
        double w1 = 2.0 * PI * f1;
        double w2 = 2.0 * PI * f2;
        double xd = a1 * Math.sin(w1 * t) + a2 * Math.sin(w2 * t);
        double vd = a1 * w1 * Math.cos(w1 * t) + a2 * w2 * Math.cos(w2 * t);
        double ad = -a1 * w1 * w1 * Math.sin(w1 * t)
                  - a2 * w2 * w2 * Math.sin(w2 * t);
        return new double[] {xd, vd, ad};
    }

    private static double dot(double[] a, double[] b) {
        double value = 0.0;
        for (int i = 0; i < PARAMETER_COUNT; ++i) {
            value += a[i] * b[i];
        }
        return value;
    }

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

    private static Result simulate(boolean adaptive, double dt, double duration) {
        int steps = (int) Math.round(duration / dt);
        double[] thetaTrue = {0.65, 18.0, 250.0, 0.55, 0.32, -0.20};
        double[] thetaHat = {0.45, 10.0, 180.0, 0.20, 0.0, 0.0};
        double[] thetaNominal = thetaHat.clone();
        double[] lower = {0.20, 0.0, 50.0, 0.0, -1.0, -1.0};
        double[] upper = {1.20, 50.0, 500.0, 1.5, 1.0, 1.0};
        double[] gamma = {2.0e2, 5.0e3, 5.0e2, 2.0e3, 4.0e2, 4.0e2};

        double sigma = 0.05;
        double lambdaE = 300.0;
        double kS = 55.0;
        double velocityScale = 2.0e-4;
        double ripplePitch = 5.0e-3;
        double controlLimit = 25.0;
        double encoderResolution = 1.0e-8;
        double alphaV = Math.exp(-2.0 * PI * 400.0 * dt);

        Result result = new Result(steps, thetaHat);
        double x = 0.0;
        double velocity = 0.0;
        double velocityEstimate = 0.0;
        double previousMeasurement = 0.0;

        for (int k = 0; k < steps; ++k) {
            double t = k * dt;
            double[] ref = referenceTrajectory(t);
            double xd = ref[0];
            double vd = ref[1];
            double ad = ref[2];

            double sensorNoise = 4.0e-9 * Math.sin(2.0 * PI * 997.0 * t)
                               + 2.0e-9 * Math.cos(2.0 * PI * 619.0 * t);
            double measuredPosition =
                Math.rint((x + sensorNoise) / encoderResolution) * encoderResolution;
            double rawVelocity = (measuredPosition - previousMeasurement) / dt;
            velocityEstimate = alphaV * velocityEstimate
                             + (1.0 - alphaV) * rawVelocity;
            previousMeasurement = measuredPosition;

            double e = measuredPosition - xd;
            double eDot = velocityEstimate - vd;
            double slidingError = eDot + lambdaE * e;
            double referenceVelocity = vd - lambdaE * e;
            double referenceAcceleration = ad - lambdaE * eDot;

            double[] regressor = {
                referenceAcceleration,
                referenceVelocity,
                measuredPosition,
                Math.tanh(velocityEstimate / velocityScale),
                Math.sin(2.0 * PI * measuredPosition / ripplePitch),
                Math.cos(2.0 * PI * measuredPosition / ripplePitch)
            };

            double unsaturatedControl = dot(thetaHat, regressor) - kS * slidingError;
            double appliedControl = clamp(unsaturatedControl, -controlLimit, controlLimit);

            if (adaptive && Math.abs(unsaturatedControl) <= controlLimit) {
                double normalization = 1.0 + dot(regressor, regressor);
                for (int i = 0; i < PARAMETER_COUNT; ++i) {
                    double thetaDot =
                        -gamma[i] * regressor[i] * slidingError / normalization
                        - sigma * (thetaHat[i] - thetaNominal[i]);
                    thetaHat[i] = clamp(thetaHat[i] + dt * thetaDot, lower[i], upper[i]);
                }
            }

            double periodicFriction =
                thetaTrue[3] * Math.tanh(velocity / velocityScale)
                + thetaTrue[4] * Math.sin(2.0 * PI * x / ripplePitch)
                + thetaTrue[5] * Math.cos(2.0 * PI * x / ripplePitch);
            double externalDisturbance = 0.04 * Math.sin(2.0 * PI * 73.0 * t);
            double acceleration =
                (appliedControl - thetaTrue[1] * velocity - thetaTrue[2] * x
                 - periodicFriction + externalDisturbance) / thetaTrue[0];

            velocity += dt * acceleration;
            x += dt * velocity;

            result.time[k] = t;
            result.reference[k] = xd;
            result.position[k] = x;
            result.error[k] = e;
            result.control[k] = appliedControl;
        }
        return result;
    }

    private static void report(String name, Result result, double settlingTime) {
        double squaredSum = 0.0;
        double peak = 0.0;
        double maxForce = 0.0;
        int count = 0;
        for (int i = 0; i < result.time.length; ++i) {
            maxForce = Math.max(maxForce, Math.abs(result.control[i]));
            if (result.time[i] >= settlingTime) {
                squaredSum += result.error[i] * result.error[i];
                peak = Math.max(peak, Math.abs(result.error[i]));
                ++count;
            }
        }
        if (count == 0) {
            throw new IllegalArgumentException("No samples remain after settling time.");
        }
        double rmsMicrometers = 1.0e6 * Math.sqrt(squaredSum / count);
        double peakMicrometers = 1.0e6 * peak;
        System.out.printf(
            Locale.US,
            "%8s: RMS error = %8.3f um, peak error = %8.3f um, max |u| = %6.3f N%n",
            name, rmsMicrometers, peakMicrometers, maxForce
        );
    }

    private static void saveCsv(Path path, Result fixed, Result adaptive)
            throws IOException {
        try (BufferedWriter writer = Files.newBufferedWriter(path)) {
            writer.write(
                "time_s,reference_m,fixed_position_m,adaptive_position_m,"
                + "fixed_error_m,adaptive_error_m,fixed_control_N,adaptive_control_N\n"
            );
            for (int i = 0; i < fixed.time.length; ++i) {
                writer.write(String.format(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    fixed.time[i], fixed.reference[i], fixed.position[i],
                    adaptive.position[i], fixed.error[i], adaptive.error[i],
                    fixed.control[i], adaptive.control[i]
                ));
            }
        }
    }

    public static void main(String[] args) {
        try {
            Result fixed = simulate(false, 1.0e-4, 2.0);
            Result adaptive = simulate(true, 1.0e-4, 2.0);
            report("Fixed", fixed, 0.2);
            report("Adaptive", adaptive, 0.2);
            Path output = Path.of("Chapter25_Lesson4_results_java.csv");
            saveCsv(output, fixed, adaptive);

            System.out.print("Final adaptive parameter estimate:");
            for (double value : adaptive.finalParameters) {
                System.out.printf(Locale.US, " %.8g", value);
            }
            System.out.println();
            System.out.println("Saved " + output.toAbsolutePath());
        } catch (IOException | RuntimeException error) {
            System.err.println("Error: " + error.getMessage());
            System.exit(1);
        }
    }
}

Chapter25_Lesson4.m

The script runs the numerical experiment, generates the error plot, writes a CSV file, and contains an optional programmatic Simulink scaffold. Set BUILD_SIMULINK_MODEL = true to create the model shell, then populate its MATLAB Function blocks with the equations already implemented in the script.

%% Chapter25_Lesson4.m
% Adaptive Control in High-Precision Motion Systems
% Compares a fixed nominal controller with a projection-based adaptive law.

clear; clc; close all;

fixedResult = simulatePrecisionStage(false);
adaptiveResult = simulatePrecisionStage(true);

reportResult("Fixed", fixedResult, 0.2);
reportResult("Adaptive", adaptiveResult, 0.2);

data = [fixedResult.time, fixedResult.reference, fixedResult.position, ...
    adaptiveResult.position, fixedResult.error, adaptiveResult.error, ...
    fixedResult.control, adaptiveResult.control];
headers = {"time_s", "reference_m", "fixed_position_m", ...
    "adaptive_position_m", "fixed_error_m", "adaptive_error_m", ...
    "fixed_control_N", "adaptive_control_N"};
outputTable = array2table(data, "VariableNames", headers);
writetable(outputTable, "Chapter25_Lesson4_results_matlab.csv");

disp("Final adaptive parameter estimate:");
disp(adaptiveResult.parameters(end, :));

figure;
plot(fixedResult.time, 1e6 * fixedResult.error, "LineWidth", 1.0); hold on;
plot(adaptiveResult.time, 1e6 * adaptiveResult.error, "LineWidth", 1.0);
grid on;
xlabel("Time (s)");
ylabel("Tracking error (um)");
title("High-Precision Stage Tracking Error");
legend("Fixed nominal", "Adaptive", "Location", "best");

%% Optional Simulink scaffold
% Set BUILD_SIMULINK_MODEL to true to create a lightweight model shell. The
% Adaptive Controller and Precision Stage blocks are MATLAB Function blocks
% whose equations can be copied from simulatePrecisionStage below.
BUILD_SIMULINK_MODEL = false;
if BUILD_SIMULINK_MODEL
    modelName = "Chapter25_Lesson4_Simulink";
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);
    add_block("simulink/Sources/Clock", modelName + "/Clock", ...
        "Position", [40 80 70 110]);
    add_block("simulink/User-Defined Functions/MATLAB Function", ...
        modelName + "/Reference", "Position", [110 55 230 135]);
    add_block("simulink/User-Defined Functions/MATLAB Function", ...
        modelName + "/Adaptive Controller", "Position", [290 40 450 150]);
    add_block("simulink/User-Defined Functions/MATLAB Function", ...
        modelName + "/Precision Stage", "Position", [510 40 650 150]);
    add_block("simulink/Sinks/Scope", modelName + "/Scope", ...
        "Position", [720 70 750 100]);
    add_line(modelName, "Clock/1", "Reference/1");
    add_line(modelName, "Reference/1", "Adaptive Controller/1");
    add_line(modelName, "Adaptive Controller/1", "Precision Stage/1");
    add_line(modelName, "Precision Stage/1", "Scope/1");
    save_system(modelName);
    fprintf("Created %s.slx. Populate the MATLAB Function blocks with the equations in this file.\n", modelName);
end

function result = simulatePrecisionStage(adaptive)
    dt = 1e-4;
    duration = 2.0;
    time = (0:dt:(duration - dt))';
    steps = numel(time);

    thetaTrue = [0.65; 18.0; 250.0; 0.55; 0.32; -0.20];
    thetaHat = [0.45; 10.0; 180.0; 0.20; 0.0; 0.0];
    thetaNominal = thetaHat;
    lower = [0.20; 0.0; 50.0; 0.0; -1.0; -1.0];
    upper = [1.20; 50.0; 500.0; 1.5; 1.0; 1.0];
    gamma = [2e2; 5e3; 5e2; 2e3; 4e2; 4e2];

    sigma = 0.05;
    lambdaE = 300.0;
    kS = 55.0;
    velocityScale = 2e-4;
    ripplePitch = 5e-3;
    controlLimit = 25.0;
    encoderResolution = 1e-8;
    alphaV = exp(-2 * pi * 400 * dt);

    reference = zeros(steps, 1);
    position = zeros(steps, 1);
    errorSignal = zeros(steps, 1);
    control = zeros(steps, 1);
    parameters = zeros(steps, 6);

    x = 0.0;
    velocity = 0.0;
    velocityEstimate = 0.0;
    previousMeasurement = 0.0;

    for k = 1:steps
        t = time(k);
        [xd, vd, ad] = referenceTrajectory(t);

        sensorNoise = 4e-9 * sin(2 * pi * 997 * t) ...
            + 2e-9 * cos(2 * pi * 619 * t);
        measuredPosition = round((x + sensorNoise) / encoderResolution) ...
            * encoderResolution;
        rawVelocity = (measuredPosition - previousMeasurement) / dt;
        velocityEstimate = alphaV * velocityEstimate ...
            + (1 - alphaV) * rawVelocity;
        previousMeasurement = measuredPosition;

        e = measuredPosition - xd;
        eDot = velocityEstimate - vd;
        slidingError = eDot + lambdaE * e;
        referenceVelocity = vd - lambdaE * e;
        referenceAcceleration = ad - lambdaE * eDot;

        regressor = [referenceAcceleration; referenceVelocity; measuredPosition; ...
            tanh(velocityEstimate / velocityScale); ...
            sin(2 * pi * measuredPosition / ripplePitch); ...
            cos(2 * pi * measuredPosition / ripplePitch)];

        unsaturatedControl = thetaHat' * regressor - kS * slidingError;
        appliedControl = min(max(unsaturatedControl, -controlLimit), controlLimit);

        if adaptive && abs(unsaturatedControl) <= controlLimit
            normalization = 1 + regressor' * regressor;
            thetaDot = -gamma .* regressor * slidingError / normalization ...
                - sigma * (thetaHat - thetaNominal);
            thetaHat = min(max(thetaHat + dt * thetaDot, lower), upper);
        end

        periodicFriction = thetaTrue(4) * tanh(velocity / velocityScale) ...
            + thetaTrue(5) * sin(2 * pi * x / ripplePitch) ...
            + thetaTrue(6) * cos(2 * pi * x / ripplePitch);
        externalDisturbance = 0.04 * sin(2 * pi * 73 * t);
        acceleration = (appliedControl - thetaTrue(2) * velocity ...
            - thetaTrue(3) * x - periodicFriction + externalDisturbance) ...
            / thetaTrue(1);

        velocity = velocity + dt * acceleration;
        x = x + dt * velocity;

        reference(k) = xd;
        position(k) = x;
        errorSignal(k) = e;
        control(k) = appliedControl;
        parameters(k, :) = thetaHat';
    end

    result = struct("time", time, "reference", reference, ...
        "position", position, "error", errorSignal, ...
        "control", control, "parameters", parameters);
end

function [xd, vd, ad] = referenceTrajectory(t)
    a1 = 4e-4; f1 = 2.0;
    a2 = 1.5e-4; f2 = 5.0;
    w1 = 2 * pi * f1;
    w2 = 2 * pi * f2;
    xd = a1 * sin(w1 * t) + a2 * sin(w2 * t);
    vd = a1 * w1 * cos(w1 * t) + a2 * w2 * cos(w2 * t);
    ad = -a1 * w1^2 * sin(w1 * t) - a2 * w2^2 * sin(w2 * t);
end

function reportResult(name, result, settlingTime)
    mask = result.time >= settlingTime;
    rmsUm = 1e6 * sqrt(mean(result.error(mask).^2));
    peakUm = 1e6 * max(abs(result.error(mask)));
    maxForce = max(abs(result.control));
    fprintf("%8s: RMS error = %8.3f um, peak error = %8.3f um, max |u| = %6.3f N\n", ...
        name, rmsUm, peakUm, maxForce);
end

Chapter25_Lesson4.nb

The downloadable file is a textual Mathematica notebook expression. Opening it in the Wolfram front end displays one initialization cell containing the complete simulation, plotting, and CSV-export workflow.

Notebook[{
 Cell["Chapter 25, Lesson 4", "Title"],
 Cell["Adaptive Control in High-Precision Motion Systems", "Subtitle"],
 Cell["The input cell simulates a precision stage with uncertain inertia, damping, stiffness, Coulomb friction, position-periodic force ripple, encoder quantization, and measurement noise. It compares fixed nominal and adaptive controllers.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];

referenceTrajectory[t_] := Module[{a1, f1, a2, f2, w1, w2, xd, vd, ad},
  a1 = 4.0*^-4; f1 = 2.0;
  a2 = 1.5*^-4; f2 = 5.0;
  w1 = 2 Pi f1; w2 = 2 Pi f2;
  xd = a1 Sin[w1 t] + a2 Sin[w2 t];
  vd = a1 w1 Cos[w1 t] + a2 w2 Cos[w2 t];
  ad = -a1 w1^2 Sin[w1 t] - a2 w2^2 Sin[w2 t];
  {xd, vd, ad}
];

simulate[adaptive_] := Module[
  {dt = 1.0*^-4, duration = 2.0, steps, thetaTrue, thetaHat,
   thetaNominal, lower, upper, gamma, sigma = 0.05, lambdaE = 300.0,
   kS = 55.0, velocityScale = 2.0*^-4, ripplePitch = 5.0*^-3,
   controlLimit = 25.0, encoderResolution = 1.0*^-8, alphaV,
   x = 0.0, velocity = 0.0, velocityEstimate = 0.0,
   previousMeasurement = 0.0, rows = {}, t, xd, vd, ad, sensorNoise,
   measuredPosition, rawVelocity, e, eDot, slidingError,
   referenceVelocity, referenceAcceleration, regressor,
   unsaturatedControl, appliedControl, normalization, thetaDot,
   periodicFriction, externalDisturbance, acceleration},

  steps = Round[duration/dt];
  thetaTrue = {0.65, 18.0, 250.0, 0.55, 0.32, -0.20};
  thetaHat = {0.45, 10.0, 180.0, 0.20, 0.0, 0.0};
  thetaNominal = thetaHat;
  lower = {0.20, 0.0, 50.0, 0.0, -1.0, -1.0};
  upper = {1.20, 50.0, 500.0, 1.5, 1.0, 1.0};
  gamma = {2.0*^2, 5.0*^3, 5.0*^2, 2.0*^3, 4.0*^2, 4.0*^2};
  alphaV = Exp[-2 Pi 400.0 dt];

  Do[
    t = k dt;
    {xd, vd, ad} = referenceTrajectory[t];
    sensorNoise = 4.0*^-9 Sin[2 Pi 997.0 t] + 2.0*^-9 Cos[2 Pi 619.0 t];
    measuredPosition = Round[(x + sensorNoise)/encoderResolution] encoderResolution;
    rawVelocity = (measuredPosition - previousMeasurement)/dt;
    velocityEstimate = alphaV velocityEstimate + (1 - alphaV) rawVelocity;
    previousMeasurement = measuredPosition;

    e = measuredPosition - xd;
    eDot = velocityEstimate - vd;
    slidingError = eDot + lambdaE e;
    referenceVelocity = vd - lambdaE e;
    referenceAcceleration = ad - lambdaE eDot;

    regressor = {
      referenceAcceleration,
      referenceVelocity,
      measuredPosition,
      Tanh[velocityEstimate/velocityScale],
      Sin[2 Pi measuredPosition/ripplePitch],
      Cos[2 Pi measuredPosition/ripplePitch]
    };

    unsaturatedControl = thetaHat.regressor - kS slidingError;
    appliedControl = Clip[unsaturatedControl, {-controlLimit, controlLimit}];

    If[TrueQ[adaptive] && Abs[unsaturatedControl] <= controlLimit,
      normalization = 1 + regressor.regressor;
      thetaDot = -gamma regressor slidingError/normalization
                 - sigma (thetaHat - thetaNominal);
      thetaHat = MapThread[Clip[#1, { #2, #3 }] &,
                           {thetaHat + dt thetaDot, lower, upper}];
    ];

    periodicFriction = thetaTrue[[4]] Tanh[velocity/velocityScale]
      + thetaTrue[[5]] Sin[2 Pi x/ripplePitch]
      + thetaTrue[[6]] Cos[2 Pi x/ripplePitch];
    externalDisturbance = 0.04 Sin[2 Pi 73.0 t];
    acceleration = (appliedControl - thetaTrue[[2]] velocity
      - thetaTrue[[3]] x - periodicFriction + externalDisturbance)/thetaTrue[[1]];

    velocity = velocity + dt acceleration;
    x = x + dt velocity;
    AppendTo[rows, Join[{t, xd, x, e, appliedControl}, thetaHat]],
    {k, 0, steps - 1}
  ];

  <|\"Data\" -> rows, \"FinalParameters\" -> thetaHat|>
];

fixed = simulate[False];
adaptive = simulate[True];
fixedData = fixed[\"Data\"];
adaptiveData = adaptive[\"Data\"];
startIndex = First@FirstPosition[fixedData[[All, 1]], value_ /; value >= 0.2];
fixedRMS = 10^6 Sqrt[Mean[fixedData[[startIndex ;; All, 4]]^2]];
adaptiveRMS = 10^6 Sqrt[Mean[adaptiveData[[startIndex ;; All, 4]]^2]];

Print[\"Fixed RMS error (um): \", N[fixedRMS]];
Print[\"Adaptive RMS error (um): \", N[adaptiveRMS]];
Print[\"Final adaptive parameter estimate: \", adaptive[\"FinalParameters\"]];

errorPlot = ListLinePlot[
  {
    Transpose[{fixedData[[All, 1]], 10^6 fixedData[[All, 4]]}],
    Transpose[{adaptiveData[[All, 1]], 10^6 adaptiveData[[All, 4]]}]
  },
  PlotLegends -> {\"Fixed nominal\", \"Adaptive\"},
  Frame -> True,
  FrameLabel -> {\"Time (s)\", \"Tracking error (um)\"},
  PlotRange -> All,
  ImageSize -> Large
];
Print[errorPlot];

csvRows = MapThread[
  Join,
  {
    fixedData[[All, {1, 2, 3, 4, 5}]],
    adaptiveData[[All, {3, 4, 5}]]
  }
];
Export[
  \"Chapter25_Lesson4_results_mathematica.csv\",
  Prepend[csvRows, {
    \"time_s\", \"reference_m\", \"fixed_position_m\", \"fixed_error_m\",
    \"fixed_control_N\", \"adaptive_position_m\", \"adaptive_error_m\",
    \"adaptive_control_N\"
  }],
  \"CSV\"
];"], "Input", InitializationCell -> True]
}, WindowSize -> {1200, 850}, StyleDefinitions -> "Default.nb"]

Model structure: Include only uncertainty terms whose basis functions are physically justified and measurable in real time.

Bounds: Derive projection limits from payload, motor, transmission, friction, and stiffness data rather than choosing arbitrary wide intervals.

Reference conditioning: Ensure desired velocity and acceleration are bounded and compatible with force, jerk, and travel limits.

Sampling: Measure worst-case execution time and jitter. Use the actual encoder resolution and current-loop delay in simulation.

Robustification: Add normalization, leakage or a dead zone, projection, and saturation-aware update logic before hardware testing.

Frequency-domain verification: Confirm that adaptation does not reduce phase or gain margins near unmodeled flexible modes.

Commissioning sequence: Validate the fixed nominal loop first, enable adaptation at low gain, inspect estimate trajectories, then expand the operating envelope gradually.

Acceptance metrics: Specify RMS and peak error, settling time, spectral error, maximum force, saturation duty cycle, thermal limits, and safe parameter bounds.

Problem 1 (Filtered-Error Dynamics): Starting from \(s=\dot e+\lambda e\), prove that bounded \(s\) implies bounded \(e\), and that \(s(t)\to0\) implies \(e(t)\to0\).

Solution: Solve the stable first-order equation:

\[ e(t)=e^{-\lambda t}e(0)+ \int_0^t e^{-\lambda(t-\tau)}s(\tau)\,d\tau. \]

If \(|s(t)|\leq\bar s\), then

\[ |e(t)|\leq e^{-\lambda t}|e(0)|+ \frac{\bar s}{\lambda}(1-e^{-\lambda t}), \]

so \(e\) is bounded. If \(s(t)\to0\), split the convolution at a sufficiently large time and use exponential decay of the old part plus smallness of the recent input to conclude \(e(t)\to0\).

Problem 2 (Derive the Adaptive Error Model): Derive the filtered-error dynamics for the plant in Section 2 and the control law in Section 4.

Solution: Since

\[ \dot q=s+\dot q_r,\qquad \ddot q=\dot s+\ddot q_r, \]

the plant becomes

\[ m\dot s+bs+\boldsymbol\theta^T\mathbf Y=u+d. \]

Substitute \(u=\hat{\boldsymbol\theta}^T\mathbf Y-k_s s\):

\[ m\dot s=-(b+k_s)s+ (\hat{\boldsymbol\theta}-\boldsymbol\theta)^T\mathbf Y+d, \]

which is the required model.

Problem 3 (Nominal Stability): Using the Lyapunov function in Section 5, prove asymptotic tracking for \(d=0\) under bounded reference signals.

Solution: Differentiate \(V\) and insert

\[ \dot{\hat{\boldsymbol\theta} }=-\mathbf\Gamma\mathbf Y s. \]

The cross terms cancel exactly:

\[ \dot V=-(b+k_s)s^2. \]

Thus all signals represented in \(V\) are bounded and \(s\in L_2\). Bounded regressors imply bounded \(\dot s\); Barbalat's lemma gives \(s\to0\), and Problem 1 gives \(e\to0\).

Problem 4 (Ultimate Bound): Assume \(|d(t)|\leq\bar d\). Give a sufficient condition on \(\varepsilon\) for the coefficient of \(s^2\) in the Lyapunov inequality to remain negative.

Solution: Young's inequality gives

\[ \dot V\leq- \left(b+k_s-\frac{\varepsilon}{2}\right)s^2+ \frac{\bar d^2}{2\varepsilon}. \]

Therefore choose \(0<\varepsilon<2(b+k_s)\). Outside the set

\[ |s|\leq \frac{\bar d}{\sqrt{2\varepsilon(b+k_s-\varepsilon/2)} }, \]

the derivative is negative. This gives a conservative ultimate-error estimate; leakage and projection add bounded terms but preserve UUB under their standard properties.

Problem 5 (Ripple Basis): A motor has dominant force-ripple pitch \(p\) and a measurable second harmonic. Extend the regressor and parameter vector.

Solution: Use

\[ \boldsymbol\psi_p(q)= \begin{bmatrix} \sin(2\pi q/p)&\cos(2\pi q/p)& \sin(4\pi q/p)&\cos(4\pi q/p) \end{bmatrix}^{T}, \]

with four unknown coefficients. Append these basis functions to \(\mathbf Y\) and their estimates to \(\hat{\boldsymbol\theta}\). The same Lyapunov cancellation applies because the uncertainty remains linear in the new coefficients.

Problem 6 (Quantization and Velocity Estimation): An encoder has resolution \(\Delta_q\) and velocity is computed by one-sample differencing with period \(T_s\). Estimate the worst-case quantization-induced difference magnitude.

Solution: Two consecutive quantization errors can differ by approximately one quantization interval, so

\[ |\delta\dot q|\lesssim\frac{\Delta_q}{T_s}. \]

Reducing \(T_s\) improves control update rate but increases this raw differentiated-noise scale. A low-pass differentiator, observer, or encoder-derived velocity estimate is therefore required. Its cutoff must retain the reference bandwidth while attenuating quantization and avoiding excitation of flexible modes.

Problem 7 (Why Parameters Need Not Converge): Explain why the benchmark's mass and stiffness estimates can remain close to their initial values while tracking improves substantially.

Solution: The Lyapunov proof requires bounded parameter error and drives the scalar filtered error to zero. It does not require each regressor direction to be persistently excited. On the selected small-amplitude trajectory, friction and periodic-ripple terms dominate the error and receive informative excitation, while mass and stiffness effects can be correlated with other terms. The controller can therefore find an equivalent compensation vector for the executed trajectory without uniquely identifying every physical coefficient.

High-precision motion control converts small structured uncertainties into decisive tracking errors. By expressing inertia, damping, stiffness, friction, and spatial force ripple with known regressors and unknown coefficients, a filtered-error adaptive controller can cancel dominant uncertainty while preserving a transparent Lyapunov proof. The nominal law gives asymptotic tracking; bounded disturbances and robust modifications lead to uniform ultimate boundedness. Successful deployment additionally requires quantization-aware velocity estimation, normalization, projection, saturation gating, bandwidth separation, and systematic logging. The supplied cross-language benchmark demonstrates these principles and also illustrates the important distinction between tracking convergence and full parameter identification.