Gauss-Seidel Load Flow

Chapter 18 ended with a set of nonlinear algebraic equations and the plain statement that they cannot be solved in closed form beyond two buses. What remains is to guess an answer, use the equations to produce a better guess, and repeat. Every load flow program ever written does exactly this; the methods differ only in how the improvement is manufactured.

The oldest device is the simplest. Take a single equation \(f(x)=0\) and rearrange it, by whatever algebra is convenient, into the form

A solution \(x^{*}\) is now a value that the function \(g\) returns unchanged — a fixed point. Start from any estimate \(x^{(0)}\) and compute \(x^{(1)}=g(x^{(0)})\), then \(x^{(2)}=g(x^{(1)})\), and so on. This is successive substitution, or the Gauss iterative method. If the sequence settles down at all, it settles at a fixed point, because the moment \(x^{(k+1)}\) equals \(x^{(k)}\) the defining equation is satisfied exactly.

Whether it settles is a separate question, and it has a clean answer. Let \(e^{(k)} = x^{(k)}-x^{*}\) be the error. Then

Each pass multiplies the error by the slope of \(g\) at the solution. The iteration converges if and only if that slope is smaller than one in magnitude, and the smaller it is the faster the error dies. The decay is geometric: the error is multiplied by a constant factor every pass, so it takes a fixed number of iterations to gain each decimal place. This is called linear convergence, and it is the essential character of the entire Gauss-Seidel family.

The geometry makes the criterion obvious. Plot \(y=g(x)\) and the line \(y=x\) on the same axes; the solution is where they cross. Starting from \(x^{(0)}\), move vertically to the curve to get \(g(x^{(0)})\), then horizontally to the line \(y=x\) to make that value the new \(x\), and repeat. The path is a staircase. If the curve is flatter than the \(45^\circ\) line the staircase closes in on the crossing; if it is steeper the staircase walks away.

Move from one equation to \(n\) of them. Suppose each equation has been solved for one variable, so that we have a set of update rules \(x_i = g_i(x_1,\dots,x_n)\), one per variable. There are two ways to run them.

The Gauss (or Jacobi) method computes every new value from the old ones, keeping the whole vector frozen until all \(n\) updates are complete:

The Gauss-Seidel method does the obvious thing instead. By the time variable \(i\) is being updated, variables \(1\) through \(i-1\) have already been improved in this very sweep. Throwing that improvement away is wasteful, so use it:

The gain is real and can be quantified for a linear system. Write \(\mathbf{A}\mathbf{x}=\mathbf{b}\) and split \(\mathbf{A}=\mathbf{L}+\mathbf{D}+\mathbf{U}\) into its strictly lower, diagonal and strictly upper parts. Solving the \(i\)-th equation for \(x_i\) gives \(\mathbf{D}\mathbf{x} = \mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}\), and the two schemes become

The error is multiplied at each pass by the largest eigenvalue in magnitude — the spectral radius \(\rho\) — of the corresponding matrix, which is the multivariable version of the slope \(g'\) of Section 19-1. For a large class of matrices, including those arising from networks, the Gauss-Seidel spectral radius is the square of the Gauss one. Squaring a number smaller than one halves the number of iterations needed. Example 1 shows this happening digit for digit.

Both schemes converge for certain if \(\mathbf{A}\) is strictly diagonally dominant, meaning that in every row the diagonal entry outweighs the sum of the others:

This is where the load flow problem hands the method a gift. Chapter 16 built \(\mathbf{Y}_{bus}\) by the rule that the diagonal entry is the sum of every admittance touching the bus and the off-diagonal entry is the negative of the admittance joining two buses. Therefore

with strict inequality wherever a shunt \(y_{i0}\) is present — line charging, a shunt capacitor, a transformer magnetising branch. \(\mathbf{Y}_{bus}\) is diagonally dominant by construction, not by luck, and that single structural fact is the reason a method as unsophisticated as Gauss-Seidel works on power networks at all. It is also the reason the method is built on \(\mathbf{Y}_{bus}\) rather than on the \(\mathbf{Z}_{bus}\) of Chapter 17, which is full and has no such structure.

Chapter 18 wrote the complex power injected at bus \(i\) in conjugated form so that the admittance sum was left intact:

The variable we want to isolate, \(V_i\), sits inside the sum multiplied by \(Y_{ii}\) — the largest coefficient in the row, by the dominance argument just made. Divide through by \(V_i^{*}\), separate the \(k=i\) term, and solve for \(V_i\):

This is the required \(x=g(x)\). The right-hand side still contains \(V_i\) through the conjugate \(V_i^{*}\), which is precisely what makes it a fixed-point equation rather than a formula. Read it physically and it is nothing but Kirchhoff's current law at bus \(i\), solved for the voltage that would make the current drawn by the specified power balance the current arriving through the network.

Two points of bookkeeping decide whether a student's arithmetic comes out right, and both were settled in Chapter 18. First, \(P_{i,\text{inj}}\) and \(Q_{i,\text{inj}}\) are net injections: generation minus load at that bus. A pure load bus injects negative power, so a load of \(4.0+j2.5\) per unit enters the formula as \(P_i=-4.0\), \(Q_i=-2.5\). A bus with both a machine and a local load carries the difference. Second, everything is in per unit on the common base of Chapter 4, so a \(400\) MW load on a \(100\) MVA base is \(4.0\), not \(400\).

An equivalent form is sometimes more convenient for hand calculation because it avoids forming \(\mathbf{Y}_{bus}\) explicitly. Writing \(Y_{ii}=y_{i0}+\sum_{k\neq i}y_{ik}\) and \(Y_{ik}=-y_{ik}\),

The two are the same equation written twice; the second makes the physical reading immediate. The numerator is the total current entering bus \(i\) — the load current plus everything arriving along the branches — and the denominator is the total admittance seen at the bus. The quotient is a voltage.

At a PQ bus both \(P_i\) and \(Q_i\) are known numbers, so the update of Section 19-3 can be applied directly with nothing left to estimate. This is the easy case and it is worth writing out in full for a small system before adding the complication of generator buses.

Number the buses so that bus \(1\) is the slack. One iteration means one complete sweep over buses \(2\) through \(n\), taken in order. For a five-bus system with buses \(2,3,4\) of PQ type, the first sweep from a flat start reads

The underlined entries are the Seidel refinement in action: by the time bus \(4\) is reached, the improved voltages at buses \(2\) and \(3\) are already in memory and are used. The slack voltage \(V_1\) never carries an iteration superscript because it is fixed data. And the conjugate \(V_i^{*(k)}\) in the numerator always uses the previous value of the bus's own voltage, since the new one does not exist yet — that is the fixed-point structure of Section 19-1.

Everything in this chapter is worked on the three-bus system that Chapter 18 introduced and solved, so that every intermediate number can be checked against a known answer.

Its admittance matrix was formed in Chapter 18 and is repeated here because every worked example needs it:

Every row obeys the dominance inequality of Section 19-2 with equality, since no shunts are present. Row two, for instance, gives \(|Y_{22}| = |26-j52| = 58.14\) against \(|Y_{21}|+|Y_{23}| = 22.36+35.78 = 58.14\). The margin is exactly zero, which warns in advance that convergence will be steady rather than quick.

At a generator bus under automatic voltage regulator control, Chapter 18 established that \(P_i\) and \(|V_i|\) are specified and \(Q_i\) and \(\delta_i\) are not. The update of Section 19-3 demands \(Q_i\), which we do not have. It also produces a complete complex number, magnitude included, when only the angle is wanted. Both difficulties are handled inside the sweep, and the two fixes together are the only structural difference between a PQ bus and a PV bus in Gauss-Seidel.

Step one — estimate the reactive injection. The reactive power flow equation of Chapter 18 is an identity that holds at the solution. Evaluate it at the present iterate and it returns the reactive power that the current voltages imply:

This is not the true \(Q_i\) — the voltages are not yet the true voltages — but it is the best estimate available, and it becomes exact as the iteration converges. The bus's own voltage enters at its old value while the neighbours already visited enter at their new ones, keeping faith with the Seidel principle.

Step two — correct the magnitude. Feed that \(Q\) into the ordinary voltage update. The result is a complex number \(V_i^{(k+1)}\) whose magnitude will generally miss the scheduled value, because nothing in the update was told about the schedule. The angle, however, is genuine new information. Keep the angle and impose the magnitude:

The correction is legitimate for a reason worth stating plainly. The load flow at a PV bus has exactly one unknown, \(\delta_i\). The update formula, being a complex equation, delivers two real numbers. Discarding the magnitude discards a quantity that was never an unknown; retaining the angle keeps the quantity that was. What the discarded magnitude actually measures is the residual error in the reactive estimate, and watching it approach the scheduled value is a useful convergence check in its own right.

Step three — check the reactive limit. Chapter 18 closed with the physical fact that a machine's reactive output is bounded above by field heating and below by end-region heating and steady-state stability. The estimate \(Q_{i}^{(k+1)}\) must be tested against those bounds at every iteration:

When a limit binds, the clamped value is used and — this is the part students forget — the magnitude correction is skipped. The bus has become a PQ bus with its reactive injection pinned at the limit, and its voltage magnitude is now an unknown free to settle wherever the network puts it. The scheduled magnitude is remembered so that the bus can be restored to PV type if a later iteration finds the constraint no longer binding. Example 6 carries this through on the test system.

The linear convergence of Section 19-1 has one redeeming feature: it is predictable. If each pass multiplies the error by roughly the same factor \(\rho\), then successive iterates march toward the solution along a nearly straight path, each step a fixed fraction of the one before. A method that knows it will keep travelling in the same direction should travel further while it is there.

That is the whole of the acceleration idea. Let \(V_i^{(k)}\) be the value the update formula produces and \(V_{i,\text{acc}}^{(k-1)}\) the accelerated value currently held. Instead of accepting the step, stretch it by a factor \(\lambda\):

The arithmetic of why over-shooting helps is one line. Suppose the unaccelerated iteration reduces the error by a factor \(\rho\), so that \(V^{(k)}-V^{*} = \rho\,(V_{\text{acc}}^{(k-1)}-V^{*})\). Substituting into the acceleration formula,

With \(\lambda=1\) the factor is \(\rho\), as it must be. Increasing \(\lambda\) drives the bracket toward zero and then past it; the value \(\lambda=1/(1-\rho)\) would annihilate the error in a single step. That ideal is unattainable in practice for two reasons, and both are instructive. The convergence ratio \(\rho\) is not known before the study is run. More fundamentally, a system with many buses has many error modes, each decaying at its own rate; a \(\lambda\) tuned to eliminate the slowest mode multiplies the bracket for the faster modes by something larger than one and amplifies them. The practical optimum therefore sits well below \(1/(1-\rho)\).

Past \(\lambda = 2\) the bracket exceeds one in magnitude for every \(\rho\) between \(0\) and \(1\), and the iteration cannot converge at all. This is the origin of the standard rule.

On the three-bus system of this chapter, \(\lambda=1.2\) reduces the iteration count from twelve to seven, \(\lambda=1.6\) is worse than no acceleration at all, and at \(\lambda=2.1\) the iteration diverges outright. A three-bus network is far stiffer than a real one and its optimum is unusually low, but the shape of the result — a modest gain, a clear optimum, rapid deterioration on the high side — is what every system shows.

The straightness of both traces is the visible signature of linear convergence. Each iteration multiplies the mismatch by the same factor — \(0.411\) without acceleration on this system — so each buys the same number of decimal places, and the cost of one more digit is always the same handful of sweeps. Chapter 20 shows what a curve that bends downward looks like.

Two stopping tests are in circulation and they are not equally good. The weaker one watches the voltages and stops when they cease to move:

It is cheap, since the differences are already computed, but it can lie. An iteration crawling along with a convergence ratio close to one produces tiny voltage changes while the equations are still badly violated, and the program reports success on a wrong answer. Worse, the danger is greatest exactly where it matters — a heavily loaded system, where the ratio is worst.

The honest test is the mismatch test of Chapter 18, applied to the equations themselves. After each sweep, evaluate the power flow equations at the current voltages, subtract from the scheduled values, and examine the residuals:

How many iterations should be expected? Gauss-Seidel has an unwelcome scaling property. The error modes of the sweep correspond to patterns of voltage spread over the network, and the slowest mode is the one that has to propagate information from one end of the system to the other. A single sweep moves information only one bus along in the direction opposite to the ordering, so the number of iterations rises roughly in proportion to the number of buses.

Compare this with the Newton-Raphson result of Chapter 20, where the iteration count is three to five almost regardless of size. That contrast — not the cost of a single iteration — is why Gauss-Seidel was displaced as the production method in the late nineteen-sixties.

When the iteration does not converge, three quite different things may be happening, and telling them apart saves a great deal of wasted effort.

Everything of Sections 19-3 to 19-7 assembles into a procedure short enough to program in an afternoon, which was a large part of its original appeal.

Three implementation habits separate a program that works from one that only sometimes works. Store \(\mathbf{Y}_{bus}\) sparsely, keeping only the non-zero entries, since a real network has three or four branches per bus and the matrix is more than \(99\%\) empty. Precompute the constants \(1/Y_{ii}\) once rather than dividing inside the loop. And carry the voltages in rectangular form throughout, converting to polar only for output — the update formula is entirely rectangular, and repeated polar conversion is both slow and a source of angle-wrapping bugs.

The method's reputation has suffered from being compared only against Newton-Raphson on large transmission cases, which is the one contest it is certain to lose. Set against the criteria that actually matter in a given application, it retains a clear set of strengths.

The last two rows explain why the method never quite disappeared. Its insensitivity to the initial estimate makes it the standard way to begin a difficult case: run one or two Gauss-Seidel sweeps from a flat start to pull the voltages into the neighbourhood of the solution, then hand over to Newton-Raphson, which converges quadratically once it is close. Almost every serious load flow program offers this hybrid.

Two further niches are worth knowing. Distribution feeders — long, radial, with \(R/X\) ratios near unity — defeat the fast decoupled method of Chapter 20, whose central approximation requires \(X \gg R\), while Gauss-Seidel and its relatives handle them without complaint. And in a system with only a handful of buses, twelve trivial sweeps beat three expensive Jacobian factorisations on any measure, which is why every examination question on hand-computed load flow is a Gauss-Seidel question.

A refinement worth naming is the \(\mathbf{Z}_{bus}\) Gauss-Seidel method, which uses the impedance matrix of Chapter 17 in place of \(\mathbf{Y}_{bus}\). Because \(\mathbf{Z}_{bus}\) is full, every bus feels every other bus in a single sweep and the iteration count collapses to a handful. The price is the storage and the formation cost of a full \(n\times n\) matrix, which is why the idea belongs to fault analysis in Chapter 21, where \(\mathbf{Z}_{bus}\) has to be built anyway, rather than to routine load flow.

Work in per unit on a \(100\) MVA base, take bus \(1\) as the slack, and start flat unless told otherwise. Carry five decimal places in intermediate results; rounding early destroys the later iterations.