Continuous Simulation

ODE Simulation

Continuous simulation advances a state vector x(t) governed by ordinary differential equations:

dx/dt = f(t, x)

The state evolves continuously. Numerical integration approximates the solution at discrete time steps.

Explicit Methods

Euler: x(t+h) = x(t) + h·f(t, x(t))

First-order accurate, O(h) local error
Simple but requires small h for stability

RK4 (Classical Runge-Kutta):

Fourth-order accurate, O(h⁴) local error
Workhorse of non-stiff ODE simulation
Four function evaluations per step

STRUCTURE OdeState
    time: real
    values: array of real

PROCEDURE RK4(f, state, h) → OdeState
    n ← LENGTH(state.values)
    t ← state.time
    y ← state.values

    k1 ← f(t, y)
    y2 ← array: FOR i ← 0..n-1: y[i] + 0.5 * h * k1[i]
    k2 ← f(t + 0.5 * h, y2)
    y3 ← array: FOR i ← 0..n-1: y[i] + 0.5 * h * k2[i]
    k3 ← f(t + 0.5 * h, y3)
    y4 ← array: FOR i ← 0..n-1: y[i] + h * k3[i]
    k4 ← f(t + h, y4)

    new_values ← array: FOR i ← 0..n-1:
        y[i] + h/6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])

    RETURN OdeState(time ← t + h, values ← new_values)

Adaptive Step Size

Embedded Runge-Kutta pairs (e.g., Dormand-Prince RK45) compute two estimates of different order. The difference estimates local error:

error ≈ |x₅ - x₄|

If error > tolerance, reject step and reduce h. If error << tolerance, increase h. This automatically uses small steps during rapid transients and large steps during smooth regions.

A system is stiff when it contains dynamics on vastly different time scales — fast transients coexist with slow evolution. Explicit methods must use tiny steps dictated by the fastest mode, even when tracking the slow dynamics.

Example: chemical kinetics where some reactions happen in microseconds and others in minutes.

Implicit Methods

Stiff solvers use implicit formulas requiring the solution of algebraic equations at each step:

Backward Euler: x(t+h) = x(t) + h·f(t+h, x(t+h))

This is a nonlinear equation in x(t+h), solved by Newton iteration. The Jacobian ∂f/∂x is needed.

BDF (Backward Differentiation Formulas): Multi-step implicit methods. BDF1 = backward Euler, BDF2-BDF5 for higher order. Standard in stiff solvers (LSODE, CVODE, SUNDIALS).

Implicit Runge-Kutta: RADAU IIA and Gauss methods offer high order with excellent stability for stiff problems.

Stiffness Detection

Eigenvalues of the Jacobian span many orders of magnitude
Explicit methods require h << 1/|λ_max| for stability, not accuracy
Ratio |λ_max / λ_min| >> 1 indicates stiffness

Differential-Algebraic Equations (DAE)

Many physical systems include algebraic constraints alongside differential equations:

dx/dt = f(t, x, y)
0 = g(t, x, y)

where x are differential variables and y are algebraic variables.

DAE Index

The index measures how many times the algebraic constraint must be differentiated to obtain an ODE system:

Index-0: Already an ODE (trivial)
Index-1: One differentiation — well handled by BDF methods (DASSL, IDA)
Index-2+: Requires index reduction (Pantelides algorithm) or specialized solvers

Applications

Electrical circuits (Kirchhoff's laws impose algebraic constraints)
Multibody mechanics (rigid joint constraints)
Chemical equilibrium coupled with reaction kinetics

Bond Graphs

Bond graphs model physical systems as networks of energy exchange. Every bond carries two variables — effort (e) and flow (f) — whose product is power.

| Domain | Effort | Flow | |--------|--------|------| | Mechanical (translational) | Force | Velocity | | Mechanical (rotational) | Torque | Angular velocity | | Electrical | Voltage | Current | | Hydraulic | Pressure | Volume flow rate | | Thermal | Temperature | Entropy flow |

Core Elements

Sources (Se, Sf): Impose effort or flow
Storage (I, C): Inertance (kinetic energy) and capacitance (potential energy)
Dissipation (R): Resistive element
Junctions (0, 1): 0-junction = common effort, 1-junction = common flow
Transformer (TF), Gyrator (GY): Convert between domains

Bond graphs automatically generate state equations by assigning causality (which variable is input vs output at each bond).

Modelica

Modelica is an equation-based, object-oriented language for multi-domain physical modeling. Instead of writing procedural simulation code, you declare equations and the compiler handles causality assignment and solver selection.

model SpringDamper
  Real x(start=1.0) "displacement";
  Real v(start=0.0) "velocity";
  parameter Real m = 1.0 "mass";
  parameter Real k = 10.0 "spring constant";
  parameter Real c = 0.5 "damping coefficient";
equation
  der(x) = v;
  m * der(v) = -k * x - c * v;
end SpringDamper;

Key Modelica features:

Acausal connectors (equations, not assignments)
Component libraries (mechanical, electrical, thermal, fluid)
Automatic symbolic manipulation (index reduction, tearing)
Tools: OpenModelica (open-source), Dymola, MapleSim

Co-Simulation and FMI

Functional Mock-up Interface (FMI)

FMI is a standard for exchanging simulation models between tools. A Functional Mock-up Unit (FMU) packages:

Model equations (compiled C code or source)
XML description of variables and parameters
Solver (for Model Exchange, the master provides the solver; for Co-Simulation, the FMU includes its own)

Co-Simulation Architecture

┌──────────────┐     ┌──────────────┐
│   FMU A      │     │   FMU B      │
│ (Mechanical) │◄───►│ (Electrical) │
└──────┬───────┘     └──────┬───────┘
       │                    │
       └────────┬───────────┘
           ┌────▼────┐
           │  Master  │
           │Algorithm │
           └──────────┘

The master orchestrates data exchange at communication points. Challenges:

Coupling stability — Explicit exchange can be unstable; implicit or iterative coupling improves stability
Step size coordination — FMUs may require different time steps
Algebraic loops — Direct feedthrough between FMUs requires iteration

Real-Time Simulation

Real-time simulation must complete each integration step within the wall-clock time it represents. Used for:

Pilot/driver training simulators — Visual, motion, and model updates at 60-1000 Hz
Hardware-in-the-loop (HIL) — Real controller hardware connected to simulated plant
Software-in-the-loop (SIL) — Embedded software tested against plant model

Constraints

Fixed time step (no adaptive stepping)
Worst-case execution time must be < step size
Deterministic execution (no garbage collection pauses, no dynamic allocation)
Model fidelity traded against computational budget

Hardware-in-the-Loop (HIL)

┌──────────┐    analog/digital    ┌──────────────┐
│   Real   │◄────────────────────►│  Real-Time   │
│Controller│    I/O interface     │  Simulator   │
│  (ECU)   │                      │ (Plant Model)│
└──────────┘                      └──────────────┘

HIL testing validates embedded controllers against plant models before physical prototypes exist. Critical in automotive (engine ECU, ABS), aerospace (flight control), and power systems (grid protection relays).

Requirements:

Nanosecond-level I/O timing for power electronics HIL
Microsecond-level for automotive powertrain
Millisecond-level for thermal and mechanical systems

Numerical Stability Summary

| Method | Order | Stability Region | Stiff-Suitable | |--------|-------|-------------------|----------------| | Forward Euler | 1 | Small disk | No | | RK4 | 4 | Moderate | No | | Backward Euler | 1 | Entire left half-plane | Yes | | Trapezoidal | 2 | Entire left half-plane | Yes | | BDF-2 | 2 | Large region | Yes | | RADAU IIA | 2s-1 | L-stable | Yes |

For non-stiff problems, explicit RK methods with adaptive stepping are efficient. For stiff problems, implicit methods (BDF, implicit RK) are essential — the cost of solving nonlinear systems at each step is repaid by enabling large stable step sizes.