Powerflow

The powerflow problem is about the calculation of voltage magnitudes and angles for all network nodes. The solution is obtained from a subset of voltages and power injections.

Power System Model

Power systems are modeled as a network of nodes (buses) and branches (lines and transformers). Power sources (generators) and sinks (loads) can be connected to the nodes. Each node in the network is fully described by the following four electrical quantities:

• \(\vert V_k \vert\): voltage magnitude
• \(\theta_k\): voltage phase angle
• \(P_k\): active power
• \(Q_k\): reactive power

There are three types of network nodes: VD (also called the slack or reference bus — types::BusType::Slack in gridoxide's own code), PV and PQ. Depending on the node type, two of the four electrical quantities are specified.

Newton Raphson

The goal is to bring a nonlinear mismatch function \(f\) to zero. The value of the mismatch function depends on a solution vector \(x\):

\[ f(x) = 0 \]

As \(f(x)\) is nonlinear, the equation system is solved iteratively using Newton-Raphson:

\[ x_{i+1} = x_i + \Delta x_i = x_i - \textbf{J}_f(x_i)^{-1} f(x_i) \]

where \(\Delta x\) is the correction of the solution vector and \(\textbf{J}_f\) is the Jacobian matrix.

Instead of computing \(\Delta x_i = - \textbf{J}_f(x_i)^{-1} f(x_i)\), the linear equation set

\[ - \textbf{J}_f(x_i) \Delta x_i = f(x_i) \]

is solved for \(\Delta x_i\).

Iterations are stopped when the mismatch is sufficiently small:

\[ f(x_i) < \epsilon \]

Powerflow Solution

The solution vector \(x\) represents the voltage \(V\) either in polar coordinates

\[ \left [ \begin{array}{c} \delta \ \vert V \vert \end{array} \right ] \]

or rectangular coordinates

\[ \left [ \begin{array}{c} V_{real} \ V_{imag} \end{array} \right ] \]

The mismatch function \(f\) represents the power mismatch

\[ \Delta S = \left [ \begin{array}{c} \Delta P \ \Delta Q \end{array} \right ] \]

or the current mismatch

\[ \Delta I = \left [ \begin{array}{c} \Delta I_{real} \ \Delta I_{imag} \end{array} \right ] \]

This results in four different formulations of the powerflow problem:

• power mismatch function and polar coordinates
• power mismatch function and rectangular coordinates
• current mismatch function and polar coordinates
• current mismatch function and rectangular coordinates

To solve the problem using Newton-Raphson, we need to formulate \(\textbf{J}_f\) and \(f\) for each powerflow problem formulation.

Powerflow with Power Mismatch Function and Polar Coordinates

The injected power at a node \(k\) is given by

\[ S_k = V_k I_k^* \]

The current injection into any node \(k\) is

\[ I_k = \sum_{j=1}^N Y_{kj} V_j \]

Substitution yields

\[ \begin{align*} S_k &= V_k \left ( \sum_{j=1}^N Y_{kj} V_j \right )^* \\ &= V_k \sum_{j=1}^N Y_{kj}^* V_j^* \end{align*} \]

\(G_{kj}\) and \(B_{kj}\) are defined as the real and imaginary part of the admittance matrix element \(Y_{kj}\), so that \(Y_{kj} = G_{kj} + jB_{kj}\). This results in

\[ \begin{align*} S_k &= V_k \sum_{j=1}^N Y_{kj}^* V_j^* \\ &= \vert V_k \vert \angle \theta_k \sum_{j=1}^N (G_{kj} + jB_{kj})^* ( \vert V_j \vert \angle \theta_j)^* \\ &= \vert V_k \vert \angle \theta_k \sum_{j=1}^N (G_{kj} - jB_{kj}) ( \vert V_j \vert \angle - \theta_j) \\ &= \sum_{j=1}^N \left \vert V_k \vert \vert V_j \vert \angle (\theta_k - \theta_j) \right (G_{kj} - jB_{kj}) \\ &= \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( cos(\theta_k - \theta_j) + jsin(\theta_k - \theta_j) \right ) (G_{kj} - jB_{kj}) \end{align*} \]

If we perform the algebraic multiplication of the two terms inside the parentheses, and collect real and imaginary parts, and recall that \(S_k = P_k + jQ_k\), we can split this into two equations: one for the real part, and one for the imaginary part.

\[ \theta_{kj} = \theta_k - \theta_j \\ P_k = \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( G_{kj}cos(\theta_{kj}) + B_{kj} sin(\theta_{kj}) \right ) \\ Q_k = \sum_{j=1}^N \vert V_k \vert \vert V_j \vert \left ( G_{kj}sin(\theta_{kj}) - B_{kj} cos(\theta_{kj}) \right ) \]

These are called the power flow equations.

We consider a power system network having \(N\) buses. We assume one VD bus, \(N_{PV}-1\) PV buses and \(N-N_{PV}\) PQ buses. We assume that the VD bus is numbered bus \(1\), the PV buses are numbered \(2,...,N_{PV}\), and the PQ buses are numbered \(N_{PV}+1,...,N\). We define the vector of unknowns as the composite vector of unknown angles \(\theta\) and voltage magnitudes \(\vert V \vert\):

\[ x = \left[ \begin{array}{c} \theta \\ \vert V \vert \\ \end{array} \right ] = \left[ \begin{array}{c} \theta_2 \\ \theta_{3} \\ \vdots \\ \theta_N \\ \vert V_{N_{PV+1}} \vert \\ \vert V_{N_{PV+2}} \vert \\ \vdots \\ \vert V_N \vert \end{array} \right] \]

The right-hand sides of the powerflow equations for \(P_k\) and \(Q_k\) depend on the elements of the unknown vector \(x\).

Expressing this dependency more explicitly, we rewrite these equations as

\[ \begin{align*} P_k^{spec} = P_k^{calc} (x) \Rightarrow P_k^{calc} (x) - P_k^{spec} &= 0 \quad \quad k = 2,...,N \\ Q_k^{spec} = Q_k^{calc} (x) \Rightarrow Q_k^{calc} (x) - Q_k^{spec} &= 0 \quad \quad k = N_{PV}+1,...,N \end{align*} \]

We define the mismatch \({f} (x)\) as

\[ \begin{align*} f(x) = \left [ \begin{array}{c} f_1(x) \\ \vdots \\ f_{N-1}(x) \\ ------ \\ f_N(x) \\ \vdots \\ f_{2N-N_{PV} -1}(x) \end{array} \right ] = \left [ \begin{array}{c} P_2(x) - P_2 \\ \vdots \\ P_N(x) - P_N \\ --------- \\ Q_{N_{PV}+1}(x) - Q_{N_{PV}+1} \\ \vdots \\ Q_N(x) - Q_N \end{array} \right] = \left [ \begin{array}{c} \Delta P_2 \\ \vdots \\ \Delta P_N \\ ------ \\ \Delta Q_{N_{PV}+1} \\ \vdots \\ \Delta Q_N \end{array} \right ] = 0 \end{align*} \]

That is a system of nonlinear equations. The nonlinearity stems from the fact that \(P_k\) and \(Q_k\) have terms containing products of unknowns and also terms containing trigonometric functions of unknowns.

The Jacobian matrix is obtained by taking all first-order partial derivates of the power mismatch function with respect to the voltage angles \(\theta_k\) and magnitudes \(\vert V_k \vert\):

\[ \theta_{jk} = \theta_j - \theta_k \\ \begin{align*} J_{jk}^{P \theta} &= \frac{\partial P_j (x ) } {\partial \theta_k} = \vert V_j \vert \vert V_k \vert \left ( G_{jk} sin(\theta_{jk}) - B_{jk} cos(\theta_{jk} ) \right ) \\ J_{jj}^{P \theta} &= \frac{\partial P_j(x)}{\partial \theta_j} = -Q_j (x ) - B_{jj} \vert V_j \vert ^{2} \\ J_{jk}^{Q \theta} &= \frac{\partial Q_j(x)}{\partial \theta_k} = - \vert V_j \vert \vert V_k \vert \left ( G_{jk} cos(\theta_{jk}) + B_{jk} sin(\theta_{jk}) \right ) \\ J_{jj}^{Q \theta} &= \frac{\partial Q_j(x)}{\partial \theta_j} = P_j (x ) - G_{jj} \vert V_j \vert ^{2} \\ J_{jk}^{PV} &= \frac{\partial P_j (x ) } {\partial \vert V_k \vert } = \vert V_j \vert \left ( G_{jk} cos(\theta_{jk}) + B_{jk} sin(\theta_{jk}) \right ) \\ J_{jj}^{PV} &= \frac{\partial P_j(x)}{\partial \vert V_j \vert } = \frac{P_j (x )}{\vert V_j \vert} + G_{jj} \vert V_j \vert \\ J_{jk}^{QV} &= \frac{\partial Q_j (x ) } {\partial \vert V_k \vert } = \vert V_j \vert \left ( G_{jk} sin(\theta_{jk}) - B_{jk} cos(\theta_{jk}) \right ) \\ J_{jj}^{QV} &= \frac{\partial Q_j(x)}{\partial \vert V_j \vert } = \frac{Q_j (x )}{\vert V_j \vert} - B_{jj} \vert V_j \vert \\ \end{align*} \]

The linear system of equations that is solved in every Newton iteration can be written in matrix form as follows

\[ \begin{align*} -J(x) \left [ \begin{array}{c} \Delta \theta \\ \Delta \vert V \vert \end{array} \right ] &= \left [ \begin{array}{c} \Delta P \\ \Delta Q \end{array} \right ] \\ \Rightarrow J(x) \left [ \begin{array}{c} \Delta \theta \\ \Delta \vert V \vert \end{array} \right ] &= \left [ \begin{array}{c} -\Delta P \\ -\Delta Q \end{array} \right ] \end{align*} \]

\[ \begin{align*} \left [ \begin{array}{cccccc} \frac{\partial \Delta P_2 }{\partial \theta_2} & \cdots & \frac{\partial \Delta P_2 }{\partial \theta_N} & \frac{\partial \Delta P_2 }{\partial \vert V_{N_{G+1}} \vert} & \cdots & \frac{\partial \Delta P_2 }{\partial \vert V_N \vert} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \Delta P_N }{\partial \theta_2} & \cdots & \frac{\partial \Delta P_N}{\partial \theta_N} & \frac{\partial \Delta P_N}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta P_N}{\partial \vert V_N \vert} \\ \frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_2} & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \theta_N} & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_{N_{G+1}} }{\partial \vert V_N \vert} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \Delta Q_N}{\partial \theta_2} & \cdots & \frac{\partial \Delta Q_N}{\partial \theta_N} & \frac{\partial \Delta Q_N}{\partial \vert V_{N_{G+1}} \vert } & \cdots & \frac{\partial \Delta Q_N}{\partial \vert V_N \vert} \end{array} \right ] \left [ \begin{array}{c} \Delta \theta_2 \\ \vdots \\ \Delta \theta_N \\ \Delta \vert V_{N_{G+1}} \vert \\ \vdots \\ \Delta \vert V_N \vert \end{array} \right ] = \left [ \begin{array}{c} -\Delta P_2 \\ \vdots \\ -\Delta P_N \\ -\Delta Q_{N_{G+1}} \\ \vdots \\ -\Delta Q_N \end{array} \right ] \end{align*} \]

Solution Steps

1. Set the iteration counter to \(i=1\). Use the initial solution \(V_{i} = 1 \angle 0^{\circ}\)
2. Compute the mismatch vector \(f({x_i})\) using the power flow equations
3. Check the stopping criterion
   • If \(\vert \Delta P_{i} \vert < \epsilon_{P}\) for all type PQ and PV buses and
   • If \(\vert \Delta Q_{i} \vert < \epsilon_{Q}\) for all type PQ
   • Then go to step 6
   • Else, go to step 4
4. Evaluate the Jacobian matrix \(\textbf{J}_f(x_i)\) and compute \(\Delta x_i\).
5. Compute the new solution vector \(x_{i+1}\) and return to step 3.
6. Stop.

Reactive Power Limits (PV → PQ Switching)

Motivation

The plain Newton-Raphson formulation on the Powerflow page treats every PV bus as if its generator could supply or absorb any amount of reactive power needed to hold \(\vert V_k \vert\) at its setpoint. Real generators can't: each has a nameplate reactive capability \(Q_k^{min} \le Q_k \le Q_k^{max}\), usually tightest for round-rotor machines operating near their active-power limit. If the voltage setpoint on a PV bus would require \(Q_k\) outside that range, the unconstrained solution is not physically achievable — the generator saturates at its limit and the bus can no longer hold its voltage.

The standard fix, used by MATPOWER's runpf (enforce_q_lims) and implemented the same way in gridoxide, is PV → PQ switching: solve normally, check every PV bus's computed \(Q_k\) against its limits, and for any bus that violates one, convert it to a PQ bus with \(Q_k\) pinned at the violated limit — trading the voltage-magnitude equation for a reactive-power equation at that bus — then re-solve. Repeat until no bus violates its limits.

What changes in the equation system

Recall from the Powerflow page that a PV bus contributes only a \(P\)-mismatch row to \(f(x)\) — its voltage magnitude \(\vert V_k \vert\) is known (the setpoint), so it isn't part of the unknown vector \(x\) and there's no \(Q\)-mismatch row for it. A PQ bus contributes both a \(P\)- and a \(Q\)-mismatch row, and its \(\vert V_k \vert\) is part of \(x\).

Switching bus \(k\) from PV to PQ therefore:

• adds \(\vert V_k \vert\) to the unknown vector \(x\) (it's no longer fixed at the setpoint — the voltage is now free to float),
• adds a \(Q\)-mismatch row \(Q_k^{calc}(x) - Q_k^{spec}\) to \(f(x)\), with \(Q_k^{spec}\) pinned at whichever limit was violated (\(Q_k^{min}\) or \(Q_k^{max}\)) instead of the bus's original \(q\_spec\),
• leaves that bus's \(P\)-mismatch row untouched — the same equation, unaffected by bus type.

This changes \(n_{unknowns}\) itself, so the Jacobian's dimensions grow by one row and one column per switched bus. Any cached symbolic factorization (fill-reducing ordering) computed for the old sparsity pattern is invalid once a switch happens and must be redone from scratch.

The outer-loop algorithm

solver::newton_raphson_enforcing_q_limits(buses, ybus, tol, max_iter, backend, max_outer_iter) implements this as an outer loop around a solver::PersistentSolver:

1. Solve to convergence with every PV bus free (PersistentSolver::solve). If this doesn't converge, stop and return that status — there's no point checking limits on a non-converged solution.
2. Compute every bus's actual \(Q_k\) (network::power_injections).
3. For each PV bus whose \(Q_k\) violates q_min/q_max: switch bus_type to PQ and pin q_spec to the violated limit.
4. If no bus was switched this pass, the solution is self-consistent — stop and return Converged.
5. Otherwise, reset the cached factorization (PersistentSolver::reset) and go back to step 1.
6. If max_outer_iter outer passes pass without stabilizing, stop and return MaxIterationsReached.

Two simplifications versus a fully general implementation (deliberate, and both match MATPOWER's own default runpf behavior rather than being a gridoxide-specific shortcut):

• Switching is one-directional. A bus switched to PQ never switches back to PV within the same call, even if a later pass's solution would put it back within limits. A bidirectional scheme needs real anti-oscillation logic (a bus can otherwise flip back and forth indefinitely across outer passes); this avoids that whole problem by not attempting it.
• Bus voltages carry over between outer passes rather than resetting to a flat start. Since only one bus's type changes per pass, the previous pass's converged state is normally very close to the next equilibrium, so re-solving after a switch typically takes only a few extra Newton iterations, not a full fresh solve.

Alternative switching strategies

One-directional switching is a deliberate simplification, not the only viable design. Two genuinely different strategies exist for the same underlying problem:

• Bidirectional switching with a hard switch-count cap. A bus moved to PQ is allowed to move back to PV in a later pass if its computed \(Q\) returns within limits — recovering a bus that was only transiently out of range (e.g. during an early, poorly-conditioned iteration far from the solution) instead of committing it to PQ for the rest of the solve. The risk this reintroduces is a bus oscillating between PV and PQ indefinitely across passes; the standard mitigation (used this way by powsybl-open-loadflow's ReactiveLimitsOuterLoop) is a hard cap on how many times any one bus may switch, after which it's forced to stay wherever it last was — real bookkeeping gridoxide's one-directional rule avoids needing at all.
• Checking mid-iteration instead of only after full convergence. Rather than a separate outer loop that fully re-solves after every switch, the check-and-switch step can instead be folded directly into the Newton iteration itself once the mismatch drops below some threshold, switching bus types and continuing the same iteration sequence rather than starting a fresh solve. One published version of this technique ("On PV-PQ Bus Type Switching Logic in Power Flow Computation", Jinquan Zhao — cited directly in VeraGrid's own source as the basis for its implementation) also bases the switch-back decision on comparing the bus's voltage against its setpoint in addition to its \(Q\) against its limits, not \(Q\) alone.

Both alternatives add real complexity (anti-oscillation bookkeeping, or reworking the iteration loop itself) to recover a case gridoxide's simpler design just accepts as a one-way commitment — a deliberate scope trade-off matching MATPOWER's own default behavior, not an oversight.

Scope: net injection, not per-device limits

Bus::q_min/q_max bound the bus's net reactive injection — the same aggregate quantity q_spec already represents when multiple loads/generators share a node (see pgm::PgmVoltageRegulator's own q_min/q_max parsing, which sums each active generator's limits at a bus exactly the way p_spec is already summed). Pinning q_spec to a violated limit only achieves that limit exactly if the bus carries no voltage-dependent ZIP-model load terms (Bus::zip_terms) — true for the common case of a bus that's purely a PV/generator connection, not true in general if a voltage-dependent load is co-located there.

Validated against real data

tests/q_limits_test.rs checks the mechanics directly: a 3-bus fixture with one PV bus whose q_min is set tight enough to force a violation, confirmed to switch to PQ, pin q_spec at exactly q_min, and converge to the same voltages across every Jacobian backend (Scalar, Block, Klu) — as well as a control case with a loose q_min that shouldn't trigger any switch at all, converging to the same answer as plain unconstrained newton_raphson.

Beyond the unit fixture, this has been exercised against all 12 real MATPOWER benchmark cases (scripts/bench/matpower_to_pgm.py populates voltage_regulator.q_min/q_max from each case's own gen matrix): 11 of the 12 have at least one PV bus whose unconstrained \(Q\) genuinely exceeds its nameplate limit — from 4 violations on the smallest case up to 166 simultaneous violations on case3120sp.