Differential equations and first- and second-order responses

Convolution tells us how to compute an LTI system's output once we know its impulse response, but where do impulse responses come from? For physical systems, the answer is physics: Newton's laws, Kirchhoff's laws and the motor equations of ROB-220 produce linear ordinary differential equations, and solving them gives the impulse and step responses. This lesson models three components found in every robot, an RC filter, a mass-spring-damper and a DC motor, and shows that their behaviour falls into two families. First-order systems respond with a single exponential characterized by a time constant. Second-order systems can overshoot and oscillate, characterized by a natural frequency and a damping ratio. These few parameters, and the specifications built on them such as rise time, overshoot and settling time, are the everyday language in which robot engineers describe and tune their machines.

From physics to differential equations

An RC low-pass filter, the anti-aliasing filter in front of an ADC (ROB-200), has a resistor R in series and a capacitor C to ground. The capacitor current is \( C\,\dot{v}_{\text{out}} \) and equals the resistor current \( (v_{\text{in}} - v_{\text{out}})/R \), so

A mass m on a spring of stiffness k with viscous damping b, pushed by a force F, obeys Newton's second law, \( m\ddot{x} = F - kx - b\dot{x} \), or

This one equation models a surprising range of robot behaviour: a joint under a proportional-derivative controller (where the controller supplies the spring and damper), the flexible coupling of a series elastic actuator (ROB-220), the bounce of a legged robot's foot, and the wobble of a camera on a mast.

A DC motor (ROB-220) with armature resistance R, inductance L, torque constant \( k_t \), back-EMF constant \( k_e \), rotor inertia J and viscous friction b obeys an electrical and a mechanical equation:

The two equations are coupled through the current and the speed. When the electrical time constant L/R is much shorter than the mechanical one, as it usually is, the inductance can be neglected, \( i = (v - k_e\omega)/R \), and the motor reduces to a single first-order equation for speed:

Back-EMF acts like extra damping, \( k_t k_e/R \), which is usually much larger than mechanical friction. This electrical damping is why a DC motor settles to a steady speed for a given voltage.

First-order systems

Every first-order LTI system can be written in the standard form

with time constant \( \tau \) and steady-state gain K. For the RC filter, \( \tau = RC \) and K = 1; for the simplified motor, \( \tau = J/(b + k_tk_e/R) \) and \( K = (k_t/R)/(b + k_tk_e/R) \). The response to a step of height A from rest is

The time constant has a memorable meaning: after one \( \tau \) the output has covered 63.2 % of its final change, after 2\( \tau \) 86.5 %, after 3\( \tau \) 95.0 % and after 4\( \tau \) 98.2 %. The initial slope is \( KA/\tau \): if the output kept rising at its initial rate it would arrive in exactly one time constant. The 10 to 90 % rise time is \( 2.2\tau \), and the 2 % settling time is about \( 4\tau \). A first-order system never overshoots. Its impulse response is \( (K/\tau)\,e^{-t/\tau} \), as lesson 3 used for the temperature sensor.

Second-order systems

The mass-spring-damper is the prototype second-order system. Dividing by m and introducing two new parameters puts it in standard form:

The undamped natural frequency \( \omega_n \) is the frequency at which the system would oscillate with no damping. The damping ratio \( \zeta \) (zeta) is the ratio of the actual damping to the critical damping \( 2\sqrt{km} \), the smallest damping that prevents oscillation. The damping ratio divides second-order behaviour into regimes:

• Underdamped, \( 0 < \zeta < 1 \): the response overshoots and oscillates with a decaying envelope at the damped frequency \( \omega_d = \omega_n\sqrt{1 - \zeta^2} \).
• Critically damped, \( \zeta = 1 \): the fastest response without overshoot.
• Overdamped, \( \zeta > 1 \): sluggish, no overshoot, behaving like two first-order systems in series.
• Undamped, \( \zeta = 0 \): oscillation forever at \( \omega_n \). Negative damping means instability.

For an underdamped system, the step response from rest is

Its shape is summarized by a few specifications that appear on every control datasheet and in every tuning discussion:

Overshoot depends only on \( \zeta \): 73 % at \( \zeta = 0.1 \), 37 % at 0.3, 16 % at 0.5, 4.6 % at 0.7 and none at 1. A damping ratio around 0.7 is a popular compromise, fast with little overshoot. Settling time depends on \( \zeta\omega_n \), the rate at which the envelope decays. Raising stiffness raises \( \omega_n \) and speeds everything up; adding damping reduces overshoot but, past critical damping, slows the response.

Natural response, forced response and the characteristic equation

Why exponentials? Try \( y = e^{st} \) in the homogeneous equation (input zero). Every derivative brings down a factor s, so \( s^2 + 2\zeta\omega_n s + \omega_n^2 = 0 \), the characteristic equation. Its roots

determine the natural response. Real, negative roots give decaying exponentials (overdamped); a repeated root gives critical damping; complex roots \( -\zeta\omega_n \pm j\omega_d \) give decaying oscillation; roots with positive real part give growth, which is instability. The total response is the natural response, set by initial conditions, plus a forced response that has the same form as the input: a constant for a step, a sinusoid at the input frequency for a sinusoidal input. These roots reappear in lesson 5 as the poles of the transfer function, and their location in the complex plane becomes the main way of thinking about system behaviour.

Simulating responses numerically

Analytical solutions exist for first- and second-order systems, but real robot models quickly get more complicated, so numerical simulation is essential. The simplest method, forward Euler, replaces derivatives by differences over a small step h: \( y_{k+1} = y_k + h\,\dot{y}_k \). The step must be well below the fastest time constant, or the simulation becomes inaccurate or even unstable (lesson 10 explains why). The code below simulates the PD-controlled joint from the example and measures its overshoot and settling time directly.

import numpy as np
J, Kp, Kd = 0.02, 8.0, 0.4
h, T = 1e-4, 1.0
n = int(T/h)
th, w = 0.0, 0.0
traj = np.empty(n)
for k in range(n):
    tau = Kp*(1.0 - th) - Kd*w        # PD torque for a unit step command
    w += h * tau / J
    th += h * w
    traj[k] = th
t = np.arange(n)*h
over = (traj.max() - 1.0)*100
outside = np.where(np.abs(traj - 1.0) > 0.02)[0]
print(f"overshoot      = {over:.1f} %  (formula 16.3 %)")
print(f"peak time      = {t[traj.argmax()]:.3f} s  (formula 0.181 s)")
print(f"2% settling    = {t[outside[-1]]:.3f} s  (rule of thumb 0.4 s)")

overshoot      = 16.3 %  (formula 16.3 %)
peak time      = 0.181 s  (formula 0.181 s)
2% settling    = 0.403 s  (rule of thumb 0.4 s)

The simulated overshoot and peak time match the formulas closely. The measured settling time also lands close to the 4/(ζω_n) rule of thumb, although that rule is only approximate: it follows the decay envelope rather than the actual curve, so for other damping ratios it can be off by 10 to 20 %.