Differential Models

Chapter 1: Velocity Constraints — &ẋ; = f(x, u)

The central equation of Chapter 13 is deceptively simple:

Read this as: "the rate of change of the state at time t equals some function f of the current state and the current input." This is a first-order ordinary differential equation (ODE) — it tells you the direction and speed at which the state moves, not the state itself. Integrating it forward in time gives you a trajectory.

The function f encodes everything about the system's physics. For a point mass sliding on a frictionless surface, f is simple. For a six-degree-of-freedom robot arm with joint limits and gear ratios, f is complex. But the form &ẋ; = f(x, u) is universal — it covers every physical system we'll study in this chapter.

A concrete example: the unicycle

A unicycle has state x = (p_x, p_y, θ) — position and heading angle. Its controls are u = (v, ω) — forward speed and turning rate. The differential equation is:

This says: the unicycle moves in the direction it's currently pointing (cosθ, sinθ) at speed v, and rotates at rate ω. There's no equation saying it can move sideways — that direction simply isn't available. The physics rules it out, not a software constraint.

Worked example: integrating the unicycle

Suppose the unicycle starts at (0, 0, 0) — origin, facing right. Apply constant controls v = 1, ω = π/4 for 2 seconds. What state do we reach?

With constant ω, the heading at time t is θ(t) = ωt = (π/4)t. The position integrates:

At t = 2: θ = π/2, p_x = sin(π/2)/(π/4) = 4/π ≈ 1.27, p_y = (1 − cos(π/2))/(π/4) = 4/π ≈ 1.27. The unicycle has turned 90° and traced a circular arc.

Chapter 2: Nonholonomic Constraints

You want to parallel park. You maneuver — forward, turn, backward, turn — in a sequence of small moves that no single straight path could achieve. You're exploiting something subtle: although your car cannot slide sideways at any instant, it can reach any nearby configuration by combining forward and turning motions cleverly. This is the heart of nonholonomy.

A holonomic constraint restricts the configurations a system can occupy — it reduces the dimension of the configuration space itself. A nonholonomic constraint restricts the velocities available at each configuration, but does not reduce the reachable configurations. The car's no-slip condition is nonholonomic: it says the wheel's velocity must be tangent to the wheel's rolling direction, but with enough maneuvering, you can reach any configuration in the plane.

The no-slip condition for a car

A simplified car (the "kinematic car") has state x = (p_x, p_y, θ) — rear-axle position and heading. The no-slip condition says the rear axle cannot move sideways. Mathematically:

This is a Pfaffian constraint — a linear constraint on the velocity vector (ṗ_x, ṗ_y). It's not integrable to a constraint on position alone, which is precisely what makes it nonholonomic.

Lie brackets and reachability

Here's what's remarkable: even though the car can't move sideways directly, it can move sideways indirectly. The Lie bracket of two vector fields f₁ and f₂ — written [f₁, f₂] — gives the "new direction" generated by alternating between them. For the car, the forward/backward motion vector field f₁ and the turning vector field f₂ have a Lie bracket that points sideways. This is why parallel parking works.

The Lie Algebra Rank Condition (LARC) says: a system is small-time locally controllable (STLC) from x if the set {f₁, f₂, [f₁,f₂], [f₁,[f₁,f₂]], ...} spans the full tangent space at x. For the kinematic car in the plane (3D state space), two vector fields plus one Lie bracket bracket gives three independent directions — full rank — so the car is controllable.

Chapter 3: Phase Space — Position Isn't Enough

Drop a ball. To predict its future trajectory, knowing where it is right now isn't enough. You also need to know how fast it's moving. A ball at position y = 5m moving upward at 3 m/s behaves completely differently from a ball at y = 5m moving downward at 3 m/s — even though they occupy the same geometric position.

For systems governed by Newton's second law (F = ma), the equation of motion is second-order in position: ẍ = f(x, ẋ, u). To convert this into the first-order form &ẋ; = f(x, u) that planning algorithms need, we double the state. The phase space (also called state space in control theory) includes both positions and velocities:

Now the system is first-order in q. The generalized coordinates q₁, ..., q_n describe configuration (positions and angles). Their time derivatives q̇₁, ..., q̇_n are the generalized velocities. Together they form the phase state.

Phase portrait: the geometry of dynamics

In the phase plane (q, q̇), every trajectory of the system traces a curve. The collection of all such curves — one for each initial condition — is the phase portrait. For a simple harmonic oscillator (spring-mass system), the phase portrait consists of concentric ellipses. For a pendulum, it contains closed loops (oscillation) and separatrices (the boundary between oscillation and full rotation).

Chapter 4: Newton-Euler Mechanics

F = ma. You've known this since high school. But what does it mean for a robot with rotating links, joints, and distributed mass? The Newton-Euler formulation extends F = ma to rigid bodies by treating translational and rotational motion separately, then coupling them through constraints at the joints.

Newton's second law: linear momentum

For the center of mass of a rigid body:

Here F is the net external force vector (3D), m is the total mass, and a_cm is the acceleration of the center of mass. This handles translation.

Euler's rotational equation

Rotation is governed by the rotational analogue: torque equals rate of change of angular momentum.

Here τ is the net torque, I is the inertia tensor (a 3×3 matrix encoding how mass is distributed around the center of mass), ω is the angular velocity vector, and α = dω/dt is the angular acceleration. The cross-product term ω × (I·ω) is the gyroscopic term — it appears because we're writing the equation in the body frame, which rotates.

Coupling through joints

A robot arm has n links connected by n joints. Each link satisfies its own Newton-Euler equations, but the forces and torques at each joint connect neighboring links. The Newton-Euler method processes this by recursive outward-inward sweeps:

The result is a set of equations of the form τ = M(q)q̈ + C(q, q̇)q̇ + g(q), where M is the mass matrix, C encodes Coriolis and centrifugal effects, and g is the gravity vector. This is the manipulator equation — the equation of motion for any serial robot arm.

python
import numpy as np

def manipulator_dynamics(q, qdot, tau, M, C_func, g_func):
    """Solve manipulator equation for qddot.

    tau = M(q) @ qddot + C(q,qdot) @ qdot + g(q)
    => qddot = M(q)^{-1} @ (tau - C(q,qdot)@qdot - g(q))

    Args:
        q:      joint angles (n,)
        qdot:   joint velocities (n,)
        tau:    applied joint torques (n,)
        M:      mass matrix (n,n) — symmetric positive definite
        C_func: Coriolis matrix function C(q, qdot) -> (n,n)
        g_func: gravity vector function g(q) -> (n,)
    Returns:
        qddot:  joint accelerations (n,)
    """
    C = C_func(q, qdot)   # Coriolis/centrifugal (n,n)
    g = g_func(q)          # gravity torques (n,)
    rhs = tau - C @ qdot - g
    qddot = np.linalg.solve(M, rhs)  # M^{-1} @ rhs
    return qddot

# Euler integration step
def step(q, qdot, tau, M, C_func, g_func, dt=0.01):
    qddot = manipulator_dynamics(q, qdot, tau, M, C_func, g_func)
    q_new    = q    + qdot  * dt
    qdot_new = qdot + qddot * dt
    return q_new, qdot_new

Chapter 5: Simple Pendulum — Phase Portrait Showcase

The simple pendulum is the most important example in all of classical mechanics. It's nonlinear, exactly solvable in terms of elliptic integrals, and shows every key feature of phase-space dynamics: stable equilibria, unstable equilibria, periodic orbits, and separatrices. It's the sandbox where intuition for differential systems is built.

The equations of motion

A mass m on a massless rod of length L, subject to gravity g. The angle θ is measured from the downward vertical. The torque from gravity is −mgL sinθ (restoring for small angles). The equation of motion is:

Dividing by mL²: θ̈ = −(g/L) sinθ − (b/mL²)θ̇ + u

Define ω₀² = g/L and γ = b/mL². The phase-space form (state = (θ, θ̇)):

Two equilibria: θ = 0 (stable, pendulum hanging down) and θ = π (unstable, pendulum balanced upright). The phase portrait tells you everything about the dynamics at a glance.

Energy interpretation

The total mechanical energy of the pendulum is E = (1/2)mL²θ̇² + mgL(1 − cosθ). With no damping and no control (u = 0, γ = 0), E is conserved — each phase portrait trajectory is a level curve of E. The separatrix corresponds to E = 2mgL, the energy needed to balance the pendulum perfectly upright.

Damping dissipates energy (dE/dt = −γθ̇² ≤ 0), spiraling the trajectory inward toward the stable equilibrium. Control torque injects or removes energy. Swing-up control — bringing the pendulum from hanging to balanced — requires carefully pumping energy to cross the separatrix, then switching to stabilization at the top.

Chapter 6: Lagrangian Mechanics

Newton-Euler is powerful but bookkeeping-heavy. For every link you need to track forces, torques, reaction forces at joints, and carefully propagate them through the chain. For a 7-DOF robot arm, that's seven recursive sweeps with 3D vector arithmetic at each step. There's a better way.

The Lagrangian formulation sidesteps all of that. Instead of tracking forces, you track energy. Define the Lagrangian as kinetic energy minus potential energy:

Here T is the total kinetic energy of the system and V is the total potential energy. Both are scalar functions of the generalized coordinates q and velocities q̇. The Lagrangian captures everything about the system's dynamics in a single scalar — no vectors, no frames, no reaction forces needed.

A systematic recipe

Worked example: single pendulum via Lagrangian

Generalized coordinate: q = θ (angle from downward vertical). The mass m is at position (L sinθ, −L cosθ).

Kinetic energy: velocity magnitude is Lθ̇, so T = (1/2)mL²θ̇².

Potential energy: height is −L cosθ, so V = −mgL cosθ (taking zero at the pivot).

Lagrangian: L = (1/2)mL²θ̇² + mgL cosθ.

Now apply the Euler-Lagrange recipe (next chapter). The result: mL²θ̈ = −mgL sinθ. Divide by mL²: θ̈ = −(g/L) sinθ. Exactly what Newton gives — but we never drew a free-body diagram.

Chapter 7: Euler-Lagrange Equations

We have the Lagrangian L(q, q̇) = T − V. How do we extract equations of motion from it? The answer is the Euler-Lagrange equations — one equation per generalized coordinate:

Here τ_i is the generalized force corresponding to coordinate q_i: a torque if q_i is an angle, a force if q_i is a position. The left side is purely a function of L and its derivatives — no physics beyond energy. The right side is the external input.

Where the Euler-Lagrange equations come from

The Euler-Lagrange equations are the optimality conditions for the action functional: S[q] = ∫_t1^t2 L(q(t), q̇(t)) dt. Hamilton's principle states: the actual trajectory q(t) that a physical system follows is the one that makes S stationary (a local extremum) over all paths with fixed endpoints. Setting δS = 0 and integrating by parts produces exactly the Euler-Lagrange equations. This variational derivation reveals their deep origin — they're not arbitrary; they're the condition for the physical path to be optimal in a precise mathematical sense.

Detailed derivation: double pendulum

Two pendulums linked: link 1 (mass m₁, length L₁, angle θ₁) and link 2 (mass m₂, length L₂, angle θ₂ relative to link 1). The position of mass 2 is:

Its velocity squared: ẋ₂² + ẏ₂² = L₁²θ̇₁² + L₂²(θ̇₁+θ̇₂)² + 2L₁L₂θ̇₁(θ̇₁+θ̇₂)cosθ₂.

Total T = (1/2)m₁L₁²θ̇₁² + (1/2)m₂[above]. Total V = −m₁gL₁cosθ₁ − m₂g[L₁cosθ₁ + L₂cos(θ₁+θ₂)]. Applying Euler-Lagrange to each θ_i yields two coupled nonlinear ODEs — the double pendulum equations. No free-body diagrams. Just differentiation.

python
import numpy as np
from scipy.integrate import solve_ivp

def double_pendulum_rhs(t, y, m1, m2, L1, L2, g=9.81):
    """Euler-Lagrange equations for double pendulum.
    y = [theta1, dtheta1, theta2, dtheta2]
    """
    th1, w1, th2, w2 = y
    d = th2 - th1   # relative angle
    cos_d, sin_d = np.cos(d), np.sin(d)

    # Mass matrix elements (from Euler-Lagrange)
    M11 = (m1 + m2) * L1**2
    M12 = m2 * L1 * L2 * cos_d
    M21 = M12
    M22 = m2 * L2**2

    # Right-hand side (gravity + Coriolis)
    rhs1 = -m2 * L1 * L2 * w2**2 * sin_d - (m1+m2)*g*L1*np.sin(th1)
    rhs2 =  m2 * L1 * L2 * w1**2 * sin_d - m2*g*L2*np.sin(th2)

    # Solve 2x2 linear system for angular accelerations
    M = np.array([[M11, M12],[M21, M22]])
    rhs = np.array([rhs1, rhs2])
    alpha = np.linalg.solve(M, rhs)

    return [w1, alpha[0], w2, alpha[1]]

# Integrate for 5 seconds from near-vertical start
sol = solve_ivp(double_pendulum_rhs, [0, 5], [3.0, 0, 0.5, 0],
                args=(1, 1, 1, 1), max_step=0.01)

Chapter 8: Multiple Decision Makers — Differential Games

Everything so far assumed one player: a single robot with a single control input. But what about two robots competing for a target? A drone trying to evade a pursuer? An autonomous car sharing a road with a human driver who has their own objectives? These are differential games — dynamic systems with multiple decision makers, each optimizing their own objective.

A two-player differential game combines two control inputs in a single differential equation:

Player 1 (minimizer) chooses u₁ ∈ U₁. Player 2 (maximizer) chooses u₂ ∈ U₂. The classic pursuit-evasion game has the pursuer minimizing time-to-capture while the evader maximizes it. Neither player can observe the other's control directly — they must reason about each other's best strategies.

The pursuit-evasion game on the plane

Pursuer P at position p ∈ ℝ², evader E at position e ∈ ℝ². Both move at fixed speeds: ‖ṗ‖ = v_P, ‖ė‖ = v_E. The state is x = p − e (relative position). The game ends when ‖x‖ ≤ r (capture radius).

The key quantity is the value function V*(x) — the optimal time-to-capture when both players play optimally. It satisfies the Hamilton-Jacobi-Isaacs (HJI) equation:

The −1 on the right comes from the fact that value decreases by 1 per second of real time. Solving the HJI equation gives V* and the saddle-point strategies: the pursuer steers to minimize V (head toward capture), the evader steers to maximize V (flee optimally).

Isaacs' result: the homicidal chauffeur

A famous differential game (Isaacs, 1965): the pursuer is a car with minimum turning radius r_min (think a fast but not very maneuverable vehicle). The evader is a slower pedestrian with no turning constraint. Who wins? It depends on the geometry. When the pursuer is fast enough and the initial configuration is in the "capture zone," the pursuer always wins regardless of evasion strategy — this zone boundary is called the barrier or safe evader set.

Nash equilibria in general-sum games

Not all multi-player games are zero-sum. In a general-sum game, each player has a distinct cost functional. A Nash equilibrium is a strategy profile where no player benefits by unilaterally switching strategies. For two players with cost functionals J₁, J₂:

Finding Nash equilibria for differential games requires solving a coupled system of HJI equations — one per player. This is significantly harder than the single-player (optimal control) case and remains an active research area.

Chapter 9: Connections and What Comes Next

You've now seen the full differential-model toolkit. Let's place it in context: what came before, what comes after, and where this mathematics appears beyond robotics.

What we built in this chapter

The bridge to Chapter 14

Chapter 14 (Sampling-Based Planning with Differential Constraints) takes everything from this chapter and asks: how do we plan trajectories for systems governed by &ẋ; = f(x, u)? The key insight is that random-tree methods (RRT, EST) extend naturally: instead of sampling configurations and connecting them with straight lines, you sample controls u and integrate &ẋ; = f(x, u) forward to generate edges. The differential model you learned here becomes the edge generator for the planner.

Where differential models appear beyond robotics

Limitations of the differential-model framework

Model accuracy. The manipulator equation assumes rigid links and perfect actuators. Real robots flex, backlash exists, and motors have nonlinear saturation. Model errors propagate through integration, growing over time. Robust and adaptive control methods address this.

Contact discontinuities. When a robot foot hits the ground, the dynamics change abruptly — a hybrid system with discrete mode switches. Pure differential equations don't capture this; hybrid systems theory extends &ẋ; = f(x, u) with discrete guards and reset maps.

Integration accuracy. Euler integration (Δx ≈ f(x,u)Δt) accumulates errors. Real implementations use Runge-Kutta 4 (RK4) for ~10× better accuracy at similar cost. ODE solver libraries (scipy.integrate.solve_ivp) handle stiff systems automatically.

Related micro-lessons