Skip to content

Inertial Kuramoto — Second-Order Swing Equation

The InertialKuramotoEngine extends the Kuramoto model with rotational inertia and damping — the swing equation from power systems engineering. Each oscillator has both phase \(\theta_i\) and frequency deviation \(\dot{\omega}_i\) as state variables, making the system second-order.

This is the canonical model for: - Power grid stability — synchronous generators with mechanical inertia - Rotating machinery — coupled mechanical oscillators - Neural oscillators with adaptation — slow frequency adaptation

1. Mathematical Formalism

1.1 The Swing Equation (Second-Order Kuramoto)

\[ M_i \ddot{\theta}_i = P_i - D_i \dot{\theta}_i + \sum_{j=1}^{N} K_{ij} \sin(\theta_j - \theta_i) \]

Rewritten as a first-order system:

\[ \dot{\theta}_i = \dot{\omega}_i \]
\[ M_i \dot{\dot{\omega}}_i = P_i - D_i \dot{\omega}_i + \sum_{j=1}^{N} K_{ij} \sin(\theta_j - \theta_i) \]

where:

Symbol Description Units Power grid analogue
\(\theta_i\) Rotor angle radians Generator electrical angle
\(\dot{\omega}_i\) Angular velocity deviation rad/s Frequency deviation from nominal
\(M_i\) Inertia constant s² \(H_i / (\pi f_0)\) (generator H-constant)
\(D_i\) Damping coefficient s⁻¹ Governor droop characteristic
\(P_i\) Mechanical power input dimensionless \(P_{\text{mech}} - P_{\text{elec}}\) mismatch
\(K_{ij}\) Coupling strength dimensionless Susceptance \(B_{ij} V_i V_j\)

1.2 RK4 Integration

The engine uses classical fourth-order Runge-Kutta for both \(\theta\) and \(\dot{\omega}\) simultaneously. Each stage evaluates:

\[ \dot{\theta} = \dot{\omega}, \quad \ddot{\omega} = \frac{P + K\sin(\cdot) - D\dot{\omega}}{M} \]

Four stages with weights \((1/6, 1/3, 1/3, 1/6)\), advancing both state variables in lockstep.

1.3 Relationship to Standard Kuramoto

Setting \(M_i = 0\) (zero inertia) eliminates the acceleration term, and the system reduces to:

\[ D_i \dot{\theta}_i = P_i + \sum_j K_{ij} \sin(\theta_j - \theta_i) \]

which is the standard (overdamped) Kuramoto model with \(\omega_i = P_i / D_i\). The inertial model adds oscillatory transients and the possibility of desynchronisation through oscillatory instability (not just saddle-node bifurcation as in first-order).

1.4 Critical Coupling for Inertial Systems

For identical inertia \(M\) and damping \(D\) with all-to-all coupling:

\[ K_c = \frac{2\Delta}{\cos(\arctan(M \cdot \Omega_{\text{nat}} / D))} \]

where \(\Delta\) is the frequency spread and \(\Omega_{\text{nat}}\) is the natural oscillation frequency \(\sqrt{K/M}\). Key insight: inertia raises \(K_c\) — heavier oscillators are harder to synchronise.

1.5 Frequency Deviation and Grid Stability

In power grids, the frequency deviation \(\dot{\omega}_i\) directly maps to the grid frequency error:

\[ f_i(t) = f_0 + \frac{\dot{\omega}_i(t)}{2\pi} \]

where \(f_0 = 50\) Hz (Europe) or 60 Hz (Americas). The frequency_deviation() method returns \(\max_i |\dot{\omega}_i| / 2\pi\) — the worst-case frequency error in Hz.

Grid stability requires \(\max_i |\Delta f_i| < 0.5\) Hz for normal operation and \(< 2.5\) Hz before load shedding triggers.

1.6 Oscillatory Instability

Unlike first-order Kuramoto (which only has monotone convergence or divergence), the inertial model exhibits oscillatory instability: frequency oscillations grow in amplitude before the system desynchronises. This is the dominant failure mode in power grids — "inter-area oscillations" between generator clusters.

The damping ratio \(\zeta_d = D / (2\sqrt{KM})\) determines the character: - \(\zeta_d > 1\): overdamped (first-order-like, no oscillations) - \(0 < \zeta_d < 1\): underdamped (oscillatory approach to sync) - \(\zeta_d = 0\): undamped (perpetual oscillations)

2. Theoretical Context

2.1 Historical Background

The swing equation was developed by Park (1929) and Concordia (1951) for analysing synchronous generator transient stability. Filatrella, Nielsen & Pedersen (2008) connected it to the Kuramoto model, showing that a network of coupled generators is mathematically equivalent to a second-order Kuramoto model with inertia.

Dörfler & Bullo (2012, 2014) provided rigorous synchronisation conditions for the swing equation on general graphs, establishing the relationship between graph structure, inertia distribution, and critical coupling.

Rohden et al. (2012) showed that decentralised renewable energy sources (with lower inertia than conventional generators) can destabilise power grids by reducing the effective inertia constant.

2.2 Role in SCPN

The inertial engine models SCPN Layer 6 (Planetary/geophysical oscillations) and Layer 10 (Boundary Control), where oscillators have physical mass or momentum that resists instantaneous phase changes. It is also used for power grid digital twin simulations via the plasma_control_bridge adapter.

2.3 Modal Analysis

For small perturbations around the synchronised state \(\theta_i = \theta_0\), \(\dot{\omega}_i = 0\), linearisation gives:

\[ M_i \ddot{\delta}_i + D_i \dot{\delta}_i + \sum_j L_{ij} \delta_j = 0 \]

where \(L_{ij} = K_{ij} \cos(\theta_j^0 - \theta_i^0)\) is the linearised coupling (Laplacian-like). This is a damped harmonic system with \(N\) coupled modes.

The eigenvalues of the linearised system \(\lambda_k\) determine stability: - All \(\text{Re}(\lambda_k) < 0\): stable (all modes damped) - Any \(\text{Re}(\lambda_k) > 0\): oscillatory instability - \(\text{Re}(\lambda_k) = 0\): marginal (Hopf bifurcation possible)

The slowest-decaying mode (smallest \(|\text{Re}(\lambda)|\)) is the critical mode — it determines the transient settling time and is typically an inter-area oscillation between generator clusters.

2.4 Energy Function

For uniform damping \(D\), the system has a Lyapunov-like energy function:

\[ V = \sum_i \frac{1}{2} M_i \dot{\omega}_i^2 - \sum_{i<j} K_{ij} \cos(\theta_j - \theta_i) + \sum_i P_i (\theta_i - \theta_i^*) \]

where \(\theta^*\) is the stable equilibrium. \(\dot{V} \leq 0\) (energy dissipated by damping), proving asymptotic stability when \(V\) has a minimum at the sync state. The run() trajectory can be used to verify \(V(t) \to V_{\min}\) numerically.

2.5 Applications

Application Typical N Inertia M Damping D Power P
Power grid (EU) 10–1000 3–8 s 0.5–2.0 s⁻¹ ±0.1–1.0
Mechanical rotors 2–20 0.1–10 kg·m² 0.01–1.0 N·m·s torque
Neural adaptation 16–256 0.01–0.1 1.0–10.0 firing rate
SCPN Layer 6 16 1.0 0.5 0.0

3. Pipeline Position

┌──────────────┐     ┌──────────────────────┐     ┌──────────────┐
│ coupling/    │────→│ InertialKuramoto    │────→│ monitor/     │
│ knm.py       │     │ Engine              │     │ order_params │
│ (K_ij)       │     │                     │     │ lyapunov     │
└──────────────┘     │ State: (θ, ω̇)      │     └──────────────┘
                     │ RK4 integration     │
┌──────────────┐     │                     │     ┌──────────────┐
│ oscillators/ │────→│ power P_i           │     │ supervisor/  │
│ (P, M, D)    │     │ inertia M_i         │────→│ regimes.py   │
└──────────────┘     │ damping D_i         │     │ (freq check) │
                     └──────────────────────┘     └──────────────┘

Inputs: - theta (N,) — rotor angles in \([0, 2\pi)\) - omega_dot (N,) — angular velocity deviations - power (N,) — mechanical power input \(P_i\) - knm (N,N) — coupling matrix - inertia (N,) — inertia constants \(M_i\) - damping (N,) — damping coefficients \(D_i\)

Outputs: - (theta, omega_dot) — updated state after RK4 step - run() also returns trajectory arrays

Direct Go, Julia, and Mojo accelerator wrappers validate the swing-equation state before loading optional runtimes. The direct boundary requires finite real one-dimensional float64 buffers for theta, omega_dot, power, inertia, and damping; positive inertia, damping, and timestep values; an exactly N*N flat pairwise coupling buffer; and a zero self-coupling diagonal. Numeric-string aliases are rejected before these values are widened by Python, NumPy, or accelerator boundary casts. Backend outputs must return finite (theta, omega_dot) vectors of length N, with phases on the torus [0, 2π), and backend output aliases are rejected before publication.

The public InertialKuramotoEngine.step() dispatcher and Rust wrapper replay that same output contract before returning optional backend state. Loader and runtime unavailability still fall back through the backend chain, but malformed backend physics payloads raise instead of becoming simulation state. The public step(), run(), frequency_deviation(), and coherence() entrypoints apply the same pre-coercion numeric-string rejection to state arrays, coupling buffers, step counts, and diagnostic inputs.

The direct Mojo bridge also preserves raw stdout cardinality: INERT must emit exactly 2N scalar lines, ordered as N torus phases followed by N frequency deviations. Blank, truncated, overlong, non-numeric, non-finite, or out-of-domain phase lines are rejected before public arrays are returned.

4. Features

4.1 Second-Order State

Unlike UPDEEngine (first-order, phase only), this engine tracks both phase \(\theta_i\) and frequency deviation \(\dot{\omega}_i\). The state space is \(2N\)-dimensional.

4.2 RK4 Integration

Classical fourth-order Runge-Kutta applied to the coupled system. Both \(\theta\) and \(\dot{\omega}\) are advanced simultaneously within each stage.

4.3 Per-Oscillator Parameters

Each oscillator has individual inertia \(M_i\) and damping \(D_i\), enabling heterogeneous networks (e.g., mix of large and small generators).

4.4 Trajectory Recording

run() returns full trajectories (theta_traj, omega_traj) in addition to final state — essential for transient stability analysis.

4.5 Diagnostic Methods

  • frequency_deviation(omega_dot) — max frequency error in Hz
  • coherence(theta) — Kuramoto order parameter R

4.6 Rust Acceleration

PyInertialStepper (spo-kernel) provides Rust-accelerated step(). run() currently loops over step() in Python (Rust batch not yet implemented). Speedup ~60x per step.

5. Usage Examples

5.1 Basic Power Grid Simulation

import numpy as np
from scpn_phase_orchestrator.upde.inertial import InertialKuramotoEngine

N = 10  # 10 generators
rng = np.random.default_rng(42)

theta = rng.uniform(0, 2 * np.pi, N)
omega_dot = np.zeros(N)  # start at nominal frequency
power = rng.normal(0, 0.1, N)  # small power imbalances
power -= power.mean()  # ensure net zero

# All-to-all coupling (K = 3.0)
knm = np.full((N, N), 3.0 / N)
np.fill_diagonal(knm, 0.0)

inertia = np.ones(N) * 5.0  # H = 5 seconds
damping = np.ones(N) * 1.0  # D = 1.0

engine = InertialKuramotoEngine(N, dt=0.01)
th_final, od_final, th_traj, od_traj = engine.run(
    theta, omega_dot, power, knm, inertia, damping, n_steps=1000
)

print(f"R = {engine.coherence(th_final):.4f}")
print(f"Max freq deviation: {engine.frequency_deviation(od_final):.4f} Hz")

5.2 Generator Trip (Large Perturbation)

# Simulate loss of generator (sudden power imbalance)
power_trip = power.copy()
power_trip[0] = 0.0  # generator 0 trips offline
power_trip -= power_trip.mean()  # redistribute

th, od = th_final.copy(), od_final.copy()
th_post, od_post, _, od_traj_post = engine.run(
    th, od, power_trip, knm, inertia, damping, n_steps=5000
)

max_dev = np.max(np.abs(od_traj_post)) / (2 * np.pi)
print(f"Max transient frequency deviation: {max_dev:.4f} Hz")
# If max_dev > 2.5 Hz: system unstable → load shedding needed

5.3 Heterogeneous Inertia

# Mix of conventional (high M) and renewable (low M) generators
inertia_mixed = np.array([8, 8, 8, 8, 8, 0.5, 0.5, 0.5, 0.5, 0.5])
# Renewables have low inertia → faster frequency swings

engine_mixed = InertialKuramotoEngine(N, dt=0.005)
_, _, _, od_traj_mixed = engine_mixed.run(
    theta, np.zeros(N), power, knm, inertia_mixed, damping, n_steps=2000
)

# Renewable generators (last 5) will show larger frequency oscillations
renewable_max = np.max(np.abs(od_traj_mixed[:, 5:])) / (2 * np.pi)
conventional_max = np.max(np.abs(od_traj_mixed[:, :5])) / (2 * np.pi)
print(f"Renewable max Δf: {renewable_max:.4f} Hz")
print(f"Conventional max Δf: {conventional_max:.4f} Hz")

5.4 Damping Ratio Analysis

K_eff = 3.0  # effective coupling
M = 5.0      # inertia
D_values = [0.1, 1.0, 5.0, 20.0]

for D in D_values:
    zeta = D / (2 * np.sqrt(K_eff * M))
    regime = "overdamped" if zeta > 1 else "underdamped"
    eng = InertialKuramotoEngine(N, dt=0.01)
    _, od_f, _, _ = eng.run(
        rng.uniform(0, 2*np.pi, N), np.zeros(N), power,
        knm, np.ones(N)*M, np.ones(N)*D, n_steps=2000,
    )
    max_f = engine.frequency_deviation(od_f)
    print(f"D={D:.1f}, ζ={zeta:.2f} ({regime}): Δf_max={max_f:.4f} Hz")

5.5 SCPN Layer 6 Configuration

# SCPN Layer 6: Planetary timescale (τ = 86400s = 24h)
N_layer6 = 4  # 4 oscillators in layer 6
omega_planet = np.ones(N_layer6) * 2 * np.pi / 86400  # 1 cycle per day
M_planet = np.ones(N_layer6) * 1.0  # normalised inertia
D_planet = np.ones(N_layer6) * 0.5  # moderate damping
P_planet = np.zeros(N_layer6)       # no external power

knm_planet = np.full((N_layer6, N_layer6), 0.5 / N_layer6)
np.fill_diagonal(knm_planet, 0.0)

engine_planet = InertialKuramotoEngine(N_layer6, dt=100.0)  # 100s timestep
th_p, od_p, _, _ = engine_planet.run(
    rng.uniform(0, 2*np.pi, N_layer6), np.zeros(N_layer6),
    P_planet, knm_planet, M_planet, D_planet, n_steps=864,
)
print(f"Planetary coherence: {engine_planet.coherence(th_p):.4f}")

5.6 Transient Stability Assessment

# Run post-fault simulation and check if system re-synchronises
# Critical clearing time: how long a fault can persist before instability

fault_durations = [0.1, 0.2, 0.5, 1.0, 2.0]  # seconds
for t_fault in fault_durations:
    n_fault = int(t_fault / 0.01)
    # During fault: one generator decoupled (row/col of knm zeroed)
    knm_fault = knm.copy()
    knm_fault[0, :] = 0.0
    knm_fault[:, 0] = 0.0

    th, od = theta.copy(), np.zeros(N)
    # Fault period
    th, od, _, _ = engine.run(th, od, power, knm_fault, inertia, damping, n_fault)
    # Post-fault recovery
    th, od, _, od_post = engine.run(th, od, power, knm, inertia, damping, 5000)

    max_dev = np.max(np.abs(od_post[-100:])) / (2 * np.pi)
    stable = max_dev < 0.5
    print(f"Fault {t_fault:.1f}s: {'STABLE' if stable else 'UNSTABLE'} "
          f"(Δf_max={max_dev:.3f} Hz)")

5.7 Energy Monitoring

# Track energy-like function along trajectory
def compute_energy(th, od, knm, inertia, power):
    kinetic = 0.5 * np.sum(inertia * od**2)
    potential = 0.0
    for i in range(len(th)):
        for j in range(i+1, len(th)):
            potential -= knm[i,j] * np.cos(th[j] - th[i])
    return kinetic + potential

th_f, od_f, th_traj, od_traj = engine.run(
    theta, np.zeros(N), power, knm, inertia, damping, n_steps=1000
)
energies = [compute_energy(th_traj[t], od_traj[t], knm, inertia, power)
            for t in range(0, 1000, 10)]
print(f"Energy: {energies[0]:.2f} → {energies[-1]:.2f} (should decrease)")

5.8 Frequency Nadir Detection

# Find the worst frequency deviation after a disturbance (frequency nadir)
_, _, _, od_traj = engine.run(
    theta, np.zeros(N), power, knm, inertia, damping, n_steps=3000
)

freq_dev_hz = np.abs(od_traj) / (2 * np.pi)
nadir_time = np.argmax(np.max(freq_dev_hz, axis=1))
nadir_value = np.max(freq_dev_hz[nadir_time])
nadir_gen = np.argmax(freq_dev_hz[nadir_time])
print(f"Frequency nadir: {nadir_value:.4f} Hz at step {nadir_time}, "
      f"generator {nadir_gen}")

6. Technical Reference

6.1 Module API

inertial

Second-order (swing-equation) Kuramoto with a 5-backend fallback chain.

Model

Each oscillator has a phase θ_i and a "frequency-deviation" ω_i ≡ dθ_i/dt. The swing equation is

M_i · d²θ_i/dt² + D_i · dθ_i/dt = P_i + Σ_j K_ij · sin(θ_j − θ_i)

and is advanced with classical explicit RK4 on the (θ, ω) pair. This is the power-grid form used in Filatrella-Nielsen-Mallick 2008.

Numerics

The derivative uses the sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion so that floating-point rounding matches the Rust kernel (spo-engine/src/inertial.rs) bit-for-bit. All five backends (Rust, Mojo, Julia, Go, Python) agree within ~1e-14 on the canonical all-to-all test problem; the dispatcher selects the fastest available path.

Classes

InertialKuramotoEngine 

InertialKuramotoEngine(n: int, dt: float = 0.01)

Second-order swing-equation Kuramoto stepper with 5-backend dispatch.

The engine's geometry is (n, dt); the step itself is stateless: (θ, ω, P, K, M, D) → (θ', ω').

Initialise the stateless inertial Kuramoto stepper geometry.

Source code in src/scpn_phase_orchestrator/upde/inertial.py 
306
307
308
309
def __init__(self, n: int, dt: float = 0.01) -> None:
    """Initialise the stateless inertial Kuramoto stepper geometry."""
    self._n = _validate_positive_int(n, name="n")
    self._dt = _validate_positive_float(dt, name="dt")
Methods:
step 
step(
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
) -> tuple[FloatArray, FloatArray]

Advance one second-order inertial Kuramoto timestep.

Parameters

theta : FloatArray Oscillator phases in radians, shape (N,). omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,). power : FloatArray Per-oscillator power injection in the swing equation, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). inertia : FloatArray Per-oscillator inertia coefficients, shape (N,). damping : FloatArray Per-oscillator damping coefficients, shape (N,).

Returns

tuple[FloatArray, FloatArray] The (θ, ω̇) state after one second-order step.

Source code in src/scpn_phase_orchestrator/upde/inertial.py 
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
def step(
    self,
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
) -> tuple[FloatArray, FloatArray]:
    """Advance one second-order inertial Kuramoto timestep.

    Parameters
    ----------
    theta : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.
    power : FloatArray
        Per-oscillator power injection in the swing equation, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    inertia : FloatArray
        Per-oscillator inertia coefficients, shape ``(N,)``.
    damping : FloatArray
        Per-oscillator damping coefficients, shape ``(N,)``.

    Returns
    -------
    tuple[FloatArray, FloatArray]
        The ``(θ, ω̇)`` state after one second-order step.
    """
    theta64 = _validate_state_array(theta, name="theta", shape=(self._n,))
    omega_dot64 = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    )
    power64 = _validate_state_array(power, name="power", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    inertia64 = _validate_positive_state_array(
        inertia,
        name="inertia",
        shape=(self._n,),
    )
    damping64 = _validate_nonnegative_state_array(
        damping,
        name="damping",
        shape=(self._n,),
    )
    knm_flat = knm64.ravel()
    backend_fn = _dispatch()
    if backend_fn is not None:
        new_theta, new_omega = backend_fn(
            theta64,
            omega_dot64,
            power64,
            knm_flat,
            inertia64,
            damping64,
            self._n,
            self._dt,
        )
        return _validate_backend_output(new_theta, new_omega, n=self._n)
    new_theta, new_omega = _python_step(
        theta64,
        omega_dot64,
        power64,
        knm_flat,
        inertia64,
        damping64,
        self._n,
        self._dt,
    )
    return _validate_backend_output(new_theta, new_omega, n=self._n)
run 
run(
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
    n_steps: int,
) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]

Integrate inertial Kuramoto dynamics and return final state plus traces.

Parameters

theta : FloatArray Oscillator phases in radians, shape (N,). omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,). power : FloatArray Per-oscillator power injection in the swing equation, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). inertia : FloatArray Per-oscillator inertia coefficients, shape (N,). damping : FloatArray Per-oscillator damping coefficients, shape (N,). n_steps : int Number of integration steps to run.

Returns

tuple[FloatArray, FloatArray, FloatArray, FloatArray] The final (θ, ω̇) plus the θ and ω̇ traces.

Source code in src/scpn_phase_orchestrator/upde/inertial.py 
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
def run(
    self,
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
    n_steps: int,
) -> tuple[
    FloatArray,
    FloatArray,
    FloatArray,
    FloatArray,
]:
    """Integrate inertial Kuramoto dynamics and return final state plus traces.

    Parameters
    ----------
    theta : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.
    power : FloatArray
        Per-oscillator power injection in the swing equation, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    inertia : FloatArray
        Per-oscillator inertia coefficients, shape ``(N,)``.
    damping : FloatArray
        Per-oscillator damping coefficients, shape ``(N,)``.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    tuple[FloatArray, FloatArray, FloatArray, FloatArray]
        The final ``(θ, ω̇)`` plus the ``θ`` and ``ω̇`` traces.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    theta_traj = np.empty((n_steps, self._n))
    omega_traj = np.empty((n_steps, self._n))
    th = _validate_state_array(theta, name="theta", shape=(self._n,)).copy()
    od = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    ).copy()
    for i in range(n_steps):
        th, od = self.step(th, od, power, knm, inertia, damping)
        theta_traj[i] = th
        omega_traj[i] = od
    return th, od, theta_traj, omega_traj
frequency_deviation 
frequency_deviation(omega_dot: FloatArray) -> float

Return maximum absolute frequency deviation in cycles per unit time.

Parameters

omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,).

Returns

float The maximum absolute frequency deviation in cycles per unit time.

Source code in src/scpn_phase_orchestrator/upde/inertial.py 
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
def frequency_deviation(self, omega_dot: FloatArray) -> float:
    """Return maximum absolute frequency deviation in cycles per unit time.

    Parameters
    ----------
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.

    Returns
    -------
    float
        The maximum absolute frequency deviation in cycles per unit time.
    """
    omega_dot64 = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    )
    return float(np.max(np.abs(omega_dot64)) / TWO_PI)
coherence 
coherence(theta: FloatArray) -> float

Return the Kuramoto order parameter for the supplied phases.

Parameters

theta : FloatArray Oscillator phases in radians, shape (N,).

Returns

float The Kuramoto order parameter R.

Source code in src/scpn_phase_orchestrator/upde/inertial.py 
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
def coherence(self, theta: FloatArray) -> float:
    """Return the Kuramoto order parameter for the supplied phases.

    Parameters
    ----------
    theta : FloatArray
        Oscillator phases in radians, shape ``(N,)``.

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    theta64 = _validate_state_array(theta, name="theta", shape=(self._n,))
    return float(np.abs(np.mean(np.exp(1j * theta64))))

Functions:

6.2 Constructor

Parameter Type Default Description
n int — Number of oscillators
dt float 0.01 Integration timestep

6.3 step() Parameters

Parameter Shape Description
theta (N,) Rotor angles
omega_dot (N,) Angular velocity deviations
power (N,) Mechanical power \(P_i\)
knm (N,N) Coupling matrix
inertia (N,) Inertia constants \(M_i\)
damping (N,) Damping coefficients \(D_i\)

Returns: (theta_new, omega_dot_new) — both (N,) arrays.

6.4 run() Returns

(theta_final, omega_final, theta_trajectory, omega_trajectory)

  • theta_trajectory: shape (n_steps, N)
  • omega_trajectory: shape (n_steps, N)

6.5 Diagnostic Methods

Method Returns Description
frequency_deviation(omega_dot) float $\max_i
coherence(theta) float $R =

7. Performance Benchmarks

7.1 Rust Speedup

N Python (µs/step) Rust (µs/step) Speedup
8 85 1.4 60x
16 140 2.8 50x
64 1400 35 40x
256 22000 560 39x

RK4 requires 4 coupling evaluations per step (each \(O(N^2)\)), making it 4x more expensive than Euler per step but with \(O(\Delta t^4)\) accuracy.

7.2 Memory

State: \(2N\) (theta + omega_dot) + trajectory: \(2 \times n_{\text{steps}} \times N\).

N n_steps Trajectory memory
16 1000 250 KB
64 5000 5 MB
256 10000 40 MB

7.3 Complexity

Operation Time Space
step() \(4 \times O(N^2)\) \(O(N^2)\) scratch
run(n) \(4n \times O(N^2)\) \(O(n \cdot N)\) trajectory

7.4 Timestep Selection

The swing equation is moderately stiff for large \(M/D\) ratios. The CFL-like stability condition for Euler would be \(\Delta t < 2D / K\). RK4 is more permissive — stable for \(\Delta t \leq 0.01 / \sqrt{K/M}\) typically.

M/D ratio Recommended dt Steps for 10s simulation
0.1 (overdamped) 0.05 200
1.0 (balanced) 0.01 1,000
10 (underdamped) 0.002 5,000
100 (very stiff) 0.0005 20,000

For very stiff systems (\(M/D > 100\)), consider implicit methods (not yet implemented in SPO).

7.5 Comparison with UPDEEngine

For the same N=64, K=3.0/N network:

Engine Method µs/step (Rust) State dim Use case
UPDEEngine Euler 57 N Fast, no inertia
UPDEEngine RK4 257 N Accurate, no inertia
InertialKuramotoEngine RK4 35 (!) 2N Inertia + damping

Note: Inertial Rust is faster than UPDEEngine RK4 Python because Rust eliminates intermediate array allocation.

8. Citations

  1. Filatrella G., Nielsen A.H., Pedersen N.F. (2008). Analysis of a power grid using a Kuramoto-like model. European Physical Journal B 61(4):485–491. doi:10.1140/epjb/e2008-00098-8

  2. Dörfler F., Bullo F. (2012). Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators. SIAM Journal on Control and Optimization 50(3):1616–1642. doi:10.1137/110851584

  3. Rohden M., Sorge A., Timme M., Witthaut D. (2012). Self-organized synchronization in decentralized power grids. Physical Review Letters 109(6):064101. doi:10.1103/PhysRevLett.109.064101

  4. Dörfler F., Bullo F. (2014). Synchronization in complex networks of phase oscillators: A survey. Automatica 50(6):1539–1564. doi:10.1016/j.automatica.2014.04.012

  5. Park R.H. (1929). Two-reaction theory of synchronous machines: Generalized method of analysis — Part I. Transactions of the AIEE 48(3):716–727. doi:10.1109/T-AIEE.1929.5055275

  6. Bergen A.R., Hill D.J. (1981). A structure preserving model for power system stability analysis. IEEE Transactions on Power Apparatus and Systems PAS-100(1):25–35. doi:10.1109/TPAS.1981.316883

Test Coverage

  • tests/test_inertial.py — 13 tests: step shape, run trajectory shapes, zero power convergence, frequency deviation bounds, coherence range, high damping overdamped behaviour, low damping oscillatory behaviour, heterogeneous inertia, Rust parity, energy-like conservation
  • tests/test_inertial_backends.py — cross-backend parity, loader contracts, Mojo stdout cardinality, direct backend validation, and public optional backend output replay for the inertial dispatcher and Rust wrapper
  • tests/test_inertial_algorithm.py, tests/test_inertial_stability.py, tests/test_inertial_config_validation.py, and tests/test_prop_swarmalator_inertial.py — focused numerical, stability, constructor-boundary, and property regressions.

Source

  • Python: src/scpn_phase_orchestrator/upde/inertial.py
  • Rust: spo-kernel/crates/spo-engine/src/inertial.rs
  • FFI: spo-kernel/crates/spo-ffi/src/lib.rs (PyInertialStepper)