A common mechatronic axis consists of a motor, a gearbox, a compliant transmission (e.g., shaft/belt), and a load (table, rotor, link). In this chapter we have already established rotational dynamics and gear transformations, so we now focus on assembling a complete motor–gear–shaft–load mechanical model and a minimal servo-drive abstraction that does not require detailed electrical-circuit modeling (introduced in Chapter 5).

We use the following variables (all functions of time): motor angle \( \theta_m \), load angle \( \theta_\ell \), motor speed \( \omega_m = \dot{\theta}_m \), load speed \( \omega_\ell = \dot{\theta}_\ell \). The motor delivers commanded torque \( \tau_m \) and the environment may apply a disturbance/load torque \( \tau_d \) at the load.

Gear ratio convention. Let \( N > 0 \) denote the (ideal) speed ratio \( \omega_m = N\,\omega_\ell^{(g)} \), where \( \omega_\ell^{(g)} \) is the gearbox output speed. For an ideal gear, power conservation implies \( \tau_m \omega_m = \tau_\ell^{(g)} \omega_\ell^{(g)} \).

flowchart TD
  R["Reference (position or speed)"] --> C["Servo logic (outer loop)"]
  C --> U["Torque command tau_m"]
  U --> M["Motor inertia Jm, viscous bm"]
  M --> G["Ideal gear ratio N"]
  G --> S["Shaft: spring Ks + damper Cs (optional)"]
  S --> L["Load inertia Jl, viscous bl, disturbance tau_d"]
  L --> Y["Measured output (theta_l or omega_l)"]
  Y --> C
        

First consider a rigid drivetrain: the gearbox output is rigidly connected to the load, so \( \theta_\ell = \theta_\ell^{(g)} \) and \( \omega_\ell = \omega_\ell^{(g)} \). The ideal gear kinematics are \( \omega_m = N\,\omega_\ell \) and \( \theta_m = N\,\theta_\ell \) (up to a constant offset).

Torque transformation (power proof). For an ideal gear (lossless),

\[ \tau_m \omega_m = \tau_\ell \omega_\ell,\quad \omega_m = N\,\omega_\ell \;\;\Rightarrow\;\; \tau_\ell = N\,\tau_m. \]

Thus, a torque at the motor is amplified at the load side by \( N \), while speed is reduced by \( N \). (If you instead “reflect” the load to the motor side, the gear reduces the apparent torque needed at the motor.)

Reflected inertia (energy proof). Total kinetic energy is

\[ T = \tfrac{1}{2}J_m \omega_m^2 + \tfrac{1}{2}J_\ell \omega_\ell^2 = \tfrac{1}{2}\Big(J_m N^2 + J_\ell\Big)\omega_\ell^2 = \tfrac{1}{2}J_{\text{eq},\ell}\,\omega_\ell^2, \]

hence the equivalent inertia “seen” at the load side is \( J_{\text{eq},\ell} = J_m N^2 + J_\ell \). Conversely, reflecting the load to the motor side yields \( J_{\text{eq},m} = J_m + \frac{J_\ell}{N^2} \), because \( \omega_\ell = \frac{\omega_m}{N} \).

Rigid-axis ODE (motor-side reflection). If viscous friction coefficients are \( b_m \) (motor) and \( b_\ell \) (load), then on the motor side:

\[ \Big(J_m + \tfrac{J_\ell}{N^2}\Big)\dot{\omega}_m + \Big(b_m + \tfrac{b_\ell}{N^2}\Big)\omega_m = \tau_m - \tfrac{\tau_d}{N}. \]

This single second-order rotational model is often adequate when the transmission is stiff and backlash is negligible (backlash was treated in Lesson 4).

High-performance servo drives frequently exhibit oscillations due to elastic transmission elements (shafts, belts, couplings). A standard minimal model is the two-inertia torsional system: motor inertia and load inertia connected by a spring-damper.

Define the gear-output angle \( \theta_g = \frac{\theta_m}{N} \) and the shaft twist \( \phi = \theta_g - \theta_\ell = \frac{\theta_m}{N} - \theta_\ell \). With torsional stiffness \( K_s \) and damping \( C_s \), the transmitted shaft torque is

\[ \tau_s = K_s \phi + C_s \dot{\phi} = K_s\Big(\tfrac{\theta_m}{N} - \theta_\ell\Big) + C_s\Big(\tfrac{\omega_m}{N} - \omega_\ell\Big). \]

Equations of motion. Apply rotational Newton’s law to each inertia (torque balance about each shaft):

\[ J_m \dot{\omega}_m = \tau_m - b_m \omega_m - \tfrac{1}{N}\tau_s, \qquad J_\ell \dot{\omega}_\ell = \tau_s - b_\ell \omega_\ell - \tau_d. \]

Together with kinematics \( \dot{\theta}_m = \omega_m \), \( \dot{\theta}_\ell = \omega_\ell \), this yields a 4th-order ODE system in the coordinates \( (\theta_m,\omega_m,\theta_\ell,\omega_\ell) \).

Undamped torsional mode (closed-form natural frequency). Set \( b_m=b_\ell=C_s=0 \) and \( \tau_m=\tau_d=0 \). Consider the relative coordinate \( \phi = \frac{\theta_m}{N} - \theta_\ell \). Using \( \tau_s = K_s \phi \) and differentiating twice, one can show that \( \phi \) satisfies

\[ \ddot{\phi} + K_s\Big(\tfrac{1}{J_m N^2} + \tfrac{1}{J_\ell}\Big)\phi = 0. \]

Therefore the torsional natural frequency is

\[ \omega_t = \sqrt{K_s\Big(\tfrac{1}{J_m N^2} + \tfrac{1}{J_\ell}\Big)}. \]

This frequency often limits achievable servo bandwidth and must be captured in the model when \( \omega_t \) is not far above the commanded motion spectrum.

Practical servo drives are typically organized as nested loops: a fast inner loop shapes actuator behavior (often making torque track a command quickly), and slower outer loops regulate velocity/position. Without invoking electrical-circuit dynamics (Chapter 5), we adopt the standard abstraction:

Assumption (fast actuation): the motor torque follows a command \( u(t) \) closely, i.e. \( \tau_m \approx u \), within the mechanical time scales of interest.

flowchart TD
  A["Position reference theta_ref"] --> B["Compute error e = theta_ref - theta_l"]
  B --> C["PD law: u = Kp*e - Kd*omega_l"]
  C --> D["Actuator behaves like torque source (tau_m ~= u)"]
  D --> E["Mechanical plant (motor + gear + load)"]
  E --> F["Output theta_l, omega_l (from encoder)"]
  F --> B
        

To keep the mathematics aligned with material covered so far, we prove stability for a rigid axis model with the output angle \( \theta_\ell \) directly actuated through an equivalent inertia. Using load-side reflection, define \( J_{\text{eq}} = J_m N^2 + J_\ell \) and \( b_{\text{eq}} = b_m N^2 + b_\ell \). The rigid load-side dynamics are

\[ J_{\text{eq}} \dot{\omega}_\ell + b_{\text{eq}} \omega_\ell = \tau_\ell - \tau_d, \qquad \dot{\theta}_\ell = \omega_\ell, \qquad \tau_\ell = N\,\tau_m. \]

Consider the PD servo law (disturbance-free here for clarity, \( \tau_d=0 \)):

\[ e = \theta_{\text{ref}} - \theta_\ell, \qquad \tau_\ell = K_p e - K_d \omega_\ell, \qquad K_p > 0,\; K_d > 0. \]

For constant reference \( \theta_{\text{ref}} \), we have \( \dot{e} = -\omega_\ell \), and the closed-loop error dynamics are

\[ J_{\text{eq}} \ddot{e} + (b_{\text{eq}} + K_d)\dot{e} + K_p e = 0. \]

Lyapunov (energy) stability proof. Define the candidate function

\[ V(e,\dot{e}) = \tfrac{1}{2}J_{\text{eq}}\dot{e}^{\,2} + \tfrac{1}{2}K_p e^2. \]

Since \( J_{\text{eq}} > 0 \) and \( K_p > 0 \), we have \( V \ge 0 \) and \( V=0 \) iff \( e=0 \) and \( \dot{e}=0 \). Differentiate along trajectories using \( J_{\text{eq}}\ddot{e} = -(b_{\text{eq}}+K_d)\dot{e} - K_p e \):

\[ \dot{V} = J_{\text{eq}}\dot{e}\ddot{e} + K_p e\dot{e} = \dot{e}\Big(-(b_{\text{eq}}+K_d)\dot{e} - K_p e\Big) + K_p e\dot{e} = -(b_{\text{eq}}+K_d)\dot{e}^{\,2} \le 0. \]

Thus the mechanical energy decays monotonically; with standard arguments (e.g., LaSalle’s invariance principle), the only invariant set satisfying \( \dot{V}=0 \) is \( \dot{e}=0 \), which forces \( e=0 \). Hence the equilibrium is asymptotically stable for \( K_p > 0 \), \( K_d > 0 \).

For simulation, assemble the two-inertia model from Section 3 and apply a PD law using measured load angle/speed. We treat the commanded motor torque as \( \tau_m = \frac{1}{N}(K_p(\theta_{\text{ref}}-\theta_\ell) - K_d \omega_\ell) \) so that the load-side torque is approximately \( \tau_\ell = N\tau_m \approx K_p e - K_d \omega_\ell \).

Define the ODE state vector as \( \mathbf{x} = [\theta_m,\omega_m,\theta_\ell,\omega_\ell]^T \). Then

\[ \dot{\mathbf{x}} = \begin{bmatrix} \omega_m \\ \frac{1}{J_m}\Big(\tau_m - b_m\omega_m - \frac{1}{N}\tau_s\Big) \\ \omega_\ell \\ \frac{1}{J_\ell}\Big(\tau_s - b_\ell\omega_\ell - \tau_d\Big) \end{bmatrix}, \quad \tau_s = K_s\Big(\tfrac{\theta_m}{N}-\theta_\ell\Big) + C_s\Big(\tfrac{\omega_m}{N}-\omega_\ell\Big). \]

A standard explicit 4th-order Runge–Kutta (RK4) step with step size \( h \) is

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \tfrac{h}{6}\big(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\big), \]

where \( \mathbf{k}_1 = f(t_k,\mathbf{x}_k) \), \( \mathbf{k}_2 = f(t_k + \tfrac{h}{2}, \mathbf{x}_k + \tfrac{h}{2}\mathbf{k}_1) \), \( \mathbf{k}_3 = f(t_k + \tfrac{h}{2}, \mathbf{x}_k + \tfrac{h}{2}\mathbf{k}_2) \), \( \mathbf{k}_4 = f(t_k + h, \mathbf{x}_k + h\mathbf{k}_3) \).

(We will revisit numerical stability and step-size selection formally in Chapter 15; here we use RK4 as a reliable baseline.)

The following implementations simulate a two-inertia servo axis with optional shaft elasticity. Suggested parameter units: inertia in kg·m², viscous damping in N·m·s/rad, stiffness in N·m/rad.

6.1 Python (NumPy + SciPy)

import numpy as np
from scipy.integrate import solve_ivp

# Parameters
Jm = 2.0e-4
Jl = 1.2e-3
bm = 2.0e-4
bl = 8.0e-4
N  = 8.0

Ks = 3.0    # set Ks=0 for rigid coupling approximation
Cs = 2.0e-2 # set Cs=0 for purely elastic

Kp = 15.0
Kd = 0.6

def theta_ref(t):
    # step considered after t=0.05 s
    return 1.0 if t >= 0.05 else 0.0

def tau_d(t):
    # disturbance torque at the load side
    return 0.0

def dynamics(t, x):
    th_m, w_m, th_l, w_l = x

    # PD servo at load side, implemented via motor torque
    e = theta_ref(t) - th_l
    tau_l_cmd = Kp * e - Kd * w_l
    tau_m = tau_l_cmd / N

    # shaft torque
    phi  = th_m / N - th_l
    dphi = w_m / N - w_l
    tau_s = Ks * phi + Cs * dphi

    # ODEs
    dth_m = w_m
    dw_m  = (tau_m - bm * w_m - (tau_s / N)) / Jm

    dth_l = w_l
    dw_l  = (tau_s - bl * w_l - tau_d(t)) / Jl

    return [dth_m, dw_m, dth_l, dw_l]

x0 = [0.0, 0.0, 0.0, 0.0]
t_span = (0.0, 0.6)
t_eval = np.linspace(t_span[0], t_span[1], 2000)

sol = solve_ivp(dynamics, t_span, x0, t_eval=t_eval, rtol=1e-8, atol=1e-10)

# Example post-processing (no plotting styling assumptions here)
theta_l = sol.y[2]
omega_l = sol.y[3]
print("Final load angle:", theta_l[-1])
print("Peak load speed:", np.max(np.abs(omega_l)))
      

Relevant libraries: scipy.integrate for ODE solving; numpy for vectorized computation. Later (Chapters 6–8), you may also use python-control to represent interconnections, but here we remain ODE-based.

6.2 C++ (RK4 from scratch + Eigen for vectors)

#include <cmath>
#include <iostream>
#include <vector>
#include <Eigen/Dense>

struct Params {
  double Jm, Jl, bm, bl, N;
  double Ks, Cs;
  double Kp, Kd;
};

static double theta_ref(double t) {
  return (t >= 0.05) ? 1.0 : 0.0;
}

static double tau_d(double /*t*/) {
  return 0.0;
}

static Eigen::Vector4d f(double t, const Eigen::Vector4d& x, const Params& p) {
  const double th_m = x(0), w_m = x(1), th_l = x(2), w_l = x(3);

  const double e = theta_ref(t) - th_l;
  const double tau_l_cmd = p.Kp * e - p.Kd * w_l;
  const double tau_m = tau_l_cmd / p.N;

  const double phi  = th_m / p.N - th_l;
  const double dphi = w_m  / p.N - w_l;
  const double tau_s = p.Ks * phi + p.Cs * dphi;

  Eigen::Vector4d dx;
  dx(0) = w_m;
  dx(1) = (tau_m - p.bm * w_m - (tau_s / p.N)) / p.Jm;
  dx(2) = w_l;
  dx(3) = (tau_s - p.bl * w_l - tau_d(t)) / p.Jl;
  return dx;
}

static Eigen::Vector4d rk4_step(double t, double h, const Eigen::Vector4d& x, const Params& p) {
  const Eigen::Vector4d k1 = f(t, x, p);
  const Eigen::Vector4d k2 = f(t + 0.5*h, x + 0.5*h*k1, p);
  const Eigen::Vector4d k3 = f(t + 0.5*h, x + 0.5*h*k2, p);
  const Eigen::Vector4d k4 = f(t + h,     x + h*k3,     p);
  return x + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
}

int main() {
  Params p;
  p.Jm = 2.0e-4; p.Jl = 1.2e-3;
  p.bm = 2.0e-4; p.bl = 8.0e-4;
  p.N  = 8.0;
  p.Ks = 3.0;    p.Cs = 2.0e-2;
  p.Kp = 15.0;   p.Kd = 0.6;

  double t0 = 0.0, tf = 0.6;
  double h  = 1.0e-4;

  Eigen::Vector4d x;
  x << 0.0, 0.0, 0.0, 0.0;

  double t = t0;
  double peak_wl = 0.0;

  while (t < tf) {
    peak_wl = std::max(peak_wl, std::abs(x(3)));
    x = rk4_step(t, h, x, p);
    t += h;
  }

  std::cout << "Final load angle theta_l: " << x(2) << "\n";
  std::cout << "Peak |omega_l|: " << peak_wl << "\n";
  return 0;
}
      

Relevant libraries: Eigen for linear algebra (vectors/matrices). The RK4 implementation is written from scratch to make the numerical integration pipeline explicit.

6.3 Java (Apache Commons Math ODE Integrator)

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince54Integrator;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;

public class TwoInertiaServo {

  static class Params {
    double Jm, Jl, bm, bl, N;
    double Ks, Cs;
    double Kp, Kd;
  }

  static double thetaRef(double t) {
    return (t >= 0.05) ? 1.0 : 0.0;
  }

  static double tauDist(double t) {
    return 0.0;
  }

  static class Plant implements FirstOrderDifferentialEquations {
    private final Params p;
    Plant(Params p) { this.p = p; }

    public int getDimension() { return 4; }

    public void computeDerivatives(double t, double[] x, double[] dx) {
      double th_m = x[0], w_m = x[1], th_l = x[2], w_l = x[3];

      double e = thetaRef(t) - th_l;
      double tauLCmd = p.Kp * e - p.Kd * w_l;
      double tauM = tauLCmd / p.N;

      double phi  = th_m / p.N - th_l;
      double dphi = w_m  / p.N - w_l;
      double tauS = p.Ks * phi + p.Cs * dphi;

      dx[0] = w_m;
      dx[1] = (tauM - p.bm * w_m - (tauS / p.N)) / p.Jm;
      dx[2] = w_l;
      dx[3] = (tauS - p.bl * w_l - tauDist(t)) / p.Jl;
    }
  }

  public static void main(String[] args) {
    Params p = new Params();
    p.Jm = 2.0e-4; p.Jl = 1.2e-3;
    p.bm = 2.0e-4; p.bl = 8.0e-4;
    p.N  = 8.0;
    p.Ks = 3.0;    p.Cs = 2.0e-2;
    p.Kp = 15.0;   p.Kd = 0.6;

    FirstOrderDifferentialEquations ode = new Plant(p);

    double t0 = 0.0, tf = 0.6;
    double[] x0 = new double[] {0.0, 0.0, 0.0, 0.0};

    DormandPrince54Integrator integrator =
        new DormandPrince54Integrator(1.0e-6, 1.0e-3, 1.0e-8, 1.0e-10);

    final double[] peakWl = new double[] {0.0};

    integrator.addStepHandler(new StepHandler() {
      public void init(double t0, double[] y0, double t) {}
      public void handleStep(StepInterpolator interpolator, boolean isLast) {
        double[] y = interpolator.getInterpolatedState();
        peakWl[0] = Math.max(peakWl[0], Math.abs(y[3]));
      }
    });

    double[] x = x0.clone();
    integrator.integrate(ode, t0, x, tf, x);

    System.out.println("Final load angle theta_l: " + x[2]);
    System.out.println("Peak |omega_l|: " + peakWl[0]);
  }
}
      

Relevant libraries: Apache Commons Math ODE package for robust adaptive integrators. For matrix-heavy extensions later, EJML is a common choice for linear algebra in Java.

6.4 MATLAB (ode45) and Simulink Implementation Notes

% Parameters
Jm = 2.0e-4; Jl = 1.2e-3;
bm = 2.0e-4; bl = 8.0e-4;
N  = 8.0;
Ks = 3.0;    Cs = 2.0e-2;
Kp = 15.0;   Kd = 0.6;

theta_ref = @(t) double(t >= 0.05);  % step to 1 rad after 0.05s
tau_d     = @(t) 0.0;

f = @(t,x) [
  x(2);
  ( (Kp*(theta_ref(t)-x(3)) - Kd*x(4))/N - bm*x(2) - ( (Ks*(x(1)/N - x(3)) + Cs*(x(2)/N - x(4)))/N ) )/Jm;
  x(4);
  ( (Ks*(x(1)/N - x(3)) + Cs*(x(2)/N - x(4))) - bl*x(4) - tau_d(t) )/Jl
];

x0 = [0;0;0;0]; % [theta_m; omega_m; theta_l; omega_l]
tspan = [0 0.6];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,x] = ode45(f, tspan, x0, opts);

theta_l = x(:,3);
omega_l = x(:,4);
disp(theta_l(end));
disp(max(abs(omega_l)));
      

Simulink build (no images): Implement the same ODEs using: (i) two cascaded integrators for motor (theta_m, omega_m); (ii) two cascaded integrators for load (theta_l, omega_l); (iii) algebraic blocks computing phi = theta_m/N - theta_l, dphi = omega_m/N - omega_l, tau_s = Ks*phi + Cs*dphi; (iv) PD block computing tau_l_cmd = Kp*(theta_ref-theta_l) - Kd*omega_l; (v) map to motor torque tau_m = tau_l_cmd/N; (vi) compute accelerations and feed the integrators. (Use a step block for theta_ref and optionally a disturbance source for tau_d.)

6.5 Wolfram Mathematica (NDSolve)

(* Parameters *)
Jm = 2.0*10^-4; Jl = 1.2*10^-3;
bm = 2.0*10^-4; bl = 8.0*10^-4;
N  = 8.0;
Ks = 3.0; Cs = 2.0*10^-2;
Kp = 15.0; Kd = 0.6;

thetaRef[t_] := If[t >= 0.05, 1.0, 0.0];
tauD[t_] := 0.0;

(* States: thetaM[t], omegaM[t], thetaL[t], omegaL[t] *)
tauS[t_, thetaM_, omegaM_, thetaL_, omegaL_] :=
  Ks*(thetaM/N - thetaL) + Cs*(omegaM/N - omegaL);

tauM[t_, thetaL_, omegaL_] :=
  (Kp*(thetaRef[t] - thetaL) - Kd*omegaL)/N;

eqns = {
  thetaM'[t] == omegaM[t],
  omegaM'[t] == (tauM[t, thetaL[t], omegaL[t]] - bm*omegaM[t] - (tauS[t, thetaM[t], omegaM[t], thetaL[t], omegaL[t]]/N))/Jm,
  thetaL'[t] == omegaL[t],
  omegaL'[t] == (tauS[t, thetaM[t], omegaM[t], thetaL[t], omegaL[t]] - bl*omegaL[t] - tauD[t])/Jl,
  thetaM[0] == 0, omegaM[0] == 0, thetaL[0] == 0, omegaL[0] == 0
};

sol = NDSolve[eqns, {thetaM, omegaM, thetaL, omegaL}, {t, 0, 0.6},
  AccuracyGoal -> 12, PrecisionGoal -> 12
];

thetaLFinal = thetaL[0.6] /. sol[[1]];
peakOmegaL = NMaximize[{Abs[omegaL[t] /. sol[[1]]], 0 <= t <= 0.6}, t][[1]];

thetaLFinal
peakOmegaL
      

Mathematica’s NDSolve provides adaptive integration and is useful for quickly testing parameter regimes (e.g., varying \( K_s \) and observing torsional oscillations).

Problem 1 (Reflected inertia): For a rigid gear with ratio \( \omega_m = N\omega_\ell \), prove that reflecting the load inertia to the motor side yields \( J_{\text{eq},m} = J_m + \frac{J_\ell}{N^2} \).

Solution: Use kinetic energy expressed in motor speed:

\[ T = \tfrac{1}{2}J_m \omega_m^2 + \tfrac{1}{2}J_\ell \omega_\ell^2 = \tfrac{1}{2}J_m \omega_m^2 + \tfrac{1}{2}J_\ell\Big(\tfrac{\omega_m}{N}\Big)^2 = \tfrac{1}{2}\Big(J_m + \tfrac{J_\ell}{N^2}\Big)\omega_m^2. \]

Define \( J_{\text{eq},m} \) by \( T = \tfrac{1}{2}J_{\text{eq},m}\omega_m^2 \). Then \( J_{\text{eq},m} = J_m + \frac{J_\ell}{N^2} \).

Problem 2 (Two-inertia torsional mode): With \( b_m=b_\ell=C_s=0 \) and no inputs, show that \( \phi = \frac{\theta_m}{N} - \theta_\ell \) satisfies \( \ddot{\phi} + K_s(\frac{1}{J_m N^2} + \frac{1}{J_\ell})\phi = 0 \).

Solution: With \( \tau_s = K_s\phi \), the inertial equations reduce to \( J_m \dot{\omega}_m = -\frac{1}{N}\tau_s \), \( J_\ell \dot{\omega}_\ell = \tau_s \). Differentiate \( \dot{\phi} = \frac{\omega_m}{N} - \omega_\ell \):

\[ \ddot{\phi} = \tfrac{1}{N}\dot{\omega}_m - \dot{\omega}_\ell = \tfrac{1}{N}\Big(-\tfrac{1}{J_m N}\tau_s\Big) - \Big(\tfrac{1}{J_\ell}\tau_s\Big) = -\Big(\tfrac{1}{J_m N^2} + \tfrac{1}{J_\ell}\Big)\tau_s. \]

Substitute \( \tau_s = K_s\phi \) to obtain the claimed ODE.

Problem 3 (Rigid-axis PD closed-loop ODE): Assume the rigid load-side model \( J_{\text{eq}}\dot{\omega}_\ell + b_{\text{eq}}\omega_\ell = \tau_\ell \), \( \dot{\theta}_\ell = \omega_\ell \), and control law \( \tau_\ell = K_p(\theta_{\text{ref}}-\theta_\ell) - K_d \omega_\ell \) with constant \( \theta_{\text{ref}} \). Derive the scalar error ODE for \( e = \theta_{\text{ref}}-\theta_\ell \).

Solution: Since \( \dot{e} = -\omega_\ell \) and \( \ddot{e} = -\dot{\omega}_\ell \), we substitute into the plant:

\[ J_{\text{eq}}(-\ddot{e}) + b_{\text{eq}}(-\dot{e}) = K_p e - K_d(-\dot{e}) \;\;\Rightarrow\;\; J_{\text{eq}}\ddot{e} + (b_{\text{eq}}+K_d)\dot{e} + K_p e = 0. \]

Problem 4 (Static error under constant disturbance torque): For the rigid axis with the same PD law, suppose a constant load disturbance \( \tau_d(t)=\tau_{d0} \) enters as \( J_{\text{eq}}\dot{\omega}_\ell + b_{\text{eq}}\omega_\ell = \tau_\ell - \tau_{d0} \). Find the steady-state position error \( e_{\infty} \) assuming the system settles to rest.

Solution: At rest, \( \omega_\ell=0 \), \( \dot{\omega}_\ell=0 \), hence \( 0 = \tau_\ell - \tau_{d0} \). With PD at rest, \( \tau_\ell = K_p e_{\infty} \). Therefore

\[ e_{\infty} = \frac{\tau_{d0}}{K_p}. \]

This highlights why purely proportional position stiffness trades off between disturbance rejection and required torque.

Problem 5 (Numerical experiment design): Using the two-inertia model, explain qualitatively what happens to the load response when you (i) increase \( K_s \) by 10×, (ii) increase \( C_s \) by 10×, holding all else fixed.

Solution: (i) Increasing \( K_s \) raises the torsional frequency \( \omega_t = \sqrt{K_s(\frac{1}{J_m N^2}+\frac{1}{J_\ell})} \), which typically reduces visible oscillation at low frequencies and makes the drivetrain closer to rigid coupling (but may introduce high-frequency ringing if the controller excites that mode). (ii) Increasing \( C_s \) increases dissipation in the relative twist coordinate, reducing torsional oscillation amplitude and speeding decay of twist energy (at the cost of more torque required for the same motion profile).

We developed motor–load dynamics from first principles, proved the reflected inertia and torque relations for ideal gears, derived the two-inertia flexible-shaft ODEs, and obtained a closed-form torsional natural frequency that explains oscillatory behavior in servo axes. Using a torque-source abstraction, we derived the rigid-axis PD error dynamics and proved asymptotic stability via an energy (Lyapunov) argument. Finally, we provided multi-language simulation implementations that directly integrate the governing ODEs.