Power System Operations

Load Flow Analysis

© Aktaibi
 Lecture outlines and outcomes:

• In this lecture you will study about:

• The load Flow analysis
To other
Bus K buses

𝑃𝑘 𝑄𝑘

𝑃𝐺𝑘 𝑄𝐺𝑘 𝑃𝐿𝑘 𝑄𝐿𝑘

P Q V δ
1. Slack bus x x ✓ ✓ Gen Load
2. PQ bus ✓ ✓ x x Vk= Vk /k
3. PV bus ✓ x ✓ x
© Aktaibi
 Power Flow:

• Successful power system operation under normal balanced three-phase steady-state

conditions requires the following:

1. Generation supplies the demand (load) plus losses

2. Bus voltage magnitudes remain close to rated values
3. Generators operate within specified real and reactive power limits
4. Transmission lines and transformers are not overloaded
© Aktaibi
 The Power-flow Problem
• The power-flow problem is:
• The computation of voltage magnitude and phase
angle at each bus in a power system under
balanced three-phase steady-state conditions.
• The calculation of real and reactive power flows in
equipment such as transmission lines and
transformers, as well as equipment losses, can be
computed.
• The starting point for a power-flow problem is a
single-line diagram of the power system, from
single-line diagram
which the input data can be obtained. Input data
consist of bus data, transmission line data, and
transformer data.
© Aktaibi
 The Power-flow Problem
• As shown in the Figure, the following
To other
four variables are associated with each buses
bus 𝑘:
• voltage magnitude 𝑉𝑘 , 𝑃𝑘 𝑄𝑘
• phase angle 𝑘 , Bus K
• real power 𝑃𝑘 ,
𝑃𝐺𝑘 𝑄𝐺𝑘 𝑃𝐿𝑘 𝑄𝐿𝑘
• reactive power 𝑄𝑘
• At each bus, two of these variables are Gen Load
specified as input data, and the other
two are unknowns to be computed by Vk= Vk /k
the power-flow program.
• For convenience, the power delivered P Q V δ
to bus 𝑘 in the Figure is separated into 1. Slack bus (reference bus) x x ✓ ✓
generator and load terms. That is,
𝑃𝑘 = 𝑃𝐺𝑘 − 𝑃𝐿𝑘 2. PQ bus (Load Bus) ✓ ✓ x x
𝑄𝑘 = 𝑄𝐺𝑘 − 𝑄𝐿𝑘 3. PV bus (voltage-controlled bus) ✓ x ✓ x
© Aktaibi
 Types of buses in Power flow analysis
Load flow analysis consists of calculating voltage magnitude and phase angle at the buses of
the power system, as well as the active and reactive power flows for the specified terminal (or
bus conditions). Four variables are associated with each bus:
• voltage |V|
• phase angle |δ|
• active power |P|
• reactive power |Q|

FYI
• Maximum and minimum VAR limits 𝑄𝐺𝑘 max and
𝑄𝐺𝑘 min that this equipment can supply are also
input data.
• If no generation, 𝑃𝑘 = − 𝑃𝐿𝑘 is negative. P Q V δ
• if the load is inductive, 𝑄𝑘 = − 𝑄𝐿𝑘 is negative. 1. Slack bus x x ✓ ✓
• A bus may be classified into the above- 2. PQ bus ✓ ✓ x x
mentioned categories as in this table: 3. PV bus ✓ x ✓ x
© Aktaibi
 Types of buses in Power flow analysis
Each bus k is categorized into one of the following three bus types:
To other
1. Swing bus (or slack bus)—There should be only one swing bus, which buses
for convenience is numbered bus 1 in this text. The swing bus is a
reference bus for which V1 /1, unless otherwise stated, it is typically 1/0o 𝑃𝑘 𝑄𝑘
per unit, is input data. The power-flow program computes 𝑃1 and 𝑄1. Bus K

2. Load bus (PQ) —𝑃𝑘 and 𝑄𝑘 are input data. The power-flow program 𝑃𝐺𝑘 𝑄𝐺𝑘 𝑃𝐿𝑘 𝑄𝐿𝑘
computes 𝑉𝑘 and 𝑘 . Most buses in a typical power-flow program are
load buses. Gen Load

3. Voltage controlled bus (PV) —𝑃𝑘 and 𝑉𝑘 are input data. Vk= Vk /k
The power flow program computes 𝑄𝑘 and 𝑘 . Examples are buses to
which generators, switched shunt capacitors, Tap-Changer transformer,
or static VAR systems are connected. Maximum and minimum VAR limits P Q V δ
𝑄𝐺𝑘 𝑚𝑎𝑥 and 𝑄𝐺𝑘 𝑚𝑖𝑛 that this equipment can supply are also input data. If
1. Slack bus x x ✓ ✓
an upper or lower reactive power limit is reached, then the reactive
power output of the generator is held at the limit, and the bus is
2. PQ bus ✓ ✓ x x
modeled as a PQ bus. 3. PV bus ✓ x ✓ x
© Aktaibi
 Power flow analysis
Note that:
𝑃𝑘 = 𝑃𝐺𝑘 − 𝑃𝐿𝑘
𝑄𝑘 = 𝑄𝐺𝑘 − 𝑄𝐿𝑘
• when bus 𝑘 is a load bus with no generation, 𝑃𝑘 = − 𝑃𝐿𝑘 is negative; that is, the real power
supplied to bus k is negative, and, if the load is inductive, 𝑄𝑘 = − 𝑄𝐿𝑘 is negative.
• The net complex power flow into the network is not known in advance, and the system power losses
are unknown until the study is complete.
• It is necessary to have only one slack bus (i.e., the slack bus as a reference) at which:
• complex power is unspecified so that it supplies the difference in the total system load plus losses and the
sum of the complex powers specified at the remaining buses.
• The slack bus must also be a generator bus.
• The complex power (P, Q) allocated to this bus is computed as part of the solution.
• In order for the variations in real and reactive powers of the slack bus to be a small percentage of its
generating capacity during the solution process, the bus connected to the largest generating station is
normally selected as the slack bus.
© Aktaibi
 How to solve the power flow problem:
• The solution requires mathematical formulation and numerical solution.
• Since load flow problems generate non-linear equations that computers cannot solve
quickly, numerical methods are required.
• The following methods are commonly used algorithms:

• Gauss Iterative Method

• Fast Decoupled Load Flow Method
• Gauss-Seidel Method
• Newton-Raphson Method
© Aktaibi
 Example 1: Solving an equation using a Numerical method
• Solve 𝑥− 𝑥−1=2
• Since we have one equation, we cannot solve this equation with the regular solving techniques. Instead we
need to use numerical analysis to solve it.
• We can use Gauss Iterative Method to solve this problem numerically. That means, assume a random value for
x and solve , and then use the new value for x , and so on.
clc
clear all
• 𝑥 (𝑣+1) = 1 + 𝑥 (𝑣) disp(' -----------------')
disp(' k x ')
disp(' -----------------')
for k=1:9
x(1)=1;
x(k+1)=1+sqrt(x(k));
disp([ k x(k) ])
end

-------------------
k x
-------------------
1 1
2 2
3.0000 2.4142
• So the value of 𝑥 = 2.61798 4.0000
5.0000
2.5538
2.5981
6.0000 2.6118
• to verify 𝑥− 𝑥−1=0 7.0000 2.6161
8.0000 2.6174
• 2.61798 − 1.618017 − 1 = −0.00009 (very small value) 9.0000 2.6179
© Aktaibi
 Example 2:
• Solving equations using Gauss-Seidel Method clc
clear all
4𝑥1 + 𝑥2 − 𝑥3 = 3 disp(' -----------------------------------------')
disp(' k x1 x2 x3 ')
2𝑥1 + 7 𝑥2 + 𝑥3 = 19 disp(' -----------------------------------------')
for k=1:11
𝑥1 − 3 𝑥2 + 12 𝑥3 = 31 x1(1)=0;
x2(1)=0;
Initial values are x1(0) = x2(0) = x3(0) = 0 x3(1)=0;
x1(k+1)=(1/4)*(3-x2(k)+x3(k));
1 x2(k+1)=(1/7)*(19-2*x1(k+1)-x3(k));
• 𝑥1 = 3 − 𝑥2 + 𝑥3 x3(k+1)=(1/12)*(31-x1(k+1)+3*x2(k+1));
4 disp([ k x1(k) x2(k) x3(k)])
1 end
• 𝑥2 = (19 − 2𝑥1 − 𝑥3)
7
1
• 𝑥3 = (31 − 𝑥1 + 3𝑥2) -----------------------------------------
12
k x1 x2 x3
-----------------------------------------
1 0 0 0

2.0000 0.7500 2.5000 3.1458

3.0000 0.9115 2.0045 3.0085

4.0000 1.0010 1.9985 2.9995

5.0000 1.0003 2.0000 3.0000

6.0000 1.0000 2.0000 3.0000

© Aktaibi
 Example 3: Find Y bus for this system

𝑦1 + 𝑦12 + 𝑦13 −𝑦12 −𝑦13

−𝑦12
𝑌𝑏𝑢𝑠 = 𝑦2 + 𝑦12 + 𝑦23 −𝑦23

−𝑦13 −𝑦23 𝑦3 + 𝑦13 + 𝑦23

© Aktaibi
 Example: find 𝒀𝒃𝒖𝒔 for this system
© Aktaibi
 POWER-FLOW SOLUTION BY GAUSS–SEIDEL
• Nodal equations 𝐼 = 𝑌𝑏𝑢𝑠 𝑉 are a set of linear equations analogous to 𝑦 = 𝐴𝑥.
• Since power-flow bus data consists of 𝑃𝑘 and 𝑄𝑘 for load buses or 𝑃𝑘 and 𝑉𝑘 for voltage-
controlled buses, nodal equations do not directly fit the linear equation format; the current
source vector 𝐼 is unknown and the equations are actually nonlinear.
• The complex power 𝑆𝑘 delivered to bus 𝑘 is :
𝑆𝑘 = 𝑃𝑘 + 𝑗𝑄𝑘 = 𝑉𝑘 𝐼𝑘∗
• For each load bus, 𝐼𝑘 can be calculated from
𝑆𝑘 𝑃𝑘 + 𝑗𝑄𝑘
𝐼𝑘 = ∗ =
𝑉𝑘 𝑉𝑘∗
• Applying the Gauss–Seidel method to the nodal equations, with 𝐼𝑘 given above, we obtain

∗(𝑖)
• Equation (6.5.2) can be applied twice during each iteration for load buses, first using 𝑉𝑘 ,
∗(𝑖) ∗(𝑖+1)
then replacing 𝑉𝑘 , by 𝑉𝑘 on the right side of (6.5.2).
© Aktaibi
 Power-flow Solution By Gauss–Seidel (FYI)
• For a voltage-controlled bus, Qk is unknown, but can be calculated from
𝑁

𝑄𝑘 = 𝑉𝑘 𝑖 ෍ 𝑌𝑘𝑛 𝑉𝑛 𝑖 sin 𝛿𝑘 𝑖 − 𝛿𝑛 𝑖 − 𝜃𝑘𝑛 𝑎𝑛𝑑 𝑄𝐺𝑘 = 𝑄𝑘 + 𝑄𝐿𝑘

𝑛=1

• If the calculated value of QGk does not exceed its limits, then Qk is used in (6.5.2) to calculate Vk(i + 1) =
Vk(i + 1)/δk(i+1). Then the magnitude Vk(i + 1) is changed to Vk, which is input data for the voltage-controlled
bus. Thus we use (6.5.2) to compute only the angle δk(i + 1) for voltage-controlled buses.
• If the calculated value exceeds its limits QGkmax or QGkmin during any iteration, then the bus type is changed
from a voltage-controlled bus to a load bus, which QGK set to its limit value. Under this condition, the voltage-
controlled device (capacitor bank, Static Var system and so on) is not capable of maintaining Vk as specified by
the input data. The power flow program then calculates a new value of Vk.
• For the swing bus, denoted bus 1, V1 and δ1 are input data. As such, no iterations are required for bus 1.
after the iteration process has converged, one pass through (6.4.10) and (6.4.11) can be made to compute
P1 and Q1.
© Aktaibi
 Example 4:
• Assume a 1 + j0.5 p.u load is connected to a generator through a line with z = 0.02 + j0.06 pu.
Also, there is a -j0.25 capacitance at bus 2. If the generator voltage is 1.0 pu, what is V2?

𝑷 + 𝒋𝑸
© Aktaibi
 Example 4: cont.
© Aktaibi
 A Matlab code to solve example 4
clc 5.0000 -11.0000i
clear all -5.0000 +15.0000i
Y22 = 5-11i; -1.0000
Y21 = -5+15i; 0 - 0.5000i
Y11 = 5-15i;
Y12 = -5+15i; ---------------------------------------------------------------
V1 = 1; k V / Theta S1
P2 = -1; ---------------------------------------------------------------
Q2 = -0.5i; 1.0000 1.2449 -9.0218 0
disp([ Y22; Y21; P2; Q2])
disp(' --------------------------------------------------------------') 2.0000 1.2492 -8.0123 1.7808 - 4.4178i
disp(' k V / Theta S1')
disp(' --------------------------------------------------------------') 3.0000 1.2502 -8.0518 1.4269 - 4.4253i

for k=1:5 4.0000 1.2502 -8.0483 1.4374 - 4.4435i

V2(1)=1;
V2(k+1)=(1/Y22)*(((P2-Q2)/conj(V2(k)))-Y21*V1); 5.0000 1.2502 -8.0484 1.4361 - 4.4435i
X(k+1) = real(V2(k+1));
Y(k+1) = imag(V2(k+1));
[THETA,V] = cart2pol(X(k+1),Y(k+1)); S1 = 1.4361 - 4.4435i
Theta=THETA*180/pi; P1 = 1.4361
S1(k)=V1*(Y11*V1+Y12*V2(k)); Q1 = - 4.4435i
disp([ k V Theta conj(S1(k))])

end
© Aktaibi
 Example 5:
• Assume a 0.8 + j0.4 p.u load at bus 2 is being supplied by a generator at bus 1 through a transmission line with
series impedance of 0.05 + j0.1 p.u. Assuming bus 1 is the swing bus with a fixed per unit voltage of 1/0, use
the Gauss-Seidel method to calculate the voltage at bus 2 after three iterations.
© Aktaibi
 A Matlab code to solve example 5
clc 4.0000 - 8.0000i
clear all
-4.0000 + 8.0000i
Y22=4-8i;
-0.8000
Y21=-4+8i;
V1=1; 0 - 0.4000i
P2=-0.8;
Q2=-0.4i; ----------------------------------
disp([ Y22; Y21; P2; Q2]) k V2 / Theta
disp(' ----------------------------------') ----------------------------------
disp(' k V / Theta') 1.0000 0.9220 -3.7314
disp(' ----------------------------------')
2.0000 0.9111 -3.7314
for k=1:6
V2(1)=1; 3.0000 0.9101 -3.7802
V2(k+1)=(1/Y22)*(((P2-Q2)/conj(V2(k)))-Y21*V1);
X(k+1) = real(V2(k+1)); 4.0000 0.9099 -3.7802
Y(k+1) = imag(V2(k+1));
[THETA,V2] = cart2pol(X(k+1),Y(k+1));
5.0000 0.9099 -3.7809
Theta=THETA*180/pi;
disp([ k V2 Theta])
end 6.0000 0.9099 -3.7809
© Aktaibi
 Example 6:
• For the three-bus system whose Ybus is given below, calculate the second iteration value of V3 using the Gauss-
Seidel method. Assume bus 1 as the slack (with V1 = 1.0 /0), and buses 2 and 3 are load buses with a per unit
load of S2 = 1 + j0.5 and S3 = 1.5 + j0.75. Use voltage guesses of 1.0 /0 at both buses 2 and 3. The bus
admittance matrix for a three-bus system is
© Aktaibi
 Example 6: cont.
© Aktaibi
 A Matlab code to solve example 6
clc 0 -10.0000i
clear all 0 + 5.0000i
0 + 5.0000i
Y22 = -10i; -1.0000
Y21 = 5i; 0 - 0.5000i
Y33 = Y22;
Y23 = Y21;
Y31 = Y21;
0 -10.0000i
Y32 = Y21; 0 + 5.0000i
0 + 5.0000i
V1 = 1; -1.5000
P2 = -1; 0 - 0.7500i
Q2 = -0.5i; -----------------------------------------------------
k V2 / Delta2 V3 / Delta3
P3 = -1.5; -----------------------------------------------------
Q3 = -0.75i;
1.0000 0.9552 -6.0090 0.9220 -12.5288
disp([ Y22 ; Y21 ; Y23 ; P2 ; Q2])
2.0000 0.9090 -12.6225 0.8630 -16.1813
disp([ Y33 ; Y31 ; Y32 ; P3 ; Q3]) 3.0000 0.8640 -14.4489 0.8254 -17.6683
disp(' -----------------------------------------------------') 4.0000 0.8385 -15.4165 0.8046 -18.6306
disp(' k V2 / Delta2 V3 / Delta3') 5.0000 0.8241 -16.0356 0.7922 -19.2126
disp(' -----------------------------------------------------') 6.0000 0.8154 -16.3999 0.7847 -19.5626
7.0000 0.8101 -16.6203 0.7801 -19.7789
for k=1:24 8.0000 0.8068 -16.7566 0.7772 -19.9132
V2(1)=1; 9.0000 0.8048 -16.8410 0.7754 -19.9967
V3(1)=1;
10.0000 0.8036 -16.8935 0.7743 -20.0489
V2(k+1)=(1/Y22)*(((P2-Q2)/conj(V2(k)))-Y21*V1-Y23*V3(k));
V3(k+1)=(1/Y33)*(((P3-Q3)/conj(V3(k)))-Y31*V1-Y32*V2(k+1));
11.0000 0.8028 -16.9263 0.7736 -20.0815
12.0000 0.8023 -16.9468 0.7732 -20.1019
X2(k+1) = real(V2(k+1)); 13.0000 0.8020 -16.9597 0.7729 -20.1147
Y2(k+1) = imag(V2(k+1)); 14.0000 0.8018 -16.9677 0.7728 -20.1227
[THETA2,V_2] = cart2pol(X2(k+1),Y2(k+1)); 15.0000 0.8017 -16.9727 0.7727 -20.1277
Delta2=THETA2*180/pi; 16.0000 0.8016 -16.9759 0.7726 -20.1309
17.0000 0.8015 -16.9779 0.7725 -20.1329
X3(k+1) = real(V3(k+1)); 18.0000 0.8015 -16.9791 0.7725 -20.1341
Y3(k+1) = imag(V3(k+1));
19.0000 0.8015 -16.9799 0.7725 -20.1349
[THETA3,V_3] = cart2pol(X3(k+1),Y3(k+1));
Delta3=THETA3*180/pi;
20.0000 0.8015 -16.9804 0.7725 -20.1354
21.0000 0.8015 -16.9807 0.7725 -20.1357
disp([ k V_2 Delta2 V_3 Delta3]) 22.0000 0.8015 -16.9809 0.7725 -20.1358
end 23.0000 0.8015 -16.9810 0.7725 -20.1360
24.0000 0.8015 -16.9810 0.7725 -20.1360
© Aktaibi