Power Flow: The Core Computation

Per-unit normalisation

Before any equations, one convention. Power system quantities are expressed per unit: each variable is divided by a base value for its type.

S_base   chosen once for the system (commonly 100 MVA)
V_base   chosen per voltage level (e.g. 400 kV, 220 kV, 20 kV)
Z_base   = V_base² / S_base
I_base   = S_base / (√3 · V_base)

The payoff is that transformer ratios disappear from the equations when bases are chosen consistently across levels, and that all quantities land in a similar numerical range, which conditions the matrices well. Voltages sit near 1.0 p.u., impedances are small fractions, and an out-of-range value is immediately visible as an error. Every serious tool works internally in per unit.¹

The network equations

Each bus i has a complex voltage V_i = |V_i| ∠ θ_i. The network’s electrical behaviour is captured by the bus admittance matrix Y, built by inspection:

Y_ii = sum of admittances of all branches connected to bus i (plus shunts)
Y_ij = − (admittance of the branch between i and j),  for i ≠ j

Y is complex, symmetric for networks without phase shifters, and extremely sparse — a bus connects to a handful of neighbours regardless of system size. Exploiting that sparsity is what makes solving 10,000-bus systems routine.²

Injected complex power at bus i is S_i = V_i (Σ_k Y_ik V_k)*, which separates into the power flow equations:

P_i = |V_i| Σ_k |V_k| ( G_ik cos(θ_i − θ_k) + B_ik sin(θ_i − θ_k) )
Q_i = |V_i| Σ_k |V_k| ( G_ik sin(θ_i − θ_k) − B_ik cos(θ_i − θ_k) )

These are nonlinear in the voltage magnitudes and angles, which is the source of every difficulty that follows.

Bus types and the unknowns

There are four quantities per bus (P, Q, |V|, θ) and two equations, so two must be specified:

The slack bus exists for two reasons: angles are only meaningful relative to a reference, and total losses are not known until the flow is solved, so some bus must take up the residual. Where the slack sits affects the reported generation at that bus but not the physics — though in a large model it is normal to distribute the slack across many generators to avoid a single fictitious swing.

Solving it

Newton–Raphson is the standard method.³ Linearise the mismatch equations around the current iterate:

[ ΔP ]   [ ∂P/∂θ   ∂P/∂|V| ] [ Δθ  ]
[    ] = [                 ] [     ]
[ ΔQ ]   [ ∂Q/∂θ   ∂Q/∂|V| ] [ Δ|V|]

Solve the sparse linear system for the update, apply, repeat. Convergence is quadratic near the solution — typically 3–5 iterations from a flat start (all angles zero, all magnitudes 1.0 p.u.) for a well-conditioned case.

Variants worth knowing:

• Fast Decoupled Load Flow (FDLF) exploits the weak coupling between P–θ and Q–|V| in high-X/R networks to use two constant, pre-factorised matrices. Slower per-iteration convergence, far cheaper per iteration; still widely used for contingency screening.⁴
• Gauss–Seidel is simple and robust from poor starting points but converges slowly; occasionally used to produce a starting point for Newton.
• Continuation power flow traces the solution as load is scaled up, following the curve around the nose point where the Jacobian becomes singular. This is how voltage stability margins are computed — the point where Newton stops converging is a physical result, not just a numerical failure.⁵

When power flow fails to converge, the usual causes are, in order of frequency: an islanded subnetwork with no slack; a data error producing an implausible impedance; a case with genuinely insufficient reactive support (no solution exists); or an ill-conditioned start. Distinguishing “no solution exists” from “the solver failed” is a real skill and a real risk — a non-converged case reported as a violation is a false alarm; a genuinely infeasible case dismissed as a numerical artefact is worse.

The DC approximation

For high-voltage meshed networks, three approximations hold well enough to be extraordinarily useful:

1. R ≪ X, so resistance (and therefore losses) is neglected.
2. Voltage magnitudes are near 1.0 p.u. everywhere.
3. Angle differences across branches are small, so sin(θ_i − θ_j) ≈ θ_i − θ_j and cos(·) ≈ 1.

The result is linear:

P_ij = (θ_i − θ_j) / x_ij            branch flow
P    = B' θ                          nodal balance

where B' is the real, sparse susceptance matrix. This is the DC power flow.⁶ Its consequences are what make it dominant in market and planning applications:

• It is a single linear solve, not an iteration, so it always “converges”.
• It can be embedded in a linear program, which is why every LMP market, every security-constrained unit commitment, and virtually every capacity expansion model uses it.
• Superposition applies, which gives us the sensitivity factors below.

Its costs are equally important. It has no reactive power, no voltages, and no losses, so it cannot detect a voltage violation and will systematically under-report flows on loaded corridors where losses matter. Its accuracy degrades on low-voltage networks (where R ≈ X or R > X), on heavily loaded systems with large angle differences, and anywhere the answer depends on reactive support. Stott, Jardim and Alsaç’s assessment is worth reading in full before relying on it: the linearisation is defensible, but the choice of reactance values and the handling of losses matter more than practitioners assume.⁷

Sensitivity factors

Linearity yields two matrices that carry an enormous amount of practical work.

PTDF — Power Transfer Distribution Factor

PTDF(l, i) is the change in flow on branch l per unit of power injected at bus i and withdrawn at the reference. From the DC equations:

PTDF = B_branch · B'⁻¹     (with the reference bus row/column removed)

PTDFs answer “if I move 100 MW from here to there, how much of it shows up on that line?” — and the answer is usually “a surprising amount of it on lines nowhere near either end”. A zonal PTDF, obtained by aggregating nodal PTDFs with a generation shift key describing how a zone’s net position is spread across its nodes, is the object at the heart of European flow-based market coupling.⁸

LODF — Line Outage Distribution Factor

LODF(l, k) is the fraction of the pre-outage flow on branch k that redistributes onto branch l when k trips. Given the PTDF matrix, LODFs follow in closed form, and post-contingency flows come from a matrix multiplication:

P_l^post = P_l + LODF(l,k) · P_k

This is why full N-1 screening over thousands of contingencies is computationally cheap in DC: no power flow is re-solved. Only the contingencies that screen as potentially violating are then re-run in AC. That two-stage structure — cheap linear screening, expensive nonlinear verification — is the standard architecture of a security analysis pipeline.⁹

Optimal power flow

Power flow evaluates a given dispatch. Optimal power flow (OPF) chooses it: minimise cost (or losses, or deviation) over controllable injections and setpoints, subject to the power flow equations and to operating limits.

min   Σ c_g(P_g)
s.t.  power flow equations
      P_g^min ≤ P_g ≤ P_g^max
      |P_l|   ≤ P_l^max
      V^min   ≤ |V| ≤ V^max

Two facts shape everything built on OPF:

The AC form is nonconvex and NP-hard in general.¹⁰ Practical solvers find local optima and are sensitive to initialisation. Convex relaxations — semidefinite, second-order cone, and the QC relaxation — are exact on some network classes (notably radial distribution feeders under conditions) and provide bounds elsewhere; this is one of the more active research areas in the field.¹¹

The DC form is a linear program, therefore solvable at scale with guaranteed global optimality and, crucially, with duals. The dual variable on each nodal balance constraint is the locational marginal price. This is not a modelling convention layered on top of the physics; it falls directly out of the optimisation, and it is the reason LMP markets and DC-OPF are inseparable.¹²

Adding contingency constraints gives SCOPF (security-constrained OPF): the dispatch must be feasible in the base case and survive every credible contingency, either without action (preventive) or with defined post-contingency actions (corrective). SCOPF is the formal statement of what a system operator is actually solving, and its size — constraints scaling as branches × contingencies — is why constraint screening and iterative constraint generation are standard practice.

Choosing a formulation

Summary

Power flow is a nonlinear system solved by Newton–Raphson in per-unit quantities on a sparse admittance matrix. The DC linearisation trades reactive power, voltages, and losses for a linear model that supports superposition, sensitivity factors, and embedding in optimisation — the trade that market and planning models almost always take. OPF turns evaluation into decision, and the dual variables of its linear form are prices.

Next: how operators use these tools to keep the system secure.

References