The DAE and Its Jacobian

The formulation, and the one decision everything follows from

gridoxide solves the differential and algebraic halves simultaneously: one Newton per step over the whole vector, with no inner iteration between a device solver and a network solver. There is no interface error to converge away and no per-step lag between the two halves.

The trapezoidal rule on the differential half, the constraint enforced at the new point:

\[ x_{k+1} - x_k - \tfrac{h}{2}\big(f(x_k,V_k) + f(x_{k+1},V_{k+1})\big) = 0 \] \[ \mathbf{Y}V_{k+1} - I_{inj}(x_{k+1},V_{k+1}) = 0 \]

giving the Newton matrix

      ⎡ I − (h/2)·∂f/∂x       −(h/2)·∂f/∂V ⎤
  J = ⎢                                     ⎥
      ⎣    −∂I_inj/∂x        Y − ∂I_inj/∂V  ⎦

Why rectangular coordinates, and why current balance

Write the network in rectangular coordinates and its block is the real form of the admittance matrix:

Y_real = ⎡ G  −B ⎤        I = YV  ⇒  i_re = G v_re − B v_im
         ⎣ B   G ⎦                   i_im = B v_re + G v_im

which is constant for the lifetime of a topology. Nothing about it depends on x or on V. Every moving part is then a device stamp, and every device stamp is local: ∂f/∂x is one dense block per device, ∂f/∂V and ∂I/∂x couple that block to its own bus's two columns and rows, and ∂I/∂V is 2×2 on one diagonal.

A polar power-mismatch Jacobian — what the power flow builds — has none of that structure: every entry moves every iteration. So dynamics borrows that module's shape (analyze once, refill values at fixed offsets, hand the values to the backend) and none of its arithmetic.

One symbolic factorization for a whole run

The sparsity pattern is fixed for the lifetime of a topology, and every event preserves it. That is deliberate. The pattern is built from a topology superset — with every branch in service and every bus diagonal stamped, whether it needs one yet or not — and an out-of-service branch has its positions re-stamped at zero rather than dropped.

So a fault, a clearing, a branch trip and a unit trip are all numeric refills against one symbolic factorization. Nothing re-analyzes for the whole run.

Tripping a generating unit looks like the exception, since its states would leave the system. Freezing it rather than removing it keeps the layout identical — its rows become x₁ − x₀ = 0 and its Jacobian block the identity the implicit rule contributes — and freezing is also the better model: nothing in the network can observe a disconnected machine's rotor.

The Norton stamp, which changes no answer

Each machine declares a constant admittance that is stamped into Y at build, and its injection returns

I_inj = I_machine(x, V) + y_norton·V

so the y_norton·V the network row adds is exactly cancelled. The stamp is purely a conditioning device: it keeps Y diagonally dominant at generator buses without moving a single result.

For a machine with no subtransient saliency the cancellation is exact and ∂I/∂V comes out zero to the last bit — which is the sharpest available check on the whole rotor-frame algebra, and is asserted as an identity rather than a tolerance. For a salient machine a residual remains: the machine is genuinely not an impedance in the network frame, because its response depends on the rotor's orientation. That residual is the saliency.

Two integration rules, one assembly

The trapezoidal rule is second-order and A-stable, which is what a stiff phasor DAE wants. It is not L-stable: its amplification factor tends to −1 rather than 0, so a discontinuity excites a numerical oscillation that alternates sign and decays only as slowly as the physical mode does.

The remedy is a couple of backward Euler steps immediately after each event, then back to trapezoidal. Both rules are the same residual with different coefficients —

x₁ − x₀ − h(a·f₁ + b·f₀) = 0     (a, b) = (½, ½) trapezoidal
                                         (1, 0)  backward Euler

— so nothing branches on which is active beyond choosing the pair. That matters because the rule switches within a run, and two separate assemblies would be two places to get the damping wrong.

The damping is not free, and its cost is known in closed form: backward Euler evaluates δ̇ at the end of the step, so on a linearly growing speed each damping step overshoots by exactly Ω_b·a·h²/2. Two steps leave a permanent offset of Ω_b·a·h², which is 2.5e-5 rad at a 1 ms step on a typical machine.

The backends

All five of gridoxide's sparse LU backends serve this, unchanged, because the DAE Jacobian satisfies the same contract the power flow's does: pattern fixed, values change every call. Block is the exception and is refused — it assumes a uniform 2×2 block per bus, and a DAE has variable-size device blocks alongside the network's.

One caveat worth knowing if you extend this: the backends are interchangeable at solve time but not at construction time. KluNative and Klu factor numerically inside new, so they must be handed real values from an actual fill, not a placeholder pattern.