CONTINGENCY ANALYSIS OF POWER SYSTEM NETWORKS

By

MOHD HAFIZ BIN ABDUL RAUF

FINAL PROJECT REPORT

Submitted to the Electrical & Electronics Engineering Programme

in Partial Fulfillment of the Requirements
for the Degree
Bachelor of Engineering (Hons)
(Electrical & Electronics Engineering)

Universiti Teknologi Petronas

Bandar Seri Iskandar
31750 Tronoh
Perak Darul Ridzuan

© Copyright 2006
by
Mohd Hafiz bin Abdul Rauf, 2006
 CERTIFICATION OF APPROVAL

CONTINGENCY ANALYSIS OF POWER SYSTEM NETWORKS

by

Mohd Hafiz bin Abdul Rauf

A project dissertation submitted to the

Electrical & Electronics Engineering Programme
Universiti Teknologi PETRONAS
in partial fulfilment of the requirement for the
Bachelor of Engineering (Hons)
(Electrical & Electronics Engineering)

Approved:

Professor Dr Ramiah Jegatheesan

Project Supervisor

UNIVERSITI TEKNOLOGI PETRONAS

TRONOH, PERAK

December 2006

HI
 CERTIFICATION OF ORIGINALITY

This is to certify that I am responsible for the work submitted in this project, that the

original work is my own except as specified in the references and

acknowledgements, and that the original work contained herein have not been

undertaken or done by unspecified sources or persons.

IV
 ABSTRACT

In this modern power system, planning and construction of the system with outage
minimization as a major design criterion is becoming more important. However,
even with the best possible construction and operating practices, outages still can
occur. Several techniques have been developed to analyze the status of the system
after the real outage occurs. This project presents an alternative way of conducting
power flow analysis in determining the new bus voltages and line currents during
outage occurrence. The algorithms have been accommodated in one complete
program making it an effective tool for contingency analysis. The program has been
tested on an IEEE bus test case and the remarkable results obtained in terms of its
reduced computation time and accuracy showed that this technique could be
implemented, particularly when involving practical work with a lot of calculations.
 ACKNOWLEDGEMENTS

This dissertation could not have been written without Dr. Ramiah Jegatheesan who
not only served as my supervisor but also encouraged and challenged me throughout
my academic program as well as this project. He patiently guided me through the
dissertation process, never accepting less than my best efforts.

Finally, I would like to thank my mother and father, who gave me the opportunity and
the spirit to educate myself. Mom, Dad, you are the best. To all my sisters and little
brother, I hope I can serve as inspiration to all of you. If I can do it, so can you.

VI
 TABLE OF CONTENTS

Certification of Approval ... iii

Certification of Originality ... iv
Abstract ... v
Acknowledgement ... vi

List of Illustration ... vii

List of Abbreviations ... ix

1.0 Introduction ... 1

2.0 Literature Review ...3
2.1 Bus Impedance Matrix ZbUS ... 3
2.2 Power Flow Analysis ... 5
2.3 Graphical User Interface ... 6
3.0 Methodology ... 7
3.1 Bus Admittance Matrix ... 7
3.1.1 Bus Admittance Matrix building algorithm ... 8
3.1.2 Factored Matrices of Ybus Matrix ... 10
3.1.3 Solving Ybus Vbm = hasfor Vbm using the factor matrices ... 12
3.1.4 Summary of formulas ... 14
3.1.5 Calculation ofthe vector Z£"n)rrom
factors of Ybus matrix ... 14
3.2 Addition of one line ... 15
3.3 Removal of one line ...19
3.4 Addition / Removal of two lines ... 19
3.5 Half Line Charging Admittance ... 22
3.6 Tap Setting for Transformers ... 23
3.7 MATLAB Programming ...24
3.7.1 Program Modules ... 25
3.7.2 Creating MS Excel Data File ... 26
3.7.3 Exporting Excel Data into MATLAB ... 27
3.7.4 Error Handling ...29
3.7.5 Counter (Stopwatch Timer) ... 29
3.8 Tools ... 30
4.0 Results and Discussion ... 31
4.1 Results and Analysis ... 31
4.1.1 Removal of 33rd element (Bus 24 - Bus 25) ... 35
4.2 Graphical User Interface ... 41
5.0 Conclusions and Recommendation ... 45

References ... 46
Appendices ... 47

vii
 LIST OF ILLUSTRATIONS

List of Figures

Figure 1: A simple GUI using MATLAB ... 6

Figure 2: Original system with voltages Vi, V2,..., Vn due to current
injectionIit I2, ..., In ... 15
Figure 3: Addition of line ... 16
Figure 4: The effect of current Ia ... 16
Figure 5: Half Line Charging Admittance ... 22
Figure 6: The tap setting for transformers (I) ... 23
Figure 7: The tap setting for transformers (II) ... 24
Figure 8: The flowchart of the MATLAB program ... 24
Figure 9: The Program Modules ... 25
Figure 10: The IEEE 30 bus line data in MS Excel under "Edata" sheet ... 26
Figure 11: The IEEE 30 bus pre-outage bus voltages
under "preout_V" sheet ... 27
Figure 12: The matrix ... 28
Figure 13: The warning dialog ... 29
Figure 14: The MATLAB software ... 30
Figure 15: Comparison of the computation time (in seconds) between
four (4) different methods of power flow, including Contingency
Analysis Program whenever line 24 - 25 was removed.
Contingency Analysis program indicates remarkable results in
term of its computation time reduction. ... 36
Figure 16: The main page of the program 41

Figure 17: The Single Contingency Section 42

Figure 18: Select the correct data 43

Figure 19: Enter the input and click CALCULATE! 43

Figure 20: The displayed results 44

vin
List of Tables

Table 1: Summary of the Bus Impedance Matrix algorithm formulas ... 4

Table 2: The Line Data for IEEE 30 Bus Systems ... 32 - 34

Table 3: The initial voltages at each bus after conducting the

power flow analysis ... 34 - 35
Table 4: Comparison of the computation time (in seconds) between
four (4) different methods of power flow, including
Contingency Analysis Program in table form ... 36

Table 5: The new voltages at each bus when a line is removed ... 37 - 38

Table 6: Comparison between Power Flow Program and

Contingency Program during the evaluation of line flows
on each line whenever line 24 - 25 was removed ... 39 - 40

IX
 LIST OF ABBREVIATIONS

IEEE Institute of Electrical and Electronics Enginners

MATLAB Matrix Laboratory (Mathworks, Inc.)
GUI Graphical User Interface
AC Alternating Current
DC Direct Current

AEP American Electric Power System

Zbus Bus Impedance Matrix
Ybus Bus Admittance Matrix

HLCA Half Line Charging Admittance

RAM Random Access Memory

x
 CHAPTER 1

INTRODUCTION

Line outages, whether it is planned or unplanned are the common issue in power
system networks. That is why planning and construction of the system with outage
minimization as a major design criterion is becoming more important and more
common. Of course, the primary requirements of a reliable power system are the
proper line construction and maintenance and proper over-current device placement
and coordination [1]. Nevertheless, even with the best possible construction and
operating practices, outages still can occur.

Whenever an outage occurs, an effective and economical way to improve the

reliability of the power system is by restoring the power temporarily using alternative
feeds [1]. However, what it more important is to determine the status of the system
when any possible line outages occur. This is the point where many techniques of
power system analysis have been introduced and developed, such as power flow
analysis, contingency analysis, system restoration, load shedding, etc [2].

A system contingency can be defined as a disturbance that can occur in the

network and can result in possible loss of parts of the network like buses, lines,
transformers or power unit in any of the network areas. Contingency analysis is an
evaluation of an existing or proposed electrical distribution system to determine the
capability of the system to restore power using alternative feeds [1]. It is performed
to assess the status of the power system network for any possible line outage before it
is actually occurs.

When a line is switched onto or off the system through the action of circuit
breakers, line currents are redistributed throughout the network and bus voltages
change. These new steady steady-state bus voltages and line currents can be
predicted quickly by what is called the contingency analysis program [6]. The large-
scale network models used for this type of evaluation, like those used for fault
calculations, do not have the be exact because the system planners and operators,
who must undertake hundreds of studies in a short time period, are more concerned
with knowing if overload levels of current and out-of-limit voltages and power flow
exist than with finding the exact values of those quantities [6].

This project aimed at developing the Contingency Analysis program, to

calculate new steady state bus voltages and line currents during the line outage and
other contingencies. The program uses the factored matrices of Bus Admittance
Matrix algorithm, which is the main advantage of this program. It is written and
demonstrated in MATLAB. In realizing the effectiveness of the developed program,
validation was performed on an IEEE Bus System and the results were compared to
the results using the traditional power flow analysis and it indicates a remarkable
advantage particularly in its computation time and the accuracy. Results obtained
from the studies indicated the proposed technique could be implemented, particularly
when involving practical work with a lot of calculations.
 CHAPTER 2

LITERATURE REVIEW

Power system network can be characterized by two methods namely, Bus Impedance
Matrix, Zbus and Bus Admittance Matrix, Ybus, which are inverse of each other.
Practically, it is much easier to construct the bus admittance matrix, rather than bus
impedance matrix.

2.1 Bus Impedance Matrix Zbus [4]

Bus Impedance Matrix, Zbus of a power network can be obtained by inverting the bus
admittance matrix, Ybm, which is easy to construct. However, when the order of
matrix is large, direct inversion requires more core storage and enormous computer
computation time. Therefore, inversion of Ybus is prohibited for large size network.

Zbus can be constructed by adding the network elements one after another.
Using impedance parameters, performance equations in bus frame of reference can
be written as:

J^bus ^bus *bus (1)

In expanded form, the above becomes

7
E, Zn ~h '
7
'21 '22 ^2N h
(2)

hN Lm ZN2 7
_h_

From this we can write

Ev - ZBlIl + Zrjlh + + Z„aL

•pq*q + + ZnNll
'spN-lN (3)
By injecting 1 p.u. current at bus q and other bus currents are set to zero, we can
obtained the Zm by measuring Ep. The indices p and q can be varied from 1 to N.

When constructing Zb„s building algorithm, elements are added one by one.
At any stage, the added element may be a branch or a link. A branch is the added
element which will create a new bus in the network. On the other hand, by adding
link, no new bus will be created.

As given in [4], the summary of the Bus Impedance Matrix algorithm

formulas are as Table 1 below:

p is not reference bus p is the reference bus

Added element p-q is

a branch
Zqi ~ Zpi Zqi-0

i = l,2, ..... m i = 1, 2, ...., m

ifq i+q

7 =7 +7 7 =7
11° PI a 99 a

Added element/?-(/ is
a link
Zj\i —£,pj—/jqj Zli =~Zqi

i - 1, 2, ...., m i = 1, 2, ..... m

Zu = Zpl ~ Zql + Za Zll =-Zgl + 2a

7 _ z,z,
^ij(modified) Zijfbefore modification)
Zll

i = 1, 2, ...., m and j = 1, 2, ...., m

Table 1: Summary of the Bus Impedance Matrix algorithm formulas

2.2 Power Flow Analysis

Power flow analysis, also known as load flow, is an important part of power system
analysis. It serves the purpose of planning, control of existing system as well as
planning its future expansion [5]. Unlike traditional circuit analysis, a power flow
study usually uses simplified notation such as a one-line diagram and per-unit
system, and focuses on various forms of AC power (i.e. reactive, real, and apparent)
rather than voltage and current.

The great importance of power flow or load-flow studies is in the planning

the future expansion of power systems as well as in determining the best operation of
existing systems. The principal information obtained from the power flow study is
the magnitude and phase angle of the voltage at each bus and the real and reactive
power flowing in each line.

When conducting a power flow analysis, the system is assumed to be

operating under balanced conditions and a single phase model is used. Four
quantities are associated with each bus. These are voltage magnitude |V|, phase angle
5, real power P, and reactive power Q. Generally, the system buses can be classified
into 3 types: [5]

1) Slack bus
It also known as slack or swing bus, which will be the reference bus. In
slack bus, voltage magnitude and phase angle are specified.

2) Load buses (P - Q buses)

In these buses the active and reactive power are specified.

3) Regulated buses (P - V buses)

These buses are the generator buses, where the real power and voltage
magnitude are specified.

There are several solutions, or methods introduced to solve for the power flow
analysis problem, namely:

1) Gauss-Seidel Method
2) Newton-Raphson Method
3) Fast Decoupled Method
2.3 Graphical User Interface (GUI)

A graphical user interface (GUI) is a ^aphical display that contains devices, or

components, that enable a user to perform interactive tasks. To perform these tasks,
the user of the GUI does not have to create a script or type commands at the
command line. Often, the user does not have to know the details ofthe task at hand.

The GUI components can be menus, toolbars, push buttons, radio buttons, list
boxes, and sliders. In MATLAB, a GUI can also display data in tabular form or as
plots, and can group related components. The following figure illustrates a simple
GUI.

This project will utilize this MATLAB's feature to make the program more
user-friendly and thus make it easier for conducting the analysis.

1••[' simple_gui 4pH

Surf |
1CU- ••"*"
^^^^k.
Mesh j
0-
Contour |

-10,
W---'~'-
40 ^\^ Select Data
40
20^ 20 peaks *
0 0

Figure 1: A simple GUI using MATLAB

 CHAPTER 3

METHODOLOGY

In contingency analysis, the network models used do not have to be exact because
the power engineers and system operators are more interested in knowing whether or
not an insecure or vulnerable condition exists in the steady state following any of the
outages.

Therefore, several approximations are made. Often, the resistance is

considered negligible and the network model becomes purely reactive. Line charging
and off-nominal tap charging of transformers are also omitted [6]. Accordingly, to
test for the effects of line and transformer outages on the bus voltages and line flows
in the network, approximate ac power-flow techniques are generally employed since
they can provide a fast solution of the many test cases which need to be run.

In many cases, linear models are considered satisfactory and the principle of
superposition is then employed. Methods of contingency analysis which used the
factored Bus Admittance Matrix (Yb,J becomes attractive from computational
viewpoint, especially if the loads can be treated as constant current injections into the
various buses of the system [6]. Hence, the concept of compensating currents is
developed, to allow the use of existing system Zbus without having to modify it
whenever a contingency occurs.

3.1 Bus Admittance Matrix

While doing symmetrical short circuit analysis, the bus impedance matrix Zbus is
employed. However, when the system is large, construction of Zbus matrix itself is
tedious. Furthermore, since Zbus matrix is always full, more computer storage is
required to store the matrix. Moreover, modification of the system such as removal
or addition of one link calls for recalculation of all the elements in Zbus matrix.
 On the other hand, the bus admittance matrix Ybus can be easily obtained
even for large power system network. The main advantage of this matrix is that, it is
always sparse (having many elements of zero) and thus requires less storage
capacity. For instance, whenever an element a - b is added, only four elements 7^,
Ybb, Yab and Yba need to be modified and the other elements in the Ybus matrix remain
unaltered.

3.1.1 Bus Admittance Matrix building algorithm [6]

This matrix can be build by adding one element at a time. To start with, all the
elements of bus admittance matrix are assumed to be zero. Then, when a network
element of admittance ya is added between buses p and q, ONLY 4 elements of bus
admittance matrix are to be modified as below.

•* pp new ~ *pp old "*" yu

%qq new ~ -* qq old ~*~ J^«

*pq new ~ *pq old ' yu (4)

Iqpnew ~ 1qpold ' yti

Even with the orientation is from bus q to /?, the result will be the same.
However, if busq happens to be reference bus,then Ypp alone getsmodified as:

Ippnew ~" fppold'ya (p)

i.e. if ypo is connected between bus^ and reference bus, then

Ypp new = Ypp 0id + ypo (6)

This method of constructing the bus admittance matrix by adding the

elements one by one and modifying the bus admittance matrix each time, is very well
suited for large network. Each time, one network element data is read and the bus
admittance matrix is updated. When all the network elements are added, the final bus
admittance matrix will be obtained.

When any fault occurs atpth bus, only the//* column ofthe Zbus matrix and
not the full Zbus matrix is needed. However, it is also possible to obtain the p'h
column of Zbus matrix, Z\P2 from the Ybus matrix. In fault study, the bus at which the
fault occurs (value of/?) may keep changing. Therefore, if the Ybus is factored and the
factored matrices are available, any required column of Zbus matrix can be obtained
easily.

The // column of the Zbus matrix, Z[p] can be obtained from the factor
matrices of Ybus matrix. Previously, we know already that

Zbus *bus 'bus (7)

Thus,

z„ 2a •- ^ ••• 7 -
A yt
7 7 • 7
h v1
" Zl.' ' ^2N

.. 7 •• 7A pN
—

(8)
zu • h y.
zr
7 7 .. 7 •• 7
^Np .V y*

Consider the bus current vector:

~v "0"
h 0

*bus ~
—

and denote this as 1, (9)

h 1

0
Jo.

Substituting Ibus = \p in the equation (8), we will get

"V >l" '^1

z?' y,
~
i e 7lp) (10)
yF VD
zr
7
^Np Y«. K,

Thus, by injecting 1 p.u current at bus/?and set the othercurrents to be zero (Ibus=lP)
Z12 is numerically equal to Vbm obtained from (7). The inverse form of (7)

ls,1busrVbus - I
isF 1bus-
 In the subsequent chapters, we will see how the Ybus is factorized and how
the factored matrices are used to obtained the//* column ofthe Zbus.

3.1.2 Factored Matrices of Ybus Matrix [6]

This is the important part of the program. Once the factored matrices of Ybus matrix
are obtained, any required columnof the Zbus can be extracted easily.

The symmetrical matrix Ybus can be factored as:

Ybus = LDLl (11)

where I is a lower triangular matrix

D is a diagonal matrix

V is the transpose of L matrix

Because of Ybus = LDL', we can state that:

Z{bll is numerically equal to Vbus obtained from LDLlVbus =Ibus, for Ibus = \p.

If the fault occurs at some other bus, say k, then ZbJk) can be obtained by
solving LDL'Vbus = Ibus for Vbus taking Ibus = \k. This calculation requires only back
substitution as the factor matrices are already available. Thus, only a repeat solution
is carried out with different values for right hand side vector.

Let us consider a 4 bus system described by:

Y
"Ml Y
l2\ Y
12\ Y
J41 \y,~ ~h
Y y Y Y
J12 J22 J32 "*42 y2 h
(12)
Y
J13 Y
J23 Y
J33 Y
J43 y, h
Y
/14
Y
J24
Y
J34
Y
J44 k. j,.

The coefficient matrix Ybus can be written as the product of three matrices L, D and

Ll. Therefore, Ybus = LDL'. (13)

10
Matrix L is a lower triangular matrix of the form

1 0 0 0

*"21 1 0 0
1 = (14)
*3. -L32 1 0

*41 ^42 ^43 1

and matrix D is a diagonal matrix of the form

dn 0 0 0

0 d22 0 0
D = (15)
0 0 ^33 0

0 0 0 dA44

The matrix L' is only the transpose of matrix L.

1 *- 21 '31 *41
0 1 c32 C42
L' = (16)
0 0 1 •°43

0 0 0 1

Thus, using equation (14) and (15) in equation (13), Ybus can also be written as:

1 0 0 0 du 0 0 0 1 **21 "^31 '4. Y

-*11 Y
J21 Y
J31 Y41
AA

-c21 1 0 0 0 ^22 0 0 0 1 c32 Y 22 Y 32 Y 42

"M2

*31 ^32 1 0 0 0 ^33 0 0 0 1 ^43 Y Y 23 Y Y 43

*41 ^42 ^43 1 0 0 0 <4 0 0 0 1 Y•14 Y 24 -34 r.44

After multiplying first 2 matrices, the result as follows:

d„ 0 0 0 1 £4 21 F Y Y Y
^31 '41 "Ml J21 J31 J441

*2i"n ^22 0 0 0 I ^32 -c42 r„

•11 E,
32 7,42
(17)
^ 3l"ll * 32^22 ^33 0 0 0 1 ^43 Y„ •23 K,
'33 K43

'4A <-42w22 •C43M33 V 0 0 0 1 7,„

-24
E,
34
K'44

11
By looking to the above equations, we can generalize the equations to obtain the
elements of the matrices. The equations are as follows:

dll=rll-fjelt2dkt (18)
k=\

j-\

^rPi-IW*)^ (19)
k=\

The elements of the factor matrices L and D can be obtained using equations (18) and
(19). For matrix of order N, the elements of factor matrices L and D are calculated in
the order:

du^M^Mi^ Jm fori = 1,2,3, ,N (20)

The L and D matrices are combined into one matrix only, to reduce the computer
storage.

dn 0 0 0
£2l d22 0 0
(21)

f
/Ml fM2 fM3 d
u44

3.1.3 Solving YbmVbus = Ibus for Vbus using the factor matrices [6]

Now the problem is to solve for vector Vbus for a given value of vector Ibus using

LDL'Vbm=Ibm. (22)

Let us define two intermediate dummy vectors V and V" as

V - V (23)

DV = V" (24)

Then LV" = Ibm (25)

In order to solve for the unknown voltage vector Vbus, the three equations above must
be solved in the order equation (25), (24) and (23).

12
For a 4 bus system, equation (25) is:

du 0 0 0
yl [7,1
Ml <M2 0 0 h
=
(26)
^33 0
M, M2 y;
M! M2 M3 *« k". j*\
The elements of V" vector are calculated using the equation below:

i-i

r."=.r,-EV» ' =1.2, AT (27)

t=l

Now equation (24) can be written as:

dn 0 0 0 y[ vl
0 ^22 0 0
=
yi (28)
0 0 d33 0
y; y,"
0 0 0 d<44 y. yt

From this above equation, it is clear that the elements of V vector are calculated
using the equation below:

V'=Vt"/dit i-1,2, ,N (29)

Finally, equation (23) can be written in the form of

1 £
Ml
£
Ml M, \y,~ y[
0 1 £ 42
M2 *- y2 yi (30)
0 0 1 M3 y, y,
0 0 0 1
[y*\ y*_

From this equation, we can generalize that

JVl

v,-v;-^kivk \ = N,N-i, ,1 (31)

t=(+i

It is to be noticed that, while calculating elements of vector Vbus, from the above
equation, V, s are computed in the order of i = N, N-\, ,1

13
 3.1.4 Summary of formulas

Elements of factor matrices are calculated from

dn =Yh-Z**2** and £n =(7, -2X-A^)/</,

*=i k=i

following the order duJl+l,£j+2l, JNi fori-1,2,3, ,N

Vbus solution is obtained from these equations:

i-i

vf=h-Y<t*v* i=1'2' ;N (32)

k=\

v;=vrid„ i = l,2, ,N (33)

K^-E1^ i=N,N-l, ,1 (34)

jt=f+i

3.1.5 Calculation of the vector Z[™;n) from factors of Ybus matrix [6]

In contingency analysis, if a contingency occur between line m and n, we need to

calculate the difference ofthe two vectors (Z£> - z££ ) where Z£> and Z£] are
the m"1 and m** column ofthe bus impedance matrix Zbus, respectively.

Let Z[m-n)=
bus
Z[m)-
bus
Z{hn)
bus (35)

Again let us define vector 1(m.n) as

1 0

m 1

(m-ri) (36)
0

n -1

N 0

14
In Vbus = ZbusIbus, if Ibus = l(m_n) (1 p.u current is injected at bus mand -1 p.u current
is injected at bus «), then the resultant vector will be the difference of the mih and rih
column of the Zbus, respectively.

Thus z£r° is numerically equals to Vbus, obtained from Vbm-=ZbJbus

for/,„s =l(m_fl). Equivalently, Z^ can be obtained solving YbusVbus =Ibus for Vbus
corresponding to Ibus=\m.n).

Knowing that Ybus can be factored as LDL', we can now state that Z(™~n) can

be obtained by solving LDL'Vbus =Ibus for Vbus corresponding to Ibus =\m.n).

3.2 Addition of one line [6]

When considering line addition to or removal from an existing system, it is not

always necessary to build a new Zbus or to calculate new triangular factors of Ybus,
especially if the only interest is to establish the impact of the changes on bus
voltages. This is where the compensating currents concept is developed as an
alternative procedure. The compensating currents are injected into the existing
system to account for the effect of the line changes.

Consider the original system with voltages Vi, V2, ..., VN due to current
injection//,/^ -, In as illustrated in figure 2 below:

Figure 2 : Original System

15
We assume that the bus voltages Vj, Vj, ..., Vn produced in the original system by the
current injections I], h, ..., In are known. Suppose that a line of impedance za is
added between buses m and n. The current injections are fixed in value and therefore
are unaffected by the addition of line of impedance za. Let the changed bus voltages
be Vi\ V2', .... Vm', ..., Vn', ..., Vn' a shown in Figure 4 on the next page.

Figure 3 Figure 4

Addition of line causes a current flow of Ia as shown in Figure 3. We have

Za*a ~ "bus "m

K
m n

[0...0 1 0...0 -1 0...0] K (37)

v'

16
Defining
m n

Ac=[Q...O 1 0...0 -1 0...0] (38)

we get

ZJa —AcVbus (39)

The effect of current Ia shown in Figure 3 can be replaced by current injection of - Ia

at bus m and Ia at bus n as shown in Figure 4. By examining Figure 4, it reveals that
the bus voltages Vj\ F/, ..., ^V are due to:

1) Current injections h, I2, ..., Inand

2) Current injections - laat bus mand Ia at bus n

Hence, to determine the bus voltages V}', V2', ..., VN' we shall now apply the
Superposition Theorem. When current injections /;, h, ..., In alone are there, bus
voltages will be the same as the bus voltages in the original system (refer Figure 2)
denoted by Vbus- When current injections -Ia at bus m and Ia at bus n are alone there,
bus voltages will be equal to ZbusICOmp where

' 0 " "0"

0 0

-h 1

0 0

comp / = -4L (40)

0 0
_1
/.
0 0

0 0

Therefore;

Miits' 'bus ~>~ Zbus^comp 'bus ZbusAcla (41)

Vbus can be obtained combining equation (39) and equation (41).

17
Substituting equation (41) into (39),

zJa=AVbus-AcZbusA'Ja (42)

Defining Z = za+ AcZbusA'c (43)

we get ZIa =Vm-Vn, thus we get

i.=<y.-y.vz (44)

Now equation (41) can be written as

Vbus- = Vbus + kV (45)

Where AV=-ZbusA'Ja (46)

Thus, new bus voltages Vbus can be written as follows:

1) Calculate Z from equation (43)

2) Calculate Ia from equation (44)
3) Calculation AVfrom equation (46)

4) Calculation Vbus from equation (45)

Note that equation (42) and (45) contain Zbus, the bus impedance matrix of the
system, which is not available. Let us now see how this could be overcome.

Calculation of Z

Let us note that ^c=[0... 0 1 0 ... 0 -1 0...0]

m n
Z = za+ AcZbvsA'c

= za+AcZ^n) (47)

Denoting m[Z^;n) ] is the mth element of Z£""} and n[Z^"n) ] is the nth element of
Z{b™~n) the above equation becomes

Z-za + m[Z<r)]-n[ZSr)] (48)

18
Calculation of/,,

Ia=(Vm-Vn)/Z (49)

Calculation of AV

AV=~ZbusA'Ja (50)

7(m-«) t

Calculation of new busvoltage V'bus

yL,=K.+*y (5i)

It isto be noted that for the above calculations we need only the vector Z£"H) and not
the bus impedance matrix ZbHS. We should use the above four equations (48 - 51) to
calculate Z, la, AVand Vbus.

3.3 Removal of one line

In most cases, line outage is more common due to fault or over loading. Hence it is
more pertinent to study removal of line rather that addition of line. Accommodating
removal of line is not at all difficult. Removal of a line m - n with impedance za is
equivalent to addition of line m - n with impedance - za.

3.4 Addition / Removal of two lines

The method discussed can be extended to addition or removal of two lines. Let us

say we have a power system for which the bus voltages are V\\ V2', ..., Vn'. Two
lines of impedance za and zb are added between buses m-n and p-q respectively. It
is required to find the new bus voltages Vj', V2', ..., Vn'. The currents in the added
lines will be Ia and Ib. Then

z.I.=Vm-K<mdzbIt=Vf-V;
i.e.

19
 V

v„

m n
V.
Za 0" • 0 1 0 ••• 0 -1 0 •• 0
=

0 ••••0 1 0 0 -1 0 • •• 0
F„

o p 1

V.

Vy

-AKus (52)

For this case the compensating current is

' 0 " "0 0"

m
-h -1 0

P -h 0 -1
~I~
a ~L~ a

comp
0 =-4 c
(53)
1 0
h. h_
n
/.
0 0 0

q h 0 1

0 0 0

Using Superposition Theorem, new bus voltages can be written as

Vbus =Vbus +7 bus I camp =Vbus - 7bus A'c (54)

Using the above equation in equation (51), we get

zna 0" ~I~ a

~I~
a

= AcVbm AcZbusAc (55)

0 zb_

i.e.

20
 Vm -Vn
+AZbuX = Avbus (56)
0 z. V -V
P M.

Defining Z = + 42*4 (57)

0 z.

The above equation becomes

V -V
i.e.
V -V
p q

V -V
= Z (58)
V -V
. p 1.

Then, the new bus voltages can be written as

Vh„,. = Vb„s + AV (59)

Where &V=-ZbusA'c (60)

Thus, new bus voltages Vbus can be determined as follows:

1) Calculate Z from equation (57)

2) Calculate from equation (58)

L

3) Calculation AV from equation (60)

4) Calculation Vbus from equation (59)

It is noted that equation (57) and (60) contain Zbus, the bus impedance matrix of the
system, which is not available. Let us now show that how this problem can be
overcome.

Calculation of Z

z =
z„a 0"
+
\m~n)[Ztn)] {m-n)[Z£*i (61)
0 V (p~q)[Z(b7n)] (p-q)[Z™]

21
Calculation of

V -V
= Z (62)
V -V
p <t.

Calculation of AV

AV=-Zbll,A'c

— 7,Cm_B) 7iP-<l)
{L
= -z bus 'bus (63)

Calculation of a new bus voltage /',bus

Vbus=Vbus+AV (64)

It is to be noted that for the above calculations, we need the vectors Z^ w) and
Z^^only and not the bus impedance matrix Zbus- We should use the above four

equations (61-64) to calculate Z, ,AVandK'

3.5 Half Line Charging Admittance [5]

Consider the diagram below:

Figure 5

22
Line current from / ->j is given by

W/+J»=*(m-m-) +*om (65)

Similarly, the line current1$ measured from busy* to bus / and defined positive in the
direction/ -> / is given by

iJI=-i,+iJo=ys<yJ-y,)+yJyJ (66)

The complex powers Syfrom bus / toy andSyfrom busy to / are

s»=y,n (67)

(68)

3.6 Tap Setting for Transformers [5]

Tap changing transformers as well as regulating transformers can be use to control

the real and reactive powers. In a tap changing transformer, when the ratio is at the
nominal value, the transformer is represented by a series admittance yt in per unit.
With off-nominal ratio, the per unit admittance is different from both sides of the
transformer, and the admittance must be modified to include the effect of the off-
nominal ratio.

1
v, yt vx !
A ) ('
) C
l:a

Figure 6

y,
a
(69)
y, y< vt
a \a\

For the case where a is real, the n model as shown in the figure below
represents the admittance matrix above. In the it model, the left side corresponds to
the non-tap side and right side corresponds to the tap side of the transformers.

23
 Figure 7

3.7 MATLAB Programming

A program is written, based on the factored matrices of bus admittance matrix

building algorithm as explained before. The programming is written using
MATLAB. The flows of the program are illustrated in the flow chart below:

Read Data from MS Excel

(Line and Bus Data)

i '

Construct Ybus and its factor matrices

' '

User input which line to be remove?

1 '

Calculate

' '

Display Results

Figure 8: The flow chart of the MATLAB program

The line data and initial bus voltages are entered in Microsoft Excel spreadsheet
format. Then, the data are imported into MATLAB workspace for further computing.

24
 3.7.1 Program Modules

This program is comprises of several modules, which means it has several

functions that working together to achieve the final results. The advantages of
writing the program in modular structure are:

1. Each module of function can be run separately, which make it reusable

and allow the programmer to work on different modules simultaneously
2. It promotes modularity which can make complex program easier to
understand

3. The easiness of doing testing and troubleshooting, when any problem

arises

Figure 9: The Program Modules

25
 3.7.2 Creating MS Excel Data File

Before starting the analysis, the line data as well as the pre-outage voltage are to be
stored inside the MS Excel. The main purpose for this is, the user will find it much
easier to manipulate the data when they are in Excel form, since they can be changed
without having to open the MATLAB.

^.,S'dl "0:' 'tf "-2) ••"' !EEE30bui„SjSdsls,ili 'iCiiiiiiJ.ilrUilUy l.!-:'J< -Mrcicioft E*£ii
• ' Home Insert Page layout Formulas Data Review View iWd-lns
|'V"i] |j Colors ' j hi \~~~i - „_J | ^SLi GfWl*h: Automatic • ,—j Gridlinss Headings
lil^Jt [gfonts " —' "—J '—' -!=*' 1 -'-:,: ~=A £J Htijjht Automalie - •J View J View
Themes r—| Margins Qnerrtation Sue Punt Break; Background Print ., . , Custom
. . 0 Effects ' - Afea „ . Title; a} Scale:. lOOfc . Vlews Print Print
Th*ntcj r-3ae j*t'Jp r-" S«lito Fit r- Sntst Optii-rii

M9 - /< '
A" " 8 c D E F G \

.. From Bus ToB Z X HLCA TAP SETTING

z 1 2 0,0192 0,0575 0,0264 1

3 1 3 0,0452 0,1852 0,0204 1

4 2 4 0,057 0S1737 0,0184 1

5 3 4 0,0132 0,0379 0,0042 1

6 2 5 0,0472 0,1983 0,0209 1

? 2 6 0,0581 0,1763 0,0187 1

3 4 6 0,0119 0,0414 0,0045 1

9' 5 7 0,046 0,116 0,0102 1

10' 6 7 0,0267 0,082 0,0085 1

IV 6 8 0.012 0.042 0,0045 1

12 i 8 9 0 0.208 0 0,978
13 6 10 0 0,556 0 0,969
ui 9 11 0 0,208 0 1

15! 9 10 0 0,11 0 1

16! 4 12 0 0,256 0 0,932

171 12 13 0 0,14 0 1

is1 12 14 0,1231 0,2559 0 1

19i 12 15 0.0662 0,1304 0 1

20' 12 16 0.0945 0.1987 0 1

21: 14 15 0.221 0,1997 • 0 1

22 16 17 0.0824 0,1923 0 1

23 15 18 0.1073 0,2185 0 1

24 18 19 0,0639 0,1292 0 1

25 19 20 0.034 0,068 0 1

26 10 20 0,0936 0,209 0 1

2? 10 17 0.0324 0.0845 0 1
28 10 21 0,0348 0.0749 0 1

29 10 22 0.0727 0,1499 0 1

30 21 22 0.0116 0,0236 0 1

31 15 23 0.1 0,202 0 1

32 f * \ 24 0,115 0,179 0 1

33 f 23 1 24 0,132 0.27 0 1
h " > • Edata • I k u •..;.«,!
'.. a.. :.:x:--...
Ready V^ J • ! :
" - " - - -

Figure 10: The IEEE 30 bus line data in MS Excel under "Edata" sheet

Figure 10 illustrate how the line data is written in MS Excel. The pre-outage
voltages are to be written, also under the same file name, but different sheet. Please
refer to Figure 11 below.

26
 %7V'^%--^V' IEfE30bijs_j,vsdatii.p(l5 ['Lompatibtht/Med?; - Miaosoft Excel

Home Insert Page layout View Add-ins

Aria! -10 -A* >- j iJ' Wrap Text General

4il
Paste , B I iff 9 9: E ^E • ^ Merge and Center - ''*jj§ • % » Conditional
Formatting -
Clipbiiar-i r' Alignment '• Number

P6 X'.""'
A D H

Vbus (real) Vbus (angle)

1,0600 0,0000
1,0430 -5,4970
1,0215 -8.0040
1,0128 -9,6620
1,0100 -14,3810
1,0121 -11,3980
1,0035 -13,1490
1.0100 -12,1150
1.0510 -14,4340
1,0444 -16,0240
1.0820 -14,4340
1,0574 -15.3030
1,0710 -15,3030
1,0425 -16.1980
1,0378 -16,2760
1,0447 -15,8810
1.0392 -16,1880
1,0280 -18.8820
1,0252 -17,0510
1,0292 -16,8520
1,0321 -16,4680
1,0327 -16,4540
1,0272 -16,6610
1,0216 -16,8290
1,0189 -16.4230
1,0013 -16,8400
1.0257 -15,9120
1,0107 -12.0570
1,0059 -17,1360
0.9945 -18.0140

i*l__^:.. vJ!iC
Ready

Figure 11: The IEEE 30 bus pre-outage bus voltages under "preoutJV" sheet

3.7.3 Exporting Excel Data into MATLAB

The line data and the bus voltages can be exported into MATLAB by using a
single command. In order to make it more flexible, the user can choose the desired
file to conduct the analysis. This is achieved by using the command "«/ge(/*/e" in the
program. By using this function, it enablesthe user to select the appropriate data file.
The syntax for this iunction is appended in the appendix.

27
 Then, by using another function "xlsread\ MATLAB will read the data in
MS Excel and stored them in the workspace in a matrix. Before the matrix can be
used for the analysis, it has to be modified slightly, as illustrated in the Figure below.

IEE E3Obus^sysdata

| Element No Fern Bus l b Bus B x hies tap 1

I 1 1 2 0.0192 0.0575 0.0264 1.0000 |

1 2 1 3 0.0452 0.1852 0.0204 1.0000 |
1 3 2 4 0.0570 0.1737 0.0184 1.0000 |
I 4 3 4 0.0132 0.0379 0.0042 1.0000 |
1 5 2 5 0.0472 0.1983 0.0209 1.0000 |
! 6 2 6 0.0581 0.1763 0.0187 1.0000 I
1 1 4 6 0.0119 0.0414 0.0045 1.0000 |
1 8 5 7 0.0460 0.1160 0.0102 1.0000 |
1 9 6 7 0.0Z67 0.0820 0.0085 1.0000 |
1 10 6 8 0.0120 0.0420 0.0045 1.0000 |
1 11 6 9 0.0000 0.2080 0.0000 0.97B0 1
[ 12 6 10 0.0000 0.5560 0.0000 0.9690 1
1 13 9 11 0.0000 0.2060 0.0000 1.0000 1
1 M 9 10 0.0000 0.1100 0.0000 1.0000 1
I IS 4 12 0.0000 0.2560 0.0000 0.9320 |
1 16 12 13 0.0000 0.140Q 0.0000 1.0000 1
I 17 12 14 0.1231 0.2559 0.0000 1.0000 |
. 1 *8 12 15 0.0662 0.1304 0.0000 1.0000 |
j 19 12 16 0.0945 0.198? 0.0000 1.0000 |
| 20 14 IS 0.2210 0.1997 0.0000 1.0000 |
1 21 16 17 0.0824 0.1923 0.0000 1.0000 |
| 22 IS 18 0.1073 0.2185 0.0000 1.0000 |
t 23 18 19 0.0639 0.1292 0.0000 1.0000 |
I 24 19 20 0.0340 0.06BO 0.0000 1.0000 |
I 25 10 20 0.0936 0.2090 0.0000 1.0000 |
i 26 10 17 0.0324 0.0845 0.0000 1.0000 |
I 27 10 21 0.0346 0.0749 0.0000 1.0000 |
I 28 10 22 0.0727 0.1499 0.0000 1.0000 |
] 29 21 22 0.0116 0.0236 0.0000 1.0000 t
1 30 15 23 0.1000 0.2020 0.0000 1.0000 |
i 31 22 24 0.1150 0.1790 0.0000 1.0000 I
i 32 23 24 0.1320 0.2700 0.0000 1.0000 i
I 33 24 25 0.18B5 0.3292 0.0000 1.0000 |
I 34 25 26 0.2544 0.3800 0.0000 1.0000 |
| 35 25 27 0.1093 0.2067 0.0000 1.0000 |
I 36 28 27 0.0000 0.3960 0.0000 0.96B0 |
I 37 27 29 0.2198 0.4153 0.0000 1.0000 |
I 30 27 30 0.3202 0.6027 0.0000 1.0000 |
[ 39 29 30 0.2399 0.4533 0.0000 1.0000 |
I 40 8 28 0.0636 0.2000 0.0214 1.0000 |
j 41 6 28 0.0169 0.0599 0.0650 1.0000 t
1 . . „„
1
»

Figure 12: The matrix

By using the Ybas building algorithm as discussed previously, the Ybas matrix
and its factor matrices are constructed. The construction of the matrix also takes the

half line charging admittance (HLCA) and tap setting for transformer into
consideration. Next, the user will be prompt to input which line to be removed. After
the program has obtained all the necessary data required, it will start calculating for
the new bus voltages, line currents as well as power flow in each line in the system.
The final results are then displayed for clearer view.

28
 3.7.4 Error Handling

In many cases, it is desirable to take specific actions when different kinds of errors
occur. For example, the program may want to prompt the user for more input, display
extended error or warning information, or repeat a calculation using default values.
The error handling capabilities in MATLAB let the program check for particular
error conditions and execute appropriate code depending on the situation.

One of the methods used in the program is warning dialog, or warndlg in

MATLAB. The example ofthe warning dialog box is illustrated below.

^ -laixj

l\ Pressing OK wit dear memoiy

£~2

Figure 13: The warning dialog

''warndlg" displays a dialog box named 'Warning Dialog' containing the string 'This
is the default warning string.' The warning dialog box will disappears after you
press the "OK" button.

3.7.5 Counter (Stopwatch Timer)

One of the main aspects upon completing this program is to determine what is the
program computation time and compare it to the existing approaches (Gauss Seidel,
Newton Raphson and Fast Decoupled). Hence, a stopwatch timer has been
implemented in this program to calculate how long it will take to complete the
calculation.

The counting is started by using command "tfc" and ended by command

"toe". The program is appended in between of the two, so that the total time elapsed
for the program to complete its computation and come to the final results is obtained.
The time will be displayed in the Command Window at the end of the calculation.
The syntax for these two functions is appended in the Appendix D.

29
 The speed of the computer computation is depends on the speed of the CPU
itself. If a fast computer is used, then computation time should be less. On the other
hand, if a slower computer is used, the increase in computation time is expected.
However, when we compares the computation time between these four methods
either by using a faster or slower computer, we should expect the ratio will be the
same for both computers.

3.8 Tools

These tools that are used throughoutthis project are listed below. They are consist if
hardware and software.

1) A PC equipped with at least Pentium III or equivalent processor, 128 MB

RAM, Microsoft 2000/ME/XP operating system
2) MATLAB software
3) MS Excel (Microsoft Office)

About MATLAB

MATLAB

.••l»iii<:ni "•! .MtMI'-MjII:',

ShowUcense OK

Figure 14: The MATXAB software

30
 CHAPTER 4

RESULTS AND DISCUSSION

The contingency analysis program has been tested on one test system, namely IEEE
30-bus Reliability Test Systems to test its practicality and reliability on real system
implementation. The IEEE 30-bus test case represents a portion of American Electric
Power System (AEP) in the Midwestern US as of December 1961. It is now being
made available to the electric utility industry as a standard case for evaluating
various analytical methods and computer programs for the solution of power system
problems [8], The one line diagram of the test case, its bus data and line data is
illustrated and appended in the Appendix B.

Before the program is run, the power flow analysis is conducted to get the
initial value of the bus voltages on each of the 30 buses. These voltages will be
stored in the MS Excel data file, and become the pre-outage bus voltages. Then, one
element of the line data is removed or in other words, it is open-circuited. Now,
instead of conducting power flow analysis again, the new bus voltages are predicted
using the contingency analysis program.

Then, the computation time for the program is determined and then compared
with the existing method available. In order to verify the results of the contingency
program, its output is compared with the output of the power flow analysis.

4.1 Results and Analysis

The program has been through series of testing of troubleshooting, to ensure the
reliability and the accuracy of the results. During one of the testing of the program
using IEEE 30 bus system, it is observed that the results are not so accurate,
compared to the power flow analysis program. At that time, the line data that are
being used are only inclusive impedance and reactance for each line in the power

31
system. However, it is noticed that they are not sufficient, and therefore the half line
charging admittance (HLCA) and the tap setting of the transformers are included in
the line data. The effects of these two added elements are discussed earlier in the

methodology section.

The Contingency Program will read the line data from MS Excel to construct
the bus admittance matrix for the networks. Once the matrix is constructed, it will
wait for further instruction. The new line data and the pre-outage bus voltages are
appended in the Table 2 and Table 3 below:

Element From To Bus Impedance Reactance HLCA Tap

No Bus No No Setting

1 1 2 0.0192 0.0575 0.0264 1

2 1 3 0.0452 0.1852 0.0204 1

3 2 4 0.0570 0.1737 0.0184 1

4 3 4 0.0132 0.0379 0.0042 1

5 2 5 0.0472 0.1983 0.0209 1

6 2 6 0.0581 0.1763 0.0187 1

7 4 6 0.0119 0.0414 0.0045 1

8 5 7 0.0460 0.1160 0.0102 1

9 6 7 0.0267 0.0820 0.0085 1

10 6 8 0.0120 0.0420 0.0045 1

11 6 9 0.0000 0.2080 0.0000 0.978

12 6 10 0.0000 0.5560 0.0000 0.969

13 9 11 0.0000 0.2080 0.0000 1

14 9 10 0.0000 0.1100 0.0000 1

32
15 4 12 0.0000 0.2560 0.0000 0.932

16 12 13 0.0000 0.1400 0.0000

17 12 14 0.1231 0.2559 0.0000

18 12 15 0.0662 0.1304 0.0000

19 12 16 0.0945 0.1987 0.0000

20 14 15 0.2210 0.1997 0.0000

21 16 17 0.0824 0.1923 0.0000

22 15 18 0.1073 0.2185 0.0000

23 18 18 0.0639 0.1292 0.0000

24 19 20 0.0340 0.0680 0.0000

25 10 20 0.0936 0.2090 0.0000

26 10 17 0.0324 0.0845 0.0000

27 10 21 0.0348 0.0749 0.0000

28 10 22 0.0727 0.1499 0.0000

29 21 22 0.0116 0.0236 0.0000

30 15 23 0.1000 0.2020 0.0000

31 22 24 0.1150 0.1790 0.0000

32 23 24 0.01320 0.2700 0.0000

33 24 25 0.1885 0.3292 0.0000

34 25 26 0.2544 0.3800 0.0000

35 25 27 0.1093 0.2087 0.0000

36 28 27 0.0000 0.3960 0.0000 0.968

33
37 27 28 0.2198 0.4153 0.0000

38 27 30 0.03202 0.6027 0.0000

39 29 30 0.2399 0.4533 0.0000

40 8 28 0.0636 0.2000 0.0214

41 6 28 0.0169 0.0599 0.0650

Table 2: The Line Data for IEEE 30 Bus Systems

Bus No Voltage magnitude Angle (°)

1 1.060 0.00

2 1.043 -5.497

3 1.022 -8.004

4 1.013 -9.662

5 1.010 -14.381

6 1.012 -11.398

7 1.003 -13.149

8 1.010 -12.115

9 1.051 -14.434

10 1.044 -16.024

11 1.082 -14.434

12 1.057 -15.303

13 1.071 -15.303

14 1.042 -16.198

15 1.038 -16.276

34
 16 1.045 -15.881

17 1.039 -16.188

18 1.028 -16.882

19 1.025 -17.051

20 1.029 -16.852

21 1.032 -16.468

22 1.033 -16.454

23 1.027 -16.661

24 1.022 -16.829

25 1.019 -16.423

26 1.001 -16.840

27 1.026 -15.912

28 1.011 -12.057

29 1.006 -17.136

30 0.995 -18.014

Table 3: The initial voltages at each bus after conducting the power flow analysis

>rd
4.1.1 Removal of 33™ element (Bus 24 - Bus 25)

When the 33rd element is removed, the new bus voltages are calculated using the
Contingency Program and tabulated in the Table 5 below. After the bus voltages are
known, the power flows in each line are also calculated and the results can be
observed from Table 6.

Then, the results are compared with the results from Power Flow Analysis.
As illustrated below, both tables show the comparison between the Contingency
Analysis and the Power Flow Analysis results.

35
 At the same time, the total time elapsed for the program to complete its
computation is recorded. The program is run for three times and the average time is
taken. Similarly, for the same contingency, three other programs for different
methods of power flow are run and the average computation time is taken. The
results are illustrated both in Table 4 and Figure 15.

Newton Fast Contingency

Method Gauss Seidcl
Raphson Decoupled Analysis

1 0,266 0,172 0,312 0,047

2 0,25 0,172 0,312 0,062

3 0,25 0,219 0,234 0,047

Average 0,255333333 0,187666667 0,286 0,052

Table 4: Comparison of the computation time (in seconds) between four (4) different methods of
power flow, including Contingency Analysis Program in table form

11

12

Average

Gauss Sei del Newton Fast Contingency

Ffephson Decoupled Analysis

Figure 15: Comparison of the computation time (in seconds) between four (4) different methods
of power flow, including Contingency Analysis Program whenever line 24 - 25 was removed.
Contingency Analysis program indicates remarkable results in term of its computation time
reduction.

36
 From the table and graph above, it is clearly seen that Contingency Analysis
gives remarkable reduced computation time compared to the existing methods for
Power Flow Analysis. This is really useful and beneficial, when dealing with a much
larger system, since it will takes less time to complete its computation.

Table 5: The new voltages at each bus when a line 33 is removed

Bus No Power Flow Analysis Contingency Analysis

1 1.0600 1.0596+ 0.0017i

2 1.0333-0.0982i 1.0378-0.0983i

3 1.0103-0.1424i 1.0111-0.1406i

4 0.9969-0.170H 0.9980-0.1683i

5 0.9781- 0.2517i 0.9779 - 0.2492i

6 0.9909 - 0.2002i 0.9917-0.1984i

7 0.9764 - 0.2287i 0.9768 - 0.2266i

8 0.9874 -0.21261 0.9871-0.2104i

9 1.0174-0.2635i 1.0173-0.260H

10 1.0037 -0.2905i 1.0033 -0.28621

11 1.0474-0.27131 1.0473-0.2678i

12 1.0194-0.2806i 1.0194 -0.277H

13 1.0326-0.28421 1.0325-0.2806i

37
14 1.0007-0.2926i 1.0006-0.2887i

15 0.9958 - 0.2932i 0.9957 - 0.2887i

16 1.0045-0.2877i 1.0043-0.2838i

17 0.9977-0.2918i 0.9974 - 0.28771

18 0.9833 - 0.3007i 0.9832 - 0.2964i

19 0.9799 - 0.3028i 0.9796 - 0.2985i

20 0.9848 - 0.3006i 0.9845 - 0.2963i

21 0.9897 - 0.2954i 0.9892 - 0.2904i

22 0.9904 - 0.2955i 0.9899 - 0.2903i

23 0.9839-0.2983i 0.9836 - 0.2922i

24 0.9778 -0.3016i 0.9774 - 0.293li

25 0.9711-0.2769i 0.9773 - 0.2887i

26 0.9520-0.279H 0.9583-0.2907i

27 0.9808 - 0.2747i 0.9864 - 0.28121

28 0.9872-0.2105i 0.9880 - 0.2097i

29 0.9556 - 0.2900i 0.9613-0.2968i

30 0.9400 - 0.3013i 0.9457 - 0.30731

38
 Table 6: Comparison between Power Flow Program and Contingency Program during the
evaluation of line flows on each line whenever line 24 - 25 was removed

POWER FLOW CONTIGENCY Deviation in (%)

ANALYSIS ANALYSIS

Using Gauss in P.U. System Real power Reactive

Seidel (P.U) power

r~ --2 1.7760-0.1314i 1.7775-0.2186i 0.15 -8.72

i 3 0.8345 + 0.0580i 0.8317+ 0.0526i -0.28 -0.54

2 4 0.4569+ 0.0062i 0.4569 + 0.028 li 0 2.19

4 0.7822-0.0257i 0.7803-0.0292i -0.19 -0.35

i

~5 0.8306-0.0085i 0.8294 + 0.0183i -0.12 2.68

ft 0.6177 -0.0319i 0.6189-0.0089i 0.12 2.3

4 ft 0.6982 -0.1843i 0.7011-0.1768i 0.29 0.75

5 -0.1419 + 0.11071 -0.1422 + 0.1043i -0.03 -0.64

~ft 0.3749-0.0243i 0.3748-0.0183i -0.01 0.6

"ft" « 0.2929-0.06321 0.2959-0.0364i 0.03 2.66

ft" 0.2822-0.078H 0.2756-0.07181 -0.66 0.63

"ft Id 0.1613+ 0.0037i 0.1574 + 0.007H -0.39 0.34

9 II 0.0000 -0.1569i 0.0003-0.15651 0.03 0.04

y 10 0.2823 + 0.0622i 0.2753 + 0.0680i -0.7 0.58

4 12 0.4458+ 0.140H 0.4398+ 0.1472i -0.6 0.71

12 P 0.0000 -0.1032i 0.0002 -0.1026i 0.02 0.06

12 ,14 0.0792 + 0.0232i 0.0786 + 0.0244i -0.06 0.12

12 15 0.1821+ 0.0653i 0.1773+ 0.0703i -0.48 0.5

39
12 Ift 0.0726+ 0.0323i 0.0719 + 0.0338i -0.07 0.15

14 i 15 0.0165+ 0.0056i 0.0154+ 0.0075i -0.11 0.19

ift r 0.0371 + 0.013H 0.0363+ 0.0145i -0.08 0.14

15 18 0.0597+ 0.0173i 0.0601+0.0174i 0.04 0.01

18 l<> 0.0274 + 0.0075i 0.0282 + 0.0084i 0.08 0.09

19 20 -0.0676-0.0266i -0.0670-0.0267i 0.06 -0.01

10 20 0.0907 + 0.0358i 0.0902+ 0.036H -0.05 0.03

10 1" 0.0533 + 0.0455i 0.0535 + 0.0438i 0.02 -0.17

10 21 0.1627+ 0.0922i 0.1560+ 0.0997i -0.67 0.75

10 22 0.0794+ 0.0407i 0.0750 + 0.0455i -0.44 0.48

2l". 22 -0.0131-0.0224i -0.0203-0.0160i -0.72 0.64

1* 23 0.0546 + 0.0242i 0.0492 + 0.03031 -0.54 0.61

0.0654+ 0.0174i 0.0546 + 0.0292i -1.08 1.18

2?; •24 0.0223+ 0.0075i 0.0167+ 0.0133i -0.56 0.58

24 25 0 0 0 0

25 26 0.0353 + 0.0237i 0.0354 + 0.0236i 0.01 -0.01

2^ 2" -0.0359-0.0234i -0.0512-0.0057i -1.53 1.77

0.1686+ 0.069H 0.1850 + 0.0529i 1.64 -1.62

r 29 0.0618+ 0.0167i 0.0627 + 0.0160i 0.09 -0.07

2" '30 0.0708+ 0.0167i 0.0708+ 0.0170i 0 0.03

2y 30 0.0369 + 0.0062i 0.0362+ 0.007H -0.07 0.09

8 28 -0.0082-0.0162i -0.0051 -0.02391 0.31 -0.77

ft 28 0.1770-0.089H 0.1910-0.0960i 1.4 -0.69

40
 The percentage of deviation from the Power Flow is calculated by finding the
different between both methods and multiplies it by 100%. This can be done simply
because all the units used have been converted to per unit system (P.U.) and
100MVA is used as the system base.

By looking to the deviation of the results of Contingency Program with

respect to the Power Flow Program, it can be seen that it does not exceed 10%,
which is a good sign. The maximum deviation is 1.38% for real power and -8.72%
for reactive power. On the hand, the program can achieve results as accurate as
0.01% for real power and 0% for reactive power. Hence, at this stage it can be
concluded that this program can give the closest results to the Power Flow Program.

4.2 Graphical User Interface

Programming for integrating Graphical User Interface with the MATLAB program is
complete but it does not stop there. It will be keep improving from time to time, as
the technology changes.

lEIgEJl

CONTINGENCY ANALYSIS
PROGRAM

MAIN MENU

LINE OUTAGE • OVERLOAD RELIEF »F

PROGRAM DETAILS

Contigency Analysisprogram is a program

that is able to predictthe steady state
bus voltages and tine currents after
each contingency is occured.

After all the NEW bus voltages and line

currents are obtained, power flows in each
line can be calculated easily- These power
will then be compared to the MAXIMUM
rating of the line. The MAIN purpose
is to determine whether the out-of-limit of
power flowexist.

Developedby Oppa. CopyrightReserved 20QB. fflhaflzjtuffiipnAeam

Figure 16: The main page of the program

41
The main page of the program consists of three sections, namely HOME, LINE
OUTAGE and OVERLOAD RELIEF. Each has its own functions. To start with

the analysis, click the LINE OUTAGE button. Then, the "CONTINGENCY
ANALYSIS" section will be displayed as illustrated in the Figure 17 below.

cm^"^""1^-
^MM

CONTINGENCY ANALYSIS
PROGRAM

lineuulaye 2
MAFN MENU

Devefeapedby Oppa. Copyright Resented 2O06.n*«fiziauJ@gn^£oin

Figure 17: The Single Contingency Section

Next, the data file is selected properly so that the program will read them.
The element number to be removed is entered in the box by clicking LOCATION
button. After clicking OK and CALCULATE button and the program will start
calculating the power flow in each line. Please refer Figure 18 and Figure 19.

42
 ^iSlxj

CONTINGENCY ANALYSIS
PROGRAM

„J=J*t
MAIN MENU

Developed by Qppe. Copyright Reserved 2006. [email protected]

Figure 18: Select the correct data

KfflSfti f% I '> i* i ^injxj

CONTINGENCY ANALYSIS
PROGRAM

MAIN MENU

Developed by Qppa. CopyrW^esened2DQB.nheftuadQgmel£(m

Figure 19: Enter the input and click CALCULATE!

43
 In order to see the results, click on the View Results! button, and another
window will appear. This is actually the one - line diagram of the IEEE 30 bus
system. Consequently, by clicking Update! button within the window, the power
flow on each line in the networks can be seen clearly. The displayed results can be
removed just by clicking the Clear! button.

mmm

Thirc Winding Transformer

1.7781-0.22251
Equivalents
13
i 9 IO
0.83231+0.0508741 C^ (-WW I
i / -2.2B47e-005-0.10273i
0.45725+0;027163i / 0.27574+0.0876131
I / 4i - v w t,12
sr
i WWH
/
/ 0.44012+0J488Bi
/
/ 14

7r3 0.07S75+0.024281iJ^'

0.83013+0.0185741
0.76067-0.(007611 13 / ^ 0.015387+0.00751721
I u \ 15
I 0.17756+0.070107i-

I 0.70138-0.17839J *Y2 \ O.Q6Q2tB+Qtf17285i

0.071965+0.0337141 16 j
I 0.049243+0.030383i 30
I
I 0.038301+0.014486i H~ll028277+T'*"
I0.61931-0.010215!II\
jr ^rt O.O063635U 8
U.UUUMU w -j-p
j j 23
*--•

QJB0377-tfLG3598i -0.067135-0.0256081 0.018677+0.013364*

JO, I ,-0.15754+0
~
'PM—I 111 19 H*= 0.070993+8.016474i
l« .00687871
-0.053675-K3.043654i
U~ /
0.037038+
^ I 0 0059273t

1 ™4*f»*:l»;075074+0.045505i I fj O.O620D7+0.018616J
sV ' ^ r • 22
\\ 9 _4.2743e-0Q5-0.158B6i -O-'oswe-oloo^oBBi 27 I
25 -H- oU 29

\\ 0.035444+0.0235451 0.t8483+0.052742i
_ 1 _2K
—:±J
0.29808-0.037042^ \ 0.19079-0.098497i . —-^

\ -0.8051563-0.023921!

w G: Generators
C: Synchronous condensers
Cbarl I Bdfl i

Figure 20: The displayed results

44
 CHAPTER 5

CONCLUSION AND RECOMMENDATION

An alternative way of doing power flow analysis when possible line outage occurs
has been presented in this project. Methods of contingency analysis which used the
factored Bus Admittance Matrix YbUS becomes attractive from computational
viewpoint, especially if the loads can be treated as constant current injections into the
various buses of the system. Even though the program does not give results as
accurate as power flow, it is not so much important because the system planners and
operators, who must undertake hundreds of studies in a short time period, are more
concerned with knowing if overload levels of current and out-of-limit voltages and
power flow exist than with finding the exact values of those quantities.

This Contingency Analysis program has been tested on an IEEE 30 buses

Test System to test its effectiveness, practicality and reliability. From the results and
findings, it is can be concluded that this program exhibits remarkable advantage in
term of its computation time reduction, as well as its accuracy to the exact solution.

Undoubtedly, it is a big advantage for conducting this analysis on a large

system. The computation time will be reduced significantly and therefore this will
benefits the electrical utility company itself.

It is therefore important to suggest some improvements and recommendations

for the benefit of the project. Since this area is growing, there are opportunities for
this approach to be widely applied in the power system engineering. The existing
system can be integrated with another system that control the generation of power, in
which the effects of both line outage and the power generation control, can be
exhibited by the program. If that can be done, this program will be a great success.

45
 REFERENCES

1. Computer aided contingency analysis [of power systems]

Carr, W.;
Rural Electric Power Conference, 1999
2-4 May 1999 Page(s):Cl/l - Cl/4

2. Fast automatic contingency analysis and ranking technique for power system
security assessment
Musirin, I.; Abdul Rahman, T.K.;
Research and Development, 2003. SCORED 2003. Proceedings. Student Conference
on 25-26 Aug. 2003

3. Comprehensive approach of power system contingency analysis

J. Deuse; K. Karoui; A. Bihain; J. Dubois
IEEE Bologna PowerTech Conference 2003
June 23-26, Bologna, Italy

4. Computer Methods in Power System Analysis, by Stagg, G.W., and EL-Ahmad,

A.H, McGraw-Hill Book Company, New York, 1968

5. Power System Analysis by Haadi Saadat, McGraw-Hill, Inc. Second Edition 2004.

6. Power System Analysis by J.J. Grainger and W.D. Stevenson, McGraw-Hill, Inc.
International Edition 1994

7. IEEE Explore
http://ieeexplore.ieee.org/Xplore/guesthome.isp

8. IEEE Power System Test Case Archive

http://www.ee.washington.edu/research/pstca/pf30/pg tca30bus.htm

46
 APPENDICES

LIST OF APPENDICES

APPENDIX A - MATLAB Code

APPENDIX B - IEEE 30 Bus One Line Diagram

APPENDIX C - IEEE 30 Bus Data

APPENDIX D - MATLAB Functions syntax

47
 APPENDIX A - The MATLAB Code

The MAIN program

function varargout = Contgency_Analysis(varargin)

% Begin initialization code - DO WOT EDIT

gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
1gui_OpeningFcn', @Contgency_Analysis_OpeningFcn, ...
'guijOutputFcn', @Contgency_Analysis_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', [] );
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{l});
end

if nargout
[varargout{1:nargout}] = gui_mainfen(gui_State, varargin{:});
else
gui_mainfcn{gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% Executes just before Contgency_Analysis is made visible.

function Contgency_Analysis_OpeningFcn(hObject, eventdata, handles,

varargin)

% This function has no output args, see OutputFcn.

% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to Contgency_Analysis (see VARARGIN)

% Choose default command line output for Contgency_Analysis

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes Contgency_Analysis wait for user response (see UIRESUME)

% uiwait(handles.figurel);

%Plot pictures on startup

scout = imreadf'gears640.jpg');
axes(handles.axesl);
image(scout);
axis off;

yaw_pedals = imread('qO.jpg');
axes(handles.axes3);
image(yaw_pedals);
axis off;

collective_lever = imread('qO.jpg');
axes(handles.axes4);
image(collective_lever);
axis off;

cyclic_stick = imread('qO.jpg');
axes(handles.axes5);
image(cyclic_stick);
axis off;

48
aa = imread('qO.jpg');
axes(handles.axes10);
image(aa);
axis off;

utp_logo = imread('ql.jpg');
axes(handles.axes7);
image (utp__logo) ;
axis off;

% Outputs from this function are returned to the command line.

function varargout = Contgency_Analysis_OutputFcn(hObject, eventdata,
handles)

% Get default command line output from handles structure

varargout{l} = handles.output;

%HOME Details
function tabl_ResizeFcn(hObject, eventdata, handles)
% hObject handle to uipanell (see GCBO)

%LINE OUTAGE Details

function tab2_ResizeFcn(hObject, eventdata, handles)
% hObject handle to uipanel2 (see GCBO)

%OVERLOAD RELIEF Details

function tab3_ResizeFcn(hObject, eventdata, handles)
% hObject handle to uipanel3 (see GCBO)

% Executes on button press in exit.

function exit_Callback(hObject, eventdata, handles}

close Contgency_Analysis;

% Executes on button press in HOME.

function scoutinfo_Callback(hObject, eventdata, handles)

handles = guidata(Contgency_Analysis);

set(handles.tabl,'Visible','on'];
set(handles.tab2,'Visible','off');
set(handles.tab3,'Visible','off');

% Executes on button press in LINE OUTAGE.

function controlinfo_Callback(hObject, eventdata, handles)

handles = guidata(Contgency_Analysis);

set(handles.tabl,'Visible','off);
set(handles.tab2,'Visible','on');
set(handles.tab3,'Visible','off);

% Executes on button press OVERLOAD RELIEF.

function program_Callback(hObject, eventdata, handles)

handles = guidata(Contgency_Analysis);

set (handles, tabl, 'Visible', 'off) ;

set (handles. tab2, 'Visible', 'off) ;
set{handles.tab3,'Visible','on');
% Executes on button press in SINGLE CONTINGENCY(under LINE OUTAGE),
function single_Callback(hObject, eventdata, handles)

lineoutage 2; %call top open a new window to analyze single contingency

49
% Executes on button press in VIEW DATA.
function data_Callback(hObject, eventdata, handles}

data; % call to open a new window with the one line diagram

% Executes on button press in ovrld_outage.

function ovrld_outage_Callback(hObject, eventdata, handles)

over_relief_lineoutage__2;

% Executes on button press in data2.

function data2jCallback(hObject, eventdata, handles)

data;

The Contingency Analysis program

function varargout = lineoutage_2(varargin)

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @lineoutage_2_0peningFcn, ...
1gui_OutputFcn', @lineoutage_2_0utputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin(l))
gui_State.gui_Callback = str2func(varargin(1}};
end

if nargout
[varargout{l:nargout}] = gui_mainfen(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% Executes just before lineoutage_2 is made visible.

function lineoutage_2_OpeningFcn(h0bject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to lineoutage_2 (see VARARGIN)

% Choose default command line output for lineoutage_2

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes lineoutage_2 wait for user response (see UIRESUME)

% uiwait(handles.figurel);

% Outputs from this function are returned to the command line.

function varargout = lineoutage_2_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

50
% Get default command line output from handles structure
varargout{1} = handles.output;

% Executes on button press in input_l.

function input_l_Callback(hObject, eventdata, handles)
% hObject handle to input_l (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

filename = uigetfile('*.xls'); %prompt to select the file name

if ~filename
msgbox('No file is selected. Please select the appropriate file','File
Not Found!','warn');
set(handles.filename,'String','No file selected')
return
else
[tmp,tmp,tmpl] = fileparts(filename);
if ~strcmp(lower(tmpl),'.xls')
msgbox('Inappropriate file selected. Please select ONLY THE EXCEL
file','Inappropriate File!','warn');
set(handles.filename, 'String', 'Wrong file type selected')
return
end
end

set(handles.filename,'String',filename}; ^display the file name

function filename_Callback(hObject, eventdata, handles)

% hObject handle to filename (see GCBO)

% Executes during object creation, after setting all properties.

function filename_CreateFcn(hObject, eventdata, handles)
% hObject handle to filename (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

if ispc
set(hObject,'BackgroundColor1,'white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor')) ;
end

% Executes on button press in read__data.

function read_data_Callback(hObject, eventdata, handles)
% hObject handle to read_data (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

filename = getappdata(0,'file')
if filename -= 1
msgbox('No file is selected. Please select the appropriate file','File
Not Found!','warn');
return

end

Edata = xlsread(filename, 'Edata' ); %read the line data from the

file
preout_V = xlsread(filename, 'preout_V'); %read the preoutage voltage from
the file

%the wait bar (to make it interesting)

h = waitbar(0,'Please wait...');
for i=l:1000, % computation here %
waitbar(i/100)
end
close(h)

51
%get the Ybus matrix,Line Data, number of element, number of bus
[Ybus,Y,Edata_new,Nele] = get_admittance_param(Edata);

%convert the voltage into real and imaginary

preout_V = read_preout_V_edited(preout_V);

fprintf('====================================================================

fprintf(['[ Element No |' ' Frm Bus 1' ' To Bus |' ' R j' '
X |* ' hlca |' ' tap |\n'])
fprintf (' =============^=============.========================:====™===

for k = l:Nele
fprintf(' \n'), fprintf(' I %2g', k), fprintf(' | %2g',
Edata_new(k,l) )
fprintf(' | %2g', Edata_new(k,2)), fprintf(' | %1.4f',
real(Edata_new(k,3)))
fprintf(' I %1.4f', imag(Edata_new(k,3))), fprintf(' | %1.4f',
imag(Edata_new(k, 4)))
fprintf(' | %1.4f', Edata_new(k, 5)), fprintf(' I')
end
fprintf (' \n=-======™==-========-====================================-=====

Edata_x = Edata_new;

^convert all variables to be used to handles

handles.Ybus - Ybus;
handles.Y = Y;
handles.Edata_new = Edata_x;
handles.preout_V = preout_V;
handles.filename = filename;
handles.Nele = Nele;
guidata(hObject,handles)

% Executes on button press in sing_cont.

function sing_cont_Callback(hObject, eventdata, handles)

set(handles.menu,'Visible','off);
set(handles.single,'Visible','on');
set(handles.double,'Visible','off');
set(handles.sing_cont,'Visible','off);
set(handles.multi_cont,'Visible','off);
cleardata;

% Executes on button press in multi_cont.

function multi_cont__Callback(hObject, eventdata, handles)

set(handles.menu,'Visible' , 'off' ) ;
set(handles.single,'Visible', 'off);
set(handles.double,'Visible','on');
set(handles.sing_cont,'Visible','off');
set(handles.multi_cont, 'Visible', 'off');
cleardata;

% Executes on button press in backl.

function backl_Callback(hObject, eventdata, handles)

set(handles.menu,'Visible','on');
set(handles.single,'Visible','off);
set(handles.double,'Visible','off');
set(handles.sing_cont,'Visible','on1);
set(handles.multi_cont,'Visible','on');
cleardata;

52
% Executes on button press in back2.
function back2_Callback(hObject, eventdata, handles}

set(handles.menu,'Visible','on');
set(handles.single,'Visible','off);
set(handles.double,'Visible','off);
set(handles.sing_cont,'Visible','on');
set(handles.multi_cont,'Visible','on');
cleardata;

% Executes on button press in location_l.

function location_l_Callback(hObject, eventdata, handles)

Nele = handles.Nele;

out = inputdlg('Which line to be remove?','Please Enter Input');

out = str2num(out{1,1});
setappdata(0,'outs',out);

out = getappdata(0,'outs');

if ((out == 0) I (out > Nele))

msgbox('The element is out of range. Please select the
appropriately','Element Not Found!','warn');
return
end

% Executes on button press in calcul__l.

function calcul_l__Callback (hObject, eventdata, handles)

Nele = handles.Nele;
Ybus = handles.Ybus;
Y = handles.Y;
Edata_x = handles.Edata_new;
preout_V = handles.preout_V;

out = getappdata(0,'outs');

%construction of the Ybus factored matrices

[L,U] = lu(Ybus);
j = diag(U);
L;
U = diag(j);

%reading the line which need to be removed (outage line)

outage = Edata_x(out,:); %read the entire line
m = outage (1); %pick up the 1st bus
n = outage(2); %pick up the 2nd bus
Za = outage(3); %pick up the impedance of the line

V = find_voltage(m,n,Y,L,U); %find the (m-n)th column of the Zbus matrix

Z = Za+V(m)-V(n); %find the total impedance

I = [preout_V(m)-preout_V(n)]/Z; %calculates the compensating current

deltaV = -V*I; %calculates the changes in voltages
Vbus = preout__V + deltaV %calculates new bus voltages

^converting the voltages to magnitudes and angles

for i = 1:30
x(i,:) = abs(Vbus(i)); %get the magnitude
y(i,:) = angle(Vbus(i))*180/pi; %get the angle
ind (i,:) = i; %generate the index (1,2,3,...)
end
Vangle = [ind x y] %display results

53
 CJP Tl (D P hh P c*P
I B I I ) J ) I D I D I B I 3 ( l ) ( D ( l ) f f l l D I I I ( l O ( D ( I ) ( I l ( I l l t U ( I ) ( D l l l l l ) ( l ) l l ( D I D ( D I I I ( I I O ( t l D ( l l l D I P ( t l l l l t to O 3 n, o M O
r t r t r t r t t t r t r t t t r t c t r t t t f t r t r t r t r t r t r t r t f t r t t t r t t t r t t t c t r t f t r t r t r t r t r t - r t r t r t r t r t r t rt x a 3 hi P- p
fupi(bpipiplUOiOiiliDi(liDJfu0ip]D>D]lDplDlD0'P'O'IlilliIliQ'IliQJDifT&iDifliQ)(b0i!liQJ o P- 3 h-<
H ii <d r+ (D p- ?«• (D O
^TiT5^^^^^^^^^^Tin3^^Ti^^Tin3^^^^^'0'0'0'0'0'0'0'D'OfOhO'0'0'0 (D p rt h-1 Hi C
aaaap.aaaaaao.aaaQ.aQ.aaao.aao.ao-ao.o-o-ao.aa.ao-ao.aa c/j o P co 11 II M
Q)[l;QJ{l]p^|ljQ;Q>|llpQj|llQ;[Lljl)Qf£l)QjQ)Q!llf(l)il]Qip]Q)P]^[lipjp]Qjp)P}p[}|[ll(U[lipi rt 3 (D ?>" P
r t t t d - f t r t r t r t r l r t r t t t r t r t t t r t r t r t f t t t r t t t r t r t r t r t r t r t r t r f r t t t r t r t r t r t r t r t f t f t r t r t p- ctP W h-t dp r t O P"u< hh a Q cn h-1 H1 rt
(D c, ~ a O P (D P1 ti — — II a> (C
— jr i-h P 13 O II ii ?r W II 3 3
ooooooooooooooooooooooooooooooooooooooooo P *• - P1 < P ll ^^ -—^
HiQ rt
3 ii O II P P3 Kl o P- rt P-
Co • • —-
C P h-> a a II II c 3 cr <d
Q)Q;[U|llQ)Qj|UQ;^[ll^Q]^|^|lj(llQ;|ll[[]pQ))l]QjQIQj|llp|l]p]Q)pQip>[llpp}QfQlfll{lipi S — r 1 p> Cd c II "^ p p rt fD ^—,

* i b U U W W U U U i i l U U M M ( O M M M M M M M P H P I - l P I - 1 l - ' l J l J l - l l C i a ) - J m ( J l i I . U J M P I l(I>
l II ^-^
p a (D Pi r t rt o o Cd 13
HowtD-Jmuiit'UiMFOBro-Joiiji^iiJMMowcc-JoiuiiiiuWHo - - - - - - - - - M II rt p a P P **. <*.
a o
•.»,-.'.i»s»,,'.,,,>,>»»>»>,,»,>,>,.l]Jti|Jlj]li|jigi(jlfJl(JI;J H- m <
S
c
p- r t
3 p
K
a
p
rt
I
>;
1
a
P
rt
s
(D
p
^XiVV^VVV^XiV^V^VX)XSX)XSXiX)^V'i3'T3X)'OrO*OX)XirO O O 0 O O O 0 O O 3 P D- co <o ! rt P ,— .—^ P hi
O00OO0000O00O0000000000000000000SSSSSSSSS tr c --^
N
P
1 pr *• t
ssssssssssssssssssssssssssssssss — CD u en a O — • X X hh
,—..—,,-,^,,-,,-,.-,,-.,-,,-,,^—.^..-,,-..-, toco -j cnui .c ui to H 1
P .—. .—. C ?r .—. K) M -^ h--
^.fcljlWWWUUIWWOJWMMMMMIOIOMMMPPPI-'MPMh'HP-'---'--"-'" 0 w a * hi - ?r ^-- ^--
•• 0
P" —.. -—- >< hj Ln
W
•*. *.. ^.
- s.
•r * n- CD — (-0 M
0 — n \ 3 '• ^— — P-
J^
3 — o rt rt ->»

CD P r 3
-
-.. LJ.
•d l-h CD
O ..—. > hi p
Hi to O 0
a 3 p-
rt hh I
j_j
p- ~—• ^-^ M
(D 1—'
< co P-
>-
tr r t op c*P CtP c*> t*> dP p
hh P •d rt TJ •a •c a (D
•1^ o w cr h- P" H- H- o CD
PJ -~~
c o <C n O s; rt
PJ dp hh co x w ?f fD CD
O <K> a -—
c P c c hi hi
s o P * rt n p. Tl •o 3
P- O h-^< •o dp 3 II H-
P 3 o rt rt •d H- r t rt P
ua •a d \ ho tr p- r t P" p- O CD
h^ p- rt 3 CD n r t (D tD **

<D P p a ?r P rt
,_,
X rt T! rt ri 3 M hh 3*
cd *—'
tr p •0 O P W O CD
•o .—•
d T3 (D a rt hi
!_,
O T3 + w
£ o - w rt O tr tr rt CD
CD S < CD P- hf, c c 3" 3
hi n> tr h- rt (D co co fD iQ
hi c 3 rt rt rt
,—. Co O P- P" O P"
3 hh p. a P" CD c
< t-1 C iQ P rt
a
> o t—'
a h-* P
Hi
— S P- — H-
p- P P P (0
M iQ — |-J (D
n h-1
H-
p I P-
3
>—i tr" 3
ID
>.
o (D
>
0
p P"
3 P
a H
ua
h3 H-
> p
fD iQ
% Executes on button press in viewdata_l.
function viewdata_l_Callback(hObject, eventdata, handles)

data;

% Executes on button press in location_2.

function location_2_Callback(hObject, eventdata, handles)

Nele = handles.Nele;

outl = inputdlg('Which line to be remove?','Please Enter Input');

outl = str2num(outlf1,1});
out2 = inputdlg('Which line to be remove?','Please Enter Input');
out2 = str2num(out2{l,l});
setappdatafO,'outl',outl);
setappdata(0,'out2',out2);

outl = getappdata(0,'outl');
out2 = getappdata(0, rout2');

if ((outl == 0) | (outl > Nele))|{(out2 == 0} I (out2 > Nele))

msgbox('The element is out of range. Please select the
appropriately','Element Not Found!','warn');
return
end

% Executes on button press in calcul_2.

function calcul_2_Callback(hObject, eventdata, handles)

Nele = handles.Nele;
Ybus = handles.Ybus;
Y = handles.Y;
Edata_x = handles.Edata_new;
preout_V = handles.preout_V;

outl = getappdata(0,'outl');
out2 = getappdata(0,'out2');

[L,U] = lu(Ybus);
j = diag(U);
L;
U = diag(j);

outage1 = Edata_x(out1, :) ;
m = outagel(1);
n = outagel (2);
z2 = outagel(3);

outage2 = Edata_x(out2,:);
p = outage2 (1);
q = outage2(2);
z3 - outage2(3);

V = find_voltage(m,n,Y,L,U);
VI = find_voltage(m,n,Y,L,U);

Z = [z2 0; 0 z3] ;
Zadd = [V(m)-V(n) Vl(m)-Vl(n); V(p)-V(q) VI(p)-VI(q)];;
Z = Z+Zadd;

Void = [preout_V(m)-preout_V(n); preout_V(p)-preout_V(q)];

I = inv(Z) * Void; %% Compensating currents, Ia and lb
deltaV - -[V VI]*I
Vbusnew = preout_V + deltaV;
X - Vbusnew;

55
for i = 1:30;
x(i,:)=abs(Vbusnew(i)),
y(i,:)=angle(Vbusnew(i[ *180/pi;
ind (i, :) = i;
end
oopa = [ind x y]

%Edata_x(out,:) = [];
nline = length(Edata_x(:,1));
for k = 1:nline
if (k =- outl)|(k == out2)
S(k) = 0;
J(k) = 0;
else
d = Edata_x(k,l);
f = Edata_x(k,2);
yt = 1/Edata_x(k,3);
hlca = Edata_x(k,4);
tap = Edata_x(k,5);
Idf = [(X(d)*yt/taprt2)-(X(f)*yt/tap)]+[X(d)*hlca; ^including HLCA
and TAP setting
S(k, :) = [X(d *conj(Idf)];
J(k,:) - [abs S(k) )];
number(k,1) = k;
end
end
power = [number S J];
pow = power(:,2);

setappdata (0, 'al',,pow(l)}j

setappdata (0, 'a2',,pow(2)};
setappdata (0, ra3',,pow(3));
setappdata (0, 'a4',,pow(4)};
setappdata (0, ra5',,pow(5)};
setappdata (0, ra6',,pow(6}) ;
setappdata (0/ 'aV,,pow(7));
setappdata (0/ 'a8'.,pow(8});
setappdata (0, 'a9\,pow(9));
setappdata (0, 'alO'',pow(10)
setappdata (0, 'all'',pow(11)
setappdata (0, 'al2'',pow(12)
setappdata (0- •al3'',pow(13)
setappdata (Of 'al4'T,pow(14)
setappdata (0, 'al5'',pow(15)
setappdata (0, 'al6'',pow(16)
setappdata (0, 'al7'',pow(17)
setappdata (0, 'al8'\pow(18)
setappdata (0, 'al9'',pow(19)
setappdata (0, 'a20'\ pow(20)
setappdata (0, 'a21'',pow(21)
setappdata (Of 'a22'',pow(22)
setappdata (0, 'a23'',pow(23)
setappdata (0, •a24'',pow(24)
setappdata (0, 'a25T',pow(25)
setappdata (0, 'a26'*,pow(26)
setappdata (0, 'a27'',pow(27)
setappdata (0, 'a28'r,pow(28)
setappdata (0, 'a29'',pow(29)
setappdata (0, 'a30'',pow(30)
setappdata (0, 'a31'',pow(31)
setappdata (0, 'a32'r,pow(32)

56
setappdatafO, 'a33',pow(33;
setappdata(0, 'a34',pow(34:
setappdata(0, ra35',pow(35:
setappdata(0, 'a36',pow(36;
setappdata(0, 'a37',pow(371
setappdata(0, 'a38',pow(38]
setappdata(0, 'a39',pow(39]
setappdata(0, 'a40',pow(40;
setappdata(0, 'a41',pow(4i;

% Executes on button press in viewdata_2.

function viewdata_2_Callback(hObject, eventdata, handles]

data;

% Executes on button press in exit.

function exit Callback(hObject, eventdata, handles)

Other functions

Calculate the difference of two columns

function V = diff_2column(p,q,Y_2,L,U)

for d=l:Y_2
Ibus(d) = 0;
end
Ibus(p) = 1;
Ibus(q) = -1;
Ibus = transpose(Ibus);

for s=l:Y_2
sum = 0;
for t = l:s-l
sum = L(s,t) * V(t) + sum;
end
V(s) = Ibus(s) - sum;
end

for s = 1:Y_2
V(s) = V(s)/U(s,s);
end

for s = Y_2:-l:l
d - s+1;
summ = 0;
if d > Y_2
summ = 0;
else
for t = d:Y_2
summ = L(t,s) * V(t) + summ;
end
end
V(s) = V(s) - summ;
end
V = transpose(V);

57
 APPENDIX B - IEEE 30 Bus One Line Diagram

6.9. DATA PREPARATION 225

Three Winding Transformer

Equivalents
13
0
HVVWH

G: Generators
C: Synchronous condensers

FIGURE 6.16
30-bus IEEE sample system.

58
 *zs
m
p
o u .^

Q) (0 CO , ,

-n o
£3
£
S3
O O O o o o o o en
tH
o o O O o o O O o CD o o o •
o o o o o o
<r
M

X
85 o o O o o o o o o <r • o sr o o o o O o o o o o o o o o o o o
M S m •r V CM CM
O c3
4J
83 rt o o O o o o o o o <o o o *o o o o o o o o o o o o o o o o o o
^^
H •H <j* <? •rt 1 1
1-1 aj s i i 1
CD a o
•=
ii o o o o o o o o o o o
O o u
a. 65 o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o
3 u K
>
rt •H o o o o o o o o o o o

u 3 o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o
u S •<*•
s aj
•H
U U
US o V- <M ID o o 01 o o o o IT) o *£} m CO CO Ol T r- CM o UJ ["- o CO o o at en

id TS o cm tH tH en o o o o CM o V- o •rt CM iH m o CO o tH o tH ID o CM o o o •rt
o CO tH rt tH CO iH
•H o
S-i -3 CM
aj 1 o o T \D CM o CO • o CO o >
CM CM m o CM to CM Lfl CM l> m * <a
1 |3 r- •rt o ,

1
•

I 1 s o • CM r- * o CM o o m o •rt (D CO CO m co CJl CM r- o CO CO o co o o CM o
m M co •rt tH
CM
H QJ
w HI cu
H -H u
(n & tn o o O o o o o o o o o
s- d m
Oi «s a o o O o o o o o o o o a o O o o o o o O o o O o o o o o o o

E- «Li
a. m ("ft
w 03 CO CM •rt
H 0 rt CO
a o o o o O o o o o o o o o
ca O <d
& !> R3 •rt •rt tH tH tH 1-1 tH tH T-t 1-1 •rt •rt •<-i •rt •rt tH rt •rt tH •rt iH tH tH •H tH rt tH tH tH •rt

i tb
o m •c
co pi o
m u tH CM o o CM o o CM o o CM a CM o o o o o o o o o O o o o o o o o
W
w «
ki ^ o
H P3 a tH CM co * m <£» r> CO en o •rt CM CO *r LO tD E^ CO CTl o iH CM CO *r m (D r~ CO en o
iH •rt tH T-t tH tH tH •rt tH •rt CM CM CM CM CM CM CM CM CM CM o
11
as
4J
115
73
tn
3
5.1' *v- s'l» J3
 APPENDIX D - MATLAB Functions syntax

uigetfile

Open standard dialog box for retrieving files

Syntax

uigetfile
uigetfile('FilterSpec')
uigetfile('FilterSpec','DialogTitle')
uigetfile('FilterSpec','DialogTitle','DefaultNarae')
uigetfile(...,'Location',[x y])
uigetfile(...,'MultiSelect',selectmode)
[FileName,PathName] = uigetfile(...)
[FileName,PathName,Filterlndex] = uigetfile(...)

\Is read

Read Microsoft Excel spreadsheet file (. xi s)

Syntax

nura = xlsread(filename)
nura = xlsreadffilename, -1)
nura = xlsreadffilename, sheet)
num = xlsread(filename, 'range')
num = xlsreadffilename, sheet, 'range')
num = xlsreadffilename, sheet, 'range', 'basic')
num = xlsreadffilename, ..., functionhandle)
[num, txt] = xlsreadffilename, ...)
[num, txt, raw] = xlsreadffilename, ...)
[num, txt, raw, X] = xlsreadffilename, ..., functionhandle)
xlsread filename sheet range basic

warndlg

Open warning dialog box

Syntax

warndlg
warndlg('warningstring')
warndlg('warningstring', 'dlgname')
h = warndlgf'warnstring','dlgname', createmode)
h = warndlg(...)

tic, toe

Measure performance using stopwatch timer

Synopsis

tic

any statements

toe

t = toe

60