The rigid-manipulator model used in the preceding lessons assigns one generalized coordinate to each actuated joint. A flexible-joint robot instead has distinct motor-side and link-side coordinates. A flexible-link robot additionally has distributed elastic deformation, theoretically described by partial differential equations and practically approximated by a finite set of vibration modes.

Let \( \boldsymbol{q} \in \mathbb{R}^n \) denote link angles and \( \boldsymbol{\theta} \in \mathbb{R}^n \) motor angles. With an elastic transmission, generally \( \boldsymbol{q} \neq \boldsymbol{\theta} \). The deflection \( \boldsymbol{\delta} = \boldsymbol{\theta}-\boldsymbol{q} \) stores potential energy and introduces lightly damped modes. High adaptation gains can therefore excite precisely the dynamics omitted by a rigid model.

flowchart TD
  A["Motor coordinate theta"] --> B["Elastic transmission: stiffness K and damping D"]
  B --> C["Link coordinate q"]
  C --> D["Flexible link modes eta_1 ... eta_r"]
  E["Motor torque u"] --> A
  D --> F["End-effector motion and vibration"]
  F --> G["Sensors: motor encoder, link encoder, strain or tip sensor"]
  G --> H["Adaptive controller"]
  H --> E
        

Flexible joint means compliance concentrated in gears, shafts, harmonic drives, tendons, or intentionally elastic actuators. Flexible link means bending or torsion distributed along a structural member. The two phenomena may coexist, but their mathematical models and sensing requirements are not identical.

A standard rigid-link, flexible-joint model is

\[ \mathbf{M}(\boldsymbol{q},\boldsymbol{a})\ddot{\boldsymbol{q}}+ \mathbf{C}(\boldsymbol{q},\dot{\boldsymbol{q}},\boldsymbol{a}) \dot{\boldsymbol{q}}+\mathbf{g}(\boldsymbol{q},\boldsymbol{a})+ \mathbf{B}_\ell\dot{\boldsymbol{q}}= \boldsymbol{\tau}_s+\boldsymbol{d}_\ell, \]

\[ \mathbf{J}_m\ddot{\boldsymbol{\theta}}+ \mathbf{B}_m\dot{\boldsymbol{\theta}}+ \boldsymbol{\tau}_s=\boldsymbol{u}+\boldsymbol{d}_m, \]

\[ \boldsymbol{\tau}_s=\mathbf{K}(\boldsymbol{\theta}-\boldsymbol{q})+ \mathbf{D}(\dot{\boldsymbol{\theta}}-\dot{\boldsymbol{q}}). \]

Here \( \mathbf{M} \) is positive definite, the usual robot identity \( \dot{\mathbf{M}}-2\mathbf{C} \) is skew-symmetric under a consistent choice of \( \mathbf{C} \), and \( \boldsymbol{a} \) collects inertial parameters. The motor inertia \( \mathbf{J}_m \), stiffness \( \mathbf{K} \), damping matrices, payload, and friction may all be uncertain.

The link equation retains the linear parameterization property

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

However, the control input \( \boldsymbol{u} \) acts on the motor equation, not directly on the link equation. Adaptive rigid-robot computed torque cannot simply be copied without an inner elastic-torque or motor-position design.

For a slender Euler-Bernoulli link with transverse deflection \( w(x,t) \), a simplified distributed model is

\[ \rho A\,\frac{\partial^2 w}{\partial t^2}+c_w\, \frac{\partial w}{\partial t}+EI\,\frac{\partial^4 w}{\partial x^4} =f(x,t), \quad 0\leq x\leq L. \]

Boundary conditions depend on whether the link is clamped, pinned, carries a tip mass, or is connected to another joint. The assumed-modes approximation uses

\[ w(x,t)\approx\sum_{i=1}^r \phi_i(x)\eta_i(t), \]

which produces a finite-dimensional Euler-Lagrange model of the form

\[ \begin{bmatrix}\mathbf{M}_{qq}&\mathbf{M}_{q\eta}\\ \mathbf{M}_{\eta q}&\mathbf{M}_{\eta\eta}\end{bmatrix} \begin{bmatrix}\ddot{\boldsymbol{q}}\\\ddot{\boldsymbol{\eta}}\end{bmatrix} +\begin{bmatrix}\boldsymbol{h}_q\\\boldsymbol{h}_\eta\end{bmatrix}+ \begin{bmatrix}\mathbf{0}&\mathbf{0}\\\mathbf{0}&\mathbf{K}_\eta\end{bmatrix} \begin{bmatrix}\boldsymbol{q}\\\boldsymbol{\eta}\end{bmatrix}+ \begin{bmatrix}\mathbf{0}&\mathbf{0}\\\mathbf{0}&\mathbf{D}_\eta\end{bmatrix} \begin{bmatrix}\dot{\boldsymbol{q}}\\\dot{\boldsymbol{\eta}}\end{bmatrix} =\begin{bmatrix}\mathbf{B}_q\\\mathbf{B}_\eta\end{bmatrix}\boldsymbol{u}. \]

Truncating after \( r \) modes creates unmodeled higher modes. A controller that suppresses the retained modes but injects energy near an omitted natural frequency can produce spillover. Thus adaptation, filtering, sampling rate, and actuator bandwidth must be selected relative to the complete structural spectrum, not only the retained model.

To isolate the adaptive idea, consider a horizontal one-link system with an elastic joint:

\[ J_\ell\ddot q+b_\ell\dot q=\tau_s, \qquad J_m\ddot\theta+b_m\dot\theta+\tau_s=u, \qquad \tau_s=K(\theta-q). \]

The true link inertia \( J_\ell \) and link damping \( b_\ell \) are treated as unknown. Stiffness \( K \) is assumed known in the executable benchmark so that the lesson can focus on adaptation of the link dynamics. Unknown stiffness requires a different parameterization because it enters both the virtual motor command and the physical elastic torque.

For a desired trajectory \( q_d(t) \), define

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

\[ \ddot q_r=\ddot q_d-\lambda\dot e, \qquad \lambda>0. \]

The filtered error \( s \) is the adaptive-control error variable. If \( s \) converges to zero, then \( \dot e+\lambda e=0 \), so the position error converges exponentially.

The desired elastic torque is chosen as

\[ \tau_{s,d}=\hat J_\ell\ddot q_r+\hat b_\ell\dot q_r-k_s s, \qquad k_s>0. \]

Because \( \tau_s=K(\theta-q) \), define a virtual motor command

\[ \theta_d=q+K^{-1}\tau_{s,d}. \]

The laboratory motor servo is

\[ u=\tau_s-k_{pm}(\theta-\theta_d)-k_{dm}(\dot\theta-\dot q_r). \]

This inner loop is deliberately simple. It does not establish exact equality \( \tau_s=\tau_{s,d} \); it makes the elastic-torque error small when motor-side gains and bandwidth are sufficiently high. Industrial designs commonly add feedforward motor acceleration, disturbance observers, torque sensing, or a rigorously backstepped inner loop.

The parameter updates are

\[ \dot{\hat J}_\ell=-\gamma_J\ddot q_r s, \qquad \dot{\hat b}_\ell=-\gamma_b\dot q_r s, \qquad \gamma_J,\gamma_b>0, \]

followed by projection onto physically meaningful intervals. Projection prevents negative inertia and limits parameter drift, but it does not create persistent excitation or guarantee convergence to the true parameters.

flowchart TD
  R["Desired link trajectory q_d"] --> F["Reference filter: q_r_dot and q_r_ddot"]
  Q["Measured link q and q_dot"] --> E["Filtered error s"]
  F --> E
  E --> A["Adaptive law: J_hat and b_hat"]
  A --> T["Desired elastic torque tau_s_d"]
  F --> T
  T --> M["Virtual motor command theta_d"]
  Q --> M
  M --> I["Fast motor-side servo"]
  TH["Measured motor theta and theta_dot"] --> I
  I --> P["Flexible-joint plant"]
  P --> Q
  P --> TH
        

First analyze the ideal reduced problem in which the inner loop realizes \( \tau_s=\tau_{s,d} \). Let \( \widetilde J=\hat J_\ell-J_\ell \) and \( \widetilde b=\hat b_\ell-b_\ell \). Since

\[ J_\ell(\dot s+\ddot q_r)+b_\ell(s+\dot q_r)=\tau_{s,d}, \]

substitution gives

\[ J_\ell\dot s=\widetilde J\ddot q_r+ \widetilde b\dot q_r-(k_s+b_\ell)s. \]

Choose

\[ V=\frac{1}{2}J_\ell s^2+ \frac{1}{2\gamma_J}\widetilde J^2+ \frac{1}{2\gamma_b}\widetilde b^2. \]

Its derivative is

\[ \dot V=s\left(\widetilde J\ddot q_r+ \widetilde b\dot q_r-(k_s+b_\ell)s\right)+ \frac{\widetilde J}{\gamma_J}\dot{\hat J}_\ell+ \frac{\widetilde b}{\gamma_b}\dot{\hat b}_\ell. \]

With the stated update laws, cross terms cancel:

\[ \dot V=-(k_s+b_\ell)s^2\leq 0. \]

Thus \( s,\widetilde J,\widetilde b \) are bounded and \( s\in L_2 \). If the desired trajectory and its derivatives are bounded, then \( \dot s \) is bounded; Barbalat's lemma yields \( s(t)\rightarrow0 \), and the stable filter \( \dot e+\lambda e=s \) yields \( e(t)\rightarrow0 \).

This proof establishes tracking, not identification. Parameter convergence requires an excitation condition on the regressor \( \boldsymbol{\phi}=[\ddot q_r\;\;\dot q_r]^T \).

Let the elastic-torque realization error be

\[ \varepsilon_\tau=\tau_s-\tau_{s,d}. \]

The actual link error dynamics become

\[ J_\ell\dot s=\widetilde J\ddot q_r+ \widetilde b\dot q_r-(k_s+b_\ell)s+\varepsilon_\tau, \]

and therefore

\[ \dot V=-(k_s+b_\ell)s^2+s\varepsilon_\tau. \]

Using Young's inequality, for any \( \alpha>0 \),

\[ s\varepsilon_\tau\leq\frac{\alpha}{2}s^2+ \frac{1}{2\alpha}\varepsilon_\tau^2. \]

If the inner loop makes \( \varepsilon_\tau \) bounded and small, the outer loop is input-to-state practically stable and the tracking error approaches a neighborhood whose size depends on the torque-tracking error. A stronger asymptotic result requires a composite Lyapunov function, backstepping, passivity, or singular-perturbation analysis of the complete motor-link dynamics.

The phrase “fast inner loop” is therefore an engineering assumption with a mathematical consequence, not a proof by itself. Its validity must be checked against the elastic natural frequency, sampling interval, delays, saturation, and unmodeled modes.

After assumed-mode truncation, a single flexible link may be represented by hub angle \( q \) and modal coordinates \( \boldsymbol{\eta} \). A control-oriented linearized model is

\[ \mathbf{M}_f(\boldsymbol{a})\ddot{\boldsymbol{z}}+ \mathbf{D}_f(\boldsymbol{a})\dot{\boldsymbol{z}}+ \mathbf{K}_f(\boldsymbol{a})\boldsymbol{z}= \mathbf{B}_f\boldsymbol{u}, \qquad \boldsymbol{z}=\begin{bmatrix}q&\boldsymbol{\eta}^T\end{bmatrix}^T. \]

Unknown payload mass changes the inertia matrix and often shifts modal frequencies. A parameter estimator can update a payload or modal parameter vector using hub, tip, strain, or acceleration measurements. Yet three difficulties arise:

1. the finite-dimensional model omits higher modes;
2. tip or strain sensing may be non-collocated with the actuator;
3. modal parameters can be weakly identifiable from ordinary trajectories.

A typical composite objective is

\[ J_c=\int_0^T\left(e_q^T\mathbf{Q}_q e_q+ \boldsymbol{\eta}^T\mathbf{Q}_\eta\boldsymbol{\eta}+ \dot{\boldsymbol{\eta}}^T\mathbf{Q}_d\dot{\boldsymbol{\eta}}+ \boldsymbol{u}^T\mathbf{R}\boldsymbol{u}\right)dt, \]

but an adaptive law must still be derived from a stability argument rather than from this performance index alone. In practice, low-pass command shaping, notch filters, strain feedback, robust leakage, and conservative adaptation bandwidth are used to prevent spillover.

Bandwidth separation. Let the dominant elastic frequency be

\[ \omega_e\approx\sqrt{K\left(\frac{1}{J_m}+ \frac{1}{J_\ell}\right)}. \]

The numerical step and sensing bandwidth must resolve this mode. Adaptation should usually be slower than the well-damped inner servo, and command spectra should avoid unmodeled modes.

Projection and leakage. A robust update can be written conceptually as

\[ \dot{\hat{\boldsymbol{a}}}= \operatorname{Proj}_\Omega\left( -\mathbf{\Gamma}\mathbf{Y}^T\boldsymbol{s}- \sigma\mathbf{\Gamma}(\hat{\boldsymbol{a}}-\boldsymbol{a}_0) \right), \quad \sigma\geq0. \]

Leakage improves boundedness under disturbance and unmodeled dynamics but generally changes exact asymptotic convergence into ultimate boundedness.

Sensing. Motor encoders alone do not reveal link deflection when compliance is significant. Link encoders, torque sensors, strain gauges, accelerometers, or observers may be required. Numerical differentiation of noisy position measurements should be replaced by filtered derivatives or observer states.

Saturation. Parameter adaptation must not interpret actuator saturation as a plant-parameter change. Freeze or modify adaptation during severe saturation, and impose explicit torque, deflection, velocity, and parameter limits.

File: Chapter25_Lesson3.py

# Chapter25_Lesson3.py
# Adaptive control of a one-DOF flexible-joint robot.
# Plant: J_l q_ddot + b_l q_dot = K(theta - q)
#        J_m theta_ddot + b_m theta_dot + K(theta - q) = u

from dataclasses import dataclass
import csv
import numpy as np
import matplotlib.pyplot as plt


@dataclass(frozen=True)
class Parameters:
    J_l: float = 0.35
    b_l: float = 0.08
    J_m: float = 0.05
    b_m: float = 0.03
    K: float = 40.0
    lam: float = 3.0
    K_s: float = 5.0
    K_pm: float = 120.0
    K_dm: float = 4.0
    gamma_J: float = 0.40
    gamma_b: float = 0.15
    J_min: float = 0.05
    J_max: float = 2.00
    b_min: float = 0.00
    b_max: float = 2.00


def reference(t: float) -> tuple[float, float, float]:
    q_d = 0.5 * np.sin(1.2 * t) + 0.2 * np.sin(0.4 * t)
    qd_d = 0.6 * np.cos(1.2 * t) + 0.08 * np.cos(0.4 * t)
    qdd_d = -0.72 * np.sin(1.2 * t) - 0.032 * np.sin(0.4 * t)
    return q_d, qd_d, qdd_d


def projection(value: float, derivative: float, low: float, high: float) -> float:
    if value <= low and derivative < 0.0:
        return 0.0
    if value >= high and derivative > 0.0:
        return 0.0
    return derivative


def dynamics(t: float, x: np.ndarray, p: Parameters) -> np.ndarray:
    q, q_dot, theta, theta_dot, J_hat, b_hat = x
    q_d, qd_d, qdd_d = reference(t)

    e = q - q_d
    e_dot = q_dot - qd_d
    q_dot_r = qd_d - p.lam * e
    q_ddot_r = qdd_d - p.lam * e_dot
    s = q_dot - q_dot_r

    tau_s = p.K * (theta - q)
    tau_s_d = J_hat * q_ddot_r + b_hat * q_dot_r - p.K_s * s
    theta_d = q + tau_s_d / p.K

    # Fast motor-side servo. It realizes the requested elastic torque approximately.
    u = tau_s - p.K_pm * (theta - theta_d) - p.K_dm * (theta_dot - q_dot_r)

    q_ddot = (tau_s - p.b_l * q_dot) / p.J_l
    theta_ddot = (u - p.b_m * theta_dot - tau_s) / p.J_m

    J_hat_dot = projection(
        J_hat, -p.gamma_J * q_ddot_r * s, p.J_min, p.J_max
    )
    b_hat_dot = projection(
        b_hat, -p.gamma_b * q_dot_r * s, p.b_min, p.b_max
    )

    return np.array(
        [q_dot, q_ddot, theta_dot, theta_ddot, J_hat_dot, b_hat_dot],
        dtype=float,
    )


def rk4_step(t: float, x: np.ndarray, dt: float, p: Parameters) -> np.ndarray:
    k1 = dynamics(t, x, p)
    k2 = dynamics(t + 0.5 * dt, x + 0.5 * dt * k1, p)
    k3 = dynamics(t + 0.5 * dt, x + 0.5 * dt * k2, p)
    k4 = dynamics(t + dt, x + dt * k3, p)
    x_next = x + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
    x_next[4] = np.clip(x_next[4], p.J_min, p.J_max)
    x_next[5] = np.clip(x_next[5], p.b_min, p.b_max)
    return x_next


def main() -> None:
    p = Parameters()
    dt, final_time = 5.0e-4, 12.0
    time = np.arange(0.0, final_time + dt, dt)
    state = np.zeros((time.size, 6), dtype=float)
    state[0] = np.array([0.0, 0.0, 0.0, 0.0, 0.18, 0.02])
    q_ref = np.zeros(time.size)

    for k, t in enumerate(time):
        q_ref[k] = reference(t)[0]
        if k + 1 < time.size:
            state[k + 1] = rk4_step(t, state[k], dt, p)

    error = state[:, 0] - q_ref
    rms_error = float(np.sqrt(np.mean(error**2)))
    print(f"RMS link tracking error: {rms_error:.6f} rad")
    print(f"Final J_hat: {state[-1, 4]:.6f} kg.m^2")
    print(f"Final b_hat: {state[-1, 5]:.6f} N.m.s/rad")

    with open("Chapter25_Lesson3_results.csv", "w", newline="", encoding="utf-8") as file:
        writer = csv.writer(file)
        writer.writerow(["t", "q", "q_ref", "theta", "J_hat", "b_hat", "error"])
        for k, t in enumerate(time):
            writer.writerow(
                [t, state[k, 0], q_ref[k], state[k, 2], state[k, 4], state[k, 5], error[k]]
            )

    plt.figure()
    plt.plot(time, q_ref, "--", label="q_d")
    plt.plot(time, state[:, 0], label="q")
    plt.plot(time, state[:, 2], label="theta", alpha=0.75)
    plt.xlabel("Time (s)")
    plt.ylabel("Angle (rad)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(time, state[:, 4], label="J_hat")
    plt.plot(time, state[:, 5], label="b_hat")
    plt.xlabel("Time (s)")
    plt.ylabel("Parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

File: Chapter25_Lesson3.cpp

// Chapter25_Lesson3.cpp
// Standard-library simulation of the same flexible-joint adaptive controller.

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

struct Parameters {
    double Jl{0.35}, bl{0.08}, Jm{0.05}, bm{0.03}, K{40.0};
    double lambda{3.0}, Ks{5.0}, Kpm{120.0}, Kdm{4.0};
    double gammaJ{0.40}, gammab{0.15};
    double Jmin{0.05}, Jmax{2.00}, bmin{0.00}, bmax{2.00};
};

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

std::array<double, 3> reference(double t) {
    return {
        0.5 * std::sin(1.2 * t) + 0.2 * std::sin(0.4 * t),
        0.6 * std::cos(1.2 * t) + 0.08 * std::cos(0.4 * t),
       -0.72 * std::sin(1.2 * t) - 0.032 * std::sin(0.4 * t)
    };
}

double projection(double value, double derivative, double low, double high) {
    if (value <= low && derivative < 0.0) return 0.0;
    if (value >= high && derivative > 0.0) return 0.0;
    return derivative;
}

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

State dynamics(double t, const State& x, const Parameters& p) {
    const double q = x[0], qdot = x[1], theta = x[2], thetadot = x[3];
    const double Jhat = x[4], bhat = x[5];
    const auto ref = reference(t);

    const double e = q - ref[0];
    const double edot = qdot - ref[1];
    const double qdotr = ref[1] - p.lambda * e;
    const double qddotr = ref[2] - p.lambda * edot;
    const double s = qdot - qdotr;

    const double tauS = p.K * (theta - q);
    const double tauSd = Jhat * qddotr + bhat * qdotr - p.Ks * s;
    const double thetaD = q + tauSd / p.K;
    const double u = tauS - p.Kpm * (theta - thetaD) - p.Kdm * (thetadot - qdotr);

    const double qddot = (tauS - p.bl * qdot) / p.Jl;
    const double thetaddot = (u - p.bm * thetadot - tauS) / p.Jm;
    const double Jhatdot = projection(Jhat, -p.gammaJ * qddotr * s, p.Jmin, p.Jmax);
    const double bhatdot = projection(bhat, -p.gammab * qdotr * s, p.bmin, p.bmax);

    return {qdot, qddot, thetadot, thetaddot, Jhatdot, bhatdot};
}

State rk4Step(double t, const State& x, double dt, const Parameters& p) {
    const State k1 = dynamics(t, x, p);
    const State k2 = dynamics(t + 0.5 * dt, add(x, k1, 0.5 * dt), p);
    const State k3 = dynamics(t + 0.5 * dt, add(x, k2, 0.5 * dt), p);
    const State k4 = dynamics(t + dt, add(x, k3, dt), p);

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = x[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
    }
    next[4] = std::clamp(next[4], p.Jmin, p.Jmax);
    next[5] = std::clamp(next[5], p.bmin, p.bmax);
    return next;
}

int main() {
    try {
        const Parameters p{};
        const double dt = 5.0e-4, finalTime = 12.0;
        State x{0.0, 0.0, 0.0, 0.0, 0.18, 0.02};
        std::ofstream csv("Chapter25_Lesson3_results_cpp.csv");
        if (!csv) throw std::runtime_error("Cannot open output CSV file.");
        csv << "t,q,q_ref,theta,J_hat,b_hat,error\n";
        csv << std::setprecision(12);

        double squaredErrorIntegral = 0.0;
        long long samples = 0;
        for (double t = 0.0; t <= finalTime + 0.5 * dt; t += dt) {
            const double qref = reference(t)[0];
            const double error = x[0] - qref;
            squaredErrorIntegral += error * error;
            ++samples;
            csv << t << ',' << x[0] << ',' << qref << ',' << x[2] << ','
                << x[4] << ',' << x[5] << ',' << error << '\n';
            x = rk4Step(t, x, dt, p);
        }

        std::cout << std::fixed << std::setprecision(6)
                  << "RMS link tracking error: "
                  << std::sqrt(squaredErrorIntegral / static_cast<double>(samples)) << " rad\n"
                  << "Final J_hat: " << x[4] << " kg.m^2\n"
                  << "Final b_hat: " << x[5] << " N.m.s/rad\n";
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }
}

File: Chapter25_Lesson3.java

// Chapter25_Lesson3.java
// Standard-library simulation of a flexible-joint adaptive controller.

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

public final class Chapter25_Lesson3 {
    private static final class Parameters {
        final double Jl = 0.35, bl = 0.08, Jm = 0.05, bm = 0.03, K = 40.0;
        final double lambda = 3.0, Ks = 5.0, Kpm = 120.0, Kdm = 4.0;
        final double gammaJ = 0.40, gammab = 0.15;
        final double Jmin = 0.05, Jmax = 2.00, bmin = 0.00, bmax = 2.00;
    }

    private static double[] reference(double t) {
        return new double[] {
            0.5 * Math.sin(1.2 * t) + 0.2 * Math.sin(0.4 * t),
            0.6 * Math.cos(1.2 * t) + 0.08 * Math.cos(0.4 * t),
           -0.72 * Math.sin(1.2 * t) - 0.032 * Math.sin(0.4 * t)
        };
    }

    private static double projection(double value, double derivative, double low, double high) {
        if (value <= low && derivative < 0.0) return 0.0;
        if (value >= high && derivative > 0.0) return 0.0;
        return derivative;
    }

    private static double[] add(double[] a, double[] b, double scale) {
        double[] out = new double[a.length];
        for (int i = 0; i < a.length; i++) out[i] = a[i] + scale * b[i];
        return out;
    }

    private static double[] dynamics(double t, double[] x, Parameters p) {
        double q = x[0], qdot = x[1], theta = x[2], thetadot = x[3];
        double Jhat = x[4], bhat = x[5];
        double[] ref = reference(t);

        double e = q - ref[0];
        double edot = qdot - ref[1];
        double qdotr = ref[1] - p.lambda * e;
        double qddotr = ref[2] - p.lambda * edot;
        double s = qdot - qdotr;

        double tauS = p.K * (theta - q);
        double tauSd = Jhat * qddotr + bhat * qdotr - p.Ks * s;
        double thetaD = q + tauSd / p.K;
        double u = tauS - p.Kpm * (theta - thetaD) - p.Kdm * (thetadot - qdotr);

        double qddot = (tauS - p.bl * qdot) / p.Jl;
        double thetaddot = (u - p.bm * thetadot - tauS) / p.Jm;
        double Jhatdot = projection(Jhat, -p.gammaJ * qddotr * s, p.Jmin, p.Jmax);
        double bhatdot = projection(bhat, -p.gammab * qdotr * s, p.bmin, p.bmax);

        return new double[] {qdot, qddot, thetadot, thetaddot, Jhatdot, bhatdot};
    }

    private static double[] rk4Step(double t, double[] x, double dt, Parameters p) {
        double[] k1 = dynamics(t, x, p);
        double[] k2 = dynamics(t + 0.5 * dt, add(x, k1, 0.5 * dt), p);
        double[] k3 = dynamics(t + 0.5 * dt, add(x, k2, 0.5 * dt), p);
        double[] k4 = dynamics(t + dt, add(x, k3, dt), p);
        double[] next = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            next[i] = x[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
        }
        next[4] = Math.max(p.Jmin, Math.min(p.Jmax, next[4]));
        next[5] = Math.max(p.bmin, Math.min(p.bmax, next[5]));
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Parameters p = new Parameters();
        double dt = 5.0e-4, finalTime = 12.0;
        double[] x = {0.0, 0.0, 0.0, 0.0, 0.18, 0.02};
        double sumSquaredError = 0.0;
        long samples = 0;

        Path output = Path.of("Chapter25_Lesson3_results_java.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(output, StandardCharsets.UTF_8)) {
            writer.write("t,q,q_ref,theta,J_hat,b_hat,error\n");
            for (double t = 0.0; t <= finalTime + 0.5 * dt; t += dt) {
                double qRef = reference(t)[0];
                double error = x[0] - qRef;
                sumSquaredError += error * error;
                samples++;
                writer.write(String.format(
                    Locale.US, "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t, x[0], qRef, x[2], x[4], x[5], error
                ));
                x = rk4Step(t, x, dt, p);
            }
        } catch (IOException ex) {
            System.err.println("Cannot write CSV: " + ex.getMessage());
            System.exit(1);
        }

        System.out.printf(Locale.US, "RMS link tracking error: %.6f rad%n",
            Math.sqrt(sumSquaredError / samples));
        System.out.printf(Locale.US, "Final J_hat: %.6f kg.m^2%n", x[4]);
        System.out.printf(Locale.US, "Final b_hat: %.6f N.m.s/rad%n", x[5]);
    }
}

File: Chapter25_Lesson3.m

% Chapter25_Lesson3.m
% Adaptive control of a one-DOF flexible-joint robot using explicit RK4.
clear; clc; close all;

p.Jl = 0.35; p.bl = 0.08; p.Jm = 0.05; p.bm = 0.03; p.K = 40.0;
p.lambda = 3.0; p.Ks = 5.0; p.Kpm = 120.0; p.Kdm = 4.0;
p.gammaJ = 0.40; p.gammab = 0.15;
p.Jmin = 0.05; p.Jmax = 2.00; p.bmin = 0.00; p.bmax = 2.00;

dt = 5e-4; tf = 12.0; t = (0:dt:tf)';
x = zeros(numel(t), 6);
x(1,:) = [0, 0, 0, 0, 0.18, 0.02];
qref = zeros(numel(t),1);

for k = 1:numel(t)
    ref = reference(t(k));
    qref(k) = ref(1);
    if k < numel(t)
        x(k+1,:) = rk4Step(t(k), x(k,:)', dt, p)';
        x(k+1,5) = min(max(x(k+1,5), p.Jmin), p.Jmax);
        x(k+1,6) = min(max(x(k+1,6), p.bmin), p.bmax);
    end
end

err = x(:,1) - qref;
fprintf('RMS link tracking error: %.6f rad\n', sqrt(mean(err.^2)));
fprintf('Final J_hat: %.6f kg.m^2\n', x(end,5));
fprintf('Final b_hat: %.6f N.m.s/rad\n', x(end,6));

results = table(t, x(:,1), qref, x(:,3), x(:,5), x(:,6), err, ...
    'VariableNames', {'t','q','q_ref','theta','J_hat','b_hat','error'});
writetable(results, 'Chapter25_Lesson3_results_matlab.csv');

figure;
plot(t, qref, '--', t, x(:,1), '-', t, x(:,3), '-');
grid on; xlabel('Time (s)'); ylabel('Angle (rad)');
legend('q_d','q','theta','Location','best');

a = figure;
plot(t, x(:,5), t, x(:,6));
grid on; xlabel('Time (s)'); ylabel('Parameter estimate');
legend('J_hat','b_hat','Location','best');

% Simulink mapping:
% 1. Two Integrator chains implement q and theta dynamics.
% 2. A MATLAB Function block implements controllerAndAdaptation().
% 3. Saturation blocks impose [Jmin,Jmax] and [bmin,bmax].
% 4. Use a fixed-step solver substantially faster than the elastic mode.

function next = rk4Step(t, x, dt, p)
    k1 = dynamics(t, x, p);
    k2 = dynamics(t + dt/2, x + dt*k1/2, p);
    k3 = dynamics(t + dt/2, x + dt*k2/2, p);
    k4 = dynamics(t + dt, x + dt*k3, p);
    next = x + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
end

function dx = dynamics(t, x, p)
    q = x(1); qdot = x(2); theta = x(3); thetadot = x(4);
    Jhat = x(5); bhat = x(6);
    ref = reference(t);
    e = q - ref(1); edot = qdot - ref(2);
    qdotr = ref(2) - p.lambda*e;
    qddotr = ref(3) - p.lambda*edot;
    s = qdot - qdotr;

    tauS = p.K*(theta - q);
    tauSd = Jhat*qddotr + bhat*qdotr - p.Ks*s;
    thetaD = q + tauSd/p.K;
    u = tauS - p.Kpm*(theta - thetaD) - p.Kdm*(thetadot - qdotr);

    qddot = (tauS - p.bl*qdot)/p.Jl;
    thetaddot = (u - p.bm*thetadot - tauS)/p.Jm;
    Jhatdot = projection(Jhat, -p.gammaJ*qddotr*s, p.Jmin, p.Jmax);
    bhatdot = projection(bhat, -p.gammab*qdotr*s, p.bmin, p.bmax);
    dx = [qdot; qddot; thetadot; thetaddot; Jhatdot; bhatdot];
end

function ref = reference(t)
    ref = [0.5*sin(1.2*t) + 0.2*sin(0.4*t); ...
           0.6*cos(1.2*t) + 0.08*cos(0.4*t); ...
          -0.72*sin(1.2*t) - 0.032*sin(0.4*t)];
end

function out = projection(value, derivative, low, high)
    if (value <= low && derivative < 0) || (value >= high && derivative > 0)
        out = 0;
    else
        out = derivative;
    end
end

In Simulink, implement the motor and link equations as separate integrator chains. Place the reference filter, desired elastic torque, projection logic, and motor servo in a MATLAB Function block. Use a fixed-step solver whose step is substantially smaller than the period of the highest retained elastic mode. Log motor-link deflection and stop the simulation if a safe deflection threshold is exceeded.

File: Chapter25_Lesson3.nb

Notebook[{
 Cell["Chapter 25, Lesson 3: Flexible-Joint Adaptive Control", "Title"],
 Cell["Chapter25_Lesson3.nb", "Subtitle"],
 Cell[BoxData@ToBoxes[
   Module[{p, dt, tf, steps, reference, projection, dynamics, rk4, x, rows,
     t, ref, error, rms, trackingPlot, parameterPlot},
    p = <|"Jl" -> 0.35, "bl" -> 0.08, "Jm" -> 0.05, "bm" -> 0.03,
      "K" -> 40., "lambda" -> 3., "Ks" -> 5., "Kpm" -> 120.,
      "Kdm" -> 4., "gammaJ" -> 0.40, "gammab" -> 0.15,
      "Jmin" -> 0.05, "Jmax" -> 2., "bmin" -> 0., "bmax" -> 2.|>;
    dt = 5.*10^-4; tf = 12.; steps = Round[tf/dt];

    reference[time_] := {0.5 Sin[1.2 time] + 0.2 Sin[0.4 time],
      0.6 Cos[1.2 time] + 0.08 Cos[0.4 time],
      -0.72 Sin[1.2 time] - 0.032 Sin[0.4 time]};

    projection[value_, derivative_, low_, high_] :=
      If[(value <= low && derivative < 0) || (value >= high && derivative > 0),
        0., derivative];

    dynamics[time_, state_] := Module[
      {q, qdot, theta, thetadot, Jhat, bhat, r, e, edot, qdotr,
       qddotr, s, tauS, tauSd, thetaD, u, qddot, thetaddot},
      {q, qdot, theta, thetadot, Jhat, bhat} = state;
      r = reference[time]; e = q - r[[1]]; edot = qdot - r[[2]];
      qdotr = r[[2]] - p["lambda"] e;
      qddotr = r[[3]] - p["lambda"] edot;
      s = qdot - qdotr;
      tauS = p["K"] (theta - q);
      tauSd = Jhat qddotr + bhat qdotr - p["Ks"] s;
      thetaD = q + tauSd/p["K"];
      u = tauS - p["Kpm"] (theta - thetaD) -
        p["Kdm"] (thetadot - qdotr);
      qddot = (tauS - p["bl"] qdot)/p["Jl"];
      thetaddot = (u - p["bm"] thetadot - tauS)/p["Jm"];
      {qdot, qddot, thetadot, thetaddot,
       projection[Jhat, -p["gammaJ"] qddotr s, p["Jmin"], p["Jmax"]],
       projection[bhat, -p["gammab"] qdotr s, p["bmin"], p["bmax"]]}
      ];

    rk4[time_, state_] := Module[{k1, k2, k3, k4, next},
      k1 = dynamics[time, state];
      k2 = dynamics[time + dt/2, state + dt k1/2];
      k3 = dynamics[time + dt/2, state + dt k2/2];
      k4 = dynamics[time + dt, state + dt k3];
      next = state + dt (k1 + 2 k2 + 2 k3 + k4)/6;
      next[[5]] = Clip[next[[5]], {p["Jmin"], p["Jmax"]}];
      next[[6]] = Clip[next[[6]], {p["bmin"], p["bmax"]}]; next];

    x = {0., 0., 0., 0., 0.18, 0.02}; rows = {};
    Do[t = k dt; ref = reference[t];
      AppendTo[rows, Join[{t}, x, {ref[[1]], x[[1]] - ref[[1]]}]];
      x = rk4[t, x], {k, 0, steps}];
    error = rows[[All, 9]]; rms = Sqrt[Mean[error^2]];
    Export["Chapter25_Lesson3_results_mathematica.csv", rows];
    Print["RMS link tracking error: ", NumberForm[rms, {8, 6}], " rad"];
    Print["Final J_hat: ", x[[5]], " kg.m^2"];
    Print["Final b_hat: ", x[[6]], " N.m.s/rad"];

    trackingPlot = ListLinePlot[
      {rows[[All, {1, 8}]], rows[[All, {1, 2}]], rows[[All, {1, 4}]]},
      PlotLegends -> {"q_d", "q", "theta"}, AxesLabel -> {"t", "angle"},
      PlotRange -> All, ImageSize -> Large];
    parameterPlot = ListLinePlot[
      {rows[[All, {1, 6}]], rows[[All, {1, 7}]]},
      PlotLegends -> {"J_hat", "b_hat"}, AxesLabel -> {"t", "estimate"},
      PlotRange -> All, ImageSize -> Large];
    Column[{trackingPlot, parameterPlot}]
   ]], "Input"]
}]

The initial estimates deliberately differ from the true values. Tracking can become accurate even while the estimates do not converge exactly. This is expected: the Lyapunov result needs bounded estimates and convergence of \( s \), whereas exact parameter recovery additionally needs persistent excitation and a model that captures the relevant dynamics.

Repeat the simulation with larger adaptation gains. Initially the response may improve, but beyond a point the parameter estimates react to the elastic transient and the motor loop must work harder. Then reduce stiffness or increase motor inertia. The elastic mode becomes slower, weakening the assumed bandwidth separation. These experiments demonstrate why “faster adaptation” is not synonymous with “better adaptive control.”

Problem 1 — Elastic natural frequency. Neglect damping and external torque. Derive the nonzero natural frequency of the one-DOF motor-link pair.

Solution. The unforced equations are

\[ J_\ell\ddot q=K(\theta-q), \qquad J_m\ddot\theta=-K(\theta-q). \]

Using a harmonic solution gives

\[ \det\begin{bmatrix}K-J_\ell\omega^2&-K\\-K&K-J_m\omega^2 \end{bmatrix}=0, \]

\[ \omega^2\left(J_\ell J_m\omega^2-K(J_\ell+J_m)\right)=0. \]

One mode is rigid-body motion with \( \omega=0 \); the elastic mode is

\[ \omega_e=\sqrt{K\left(\frac{1}{J_m}+ \frac{1}{J_\ell}\right)}. \]

Problem 2 — Lyapunov cancellation. Starting from the ideal torque-realization error dynamics, verify the parameter-update signs.

Solution. With \( \widetilde J=\hat J-J \), the Lyapunov cross term is \( s\widetilde J\ddot q_r \). The parameter-energy derivative contributes \( \widetilde J\dot{\hat J}/\gamma_J \). Choosing \( \dot{\hat J}=-\gamma_J\ddot q_r s \) cancels it. The same argument gives \( \dot{\hat b}=-\gamma_b\dot q_r s \).

Problem 3 — Practical stability with torque error. Suppose \( |\varepsilon_\tau|\leq\bar\varepsilon \). Find an ultimate bound on \( |s| \) using a direct sign argument.

Solution. Since

\[ \dot V\leq-(k_s+b_\ell)|s|^2+|s|\bar\varepsilon, \]

the derivative is negative whenever \( |s|>\bar\varepsilon/(k_s+b_\ell) \). Hence the filtered error is driven toward a neighborhood no larger than this conservative radius, subject to bounded parameter estimates and inner-loop signals.

Problem 4 — Parameter convergence. Can a constant desired position identify both \( J_\ell \) and \( b_\ell \)?

Solution. After transients, \( \dot q_r\approx0 \) and \( \ddot q_r\approx0 \), so the regressor loses rank. Tracking may be achieved, but the data contain insufficient excitation for unique identification of inertia and damping.

Problem 5 — Flexible-link truncation. Explain why increasing the number of retained modes is not automatically harmless.

Solution. Each added mode increases model order, sensing demands, numerical stiffness, and parameter dimension. Weakly observable modes can degrade estimation, and a controller designed around an inaccurate high-order model can inject excessive high-frequency control. Mode selection should reflect the command spectrum, actuator bandwidth, sensor placement, and uncertainty of omitted dynamics.

Problem 6 — Simulation study. In each implementation, set \( K=10 \) while keeping gains unchanged. Predict and then observe the effect.

Solution. Lower stiffness reduces the elastic natural frequency and increases motor-link deflection required for a given torque. The inner loop is no longer as fast relative to the outer loop, so \( \varepsilon_\tau \) grows and tracking generally deteriorates. Retuning the motor servo and reducing adaptation bandwidth may be necessary.

Flexible joints introduce separate motor and link coordinates; flexible links add distributed deformation and vibration modes. Their adaptive control requires more than substituting estimated parameters into a rigid computed-torque law. The lesson constructed a two-loop design: an adaptive link-side torque request and a fast motor-side realization loop. A Lyapunov proof establishes asymptotic link tracking only under exact elastic-torque realization; the full coupled system contains a torque-error perturbation that motivates composite, backstepping, passivity, or singular-perturbation analysis. For flexible links, assumed-mode truncation, sensor placement, excitation, and spillover are central design constraints.