Power flow

Three engines implement IPowerFlowSolver: the dense Newton-Raphson (src/GridSim.Core/PowerFlow/NewtonRaphsonPowerFlow.cs), the sparse Newton-Raphson (src/GridSim.Core/PowerFlow/SparseNewtonRaphsonPowerFlow.cs) and the backward/forward sweep for distribution feeders (src/GridSim.Core/PowerFlow/BackwardForwardSweep.cs). Everything that does not depend on the matrix representation - bus-spec assembly, the Q-limit outer loop, result assembly, line flows - is shared, as pure functions, in src/GridSim.Core/PowerFlow/NewtonShell.cs, so the dense and sparse engines produce bit-for-bit identical results on the same case.

The Newton-Raphson formulation as implemented

Bus spec and start point

NewtonShell.BuildBusSpec aggregates generators and loads onto buses:

Pspec[i] = (sum Pg - Pd) / BaseMva, Qspec[i] = (sum Qg - Qd) / BaseMva.
Per-bus aggregate Qmin/Qmax (summed over the bus's generators) and the generator voltage setpoint Vg are kept for the Q-limit loop.
Start point: a voltage-controlled bus (Slack or PV) with a generator starts at Vg; otherwise Vm = 1.0 under FlatStart, else the bus's model Vm. The slack angle is always the model's VaDeg (in radians internally); other angles start at 0 under FlatStart, else at the model values.

Mismatch and unknown ordering

Each iteration evaluates the polar power injections at the current state (NewtonRaphsonPowerFlow.ComputeInjections - internal so the state estimator's measurement model evaluates h(x) through the exact same code path):

P[i] = Vm[i] * sum_k Vm[k]*(G[i,k]*cos thetaik + B[i,k]*sin thetaik)
Q[i] = Vm[i] * sum_k Vm[k]*(G[i,k]*sin thetaik - B[i,k]*cos thetaik)     thetaik = Va[i] - Va[k]

The mismatch vector stacks real power over the pvpq buses (every non-slack bus, in bus order) then reactive power over the pq buses:

f = [ Pspec - P  over pvpq   (npvpq rows)
      Qspec - Q  over pq     (npq rows)   ]        m = npvpq + npq

Convergence is max|f| < Tolerance, checked before the linear solve, so a converged start point costs zero factorisations. Both Newton engines also call PowerFlowOptions.Cancellation.ThrowIfCancellationRequested() at the top of every inner iteration (~ms granularity), so a long deferrable solve can be aborted the instant a cycle deadline passes; the default None token never throws, so an ordinary solve is unchanged. The unknowns follow the same ordering: deltatheta for pvpq buses, then deltaVm for pq buses. The dense solver keeps explicit pvpq/pq index lists; the sparse solver keeps two position arrays pPos[i]/qPos[i] (-1 when bus i has no P/Q equation) that map a bus to its row/column in the linear system - the same ordering, indexed differently.

Jacobian blocks

Both engines build the standard polar blocks H, N, M, L (dP/dtheta, dP/dVm, dQ/dtheta, dQ/dVm), with the diagonal entries expressed through the already-computed injections:

Entry	i = j (diagonal)	i != j (off-diagonal)
H = dP/dtheta	`-Q[i] - B[i,i]*Vm[i]^2`	`Vm[i]Vm[j](Gsin theta - Bcos theta)`
N = dP/dVm	`P[i]/Vm[i] + G[i,i]*Vm[i]`	`Vm[i](Gcos theta + B*sin theta)`
M = dQ/dtheta	`P[i] - G[i,i]*Vm[i]^2`	`-Vm[i]Vm[j](Gcos theta + Bsin theta)`
L = dQ/dVm	`Q[i]/Vm[i] - B[i,i]*Vm[i]`	`Vm[i](Gsin theta - B*cos theta)`

The dense engine writes them into a reused double[m,m] (NewtonRaphsonPowerFlow.BuildJacobianInto); the sparse engine walks each bus's Y-bus row and emits only structurally non-zero entries into a CsrMatrix.Builder (SparseNewtonRaphsonPowerFlow.BuildJacobian) - the Jacobian inherits the admittance graph, four blocks interleaved.

Linear solve and update

Dense: ILinearSolver.SolveInPlace(J, f) - partial-pivoting LU; f becomes the correction. Updates are undamped: Va += delta on pvpq, Vm += delta on pq. A near-singular pivot throws.
Sparse: ISparseLinearSolver.Solve(J, f) (direct Gilbert-Peierls LU by default). If it reports singular (null), the engine falls back to ILU(0)-preconditioned BiCGSTAB with tol = min(1e-10, Tolerance*1e-3) and maxIter = 4m + 200; if that also fails, the iteration stops and the solve reports non-convergence. A NaN/Inf correction is never applied. The sparse update is adaptively damped: if the largest correction component exceeds 0.4 the whole step is scaled to 0.4 (direction preserved), and Vm is floored at 0.2 pu - guards against a large flat-start step on a stiff network overshooting into a negative-voltage region. Near the solution the cap never binds, so quadratic convergence is intact.

The reactive-limit outer loop

With EnforceQLimits (default on), after each inner Newton run NewtonShell.ApplyQLimits checks every bus that started PV:

Forward: if the bus's generators would need more than the aggregate Qmax (or less than Qmin), within QLimitTolMva = 1e-3 MVA, the output is pinned at the limit and the bus is released to PQ (Qspec fixed at the limit) - the numerical analogue of a machine hitting its excitation limit. The pinned limit is recorded for revert bookkeeping.
Revert (opt-in, RevertQLimits, default off): a pinned-PQ bus returns to PV when its solved voltage shows the limit is no longer binding - pinned at Qmax with Vm > Vg + 1e-4, or at Qmin with Vm < Vg - 1e-4 (RevertVoltTolPu). At most one revert per bus per solve (anti-oscillation); a re-violating bus is re-pinned and stays pinned. The default stays off so historical one-way switching - and every pinned published-solution regression - is unchanged.

The loop re-solves until no bus switches, capped at MaxQLimitPasses (12). The slack bus is never limited; NewtonShell.DiagnoseSlackQ reports (as SlackQNote on either Newton engine, diagnostic only) when the slack's solved reactive injection exceeds the aggregate limits of the machines on it - reassigning the slack or a distributed slack is deliberately out of scope.

Options

src/GridSim.Core/PowerFlow/PowerFlowOptions.cs:

Option	Default	Meaning
`Tolerance`	`1e-8`	Convergence bound on `max\|f\|`, per-unit power.
`MaxIterations`	30	Inner Newton cap (per Q-limit pass).
`FlatStart`	false	Start PQ voltages from 1.0 / angles from 0.
`EnforceQLimits`	true	Run the PV->PQ outer loop.
`MaxQLimitPasses`	12	Outer-loop cap.
`RevertQLimits`	false	Opt-in PQ->PV revert (MATPOWER `enforce_q_lims` style).
`Trace`	null	Opt-in `SolveTrace` capture of everything the solve computes (Y-bus, spec, per-iteration mismatch/correction/Jacobian, Q-limit events); null = zero overhead.
`Cancellation`	`None`	`CancellationToken` checked at the top of every NR inner iteration in both `NewtonRaphsonPowerFlow` and `SparseNewtonRaphsonPowerFlow`, so the parallel cycle engine can abort a deferrable solve at its deadline. `default`/`None` is a guaranteed no-op - `ThrowIfCancellationRequested` never throws for it, so every existing solve is bit-for-bit unchanged (73 solver-regression tests green; an explicit test pins exact `Vm`/`VaDeg` equality between the token and no-token paths).

The OLTC outer loop

src/GridSim.Core/PowerFlow/SteadyStateSolver.cs:SteadyStateSolver.Solve wraps any IPowerFlowSolver in tap regulation. With no OLTCs it is the plain power flow. Otherwise, after each solve, every OLTC is examined against the solved voltage of its controlled bus:

Vm < Target - Deadband -> tap down one StepPu (bounded by MinTap)
- lowering the tap raises the controlled bus;
Vm > Target + Deadband -> tap up one step (bounded by MaxTap);
inside the deadband -> no move.

If any tap moved, the branch list is cloned with the new ratios (Dispatcher.Clone) and the flow re-solved; the loop ends when no tap moves or after maxOltcPasses (default 25). One step per pass per transformer means convergence is monotone and bounded, at the cost of one full power flow per step of total tap travel.

Two properties matter to consumers:

RegulatedModel - the model actually solved, whose Branch.TapRatios match the returned result. The loop solves on internal clones, so anything needing post-solve taps (the live dashboard, SCADA) must read this, never the model it passed in.
OltcPasses - how many re-solve passes the last call took.

Both are mutable per-solve state: the class is not thread-safe. Parallel N-1 sweeps use one instance per thread; a concurrent Solve on a shared instance is detected by an interlocked guard and throws rather than serialising, because sharing is a caller bug.

The distribution BFS sweep solver

BackwardForwardSweep handles what transmission-oriented NR does badly: radial feeders with high R/X ratios. It requires a single slack (the source) and a connected, at-most-weakly-meshed topology; each iteration is two passes over the feeder tree:

Backward (leaves -> source): node currents I = conj(S_draw/V) + y_shunt*V are accumulated up the tree to give each branch's series current. Folding across a tapped branch applies the ideal-transformer ratio (I/conj(tap) when the parent is the tapped from end, conj(tap)*I when it is the to end).
Forward (source -> leaves): voltages step across the series impedance: V_to = V_from/tap - Z*J, or for a child on the tapped from side V_from = tap*V_to - |tap|^2*Z*J. Half line-charging appears as end shunts, the from-end half referred through 1/|tap|^2 - exactly the shunt part of Yff, so the converged voltages solve the same MATPOWER pi+tap model as the Newton engines, not just at unity tap.

On top of the plain radial sweep:

Weakly-meshed feeders. A spanning tree is grown from the source; the leftover branches become breakpoints. A loop-impedance matrix is built once by superposition (one linear sweep per loop, loads off, slack grounded) and each iteration solves it (small complex Gauss-Jordan with partial pivoting) to inject compensating breakpoint currents until every loop's closure residual V_f - V_t - z*I vanishes - Shirmohammadi-style multi-port compensation. At most MaxCompensatedLoops = 6 loops are accepted; anything more meshed (IEEE 14-bus already has 7 independent loops) is declined as "not a distribution feeder".
PV buses. A non-slack regulating generator holds |V| at its setpoint by adjusting a reactive correction each sweep, using the source->bus path reactance as the dV/dQ sensitivity, clamped to the aggregate Q-limits; a saturated bus acts as PQ at the limit.

Convergence requires all three at once: the largest per-bus voltage change below Tolerance, all PV setpoint errors and all loop residuals below max(Tolerance, 1e-6) (controls settle to a looser tolerance than voltages). A disconnected network, a missing slack, or too many loops returns the non-converged report (Vm = 1, MaxMismatch = inf) rather than a doubtful answer. Result assembly reuses the pi+tap stamps, computes the real power mismatch per bus as a check figure (a regulating PV bus's Q is not checked - it solves its own reactive, like the slack), and applies the same non-finite guard as the Newton engines.

Engine selection

src/GridSim.Core/PowerFlow/SolverSelector.cs:SolverSelector is the single source of truth; the CLI and both dashboards route through it so a case always solves with the same engine.

Rule	Value	Engine
`LooksRadial`: `Branches.Count == BusCount - 1` and exactly one slack	-	`DistributionSweep` (BFS)
otherwise, `BusCount > sparseThreshold`	`DefaultSparseThreshold = 3000`	`SparseNewton`
otherwise	-	`DenseNewton`

The radial test only catches genuine tree networks - any meshed transmission case (case9 upward) fails it. The 3,000-bus threshold is the dense LU memory/fill-in ceiling; the same constant gates the estimator's dense/sparse step choice (EstimationOptions.SparseBusThreshold) and the CLI's SolverBackend.Pick. Create returns the engine plus a display label ("Newton-Raphson (dense)"/"(sparse)"/"Backward-forward sweep") - Auto is always resolved, never returned.

`PowerFlowResult` and `LineFlow`

src/GridSim.Core/PowerFlow/PowerFlowResult.cs:

Field	Meaning
`Status`	`Converged` / `MaxIterations` / `NonFinite`. `NonFinite` means a NaN/Inf appeared in the solved state - the linear solve diverged - and the state is invalid even if a stale mismatch reads below tolerance. Defaults to `Converged` for presolved states loaded from disk/bridge that only carry the boolean.
`Converged`	`true` only when `Status == Converged`. Every other status means the state is untrustworthy.
`Iterations`	Total inner Newton iterations across all Q-limit passes (or sweep iterations for BFS).
`MaxMismatch`	The final `max\|f\|` in per-unit (for a non-converged BFS, the larger of mismatch and last voltage delta).
`BusIds`	Result index -> external bus id; the only way to map arrays back to the model.
`Vm`, `VaDeg`	Solved voltages, pu and degrees.
`PInjMw`, `QInjMvar`	Net injection (generation - demand) per bus, engineering units.
`TotalLossMw`, `TotalLossMvar`	Sum of net injections = network losses.
`LineFlows`	Per-branch flows (below).

LineFlow(From, To, PFromMw, QFromMvar, PToMw, QToMvar, PLossMw, QLossMvar, Name, RatingMva, LoadingPercent) conventions (NewtonShell.ComputeLineFlows):

PFrom/QFrom and PTo/QTo are the complex power flowing into the branch at each end, computed from the pi+tap admittances (S_f = V_f*conj(Yff*V_f + Yft*V_t)*BaseMva, likewise for the to end). A through-flow is therefore positive at one end and negative at the other.
PLossMw = PFromMw + PToMw (and likewise Q) - the loss in this branch; charging makes QLoss negative on lightly loaded lines.
LoadingPercent = max(|S_f|, |S_t|) / RateAMva x 100, and 0 when the branch is unrated (RateAMva = 0).
From/To are external bus ids, not indices.