Section IV: Load Flow by Gauss-Seidel Method
The basic power flow equations (4.6) and (4.7) are nonlinear. In an n -bus power system, let the number of P-Q buses be np and the number of P-V (generator) buses be ng such that n = np + ng + 1. Both voltage magnitudes and angles of the P-Q buses and voltage angles of the P-V buses are unknown making a total number of 2np + ng quantities to be determined. Amongst the known quantities are 2np numbers of real and reactive powers of the P-Q buses, 2ng numbers of real powers and voltage magnitudes of the P-V buses and voltage magnitude and angle of the slack bus. Therefore there are sufficient numbers of known quantities to obtain a solution of the load flow problem. However, it is rather difficult to obtain a set of closed form equations from (4.6) and (4.7). We therefore have to resort to obtain iterative solutions of the load flow problem.
At the beginning of an iterative method, a set of values for the unknown quantities are chosen. These are then updated at each iteration. The process continues till errors between all the known and actual quantities reduce below a pre-specified value. In the Gauss-Seidel load flow we denote the initial voltage of the i th bus by Vi(0) , i = 2, ... , n . This should read as the voltage of the i th bus at the 0th iteration, or initial guess. Similarly this voltage after the first iteration will be denoted by Vi(1) . In this Gauss-Seidel load flow the load buses and voltage controlled buses are treated differently. However in both these type of buses we use the complex power equation given in (4.5) for updating the voltages. Knowing the real and reactive power injected at any bus we can expand (4.5) as
 | (4.11) |
We can rewrite (4.11) as
 | (4.12) |
In this fashion the voltages of all the buses are updated. We shall outline this procedure with the help of the system of Fig. 4.1, with the system data given in Tables 4.1 to 4.3. It is to be noted that the real and reactive powers are given respectively in MW and MVAr. However they are converted into per unit quantities where a base of 100 MVA is chosen.
Updating Load Bus Voltages
Let us start the procedure with bus-2 of the 5 bus 7 line system given in fig: 4.1. Since this is load bus, both the real and reactive power into this bus is known. We can therefore write from (4.12)
 | (4.13) |
From the data given in Table 4.3 we can write
 |
It is to be noted that since the real and reactive power is drawn from this bus, both these quantities appear in the above equation with a negative sign. With the values of the Y bus elements given in Table 4.2 we get V21 = 0.9927 < − 2.5959° .
The first iteration voltage of bus-3 is given by
 | (4.14) |
Note that in the above equation since the update for the bus-2 voltage is already available, we used the 1st iteration value of this rather than the initial value. Substituting the numerical data we get V3(1) = 0.9883 < − 2. 8258° . Finally the bus-4 voltage is given by
 | (4.15) |
Solving we get V4(1) = 0. 9968 < −3.4849° .
Updating P-V Bus Voltages
It can be seen from Table 4.3 that even though the real power is specified for the P-V bus-5, its reactive power is unknown. Therefore to update the voltage of this bus, we must first estimate the reactive power of this bus. Note from Fig. 4.11 that
 | (4.16) |
And hence we can write the kth iteration values as
 | (4.17) |
For the system of Fig. 4.1 we have
 | (4.18) |
This is computed as 0.0899 per unit. Once the reactive power is estimated, the bus-5 voltage is updated as
 | (4.19) |
It is to be noted that even though the power generation in bus-5 is 48 MW, there is a local load that is consuming half that amount. Therefore the net power injected by this bus is 24 MW and consequently the injected power P5, inj in this case is taken as 0.24 per unit. The voltage is calculated as V5(1) = 1.0169 < − 0.8894° . Unfortunately however the magnitude of the voltage obtained above is not equal to the magnitude given in Table 4.3. We must therefore force this voltage magnitude to be equal to that specified. This is accomplished by
 | (4.20) |
This will fix the voltage magnitude to be 1.02 per unit while retaining the phase of − 0.8894 ° . The corrected voltage is used in the next iteration.
Convergence of the Algorithm
As can be seen from Table 4.3 that a total number of 4 real and 3 reactive powers are known to us. We must then calculate each of these from (4.6) and (4.7) using the values of the voltage magnitudes and their angle obtained after each iteration. The power mismatches are then calculated from (4.9) and (4.10). The process is assumed to have converged when each of ΔP2 , ΔP3, ΔP4 , ΔP5 , ΔQ2 , ΔQ3 and ΔQ4 is below a small pre-specified value. At this point the process is terminated.
Sometimes to accelerate computation in the P-Q buses the voltages obtained from (4.12) is multiplied by a constant. The voltage update of bus- i is then given by
 | (4.21) |
where λ is a constant that is known as the acceleration factor . The value of λ has to be below 2.0 for the convergence to occur. Table 4.4 lists the values of the bus voltages after the 1st iteration and number of iterations required for the algorithm to converge for different values of λ. It can be seen that the algorithm converges in the least number of iterations when λ is 1.4 and the maximum number of iterations are required when λ is 2. In fact the algorithm will start to diverge if larger values of acceleration factor are chosen. The system data after the convergence of the algorithm will be discussed later.
Table 4.4 Gauss-Seidel method: bus voltages after 1 st iteration and number of iterations required for convergence for different values of l .
| l | Bus voltages (per unit) after 1st iteration | No of iterations for convergence | |||
| V2 | V3 | V4 | V5 | ||
| 1 | 0.9927Ð- 2.6° | 0.9883Ð- 2.83° | 0.9968Ð- 3.48° | 1.02 Ð- 0.89° | 28 |
| 2 | 0.9874Ð- 5.22° | 0.9766Ð- 8.04° | 0.9918Ð- 14.02° | 1.02Ð- 4.39° | 860 |
| 1.8 | 0.9883Ð- 4.7° | 0.9785Ð- 6.8° | 0.9903Ð- 11.12° | 1.02Ð- 3.52° | 54 |
| 1.6 | 0.9893Ð- 4.17° | 0.9807Ð- 5.67° | 0.9909Ð- 8.65° | 1.02Ð- 2.74° | 24 |
| 1.4 | 0.9903Ð- 3.64° | 0.9831Ð- 4.62° | 0.9926Ð- 6.57° | 1.02Ð- 2.05° | 14 |
| 1.2 | 0.9915Ð- 3.11° | 0.9857Ð- 3.68° | 0.9947Ð- 4.87° | 1.02Ð- 1.43° | 19 |
Section V: Solution of a Set of Nonlinear Equations by Newton-Raphson Method
In this section we shall discuss the solution of a set of nonlinear equations through Newton-Raphson method. Let us consider that we have a set of n nonlinear equations of a total number of n variables x1 , x2 , ... , xn. Let these equations be given by
 | (4.22) |
where f1, ... , fn are functions of the variables x1 , x2 , ... , xn. We can then define another set of functions g1 , ... , gn as given below
 | (4.23) |
Let us assume that the initial estimates of the n variables are x1(0) , x2(0) , ... , xn(0) . Let us add corrections Δx1(0) , Δx2(0) , ... , Δxn(0) to these variables such that we get the correct solution of these variables defined by
 | (4.24) |
The functions in (4.23) then can be written in terms of the variables given in (4.24) as
 | (4.25) |
We can then expand the above equation in Taylor 's series around the nominal values of x1(0) , x2(0) , ... , xn(0) . Neglecting the second and higher order terms of the series, the expansion of gk , k = 1, ... , n is given as
 | (4.26) |
where 
is the partial derivative of gk evaluated at x2(1) , ... , xn(1) .
is the partial derivative of gk evaluated at x2(1) , ... , xn(1) . Equation (4.26) can be written in vector-matrix form as
 | (4.27) |
The square matrix of partial derivatives is called the Jacobian matrix J with J (1) indicating that the matrix is evaluated for the initial values of x2(0) , ... , xn(0) . We can then write the solution of (4.27) as
 | (4.28) |
 , k= 1, ....n |
Since the Taylor 's series is truncated by neglecting the 2nd and higher order terms, we cannot expect to find the correct solution at the end of first iteration. We shall then have
 | (4.29) |
These are then used to find J (1) and Δgk (1) , k = 1, ... , n . We can then find Δx2(1) , ... , Δxn(1) from an equation like (4.28) and subsequently calculate x2(1) , ... , xn(1). The process continues till Δgk , k = 1, ... , n becomes less than a small quantity.
Section VI: Load Flow By Newton-Raphson Method
Let us assume that an n -bus power system contains a total np number of P-Q buses while the number of P-V (generator) buses be ng such that n = np + ng + 1. Bus-1 is assumed to be the slack bus. We shall further use the mismatch equations of ΔPi and ΔQi given in (4.9) and (4.10) respectively. The approach to Newton-Raphson load flow is similar to that of solving a system of nonlinear equations using the Newton-Raphson method: At each iteration we have to form a Jacobian matrix and solve for the corrections from an equation of the type given in (4.27). For the load flow problem, this equation is of the form
 | (4.30) |
where the Jacobian matrix is divided into submatrices as
 | (4.31) |
It can be seen that the size of the Jacobian matrix is ( n + np − 1) x ( n + np −1). For example for the 5-bus problem of Fig. 4.1 this matrix will be of the size (7 x 7). The dimensions of the submatrices are as follows:
J11: (n - 1) ´ (n - 1), J12: (n - 1) ´ np, J21: np ´ (n - 1) and J22: np ´ np
The submatrices are
 | (4.32) |
 | (4.33) |
 | (4.34) |
 | (4.35) |
Load Flow Algorithm
The Newton-Raphson procedure is as follows:
Step-1: Choose the initial values of the voltage magnitudes |V| (0) of all np load buses and n − 1 angles δ (0) of the voltages of all the buses except the slack bus.
Step-2: Use the estimated |V|(0) and δ (0) to calculate a total n − 1 number of injected real power Pcalc(0) and equal number of real power mismatch ΔP (0) .
Step-3: Use the estimated |V| (0) and δ (0) to calculate a total np number of injected reactive power Qcalc(0) and equal number of reactive power mismatch ΔQ (0) .
Step-3: Use the estimated |V| (0) and δ (0) to formulate the Jacobian matrix J (0) .
Step-4: Solve (4.30) for δ (0) and Δ |V| (0) ÷ |V| (0).
Step-5 : Obtain the updates from
 | (4.36) |
 | (4.37) |
Step-6: Check if all the mismatches are below a small number. Terminate the process if yes. Otherwise go back to step-1 to start the next iteration with the updates given by (4.36) and (4.37).
Formation of the Jacobian Matrix
We shall now discuss the formation of the submatrices of the Jacobian matrix. To do that we shall use the real and reactive power equations of (4.6) and (4.7). Let us rewrite them with the help of (4.2) as
 | (4.38) |
 | (4.39) |
A. Formation of J11
Let us define J11 as
 | (4.40) |
It can be seen from (4.32) that Lik's are the partial derivatives of Pi with respect to δk. The derivative Pi (4.38) with respect to k for i ≠ k is given by
 | (4.41) |
Similarly the derivative Pi with respect to k for i = k is given by
 |
Comparing the above equation with (4.39) we can write
 | (4.42) |
B. Formation of J21
Let us define J21 as
 | (4.43) |
From (4.34) it is evident that the elements of J21 are the partial derivative of Q with respect to δ . From (4.39) we can write
 | (4.44) |
Similarly for i = k we have
 | (4.45) |
The last equality of (4.45) is evident from (4.38).
Formation of the Jacobian Matrix
C. Formation of J12
Let us define J12 as
 | (4.46) |
As evident from (4.33), the elements of J21 involve the derivatives of real power P with respect to magnitude of bus voltage |V| . For i ≠ k , we can write from (4.38)
 | (4.47) |
For i = k we have
 | (4.48) |
Formation of J22
For the formation of J22 let us define
 | (4.49) |
For i ≠ k we can write from (4.39)
 | (4.50) |
Finally for i = k we have
 | (4.51) |
We therefore see that once the submatrices J11 and J21 are computed, the formation of the submatrices J12 and J22 is fairly straightforward. For large system this will result in considerable saving in the computation time.
Solution of Newton-Raphson Load Flow
The Newton-Raphson load flow program is tested on the system of Fig. 4.1 with the system data and initial conditions given in Tables 4.1 to 4.3. From (4.41) we can write
 |
Similarly from (4.39) we have
 |
Hence from (4.42) we get
 |
In a similar way the rest of the components of the matrix J11(0) are calculated. This matrix is given by
 |
For forming the off diagonal elements of J21 we note from (4.44) that
 |
Also from (4.38) the real power injected at bus-2 is calculated as
 |
Hence from (4.45) we have
 |
Similarly the rest of the elements of the matrix J21 are calculated. This matrix is then given as
 |
For calculating the off diagonal elements of the matrix J12 we note from (4.47) that they are negative of the off diagonal elements of J21 . However the size of J21 is (3 X 4) while the size of J12 is (4 X 3). Therefore to avoid this discrepancy we first compute a matrix M that is given by
 |
The elements of the above matrix are computed in accordance with (4.44) and (4.45). We can then define
 |
Furthermore the diagonal elements of J12 are overwritten in accordance with (4.48). This matrix is then given by
 |
Finally it can be noticed from (4.50) that J22 = J11 (1:3, 1:3). However the diagonal elements of J22 are then overwritten in accordance with (4.51). This gives the following matrix
 |
From the initial conditions the power and reactive power are computed as
 |
Consequently the mismatches are found to be
 |
Then the updates at the end of the first iteration are given as
 |
The load flow converges in 7 iterations when all the power and reactive power mismatches are below 10−6 .
Section VII: Load Flow Results
In this section we shall discuss the results of the load flow. It is to be noted here that both Gauss-Seidel and Newton-Raphson methods yielded the same result. However the Newton-Raphson method converged faster than the Gauss-Seidel method. The bus voltage magnitudes, angles of each bus along with power generated and consumed at each bus are given in Table 4.4. It can be seen from this table that the total power generated is 174.6 MW whereas the total load is 171 MW. This indicates that there is a line loss of about 3.6 MW for all the lines put together. It is to be noted that the real and reactive power of the slack bus and the reactive power of the P-V bus are computed from (4.6) and (4.7) after the convergence of the load flow.
Table 4.4 Bus voltages, power generated and load after load flow convergence.
| Bus no. | Bus voltage | Power generated | Load | |||
| Magnitude (pu) | Angle (deg) | P (MW) | Q (MVAr) | P (MW) | P (MVAr) | |
| 1 | 1.05 | 0 | 126.60 | 57.11 | 0 | 0 |
| 2 | 0.9826 | - 5.0124 | 0 | 0 | 96 | 62 |
| 3 | 0.9777 | - 7.1322 | 0 | 0 | 35 | 14 |
| 4 | 0.9876 | - 7.3705 | 0 | 0 | 16 | 8 |
| 5 | 1.02 | - 3.2014 | 48 | 15.59 | 24 | 11 |
The current flowing between the buses i and k can be written as
 | (4.52) |
Therefore the complex power leaving bus- i is given by
 | (4.53) |
Similarly the complex power entering bus- k is
 | (4.54) |
Therefore the I 2 R loss in the line segment i-k is
 | (4.55) |
The real power flow over different lines is listed in Table 4.5. This table also gives the I2 R loss along various segments. It can be seen that all the losses add up to 3.6 MW, which is the net difference between power generation and load. Finally we can compute the line I2X drops in a similar fashion. This drop is given by
 | (4.56) |
However we have to consider the effect of line charging separately.
Table 4.5 Real power flow over different lines.
| Power dispatched | Power received | Line loss (MW) | ||
| from (bus) | amount (MW) | in (bus) | amount (MW) | |
| 1 | 101.0395 | 2 | 98.6494 | 2.3901 |
| 1 | 25.5561 | 5 | 25.2297 | 0.3264 |
| 2 | 17.6170 | 3 | 17.4882 | 0.1288 |
| 3 | 0.7976 | 4 | 0.7888 | 0.0089 |
| 5 | 15.1520 | 2 | 14.9676 | 0.1844 |
| 5 | 18.6212 | 3 | 18.3095 | 0.3117 |
| 5 | 15.4566 | 4 | 15.2112 | 0.2454 |
| | | | | Total = 3.5956 |
Load Flow Results
Consider the line segment 1-2. The voltage of bus-1 is V1 = 1.05 <>V2 = 0.9826 < − 5.0124° per unit. From (4.52) we then have
 per unit |
Therefore the complex power dispatched from bus-1 is
 |
where the negative signal indicates the power is leaving bus-1. The complex power received at bus-2 is
 MW |
Therefore out of a total amount of 101.0395 MW of real power is dispatched from bus-1 over the line segment 1-2, 98.6494 MW reaches bus-2. This indicates that the drop in the line segment is 2.3901 MW. Note that
 MW |
where R12 is resistance of the line segment 1-2. Therefore we can also use this method to calculate the line loss.
Now the reactive drop in the line segment 1-2 is
 MV Ar |
We also get this quantity by subtracting the reactive power absorbed by bus-2 from that supplied by bus-1. The above calculation however does not include the line charging. Note that since the line is modeled by an equivalent- p , the voltage across the shunt capacitor is the bus voltage to which the shunt capacitor is connected. Therefore the current I 12 flowing through line segment is not the current leaving bus-1 or entering bus-2 - it is the current flowing in between the two charging capacitors. Since the shunt branches are purely reactive, the real power flow does not get affected by the charging capacitors. Each charging capacitor is assumed to inject a reactive power that is the product of the half line charging admittance and square of the magnitude of the voltage of that at bus. The half line charging admittance of this line is 0.03. Therefore line charging capacitor will inject
 MV Ar |
at bus-1. Similarly the reactive injected at bus-2 will be
 MV Ar |
The power flow through the line segments 1-2 and 1-5 are shown in Fig. 4.2.
(a)
(b)
Fig. 4.2 Real and reactive power flow through (a) line segment 1-2 and (b) line segment 1-5. The thin lines indicate reactive power flow while the thick lines indicate real power flow.
Section VIII: Load Flow Programs In Matlab
The load flow programs are developed in MATLAB. Altogether there are 4 mfiles that are attached with this chapter. The program listings and descriptions of these mfiles are given below. It must however be emphasized that these are not general purpose programs and are written only for the examples of this chapter.
- Forming Ybus Matrix
- Gauss- Seidel Load Flow
- Solving Nonlinear Equations using Newton- Raphson
- Newton - Raphson Load Flow
Forming Ybus Matrix
This is a function that can be called by various programs. The function can be invoked by the statement
[yb,ych]=ybus;
where 'yb' and 'ych' are respectively the Ybus matrix and a matrix containing the line charging admittances. It is assumed that the system data of Table 4.1 are given in matrix form and the matrix that contains line impedances is 'zz', while 'ych' contains the line charging information. This program is stored in the file ybus.m. The program listing is given below.
| % Function ybus function [yb,ych]=ybus zz=[0 0.02+0.1i 0 0 0.05+0.25i 0.02+0.1i 0 0.04+0.2i 0 0.05+0.25i 0 0 0.05+0.25i 0 0.1+0.5i % The line chargings are ych=j*[0 0.03 0 0 0.02 |
Gauss-Seidel Load Flow
The Gauss-Seidel program is stored in the file loadflow_gs.m. This calls the ybus.m function discussed above. The program allows the selection of the acceleration factor. The program lists the number of iterations required to converge, bus voltages and their magnitudes and real and reactive power. The program listing is given below.
% Program loadflow_gs
d2r=pi/180;w=100*pi;
% The Y_bus matrix is
[ybus,ych]=ybus;
% The given parameters and initial conditions are
p=[0;-0.96;-0.35;-0.16;0.24];
acc=input('Enter the acceleration constant: ');
% The Gauss-Seidel iterations starts here
while del>1e-6
% P-Q buses
for i=2:4
f
tmp2=tmp2+0;
vt=(tmp1-tmp2)/ybus(i,i);
v(i)=v(i)+acc*(vt-v(i));
v(i)=v(i)+acc*(vt-v(i));
tmp1=(p(5)-j*q5)/conj(v(5));
delq=q+imag(s)';
delpq=[delp(2:5);delq(2:4)];
'GS LOAD FLOW CONVERGES IN ITERATIONS',indx,pause
'FINAL VOLTAGE MAGNITUDES ARE',abs(v)',pause
'THE REACTIVE POWERS IN EACH BUS IN MVar ARE',(-imag(s)+[0 0 0 0 0.11])*100'FINAL VOLTAGE MAGNITUDES ARE',abs(v)',pause
Solving Nonlinear Equations using Newton-Raphson
This program gives the solution of the nonlinear equations of Example 4.1 using the Newton-Raphson method. The equations and the Jacobian matrix are explicitly entered in the program itself. The program gives the number of iterations and the final values of x1 , x2 and x3 . The program listing is given below.
% Program nwtraph
clear all
x=[1;1;1];
% The Newton-Raphson iterations starts here
del=1;
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,pause
'FINAL VALUES OF x ARE',x Newton-Raphson Load Flow
The Newton-Raphson load flow program is stored in the files loadflow_nr.m. The outputs of the program can be checked by typing
indx the number of iterations v bus voltages in Cartesian form abs(v) magnitude of bus voltages angle(v)/d2r angle of bus voltage in degree preal real power in MW preac reactive power in MVAr pwr power flow in the various line segments qwr reactive power flow in the various line segments q reactive power entering or leaving a bus pl real power losses in various line segments ql reactive drops in various line segments |
It is to be noted that in calculating the power and reactive power the conventions that the power entering a node is positive and leaving it is negative are maintained. The program listing for the Newton-Raphson load flow is given below.
% Program loadflow_nr
clear all
d2r=pi/180;w=100*pi;
% The Ybus matrix is
[ybus,ych]=ybus;
g=real(ybus);b=imag(ybus);
% The given parameters and initial conditions are
p=[0;-0.96;-0.35;-0.16;0.24];
del=1;indx=0;
% The Newton-Raphson iterations starts here
d
% The Jacobian matrix
j11(i,k)=mv(ii)*mv(kk)*(g(ii,kk)*sin(th(ii)-th(kk))-b(ii,kk)*cos(th(ii)-th(kk)));
end
for i=1:3
delpq=[delp(2:5);delq(2:4)];
corr=inv(jacob)*delpq;
th=th+[0;corr(1:4)];
mv=mv+[0;mv(2:4).*corr(5:7);0];
del=max(abs(delpq));
indx=indx+1;
end
preal=(pcal+[0 0 0 0 0.24])*100;
preac=(qcal+[0 0 0 0 0.11])*100;
% Power flow calculations
for i=1:5
for i=1:4
c(k,i)=100*abs(ych(k,i))*abs(v(k))^2;
cur(i,k)=cu;cur(k,i)=-cur(i,k);
c(k,i)=100*abs(ych(k,i))*abs(v(k))^2;
cur(i,k)=cu;cur(k,i)=-cur(i,k);
pwr=real(s);
q=qwr-c;
% Power loss
ilin=abs(cur);
for i=1:4pl(i,k)=100*r*ilin(i,k)^2;pl(k,i)=pl(i,k);
