Newton-Raphson and Fast Decoupled Load Flow

Chapter 19 solved the three-bus test system in twelve sweeps and left a warning: the sweep count grows roughly in proportion to the number of buses. The reason is structural. A Gauss-Seidel update at bus \(i\) consults only the buses directly connected to it. Information about a disturbance at one end of the network therefore travels at most a few buses per sweep, and a system whose electrical diameter is thirty buses needs at least thirty sweeps before the far end has even heard about the near end.

Look at the same fact through the error. Each sweep multiplies the error by a fixed ratio \(\rho\), so the number of iterations to reach a tolerance \(\varepsilon\) is

Every extra decimal place costs the same fixed number of sweeps, and \(\rho\) itself creeps toward one as the network grows. What is wanted is a method whose ratio is not constant but improves as the solution is approached. Such a method exists, and it has been known since the seventeenth century.

The change of viewpoint is this. Gauss-Seidel treats the equations one at a time and never assembles a picture of how a change at one bus affects the mismatch at another. Newton-Raphson builds exactly that picture — the matrix of all sensitivities \(\partial(\text{mismatch})/\partial(\text{voltage})\) — and uses it to compute a correction to every unknown simultaneously. The network's connectivity is accounted for in one step rather than propagated sweep by sweep.

Nothing is free. Each Newton iteration must form and solve a linear system whose order is the number of unknowns counted in Section 18-7, and that is far more arithmetic than a Gauss-Seidel sweep. The method won the argument only when sparse matrix techniques — storing and factorising a matrix that is \(99\%\) zeros without touching the zeros — made a single iteration affordable. That happened around 1967, and the load flow has been a Newton-Raphson calculation ever since.

Take a single equation \(g(x)=0\) with a solution \(x^{*}\), and an estimate \(x^{(0)}\) which is wrong by an unknown amount \(\Delta x^{(0)}\). By definition

Expand the left side in a Taylor series about the estimate:

The series is exact but useless, because it contains every power of the unknown correction. Truncate it after the linear term. What remains is a single linear equation that can be solved for \(\Delta x^{(0)}\) at once:

Why the error squares follows directly from what was thrown away. Write \(e^{(k)} = x^{(k)}-x^{*}\) and expand \(g\) about the true root instead:

The new error is proportional to the square of the old one. The truncated second-order term is not an incidental loss; it is the entire remaining error, and squaring a small number makes it very much smaller. An error of \(10^{-2}\) becomes \(10^{-4}\), then \(10^{-8}\), then \(10^{-16}\) — at which point the arithmetic has run out of digits.

Two properties of the method must be stated alongside the good news. The derivative appears in the denominator, so a point where \(g'\) is near zero produces an enormous correction and can throw the iteration anywhere. And the quadratic behaviour is local: it is guaranteed only once the estimate is close enough. From a distant start Newton's method may converge to a different root, or oscillate, or diverge. Example 1 shows both roots of a load flow being found by the same iteration from two different starting points — which is not a defect but the two-solution structure of Section 18-2 reappearing.

Nothing in the derivation depended on there being one equation. Let a set of \(n\) nonlinear functions of \(n\) variables be specified to take given values:

Move the specified values across to define residual functions that must vanish at the solution — these are the mismatches of Chapter 18 in general dress:

Assume an estimate \(x_k^{(0)}\) and unknown corrections \(\Delta x_k^{(0)}\) such that \(g_k\big(x_1^{(0)}+\Delta x_1^{(0)},\dots\big) = 0\) for every \(k\). Expand each \(g_k\) in a Taylor series about the estimate and keep only the first-order terms — now there are \(n\) partial derivatives per function instead of one:

Writing all \(n\) of these equations together, the array of partial derivatives is the Jacobian matrix \(\mathbf{J}\), and the whole set becomes a linear system in the correction vector:

Update and repeat:

Because the series was truncated, one step does not land on the answer; the Jacobian is rebuilt at the new point and the step is taken again, until every \(|\Delta g_k|\) falls below the tolerance. The inverse is written for clarity only. No program ever forms \(\mathbf{J}^{-1}\): the linear system is solved by \(LU\) factorisation, which for a sparse Jacobian costs a small multiple of \(n\) rather than \(n^3\).

Now identify the general machinery with the load flow. The variables \(x\) are the unknown voltage angles and magnitudes; the functions \(f\) are the power flow equations of Section 18-6; the specified values \(\eta\) are the scheduled injections; and the residuals \(g\) are exactly the mismatches \(\Delta P_i\) and \(\Delta Q_i\).

Chapter 18 counted the unknowns. For a system with \(n\) buses, one slack and \(n_g\) PV buses, there are \(n-1\) unknown angles and \(n_p = n-1-n_g\) unknown magnitudes. Ordering the mismatches with all the \(\Delta P\) above all the \(\Delta Q\), and the corrections with all the \(\Delta\delta\) above all the \(\Delta|V|\), the Jacobian falls naturally into four blocks:

The block dimensions follow from the counts:

and each block collects one family of derivatives:

The total order, \((n-1)+n_p = 2n-2-n_g\), is precisely the number of unknowns counted in Section 18-7 — as it must be, or the problem would not be well posed. For a \(300\)-bus system with \(40\) generators the Jacobian is \(559\times559\); for a modern interconnection with \(50{,}000\) buses it has of the order of \(10^{5}\) rows. Neither is solvable by dense algebra, and neither needs to be. The Jacobian inherits its sparsity pattern directly from \(\mathbf{Y}_{bus}\): the entry \(\partial P_i/\partial\delta_k\) is non-zero only if buses \(i\) and \(k\) are directly connected, which for a real network means three or four entries per row.

The elements are obtained by differentiating the polar power flow equations of Chapter 18. Separating the \(k=i\) term, whose angle argument \(\theta_{ii}+\delta_i-\delta_i = \theta_{ii}\) contains no angle difference at all, gives the form that is easiest to differentiate:

using \(Y_{ii}=G_{ii}+jB_{ii}\), so that \(|Y_{ii}|\cos\theta_{ii}=G_{ii}\) and \(|Y_{ii}|\sin\theta_{ii}=B_{ii}\).

The block \(\mathbf{J}_{11}\), written \(L_{ik}=\partial P_i/\partial\delta_k\). For \(k\neq i\), only the single term containing \(\delta_k\) survives the differentiation, and \(\delta_k\) enters the cosine with a plus sign:

For \(k=i\), the angle \(\delta_i\) appears in every term of the sum, with a minus sign, so every term contributes and the signs reverse:

That is a small piece of good fortune with a large payoff. The diagonal element requires no fresh trigonometry at all: the injected power \(Q_i\) has already been computed to form the mismatch, and \(B_{ii}\) is stored data. The same collapse happens in all four blocks, and it is the reason a Jacobian can be assembled almost as cheaply as the mismatches themselves.

Working through the other three blocks in the same way — differentiating the sums, then recognising the result as a combination of \(P_i\), \(Q_i\) and a diagonal admittance — gives the complete set.

Three points of housekeeping complete the picture. A row of \(\Delta Q\) and a column of \(\Delta|V|\) exist only for PQ buses, so the rows and columns belonging to PV buses are simply deleted from \(\mathbf{J}_{12}\), \(\mathbf{J}_{21}\) and \(\mathbf{J}_{22}\). The slack bus contributes neither a row nor a column anywhere. And the whole Jacobian must be rebuilt at every iteration, because every element depends on the present voltages — this is the expensive part of the method and the target of every approximation in Sections 20-7 to 20-9.

The procedure is the general iteration of Section 20-3 with the load flow quantities substituted in.

Step 8 deserves attention because it is where the method's dimension changes during the run. Reclassifying a PV bus adds \(\Delta Q_i\) to the mismatch vector and \(\Delta|V_i|/|V_i|\) to the correction vector, so \(\mathbf{J}\) grows by one in each direction and its factorisation must be redone. A case with many limit violations therefore converges more slowly than one with none — not because Newton's method is degraded, but because the problem being solved keeps changing. Programs usually delay limit enforcement until the second or third iteration, by which time the reactive estimates are meaningful and spurious switching is avoided.

The correction vector also gives a second, sharper convergence test. Once the corrections themselves become small, they are approximately equal to the remaining error, so \(\max|\Delta\delta|\) and \(\max|\Delta|V||\) below about \(10^{-5}\) confirm the mismatch test rather than substituting for it.

Chapter 18 noticed a structural feature of the power flow equations and promised to exploit it here. In a transmission network the series reactance dominates the resistance, \(X \gg R\), so in \(Y_{ik}=G_{ik}+jB_{ik}\) the susceptance dominates: \(|B_{ik}| \gg |G_{ik}|\), and the admittance angle \(\theta_{ik}\) is close to \(\pm90^\circ\). At the same time, angle differences across a healthy line are small, a few degrees at most, so \(\delta_k-\delta_i \approx 0\).

Apply both facts to the off-diagonal Jacobian entries. With \(\alpha_{ik}=\theta_{ik}+\delta_k-\delta_i \approx \theta_{ik} \approx 90^\circ\),

But \(L_{ik}\) belongs to \(\mathbf{J}_{11}\) and \(O_{ik}\) to \(\mathbf{J}_{22}\), while \(N_{ik}\) belongs to \(\mathbf{J}_{12}\) and \(M_{ik}\) to \(\mathbf{J}_{21}\). The two diagonal blocks are large and the two off-diagonal blocks are small — which is the shading in the figure of Section 20-4. The physical statement behind the algebra is the one every operator knows: real power follows angle, reactive power follows voltage magnitude. Governors redistribute megawatts; excitation and capacitors correct volts.

Set the small blocks to zero. The single linear system splits into two independent halves:

The price is paid in convergence, not in accuracy. The equations being solved are untouched; only the matrix used to compute the corrections has been approximated. The iteration therefore still converges to the exact solution of the exact load flow equations — the mismatch test is still applied to the full nonlinear equations — but it no longer converges quadratically, because the matrix is no longer the true derivative. Decoupled Newton typically needs one or two more iterations than full Newton, each much cheaper.

This is a general and important principle. In Newton's method the Jacobian controls only the route to the answer. The answer itself is fixed entirely by the mismatch equations. Any approximation to \(\mathbf{J}\) that still points roughly downhill gives the right result, more slowly. Everything in the next section is an aggressive application of that licence.

Decoupling removed two blocks but left \(\mathbf{J}_{11}\) and \(\mathbf{J}_{22}\) to be rebuilt and refactorised at every iteration. Stott and Alsac showed in 1974 that with a few further approximations both matrices become constant — independent of the voltages — so each can be factorised once, before the iteration begins, and reused unchanged to the end. That is the fast decoupled load flow, and it is still the method behind real-time contingency analysis.

Start from the elements of \(\mathbf{J}_{11}\) and apply, in order, three approximations that the same physics justifies:

The last step deserves a comment. \(B_{ii}\) is the sum of the susceptances of every branch at the bus, typically twenty to fifty per unit, so \(|V_i|^2|B_{ii}|\) is an order of magnitude larger than any bus's reactive injection. Dropping \(Q_i\) is therefore a genuine small correction and not wishful thinking.

Divide the \(i\)-th equation of \(\Delta\mathbf{P}=\mathbf{J}_{11}\Delta\boldsymbol{\delta}\) by \(|V_i|\), take a factor \(|V_k|\approx 1\) out of the off-diagonal terms, and what is left on the matrix side contains no voltages at all:

The two matrices are not identical, and the differences are deliberate. \(\mathbf{B}'\) governs the \(P\)–\(\delta\) half, where the influences to be suppressed are those that have nothing to do with active power transfer: line resistances (so branch susceptances are taken as \(1/x_{ik}\), not \(-\operatorname{Im}(1/z_{ik})\)), line charging, and shunt reactors and capacitors. \(\mathbf{B}''\) governs the \(Q\)–\(|V|\) half, where shunts matter a great deal — they are the very devices that supply reactive power — so it is built from the imaginary part of \(\mathbf{Y}_{bus}\) exactly as it stands, and its rows and columns cover only the PQ buses.

The iteration alternates half-steps rather than taking one combined step. Solve the \(P\)–\(\delta\) system, update the angles, recompute the reactive mismatches with those new angles, solve the \(Q\)–\(|V|\) system, update the magnitudes, recompute the real mismatches, and repeat. Using the freshly updated angles in the reactive half-step is the Gauss-Seidel idea of Chapter 19 reappearing at a higher level, and it measurably improves the convergence.

One approximation remains. In the fast decoupled \(P\)–\(\delta\) equation, suppose the voltage magnitudes are not merely nearly constant but taken as exactly \(1.0\) per unit everywhere, and the reactive half of the problem is abandoned altogether. Then \(\Delta P_i/|V_i| = \Delta P_i\), the matrix \(\mathbf{B}'\) is constant, and — since the relation is now exactly linear — no iteration is needed at all:

The DC load flow is not a competitor to the methods above; it answers a different question. It is used wherever an enormous number of cases must be screened and only megawatt flows matter: security-constrained dispatch, transmission planning studies, locational pricing in electricity markets, and the first pass of a contingency analysis that hands its worst cases to a full AC solution. Its errors on real networks are typically a few per cent in flow, which is acceptable for screening and unacceptable for anything voltage-related. On the test system of this chapter it returns \(\delta_2=-3.86^\circ\) against the true \(-2.70^\circ\) — the discrepancy comes almost entirely from the \(R/X\) ratio of \(0.5\), which is far outside the range the approximation was designed for.

Setting the whole family side by side on the same three-bus case, solved to the same tolerance of \(10^{-5}\) per unit from the same flat start:

The shapes matter more than the counts. Both Gauss-Seidel traces and the fast decoupled trace are straight, because a fixed factor per iteration is a fixed drop per iteration on a logarithmic scale. The Newton trace bends downward and keeps bending, because its factor improves as it goes. Extend the horizontal axis to a three-hundred-bus system and the straight lines flatten while the bend does not move — which is the whole argument of this chapter in one picture.

Work in per unit on a \(100\) MVA base with bus \(1\) as the slack, and start flat unless told otherwise. Angles in the Jacobian are in radians; convert only when reporting.