Formation of Bus Admittance Matrix using MATLAB Software

S.NO.5 Formation of Bus Admittance Matrix using MATLAB Software Aim: To develop a computer program to form the bus admittance matrix, Ybus of a power system. Theory: The Ybus /Zbus matrix constitutes the models of the passive portions of the power network. Ybus matrix is often used in solving load flow problems. It has gained widespread applications owing to its simplicity of data preparation and the ease with which the bus admittance matrix can be formed and modified for network changes. Of course, sparsity is one of its greatest advantages as it heavily reduces computer memory and time requirements. In short circuit analysis, the generator and transformer impedances must also be taken into account. In contingency analysis, the shunt elements are neglected, while forming the Z-bus matrix, which is used to compute the outage distribution factors. This can be easily obtained by inverting the Y-bus matrix formed by inspection method or by analytical method. The impedance matrix is a full matrix and is most useful for short circuit studies. Initially, the Y-bus matrix is formed by inspection method by considering line data only. After forming the Y-bus matrix, the modified Y-bus matrix is formed by adding the generator and transformer admittances to the respective diagonal elements and is inverted to form the Zbus matrix. The performance equation for a n-bus system in terms of admittance matrix can be written as, Y .... Y   V  Y I  Y I  Y .... Y   V        .   . .  .        .  .  . .   Y  I  Y .... Y   V  (or) I = Ybus.V The admittances Y11, Y12,… Y1n are called the self-admittances at the nodes and all other admittances are called the mutual admittances of the nodes. 1

11

12

2

21

22

n

n 1

n 2

1

In

2 n

2

nn

n

Formulae Used: n

Main diagonal element in Y-bus matrix =

Y j 1

ij

 Bij

where Bij is the half line shunt admittance in mho. Yij is the series admittance in mho. Off-diagonal element in Y-bus matrix , Yij = -Yij

where Yij is the series admittance in mho. Flowchart: START

Enter the number of buses,n and lines Enter the mutual admittance between the buses

Set the bus count i =1

i = i +1

Calculate the diagonal term, Yii = sum of all admittances connected to bus i.

Is i = n

Calculate the off-diagonal term, Yij=Negative sum of the admittances connected from bus i to bus j. Compute Z bus matrix by inverting Y bus matrix Print Y bus and Z bus matrices

STOP

Algorithm: Step 1:

Read the values of number of buses and the number of lines of the given system.

Step 2:

Read the self-admittance of each bus and the mutual admittance between the buses.

Step 3:

Calculate the diagonal element term called the bus driving point admittance, Yij which is the sum of the admittance connected to bus i.

Step 4:

The off-diagonal term called the transfer admittance, Yij which is the negative of the admittance connected from bus i to bus j.

Step 5:

Check for the end of bus count and print the computed Y-bus matrix.

Step 6:

Compute the Z-bus matrix by inverting the Y-bus matrix.

Step 7:

Stop the program and print the results.

Sample Problem: The bus and branch datas for a 3 bus system is given in table below. Form Y bus matrix by inspection method. Bus Code 1-2 1–3 2-3

Impedance 0.06 + j0.18 0.02 + j0.06 0.04 + j0.12

Bus Number 1 2 3

Admittance j0.05 j0.06 j0.05

Solution: Formation of Y bus:

Ybus

  1 1   j0.05    0.06  j0.18 0.02  j 0.06  1   0.06  j0.18  1   0.02  j0.06 

1 0.06  j0.18   1 1    j 0.06  0.06  j0.18 0.04  j0.12  1 0.04  j0.12

    1  0.04  j0.12     1 1    j0.05  0.02  j0.06   0.04  j0.12 1 0.02  j0.06

Theoretical output:

Ybus

6.66  j19.95   1.66  j 5  5  j15

 1.66  j 5 4.16  j12.44  2.5  j 7.5

 5  j15

  2.5  j 7.5  7.5  j 22.45

Result: The Y bus matrix was formed for the given system by direct inspection method and the results were verified using MATLAB program.

S.NO.6 Z-bus Building Algorithm using MATLAB Software Aim: To develop a computer program to obtain the building algorithm for bus impedance matrix of the given power system. Theory: The Ybus /Zbus matrix constitutes the models of the passive portions of the power network. The impedance matrix is a full matrix and is most useful for short circuit studies. An algorithm for formulating [Zbus] is described in terms of modifying an existing bus impedance matrix designated as [Zbus]old. The modified matrix is designated as [Zbus]new. The network consists of a reference bus and a number of other buses. When a new element having self impedance Zb is added, a new bus may be created (if the new element is a tree branch) or a new bus may not be created (if the new element is a link). Each of these two cases can be subdivided into two cases so that Zb may be added in the following ways: 1. Adding Zb from a new bus to reference bus. 2. Adding Zb from a new bus to an existing bus. 3. Adding Zb from an existing bus to reference bus. 4. Adding Zb between two existing buses. Type 1 modification: In type 1 modification, an impedance Zb is added between a new bus p and the reference bus as shown in Figure 1 1 n Network

p Zb

Vp

Ref. Bus Figure 1. Type 1 modification of Zbus

Let the current through bus p be Ip, then the voltage across the bus p is given by, Vp = Ip Zb

The potential at other buses remains unaltered and the system equations can be written as, V1    0   I1        0   I 2  V2      Z bus  old  0            0       V    0  I n   n                   0 0 0 0 0  Z    b   I p  V p   Type 2 modification: In type 2 modification, an impedance Zb is added between a new bus p and an existing bus k as shown in Figure 2. The voltages across the bus k and p can be expressed as, Vk(new) = Vk + Ip Zkk Vp

= Vk(new) + Ip Zp = Vk + Ip(Zb + Zkk)

where, Vk is the voltage across bus k before the addition of impedance Zb Zkk is the sum of all impedance connected to bus k. 1 n Ik + Ip

Network

k Z

Ip

b

p Ref. Bus Figure 2.Type 2 Modification of Zbus

The system of equations can be expressed as, V1    Z 1k       Z 2k V2        Z bus  old             V       n               V   Z Z     Z  Z  kk b  p   k1 k 2

I1    I 2        I   n   I   p

Type 3 Modification: In this modification, an impedance Zb is added between a existing bus k and a reference bus. Then the following steps are to be followed: 1. Add Zb between a new bus p and the existing bus k and the modifications are done as in type 2. 2. Connect bus p to the reference bus by letting Vp = 0. To retain the symmetry of the Bus Impedance Matrix, network reduction technique can be used to remove the excess row or column. Type 4 Modification: In this type of modification, an impedance Zb is added between two existing buses j and k as shown in Figure 3. From Figure 3, the relation between the voltages of bus k and j can be written as, Vk – Vj = IbZb (3) 1 n Ij + Ib

Network

j Z

Ib

b

Ik - Ib

k Ref. Bus

Figure 3.Type 4 Modification of Zbus

The voltages across all the buses connected to the network changes due to the addition of impedance Zb and they can be expressed as, V1 = Z11I1 + Z12I2 + - - - - - - - - + Z1j(Ij + Ib) + Z1k(Ik – Ib)+- - V2 = Z21I1 + Z22I2 + - - - - - - - - + Z2j(Ij + Ib) + Z2k(Ik – Ib)+ - -     Vj = Zj1I1 + Zj2I2 + - - - - - - - - + Zjj(Ij + Ib) + Zjk(Ik – Ib) + - - (4) Vk = Zk1I1 + Zk2I2 + - - - - - - - - + Zkj(Ij + Ib) + Zkk(Ik – Ib) + - -     Vn = Zn1I1 + Zn2I2 + - - - - - - - - + Znj(Ij + Ib) + Znk(Ik – Ib) + - - -

On solving the Equations (3) and (4), the system of equations can be rewritten as,

 ( Z 1 j  Z 1k ) V1        V2      Z bus  old            V   ( Z kj  Z kk )  n                 V   Z bb  p  ( Z j1  Z k 1 )    ( Z jk  Z kk ) 

          

I1    I 2        I   n   I   p

(5)

where, Zbb = Zjj + Zkk – 2 Zjk + Zb

Procedure for formation of Zbus matrix: Step1: Number the nodes of the given network, starting with those nodes at the ends of branches connected to the reference node. Step2: Start with a network composed of all those branches connected to the reference node. Step3: Add a new node to the ith node of the existing network. Step4: Add a branch between ith and jth nodes. Continue until all the remaining branches are connected. Sample problem: Form bus impedance matrix using building algorithm:

Solution: Step1: Add an element between ref (0) bus and a new bus (1).

Z = [j0.2]

Step2: Add an element between existing bus (1) to a new bus (2).

Z=

 j 0.2  j 0.2 

j 0.2  j 0.6

Step3: Add an element between existing (2) Bus to a ref (0) Bus.

 j 0.2 j 0.2 j 0.2  Z=  j 0.2 j 0.6 j 0.6   j 0.2 j 0.6 j 0.8 New Z Bus: Z11 Z11 Z12 Z22 Z22

= Z11-(Z31*Z13)/Z33 = j0.2 – (j0.2*j0.2)/j0.8 = j0.05 =Z21= Z12-(Z32*Z13)/Z33 = j0.2 - (j0.6*j0.2)/j0.8 = j0.05 =Z22-(Z32*Z23)/Z33 =J0.6-(j0.6*j0.6)/j0.8 =j0.15 Z Bus

=

 j 0.05  j 0.05 

j 0.05  j 0.15

Result: The bus impedance matrix using building algorithm for the given system was formed and the results were verified using MATLAB program.

S.NO.7 Gauss Seidal Load flow analysis using MATLAB software Aim: To develop a computer program to solve the set of non linear load flow equations using Gaussseidal load flow algorithm. Theory: Load flow analysis is the most frequently performed system study by electric utilities. This analysis is performed on a symmetrical steady-state operating condition of a power system under ‘normal’ mode of operation and aims at obtaining bus voltages and line/transformer flows for a given load condition. This information is essential both for long term planning and next day operational planning. In long term planning, load flow analysis helps in investigating the effectiveness of alternative plans and choosing the ‘best’ plan for system expansion to meet the projected operating state. In operational planning, it helps in choosing the ‘best’ unit commitment plan and generation schedules to run the system efficiently for them next day’s load condition without violating the bus voltage and line flow operating limits. The Gauss seidal method is an iterative algorithm for solving a set of non- linear algebraic equations. The relationship between network bus voltages and currents may be represented by either loop equations or node equations. Node equations are normally preferred because the number of independent node equation is smaller than the number of independent loop equations. The network equations in terms of the bus admittance matrix can be written as, I bus  Ybus Vbus

(1)

For a n bus system, the above performance equation can be expanded as,  I 1  Y11 Y12  Y1 p  Y1n  V1     Y Y  Y2 p  Y2 n  V2  22  I 2   12                             I  Y p1 Y p 2  Y pp Y pn  V p  p                     Ynn  Vn   I n  Yn1 Yn 2  Ynp where n is the total number of nodes. Vp is the phasor voltage to ground at node p. Ip is the phasor current flowing into the network at node p.

(2)

At the pth bus, current injection: I p  Y p1V1  Y p 2V2  ............  Y ppV p  .............Y pnVn 

n

 Y pqVq

q 1

 Y ppV p 

n

 Y pqVq

(3)

q 1 q p

  n  1  Vp   I p   Y pqVq  ; p  2,....n Y pp   q 1   q p At bus p , we can write Pp – jQp = V p I p Hence, the current at any node p is related to P, Q and V as follows:

(4)

( Pp  jQ p )

( for any bus p except slack bus s) (5) V p Substituting for Ip in Equation (4),   n   P  jQ 1 p p Vp    Y pqVq  ; p  2, .....n (6)  Y pp  V p*  q 1   q p Ip has been substituted by the real and reactive powers because normally in a power system these quantities are specified.

Ip 

Algorithm: Step 1: Step 2: Step 3: Step 4: Step 5: next

Read the input data. Find out the admittance matrix. Choose the flat voltage profile 1+j0 to all buses except slack bus. Set the iteration count p = 0 and bus count i = 1. Check the slack bus, if it is the generator bus then go to the next step otherwise go to

step 7. Step 6: Before the check for the slack bus if it is slack bus then go to step 11 otherwise go to next step. Step 7: Check the reactive power of the generator bus within the given limit. Step 8: If the reactive power violates a limit then treat the bus as load bus. Step 9: Calculate the phase of the bus voltage on load bus

Step 10: Calculate the change in bus voltage of the repeat step mentioned above until all the bus voltages are calculated. Step 11: Stop the program and print the results

Flowchart:

Start

Read the input data values

Form Y Bus matrix

Set flat voltage profile 1+j0 except slack bus

Set iteration count, p=0 A

Set the bus count, i = 1 D

Yes Y

Check for slack bus

E

No

It is a load bus calculate

No

Check for Gen bus

p 1 ical

V

1  Yii

n  Pi  jQi j 1   Y V  Y V     ik k ik k * k 1 j 1  Vi 

Yes Calculate



* i 1

Q   Im Vip Y V  Y V p1

i

n

p1

k 1

ik

k

k i

ik

k

p

 C

Yes Check

Qi

p 1

Set Qi=Qi min

Qmin No

Yes Check

Qi

p 1

 Qmax No B

Set Qi=Qi max

B

Treat this as gen bus & calculate Vpi i 1 n  1  Pi  jQ p 1 Vi p 1   Y V  Yik Vkp     ik k * Yii  Vi k 1 k  i 1  C

Calculate the change in voltage

Vi p 1

E

Increment the bus count

No

D

Check i n A

Yes No

Check Vi p 1 

Increment iteration count P = P+1

Yes Print the result

Stop

Sample Problem: The load flow data for a 3 bus system is given in tables below. The voltage magnitude at bus 2 is to be maintained at 1.04 p.u. The maximum and minimum reactive power limits for bus 2 are 0.5 to 0.2 respectively. Taking bus 1 as slack bus, determine voltages of the various buses at the end of first iteration starting with flat voltage profile for all buses except slack bus using Gauss-Seidal method with acceleration factor of 1.6.

Bus Code 1–2 1–3 2–3 Bus Code 1 2 3

Impedance 0.06 + j0.18 0.02 + j0.06 0.04 + j0.12 Assumed Voltage 1.06 + j0 1 + j0 1 + j0

Bus Number 1 2 3 Generation MW MVAr 0 0 0.2 0 0 0

Admittance j0.05 j0.06 j0.05 Load MW MVAr 0 0 0 0 0.6 0.25

Solution: Formation of Ybus:

Ybus

6.66  j19.95   1.66  j 5  5  j15

 1.66  j 5 4.16  j12.44  2.5  j 7.5

 5  j15

  2.5  j 7.5  7.5  j 22.45

Calculation of Q2:  * n  Q2 =  Im V2  YpqVq  q 1   =  Im1.04(1.66  j 5)(1.06)  (4.16  j12.5)1.04  (2.5  j 7.5) =  Im1.04(1.763  j 5.30)  (4.16  j12.5)1.04  (2.5  j 7.5) =  Im0.07  j 0.14 Q2 = 0.14, it violates the limits of the reactive power. Q2 = Q min = 0.2 as Q2  Qmin [ If suppose, Q2  Qmax then Q2 = Qmax] Calculation of Bus voltages:  0 .2  j 0 .2  (1) V2 = 0.075 71.63  ((1.66  j 5)(1.06)  (2.5  j 7.5)(10)) Volts  1.04  = 0.075   71.634.452  j12.99 (1)

V2 = 1.047+j 0.007 volts

Accelerated voltage, (1) V2 = 1.04+ 1.6(1.047+j0.007 -1.04) = 1+0.048-j0.048 (1) V2 =1.0512+j0.0112 Volts

V3 =0.0423 71.49 0.6  j 0.25  ((5  j15)(1.06)  (2.5  j 7.5)(1.0512  j 0.0112)) (1)

V3

(1)

= 1.041 – j 0.017 Volts

Accelerated voltage, (1)

V3 = 1+1.6(1.041 – j 0.17- 1 ) V 3(1) = 1.0656-j0.272 Volts Theoretical Output: (1)

V1=1.06+j0 Volts, V2 =1.0512+j0.0112 Volts, V 3(1) = 1.0656-j0.272 Volts

Result: The given set of load flow equations for a given power system were solved using Gauss-Seidal method.

S.NO.9 Fast Decoupled Load Flow Analysis using MATLAB Software Aim: To become proficient in the usage of software in solving load flow problems using Fast decoupled load flow method. Theory: Load flow study is useful in planning the expansion of power system as well as determining best operation of the system. The principle obtained from load flow study is the magnitude and phase angle of the voltage at each bus and real and reactive power flowing in each line. Load flow analysis may be performed using A.C. network analyzer and also by digital computer. But nowa-days digital computer oriented load flow analysis is a standard practice. The fast decoupled load flow method is a very fast method of obtaining load flow solutions. This method requires less number of arithmetic operations to complete an iteration consequently. This method requires less time per iterations. In N-R method, the elements of Jacobian are to be computed in each iteration .So the time per iteration is considerably more in N-R method than in FDLF. The rate of convergence in FDLF method is slow requiring considerably more number of iterations to obtain a solution than in the case of N-R method. However accuracy is same in both the cases. In this method both the speeds as well as the sparsity are exploited. This is an extension of N-R method formulated in polar co-ordinates with certain approximation which results into a fast algorithm for load flow solution. In practice, transmission system operating under steady state possesses strong interdependence between active powers and bus voltages, angles, similarly there is strong interdependence between bus voltage and reactive power

H km 

Pk  m

;

N km 

J km 

Qk  m

;

Lkm 

H pq 

Pp

Pk E m E m Qk E m

L pq 

;

E m QP E q

 q E q The equation for power flow are again expressed below for calculating elements of Jacobian (ie H & L)

Pp  E p E qY pp cos  pp 

n

E

q 1,  p

p

E q Y pq cos pq   p   q 

Q p  E p E q Y pp sin  pp   E p E q Y pq sin  pq   p   q 

Therefore the elements of Jacobian (ie H & L) can be calculated as from the equations above of power. OFF diagonal element of H is Pp H PQ  sin  pq   p   q   q

Flowchart: Start

Read the input data

Form the Y bus matrix

Form B’ and B” matrix

Set flat voltage profile except for slack bus

Set convergence criterion є C

Set iteration count p=0

Calculate real & reactive power

Calculate

 P   H Q    N   

 P   Q   ,    V   V 

M     L  V 

No

Are

 P   Q   ,    V   V  Yes

Calculate real & reactive line flow, bus powers

Find   & V  by solving the equations:      1  2    V     V        3  4  Q  Uptate voltage magnitude t phase angles

V   V   V          new

Print the results

new

old

old

Stop C

Algorithm: Step 1: Read the slack bus voltages, real bus powers and reactive bus powers, bus voltage magnitudes and reactive power limits. Step 2: Form the Y bus matrix without line charging admittance and shunt admittance. Step 3: Form B matrix, form Y bus matrix obtained in step 2. Step 4: Form Y bus matrix with double the line charging admittance. Step 5: Form B” matrix from Y bus matrix obtained in step 4. Step 6: Calculate the inverse of B’ & B” matrices. Step 7: Initialize the bus voltage. Step 8: Calculate [∆P/|V|] , [ΔQ/|V|] Step 9: If ΔP/ |V| & ΔQ/|V| are less than or equal to tolerance limit, solution has convergence and go to step 12 otherwise increase iteration count and go to step 10. Step 10: Calculate [Δδ] = [B’]-1 [ΔP/ |V|] [Δ|V|] = [B’’]-1[ΔQ/|V|] Step 11: Update [δ] & [|V|] for all buses except slack bus. [δ]new = [δ]old + [Δ δ]; [|V|]new = [|V|] old + [Δ|V|] Step 12: Compute slack bus power, line flows, real power loss, reactive power loss etc. Sample Problem: For the system shown in Figure 4.4 determine the voltage at the end of the 1st iteration by FDLF method. The line reactances are marked in the figure. 1

j 0.1

2 j 0.2

j 0.2 Figure 4.4

3

Bus specifications: Bus code 1 2 3

Assumed volt. Generation Load P Q P 1 1.1 5.52 0≤QG2≤5.5 0.7 3.65

Q 0.1 0.53

Solution: P2 = PG2-PD2 = 4.82 p.u P3=PG3-PD3 = -3.65 p.u Q3= QG3-QD3= -0.53 p.u  j15 j10 j 5  15  90   Ybus   j10  j15 j 5   1090  j 5 590 j 5  j10 15  5 B'     5 10

1090 590  15  90 590  590 10  90

B''  10

Flat voltage profile: V10  1  j 0  1.00 0 V20  1.1  j 0  1.10 0 V30  1  j 0  1.00 0

Calculation of P and Q: Pp 

n

 V p Vq

q 1

Y pq cos( p   pq   q )

P20,cal  0 P30,cal  0 n

Q p   V p Vq Y pq sin( p   pq   q ) q 1

= V2 .V1 .Y21 sin  21   1   2   V2 .V2 .Y22 sin  22   V2V3Y33 sin  23   3   2 

= -1.65 As 0≤QG2≤5.5 Q2= QG2-QD2 QG2=-1.65+0.1 = -1.55 Hence it is not within the specified limits. Q3= V3 .V1 .Y31 sin  31   1   3   V3 .V2 .Y32 sin  32   2   3   V3V3Y33 sin  33  = 0.5 Calculation of change in power: P2  P2.spec  P2.cal  4.82  0  4.82 P3  P3.spec  P3.cal  3.65  0  3.65 Q3  Q3.spec  Q3.cal  0.53  0.5  1.03

Find the largest value of ΔP2, ΔP3, ΔQ3 Let the largest change of ΔE= 4.82 ΔE≤ E; 4.82≤0.01 Find Δδ and V : 



   B ' 1  P   V 

 V   B 

'' 1

B 

' 1

B 

'' 1

 Q     V 

15   5 0.08 =  0.04

5 10

1

0.04 0.12

 1 / 10  0.1

   

0.08

0.04

0.04   P    0.12  V 

 2  0.08 0.04 4.38     0.04 0.12  3.65    3  Δδ2=0.2045

Δδ3=-0.2627 



 V   0.1 Q   V 

V3  0.103 Find the new values of phase angle and magnitude of the voltage:

 i p 1   i p   i p p 1

p

p

Vi  Vi   Vi p =0; i=2,3  21   20   20  0.2045 rad  11.710

 31   30   30  0.2627 rad  15.05 0 p =0; i=3 1 0 0 V3  V3   V3  1  0.103  0.897

New values are: V1'  V1'  1'  10 0 ( Slack bus ) V2'  V2'  2'  1.10 0 (Generator bus ) V3'  V3'  3'  0.897  15.05 0 ( Load bus )

Result: The load flow study on the given power system using Fast decoupled method was conducted using MATLAB and results was verified.