Power System Dynamic Security Assignment via Prony analysis by qjz15341

VIEWS: 251 PAGES: 99

									      Power System Dynamic Security
       Assignment via Prony analysis
                                           By

                          Ek Wee. Tan (s40205388)
                                     B.Engg (Hons)




     A thesis submitted at the University of Queensland in partial fulfillment of the
                              requirement of the degree of



                        Bachelor of Engineering

   In the Department of Information Technology and Electrical Engineering (ITEE)
                                   October 2003

                           Supervisor: Dr Zhao Yang. Dong




Written by Tan Ek Wee                      -i-                                 s40205388
                                                           STATEMENT OF ORIGINALITY


20th October 2003


3/77 Warren Street
St Lucia
Queensland 4067
Australia


Professor Simon Kaplan
Head of School
School of Information Technology and Electrical Engineering
University of Queensland
St Lucia, Queensland 4072


Dear Sir,
As stated in the fulfillment of the requirement for the Bachelor of Electrical Engineering
Degree (Honors), The work presented here was performed under the supervision of Dr
Zhao Yang Dong. I hereby present the thesis entitle:


      “Power system Dynamic Stability via Prony Analysis”

I declare the research and work in this thesis submitted is my own, except as
acknowledged in the text, and has not been previously submitted for a degree in the
University of Queensland or any other institution.


Yours Sincerely


Ek Wee Tan




Written by Tan Ek Wee                      - ii -                              s40205388
                                                                   ACKNOWLEDGEMENTS


                             Acknowledgements
I would like to take this opportunity to thank the following people for their support and
guidance throughout the duration of my course.


Dr Zhao Yang. Dong for is advice and willingness to assist even outside the consultation
times.


Bertrem Koh, and Wei Jian Lim, for their patient and helpful guidance in the software
programming (MATLAB)


Specially deserving of my thanks is my loving wife, Janet Koh, for her constant support
and encouragements that made this thesis possible.




Written by Tan Ek Wee                     - iii -                             s40205388
                                                                                        ABSTRACT


                                       Abstract
This project describes one of many methods of system identification techniques, Prony
analysis and its application to power system stability assessment. Prony method provides
an efficient approach toward obtaining power system dynamic stability characteristics
from the network's response to disturbances. Prony analysis directly estimates the
strength, phase, damping and frequency of the system under study. The propose method
uses Prony analysis to identify the network function of the system and to decompose the
large system into a parallel combination of simple first order systems. With the network
identified, modal analysis method is used to determine the characteristic mode of the
system      model       linearized     about        a    specific     operating        point.


In this thesis, after power system stability fundamentals, Prony's method is presented,
then power system models are developed using power system block-set (PSB) in
MATLAB software for Prony analysis on stability characteristics. Observations were
made both before and after fault applied, results were than compared with the modes of
the system. The research outcome has shown that Prony analysis is robust and efficient
for power system stability assessment, especially in view of its modal analysis.




Written by Tan Ek Wee                      - iv -                                  s40205388
                                                                   TABLE OF CONTENT


                             Table of Contents
Statement of Originality
Acknowledgement
Abstract
Table Of Content                                                   Page no.
Chapter 1     Introduction                                            2
              Classification Of Transient                             4
              Transient simulation and Digital Simulation             6
              Useful Website Related To Transient Simulations         8
              Reference                                               10

Chapter 2     Transient Phenomena In Electrical System (General)      11
              Mathematical Approach In Linear Dynamic System          12
              Reference                                               15

Chapter 3     Power System Stability                                  16
              Steady State Stability                                  17
              Transient Stability                                     20
              Reference                                               24

Chapter 4     Modal Analysis                                          25
              Introduction To Modal Analysis                          26
              Application Of Modal Analysis To Power System           29
              Case Study                                              30
              Results Obtained                                        31

Chapter 5     System Identification Techniques                        37
              S-Domain                                                38
              Z-Domain                                                40
              Recursive Least Square Curve Fitting                    42
              Prony’s Method                                          45
              Reference                                               47

Chapter 6     Prony Analysis                                          48
              Prony Analysis Concept                                  49
              Prony Program                                           52
              Prony Vs Eigen Analysis                                 56
              Reference                                               58

Chapter 7     Case Study                                              59
              Power System Blockset (PSB)                             60
              Simulation                                              61
              Results                                                 61

Chapter 8     Conclusion and Future Considerations                    64
Appendix


Written by Tan Ek Wee                   -v-                                s40205388
Written by Tan Ek Wee   - vi -   s40205388
                                                              CHAPTER 1: INTRODUCTION


Chapter 1 Introduction


Introduction:
The operation of an electrical power system involves continuous electromechanical and
electromagnetic distribution of energy among the system components. During normal
operation, under constant load and topology, these energy exchanges are not modelled
explicitly and the system behaviors can be represented by voltage and current phasers in
frequency domain.


However, the switching events and system disturbance causes the energy exchanges
subject the circuit components to higher stresses, resulting from excessive current or
voltage variations, the prediction of which is the main objective of the power system
transient simulation.


Chapter 2: Describes different types of transients and the effects of each different
transient. Transient processes can be distinguished by their characteristic features, such as
the origin of the transient, the rate of change in the system, and the change in system
structure.




Written by Tan Ek Wee                       -1-                                  s40205388
Chapter 3: Power system stability, the ability to remain stable after the occurrence of
external or internal interferences. It briefly introduces the different types of stabilities and
the concept of maintains a stable power system. Mathematical approach to solve stability
problems of the power system, hence able to determine the critical power angle and the
amount of damp introduced when a transient occurs.


Chapter 4: In this chapter, it gives a brief introduction of modal analysis, in the case
where the system is determined by it power system oscillation behavior. A case study
were conducted using inter-area power system using matlab to further understand the use
of modal analysis and modal identifications, since Prony analysis is one of the modal
identification method.

Chapter 5: System identification techniques introduce the different kinds of
identification techniques that are currently used in identifying power system modes. It
gives a mathematical approach for each of the different approach.


Chapter 6: Prony analysis concept, it gives a detail mathematical explanation and tactic
of retrieving the exponential algorithm signal. A Prony program were written in matlab to
simulate input signals and returns a set of data which includes the amplitude, damping
factor, phase and frequency of the power system signal.


Chapter 7: A case study was conducted using Power System Blockset provided by the
software “MatLab” to simulate a single phase power system. The simulated power
system was set in a way to input a fault at bus 2 so as to create the transient liked signal.
Hence the signal was analysis using Prony analysis and the results are recorded in this
chapter.




Written by Tan Ek Wee                        -1-                                    s40205388
                                            CHAPTER 1: CLASSIFICATION OF TRANSIENT



   Power system control
                                                                                   Operator Actions

                                                                             LFC

                                                       Prime Control Mover

                                              Protection

                                           Generator Control

                                      HVDC, SVC, etc


   Power System phenomena
                                                                                      Daily
                                                                                      Load
                                                                                     Followin
                                                                      Tie-Line
                                                                     regulation

                                                           Long Term Dynamics

                                            Transient Stability

                                 Sub-synchronous

                          Switching

           Lighting

10-7            10-5           10-3             10-1              101              103                105


                                          1 cycle      1 second     1minute           1 hour 1 day
   Figure 1.1 Time frame of various transient phenomena




       Written by Tan Ek Wee                        -2-                                    s40205388
                                     CHAPTER 1: CLASSIFICATION OF TRANSIENT


Figure 1.1 shows the typical time frame of a full range of power system transient. The
transient on the left of the figure involves predominantly interaction between the
magnetic fields of the inductances and the electric fields of the capacitance in the system,
they are often referred to as electromagnetic transient, the transient on the right of the
figure are mainly involves by the interactions between the mechanical energy stored in
the rotating machine and the electrical energy stored in the network, they are often
referred as electromechanical transients. There is a grey area in the middle, namely the
transient stability region, where both effect plays a part and need adequate representation
[7]


Take for example, the lighting stroke produces the highest voltages surges and thus
determine the insulation levels, However at operating voltage of 400kV and above,
system generated over voltage, such as those caused by the energisation of transmission
lines, can often be the determining factor for insulation coordination.


The main aim of this thesis project is the use of efficient use of computation techniques to
the solution of the electromagnetic transient problems in the system of a small size
topology involving linear components. This is an essential part in power system design to
ensure satisfactory operation, derive the components ratings and optimize controller and
protection settings. It is also an important diagnostic tool to provide post-mortem
information with the following incidents.


Classification of electromagnetic transient
Transient waveforms contain one or more oscillatory components and can thus be
characterized by the natural frequencies of these oscillations. However the simulation of
the process, the accurate determination of the system oscillations is closely related to the
equivalent circuit used to represent the system components.


Hence from the modeling point of view, therefore it is more appropriate to classify
transient by the time range study, as shown in Figure 1.1.




Written by Tan Ek Wee                         -3-                               s40205388
         CHAPTER 1: TRANSIENT SIMULATION AND DIGITAL SIMULATION


Lighting, the fastest-acting disturbance, requires simulation in the region of nano top
micro-seconds, of course in this time frame, the variation if the power frequency voltage
and current levels will negligible and the electronic controllers will not respond, on the
other hand, the stray capacitance and inductance of the system components will exercise
the greatest influence in the response.


Switch of phenomena may require simulations on different time frames with
corresponding components models, i.e. either a fast transient model uses a stray
parameters or one based on simpler equivalent circuit but including the dynamics of the
power electronic controllers. In each case, the simulation step size will need to be at least
one tenth of the smallest time constant of the system represented.


Power systems are of two types, i.e. Those with essentially lumped parameters, such as
electrical machines and capacitor or reactors banks, and those with distributed
parameters, including overhead lines and underground or submarine cables. Following a
switching event, these circuit elements are subjected to voltage and current involving
frequencies between 50 Hz to 100kHz. Obviously within such a vast range of the values
of the components parameters and of the earth path will vary greatly with frequency. The
simulation process therefore must be capable of reproducing the frequency variation of
both the lumped and the distributed parameters. The simulation must also represent such
non-linearities as magnetic saturation, surge diverter characteristic and circuit breaker
arcs.


Transient simulation
In the past, transient simulations for the power system are used using electronic analogue
computer, the transient network analyzer (TNA) and the HVDC simulator.


The analogue computer basically solved ordinary differential equations by means of
several units designed to perform specific function, such as adder, multiplier and
integrator as well as signal generators and multi-channel cathode-ray oscilloscope [1]




Written by Tan Ek Wee                       -4-                                  s40205388
          CHAPTER 1: TRANSIENT SIMULATION AND DIGITAL SIMULATION


Greater versatility was achieved with the use of scaled down models and in particular the
TNA [1].     It is capable of emulating the behaviors of the actual power system
components using only low voltage and current levels. Early limitations include the use
of lumped parameters to represent transmission lines, unrealistic modeling of losses,
ground mode of transmission lines and magnetic non-linearities. However all these were
largely overcome [2] and TNAs are still in use for their advantage of operating in real
time, thus allowing many runs to be performed quickly and statistical data obtained, by
varying the instant of switching. The real-time nature of the TNA permits the connection
of the actual control hardware and its performance validated, prior to their
commissioning in the actual software associated with FACTS and HVDC transmission.
However due to their cost and maintenance requirements TNAs and HVDC models are
being gradually displaced by real time digital simulators.


Digital simulation
Owning to the complexity of the modern power systems, the simulators described above
could only be relied upon to solve simple problems. The advent of the digital computer
provides the simulation to the development of more accurate and general solutions.


While the electrical power system variables are continuous, digital simulation is by
nature discrete. The main task in digital simulation has therefore been the development of
suitable methods for the solution of the differential and algebraic equations at discrete
points.


The two broad classes of method used in the digital simulation of the differential
equations representing continuous systems are numerical integration and difference
equation. Although the numerical integration method does not produce an explicit
difference equation to be simulated, each step of the solution can be characterized by a
difference equation. [3]




Written by Tan Ek Wee                      -5-                                 s40205388
 CHAPTER 1: USEFUL WEBSITE RELATED TO TRANSIENT SIMULATIONS


  The EMTP (Electromagnetic Transient Program) is a computer algorithm for the efficient
  analysis of power system electromagnetic transient. It has become an industrial standard
  and many people have contributed to enhance its capability, with the rapid increase in
  size and complexity, maintenance and support become a major problem.
  With the development and availability of personal computer (PCs), programs are created
  to run simulations of power system transient. There is a selection of EMTP-type
  programs and their corresponding website are listed below


  EMTP (Electromagnetic Transient Program) – type programs


Program                 Organization                      Website address
EPRI/DCG EMTP           EPRI                              www.emtp96.com /
ATP program                                               www.empt.org /
MicroTran               MicroTran      Power     System www.microtran.com /
                        Analysis Corporation
PSCAD / EMTDC           Manitoba       HVDC      research www.hvdc.ca /
                        center
NETOMAC                 Simens                            www.ev.siemen.de/en/pages
                                                          /
NPLAN                   BCP Busarello + Cott + Partner
                        Inc
EMTAP                   EDSA                              www.edsa.com /
PowerFactory            DlgSILENT                         www.digsilent.de /
Arene                   Anhelco                           www.anhelco.com /
Hypersim                IREQ (Real – Time Simulator)      www.ireq.ca /
RTDS                    RTDS Technologies                 www.rtds.ca /
PowerSystemToolbox      Cherry Tree (MATLAB based)        www.eagle.ca/cherry /
Transient Performance MPR (MATLAB based)                  www.mpr.com /
Advisor (TPA)




  Written by Tan Ek Wee                    -6-                                 s40205388
  CHAPTER 1: USEFUL WEBSITE RELATED TO TRANSIENT SIMULATIONS


  Other transient simulation programs


Program                Organisation                       Website address
ATOSEC5                University of Quebec at Trios Cpee.uqtr.uquebec.ca/dcdtod
                       Riveres                            c/ato5_1htm
Xtrans                 Delft University of Technology     Eps.et.tudelft.nl
KREAN                  The Norwegain University of www.elkraft.ntnu.no/sie10aj/
                       Science and Technology             Krean1990.pdf /
Power System           MATHworks (MATLAB based)           www.mathworks.com/produ
                                                          ct/
Blockset               TransEnergie Technologies          www.transergie-tech.com/en/
SABER                  Avant (formerly Analogy Inc.)      www.analogy.com/
SIMSEM                 Swiss     Federal   Institute   of Simsen.epfl.ch/
                       Technology




  Written by Tan Ek Wee                     -7-                               s40205388
                                                           CHAPTER 1: REFERENCES


References
   1. PETERSON, H.A: “An electric circuit transient analyser”, General electric
      Review, 1939, p.394
   2. BORGONOVO,G., CAZZANI, M.,CLERICI, A., LUCCHINI, G. and VIDONI,
      G.: ”Fives years of experience with the new C.E.S.I.TNA”, IEEE Canadian
      Communication and power conference, Montreal,
   3. BICKFORD, J.P., MULLINEXUX, N. and REED J.R.: “Computation of power
      system transient” (IEE Monograph Series 18, Peter Peregrinus Ltd, London,
      1976)
   4. DERUSSO, P.M, ROY,R.J,CLOSE,C.M and DESROCHERS,A.A.: “State
      variables for Engineers”(John Wiley, New York, 2nd Edition)
   5. DOMMEL, H.W.: “Digital computer solution of electromagnetic transients in
      single and multi-phase networks”, IEEE Transactions on power Apparatus and
      systems, 1969,88(2), pp.734-71
   6. MOHAN,N.,ROBBINS, W.P., UNDELAND, T.M,NILSSEN,R. and MO,O.:
      Simulation of power electronic and motion control systems –an overview”,
      Proceeding of the IEEE, 1994, pp.1287-1302.
   7. NEVILLE,W. and JOS,.A. “Power systems electromagnetic transient simulation”,
      IEE Power and Energy series, 2003




Written by Tan Ek Wee                   -8-                            s40205388
           CHAPTER 2: TRANSIENT PHENOMENA IN ELECTRICAL SYSTEM


Chapter 2 Transient Phenomena in Electrical System (General)
Types of transient processes
A transient of a given system from one operating state (condition) to another state of
energy in the electrical or electromechanical circuit changes. Most importantly, this
complex phenomenon should be considered simultaneously in time and space, hence it is
a important factor that should be taken into account when solving a particular engineering
problem. The consideration of the factors, which are critical for a given process, makes it
possible to simplify the solution of the problem. Therefore the cases of transient
processes can be distinguished by their characteristic features as listed [1]:
    •   By their origin or the kinds of disturbing actions and the types of disturbance
    •   By assumptions made when deriving and solving differential equations
    •   By the rate of changes (the speed of the processes) in the system
    •   By the structure of the system


Quality of the transient process is the operating condition immediately following the
transient process. It must have a sufficient stability margin, the maximum change in the
parameters which still the system still continue to operate with stability determines the
value of this margin. it is important to have this “quality” (the values of the voltage and
frequency, degree of voltage curve distortion), which can be use to estimate the new
operating condition. Another condition that relates directly to the transient behavior of
the system, the time required to return to it steady-state condition. It can be aperiodic,
oscillatory, monotonic or intermediate between oscillatory or monotonic. Besides
estimating the transient behavior in an electrical power system, it is more important to
know the effect of the transient on the operating condition of the entire system. For
example, it is always required that during a “well behave” transient process in the voltage
V(t) will not drop down to the values dangerous from the point of view if the stability of
the generators and motors of the entire system. However the voltage deviation during the
transient process cannot be estimated only by the instantaneous change in its absolute
value. It is essential to know how long this deviation existed. Therefore, determining the
standard deviation from the voltage, which is optimal for the given part of the system,
often makes such estimation.


Written by Tan Ek Wee                        -9-                                 s40205388
CHAPTER 2: MATHEMATICAL APPROACH IN LINEAR DYNAMIC SYSTEM


Block diagram of the system




Mathematical approach in linear dynamic system
The equivalent circuits of the system element relate the process variable to the parameters
of the element occurs. The above block diagram shows the connection between the
elements in the system and shows the changes undergo by the operating variables in these
element, where the output of each element y(t) is expressed as a function of the variable
at its input x(t).it is characterized by the ratio of the output to the input, also known as the
transfer function, TR = y(t) / x(t).
Linear algebra and circuit theory and concepts are used in power system analysis, as so to
describe the formulation of state equations and linear dynamic systems.
An electrical Power System is basically a continuous system. And with the advance
technology in the modern age, Digital simulation and many other computations are
simple if the correct analytical tools were used. Digital simulation is by nature a discrete
time process and can only provide solution for the differential and algebraic equation at
discrete points in time.




Written by Tan Ek Wee                        - 10 -                                 s40205388
CHAPTER 2: MATHEMATICAL APPROACH IN LINEAR DYNAMIC SYSTEM


Continuous system
An nth order linear dynamic system is described by an nth order linear differential
equation, which can be written as n first order linear differential equations [1]


       x1 ( t ) = a11 x1 ( t ) + a11 x2 ( t ) + K + a1n xn ( t ) + b11u1 ( t ) + b12u2 ( t ) + K + b1mum ( t )
       &
       x2 ( t ) = a21 x1 ( t ) + a22 x2 ( t ) + K + a2 n xn ( t ) + b21u1 ( t ) + b22u2 ( t ) + K + b2 mum ( t )
       &
       M
       xn ( t ) = an1 x1 ( t ) + an 2 x2 ( t ) + K + ann xn ( t ) + bn1u1 ( t ) + bn 2u2 ( t ) + K + bnmum ( t )
       &
                                                                                                                2.1
Expressing the above equation in matrix form
                   &
                   x1        a11     a12    L      a1n      x1         b11   b12 L        b1m      u1
                   &
                   x2        a21     a22 L         a2 n     x2         b21   b22 L        b2 m     u2
                         =                                         +
                   M           M      M O           M       M           M     M O           M       M
                   &
                   xn        an1     an 2 L        ann      xn         bn1 bn 2 L         bnm      um
                                                                                                                     2.2
Or in short notation
                                                  x = [ A] x + [ B ] u
                                                  &
                                                                                                                     2.3
Which is also referring to as state equation.
Since there is an equation for the input of the power system, there should be another
equation for the system output state vector, which relates to the input vector
      y1 ( t ) = c11 x1 ( t ) + c11 x2 ( t ) + K + c1n xn ( t ) + d11u1 ( t ) + d12u2 ( t ) + K + d1mum ( t )
      y2 ( t ) = c21 x1 ( t ) + c22 x2 ( t ) + K + c2 n xn ( t ) + d 21u1 ( t ) + d 22u2 ( t ) + K + d 2 mum ( t )
      M
      yn ( t ) = cO1 x1 ( t ) + cO 2 x2 ( t ) + K + cOn xn ( t ) + dO1u1 ( t ) + d O 2u2 ( t ) + K + d Omum ( t )
                                                                                                               2.4
Writing it into matrix form again we have,
                  y1   c11 c12             L      c1n      x1   d11 d12             L      d1m      u1
                  y2   c   c22             L      c2 n     x2   d    d 22           L      d2m      u2
                     = 21                                     + 21
                  M     M M                O       M       M     M M                O       M        M
                  yn   cO1 cO 2            L      cOn      xn   d O1 d O 2          L      dOm      um
                                                                                                                     2.5
Or in compact equation


Written by Tan Ek Wee                                     - 11 -                                             s40205388
CHAPTER 2: MATHEMATICAL APPROACH IN LINEAR DYNAMIC SYSTEM


                                         y = [C ] x + [ D ] u
                                                                                    2.6
Which is also call the output equation
Therefore the equation in 2.3 and 2.6 can be solved by using transformation methods, the
convolution or numerically in an iterative procedure. The form of state variable equations
is not unique and mainly depends on the choice of state variable [2].


State variable formulation can be expressed by the following transfer function
                               a0 + a1s + a2 s 2 + a3 s 3 + L + aN s N Y ( s )
                     G( s) =                                            =
                                b0 + b1s + b2 s 2 + b3 s 3 + L + bn s n   U (s)
                                                                                    2.7
Where n ≥ N
By dividing the numerator and denominator by bn provides the standard form, such that
the s term appears in the denominator with unit coefficient, the rational function in s in
time domain
                               a0 + a1s + a2 s 2 + a3 s 3 + L + aN s N Y ( s )
                     G( s) =                                           =
                                1 + b1s + b2 s 2 + b3 s 3 + L + bn s n   U (s)
                                                                                    2.8




Written by Tan Ek Wee                          - 12 -                             s40205388
                                                            CHAPTER 2: REFERENCES


References
   1. Venikov, V.A. “Transient Processes in Electrical Power System”, 1977 pp19 – 30


   2. Marple.S, Lawrence.Jr “Digital spectral Analysis” Prentice-Hall, 1987




Written by Tan Ek Wee                   - 13 -                                s40205388
                                               CHAPTER 3: POWER SYSTEM STABILITY


Chapter 3 Power System Stability


Power System Stability
There is a need for the research of the behavior of an electric power system when subject
to small disturbances. It is assume the system under study has been disturbed from the
steady-state condition. This small disturbance may be temporary or permanent. If the
system is stable, it would than be expected to be temporary disturbance to the system then
back to the initial condition, while permanent disturbance would cause the system to
acquire a new operating state after a transient period. In either case synchronism should
not be lost.
Stability problems are generally divided into two major categories, the steady-state
stability and transient stability. Steady-state stability refers to the ability for the power
system to regain synchronism after a small disturbance. Transient stability deals with the
outage of a line or the sudden application or removal of loads. Therefore transient
stability studies are required to ensure that the system can withstand the transient
condition generating and transmitting facilities are planned. Hence they can be classified
as the followings:
    •   Transient Stability
        Ability of the power system to maintain synchronism during and immediately
        following a major disturbance such as a transmission line fault or the loss of a
        large generating unit


    •   Steady-state Stability
        Ability of the power system to maintain synchronism at all points for incremental
        slow-moving changes in power outputs of units or power transfer over
        transmission facilities


    •   Small Signal (Oscillatory Stability)
        Ability of the power system to maintain synchronism during small changes in
        operating conditions, which produce small changes in generator angle, speed and
        power


Written by Tan Ek Wee                      - 14 -                                s40205388
                                                CHAPTER 3: STEADY-STATE STABILITY


Steady-state Stability
Refers to the ability of the power system to remain in synchronism when subjected to
small disturbances. The motion is free and the stability is assured if the system return to
its original state, such behavior can be determined in a linear system by examining the
characteristic equation of the system, assume that automatic controls, like the governor
and voltage regulator are not active.


To illustrate the stability problem, we will consider the behavior of a one-machine system
connected to an infinite bus.




Swing equation
Swing equation describes the relative motion of the machine, i.e. the relative position of
the rotor axis and the resultant magnetic field axis is fixed, and the angle between the two
is known as the power angle or torque angle. The equation is given by
                                   H d 2δ
                                             = Pm − Pmax sin δ
                                  π f 0 dt 2
                                                                                     3.1
the swing equation is a non-linear function of a power angle, however the swing equation
may be linearized with a little loss of accuracy if we consider a small deviation ∆δ in the
power angle from the initial operation point δ 0 . i.e. δ = δ 0 + ∆δ


Written by Tan Ek Wee                        - 15 -                              s40205388
                                               CHAPTER 3: STEADY-STATE STABILITY


As long as there is a difference in angular velocity between the rotor and the resultant air
gap field, induction motor action will take place between them, and the torque will be set
up on the rotor tending to minimize the difference between the two angular velocities,
this is known as the damping torque, the damping power is described as:
                                                  dδ
                                         Pd = D
                                                  dt
                                                                                          3.2
When the synchronizing power coefficient Ps is positive, because of the damping power,
oscillations will damp out eventually, and the operation at the equilibrium angle will be
restored. No loss of synchronism occurs and the system is stable. If the damping is
accounted for, the swing equation becomes:
                               H d 2δ       d ∆δ
                                         +D      + Ps ∆δ = 0
                              π f 0 dt 2
                                             dt
                                                                                          3.3
In terms of second order differential equation, we have
                              d 2 ∆δ         d ∆δ
                                   2
                                     + 2ζω n      + ω n 2 ∆δ = 0
                               dt             dt
                                                                                          3.4
Where ω n is the natural frequency of the oscillation and ζ is define as the damping ratio,

                                               π f0
for normal operating conditions, ζ = D / 2             < 1 and the roots of the characteristic
                                                HPs

equation are complex

                                s1 , s2 = −ζω n ± jω n 1 − ζ 2

                                    s1 , s2 = −ζω n ± jω d
                                                                                          3.5
It is clear that for positive damping, roots of the characteristic equation have a negative
real part if synchronizing power coefficient Ps is positive. The response is bounded and
the system is stable.[1]




Written by Tan Ek Wee                       - 16 -                                 s40205388
                                             CHAPTER 3: STEADY-STATE STABILITY




Ref [2]
The above figure shows the effect of the machine at different conditions; note that
automatic controls, like the governor and voltage regulator are not active.




Written by Tan Ek Wee                      - 17 -                              s40205388
                                                      CHAPTER 3: TRANSIENT STABILITY


Transient Stability
The transient stability is the determination of whether or not synchronism is maintain
after the machine has been subjected to severe disturbance. This might happen when
there is a sudden increase in load, loss of generation, loss of large load or a fault in the
system. In most cases, oscillations are of such magnitude that linearization is not
permissible and the non-linear equation must be solved.
Therefore a method was introduced which is known as equal-area criterion, which can be
used for a quick prediction of stability. The method is only applicable to a one-machine
system connected to an infinite bus or a two-machine system.




ref[2]
The equation derived from the swing equation, with the consideration of the damping
neglected
         2             δ
 dδ              2π f 0
             =          ( Pm − Pe ) dδ
 dt               H δ0

or

dδ   2π f 0 δ
   =          ( Pm − Pe ) dδ
dt    H δ0

                                                                                               3.6



Written by Tan Ek Wee                        - 18 -                                 s40205388
                                                         CHAPTER 3: TRANSIENT STABILITY


The equation 3.6 gives the relative speed of the machine with respect to the
synchronously revolving reference frame. For stability, this speed must become zero at
some time after the disturbance, therefore we have the stability criterion as,
δ
     ( Pm − Pe ) dδ   =0
δ0

                                                                                          3.7
Consider the machine operating at the equilibrium point δ o , corresponding to the

mechanical power input Pm 0 = Pe 0 as shown in the figure below. Consider a sudden step

increase in input power represented by the horizontal line Pm1 . Since Pm1 > Pe 0 , the

accelerating power on the rotor is positive and the power angle δ increases. The excess
energy stored in the rotor during the initial acceleration is


δ1
     ( Pm1 − Pe )dδ = area abc = area A1
δ0

                                                                                          3.8
With increase in δ , the electrical power increases, and when δ = δ 1 , the electrical power

matches the new input power Pm1 . Even through the accelerating power is zero at this

point, the rotor is running above synchronous speed; hence, δ and electrical power Pe

will continue to increase. Now Pm < Pe , causing the rotor to decelerate toward

synchronous speed until δ = δ max . According to 3.7, the rotor must swing past point b
until an equal amount of energy is given up by the rotating masses. The energy given up
by the rotor as it decelerates back to synchronous speed is


δ max
         ( Pm1 − Pe )d δ = area bde = area A2
    δ1

                                                                                          3.9




Written by Tan Ek Wee                           - 19 -                            s40205388
                                                             CHAPTER 3: TRANSIENT STABILITY


The result is that the rotor swings to point b and the angle δ max , at which point


area A1 = area A 2

                                                                                           3.10
This is known as the equal-area criterion. The rotor angle would then oscillate back and
forth between δ 0 and δ max at its natural frequency. The damping present in the machine
will cause these oscillations to subside and the new steady state operation would be
established at point b.


                                   Equal-area Criterion
         Pe


                                        d

                                   A2
                    c     b
         Pm1                                e
                     A1


         Pm0        a

                                                                               δ
                   δ0         δ1            δ max                  π


Application to three-phase fault
Consider the fault at some distance between the sending end and the receiving end. With
the fault location located far away from the sending end, the equivalent transfer reactance
between bus bar increased, lowing the power transfer capability, where finally a post
fault power angle curve were formed, assuming the faulted line is removed. When a fault
occurs, the operating point shifts immediately. An excess of mechanical input over its
electrical output accelerates the rotor, and thereby storing kinetic energy, and the angle
δ . Assume that the fault is cleared at time δ1 by isolating the faulted line. The net power
will decrease and previously stored kinetic energy will be reduce to zero. Since Pe is still
greater than of Pm, the rotor will continue to decelerate, the rotor angle will then oscillate


Written by Tan Ek Wee                               - 20 -                            s40205388
                                                     CHAPTER 3: TRANSIENT STABILITY


back and forth at its natural frequency. The damping in the machine will cause this
oscillation to subside and a new operation point will be established. [2]




Ref [2]
The application of equal-area criteria gives the critical angle to maintain stability, hence
at all time, it is a must to keep the machine running within the boundary limit.




Written by Tan Ek Wee                       - 21 -                                 s40205388
                                        CHAPTER 3: POWER SYSTEM STABILITY


References
[1]   Hadi Saadat “Power System Analysis” Series in Electrical and Computer
      Engineering, chapter 11, pg 460-256


[2]   Les Hajagos “Power system stability Concept”, Kestrel Power Engineering Ltd.




Written by Tan Ek Wee                  - 22 -                             s40205388
                                                          CHAPTER 4: MODAL ANALYSIS


Chapter 4 Modal Analysis
Introduction
To really understand the stability of the power system, one must understand the power
system oscillation, which generally requires a combination of analytical tools.
Oscillations are often observed in transient stability simulations performed for the system
planning or other operations studies. In most cases, operators will pay close attention to
disturbances and when system oscillation stability limit has been exceeded. Hence, in
such cases, specialized analytical tools are necessary so that the nature of the oscillation
could be understood, so system’s control and operating procedures can be modified to
ensure system security. It is also important to make use of all analytical tools to
understand all aspects of the system dynamic performance.


Specialized tools for the study of power system oscillations can be categorized into two
basic groups, namely
       •   Modal Analysis
       •   Modal Identification
Where, modal analysis involves in determination of the characteristic modes of the system
model linearized about a specific operating point. And Modal Identification involves the
determination of characteristic modes from the dynamic behavior obtain either from
system measurement, or from transient stability simulations using a non-linear model. [ref:
Analysis and control of power system oscillation]


In this chapter, we will discuss a fair bit of modal analysis, and simulation of a modal
analysis software program using “MatLab”. The simulation uses a simple 4 machines, 2-
area system. Modes will then be compared with the application Power stability system
(PSS), and power system without the use of PSS.




Written by Tan Ek Wee                      - 23 -                               s40205388
                              CHAPTER 4: INTRODUCTION TO MODAL ANALYSIS


Basic Modal analysis and formulation
In the modern days, as power has become more and more important to the society, it is
important that simulation of power system has equations to represent the transmission
network and power systems dynamic devices, such as generators and their controls. The
equations are in the general form: to be made into a routine procedure for the system
planning and operation. Most simulators use non-linear algebraic-differential
                                       dy
                                          = f ( x, v, d )
                                       dt
                                         0 = g ( x, v )
                                         y = h ( x, v )
                                                                                     4-1
where
x       is a vector of dynamic variables, referred to as “states”
v       is a vector network voltage
d       is a vector of disturbances or control inputs
y       is a vector of monitored outputs for control or information
f       represents the nonlinear dynamic characteristic of the system dynamic
        components
g       represents the nonlinear network equations
h       represents the nonlinear output equation


In the modal analysis, the nonlinear equations are linearized about a operating point by
truncating a Taylor series expansion of f and g to the first term. Therefore at steady state,
the operating point is defined by x0 and v0, then f(x0 , v0) and g(x0 , v0) are both zero.
Then the linearized equations are:
                               d ∆x ∂f       ∂f       ∂f
                                   = ∆x + ∆v +           ∆d
                                dt   ∂x      ∂v       ∂d
                                        ∂g      ∂g
                                    0=     ∆x +     ∆v
                                        ∂x      ∂v
                                         ∂h      ∂h
                                   ∆y = ∆x + ∆v
                                         ∂x      ∂v
                                                                                     4-2
Or




Written by Tan Ek Wee                       - 24 -                               s40205388
                             CHAPTER 4: INTRODUCTION TO MODAL ANALYSIS


                                 d ∆x
                                      = Ad ∆x + bdv ∆v + bdd ∆d
                                  dt
                                       0 = Cd ∆x + Yn ∆v
                                      ∆y = C0 ∆x + K n ∆v
                                                                                    4-3
Of the above-linearized equations, such as this are the starting points of all power system
modal analysis programs. But the way they present usually depends on the size of the
power system model being analyzed. Since the equation includes governing the system’s
dynamic devices and static transmission network equations, equation 3-3 is often referred
to as being in augmented form.


The system state equations
The state equation of the system is given by the following where the algebraic equations
are eliminated,
                                     d ∆x
                                          = A∆x + Bdd ∆d
                                      dt
                                                                                    4-4
With

                                    A = Ad − Bdv (Yn ) Cd
                                                      −1



                                    C = C0 − K 0 ( Yn ) Cd
                                                      −1


                                                                            4-5
Equation 3-3 produces an sparse structure (because the device dynamic equations are
block diagonal and the network matrix is by nature sparse), once it has been reduced to
state level equation, sparsity were greatly reduced, this lack of sparsity was one of the
factor that results the restricted model analysis programs to relative small power system
models.
                                                                      λ it
The characteristic mode of the equation 3-4 has a general form of ui e . The number of
characteristic modes is equal to the amount of states in a given power system model. The
vector ui is called a characteristic vector (or eigen vector) and the value λi is called the
characteristic value (or Eigenvalue). Real Eigenvalues indicate modes, which are
aperiodic. Complex Eigenvalue indicate modes, which are oscillatory. Let say for a




Written by Tan Ek Wee                        - 25 -                             s40205388
                             CHAPTER 4: INTRODUCTION TO MODAL ANALYSIS


complex eigenvalue α ± jω , the amplitude of the mode varies with as e and frequency
                                                                      αt




of the oscillation will be ω / 2π . The damping ratio
                                                      (ζ ) can be define as:
                                                 −α
                                      (ζ ) =
                                               α2 +ω2
The Eigenvalues satisfy:
                                      det( A − λi I ) = 0
                                                                                   4-6
The right eigenvectors are column which satisfy:
                                         Aui = λi ui
                                                                                   4-7
The left eigenvectors are row vectors which satisfy:
                                         vi A = λi ví
                                                                                    4-8
The left and right eigenvectors are orthogonal, and are usually scaled to be orthonormal.
The state equation given by 3-4 can be expressed in terms of modal variables by using
the transformation x = Uz , which leads to:
                            dz
                               = VAUz + VBdd ∆d = Λz + VBdd ∆d
                            dt
                            y = CUz
                                                                                   4-9
The matrix Λ is a diagonal matrix of the Eigenvalues.




Written by Tan Ek Wee                       - 26 -                             s40205388
    CHAPTER 4: USING MODAL ANALYSIS FOR SMALL POWER SYSTEM


Using Modal analysis for small power system.
Eigen analysis is being discussed in the previous section, which normally referred to as
QR method 3-1. The Q is straightforward to calculate, because of its structure. From
these, the eigenvectors of the original matrix can be obtained by reversing the
transformations.
Practical computer limitation has restrict the size of the state matrix to which full modal
analysis using the QR method. Hence in general, eigen analysis restrict the size of the
system for QR Eigenvalue calculation up to about 500 states, this implies that only 50
generators may be modeled since about 10 states are required for a generator model and
together with its controls. Currently, the power systems available uses model of at least
400 states.
However if the states level falls with the calculation limit, it is best that it uses QR
calculation to analysis the full power system model, since all the system Eigenvalues and
vectors are known. Also with this type of analysis, it is easy to provide the performance
measures of participation and residues.
As the power system grew, and the modeling requirements also expand, hence more
analysis techniques were required to continue with the modal analysis of the power
systems. Prony analysis is one of the identification techniques that will be discussed later
chapters.




Written by Tan Ek Wee                      - 27 -                               s40205388
                                                                 CHAPTER 4: CASE STUDY


Case Study
The single-line diagram of the power system is shown in the following figure (4-1)
The case study focuses on the inter-area system that is supported by 4 generators. SVC is
fitted with the controller to regulate their respective bus voltages. With the use of the
MATLAB software, modal analysis were used to analysis the following inter-area power
system. The results will compare the modes of each system prior to and after installation
of PSS system on each generator.




                  Figure 4-1: 2 Area inter-area power system


The simulation uses Power System toolbox program that is installed within the MATLAB
application program. svm_mgen.m is a driver for small signal stability analysis in the
Power System Toolbox, it requires a input data file as to perform a load flow analysis. It
forms the state matrix of the power system model, linearized about an operating point set
by a load flow and finally performs a modal analysis of the power system. The source
code is attach to the appendix, together with the data of the inter-area power system
model, data2.m.


Using the svm_mgen.m driver to run a modal analysis on the two-area test case (fig4-1),
data2.m obtains the following results.



Written by Tan Ek Wee                      - 28 -                               s40205388
                                                                 CHAPTER 4: CASE STUDY



The total number of states within the system in two-area test case without the installation
of Power Stability System (PSS) :
States =       13 (Generator 1)
               12 (Generator 2)
               10 (Generator 3)
               13 (Generator 4)
               1 (SVC)
Total states = 49




Written by Tan Ek Wee                      - 29 -                               s40205388
                                                              CHAPTER 4: CASE STUDY


The followings are the readings of the Eigenvalues calculated, Real Eigenvalues indicate
modes that are aperiodic. Complex Eigenvalue indicate modes that are oscillatory. As
you can see that there are 10 pairs of complex values, each representing the oscillatory
modes.
l=
 1.0e+002 *
         -0.0000
         -0.0011 - 0.0007i
         -0.0011 + 0.0007i
         -0.0020
         -0.0020
         -0.0020
         -0.0026 - 0.0064i
         -0.0026 + 0.0064i
         -0.0032 - 0.0063i
         -0.0032 + 0.0063i
         -0.0152 - 0.0095i
         -0.0152 + 0.0095i
         -0.0198
         -0.0199
         -0.0200
         -0.0200
         -0.0251
         -0.0351
          0.0012 - 0.0360i
          0.0012 + 0.0360i
         -0.0432
         -0.0498
         -0.0056 - 0.0687i
         -0.0056 + 0.0687i
         -0.0045 - 0.0716i
         -0.0045 + 0.0716i
         -0.1000
         -0.1000
         -0.1000
         -0.1000
         -0.0863 - 0.0990i
         -0.0863 + 0.0990i
         -0.1925 - 0.1027i
         -0.1925 + 0.1027i
         -0.2834
         -0.3071
         -0.3192
         -0.3570
         -0.3581
         -0.3629
         -0.3658
         -0.4798
         -0.5140 - 0.3070i
         -0.5140 + 0.3070i
         -0.9711
         -0.9998
         -1.0000
         -1.0023
         -1.1217



Written by Tan Ek Wee                    - 30 -                              s40205388
                                                                      CHAPTER 4: CASE STUDY


Hence with the Eigenvalues, we are able to get the damp and freq of each machine and
the modes that they are operating in:
>> [mac_state freq(1:48) damp(1:48)]

  1.0000     1.0000    1.0000      0         1.0000      1st column represents the n-1 states
  2.0000     2.0000    1.0000      0.0105    0.8495      2nd column represent the control state number
  3.0000     3.0000    1.0000      0.0105    0.8495      3rd column represent the generator number
  4.0000     4.0000    1.0000      0         1.0000
  5.0000     5.0000    1.0000      0         1.0000      4th column represent the freq
  6.0000     6.0000    1.0000      0         1.0000      5th column represent the damping factor
  7.0000     7.0000    1.0000      0.1021    0.3726
  8.0000     9.0000    1.0000      0.1021    0.3726
  9.0000     10.0000   1.0000      0.1000    0.4502
 10.0000     11.0000   1.0000      0.1000    0.4502
 11.0000     15.0000   1.0000      0.1513    0.8476
 12.0000     16.0000   1.0000      0.1513    0.8476
 13.0000     17.0000   1.0000      0         1.0000
 14.0000     1.0000    2.0000      0         1.0000                  Damping ratio less then 0, at
 15.0000     2.0000    2.0000      0         1.0000                  state 19,20.
 16.0000     3.0000    2.0000      0         1.0000
 17.0000     4.0000    2.0000      0         1.0000
 18.0000     5.0000    2.0000      0         1.0000
 19.0000     6.0000    2.0000      0.5729   -0.0323
 20.0000     7.0000    2.0000      0.5729   -0.0323
 21.0000     8.0000    2.0000      0         1.0000
 22.0000     9.0000    2.0000      0         1.0000
 23.0000     15.0000   2.0000      1.0941    0.0810
 24.0000     16.0000   2.0000      1.0941    0.0810
 25.0000     17.0000   2.0000      1.1402    0.0622
 26.0000     1.0000    3.0000      1.1402    0.0622
 27.0000     2.0000    3.0000      0         1.0000
 28.0000     3.0000    3.0000      0         1.0000
 29.0000     4.0000    3.0000      0         1.0000
 30.0000     5.0000    3.0000      0         1.0000
 31.0000     6.0000    3.0000      1.5749    0.6572
 32.0000     7.0000    3.0000      1.5749    0.6572
 33.0000    15.0000    3.0000      1.6345    0.8823
 34.0000    16.0000    3.0000      1.6345    0.8823
 35.0000    17.0000    3.0000      0         1.0000
 36.0000    1.0000     4.0000      0         1.0000
 37.0000    2.0000     4.0000      0         1.0000
 38.0000    3.0000     4.0000      0         1.0000
 39.0000    4.0000     4.0000      0         1.0000
 40.0000    5.0000     4.0000      0         1.0000
 41.0000    6.0000     4.0000      0         1.0000
 42.0000    7.0000     4.0000      0         1.0000
 43.0000    9.0000     4.0000      4.8855    0.8585
 44.0000    10.0000    4.0000      4.8855    0.8585
 45.0000    11.0000    4.0000      0         1.0000
 46.0000    15.0000    4.0000      0         1.0000
 47.0000    16.0000    4.0000      0         1.0000
 48.0000    17.0000    4.0000      0         1.0000




Written by Tan Ek Wee                           - 31 -                               s40205388
                                                            CHAPTER 4: CASE STUDY


l=
 1.0e+002 *             The following data are obtained from data2_pss. All 4 PSS
                        were installed onto each generator, and by using the
 -0.0000
 -0.0001                svc_mgen.m driver file to simulate the data, the results are
 -0.0010                as shown:
 -0.0010
 -0.0010                The number of total states has increase due to the PSS was
 -0.0020                introduced into the system.
 -0.0020
 -0.0020                system.
 -0.0021                state = 15
 -0.0065
 -0.0032 - 0.0061i              15
 -0.0032 + 0.0061i              13
 -0.0026 - 0.0068i
 -0.0026 + 0.0068i              16
 -0.0126 - 0.0122i              1
 -0.0126 + 0.0122i
 -0.0199                Total = 60 states
 -0.0200
 -0.0200
 -0.0200
 -0.0269
 -0.0363
 -0.0013 - 0.0370i
 -0.0013 + 0.0370i
 -0.0410
 -0.0499
 -0.0061 - 0.0710i
 -0.0061 + 0.0710i
 -0.0097 - 0.0707i
 -0.0097 + 0.0707i
 -0.1000
 -0.1000
 -0.1000                       There are 12 pairs of complex values, each pair
 -0.1000                       representing the oscillatory system modes
 -0.0867 - 0.0988i
 -0.0867 + 0.0988i
 -0.1942 - 0.1074i             The states associated with this mode are
 -0.1942 + 0.1074i
 -0.2838                        sparse(abs(p_norm(:,19)))
 -0.3072
 -0.3198
 -0.3571                                  (42,1)   0.4137
 -0.3585                                  (58,1)   1.0000
 -0.3634
 -0.3661
 -0.4770
 -0.5134 - 0.3062i
 -0.5134 + 0.3062i
 -0.9768 - 0.0172i
 -0.9768 + 0.0172i
 -0.9851 - 0.0909i
 -0.9851 + 0.0909i
 -1.0000
 -1.0000
 -1.0000
 -1.0000
 -1.0000
 -1.0185 - 0.0255i
 -1.0185 + 0.0255i
 -1.1134




Written by Tan Ek Wee            - 32 -                                     s40205388
                                                                      CHAPTER 4: CASE STUDY


>> [mac_state freq(1:59) damp(1:59)]

ans =

 1.0000      1.0000   1.0000      0        1.0000
 2.0000      2.0000   1.0000      0        1.0000
 3.0000      3.0000   1.0000      0        1.0000        1st column represents the n-1 states
 4.0000      4.0000   1.0000      0        1.0000        2nd column represent the control state number
 5.0000      5.0000   1.0000      0        1.0000
 6.0000      6.0000   1.0000      0        1.0000        3rd column represent the generator number
 7.0000      7.0000   1.0000      0        1.0000        4th column represent the freq
 8.0000      9.0000   1.0000      0        1.0000
 9.0000     10.0000   1.0000      0        1.0000        5th column represent the damping factor
 10.0000    11.0000   1.0000      0        1.0000
 11.0000    12.0000   1.0000      0.0975   0.4682
 12.0000    13.0000   1.0000      0.0975   0.4682
 13.0000    15.0000   1.0000      0.1088   0.3591
 14.0000    16.0000   1.0000      0.1088   0.3591
 15.0000    17.0000   1.0000      0.1944   0.7169
 16.0000    1.0000    2.0000      0.1944   0.7169
 17.0000    2.0000    2.0000      0        1.0000                  Damping ratio less then 0.05,
 18.0000    3.0000    2.0000      0        1.0000
 19.0000    4.0000    2.0000      0        1.0000                  at state 23, 24.
 20.0000    5.0000    2.0000      0        1.0000
 21.0000    6.0000    2.0000      0        1.0000
 22.0000    7.0000    2.0000      0        1.0000
 23.0000    8.0000    2.0000      0.5883   0.0346
 24.0000    9.0000    2.0000      0.5883   0.0346
 25.0000    12.0000   2.0000      0        1.0000
 26.0000    13.0000   2.0000      0        1.0000
 27.0000    14.0000   2.0000      1.1297   0.0849
 28.0000    15.0000   2.0000      1.1297   0.0849
 29.0000    16.0000   2.0000      1.1256   0.1363
 30.0000    17.0000   2.0000      1.1256   0.1363
 31.0000    1.0000    3.0000      0        1.0000
 32.0000    2.0000    3.0000      0        1.0000
 33.0000    3.0000    3.0000      0        1.0000
 34.0000    4.0000    3.0000      0        1.0000
 35.0000    5.0000    3.0000      1.5725   0.6596
 36.0000    6.0000    3.0000      1.5725   0.6596
 37.0000    7.0000    3.0000      1.7089   0.8751
 38.0000   12.0000    3.0000      1.7089   0.8751
 39.0000   13.0000    3.0000      0        1.0000
 40.0000   14.0000    3.0000      0        1.0000
 41.0000   15.0000    3.0000      0        1.0000
 42.0000   16.0000    3.0000      0        1.0000
 43.0000   17.0000    3.0000      0        1.0000
 44.0000   1.0000     4.0000      0        1.0000
 45.0000   2.0000     4.0000      0        1.0000
 46.0000   3.0000     4.0000      0        1.0000
 47.0000   4.0000     4.0000      4.8735   0.8588
 48.0000   5.0000     4.0000      4.8735   0.8588
 49.0000   6.0000     4.0000      0.2735   0.9998
 50.0000   7.0000     4.0000      0.2735   0.9998
 51.0000   9.0000     4.0000      1.4474   0.9958
 52.0000   10.0000    4.0000      1.4474   0.9958
 53.0000   11.0000    4.0000      0        1.0000
 54.0000   12.0000    4.0000      0        1.0000
 55.0000   13.0000    4.0000      0        1.0000
 56.0000   14.0000    4.0000      0        1.0000
 57.0000   15.0000    4.0000      0        1.0000
 58.0000   16.0000    4.0000      0.4058   0.9997
 59.0000   17.0000    4.0000      0.4058   0.9997




Written by Tan Ek Wee                           - 33 -                               s40205388
                                                                CHAPTER 4: CASE STUDY


As seen from the above simulation results, the modes of the systems are as listed below.
The table shows the modes of the rotor oscillation prior to and after installation of PSSs
on the four generators




 Mode    Without PSS       Damping        With PSS        Damping         Mode Shift
                            Factor                         Factor
   1     -0.26 ± j0.64      0.3726       -0.32 ± j0.61     0.4682        -0.06 ± j0.03
   2      0.12 ± j0.36     -0.0323       -0.13 ± j0.37     0.0346        -0.25 ± j0.01
   3     -1.82 ± j1.03      0.6572       -1.94 ± j1.07     0.6596        -0.12 ± j0.04
   4     -5.01 ± j3.07      0.8585       -5.14 ± j3.06     0.8588        -0.13 ± j0.01

From the above results, it shows that the installation of the Power Stability System (PSS)
helps to improve system damping.
Hence to conclude, the result uses QR eigen-analysis routine to perform as the system has
less then 500 states. Since all the system Eigenvalues and vectors are then know, it is
easy to provide the performance measures of participation and residue. The PSS approach
is based on the damping and synchronizing torque concept, which provides the basic for
tuning of both these stabilizers in the multi-machine power system.




Written by Tan Ek Wee                     - 34 -                               s40205388
                           CHAPTER 5: SYSTEM IDENTIFICATION TECHNIQUIES


Chapter 5 System identification techniques
There are several types of system identification techniques available in identifying
signals of a power system. In the modern time such as now, transient phenomena has
been carried out in real time and with the availability of cheap computing power, this has
further develop the ease of identification of transients, greatly reduce the modeling
limitations and costly maintenance. These are a few of the identification methods:


1. s-domain identification (frequency domain)
2. z-domain identification (time domain)
3. Recursive least-squares curves fitting algorithm
4. Prony’s method
(Reference from: Power system electromagnetic transient simulation)


s-domain identification (frequency domain)
The rational function in s to be fitted to the frequency domain data is:


                                       a 0 + a1s1 + a 2 s 2 + L + ans N
                            H(s) =
                                     1 + b 0 + b1s1 + b 2 s 2 + L + bns n
                                                                                      A.1

where N <= n
the frequency response of equation A.1 is:


                                                 n
                                                       (ak ( jω ) k )
                                  H ( jω ) =    k =0
                                                  n
                                                       (bk ( jω )k )
                                                k =0
                                                                                      A.2

where b0 = 1
letting the sample data be c(j ) + jd(j ), and equating it to equation A.2 yields




Written by Tan Ek Wee                         - 35 -                                s40205388
                                                      CHAPTER 5: S-DOMAIN IDENTIFICATION



          c( jω ) + jd ( jω ) =
                                ( a − a w − a w −L ) + j ( a w − a w − a w −L )
                                   0     2
                                             2
                                             k    4
                                                      4
                                                      k                 1       1       3
                                                                                                3
                                                                                                k       5
                                                                                                                5
                                                                                                                k

                                  (1 − b w − b w − L ) + j ( b w − b w − b w − L )
                                        2
                                             2
                                             k   4
                                                      4
                                                      k             1       1       3
                                                                                            3
                                                                                            k       5
                                                                                                            5
                                                                                                            k

                                                                                                                              A.3

or

            (c( jω ) + jd ( jω )) (1 − b2 wk − b4 wk − L ) + j ( b1w1 − b3 wk − b5 wk − L
                                           2       4                        3       5
                                                                                                                        )
            = ( a0 − a2 wk2 − a4 wk − L ) + j ( a1w1 − a3 wk − a5 wk − L
                                  4                        3       5
                                                                                        )
                                                                                                                              A.4

Splitting into real and imaginary part yields:


               −d k ( jω ) • ( b1wk − b3 wk − b5 wk − L ) + c( jω ) ( b2 wk2 + b4 wk4 − L
                                          3       5
                                                                                                                    )
               − ( a0 − a2 wk2 + a4 wk − L ) = −c ( jω )
                                     4



               d k ( jω ) • ( b2 wk2 − b4 wk4 − L ) + c( jω ) ( b1w1 − b3 wk − b5 wk − L
                                                                   k
                                                                           3       5
                                                                                                                )
               − ( a1wk − a3 wk + a5 wk − L ) = −d ( jω )
                              3       5




This must hold each sample points and therefore assembling into a matrix equation gives


++++ A.5


Where the terms t1, t2, t3, and t4 are
                   −1π                      tπ 1
       t1 = sin        ω k d ( jω k ) + cos
                         1
                                               ω k c( jω k )
                    2                        2
                  −lπ
       t2 = cos       ω1
                       k
                   2
                   lπ                      tπ
       t3 = cos       ω k d ( jω k ) + sin
                        1
                                              ω 1 c( jω k )
                                                k
                    2                       2
                −lπ
       t4 = sin     ωk
                     1

                 2
       l = column number
       k = row or sample number




Written by Tan Ek Wee                             - 36 -                                                                    s40205388
                                                CHAPTER 5: S-DOMAIN IDENTIFICATION


equation A.5 is of the form:



                                  [ A11 ] [ A12 ]     a
                                                          =
                                                              c
                                  [ A 21 ] [ A 22 ]   b       d
                                                                       A.6

Where
AT = (a0, a1, a2, a3, ….. an)
BT = (b1, b2, b3, ….. bn)
CT = ( -c (j 1), -c (j 2), -c (j 3), … , -c (j k))
DT = ( -d (j 1), -d (j 2), -d (j 3), … , -d (j k))


Equation B.6 is solved for the required coefficient (a’s and b’s).




Written by Tan Ek Wee                        - 37 -                  s40205388
                                                         CHAPTER 5: Z-DOMAIN IDENTIFICATION


z-domain identification (time domain)
When the sampled data consists of samples in time rather than frequency, a rational
function in the z-domain can be identified, provided the system has been excited by
waveform that contain frequency components at which the matching is required. This can
be achieved by multi-sine injection


                          I ( z ) (1 + b1 z −1 + b2 z −2 + K + bn z − n )
                          = ( a0 + a1 z −1 + a2 z −2 + K + an z − n ) V ( z )

                                                                                                B.1


or


                          I ( z ) = − I ( z ) ( b1 z −1 + b2 z −2 + K + bn z − n )
                         + ( a0 + a1 z −1 + a2 z −2 + K + an z − n ) V ( z )

                                                                                                B.2


Taking the inverse transform gives:


                                       n                            n
                         i (t ) = −          bk i ( t − k ∆t ) +          ak v ( t − k ∆t )
                                      k =1                         k =0

                                                                                                B.3


And in matrix form




Written by Tan Ek Wee                                 - 38 -                                  s40205388
                                              CHAPTER 5: Z-DOMAIN IDENTIFICATION


 v(t1 ) v(t1 − ∆t ) v(t1 − 2∆t ) L v(t1 − n∆t ) −i(t1 − ∆t ) −i (t1 − 2∆t ) −i(t1 − n∆t )
 v(t2 ) v(t2 − ∆t ) v(t2 − 2∆t ) L v(t2 − n∆t ) −i (t2 − ∆t ) −i (t2 − 2∆t ) −i (t2 − n∆t )
   M         M            M       M      M            M              M              M
 v(tk ) v(tk − ∆t ) v(tk − 2∆t ) L v(tk − n∆t ) −i (tk − ∆t ) −i (tk − 2∆t ) −i (tk − n∆t )




  a0
          i (t1 )
   a1
          i (t2 )
* M =
            M
  an
          i (tk )
   b1
  b2
   M
  bn
                                                                                     B.4


Where
        .k = time sample number
        .n = order of the rational function (k>p, i.e. over-sampled).
The time step must be choose sufficiently small to avoid aliasing, i.e. t = 1/(K1 fmax),
where K1 > 2 (Nyquist criteria) and fmax is the highest frequency injected.




Written by Tan Ek Wee                       - 39 -                                s40205388
CHAPTER 5: RECURSIVE LEAST SQUARE CURVE-FITTING ALGORITHM


Recursive least square curve-fitting algorithm
The least squares curve-fitting algorithm is describe here to extract the fundamental
frequency data based on least squared error technique. We assume a sinusoidal signal
with a frequency of    radians/sec and a phase shift of             relative to some arbitary time T0
i.e.


                                    y(t) = A sin ( t – )
                                                                                             C.1
where    = T0
This can be rewritten as


                      y(t) = A sin( t) cos( T0) – Acos( t) sin( T0)
                                                                                             C.2
Letting C1 = A cos( T0) and C2 = A sin( T0) and representing sin( t) and cos( t) by
function F1 (t) and F2 (t) respectively, then:


                                 y(t) = C1F1 (t) + C2F2 (t)
                                                                                             C.3
F1 (t) and F2 (t) are known if the fundamental frequency                   is known. However, the
amplitude and phase of the fundamental frequency need to be found, so equation C.3 has
to solve for C1 and C2. if the signal y(t) is distorted, then its deviation from a sinusoid
can be described by an error function E, i.e.


                                       x(t) = y(t) + E
                                                                                             C.4
For the least square method of curve fitting, the size of the error function is measured by
the sum if the individual residual squared values, such that:


                                            n                   2

                                      E=          { xi − yi }
                                           i =1

                                                                                             C.5


Written by Tan Ek Wee                       - 40 -                                        s40205388
CHAPTER 5: RECURSIVE LEAST SQUARE CURVE-FITTING ALGORITHM


where xi = x(t0 + i t) and yi = y(t0 + i t). from the equation C.3


                                         n                                     2

                              E=               { xi − C1 F1 (ti ) − C2 F2 (ti )}
                                        i =1

                                                                                     C.6
where the residual value r at each discrete step is define as:


                                   ri = xi - C1F1 (ti) - C2F2 (ti)
                                                                                     C.7
In matrix form:


                              r1             x1        F1 (t1 )   F2 (t1 )
                              r2             x2        F1 (t2 ) F2 (t2 )     C1
                                    =              −
                              M                M          M          M       C2
                              rn             xn        F1 (tn ) F2 (tn )

                                                                                     C.8
or
                                               [r] = [X] – [F][C]
                                                                                     C.9
The error component can be described in terms of the residual matrix as follows:
E = [r]T[r] = r12 + r22 + … + rn2
 = [[X] – [F][C]]T [[X] – [F][C]]
 = [X]T[X] – [C]T[F]T[X] - [X]T[F][C] + [C]T[F]T[F][C]
                                                                                     C.10
This error then needs to be minimized, i.e.


                            ∂E
                               = −2[ F ]T [ X ] + 2[ F ]T [ F ][C ] = 0
                            ∂C
                                      [ F ]T [ F ][C ] = [ F ]T [ X ]

                                    [C ] = [[ F ]T [ F ]]−1[ F ]T [ X ]
                                                                                     C.11


Written by Tan Ek Wee                                  - 41 -                      s40205388
CHAPTER 5: RECURSIVE LEAST SQUARE CURVE-FITTING ALGORITHM



if [A] = [F]T [F] and [B] = [F]T [X] then
                                               [C] = [A]-1[B]
                                                                                                        C.12
therefore
                            F1                          F1 F1 (t1 )    F1 F2 (t1 )       a11   a12
                   [ A] =
                            F2
                                 [ F1    F2 ] =
                                                        F2 F1 (t1 ) F2 F2 (t1 )
                                                                                     =
                                                                                         a21   a22

elements of matrix[A] can then be derived as follows
                                                                   T
                                                 F1 (t1 )              F1 (t1 )
                                        a11n   =    M                     M
                                                        F1 (tn )       F1 (tn )
                                                 n −1
                                           =            F12 (ti ) +F12 (tn )
                                                 i =1


                                               = a11n−1 + F12 (tn )

                                                                                                        C.13
etc
similarly
                                                    F1 (ti ) x(ti )   b
                                        [b ] =      F2 (ti ) x(ti )
                                                                    = 1
                                                                      b2

and
                                          b1n = b1n−1 + F1 (tn ) x(tn )

                                                                                                        C.14
                                         b2n = b2n−1 + F2 (tn ) x(tn )

                                                                                                        C.15


From the matrix element equation, C1 and C2 can be calculated recursively using
sequential data.




Written by Tan Ek Wee                                     - 42 -                                     s40205388
                                                                          CHAPTER 5: PRONY ANALYSIS


Prony Analysis
Prony analysis identifies a rational function that will have a prescribed time-domain
response[3] PARK, T.W. and BURRUS, C.S. : ”Digital filter design” (John Wiley
Interscience, New York, 1987)
Given the rational Function
                                  Y ( z ) a0 + a1 z −1 + a2 z −2 + L + an z − n
                       H ( z) =          =
                                  U ( z ) 1 + b1 z −1 + b2 z −2 + L + bn z − n
                                                                                        (B.1)
The impulse response of h(k) is related to H(z) by z-transform, i.e.
                                                      ∞
                                         H ( z) =          h ( k ) z −1
                                                    k =0

                                                                                        (B.2)
Which can be written as

                              Y ( z ) (1 + b1 z −1 + b2 z −2 + L + bn z − n )
                              = ( a0 + a1 z −1 + a2 z −2 + L + an z − n )

                                                                                         (B.3)


This is the time domain equivalent of a convolution in the time domain. Using the first L
term of the impulse response the convolution can be expressed as:
                           h0 0           0       L           0             −a0
                           h1 h0          0       L           0             −a1
                                                                      1
                           h2 h1          h0      L           0             M
                                                                      b1
                           M   M           M      O            M           −aN
                                                                      b2 =
                          hn hn −1       hn − 2   L           h0            0
                                                                       M
                          hn+1 hn hn −1 L                     h1            0
                                                                      bn
                            M   M    M O                       M            M
                           hL hL −1 hL − 2 L                 hL −3          0

                                                                                         (B.4)
Partitioning gives:

                                   a
                                     =
                                               [ H1 ]           1
                                   b        [ h1 ] [ H 2 ]      b


Written by Tan Ek Wee                             - 43 -                              s40205388
                                                                    CHAPTER 5: PRONY ANALYSIS


                                                                                   (B.5)


The dimension of the vectors and matrices are:
        a       (N+1) vector
        b       (n+1) vector
        [H1]    (N+1) x (n+1) matrix
        [h1]    vector of the last (L-N) terms of impulse response
        [H2]    (L-N) x (n) matrix


The b coefficient are determined by using the sample points more than n time steps after
the input has been removed. When this occurs the output is no longer a function of the
input (equation B.3) but only depends on the b coefficients and previous output values
(lower partition of equation B.5) i.e.
                                         0 = [ h1 ] + [ H 2 ] b

or
                                             [ h1 ] = − [ H 2 ] b
                                                                                  (B.6)
Once the b coefficients are determined the a coefficients are obtained from the upper
partition of equation B.5, i.e.
                                               b = [ H1 ] a

                                                                                  (B.7)
when L = N + n then H2 is square and , if non-singular, a unique solution for b is
obtained. If H2 is singular many solutions exist, in which case h(k) can be generated by
lower order system.
When m > n + N the system is over-determined and the task is to find a and b coefficients
that gives the best fit (minimized the error). This can be obtained solving equation B.6
using the SVD or the normal equation approach. I.e.

                                  [ H 2 ] [ h1 ] = − [ H 2 ] [ H 2 ] b
                                         T                    T




Written by Tan Ek Wee                              - 44 -                       s40205388
                                                          CHAPTER 5: REFERENCES


References
   1. ABUR, A. and SINGH, H.: “Time domain modeling of external systems for
      electromagnetic transients programs”, IEEE Transients on Power Systems, 1993
   2. PARK, T.W. and BURRUS, C.S.: “Digital filter design” (John Weley
      Interscience, New York, 1987)
   3. NEVILLE,W. and JOS,.A. “Power systems electromagnetic transient simulation”,
      IEE Power and Energy series, 2003




Written by Tan Ek Wee                 - 45 -                             s40205388
                                                            CHAPTER 6: PRONY ANALYSIS


Chapter 6        Prony Analysis


Introduction of Prony analysis
For the past 40 years or so, many power system have become interconnected so as to
exchange power so as to exchange power and at the same time to keep the operating cost
to its minimum. However, due to the interconnection between neighboring power system,
although most of the time it is not over loaded, are often relatively weak when compared
to the connections within the system. The concern were the damping, the damping
torques are lower across the weak ties, and whenever a small disturbances occurred, low
frequency oscillation may happen on the transmission lines. But with the help of modern
technology, computer software and programs are created to help the studies of power
system oscillation, stability and other complex calculations. The necessary complex
mathematics is made transparent to the program user, and the emphasis studies may be
place on the physical nature of the oscillatory phenomena and the accuracy of the model
data.


Prony analysis is a signal processing method that extends Fourier analysis. It is a
technique of analyzing signals to determine model, damping, phase, frequency and
magnitude information contain within the signal


Prony analysis
Prony method is a technique for modeling sampled data as a linear combination of
exponential, it has a close relationship to least squares linear prediction algorithm used
for AR (Autoregressive) and ARMA (Autoregressive moving average) parameter
estimation. A spectral interpretation of Prony’s method may be obtained by computing
the energy spectral density (ESD) of the deterministic exponential model.


There are three basic steps in the Prony method:
Step 1 :         Determine the linear prediction parameters that fit the available data.




Written by Tan Ek Wee                        - 46 -                                s40205388
                                                    CHAPTER 6: PRONY ANALYSIS CONCEPT


Step 2 :                Find the roots of the polynomial formed in step 1, find the prediction
                        coefficient that will yield the estimates of damping and sinusoidal
                        frequencies of each of the exponential term.
Step 3 :                With the solution obtain in step 2, we obtain a second linear equation.
                        With the equation, we yield the estimate of the exponential amplitude and
                        sinusoidal initial phase.


Prony analysis concept
Ref: Digital spectral Analysis by Marple.S, Lawrence.Jr

Since there are many data sample that are being used, such as the exponential parameters,
then an exact exponential fit to the data may be made. Consider the p-exponent discrete-
time function
            p
x ( n) =          hk zkn−1
           k =1

                                                                                              P.1
2p complex samples x[1],…,x[2p] are used to fit an exact exponential model to the 2p
complex parameters h1 ,…,hp , z1,…,zp. The p equation of P.1 for 1<n<p may be
expressed in matrix form as :
  z10       z20
                       L      z0p
                                      h1   x[1]
    1
  z1        z12        L      zp1
                                      h2   x[2]
                                         =
   M         M         O       M       M     M
   p −1
 z1        z2p −1      L     z p −1
                               p
                                      hp   x[ p ]

                                                                                              P.2
The matrix of the time-index z element has a Vandermonde structure. If a method can be
found to separately determine the z-elements, then equation P.2 represents a set of linear
simultaneous equations that can be solved for the unknown vector of the complex
amplitudes. Prony’s contribution was the discovery of such a method.
The key to the separation is to recognize that equation P.1 is the solution to some
homogenous linear constant-coefficient difference equation. In order to find the form of
this difference equation, first we must define the polynomial that has the zk exponents as
its roots,


Written by Tan Ek Wee                               - 47 -                               s40205388
                                                                         CHAPTER 6: PRONY ANALYSIS CONCEPT


             p
φ ( z ) = ∏ z − zk
            k =1

                                                                                               P.3
If the product of P.3 are expand into power series, then polynomial may be represented as
the summation
             p
φ ( z ) = ∏ a [ m] z p −m
            m =1

                                                                                               P.4
With the complex coefficient a[m] such that a[0] = 1. Shifting the index on equation P.1
from n to n-m and multiplying by the parameter a[m] results,
                                      p
a [ m] x [ n − m ] = a [ m ]                hk zk − m −1
                                                n

                                     k =1

                                                                                               P.5
Forming similar products a[0]x[n], … , a[m-1]x[n-m-1] and summing produces
 p                              p           p
       a [ m ] x [ n − m] =          hi          a [ m ]zin − m −1
m =1                          i =0        m= 0

                                                                                               P.6
Which is valid for p+1 <= n <= 2p. making the substitution zin − m −1 = zin− p zip − m −1 ,
 p                             p                   p
       a [ m]x [ n − m ] =           hi zin − p          a [ m ] zip − m −1 = 0
m=0                           i =0                m =0

                                                                                               P.7
The right hand equation then sums up in equation P.7 may be recognized as the
polynomial defines by the equation P.4, evaluated at each of its roots zi, yielding the zero
result indicated. Hence the equation P>7 is the linear difference equation whose
homogeneous solution is given by equation P.1, the polynomial P.4 is the characteristic
equation associated with this linear difference equation,
The p equations representing the valid values of a[n] the satisfy the equation P.7 may be
expressed as the p x p matrix equation




Written by Tan Ek Wee                                                    - 48 -              s40205388
                                                  CHAPTER 6: PRONY ANALYSIS CONCEPT


    x[ p]    x[ p − 1]       L     x[1]     a[1]    x[ p + 1]
  x[ p + 1]    x[ p ]        L     x[2]     a[2]    x[ p + 2]
                                                  =
      M          M           O       M        M         M
 x[2 p − 1] x[2 p − 2]       L     x[ p ]   a[ p]    x[2 p ]
                                                                                P.8
Equation P.8 demonstrates that, with 2p complex data samples, it is possible to decouple
the hk and zk parameters. The complex polynomial coefficients a[1], …, a[p], which are
functions of only the time dependent components zk of the exponential model, from a
linear predictive relationship among the time samples.
The Prony procedure to fit p exponentials to 2p data samples may be summarized in three
major steps as explain earlier in this chapter. The damping α i and sinusoidal frequency

f i may be determined from the root zi using the relationships

f i = tan −1 Im { zi } / Re { zi } / 2π T

                                                                                P.9
α i = In zi / T
                                                                                P.10
To complete the prony procedure, the roots computed in the second step are then used to
construct the matrix element of P.2, which is then solved for the p complex parameters
h[1], …, h[p]. The amplitude Ai and the initial phase θ i may be determined from each hi
parameter with the relationships
Ai = hi

                                                                                P.11
θ i = tan −1 Im {hi } / Re {hi }

                                                                                P.12




Written by Tan Ek Wee                             - 49 -                     s40205388
                                                             CHAPTER 6: PRONY PROGRAM


Prony Program
With the above description of the mathematical approach for prony analysis, it is applied
to program in “Matlab”. It is programmed to model sampled data of a linear combination
of exponentials. An input signal is created to run the simulation, and then it seeks to fit a
deterministic exponential model to data.
The following are the steps for Prony Program
       •   Acquire Data
               o N sample
               o T seconds/sample
       •   Select Input Parameter
               o Number of complex exponentials,
       •   Estimate Complex Exponential
       •   Compute real Exponential parameters


%%% ---------------------------------------------------------------------------------- %%%
clear all
% Declare the Time intervals%
t = [0:0.001:1];

% Input signal given as y(t) = 150*exp(-3*t).*sin(2*pi*51*t+(pi/4)) +
% 160*sin(2*pi*49*t+(pi/5))
Xx = -3*t;
Yy = (2*pi*51*t+(pi/4));
Zz = (2*pi*49*t+(pi/5));
input_signal = 150*exp(Xx).*sin(Yy) + 160*sin(Zz);

% Plot out the input signal waveform %
subplot(211),plot(t,input_signal)
title('Input Signal')
xlabel('Time')

% Fourier Transform of the input %
FY = fft(input_signal,512);
Pyy = FY.*conj(FY) / 512;
f = 1000*(0:256)/512;
% Plot out the signal wave form in frequency domain %
subplot(212),plot(f,Pyy(1:257))




Written by Tan Ek Wee                        - 50 -                                  s40205388
                                                          CHAPTER 6: PRONY PROGRAM


%%% ------------------------- Step 1 ------------------------------ %%%

%order = input('Please Enter the number of Orders required :       ')
order = 4;
%Sampling time = 0.001
T = 0.001;

% Step 1, find the linear prediction A y(k) = a(1)y(k-1) + a(2)y(k-2) + ..... + a(n)y(k-n)
% Then applied repleatedly to form matrix.

% [n x 2] matrix of the system in second order, X
X=[];
m = length(input_signal) - order;
step = order;
for i = 1:order
   for j = 1:m
      X(i,j) = input_signal(step-1 + j);
   end
   step = step-1;
end

z = [];
for l = 1:(length(input_signal)-order)
   z(l,1) = input_signal(l+order);
end

% Find Theta using Moore-Penrose pseudo-inverse of matrix X
% theta = (X)++ * z

%%% ------------------------- Step 2 ------------------------------ %%%
X=X';
theta = pinv(X)*z;
% Roots of the LPM
% z^n - [a(1)z^n-1 + a(2)z^n-2 + ... + a(n)] = 0
LPM = [1 -theta']
% Roots of the LPM
rootz = roots(LPM);
Freq = [];
Freqe = [];
Frequency = atan(imag(rootz)/real(rootz))./(2*pi*T);
for i = 1:1:order
   if Frequency(:,i) ~ 0
      Freq = Frequency(:,i);
   end
end
Freq
damping_factor = real(log(rootz)/T)


Written by Tan Ek Wee                      - 51 -                                s40205388
                                                          CHAPTER 6: PRONY PROGRAM


 %%% ------------------------- Step 3 ------------------------------ %%%

 Z=[];
 for k = 1:order
    for m = 1:order
       Z(k,m) = rootz(m)^(k-1);
    end
 end
 V = [];
 for n = 1:order
    V(n) = input_signal(n);
 end
 H = inv(Z)*V';
 phase_rad = atan(imag(H)/real(H));
 phase_rad = phase_rad(:,2)

 amplitude = sqrt(imag(H).^2 + real(H).^2)

 %%% --------------------------------- -------------------------------- %%%




Written by Tan Ek Wee                      - 52 -                             s40205388
                                                       CHAPTER 6: PRONY PROGRAM


The figure on the right shows the
output results of the prony
program. The first waveform
shows the input signal

y(t) = 150*exp(-3*t).*sin
       (2*pi*51*t+(pi/4)) +
        160*sin
       (2*pi*49*t+(pi/5))

The second waveform on the
figure shows the Fourier
transform signal, as shown, the
main shoots is at approximate
near 50Hz




The Prony program has the following values calculated out via prony.
     Frequency            Damping Factor           Amplitude            Phase
       51.000                -3.000                  75                -0.7854
      -51.000                -3.000                  75                 0.7854
       49.000                 0.000                  80                -0.9425
      -49.000                 0.000                  80                 0.9425

Given the input signal

                         Amplitude                               Frequency



y(t) = 150*exp(-3*t).*sin (2*pi*51*t+(pi/4)) + 160*sin (2*pi*49*t+(pi/5))


          Damping

                                  Phase



Written by Tan Ek Wee                     - 53 -                             s40205388
                                    CHAPTER 6: PRONY VERSUS EIGEN ANALYSIS


Prony versus Eigen analysis
Both of the analysis method has its own advantages and disadvantages. It is essential to
consider some of the important factors that includes: the system non-linearity, the size of
the model that is going to be analysis, the difficulty in the selection of suitable processing
parameters, and the reliability of the results.
In the full state eigen-analysis, the size of the model is limited to 500 states, which means
the max data set with 2000 buses and 300 machines must be reduced using dynamic
equivalencing to about 30 machine, not to forget the other dynamic devices such as
HDVC lines, PSS and SVC. As for prony analysis, there is no limit to the size of the
model, because only the output is being analyzed. This is the main advantage over Eigen
analysis, as the standard stability study results are directly useable. This eliminates the
task of deriving a medium-scale model (500 states) that is required for the linearization
method.
The selection of the processing parameters is dependent on the robustness of the
numerical algorithm. Sensibility should be determined for each power system model to
be analyzed. Base on the numerical results documented in [3], the results shows that
Prony analysis appears to be more sensitive to changes of the processing parameters.
The reliability of the results is directly related to the algorithm robustness.


Overall, Eigenvalue methods have been shown to be effective in the study of small
disturbance stability of power system. However the method involves the use of system
matrices and often requires a large amount of computer storage and CPU time


The Prony Analysis Method uses time domain information, such as swing curves, to
extract the dominant system modes, it is suitable for multi-frequency problems and can
be generally applied to power system damping problems. Prony Analysis requires the
mode to be excited with pulse-type function. The method construct a discrete linear
autogressive model to fit the time domain data recorded. From the model the
characteristic equation are then found to give the amplitude and the phase for each mode
in the signal.




Written by Tan Ek Wee                        - 54 -                               s40205388
                                   CHAPTER 6: PRONY VERSUS EIGEN ANALYSIS


In general multi-output swing curves are available form the power system but some
Prony methods use single output signal. Individual signals that are analyzed
independently often results in conflicting estimates due to noise effects. There is another
article stating that the work of Trudnowski et al [1] consider a simple extension of the
Prony Analysis that allows multiple signals to be analyzed simultaneously to give one set




Written by Tan Ek Wee                      - 55 -                                s40205388
                                                             CHAPTER 6: REFERENCES


   References
   [1] Trudnowski D .J, C .J Smith And D. A. Pierre, “An Application Of Prony
   Methods In PSS Design For Multi-Machine Systems” IEEE Trans. Power System,
   Vol 6, No 1, 1991, Pp118-124


   [2] Marple.S, Lawrence.Jr “Digital spectral Analysis” Prentice-Hall, 1987


   [3] C. E. Grund, . J. J. Paserba, J. F. Hauer, S. Nelsson “ Comparison Of Prony And
   Eigenanalysis For Power System Control Design” Ge Industrial And Power System
   Engineering Department Schenectady.


   [4] Trudnowski D .J, J.M. Johnson, J . F. Hauer “Making Prony Analysis More
   Accurate Using Multiple Signals” IEEE Trans. Power Systems, Vol 14, No.1, 1999




Written by Tan Ek Wee                    - 56 -                                s40205388
                                                                     CHAPTER 7: CASE STUDY


Chapter 7 Case study


Simulation using Power System Block-set
The Power System Block-set (PSB) provided by the software “MatLab” application was
design to provide a modern design tool that will allow Engineers to rapidly and easily
build model that will simulate the design power system. The power system block-set uses
SimuLink environment, allowing the user to create the power system model using just
some click and drag procedure. Not only it is easy to create model, but also the analysis
of the circuit can be included in its interactions with mechanical, thermal, control, and
many other disciplines. Due to the fact that Simulink has a extensive modeling library, it
is possible to simulate all the electrical parts of the simulation, hence in this case study, it
uses Simulink and Power System Block-set to test the Prony program and research on the
different results obtain for different scenario.


Description of the Power system
The Power system model that is going to be used in this case study uses single
compensated transmission system (single phase). It consist of a 735 kv, 300 km line is
used to transmit the power across from Bus B1 (735kv equivalent system) to the next bus
B2 , followed by a transformer to reduce the equivalent kv to 315lv load. In order to
simplify the system, only one phase of the system were used here. In this simulation, in
order to make it more practical, the transmission line is series compensated at its center
by a capacitor representing 40% of the line reactance, and therefore increases the
transmission capacity. The line is also shunt compensate at both ends by 330 Mvar shunt
reactance (110 Mvar / Phase). The series capacitor is protected by a metal oxide varistor
(MOV) simulated by the Surge arrester block. The 250 MVA, 735 kV/ 315 kV
transformer is a Saturable Transformer block simulating one phase of the three-phase
750 MVA transformer. A Multimeter block is used to monitor the fault current as well as
the flux and magnetizing current of the transformer.




Written by Tan Ek Wee                        - 57 -                                 s40205388
                                            CHAPTER 7: POWER SYSTEM BLOCKSET




The above diagram is the model used to test the prony program. As seen in the diagram,
faults were applied at bus 2 with the switching time at 10th cycle and recovered at 10.01
cycle. Transient fault normally last a lot faster then that. There are three scopes in
particular. And the following are the results of the scope 1 and 2




Written by Tan Ek Wee                       - 58 -                                 s40205388
                                                                    CHAPTER 7: SIMULATION




Observation of the scopes
Start the simulation, and the waveforms were formed at the scopes. Fault is being applied
at t = 10 cycles, a small dc voltage is being applied to bus B2. As observe, the fault
current reaches 10kA (trace 1 scope 2) during the fault, the MOV conduct at every half
cycle (trace 2 scope 1), but due to the short transient applied, it was not able to see the
effect clearly. And the voltage across the capacitor (trace 1 scope 1) is limited to 260kv.
At t = 10.01 cycles, the fault is cleared. The 15 Hz modes are clearly seen on the
capacitor voltage (trace 1 of Scope1) and bus B2 voltage (Scope3). During fault the flux
in the transformer is trapped to around 1 pu. At fault clearing the flux offset and 15 Hz
components cause transformer saturation (flux>1.2 pu, trace 3 of Scope2), producing
magnetizing current pulses (trace 2 of Scope2).




Written by Tan Ek Wee                       - 59 -                                 s40205388
                                                                  CHAPTER 7: CASE STUDY


But the main observation will be focusing on is the scopedata on bus B2 voltage
(Scope3), the data obtain before and after the transient being applied and the results using
prony program. With the results from the scope3 , it is analysis and hence determines the
phase, Amplitude, damping factor and lastly the frequency of the signal given.


Results
Using the Prony program, but the signals are analyzed, first the transient signal. The
scope shown below is the result obtained, and the results after applying the prony
program.




     Frequency           Damping Factor             Amplitude                Phase
      121.5482             -113.3798                 0.0094                  0.5524
     -121.5482             -113.3798                 0.0094                 -0.5524
      59.9773               -0.5534                  0.0204                  0.7790
      -59.9773              -0.5534                  0.0204                 -0.7790




Written by Tan Ek Wee                      - 60 -                                s40205388
                                                                 CHAPTER 7: CASE STUDY


The scope data below shows the power system without fault applied, and the results are
as shown




    Frequency           Damping Factor             Amplitude               Phase
     101.0087             -205.6958                 0.0001                0.0073
    -101.0087             -205.6958                 0.0001                -0.0073
      60.0000               0.3726                  0.5057                 1.5363
     -60.0000               0.3726                  0.5057                -1.5363

The results obtain tally with the result found using the Power System Blockset. The
Power GUI, part of the PSB provides results such as measurements of frequencies,
amplitude and states of the machines etc. the results shown are very close to the
calculated values via Prony program. The program help analyzing the oscillatory
behavior of power system, using prony as a identification tool. Hence with the results
obtain, Engineers could further use is to identify modes of the system, which will be use
in the study and control of power system oscillations.




Written by Tan Ek Wee                     - 61 -                                s40205388
                 CHAPTER 8: CONCLUSION AND FUTURE CONSIDERATIONS


Chapter 8      Conclusion and Future Considerations

To conclude, the thesis of “Power System Dynamic Stability via Prony Analysis”
introduces several identification techniques to clearly identify the modes of oscillations in
a power system model. It is essential to have special model identification tools, as it can
be used with either simulated or measured responses, they are also very useful for
validating linear and non-linear simulation models against observed system behavior


Power system engineers have many signal analysis tools to select from, and Prony
analysis uses time domain correlation through parametric spectrum analysis, which has a
number of advantages over Eigen analysis. Both Prony and Eigen analysis can produce
similar results in the determination of the modal parameters and the construction of the
controller design model, however Prony analysis requires care in the choice of processing
parameters and should be used in combination with other methods. This thesis paper only
describes the concept of Prony, hence more development to the Prony analysis is possible
and it will become a very promising tool.


In this thesis, after power system stability fundamentals, Prony's method is presented,
then power system models are developed using power system block-set (PSB) in
MATLAB software for Prony analysis on stability characteristics. Observations were
made both before and after fault applied, results were than compared with the modes of
the system. The research outcome has shown that Prony analysis is robust and efficient
for power system stability assessment, especially in view of its modal analysis.


Prony analysis are conducted when new generating and transmitting facilities are
planned, the studies are helpful in determining such things as the nature of the network
system, modal damping controllers and other devices and controls.




Written by Tan Ek Wee                       - 62 -                                 s40205388
                 CHAPTER 8: CONCLUSION AND FUTURE CONSIDERATIONS


Future considerations
In this thesis, it only touches on the Prony analysis concept as one of the modal
identification techniques. And the simulation is mainly base on a single-phase system,
which in reality, actual power system model mainly comprises of three-phase system and
in a large power system model. Hence it is possible to further develop the Prony Program
with other combinations to increase the reliability of the analysis.


The development of both the modal analysis and modal identification tools have been
accelerated by the occurrence of lightly-damped and the unstable oscillations that were
observed in practice or in planning simulations. Hence it is necessary to have knowledge
of modal analysis and its application to the power system. A complete understanding of
power system oscillations generally requires a combination of analytical tools. Therefore
for future consideration, a thorough study on modal analysis and modal identification will
be appropriate, and a case study, which is related to the real world of power system
industry, helps to further prove that Prony analysis is a better approach of most
electromechanical transient analysis, especially in damping mode identification.




Written by Tan Ek Wee                       - 63 -                             s40205388
                                                                  APPENDIX : svm_mgen.m


Appendix A

% svm_mgen.m
% 6:27 PM 18/8/97
% m.file to generate state variable models
% This m file takes the load flow and dynamic data from a data m file
% and calulates the state space matrices:
%             A matrix in a_mat
%             B matrices
%             for a change in Exciter Vref in b_vr
%             for a change in Turbine governer Pref b_pr
%             C matrices
%             change in generator speed c_sp
%             change in generator electrical torque c_t on generator base
%             change in generator electrical power c_p on generator base
%             change in bus voltage angles c_ang
%             change in bus voltage magnitude c_v
%             change in from_bus line active power c_pf1
%             change in from_bus line reactive power c_qf1
%             change in to_bus line active power c_pf2
%             change in to_bus line reactive power c_qf2
%             change in from bus current magnitude c_ilmf
%             change in to bus current magnitude c_ilmt
%             change in from bus real current c_ilrf
%             change in to bus real current c_ilrt
%             change in from bus imaginary current c_ilif
%             change in to bus imaginary current c_ilit
%             change in field voltage c_Efd
%             D matrices
%             combination of output and input, e.g.,
%             for power out and p_ref in --- d_ppr
%             change in Field voltage for a change in Vref d_Efdvr
% l is the eigenvalue vector of a_mat
% u is the right eigenvector matrix of a_mat
% v is the left eigenvector matrix of a_mat (vu = I)
% p is the unscaled participation matrix
% p_norm is the scaled participation matrix
% the maximum value of each column in p_norm is unity
% all scaled participations less than 0.1 are set to zero
% to find the states associated with the jth eigenvalue
% use sparse(abs(p_norm(:,j)))
% pr gives the participation factors of the eigenvalues
% associated with the rotor angle states of each generator
% use sparse(abs(pr(k,:))) to find the modes associated
% with the rotor angle of the kth generator
%



Written by Tan Ek Wee                      - 64 -                           s40205388
                                                                  APPENDIX : svm_mgen.m


% Author: Graham Rogers
% Modified: August 1997
% Induction Generator added
% Modified: August 1997
% load modulation and output matrices for line flow added
% Modified: April 1997
% HVDC added
% Version 1.0
% Date: September 1996
% Copyright - Joe Chow/Cherry Tree Scientific Software 1991-1997
%

clear
clear global
close %graphics windows
% set up global variables
pst_var

jay = sqrt(-1);
%
%
% load input data from m.file
disp('linearized model development by perturbation of the non-linear model')
[dfile,pathname]=uigetfile('d*.m','Select Data File');
if pathname == 0
     error(' you must select a valid data file')
else
     lfile =length(dfile);
     % strip off .m and convert to lower case
     dfile = lower(dfile(1:lfile-2));
     eval(dfile);
end
% check for valid dynamic data file
if isempty(mac_con)
   error(' the selected file is not a valid data file')
end
sys_freq = input('enter the base system frequency in Hz - [60]');
if isempty(sys_freq);sys_freq = 60;end
basrad = 2*pi*sys_freq; % default system frequency is 60 Hz
basmva = input('enter system base MVA - [100]');
if isempty(basmva);basmva = 100;end % default 100 MVA base
syn_ref = 1 ; % machine 1 is reference by default

disp(' ')
lf = input('do you wish to perform a post fault load flow?Y/N[N]','s');
if isempty(lf); lf = 'N';end



Written by Tan Ek Wee                      - 65 -                              s40205388
                                                                     APPENDIX : svm_mgen.m


if lf =='y'; lf = 'Y'; end
if lf == 'Y'
     disp('enter the changes to bus anf line required to give the post fault condition')
     disp('when you have finished, type return and press enter')
     keyboard
end

% solve for loadflow - loadflow parameter
ans = input('Do you want to solve loadflow > (y/n)[y] ','s');
if isempty(ans);ans='y';end
if ans=='Y';ans='y';end
if ans == 'y'
   if isempty(dcsp_con)
        n_conv = 0;
        n_dcl = 0;
        tol = 1e-9; % tolerance for convergence
        iter_max = 30; % maximum number of iterations
        acc = 1.0; % acceleration factor
        [bus_sol,line,line_flw] = ...
        loadflow(bus,line,tol,iter_max,acc,'n',2);
        bus = bus_sol; % solved loadflow solution needed for
        % initialization
        save sim_fle bus line n_conv n_dcl
   else
       [bus_sol,line,line_flw,rec_par,inv_par, line_par] = lfdcs(bus,line);
       bus = bus_sol;
       save sim_fle bus line rec_par inv_par line_par
   end
else
        load sim_fle
end

n_exc= 0;
n_dc =    0;
n_smp = 0;
n_st3 = 0;
n_pss=    0;
n_tg =    0;
n_svc = 0;
n_lmod = 0;
n_rlmod = 0;
% note dc_indx called in load flow
mac_indx;% identifies generators
ntot = n_mac+n_mot+n_ig;
ngm = n_mac+n_mot;
pm_sig = zeros(n_mac,2);



Written by Tan Ek Wee                       - 66 -                                 s40205388
                                                                 APPENDIX : svm_mgen.m


mac_exc=0;
% check for infinite buses
if n_ib~=0
   %remove controls associated with infinite bus generators
   %remove exciters
   if ~isempty(exc_con)
      n_exc = length(exc_con(:,1));
      net_idx = zeros(n_exc,1);
      for j = 1:n_ib
         net_idx = net_idx|exc_con(:,2)==mac_con(mac_ib_idx(j),1);
      end
      if length(net_idx)==1;
         if net_idx ==1;exc_con=[];end
      else
         perm = diag(~net_idx);
         perm = perm(~net_idx,:);
         exc_con = perm*exc_con;
      end
   end
   % remove pss
   if ~isempty(pss_con)
      n_pss = length(pss_con(:,1));
      net_idx = zeros(n_pss,1);
      for j = 1:n_ib
         net_idx = net_idx|pss_con(:,2) == mac_con(mac_ib_idx(j),1);
      end
      if length(net_idx)==1
         if net_idx == 1;pss_con = [];end
      else
         perm = diag(~net_idx);
         perm = perm(~net_idx,:);
         pss_con = perm*pss_con;
      end
   end
   % remove turbine/governos
   if ~isempty(tg_con)
         n_tg = length(tg_con(:,1));
         net_idx= zeros(n_tg,1);
         for j=1:n_ib
         net_idx =net_idx| tg_con(:,2) == mac_con(mac_ib_idx(j),1);
      end
      if length(net_idx)==1
            if net_idx==1;tg_con = [];end
      else
            perm = diag(~net_idx);
            perm = perm(~net_idx,:);



Written by Tan Ek Wee                     - 67 -                        s40205388
                                                      APPENDIX : svm_mgen.m


      tg_con = perm*tg_con;
    end
  end
end

if ~isempty(exc_con)
   exc_indx;%identifies exciters
   mac_exc = mac_int(exc_con(:,2));
else
   n_exc=0;
end
mac_pss=0;
if ~isempty(pss_con)
     pss_indx;%identifies power system stabilizers
     mac_pss=mac_int(pss_con(:,2));
else
   n_pss=0;
end
mac_tg=0;
if ~isempty(tg_con)
     tg_indx;%identifies turbine/governor
     mac_tg = mac_int(tg_con(:,2));
else
   n_tg =0;
end
if ~isempty(svc_con)~=0
   f = svc_indx; %identifies svc
else
   n_svc = 0;
end
if ~isempty(lmod_con)~=0
   f = lm_indx; % identifies load modulation buses
else
   n_lmod = 0;
end
if ~isempty(rlmod_con)~=0
   f = rlm_indx; % identifies load modulation buses
else
   n_rlmod = 0;
end

%initialize induction motor
if n_mot~=0
    vdp = zeros(n_mot,2);
    vqp = zeros(n_mot,2);
    slip = zeros(n_mot,2);



Written by Tan Ek Wee                     - 68 -             s40205388
                                                                 APPENDIX : svm_mgen.m


  dvdp = zeros(n_mot,2);
  dvqp = zeros(n_mot,2);
  dslip = zeros(n_mot,2);
end
bus = mac_ind(0,1,bus,0);

%initialize induction generator
if n_ig~=0
    vdpig = zeros(n_ig,2);
    vqpig = zeros(n_ig,2);
    slig = zeros(n_ig,2);
    dvdpig = zeros(n_ig,2);
    dvqpig = zeros(n_ig,2);
    dslig = zeros(n_ig,2);
    tmig = zeros(n_ig,2);
end
bus = mac_igen(0,1,bus,0);

%initialize svc
if n_svc ~=0
    B_cv = zeros(n_svc,2);
    dB_cv = zeros(n_svc,2);
    B_con = zeros(n_svc,2);
    dB_con = zeros(n_svc,2);
end
bus = svc(0,1,bus,0);

if n_conv~=0
    % pick up HVDC initial variables from load flow
    Vdc(r_idx,1) = rec_par(:,2); Vdc(i_idx,1) = inv_par(:,2);
    i_dc(r_idx,1) = line_par; i_dc(i_idx,1) = line_par;
    i_dcr(:,1) = i_dc(r_idx,1); i_dci(:,1) = i_dc(i_idx,1);
    alpha(:,1) = rec_par(:,1)*pi/180;
    gamma(:,1) = inv_par(:,1)*pi/180;
end
f = dc_cont(0,1,bus,0); % initialize the dc controls - sets up data for red_ybus
% this has to be done before red_ybus is used since the induction motor,svc and
% dc link initialization alters the bus matrix
v = ones(length(bus(:,1)),2);
bus_v=v;
theta = zeros(length(bus(:,1)),2);
disp(' ')
disp('Performing linearization')
% set line parameters
if ~isempty(lmon_con)
  R = line(lmon_con,3); X = line(lmon_con,4); B = line(lmon_con,5);



Written by Tan Ek Wee                     - 69 -                               s40205388
                                                          APPENDIX : svm_mgen.m


  tap = line(lmon_con,6); phi = line(lmon_con,7);
end
% step 1: construct reduced Y matrix
[Y_gprf,Y_gncprf,Y_ncgprf,Y_ncprf,V_rgprf,V_rncprf,boprf] = red_ybus(bus,line);
bus_intprf = bus_int;% store the internal bus numbers for the pre_fault system
nbus = length(bus(:,1));
if isempty(load_con)
   nload = 0;
else
   nload = length(load_con(:,1));
end
state = zeros(n_mac,1);
gen_state = state;
TR_state = state;
TB_state = state;
TA_state = state;
Efd_state = state;
R_f_state = state;
pss1_state = state;
pss2_state = state;
pss3_state = state;
tg_state = state;
state = zeros(n_mac+n_mot+n_ig+n_svc+n_lmod+n_rlmod+n_dcl,1);
max_state = 6*n_mac+5*n_exc+3*n_pss+3*n_tg+3*n_mot+3*n_ig+2*n_svc +n_lmod
+n_rlmod+ 5*n_dcl;
%17 states per generator,3 per motor, 3 per ind. generator,
% 2 per SVC,1 per lmod,1 per rlmod, 5 per dc line
theta(:,1) = bus(:,3)*pi/180;
v(:,1) = bus(:,2).*exp(jay*theta(:,1));
if n_conv ~= 0
     % convert dc LT to Equ HT bus
     Pr = bus(rec_ac_bus,6);
     Pi = bus(inv_ac_bus,6);
     Qr = bus(rec_ac_bus,7);
     Qi = bus(inv_ac_bus,7);
     VLT= v(ac_bus,1);
     i_acr = (Pr-jay*Qr)./conj(VLT(r_idx));
     i_aci = (Pi - jay*Qi)./conj(VLT(i_idx));
     v(rec_ac_bus,1) = VLT(r_idx) + jay*dcc_pot(:,2).*i_acr;
     v(inv_ac_bus,1) = VLT(i_idx) + jay*dcc_pot(:,4).*i_aci;
     theta(ac_bus,1) = angle(v(ac_bus,1));
end
bus_v(:,1) = v(:,1);
v(:,2) = v(:,1);
bus_v(:,2)=v(:,1);
theta(:,2) = theta(:,1);



Written by Tan Ek Wee                 - 70 -                           s40205388
                                                     APPENDIX : svm_mgen.m


% find total number of states
ns_file
NumStates = sum(state);
exc_sig = zeros(n_mac,2);
if n_tg ~=0
    tg_sig = zeros(n_tg,2);
else
    tg_sig = zeros(1,2);
end
if n_svc ~=0
    svc_sig = zeros(n_svc,2);
else
    svc_sig = zeros(1,2);
end
if n_lmod ~= 0
    lmod_sig = zeros(n_lmod,2);
else
    lmod_sig = zeros(1,2);
end
if n_rlmod ~= 0
    rlmod_sig = zeros(n_rlmod,2);
else
    rlmod_sig = zeros(1,2);
end

if n_conv ~= 0
    dc_sig = zeros(n_conv,2);
else
    dc_sig = zeros(2,2);
end
% set initial state and rate matrices to zero

eterm = zeros(n_mac,2);
pelect = zeros(n_mac,2);
qelect = zeros(n_mac,2);
psi_re = zeros(n_mac,2);
psi_im = zeros(n_mac,2);
psi = zeros(n_mac,2);
mac_ang = zeros(n_mac,2);
mac_spd = zeros(n_mac,2);
edprime = zeros(n_mac,2);
eqprime = zeros(n_mac,2);
psikd = zeros(n_mac,2);
psikq = zeros(n_mac,2);
dmac_ang = zeros(n_mac,2);
dmac_spd = zeros(n_mac,2);



Written by Tan Ek Wee                       - 71 -          s40205388
                                             APPENDIX : svm_mgen.m


dedprime = zeros(n_mac,2);
deqprime = zeros(n_mac,2);
dpsikd = zeros(n_mac,2);
dpsikq = zeros(n_mac,2);
if n_exc~=0
    V_TR = zeros(n_exc,2);
    V_As = zeros(n_exc,2);
    V_A = zeros(n_exc,2);
    V_R =zeros(n_exc,2);
    Efd = zeros(n_exc,2);
    R_f = zeros(n_exc,2);
    dV_TR = zeros(n_exc,2);
    dV_As = zeros(n_exc,2);
    dV_R =zeros(n_exc,2);
    dEfd = zeros(n_exc,2);
    dR_f = zeros(n_exc,2);
    pss_out = zeros(n_exc,2);
end
if n_pss~=0
    pss1 = zeros(n_pss,2);
    pss2 = zeros(n_pss,2);
    pss3 = zeros(n_pss,2);
    dpss1 = zeros(n_pss,2);
    dpss2 =zeros(n_pss,2);
    dpss3 =zeros(n_pss,2);
end
if n_tg~=0
    tg1 = zeros(n_tg,2);
    tg2 = zeros(n_tg,2);
    tg3 = zeros(n_tg,2);
    dtg1 =zeros(n_tg,2);
    dtg2 = zeros(n_tg,2);
    dtg3 = zeros(n_tg,2);
end
if n_lmod~=0
     lmod_st = zeros(n_lmod,2);
    dlmod_st = zeros(n_lmod,2);
end
if n_rlmod~=0
    rlmod_st = zeros(n_rlmod,2);
    drlmod_st = zeros(n_rlmod,2);
end

%HVDC links
if n_conv~= 0
    dv_conr = zeros(n_dcl,2);



Written by Tan Ek Wee               - 72 -          s40205388
                                                                APPENDIX : svm_mgen.m


   dv_coni = zeros(n_dcl,2);
   di_dcr = zeros(n_dcl,2);
   di_dci = zeros(n_dcl,2);
   dv_dcc = zeros(n_dcl,2);
end
% set dimensions for A matrix and permutation matrix
a_mat = zeros(NumStates);
p_mat = sparse(zeros(NumStates,max_state));
c_spd = zeros(length(not_ib_idx),NumStates);
c_p = zeros(length(not_ib_idx),NumStates);
c_t = zeros(length(not_ib_idx),NumStates);

%determine p_mat: converts the vector of length max_states to
%a column of a_mat or b

p_m_file

% step 2: initialization
flag = 0;
%machines
f = mac_em(0,1,bus,flag);
f= mac_tra(0,1,bus,flag);
f = mac_sub(0,1,bus,flag);
f = mac_ib(0,1,bus,flag);
%calculate initial electrical torque
psi = psi_re(:,1)+jay*psi_im(:,1);
if n_mot~=0&n_ig==0
    vmp = vdp(:,1) + jay*vqp(:,1);
    int_volt=[psi; vmp]; % internal voltages of generators and motors
elseif n_mot==0&n_ig~=0
    vmpig = vdpig(:,1) + jay*vqpig(:,1);
    int_volt = [psi; vmpig]; % int volt of synch and ind generators
elseif n_mot~=0&n_ig~=0
    vmp = vdp(:,1) + jay*vqp(:,1);
    vmpig = vdpig(:,1) + jay*vqpig(:,1);
     int_volt = [psi; vmp; vmpig];
else
     int_volt = psi;
end
cur(:,1) = Y_gprf*int_volt; % network solution currents into generators
b_v(boprf(nload+1:nbus),1) = V_rgprf*int_volt; % bus voltage reconstruction
if nload~=0
     vnc = nc_load(bus,flag,Y_ncprf,Y_ncgprf);%vnc is a dummy variable
    cur(:,1) = cur(:,1) + Y_gncprf*v(bus_intprf(load_con(:,1)),1);% modify currents for
nc loads
end



Written by Tan Ek Wee                    - 73 -                               s40205388
                                                                APPENDIX : svm_mgen.m


cur_re(1:n_mac,1) = real(cur(1:n_mac,1));
cur_im(1:n_mac,1) = imag(cur(1:n_mac,1));
cur_mag(1:n_mac,1) = abs(cur(1:n_mac,1)).*mac_pot(:,1);
if n_mot~=0
    idmot(:,1) = -real(cur(n_mac+1:ngm,1));%induction motor currents
    iqmot(:,1) = -imag(cur(n_mac+1:ngm,1));%current out of network
end
if n_ig~=0
    idig(:,1) = -real(cur(ngm+1:ntot,1));%induction generator currents
    iqig(:,1) = -imag(cur(ngm+1:ntot,1));%current out of network
end

if n_conv ~=0
    % calculate dc voltage and current
    V0(r_idx,1) = abs(v(rec_ac_bus,1)).*dcc_pot(:,7);
    V0(i_idx,1) = abs(v(inv_ac_bus,1)).*dcc_pot(:,8);
    Vdc(r_idx,1) = V0(r_idx,1).*cos(alpha(:,1)) - i_dcr(:,1).*dcc_pot(:,3);
    Vdc(i_idx,1) = V0(i_idx,1).*cos(gamma(:,1)) - i_dci(:,1).*dcc_pot(:,5);
    Vdc_ref = Vdc(i_idx,1);
    i_dc(r_idx,1) = i_dcr(:,1);
    i_dc(i_idx,1) = i_dci(:,1);
end

telect(:,1) = pelect(:,1).*mac_pot(:,1)...
   + cur_mag(:,1).*cur_mag(:,1).*mac_con(:,5);
%pss
f = pss(0,1,bus,flag);
%exciters
f = smpexc(0,1,bus,flag);
f = exc_dc12(0,1,bus,flag);
f = exc_st3(0,1,bus,flag);
% turbine governors
f = tg(0,1,bus,flag);
if ~isempty(lmod_con)
   disp('load modulation')
   f = lmod(0,1,bus,flag);
end
if ~isempty(rlmod_con)
   disp('reactive load modulation')
   f = rlmod(0,1,bus,flag);
end

% initialize non-linear loads
if ~isempty(load_con)
   disp('non-linear loads')
   vnc = nc_load(bus,flag,Y_ncprf,Y_ncgprf);



Written by Tan Ek Wee                     - 74 -                              s40205388
                                        APPENDIX : svm_mgen.m


else
  nload = 0;
end
% hvdc lines
f = dc_line(0,1,bus,flag);

mach_ref(1) = 0;
mach_ref(2) = 0;
sys_freq(1) = 1;
sys_freq(2) = 1;

%set states
%generators
mac_ang(:,2) = mac_ang(:,1);
mac_spd(:,2) = mac_spd(:,1);
eqprime(:,2) = eqprime(:,1);
psikd(:,2) = psikd(:,1);
edprime(:,2) = edprime(:,1);
psikq(:,2)=psikq(:,1);
%exciters
if n_exc~=0
     V_TR(:,2)=V_TR(:,1);
     V_As(:,2) = V_As(:,1);
     V_A(:,2) = V_A(:,1);
     V_R(:,2)=V_R(:,1);
     Efd(:,2)=Efd(:,1);
     R_f(:,2)=R_f(:,1);
end
%pss
if n_pss~=0
    pss1(:,2)=pss1(:,1);
    pss2(:,2)=pss2(:,1);
    pss3(:,2)=pss3(:,1);
end
%turbine governor
if n_tg~=0
    tg1(:,2)=tg1(:,1);
    tg2(:,2)=tg2(:,1);
    tg3(:,2)=tg3(:,1);
end
telect(:,2) =telect(:,1);
if n_mot~=0
    vdp(:,2) = vdp(:,1);
    vqp(:,2) = vqp(:,1);
    slip(:,2) = slip(:,1);
end



Written by Tan Ek Wee          - 75 -          s40205388
                                                                    APPENDIX : svm_mgen.m


if n_ig~=0
    vdpig(:,2) = vdpig(:,1);
    vqpig(:,2) = vqpig(:,1);
    slig(:,2) = slig(:,1);
    tmig(:,2) = tmig(:,1);
end
if n_svc ~= 0
     B_cv(:,2) = B_cv(:,1);
     B_con(:,2) = B_con(:,1);
end
if n_lmod ~=0
     lmod_st(:,2) = lmod_st(:,1);
end
if n_rlmod ~=0
     rlmod_st(:,2) = rlmod_st(:,1);
end
if n_conv~=0
     v_conr(:,2) = v_conr(:,1);
     v_coni(:,2) = v_coni(:,1);
     i_dcr(:,2) = i_dcr(:,1);
     i_dci(:,2) = i_dci(:,1);
     v_dcc(:,2) = v_dcc(:,1);
end
% set interconnection variables in perturbation stage to defaults
% this accounts for any generators which do not have the
% corresponding controls
eterm(:,2) = eterm(:,1);
pmech(:,2) = pmech(:,1);
vex(:,2) = vex(:,1);
exc_sig(:,2) = exc_sig(:,1);
tg_sig(:,2) = tg_sig(:,1);
svc_sig(:,2) = svc_sig(:,1);
lmod_sig(:,2) = lmod_sig(:,1);
rlmod_sig(:,2) = rlmod_sig(:,1);

if n_conv~=0
    Vdc(:,2) = Vdc(:,1);
    i_dc(:,2) = i_dc(:,1);
    dc_sig(:,2) = dc_sig(:,1);
    cur_ord(:,2) = cur_ord(:,1);
    alpha(:,2) = alpha(:,1);
    gamma(:,2) = gamma(:,1);
end
%perform perturbation of state variables

p_cont



Written by Tan Ek Wee                      - 76 -                          s40205388
                                                                    APPENDIX : svm_mgen.m



% setup matrix giving state numbers for generators
mac_state = zeros(sum(state(1:n_mac)),3);
for k = 1:n_mac
   if state(k)~=0
      if k == 1
         j = 1;
      else
         j = 1+sum(state(1:k-1));
      end
      jj = sum(state(1:k));
      mac_state(j:jj,1) = (j:jj)';
      mac_state(j:jj,2) = st_name(k,1:state(k))';
      mac_state(j:jj,3) = k*ones(state(k),1);
   end
end
ang_idx = find(mac_state(:,2)==1);
% Form transformation matrix to get rid of zero eigenvalue
% Use generator 1 as reference
% check for infinite buses
if isempty(ibus_con)
     ref_gen = input('Do you want to set gen 1 as reference, Y/N(N)','s');
     if isempty(ref_gen);ref_gen='N';end
     if ref_gen=='n';ref_gen='N';end
     if ref_gen ~= 'N'
          p_ang = eye(NumStates);
          p_ang(ang_idx,1) = -ones(length(ang_idx),1);
          p_ang(1,1) = 1;
          p_angi = inv(p_ang);
         %transform state matrix
         a_mat = p_ang*a_mat*p_angi;
         %transform the c matrices
         c_v = c_v*p_angi;
         c_ang = c_ang*p_angi;
         c_spd = c_spd*p_angi;
      if ~isempty(lmon_con)
         c_pf1 = c_pf1*p_angi;
         c_qf1 = c_qf1*p_angi;
         c_pf2 = c_pf2*p_angi;
         c_qf2 = c_qf2*p_angi;
      end
      %transform the b matrices
      if n_tg~=0;b_pr = p_ang*b_pr;end
      if n_exc~=0;b_vr = p_ang*b_vr;end
      if n_svc~=0;b_svc = p_ang*b_svc;end
      if n_lmod~=0;b_lmod = p_ang*b_lmod;end



Written by Tan Ek Wee                       - 77 -                           s40205388
                                                             APPENDIX : svm_mgen.m


  end
end
disp('calculating eigenvalues and eigenvectors')
%eigenvectors and eigenvalues of a_mat
[u l] = eig(a_mat); % u is the right eigenvector

% sort the eigenvalues
[l l_idx] =sort( diag(l));


%reorder the eigenvector matrix
u = u(:,l_idx);

for j = 1:NumStates
  if imag(l(j))~=0
     %scale the complex eigenvectors so that the maximum element is unity modulus
     [maxu,mu_idx] = max(abs(u(:,j)));
     u(:,j) = u(:,j)/maxu;
  end
end
v = inv(u); % left eigenvectors
% find the participation factors
for j = 1:NumStates
  p(:,j) = (conj(v(j,:)))'.*u(:,j);% p are the unnormalized participation vectors
  [p_max,p_max_idx] = max((p(:,j)));
  p_norm(:,j) = p(:,j)/p(p_max_idx,j);% p_norm has biggest element = 1
  p_big = abs(p_norm(:,j))>0.1;%big sorts out normalized participation > 1
  p_norm(:,j) = p_big.*p_norm(:,j);% p_norm now contains only values of p_norm > 0.1
end

% find states associated with the generator angles
pr = p_norm(ang_idx,:);
% frequency and damping ratio
freq = abs(imag(l))/2/pi;
zero_eig = find(abs(l)==0);
if ~isempty(zero_eig)
    damp(zero_eig)= ones(length(zero_eig),1);
else
    nz_eig = find(abs(l)~=0);
    damp(nz_eig,1) = -real(l(nz_eig))./abs(l(nz_eig));
end


disp(' the a matrix is contained in a_mat')
disp(' l is the vector of eigenvalues')
disp(' u is the matrix of right eigenvectors')



Written by Tan Ek Wee                        - 78 -                       s40205388
                                                                    APPENDIX : svm_mgen.m


disp(' v is the matrix of left eigenvectors')
disp(' p is the matrix of paticipation factors(sensitivity to changes in diagonal of a_mat)')
disp(' p_norm is the matrix of normalized participation factors')
disp(' greater than 0.1')
disp(' pr shows the normalized participation factors associated')
disp(' with the generator rotor angles')
disp(' the damping ratios are given in damp')
disp(' the frequencies are given in freq')
disp(' ')
disp(' c_p is the output matrix corresponding to generator powers')
disp(' c_v is the output matrix corresponding to bus voltage magnitudes')
disp(' c_ang is the output matrix corresponding to bus voltage angles')
disp(' c_pf1 is the output matrix for the line power flow at the from buses')
disp(' c_pf2 is the output matrix for the line power flow at the to buses')
disp(' c_qf1 is the output matrix for the line reactive power flow at the from buses')
disp(' c_qf2 is the output matrix for the line reactive power flow at the to buses')
disp(' c_spd is the output matrix corresponding to the generator speeds')
disp(' c_t is the output matrix corresponding to generator electrical torque')
disp(' c_Efd is the output matrix corresponding to field voltage output from the exciter')
disp(' ')
disp(' b_pr is the input matrix corresponding to governor p_ref')
disp(' b_vr is the input matrix corresponding to the exciter v_ref')
disp(' b_svc is the input matrix corresponding to the svc reference input')
disp(' b_lmod is the input matrix corresponding to the lmod input signal')
disp(' b_rlmod is the input matrix corresponding to the rlmod input signal')

disp(' ')
disp(' d_ppr is the d matrix for power out and governor p_ref change in')
disp(' d_pvr is the d matrix for power out and the exciter v_ref change in')
disp(' d_vpr is the d matrix for bus voltage magnitude out and governor p_ref change in')
disp(' d_vvr is the d matrix for bus voltage magnitude out and exciter v_ref change in')
disp(' d_angpr is the d matrix for bus voltage angle out and governor p_ref change in')
disp(' d_angvr is the d matrix for bus voltage angle out and exciter_ref change in')
disp(' d_Efdvr is the d matrix for field voltage out and exciter_ref change in')
disp(' ')
disp('speed is always the second state in each generator')
disp(' ')
disp(' the state variable numbers for the kth generator')
disp(' are contained in st_name(k,:)')
disp(' mac_state(j,2) gives the control state number corresponding')
disp(' to state j')
disp(' mac_state(j,3) gives the generator corresponding to state j')
disp(' ')
us_idx = find(damp<0);
ud_idx = find(damp>0&damp<0.05);
s_idx = find(damp>0.05);



Written by Tan Ek Wee                       - 79 -                                s40205388
                                                                APPENDIX : svm_mgen.m


flag = 0;
while(flag == 0)
  disp('You can examine the system modal characteristics')
  disp('Type 1 to see all unstable eigenvalues')
  disp(' 2 to see all eigenvalues with damping ratio less than 0.05')
  disp(' 3 to see all eigenvalues with damping ratio greater than 0.05')
  disp(' 4 to see participation factors associated with the unstable eigenvalues')
  disp(' 5 to see paticipation factors associated with the modes having damping ratio
          less than 0.05')
  disp(' 6 to see participation factors associated with rotor angle modes')
  disp(' 0 to quit and plot your own curves')
  sel = input('enter selection >> ');
  if isempty(sel);sel = 0;end
  if sel == 1
     if isempty(us_idx)
         disp('there are no unstable modes')
     else
         plot(real(l(us_idx)),imag(l(us_idx)),'rx')
         title('unstable modes')
         xlabel('real part')
         ylabel('imaginary part')
     end
     disp('paused: press any key to continue')
     pause
  elseif sel == 2
     if isempty(ud_idx)
         disp('there are no stable modes modes with damping ratio < 0.05')
     else
         plot(real(l(ud_idx)),imag(l(ud_idx)),'go')
         title('stable modes with damping < 0.05')
         xlabel('real part')
         ylabel('imaginary part')
     end
     disp('paused: press any key to continue')
     pause
  elseif sel == 3
     if isempty(s_idx)
         disp('there are no modes with damping ratio > 0.05')
     else
         plot(real(l(s_idx)),imag(l(s_idx)),'b+')
         title('modes with damping > 0.05')
         xlabel('real part')
         ylabel('imaginary part')
     end
     disp('paused: press any key to continue')
     pause



Written by Tan Ek Wee                    - 80 -                              s40205388
                                                                    APPENDIX : svm_mgen.m


 elseif sel == 4
   if isempty(us_idx)
        disp('there are no unstable modes')
   else
        ml = length(us_idx);
        str_ml = num2str(ml);
        mode = input(['select an index number between 1 and ',str_ml]);
        if ~isempty(mode)
            bar(sparse(abs(p_norm(:,us_idx(mode)))))
            mode_num = num2str(us_idx(mode));
            str_damp = num2str(damp(us_idx(mode)));
            str_freq = num2str(freq(us_idx(mode)));
            title(['participation factors of mode ',mode_num,'; damping ',str_damp,';
            frequency ',str_freq])
            ylabel(' participation magnitude')
            xlabel(' state number')
      end
   end
   disp('paused: press any key to continue')
   pause
 elseif sel == 5
   if isempty(ud_idx)
         disp('there are no stable modes with damping ratio < 0.05')
   else
         ml = length(ud_idx);
         str_ml = num2str(ml);
         mode = input(['select an index number between 1 and ',str_ml]);
       If ~isempty(mode)
           bar(sparse(abs(p_norm(:,ud_idx(mode)))))
           mode_num = num2str(ud_idx(mode));
           str_damp = num2str(damp(ud_idx(mode)));
           str_freq = num2str(freq(ud_idx(mode)));
           title(['participation factors of mode ',mode_num,'; damping ',str_damp,';
                 frequency ',str_freq])
            ylabel(' participation magnitude')
            xlabel(' state number')
      end
   end
   disp('paused: press any key to continue')
   pause
 elseif sel == 6
      mn = input('input the generator number');
   if isempty(mn);mn = 0;end;
   if mn >0&mn<n_mac
        bar(sparse(abs(pr(mn,:))))
        str_mn = num2str(mn);



Written by Tan Ek Wee                      - 81 -                                 s40205388
                                                                     APPENDIX : svm_mgen.m


      title(['participation factors associated with number ',str_mn,' generator rotor'])
      ylabel(' participation magnitude')
      xlabel(' mode number')
      disp('paused: press any key to continue')
      pause
    end
  elseif sel == 0
    flag = 1;
  else
    error('invalid selection')
  end
end
clear global




Written by Tan Ek Wee                       - 82 -                                 s40205388
                                                                       APPENDIX : data2.m



% Two Area Test Case
% sub transient generators with static exciters, turbine/governors
% 50% constant current active loads
% load modulation

disp('Two-area test case with subtransient generator models')
disp('Static exciters')
disp('turbine/governors')
% bus data format
% bus:
% col1 number
% col2 voltage magnitude(pu)
% col3 voltage angle(degree)
% col4 p_gen(pu)
% col5 q_gen(pu),
% col6 p_load(pu)
% col7 q_load(pu)
% col8 G shunt(pu)
% col9 B shunt(pu)
% col10 bus_type
%      bus_type         - 1, swing bus
%                       - 2, generator bus (PV bus)
%                       - 3, load bus (PQ bus)
% col11 q_gen_max(pu)
% col12 q_gen_min(pu)
% col13 v_rated (kV)
% col14 v_max pu
% col15 v_min pu


bus = [...
  1 1.03       18.5    7.00    1.61    0.00 0.00 0.00 0.00 1 5.0 -1.0 22.0    1.1 .9;
  2 1.01       8.80    7.00    1.76    0.00 0.00 0.00 0.00 2 5.0 -1.0 22.0    1.1 .9;
  3 0.9781     -6.1    0.00     0.00   0.00 0.00 0.00 3.00 3 0.0 0.0 230.0    1.5 .5;
  4 0.95       -10     0.00     0.00    9.76 1.00 0.00 0.00 3 0.0 0.0 115.0   1.05 .95;
  10 1.0103    12.1    0.00    0.00    0.00 0.00 0.00 0.00 3 0.0 0.0 230.0    1.5 .5;
  11 1.03      -6.8    7.16    1.49    0.00 0.00 0.00 0.00 2 5.0 -1.0 22.0    1.1 .9;
  12 1.01      -16.9   7.00    1.39    0.00 0.00 0.00 0.00 2 5.0 -1.0 22.0    1.1 .9;
  13 0.9899    -31.8   0.00     0.00   0.00 0.00 0.00 5.00 3 0.0 0.0 230.0     1.5 .5;
  14 0.95      -35     0.00     0.00   17.65 1.00 0.00 0.00 3 0.0 0.0 115.0    1.05 .95;
  20 0.9876     2.1     0.00    0.00   0.00 0.00 0.00 0.00 3 0.0 0.0 230.0     1.5 .5;
  101 1.00     -19.3   0.00     1.09   0.00 0.00 0.00 0.00 2 2.0 0.0 500.0     1.5 .5;
  110 1.0125   -13.4   0.00     0.00   0.00 0.00 0.00 0.00 3 0.0 0.0 230.0     1.5 .5;
  120 0.9938   -23.6   0.00     0.00   0.00 0.00 0.00 0.00 3 0.0 0.0 230.0     1.5 .5 ];




Written by Tan Ek Wee                        - 83 -                            s40205388
                                                                         APPENDIX : data2.m




% line data format
% line: from bus, to bus, resistance(pu), reactance(pu),
% line charging(pu), tap ratio, tap phase, tapmax, tapmin, tapsize

line = [...
1    10       0.0      0.0167   0.00     1.0    0.   0.    0.    0.;
2    20       0.0      0.0167   0.00     1.0    0.   0.    0.    0.;
3    4        0.0      0.005    0.00     1.0    0.   1.2   0.8   0.05;
3    20       0.001    0.0100   0.0175   1.0    0.   0.    0.    0.;
3    101      0.011    0.110    0.1925   1.0    0.   0.    0.    0.;
3    101      0.011    0.110    0.1925   1.0    0.   0.    0.    0.;
10 20         0.0025   0.025    0.0437   1.0    0.   0.    0.    0.;
11 110        0.0      0.0167   0.0      1.0    0.   0.    0.    0.;
12 120        0.0      0.0167   0.0      1.0    0.   0.    0.    0.;
13 14         0.0      0.005    0.00     1.0    0.   1.2   0.8   0.05;
13 101        0.011    0.11     0.1925   1.0    0.   0.    0.    0.;
13 101        0.011    0.11     0.1925   1.0    0.   0.    0.    0.;
13 120        0.001    0.01     0.0175   1.0    0.   0.    0.    0.;
110 120       0.0025   0.025    0.0437   1.0    0.   0.    0.    0.];

% Machine data format
% Machine data format
%   1. machine number,
%   2. bus number,
%   3. base mva,
%   4. leakage reactance x_l(pu),
%   5. resistance r_a(pu),
%   6. d-axis sychronous reactance x_d(pu),
%   7. d-axis transient reactance x'_d(pu),
%   8. d-axis subtransient reactance x"_d(pu),
%   9. d-axis open-circuit time constant T'_do(sec),
% 10. d-axis open-circuit subtransient time constant
%         T"_do(sec),
% 11. q-axis sychronous reactance x_q(pu),
% 12. q-axis transient reactance x'_q(pu),
% 13. q-axis subtransient reactance x"_q(pu),
% 14. q-axis open-circuit time constant T'_qo(sec),
% 15. q-axis open circuit subtransient time constant
%         T"_qo(sec),
% 16. inertia constant H(sec),
% 17. damping coefficient d_o(pu),
% 18. dampling coefficient d_1(pu),
% 19. bus number
%



Written by Tan Ek Wee                          - 84 -                       s40205388
                                                                        APPENDIX : data2.m


% note: all the following machines use sub-transient model
mac_con = [ ...

1 1 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.24 0.4 0.05...
              6.5 0 0 1;
2 2 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
              6.5 0 0 2;
3 11 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
              6.5 0 0 11;
4 12 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
              6.5 0 0 12];

exc_con = [...
1   1     0.01   46.0  0.06      0        0        1.0 -0.9...
0.0 0.46 3.1     0.33  2.3       0.1      0.1      1.0 0        0 0 ;
3   2     0.01   7.04  1.0       6.67     1.0      10.0 -10.0...
0.2 -0.2 200.0   4.37  20        4.83     0.09     1.1 8.63 1.0 6.53;
0   3     0.01   200.0 0          0        0       5.0 -5.0...
0   0     0      0     0          0        0       0    0      0   0;
2   4     0.01   300.0 0.01       0        0       4.95 -4.9...
1.0 1.33 3.05    0.279 2.29       0.117    0.1     0.675 0 0 0];


pss_con = [];
%1 3 300.0 20.0 0.06 0.04 0.08 0.04 0.2 -0.05];

%pss_con = [...
%1 1 100.0 10.0 0.05 0.01 0.00 0.00 0.2 -0.05;
%1 2 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;
%1 3 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;
%%1 4 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;];

% governor model
% tg_con matrix format
%column           data             unit
% 1 turbine model number (=1)
% 2 machine number
% 3 speed set point wf        pu
% 4 steady state gain 1/R     pu
% 5 maximum power order Tmaxpu on generator base
% 6 servo time constant Ts    sec
% 7 governor time constant Tc sec



Written by Tan Ek Wee                     - 85 -                           s40205388
                                                                        APPENDIX : data2.m


% 8 transient gain time constant T3           sec
% 9 HP section time constant T4 sec
% 10 reheater time constant T5 sec

tg_con = [...
1 1 1 1.0 25.0   0.1   0.5 0.0 1.25 5.0;
1 2 1 1.0 25.0   0.1   0.5 0.0 1.25 5.0;
1 3 1 1.0 25.0   0.1   0.5 0.0 1.25 5.0;
1 4 1 1.0 25.0   0.1   0.5 0.0 1.25 5.0];



load_con = [4 0 0 .5 0;
       14 0 0 .5 0;
       101 0 0 0 0];


svc_con = [ 1 101 100 10 -10 100 0.05];


%Switching file defines the simulation control
%     row 1 col1 simulation start time (s) (cols 2 to 6 zeros)
%     col7 initial time step (s)
%     row 2 col1 fault application time (s)
%     col2 bus number at which fault is applied
%     col3 bus number defining far end of faulted line
%     col4 zero sequence impedance in pu on system base
%     col5 negative sequence impedance in pu on system base
%     col6 type of fault      - 0 three phase
%                             - 1 line to ground
%                             - 2 line-to-line to ground
%                             - 3 line-to-line
%                             - 4 loss of line with no fault
%                             - 5 loss of load at bus
%                             - 6 no action
%     col7 time step for fault period (s)
%     row 3 col1 near end fault clearing time (s) (cols 2 to 6 zeros)
%     col7 time step for second part of fault (s)
%     row 4 col1 far end fault clearing time (s) (cols 2 to 6 zeros)
%     col7 time step for fault cleared simulation (s)
%     row 5 col1 time to change step length (s)
%     col7 time step (s)
%
%
%
% row n col1 finishing time (s) (n indicates that intermediate rows may be inserted)



Written by Tan Ek Wee                       - 86 -                            s40205388
                                               APPENDIX : data2.m



sw_con = [...
0    0 0 0 0 0 0.01;%sets intitial time step
0.1 3 101 0 0 0 0.01; %no fault
0.15 0 0 0 0 0 0.01; %clear near end
0.20 0 0 0 0 0 0.01; %clear remote end
%0.50 0 0 0 0 0 0.01; % increase time step
%1.0 0 0 0 0 0 0.01; % increase time step
10.0 0 0 0 0 0 0]; % end simulation

%ibus_con = [0 1 1 1];




Written by Tan Ek Wee              - 87 -         s40205388
                                                                 APPENDIX : data2_pss.m


data2_pss.m

% Two Area Test Case
% Generator plant replaces generator 1
% with simple exciters and turbine governors on all generators
% pss on generators in plant (1 to 4)


% bus data format
% bus:
% col1 number
% col2 voltage magnitude(pu)
% col3 voltage angle(degree)
% col4 p_gen(pu)
% col5 q_gen(pu),
% col6 p_load(pu)
% col7 q_load(pu)
% col8 G shunt(pu)
% col9 B shunt(pu)
% col10 bus_type
%     bus_type - 1, swing bus
%          - 2, generator bus (PV bus)
%          - 3, load bus (PQ bus)
% col11 q_gen_max(pu)
% col12 q_gen_min(pu)

bus = [...
    1 1.03 18.5 1.75 0.4 0.00 0.00 0.00 0.00 2 10.0 -10.0;
    102 1.03 18.5 1.75 0.4 0.00 0.00 0.00 0.00 2 10.0 -10.0;
    103 1.03 18.5 1.75 0.4 0.00 0.00 0.00 0.00 2 10.0 -10.0;
    104 1.03 18.5 1.75 0.4 0.00 0.00 0.00 0.00 2 10.0 -10.0
           2 1.01 8.80 7.00 1.76 0.00 0.00 0.00 0.00 2 10.0 -10.0;
           3 0.9781 -6.1 0.00 0.00 9.76 1.00 0.00 0.00 3 0.0 0.0;
           10 1.0103 12.1 0.00 0.00 0.00 0.00 0.00 0.00 3 0.0 0.0;
           11 1.03 -6.8 7.16 1.49 0.00 0.00 0.00 0.00 2 10.0 -10.0;
           12 1.01 -16.9 7.00 1.39 0.00 0.00 0.00 0.00 2 10.0 -10.0;
           13 0.9899 -31.8 0.00 0.00 17.67 1.00 0.00 0.00 3 0.0 0.0;
           20 0.9876 2.1 0.00 0.00 0.00 0.00 0.00 0.00 3 0.0 0.0;
   101 1.00 -19.3 0.00 1.09 0.00 0.00 0.00 0.00 1 0.0 0.0;
   110 1.0125 -13.4 0.00 0.00 0.00 0.00 0.00 0.00 3 0.0 0.0;
   120 0.9938 -23.6 0.00 0.00 0.00 0.00 0.00 0.00 3 0.0 0.0];

% line data format
% line: from bus, to bus, resistance(pu), reactance(pu),
% line charging(pu), tap ratio




Written by Tan Ek Wee                      - 88 -                       s40205388
                                                             APPENDIX : data2_pss.m


line = [...
    1     10   0.0      0.0668   0.00     1.   0.;
    102 10     0.0      0.0668   0.00     1.   0.;
    103 10     0.0      0.0668   0.00     1.   0.;
    104 10     0.0      0.0668   0.00     1.   0.;
     2 20      0.0      0.0167   0.00     1.   0.;
     3 20      0.001    0.0100   0.0175   1.   0.;
     3 101     0.011    0.110    0.1925   1.   0.;
     3 101     0.011    0.110    0.1925   1.   0.;
     10 20     0.0025   0.025    0.0437   1.   0. ;
     11 110    0.0      0.0167   0.0      1.   0. ;
     12 120    0.0      0.0167   0.0      1.   0. ;
     13 101    0.011    0.11     0.1925   1.   0. ;
     13 101    0.011    0.11     0.1925   1.   0. ;
     13 120    0.001    0.01     0.0175   1.   0. ;
     110 120   0.0025   0.025    0.0437   1.   0. ;];


% Machine data format
% Machine data format
%    1. machine number,
%    2. bus number,
%    3. base mva,
%    4. leakage reactance x_l(pu),
%    5. resistance r_a(pu),
%    6. d-axis sychronous reactance x_d(pu),
%    7. d-axis transient reactance x'_d(pu),
%    8. d-axis subtransient reactance x"_d(pu),
%    9. d-axis open-circuit time constant T'_do(sec),
% 10. d-axis open-circuit subtransient time constant
%            T"_do(sec),
% 11. q-axis sychronous reactance x_q(pu),
% 12. q-axis transient reactance x'_q(pu),
% 13. q-axis subtransient reactance x"_q(pu),
% 14. q-axis open-circuit time constant T'_qo(sec),
% 15. q-axis open circuit subtransient time constant
%            T"_qo(sec),
% 16. inertia constant H(sec),
% 17. damping coefficient d_o(pu),
% 18. dampling coefficient d_1(pu),
% 19. bus number
%
% note: all the following machines use sub-transient model
mac_con = [ ...

1 1 225 0.16 .00234 1.81 0.30 0.217 7.80 0.022...



Written by Tan Ek Wee                      - 89 -                   s40205388
                                                                APPENDIX : data2_pss.m


              1.76 0.61 0.254 0.9 0.074...
              3.53 0.00 0.00 1 0 0;
2 102 225 0.16 .00234 1.81 0.30 0.217 7.80 0.022...
              1.76 0.61 0.2540 0.9 0.074...
              3.53 0.00 0.00 102 0 0;
3 103 225 0.16 .00234 1.81 0.30 0.217 7.80 0.022...
              1.76 0.61 0.2540 0.9 0.074...
              3.53 0.00 0.00 103 0 0;
4 104 225 0.16 .00234 1.81 0.30 0.217 7.80 0.022...
              1.76 0.61 0.2540 0.9 0.074...
               3.53 0.00 0.00 104 0 0;
5 2 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
 6.5 6.5 0 2 0 0;
6 11 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
 6.5 6.5 0 11 0 0;
7 12 900 0.200 0.0025 1.8 0.30 0.25 8.00 0.03...
              1.7 0.55 0.25 0.4 0.05...
 6.5 6.5 0 12 0 0];

%exciter models all smpexc - additional zeros provided to match other exciter models


exc_con = [...
0 1     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 2     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 3     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 4     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 5     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 6     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0;
0 7     0.02     200.0   0.05   0   0   5.0   -5.0...
0 0     0        0 0     0      0   0   0      0 0];

% pss model on generators 1 to 4

pss_con = [...
1 1 100.0 10.0 0.05 0.01 0.00 0.00 0.2 -0.05;
1 2 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;
1 3 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;



Written by Tan Ek Wee                    - 90 -                              s40205388
                                                                  APPENDIX : data2_pss.m


1 4 100.0 10.0 0.05 0.01 0.05 0.01 0.2 -0.05;];

% governor model
% tg_con matrix format
%column           data                unit
% 1 turbine model number (=1)
% 2 machine number
% 3 speed set point wf            pu
% 4 steady state gain 1/R         pu
% 5 maximum power order Tmaxpu on generator base
% 6 servo time constant Ts        sec
% 7 governor time constant Tc sec
% 8 transient gain time constant T3   sec
% 9 HP section time constant T4 sec
% 10 reheater time constant T5 sec

tg_con = [...
1 1 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 2 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 3 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 4 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 5 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 6 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0;
1 7 1 25.0      1.0   0.1   0.5 0.0 1.25 5.0];




%Switching file defines the simulation control
%     row 1 col1 simulation start time (s) (cols 2 to 6 zeros)
%     col7 initial time step (s)
%     row 2 col1 fault application time (s)
%     col2 bus number at which fault is applied
%     col3 bus number defining far end of faulted line
%     col4 zero sequence impedance in pu on system base
%     col5 negative sequence impedance in pu on system base
%     col6 type of fault      - 0 three phase
%                             - 1 line to ground
%                             - 2 line-to-line to ground
%                             - 3 line-to-line
%                             - 4 loss of line with no fault
%                             - 5 loss of load at bus
%     col7 time step for fault period (s)
%     row 3 col1 near end fault clearing time (s) (cols 2 to 6 zeros)
%     col7 time step for second part of fault (s)
%     row 4 col1 far end fault clearing time (s) (cols 2 to 6 zeros)



Written by Tan Ek Wee                            - 91 -                  s40205388
                                                                 APPENDIX : data2_pss.m


%     col7 time step for fault cleared simulation (s)
%     row 5 col1 time to change step length (s)
%     col7 time step (s)
%
%
%
%     row n col1 finishing time (s) (n indicates that intermediate rows may be inserted)

sw_con = [...
0    0 0 0      0   0   0.01;%sets intitial time step
0.1 3 101 0     0   0   0.005; %apply three phase fault at bus 3, on line 3-101
0.15 0 0 0      0   0   0.005556; %clear fault at bus 3
0.20 0 0 0      0   0   0.005556; %clear remote end
0.50 0 0 0      0   0   0.01; % increase time step
1.0 0 0 0       0   0   0.02; % increase time step
5.0 0 0 0       0   0   0]; % end simulation




Written by Tan Ek Wee                    - 92 -                               s40205388

								
To top