VIEWS: 251 PAGES: 99 CATEGORY: Education POSTED ON: 6/22/2010
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