Many actuators and motion axes in control and robotics (robot joints, DC motor drives, compliant links, camera gimbals) are well approximated by a single dominant second-order mode. At moderate bandwidths and within a linear operating region, their input–output dynamics can be described by a scalar output \( y(t) \) and scalar input \( u(t) \) obeying a constant-coefficient second-order linear differential equation

\[ a_2 \ddot{y}(t) + a_1 \dot{y}(t) + a_0 y(t) = b_0 u(t), \quad a_2 \neq 0. \]

Here \( a_2, a_1, a_0, b_0 \) are real constants determined by the physical parameters (mass, damping, stiffness, inductance, resistance, etc.). The goal of this lesson is to show that any such model can be written in the standard second-order form

\[ \ddot{y}(t) + 2 \zeta \omega_n \dot{y}(t) + \omega_n^2 y(t) = K \,\omega_n^2 u(t), \]

where:

• \( \omega_n > 0 \) is the undamped natural frequency,
• \( \zeta \) is the damping ratio, and
• \( K \) is the low-frequency (static) gain.

These three parameters \( (\omega_n, \zeta, K) \) fully determine the step response and transient characteristics, which will be studied in later lessons.

flowchart TD
  P["Physical system (mass-spring, DC motor, robot joint)"]
    --> D["Second-order ODE a2 y'' + a1 y' + a0 y = b0 u"]
  D --> N["Normalize coefficients"]
  N --> S["Standard parameters (omega_n, zeta, K)"]
  S --> G["Transfer function G(s)"]
  G --> T["Time-response metrics (rise time, overshoot, settling)"]
        

Consider again the general second-order LTI model with scalar input and output:

\[ a_2 \ddot{y}(t) + a_1 \dot{y}(t) + a_0 y(t) = b_0 u(t), \quad a_2 \neq 0. \]

We can always normalize the leading coefficient by dividing both sides by \( a_2 \):

\[ \ddot{y}(t) + \frac{a_1}{a_2}\dot{y}(t) + \frac{a_0}{a_2} y(t) = \frac{b_0}{a_2} u(t). \]

We now define the standard second-order parameters:

\[ \omega_n := \sqrt{\frac{a_0}{a_2}}, \qquad 2 \zeta \omega_n := \frac{a_1}{a_2}, \qquad K := \frac{b_0}{a_0}. \]

Substituting these definitions into the normalized equation gives

\[ \ddot{y}(t) + 2 \zeta \omega_n \dot{y}(t) + \omega_n^2 y(t) = K \,\omega_n^2 u(t), \]

because \( \frac{a_1}{a_2} = 2 \zeta \omega_n \), \( \frac{a_0}{a_2} = \omega_n^2 \), and \( \frac{b_0}{a_2} = \frac{b_0}{a_0}\frac{a_0}{a_2} = K \omega_n^2 \). This proves that any scalar second-order LTI model with \( a_2 \neq 0 \) can be written in the standard form with uniquely determined \( \omega_n > 0 \), \( \zeta \), and \( K \).

Note that:

• \( \omega_n \) depends only on the ratio \( a_0 / a_2 \) and has units of rad/s.
• \( \zeta \) is dimensionless and involves both \( a_1 \) and \( a_0, a_2 \).
• \( K \) reflects the steady-state gain of the system, provided the closed-loop configuration is appropriate and the system is stable.

A standard mechanical example is the mass–spring–damper system used extensively in robotics to model compliant joints or end-effector dynamics. Let \( x(t) \) be the displacement of a mass \( m \), subject to viscous damping \( c \) and spring stiffness \( k \), with applied force \( f(t) \). Newton's law yields

\[ m \ddot{x}(t) + c \dot{x}(t) + k x(t) = f(t). \]

Identifying \( y(t) = x(t) \) and \( u(t) = f(t) \), the coefficients are \( a_2 = m \), \( a_1 = c \), \( a_0 = k \), \( b_0 = 1 \). Thus

\[ \omega_n = \sqrt{\frac{k}{m}}, \qquad \zeta = \frac{c}{2\sqrt{km}}, \qquad K = \frac{b_0}{a_0} = \frac{1}{k}. \]

Substituting these into the standard form gives

\[ \ddot{x}(t) + 2 \zeta \omega_n \dot{x}(t) + \omega_n^2 x(t) = K \,\omega_n^2 f(t) = \frac{\omega_n^2}{k} f(t) = \frac{1}{m} f(t), \]

so the original mass–spring–damper model is indeed a special case of the standard second-order system. In robot joint modeling, \( x(t) \) can represent a link angle, the spring stiffness can represent joint elasticity, and \( f(t) \) the torque generated by a motor plus transmission.

For control design we usually work with the Laplace-domain transfer function. Taking the Laplace transform of the standard form under zero initial conditions, with \( Y(s) \) and \( U(s) \) denoting the Laplace transforms of \( y(t) \) and \( u(t) \), we obtain

\[ \begin{aligned} &\mathcal{L}\{\ddot{y}(t)\} + 2 \zeta \omega_n \mathcal{L}\{\dot{y}(t)\} + \omega_n^2 \mathcal{L}\{y(t)\} \\ &= K \,\omega_n^2 \mathcal{L}\{u(t)\}. \end{aligned} \]

Zero initial conditions give \( \mathcal{L}\{\ddot{y}(t)\} = s^2 Y(s) \), \( \mathcal{L}\{\dot{y}(t)\} = s Y(s) \), \( \mathcal{L}\{y(t)\} = Y(s) \), so

\[ \bigl(s^2 + 2 \zeta \omega_n s + \omega_n^2\bigr) Y(s) = K \,\omega_n^2 U(s). \]

The transfer function from input \( U(s) \) to output \( Y(s) \) is therefore

\[ G(s) := \frac{Y(s)}{U(s)} = \frac{K \,\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}. \]

The denominator polynomial \( s^2 + 2 \zeta \omega_n s + \omega_n^2 \) is called the characteristic polynomial of the system and fully determines the homogeneous response. The numerator \( K \omega_n^2 \) determines the low-frequency gain and scaling of the forced response.

Conversely, any strictly proper second-order transfer function with real coefficients

\[ G(s) = \frac{b_0}{a_2 s^2 + a_1 s + a_0} \]

can be written in the standard form by factoring out \( a_2 \) from the denominator and applying the definitions of \( \omega_n, \zeta, K \) given earlier.

The homogeneous dynamics of the standard second-order system are obtained by setting \( u(t) = 0 \):

\[ \ddot{y}(t) + 2 \zeta \omega_n \dot{y}(t) + \omega_n^2 y(t) = 0. \]

Seeking exponential solutions of the form \( y(t) = e^{s t} \) leads to the characteristic equation

\[ s^2 + 2 \zeta \omega_n s + \omega_n^2 = 0. \]

Applying the quadratic formula,

\[ s_{1,2} = -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1}. \]

These roots \( s_1, s_2 \) are exactly the poles of \( G(s) \). Depending on the value of \( \zeta \), the roots are:

• real and distinct if \( \zeta^2 - 1 > 0 \),
• real and repeated if \( \zeta^2 - 1 = 0 \),
• complex conjugate if \( \zeta^2 - 1 < 0 \).

The general homogeneous solution is

\[ y_h(t) = \begin{cases} C_1 e^{s_1 t} + C_2 e^{s_2 t}, & \text{if } s_1 \neq s_2, \\ (C_1 + C_2 t)e^{s_1 t}, & \text{if } s_1 = s_2, \end{cases} \]

where \( C_1, C_2 \) depend on the initial position and velocity. In the next lesson we will analyze how the value of \( \zeta \) controls oscillations, overshoot, and settling behavior of this solution when combined with particular inputs such as step or ramp.

In robotics, a single revolute joint driven by a motor and gear train can often be approximated by a second-order model in joint angle. A typical workflow is:

1. Identify physical parameters (equivalent inertia, damping, stiffness).
2. Form the torque balance equation as a second-order ODE.
3. Normalize to standard form and compute \( \omega_n, \zeta, K \).
4. Build a transfer function model for analysis and simulation.
5. Embed the model in a controller design toolchain (MATLAB, python-control, ROS control, WPILib, etc.).

flowchart TD
  A["Robot joint parameters (inertia, damping, stiffness)"]
    --> B["Torque balance: J theta'' + B theta' + K theta = tau"]
  B --> C["Normalize coefficients"]
  C --> D["Compute omega_n, zeta, K"]
  D --> E["Construct G(s)"]
  E --> F["Use in control design & simulation"]
        

In Python, the python-control library (often combined with robotics stacks such as ROS for higher-level integration) provides transfer function and simulation primitives suitable for second-order systems.

import numpy as np
import matplotlib.pyplot as plt

# python-control: pip install control
import control as ctrl

# Standard second-order parameters
wn = 4.0      # natural frequency (rad/s)
zeta = 0.7    # damping ratio
K = 1.0       # DC gain

# Transfer function G(s) = K*wn^2 / (s^2 + 2*zeta*wn*s + wn^2)
num = [K * wn**2]
den = [1.0, 2.0 * zeta * wn, wn**2]
G = ctrl.TransferFunction(num, den)

# Step response
t, y = ctrl.step_response(G)

plt.figure()
plt.plot(t, y)
plt.xlabel("time (s)")
plt.ylabel("output y(t)")
plt.title("Standard second-order step response (Python)")
plt.grid(True)
plt.show()

# Remark:
# In robotics, G can model a single joint; G is then combined with
# trajectory generators and controllers (e.g. PID, state feedback)
# within a ROS or ROS2 control stack.
      

For more explicit state integration (e.g. embedding in a custom simulator), one can write the second-order ODE as a first-order system in velocity and integrate with scipy.integrate.solve_ivp, using the same \( \omega_n, \zeta, K \) parameters.

In C++, second-order joint models are often integrated in real time inside a robot controller loop (e.g. within ros_control or similar frameworks). Below is a simple explicit Euler integrator for the standard second-order ODE with unit step input \( u(t) = 1 \).

#include <iostream>
#include <vector>

// Simple numerical simulation of y'' + 2*zeta*wn*y' + wn^2*y = K*wn^2*u
int main() {
    double wn = 4.0;
    double zeta = 0.7;
    double K = 1.0;

    double dt = 0.001;        // integration step (s)
    int N = 5000;             // number of steps (5 seconds)

    double y = 0.0;           // position
    double ydot = 0.0;        // velocity
    double u = 1.0;           // unit step input

    for (int k = 0; k < N; ++k) {
        double t = k * dt;

        // Second-order dynamics
        double yddot =
            -2.0 * zeta * wn * ydot
            - wn * wn * y
            + K * wn * wn * u;

        // Explicit Euler integration
        ydot += dt * yddot;
        y    += dt * ydot;

        // Log a coarse sample for plotting in an external tool
        if (k % 1000 == 0) {
            std::cout << t << " " << y << std::endl;
        }
    }

    return 0;
}

// Note: In a robotics application, this loop would typically run inside a
// real-time control thread, and the state y, ydot would be coupled with
// joint-level feedback, possibly using Eigen for vectorized multi-joint
// models and ros_control for interfacing with the hardware.
      

This pattern extends naturally to vector systems \( \mathbf{y}(t) \) using matrix-based formulations with the Eigen library in multi-degree-of-freedom robots.

In educational and industrial robotics that use Java (for example, FRC robots), the WPILib library provides linear system utilities that can represent standard second-order plants in state-space. Here we show both a bare-bones numerical simulation and a sketch of how WPILib would be used.

public class SecondOrderSimulation {

    // Standard second-order parameters
    private final double wn;    // natural frequency
    private final double zeta;  // damping ratio
    private final double K;     // DC gain

    private double y;           // output position
    private double ydot;        // output velocity

    public SecondOrderSimulation(double wn, double zeta, double K) {
        this.wn = wn;
        this.zeta = zeta;
        this.K = K;
        this.y = 0.0;
        this.ydot = 0.0;
    }

    // Single simulation step with explicit Euler
    public void step(double u, double dt) {
        double yddot =
            -2.0 * zeta * wn * ydot
            - wn * wn * y
            + K * wn * wn * u;

        ydot += dt * yddot;
        y    += dt * ydot;
    }

    public double getY() {
        return y;
    }

    public static void main(String[] args) {
        SecondOrderSimulation sim =
            new SecondOrderSimulation(4.0, 0.7, 1.0);

        double dt = 0.001;
        double u = 1.0; // step input

        for (int k = 0; k < 5000; ++k) {
            sim.step(u, dt);
            if (k % 1000 == 0) {
                double t = k * dt;
                System.out.println(t + " " + sim.getY());
            }
        }
    }
}
      

In WPILib, one can instead instantiate a second-order plant via predefined utilities (e.g., for arms or elevators) and simulate it within the edu.wpi.first.math.system package:

// Sketch (requires appropriate WPILib dependencies):
//
// import edu.wpi.first.math.system.LinearSystem;
// import edu.wpi.first.math.system.plant.LinearSystemId;
// import edu.wpi.first.math.numbers.N1;
//
// LinearSystem<N1, N1, N1> armPlant =
//     LinearSystemId.createSingleJointedArmSystem(
//         /* armMomentOfInertia */ J,
//         /* armLength */ L,
//         /* armMass */ m,
//         /* motor */ motor,        // motor model
//         /* gearing */ G);
//
// // This armPlant behaves approximately as a second-order system in arm angle.
// // WPILib then provides simulation and controller classes built on top of it.
      

Such abstractions hide the explicit ODE but internally rely on the same second-order structure discussed in this lesson.

MATLAB's Control System Toolbox and Simulink provide built-in tools for defining and simulating standard second-order systems. Below is a script to define the transfer function and plot its step response.

wn   = 4.0;   % natural frequency (rad/s)
zeta = 0.7;   % damping ratio
K    = 1.0;   % DC gain

num = K * wn^2;
den = [1, 2*zeta*wn, wn^2];

G = tf(num, den);

figure;
step(G);
grid on;
title('Standard second-order step response (MATLAB)');
xlabel('time (s)');
ylabel('output y(t)');

% In Robotics System Toolbox, similar transfer functions are used as
% components inside manipulator and mobile robot controllers.
      

The same system can be created in Simulink programmatically:

model = 'second_order_model';
new_system(model);
open_system(model);

add_block('simulink/Sources/Step',        [model '/Step']);
add_block('simulink/Continuous/Transfer Fcn', [model '/G']);
add_block('simulink/Sinks/Scope',         [model '/Scope']);

set_param([model '/G'], 'Numerator',   'K*wn^2', ...
                         'Denominator', '[1 2*zeta*wn wn^2]');

add_line(model, 'Step/1', 'G/1');
add_line(model, 'G/1',    'Scope/1');

set_param(model, 'StopTime', '5');
sim(model);
      

In robotics, such Simulink models are frequently embedded in more complex diagrams including multibody dynamics, sensors, and controllers.

Wolfram Mathematica offers symbolic and numeric tools for analyzing second-order systems. The code below constructs the transfer function and computes the response to a unit step input.

wn   = 4.0;
zeta = 0.7;
K    = 1.0;

s = LaplaceTransformVariable["s"];

G = TransferFunctionModel[
      K*wn^2/(s^2 + 2*zeta*wn*s + wn^2),
      s
    ];

y[t_] = OutputResponse[G, UnitStep[t], {t, 0, 5}];

Plot[y[t], {t, 0, 5},
     AxesLabel -> {"t", "y(t)"},
     PlotLabel -> "Standard second-order step response (Mathematica)"]
      

For robotic applications, Mathematica can symbolically derive second-order models from multibody dynamics and then export them to code or controller design tools.

Problem 1 (Identification of Standard Parameters): A second-order system is described by \( 0.5 \ddot{y}(t) + 1.0 \dot{y}(t) + 8 y(t) = 4 u(t) \). Express this system in standard form and determine \( \omega_n \), \( \zeta \), and \( K \).

Solution: Here \( a_2 = 0.5 \), \( a_1 = 1.0 \), \( a_0 = 8 \), \( b_0 = 4 \). Then

\[ \omega_n = \sqrt{\frac{a_0}{a_2}} = \sqrt{\frac{8}{0.5}} = \sqrt{16} = 4 \text{ rad/s}, \]

\[ 2 \zeta \omega_n = \frac{a_1}{a_2} = \frac{1.0}{0.5} = 2 \quad \Rightarrow \quad \zeta = \frac{2}{2 \cdot 4} = \frac{1}{4} = 0.25, \]

\[ K = \frac{b_0}{a_0} = \frac{4}{8} = 0.5. \]

The standard form is

\[ \ddot{y}(t) + 2 (0.25)(4) \dot{y}(t) + 4^2 y(t) = 0.5 \cdot 4^2 u(t), \]

or simply \( \ddot{y}(t) + 2 \cdot 0.25 \cdot 4 \dot{y}(t) + 16 y(t) = 8 u(t) \), which matches the original equation after multiplying through by \( 0.5 \).

Problem 2 (Mass–Spring–Damper Parameters): For a robotic joint modeled as a mass–spring–damper with \( m = 2 \,\mathrm{kg} \), \( c = 1 \,\mathrm{N\,s/m} \), \( k = 18 \,\mathrm{N/m} \), find \( \omega_n \), \( \zeta \), and the static gain \( K \) from torque input to displacement output.

Solution: Using the formulas \( \omega_n = \sqrt{k/m} \), \( \zeta = c / (2 \sqrt{km}) \), \( K = 1/k \), we obtain

\[ \omega_n = \sqrt{\frac{18}{2}} = \sqrt{9} = 3 \text{ rad/s}, \]

\[ \zeta = \frac{c}{2 \sqrt{km}} = \frac{1}{2 \sqrt{18 \cdot 2}} = \frac{1}{2 \sqrt{36}} = \frac{1}{12} \approx 0.0833, \]

\[ K = \frac{1}{k} = \frac{1}{18} \approx 0.0556 \text{ m/N}. \]

The joint is lightly damped (small \( \zeta \)), with natural frequency \( 3 \,\mathrm{rad/s} \).

Problem 3 (Poles of the Standard Second-Order Transfer Function): Consider the standard transfer function \( G(s) = \dfrac{K \omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2} \). Show that the poles are given by \( s_{1,2} = -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1} \).

Solution: The poles are the roots of the denominator \( s^2 + 2 \zeta \omega_n s + \omega_n^2 = 0 \). Applying the quadratic formula,

\[ s_{1,2} = \frac{-2 \zeta \omega_n \pm \sqrt{(2 \zeta \omega_n)^2 - 4 \cdot 1 \cdot \omega_n^2}}{2}. \]

Simplify the discriminant:

\[ (2 \zeta \omega_n)^2 - 4 \omega_n^2 = 4 \omega_n^2 (\zeta^2 - 1). \]

Therefore

\[ s_{1,2} = \frac{-2 \zeta \omega_n \pm \sqrt{4 \omega_n^2 (\zeta^2 - 1)}}{2} = -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1}, \]

which proves the desired expression.

Problem 4 (Stability Condition from Coefficients): Consider the general second-order transfer function \( G(s) = \dfrac{b_0}{a_2 s^2 + a_1 s + a_0} \) with \( a_2 > 0 \). Give conditions on \( a_1 \) and \( a_0 \) that guarantee both poles have negative real parts (asymptotic stability), and express these conditions in terms of \( \omega_n \) and \( \zeta \).

Solution: For a second-order polynomial \( a_2 s^2 + a_1 s + a_0 \) with \( a_2 > 0 \), classical results (or the Routh–Hurwitz criterion) show that asymptotic stability is guaranteed if

• \( a_1 > 0 \),
• \( a_0 > 0 \).

Indeed, if \( a_0 < 0 \) or \( a_1 < 0 \), there is at least one root with nonnegative real part. In terms of standard parameters,

\[ a_0 = a_2 \omega_n^2, \qquad a_1 = 2 a_2 \zeta \omega_n. \]

With \( a_2 > 0 \) and \( \omega_n^2 > 0 \), we automatically have \( a_0 > 0 \). The condition \( a_1 > 0 \) becomes \( 2 a_2 \zeta \omega_n > 0 \), which implies \( \zeta > 0 \), since \( a_2 > 0 \) and \( \omega_n > 0 \). Thus, in the standard second-order parameterization, asymptotic stability is equivalent to \( \omega_n > 0 \) and \( \zeta > 0 \).

In this lesson we established the standard second-order system model as the canonical representation of many physical and robotic actuators. Starting from a general second-order LTI differential equation, we introduced the natural frequency \( \omega_n \), damping ratio \( \zeta \), and static gain \( K \), and showed how these parameters determine both the differential equation and the transfer function denominator. We derived the characteristic equation, expressed the poles in closed form, and connected these ideas to numerical simulation across several programming languages and robotics-oriented software stacks. Subsequent lessons will use this model to relate \( \omega_n \) and \( \zeta \) to time-domain performance measures such as rise time, overshoot, and settling time.