16-Jul-20

EEE 471 Power System Analysis-I

Chapter 6: Power Flow Analysis

EEE 471
POWER SYSTEM ANALYSIS - I
Assoc. Prof. Dr. A. Mete VURAL
E-mail: [email protected]

Chapter 6: Power Flow Analysis

1
 16-Jul-20

Introduction to Power Flow Analysis

 In a three phase ac power system, active and reactive power flows from the
generating stations to the loads through transmission line/network.
 The flow of active and reactive powers on a transmission line of a power
system/network is called power flow or load flow.
 Power flow analysis is very important in planning stages of power network to
add new transmission lines and new generation sides.
 The results of power flow analysis is widely used by electrical power
engineers during the planning stages and operation of power systems.
 In power flow studies the power system is represented by one-line (single-
line) diagram.
 In power studies, the analyzed is assumed to be balanced and operating
under steady-state and normal conditions.

Objective of Power Flow Study

 To find bus voltage magnitude and bus voltage phase angle of each bus
in a given power system.
 To find real and reactive power flows on each line of a given power
system.
 To find real and reactive power losses of each line of a given power
system.
 The solution of the power flow study is used to determine
under/overloaded transmission lines so that suitable line compensation
methods can be designed and applied.
 For example to add a FACTS device in a power system, before studying
power flow analysis, it is impossible to define the type, parameters, and
the best location of the FACTS device which is considered for
compensation.

2
 16-Jul-20

Properties of Power Flow Study

 The bus admittance matrix of the power system is required for power flow
study giving impedance/admittance information of each
line/generator/transformer.
 Non-linear power flow equations for each bus should be described,
linearized, and solved by iterative methods.
 Commonly used iterative methods are Gauss-Seidel, Newton-Raphson,
and fast decoupled methods.
 For a power system having many buses and branches (lines) the
equations are generally put in matrix form for calculation efficiency.
 For example in a 30-bus power system, there are 2*30=60 unknows to be
solved by power flow study.
 Hand-made calculations are very hard to implement so computer software
is used for power flow (load flow) analysis.

Some Examples to Power Flow Analysis Software

 PSS/E (commercial)
 DigSILENT (commercial)
 EasyPower (commercial)
 Eurostag (commercial)
 ETAP (commercial)
 PSASP (power system analysis software package)
 Powerworld (commercial, demo is available)
 Matpower (non-commercial, matlab based free program)
 PSAT (non-commercial matlab based free program)

3
 16-Jul-20

Bus Admittance Matrix

Generator

Generator impedance

Bus (node)

Line reactance
PU admittance
between
bus i & j
PU impedance

reactance
resistance
Bus-0 is ground node taken as between
Reference (not shown) bus i & j

Impedance diagram of a simple power system

Transformations Made to Obtain Admittance Diagram:

 All line impedances are converted into line admittances
 All voltage source+impedance are converted into a current source+parallel admittance

admittance

Admittance diagram of a simple power system

4
 16-Jul-20

Applying KCL to Buses 1-4:

Rearranging the above equations:

We introduce the following admittances:

The node equations reduces to:

Y14=Y41=0 and Y24=Y42=0

Since there is no connection between Bus 1 and Bus 4
and there is no connection between Bus 2 and Bus 4

10

5
 16-Jul-20

For a general n-bus power system the following matrix equation can be written:

or

Injected Bus voltage

Bus admittance matrix
current vector vector
Positive if injected into bus
Negative if outgoing from bus

mutual or transfer admittance self-admittance

or driving-point admittance
11

Since

If injected current vector is known

voltage bus vector can be calculated

The inverse of bus admittance matrix is called Bus Impedance Matrix (Zbus)

Notes on Bus Admittance Matrix (Ybus):

 Zbus can be calculated only if Ybus is inversible (non-singular)

 Ybus is a symmetrical matrix since Yij=Yji (there is only one connection configuration between buses i & j)
 In real, many off-diagonal elements are zero, since each bus is connected to only a few nearby buses in a power
system
 Ybus is a sparse matrix because of many zeros and there are efficient numerical methods to find its inverse
rather than calculating its inverse directly. Since for example for a 300-bus system to find Zbus we need to find the
inverse of a 300x300 matrix, which is not so efficient and slow even with a fast computer.
 In real, Zbus which is required for short-circuit analysis can be obtained more efficiently and fast by applying
building algorithm without the need for matrix inversion.

12

6
 16-Jul-20

Example: Find the bus admittance matrix of the following simple power system with the given
admittance diagram

Impedance diagram of the power system

13

Example:

E1 E2

14

7
 16-Jul-20

The code of Matlab function «ybus1.m»

% Bus Admittance Matrix
% Copyright (c) 1998-2010 by H. Saadat.

function[Ybus] = ybus1(zdata)
nl=zdata(:,1); nr=zdata(:,2); R=zdata(:,3); X=zdata(:,4);
nbr=length(zdata(:,1)); nbus = max(max(nl), max(nr));
Z = R + j*X; %branch impedance
y= ones(nbr,1)./Z; %branch admittance
Ybus=zeros(nbus,nbus); % initialize Ybus to zero
for k = 1:nbr; % formation of the off diagonal elements
if nl(k) > 0 && nr(k) > 0
Ybus(nl(k),nr(k)) = Ybus(nl(k),nr(k)) - y(k);
Ybus(nr(k),nl(k)) = Ybus(nl(k),nr(k));
end
end
for n = 1:nbus % formation of the diagonal elements
for k = 1:nbr
if nl(k) == n || nr(k) == n
Ybus(n,n) = Ybus(n,n) + y(k);
else, end
end
end 15

Let’s execute the following commands in Matlab’s workplace

Put a semicolon ; after writing each X

16

8
 16-Jul-20

17

18

9
 16-Jul-20

19

20

10
 16-Jul-20

21

Notes on the Previous Solution Technique:

 The solution of the equation Vbus = Zbus.Ibus by matrix inversion is very inefficient.
 Actually, it is not necessary to obtain the inverse of Ybus. Instead, direct solution is obtained by optimally
ordered triangular factorization.
 In Matlab, the solution of the linear equation AX=B can be obtained by using the matrix division operator
(X=A\B), which is based on the triangular factorization and Gaussian elimination.
 This technique is superior in both execution time (two to three times faster) and numerical accuracy.

We can obtain the direct solution of the previous example by executing the following command in
Matlab.

Vbus = Y \ Ibus

22

11
 16-Jul-20

Example: Find the voltages of the buses of the following power system with the given impedance diagram
using Matlab function «ybus1.m»

23

Solution:

First we write the following Matlab command to enter power system data:

>> z= [0 1 0 0.8;0 2 0 1.0;0 3 0 2.0;1 2 0 0.5; 1 3 0 0.2;2 3 0 0.4];

z=

0 1.0000 0 0.8000
0 2.0000 0 1.0000
0 3.0000 0 2.0000
1.0000 2.0000 0 0.5000
1.0000 3.0000 0 0.2000
2.0000 3.0000 0 0.4000

24

12
 16-Jul-20

>> Y = ybus1(z)

Y=

0.0000 - 8.2500i 0.0000 + 2.0000i 0.0000 + 5.0000i Bus Admittance Matrix

0.0000 + 2.0000i 0.0000 - 5.5000i 0.0000 + 2.5000i of the power system
0.0000 + 5.0000i 0.0000 + 2.5000i 0.0000 - 8.0000i
No injected current at Bus 3
(no Gen. is connected)
>> 1.28/(j*0.8)

ans = >> Ibus=[-j*1.6; 0.6-j*1.0392; 0]

0.0000 - 1.6000i Ibus =

>> (1.2*cos(pi/6)+j*1.2*sin(pi/6))/j 0.0000 - 1.6000i

0.6000 - 1.0392i
ans = 0.0000 + 0.0000i

0.6000 - 1.0392i 25

>> Vbus=inv(Y)*Ibus

Vbus =

0.9791 + 0.1860i
0.9594 + 0.2676i
0.9118 + 0.1999i

or

>> Vbus = Y \ Ibus

Vbus =

0.9791 + 0.1860i
0.9594 + 0.2676i
0.9118 + 0.1999i
26

13
 16-Jul-20

How Ybus matrix is effected if

• Adding transformer at a bus ?

• Adding new lines ?
• Removing existing lines ?
• Adding shunt capacitor at a bus ?
• Adding shunt reactor at a bus ?

27

Let’s extend the previous example to find other power system parameters:

Q12 Real power flow from Bus 1 to Bus 2

P12 P12=((V1*V2)/X12)*sin(theta1-theta2)

Reactive power flow from Bus 1 to Bus 2

Q12=(V12 - V1*V2cos(theta1-theta2))/X12

28

14
 16-Jul-20

>> Vbus1= 0.9791 + 0.1860i

Vbus1 =

0.9791 + 0.1860i

>> Vbus2= 0.9594 + 0.2676i

Vbus2 =

0.9594 + 0.2676i

>> Vbus1_mag=sqrt(0.9791^2+0.1860^2)

Vbus1_mag =

0.9966

>> Vbus1_ang=atan2(0.1860,0.9791)

Vbus1_ang =
29
0.1877

>> Vbus2_mag=sqrt(0.9594^2+0.2676^2)

Vbus2_mag =

0.9960

>> Vbus2_ang=atan2(0.2676,0.9594)

Vbus2_ang =

0.2720

>> P12=((Vbus1_mag*Vbus2_mag)/(0.5))*sin(Vbus1_ang-Vbus2_ang)

P12 =
-0.1671

>> Q12=(Vbus1_mag*Vbus2_mag*cos(Vbus1_ang-Vbus2_ang)-Vbus2_mag^2)/0.5

Q12 =
-0.0059

30

15
 16-Jul-20

Solution of Non-Linear Algebraic Equations

 The previous method is very practical if all voltage magnitudes of generators are known
 Actually, only one bus voltage magnitude and phase angle is known to solve power flow
problem (slack bus)
 In the previous method, the load information is not available, however should be
considered in an actual power flow problem.
 Actually a power flow problem is a non-linear multi-variable problem that should be
solved by computer software iteratively.
 The followings are the common methods used for iterative solution of non-linear power
flow equations:

 Gauss-Seidel Method
 Newton-Raphson Method (powerful and popular)
 Quasi-Newton Method

31

Gauss-Seidel Method:

32

16
 16-Jul-20

Solution:

initial estimate of x

33

4.5

3.5
3 2 
g(x) =-1/9x +6/9x +4/9 x
3

2.5

The root of the polynomial 1.5

after 9 iterations 1

0.5

0
The Limitations of Gauss-Seidel Method: 0 1 2 3 4 5
x

 This method needs many iterations to achieve the desired accuracy.

 There is no guarantee for the convergence.
 The solution is highly dependent on initial estimate of x. For example if initial estimate of x were
chosen as 6, the process would diverge.
 For a power system having many buses, convergence test is difficult to apply and there is no general
method for convergence test.

34

17
 16-Jul-20

Newton-Raphson Method:

35

36

18
 16-Jul-20

Equations of NR algorithm
k: iteration number

37

Solution:

38

19
 16-Jul-20

The root of the polynomial after 6 iterations

50

40

30
f(x) = x -6x +9x-4
3 2

20

10

-10
0 1 2 3 4 5 6
39
x

Type the following command in Matlab Workspace to examine the same example with different
initial estimates
>> chp6ex4
>> chp6ex4 Enter the initial estimate -> 21
Enter the initial estimate -> 9 iter Dc J dx x
iter Dc J dx x 1 1.0e+03 *
1 -320.0000 144.0000 -2.2222 6.7778
-6.8000 1.0800 -0.0063 0.0147
2 -92.7298 65.4815 -1.4161 5.3617
2 1.0e+03 *
3 -25.9042 30.9022 -0.8383 4.5234
-2.0101 0.4812 -0.0042 0.0105
4 -6.4975 16.1025 -0.4035 4.1199
3 -592.2208 215.0830 -2.7535 7.7726
5 -1.1669 10.4817 -0.1113 4.0086
4 -173.0465 96.9703 -1.7845 5.9881
6 -0.0774 9.1029 -0.0085 4.0000
5 -49.4669 44.7152 -1.1063 4.8818
7 -0.0004 9.0006 -0.0000 4.0000
6 -13.2884 21.9152 -0.6064 4.2755

7 -2.9557 12.5336 -0.2358 4.0397

8 -0.3665 9.4808 -0.0387 4.0010

9 -0.0091 9.0121 -0.0010 4.0000

10 -0.0000 9.0000 -0.0000 4.0000 40

20
 16-Jul-20

The Comments on NR Solution:

 Newton-Raphson (NR) method converges considerably more rapidly than the Gauss-Seidel method.
 Although the initial estimate of x=6 makes Gauss-Seidel method diverges, the same initial estimate
of x makes NR method converges.
 NR method is generally more tolerant when choosing the initial estimate of x, even initial estimate
is not close enough to the root of the equation.
 The limitation of NR method is that it requires the first order derivative of the function whose root is
under consideration.

41

N-dimensional Newton-Raphson Method

42

21
 16-Jul-20

Jacobian Matrix consisting of partial derivatives

43

The Comments on N-Dimensional NR Algorithm:

 Elements of Jacobian matrix are the partial derivatives evaluated at (k)

 It is assumed that Jacobian matrix has an inverse during each iteration.
 N-dimensional NR algorithm reduces the non-linear equations to a N-dimensional linear equations
to be solved by iteration.
 Actually the inversion of Jacobian matrix is inefficient and not necessary. Instead, a direct solution is
obtained by optimally ordered triangular factorization.
 In Matlab, the solution of linear equations ΔC=JΔX is obtained by using the matrix division operator
which is based on the triangular factorization and Gaussian elimination.

For example,

ΔX = J \ ΔC

44

22
 16-Jul-20

Solution:
3

(Jacobian matrix ) 2 2 2
x +y =4

-1
solution1
-2
x solution2
e +y=1

-3
-3 -2 -1 0 1 2 3
x

45

>> chp6ex5 >> chp6ex5

Enter initial estimates, col. vector [x1; x2] -> [1;1] Enter initial estimates, col. vector [x1; x2] -> [0.5;-1]
Iter DC Jacobian matrix Dx x Iter DC Jacobian matrix Dx x
1 2.0000 2.0000 2.0000 -2.1640 -1.1640 1 2.7500 1.0000 -2.0000 0.8034 1.3034
-2.7183 2.7183 1.0000 3.1640 4.1640 0.3513 1.6487 1.0000 -0.9733 -1.9733

2 -14.6933 -2.3279 8.3279 -2.8927 -4.0567 2 -1.5928 2.6068 -3.9466 -0.2561 1.0473
-3.4762 0.3122 1.0000 -2.5730 1.5910 -0.7085 3.6818 1.0000 0.2344 -1.7389

3 -14.9879 -8.1134 3.1820 1.5979 -2.4588 3 -0.1205 2.0946 -3.4778 -0.0422 1.0051
-0.6083 0.0173 1.0000 -0.6360 0.9550 -0.1111 2.8499 1.0000 0.0092 -1.7296

4 -2.9577 -4.9176 1.9101 0.5669 -1.8919 4 -0.0019 2.0102 -3.4593 -0.0009 1.0042
-0.0406 0.0855 1.0000 -0.0891 0.8660 -0.0025 2.7321 1.0000 0.0000 -1.7296

5 -0.3293 -3.7838 1.7319 0.0742 -1.8177 5 -0.0000 2.0083 -3.4593 -0.0000 1.0042
-0.0168 0.1508 1.0000 -0.0279 0.8380 -0.0000 2.7296 1.0000 -0.0000 -1.7296

6 -0.0063 -3.6354 1.6761 0.0014 -1.8163

-0.0004 0.1624 1.0000 -0.0007 0.8374

7 -0.0000 -3.6325 1.6747 0.0000 -1.8163

-0.0000 0.1626 1.0000 -0.0000 0.8374 46

23
 16-Jul-20

iter =

DC =

1.0e-04 *

-0.3669
-0.1784
0.0957

Solution: J=

4.0000 -6.0000 8.0000 only the

3.0000 8.0000 -3.0000 results of
-3.0000 4.0000 1.0000
6. iteration
(Jacobian matrix )
are shown
Dx =

1.0e-05 *

-0.5577
-0.1130
-0.2645

x=

2.0000
3.0000
4.0000 47

Power Flow Solution

In a power flow study the power system buses are classified into three-types:

or known as PQ Bus
P,Q are known
V and angle of the bus
are calculated

or known as PV Bus
P, V are known
Q and angle of the bus
are calculated

V, angle of the bus are

known
P, Q are calculated

48

24
 16-Jul-20

Power Flow Equation

 Transmission lines are represented by their equivalent pi models
 All the quantities are converted into pu values

Appliying KCL at Bus i

or
A typical bus (Bus-i) of a power system

(Total sum of entering currents = Total sum of outgoing currents)

49

Power Flow Equation

This equation is the non-linear power flow equation for Bus i

that should be solved using numerical techniques

50

25
 16-Jul-20

Gauss-Seidel Power Flow Solution

 yij is the pu admittance of the branch (or line) ij

 Pisch is the net injected real power at Bus i expressed in pu
 Qisch is the net injected reactive power at Bus i expressed in pu
 For a generator Bus, Pisch and Qisch are positive
 For a load Bus, Pisch and Qisch are negative
 k is the iteration number

51

Gauss-Seidel Power Flow Solution

 Power flow equation is usually expressed in terms of the elements of bus admittance matrix Ybus.
 Off-diagonal elements of Ybus are Yij = -yij
 Diagonal elements of Ybus are Yii = sum(yij)

53

26
 16-Jul-20

Gauss-Seidel Power Flow Solution

 Since voltage magnitude and phase angle of slack bus are known, there are 2(n-1) equations which must be
solved iteratively. For example if bus number is 4, there are 2(4-1) = 6 equations to be solved iteratively.

 Under normal operating conditions, the voltage magnitude of the buses are in the neighbourhood of 1.0 pu.

 Voltage magnitudes of the load buses are generally lower than the slack bus value.

 Voltage magnitudes of the generator buses are generally higher than the slack bus value.

 Phase angle of the load buses are generally below the reference angle of the slack bus.

 Phase angle of the generator buses are generally higher than the reference angle of the slack bus.

 Thus for Gauss-Seidel method 1.0 + j0.0 is a satisfactory initial voltage estimate for the unknown bus voltage

54

Gauss-Seidel Power Flow Solution

Pisch and Qisch are
For PQ Buses Use:
known

Pisch and magnitude of Use first:

For PV Buses bus voltage (|Vi|) are
known
second

Since bus voltage magnitude is known, to find the phase angle the following equation is used:

55

27
 16-Jul-20

Gauss-Seidel Power Flow Solution

The voltage magnitude and the

For Slack Bus phase angle of the slack bus are Use:
known

56

Line Flows and Line Losses

After finding all bus voltages with the Gauss-Seidel method, line flows and line losses can be computed:

(Complex power flow from Bus i to j )

(Complex power flow from Bus j to i )

(The power loss in the line ij)

58

28
 16-Jul-20

59

60

29
 16-Jul-20

61

62

30
 16-Jul-20

63

64

31
 16-Jul-20

65

66

32
 16-Jul-20

67

End of Chapter 6

73

33