FEASIBILITY O F S E C U R I T Y ASSESSMENT IN P O W E R SYSTEMS U S I N G F U L L T I M E - D O M A I N S O L U T I O N S IN T H E E M T P by Mazana Lukic B.Sc, Fakultet Elektrotehnike i Racunarstva, Zagreb 1995 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS F O R T H E D E G R E E O F M A S T E R O F APPLIED SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard loer ZO , £ , O O Q T H E UNIVERSITY O F BRITISH C O L U M B I A Abstract FEASIBILITY O F SECURITY ASSESSMENT IN POWER SYSTEMS USING F U L L T I M E - D O M A I N SOLUTIONS IN T H E E M T P by Mazana Lukic Chairperson of the Supervisory Committee: Professor Jose R. Marti Professor Takahide Niimura Department of Electrical and Computer Engineering Security assessment is one of the key functions in planning and operation of power systems. It has become especially important as power systems are now operating under more stressed conditions due to increase in load, lack of initiative for developing the transmission network, and uncertainties in power transfer schedules brought by deregulation. Static security assessment uses well established tools for power system analysis. However, increased complexity of the system and its operation close to the operating limits has brought the attention to the forms of instability that fail to be recognized by static tools. Therefore dynamic security assessment is becoming an increasingly important area of research. In this work the Electromagnetic Transients Program (EMTP), traditionally used for analysis of fast transients, has been considered as a tool for dynamic security analysis. The classical generator model and exponential load model have been coded in the C++ programming language and adapted for full time-domain simulations in the EMTP. The Dynamic Data Exchange ii facility in the Microsoft Windows operating system environment is used to enable on-line switches manipulation and a user interface is developed for simulation monitoring. The proposed models were tested for different systems and compared with the results from the literature obtained by conventional stability programs. The results presented in this work show that the proposed models are applicable for stability analysis using full time-domain solution in the EMTP. Emphasis is made on the validity of the classical generator model restricted to the first second after initiating a disturbance. User interaction with the running simulation is made possible through switches manipulation and visual monitoring that makes the EMTP more suitable for stability type of analysis. The EMTP was successfully tested for simulation of a large power system with 3600 buses. in T A B L E O F C O N T E N T S Abstract ii Table of Contents iv List of figures vii List of tables x List of symbols xi List of abbreviations xii Acknowledgments xiii Chapter 1: INTRODUCTION 1 Chapter 2: POWER SYSTEM STABILITY PROBLEM 5 2.1 Classification of stability - transient vs. long-term stability 5 2.2 Rotor angle stability 7 2.3 Voltage stability 10 2.4 Traditional approaches to stability studies 11 2.4.1 Rotor angle stability analyses 11 2.4.1.1 Modal analysis for small signal stability evaluation 11 2.4.1.2 Equal area criterion 12 2.4.1.3 Time-domain simulation for evaluation of transient rotor angle stability 14 2.4.2 Voltage stability analyses 15 2.4.2.1 Power flow methods, P-V and Q-V curves 15 2.4.2.2 Modal analysis for voltage stability evaluation 17 2.4.2.3 Time-domain simulation in voltage stability studies 18 Chapter 3: POWER SYSTEM SECURITY ASSESSMENT. 20 3.1 Definitions and Criteria 20 3.2 Available Transfer Capability determination 22 3.2.1 Example: A T C computation and static voltage security analysis of the 30 bus test system 28 IV Chapter 4: POWER SYSTEM MODELLING FOR STABILITY ANALYSIS USING THEEMTP 35 4.1 Synchronous machine modelling 35 4.1.1 Classical model of a synchronous machine 36 4.1.2 The swing equation 39 4.1.3 Numerical integration of the swing equation 43 4.1.4 Calculation of active power in single-phase system representation using full time-domain solution 47 4.1.5 Prediction scheme 49 4.1.6 Block diagram of a classical generator model implemented in CONNEC 50 4.2 Load modelling 52 4.2.1 Exponential load 52 4.2.2 Extracting time-varying RMS voltage response from the instantaneous values 54 4.2.3 Representation of PQ loads in the EMTP 59 4.2.4 Block diagram of an exponential load model implemented in CONNEC 61 4.3 Transmission lines and transformers 63 4.3.1 Transmission lines in stability studies 63 4.3.2 Three-phase switching in single-phase time-domain system simulations 66 4.3.3 Representation of transformers in the EMTP for stability studies. 67 Chapter 5: ADAPT A TION OF THE EMTP FOR ON-LINE STABILITY SIMULATIONS. 71 5.1 User-supplied subroutines in the EMTP version MicroTran... 72 5.1.1 ALPHA subroutine 72 5.1.2 CONNEC subroutine 76 5.2 On-line switches manipulation through Dynamic Data Exchange in Microsoft Windows 77 5.2.1 Dynamic Data Exchange 77 5.2.2 Implementation of the DDE for on-line switches manipulation in the EMTP 78 5.3 User interface for on-line simulation monitoring and switches manipulation 82 Chapter 6: CASE STUDIES 84 6.1 Generator-infinite bus case 84 6.2 Three-generator case 91 6.3 Load connected to an infinite bus case 96 6.4 Full time-domain simulation of a large power system 100 Chapter 7: CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 105 Chapter 8: REFERENCES. 108 VI LIST O F F I G U R E S Number Page Figure 1. Time responses of system components and controls 6 Figure 2. Illustration of equal area criterion 13 Figure 3. Framework for on-line A T C calculation 24 Figure 4. Modified IEEE 30-bus test system 29 Figure 5. Transmission line congestion in the base case 30 Figure 6. Reactive power loading of the pool generation units in the base case 30 Figure 7. Algorithm for T T C calculation 31 Figure 8. Changes in the voltage profile due to the bilateral transaction 24-14....32 Figure 9. Changes in the voltage profile due to the bilateral transaction 17-27....33 Figure 10. Changes in the voltage profile due to bilateral transactions 24-14 and 17-27 34 Figure 11. Simplified transient model of a synchronous machine 37 Figure 12. Thevenin equivalent of the network as seen from the C O N N E C subroutine 38 Figure 13. Illustration of integration with the trapezoidal rule 44 Figure 14. Equivalent circuits for differential equations in (a) continuous and (b) discrete time representation 46 Figure 15. Calculating electrical power output Peftk) 49 Figure 16. Classical generator model block diagram 51 Figure 17. Circuit configuration for filtering 60 H z component of viN(t) 55 Figure 18. Variation of filtered voltage magnitude as a function of frequency ....56 Figure 19. Equivalent circuit of the filter discretized with the trapezoidal rule 57 Figure 20. Equivalent circuit for filtering the voltage as seen from the C O N N E C subroutine 59 Figure 21. Load representation by passive admittances 60 Figure 22. Equivalent circuit as seen from the C O N N E C subroutine 61 Figure 23. Block diagram of an exponential load model 62 Figure 24. Nominal n equivalent circuit of a transmission line 64 vii Figure 25. Dommel's equivalent circuit of a constant-parameter ideal line model as seen from the line ends 65 Figure 26. CP line model 65 Figure 27. Equivalent circuit of a single-phase transformer 68 Figure 28. Equivalent Tt-circuit of a transformer with off-nominal turns ratio t 69 Figure 29. Representation of a thyristor in MicroTran 72 Figure 30. Calling C / C + + subroutine (ALPHAC) from the F O R T R A N code 75 Figure 31. Block diagram of the interaction of subroutine modules with the main program through a dynamic link library 76 Figure 32. Nonlinear element connected to a linear network 76 Figure 33. Block diagram of the communication between the programs 79 Figure 34. Block diagram of a server-client communication 81 Figure 35. Circuit for testing the switches manipulation and on-line monitoring 82 Figure 36. Interface for switches manipulation and on-line simulation monitoring of a test system 83 Figure 37. Equivalent circuit of a generator-infinite bus case 85 Figure 38. Rotor angle response to a disturbance with clearing time tc=0.067 s 87 Figure 39. Generator active power after the disturbance with clearing time tc=0.067 s 87 Figure 40. Fault current and voltage for £c=0.067 s 88 Figure 41. Rotor angle response to a disturbance with clearing time Cc-0.083 s 88 Figure 42. Generator active power after the disturbance with clearing time tc=0.083 s 89 Figure 43. Rotor angle response to a disturbance with clearing time fc=0.092s 89 Figure 44. Generator active power after the disturbance with clearing time fc=0.092 s 90 Figure 45. Three-generator nine-bus system 92 Figure 46. Relative rotor angle differences for three machines (a) S2-S1, (b) &-S, forD=0 94 Figure 47. Electrical power outputs of generators 1, 2 and 3 for D=Q 95 V l l l Figure 48. Relative rotor angle differences for three machines (a) S2-S1, (b) &-& for D= 1.5 95 Figure 49. Electrical power outputs of generators 1, 2 and 3 for D= 1.5 96 Figure 50. Load connected to an infinite bus 97 Figure 51. Bus voltages for the initial operating state with constant power l o a d ( t f = £ = 0 ) 97 Figure 52. RMS value of the load voltage at bus B over a longer period of time for 0) 98 Figure 53. Voltage collapse due to a high consumption of a constant power load(a=/?=0)atbusB 99 Figure 54. RMS voltage at bus B for cases of (a) constant power (a=6=0), (b) constant current (ct=8=V) and (c) constant impedance («=/?= 2) load model for power consumption of P= 1.8430 and Q=0.15 at the voltage level of 1 p.u 100 Figure 55. 500 kV transmission network of the reduced WSCC system 102 Figure 56. Comparison of voltages for simulation time of 0.3 s 103 Figure 57. Comparison of voltages for simulation time of 1.2 s 103 I X LIST O F T A B L E S Number Page Table 1. Generator data 92 Table 2. Transmission lines and load data 93 x LIST O F S Y M B O L S Ts - synchronizing torque [Nm] T D - damping torque [Nm] T E - electromagnetic torque [Nm] co - angular speed [rad/s] § - rotor angle [rad] J - inertia constant [kg m2] T M - mechanical torque [Nm] D - damping constant [Js/rad2] H - inertia constant [Ws/VA] p(t) - instantaneous power [VA] XI LIST O F A B B R E V I A T I O N S EMTP - The Electromagnetic Transients Program A T C - Available Transfer Capability NERC - North American Reliability Council DDE - Dynamic Data Exchange WSCC - Western Systems Coordinating Council T T C - Total Transfer Capability T R M - Transmission Reliability Margin E T C - Existing Transmission Commitments C B M - Capacity Benefit Margin EMS - Energy Management System SE - State Estimator SA - Security Analysis COP - Current Operating Plan OASIS - Open Access Same-Time Information System xii A C K N O W L E D G M E N T S I wish to express my sincere gratitude to all the people of the Power Group here at UBC, who have helped me to accomplish this work, accepted me as a colleague and a friend and gave me the support I needed while being far away from home. In particular I would like to thank to: Dr. Jose Marti, for making it possible for me to study here, for dedicated supervision of my work, for his encouragement and invaluable guidance. Dr. Tak Niimura, for his constant support, for his initiative that led me to expand my knowledge on deregulation issues, for his patience and concern. Dr. Hermann Dommel, for introducing me to the world of the EMTP, for his valuable time and all the support and patience that helped me realize this work. The financial assistance of the Natural Science and Engineering Council of Canada and BC Hydro & Power Authority is gratefully acknowledged. I also wish to thank all my friends, here and at home, for being there for me. Finally, I would like to thank my family, for their love and support that encouraged me to follow my dreams. X I U C h a p t e r 1 INTRODUCTION Energy management systems (EMS) used to base their on-line security assessment solely on steady-state analysis, and dynamic security limits were determined in off-line studies. Dynamic security analysis was associated with transient rotor angle stability only and the research was mainly focused on accelerating the process of calculating transient stability limits for application in on-line environments. Reports on slowly developing voltage instabilities in the late 1980s stimulated research on voltage stability phenomena. What has been realized in this more recent work is that for many power systems the security limits can depend on different criteria (steady state adequacy, transient stability, voltage stability or frequency stability), but studying each one of them requires the use of specific software. Transient rotor angle and voltage stability are recognized as the most constraining criteria and growing emphasis is given to time-domain simulations for their analysis. The need to cover both phenomena in on-line and off-line environment has motivated work on a unified approach for security assessment [16],[19]. The motivation of the present work originates in the idea of designing a single tool capable of analysing the overall stability state of the system. The Electromagnetic Transients Program (EMTP) traditionally used for analysis of fast system transients can be considered for that purpose. It offers detailed 1 models of power system components and methodology that is proven to give reliable and accurate results. The objective of this thesis is to determine the feasibility of stability analysis using the EMTP type of solution. Because the models that exist in the EMTP require a considerable amount of data which are often not available for large scale system studies, and which exceeds the needs of stability requirements (i.e. three-phase synchronous machine representation), two new models for the EMTP were developed to allow the basic stability aspects to be examined. The models are tested using the examples from the literature with satisfactory results. Further on, the possibility of adapting the EMTP for on-line simulations (which does not imply the real-time issues) has been analysed. This is a very important feature as security assessment involves studying of a large number of contingencies. The ability to interact with the system model during the simulation is essential as well as monitoring of the system variables and gathering and evaluation of the most relevant information. Therefore, the interface for simulation monitoring and switches manipulation was developed and tested on a simple model. A significant contribution of an attempt to use stability models in the full time-domain environment is testing their validity. While stability models are designed with many approximations that have been accepted over the years for systems operating with large security margins, today's environment may not allow the space for those approximations. It can be expected that more and more accurate analysis will be required in the future, which will include a three-phase representation and all the systems controls and dynamics. 2 Chapter 2 of this thesis discusses power system stability. Stability is classified according to time frames and driving force that causes instability. It is however emphasised that instability is a single phenomenon, the dynamics of different components are interlinked and it is difficult to make a clear separation between, for example, transient rotor angle and voltage stability. More appropriate separation would be according to time frames influenced by the dynamics of power system components. Further on, an overview of traditional approaches to stability analysis is presented. Chapter 3 gives an overview of power system security assessment, its definitions and criteria. Security is classified as static or dynamic depending on considered operating states. The key factor for power system operation obtained from security assessment is the Available Transfer Capability (ATC) of a transmission network. Guidelines for A T C calculation based on the North American Reliability Council (NERC) recommendations are given and an example of static security analysis and A T C determination presented. Chapter 4 explains the modelling of a synchronous machine and exponential load for stability studies in the EMTP. The method of calculating active power and RMS value of voltage from instantaneous values is described. EMTP models for transmission lines and transformers are presented and their application for stability studies explained. Chapter 5 describes the subroutines in MicroTran version of the EMTP used for implementation of the developed stability models and for on-line switches manipulations. The Dynamic Data Exchange (DDE) option in the Microsoft Windows operating system has been implemented for communication between independently running programs. The user interface for on-line switches manipulation and simulation monitoring is presented. 3 Chapter 6 deals with the case studies performed to validate the proposed models. A generator-infinite bus case and the three-generator case are used to test the synchronous machine model. The load model is tested on a single load-infinite bus case. Comparison is made between constant impedance, constant current and constant power load model. The results are presented graphically and the comparison is made with the examples from the literature. Application of the EMTP for simulation of large power systems is tested with the reduced Western Systems Coordinating Council (WSCC) system model. This system includes 3615 buses and 362 generators. The comparison of results obtained with the full time-domain solution and the power flow steady-state solution is presented. Chapter 7 presents the contributions made with this work and concludes the thesis. The recommendations for future work are given at the end. 4 Chapter 2 POWER SYSTEM STABILITY PROBLEM Large interconnected power systems of North America and elsewhere are facing new challenges in system operation caused by the recent deregulation of the power industry. Stability problems are known to engineers from the early days, but the complexity of today's systems requires new approaches and new tools for their analysis. This chapter provides an overview of the power system stability issues. Stability is classified according to time frames and type of phenomena that drives instability. Summary of methods traditionally accepted and used for rotor angle and voltage stability analyses is given in the last section of this chapter. 2.1 Classification of stability - transient vs. long-term stability A power system is considered to be stable if it is in a state of equilibrium under normal operating conditions and it remains in a state of acceptable equilibrium after being subjected to a disturbance [15]. Instability phenomena time frames range from a fraction of a second to tens of minutes, and many power system components and controls play a role in its development. However, depending on the system characteristics and the disturbance, only some of them will significantly participate in a particular 5 instability scenario. Figure 1 shows response time frames for some of the system components and controls. In general, stability can be classified into transient and long-term time frames. Transient stability time span is from zero to about ten seconds. Long-term stability time frame of typically two to three minutes is also addressed in the literature as "mid-term", "post-transient" or "post-disturbance" stability. Instability beyond these time frames is usually associated with major disturbances that are exceeding the normal system design criteria, resulting in separation of the power system into a number of islands. Transient stability Long-term stability Induction motor dynamics Generator/excitation dynamics Mech, switched capacitors/ reactors LTC transformers and distribution voltage regulators Prime mover control Load/power transfer increase Load diversity/ thermostat Excitation limiting , G a s turbine start-up Undervoltage load shedding Powerplant operator Static Var compensators Ge neration change / A G C Generator inertial dynamics Boiler dynamics Line/transf. overload DC DC converter LTCs System operator Protective Relaying L L 1 10000 1 hour 0.1 I 10 1 100 1 minute '1000 10 minutes Time (s) Figure 1. T ime responses of system components and controls Even though the power system instability is a single phenomenon, traditionally it has been studied as one of three distinctive forms of instability: 6 rotor angle, voltage and frequency instability. By this classification complexity of the stability analysis has been reduced, and different methods of calculation and prediction of stability have been introduced. However, one has to keep in mind that for the overall system stability evaluation all aspects have to be properly addressed and analyzed. Since the different forms of instability are interlinked and often happen within the same time frame, it is not always easy to separate them. For example, transient voltage stability is often interlinked with transient rotor angle stability1, and slower forms of voltage stability are interlinked with small-signal rotor angle stability. Only in extreme cases as stated in [23], it is possible to differentiate between pure angle stability and pure voltage stability phenomena. Frequency stability in the interconnected power systems is associated with long-term large-scale system disturbances where the main concern is survival of islanded systems. Up to this date this type of stability has not been considered for on-line security assessment of power systems, and as such is beyond the scope of this work. 2.2 Rotor angle stability The first type of stability problem recognized by power system engineers was related to synchronous machines. During normal operation the mechanical power input is balanced with the electrical output and the generator rotates with a constant synchronous speed dictating the frequency of the system to either 50 or 60 Hz. If the power balance is upset by a disturbance 1 Transient rotor angle stability is often referred in the literature as "transient stability". 7 in the network, such as opening of a transmission line due to a fault, this will result in deviation of the rotating speed of a machine from its nominal value. Since, unless corrective measures are taken, this situation may lead to a blackout and machine damage, the concern is whether or not the rotating machines are able to recover and reach a new steady state following a disturbance. Rotor angle stability is the ability of interconnected synchronous machines to retain their constant speed and remain in synchronism. It is usually categorized either as a small-signal (small-disturbance) or as transient rotor angle stability. The system is small-signal stable if it is able to maintain synchronism under small disturbances due to continuous variations in loads and generation. The disturbances are considered sufficiently small if the system equations can be linearized. Instability caused by this type of disturbances may result either in a steady increase in rotor angle due to lack of sufficient synchronizing torque or in rotor oscillations of increasing amplitude due to insufficient damping torque. The change in electrical torque of a synchronous machine following a disturbance can be formulated as the combination of two torque components: ATe=TsAS + TDAco (2.1) where TsAS is the synchronizing torque component and TDACO is the damping torque component. Transient rotor angle stability represents the ability of the power system to maintain synchronism when subjected to a severe transient disturbance. It results in large deviations of generator rotor angles and it is influenced by the well known nonlinear power-angle relationship [5]: 8 p = EV sin 8 (2.2) X where P is the power output of the synchronous machine and 8 is the angle of its rotor when the machine with its internal voltage E\5_ behind the reactance Xis connected for example to the infinite bus of V\0°_. Stability in this case depends on the initial state of the system as well as on the severity of the disturbance. In the stable case the rotor angle increases to a maximum, then decreases and oscillates with decreasing amplitude until it reaches a steady state (usually different to the one prior to the disturbance). In unstable cases there are two possible ways in which the angle may respond, one is called the first-swing instability and it occurs when the rotor angle increases steadily until synchronism is lost. This is usually due to insufficient synchronizing torque. In the second case, the system is stable for the first swing, but becomes unstable due to growing oscillations. This case can be considered as if the initially stable rotor angle response to a transient disturbance reaches a new post-disturbance steady state which is small-signal unstable. The fundamental equation in stability studies describes the effect of imbalance between the electromagnetic and mechanical torque of the synchronous machine, and is often referred to as the swing equation (2.3). (mech. rad/s), and Tm and Te are mechanical and electromagnetic torque respectively in (N-m). (2.3) where J is a moment of inertia (kg-m2), co is angular velocity of the rotor 9 2.3 Voltage stability Voltage stability has traditionally been analyzed as a steady-state problem suitable for static power flow analyses. This approach is based on the fact that voltage instability is related to the ability of the system to transfer reactive power from generation to loads. However, the network maximum power transfer limit does not necessarily correspond to the voltage stability limit. Since a power system is a dynamic system, voltage instability has to be considered as a dynamic process. This dynamics are usually associated with the loads and the means for their voltage control. The most general definition of voltage stability is that the power system is voltage stable if it is able to maintain steady acceptable voltages at all buses in the system in its normal operating state and after being subjected to a disturbance. Voltage instability is the absence of voltage stability resulting in the progressive decrease (or increase) in voltage. It normally involves large disturbances such as rapid load increase or power transfers. Alike to the case of rotor angle stability, the state of the system can be classified into categories such as small-disturbance voltage stable, voltage stable, and undergoing voltage collapse [23]. The system is considered to be small-disturbance voltage stable if following a small disturbance the voltages near loads remain identical or close to the pre-disturbance values. A power system is voltage stable if after subjected to a disturbance the voltages near loads reach post-disturbance equilibrium values. Finally, the system is undergoing voltage collapse if post-disturbance voltages are below acceptable limits and are progressively approaching lower values. 10 Voltage stability phenomenon can be separated into two time frames, transient and long term. In the case of transient voltage instability, voltage collapses quickly, within the time period of ten seconds. This collapse is usually caused by fast acting load components such as induction motors and dc converters. It is often difficult to distinguish between transient rotor angle instability and transient voltage instability as both phenomena may co-exist. In his book [23] Taylor raises an important question: "Does voltage collapse cause loss of synchronism, or does loss of synchronism cause voltage collapse?" Long-term voltage instability usually involves high loads, high power imports from remote generation and a sudden large disturbance. Instability can be driven by tap changers and distribution voltage regulators that try to restore the load power levels therefore causing further decrease in voltage. The instability may also evolve over the longer period of time caused by large load buildup or increase in power transfers. 2.4 Traditional approaches to stability studies 2.4.1 Rotor angle stability analyses 2.4.1.1 Modal analysis for small signal stability evaluation Modal analysis with eigenvalue approach can be used to analyze small signal rotor angle stability problems. As previously mentioned, if the disturbance is small enough it is possible to linearize the non-linear network equations around the operating point using for example Taylor's series expansion. The linearized form of the system equations can be expressed as: 11 A[x] = [A][Ax] + [B][Au] A[y] = [C][Ax] + [D][Au] (2.4) where [Ax], [Ay] and [Au] are the state, output and input vectors respectively, and [A], [B], [C] and [D] are the state, input, output and feedforward matrices as explained in [15]. Eigenvalues of the state matrix [A], better known as the Jacobian matrix, can give insight into the small signal stability characteristics of the system. For example, if the eigenvalues have negative real parts the system is asymptotically stable, and if at least one of the eigenvalues has the positive real part the system is unstable. Eigenvectors associated with eigenvalues can provide information on relative activity of the state variables Axt when particular "mode" is excited. In other words it can give us an idea what elements in the system are participating in different modes of the system response to a disturbance. 2.4.1.2 Equal area criterion It is not essential to solve the swing equation to determine if the response of a rotor angle to a system disturbance increases indefinitely or oscillates around the equilibrium point. The fast screening of critical contingencies for stability evaluation can be done by transforming the multi-machine power system parameters into a one-machine infinite bus system [25], and applying the equal area criterion method. This method has traditionally been used to demonstrate the basic concepts of rotor angle transient instability in power systems. 12 m \ prefault " "^N^X*^^^ postfault / LXy^ \ \ during fault • • fault clearing Figure 2. Illustration of equal area criterion The simple idea behind this concept is that for the rotor angle to remain stable kinetic energy Ej gained during the rotor acceleration due to a fault has to be lost during deceleration after the fault is cleared (£2). Gained and lost kinetic energy can be expressed as "c7 E] = area Ax = j(Pm - Pe) dS s„, E2 = area A2 = J(Pe -P m)dS (2.5) where Pm and Pe are mechanical and electrical power output, and 8 is the rotor angle. If the clearing time is delayed, gained kinetic energy will be greater than the one transferred back to the system, which will result in the continuous increase of the rotor angle and the loss of synchronism. 13 2.4.1.3 Time-domain simulation for evaluation of transient rotor angle stability The widely used method of analyzing transient rotor angle stability is step-by-step time-domain simulation in which the nonlinear differential equations, i.e. the swing equations, are solved using numerical integration techniques. The system of differential equations with known initial values can be written formally as = / ( M , 0 (2-6) dt where [x] is the state vector of n independent variables, and t is the time. Integration techniques that can be used are numerous and well described in the literature, therefore they will not be discussed here. In the traditional step-by-step method the three-phase network is represented by its single-phase equivalent (balanced operation of the system is assumed). Unbalanced conditions such as the case of a single line-to-ground fault are represented by positive, negative and zero symmetrical components and appropriate connection of the sequence networks [15]. Difference and algebraic equations representing the state of the system prior to, during, and after a considered disturbance are solved step-by-step giving the time response of the calculated variables. It is important to emphasize here that this method works in the frequency domain where all the time varying quantities are represented with phasors, and it differs from the full time-domain solutions provided for instance by the EMTP program. 14 2.4.2 Voltage stability analyses 2.4.2.1 Power flow methods, P- V and Q- V curves Power flow methods are traditionally most widely used when dealing with voltage stability. The power flow problem basically solves the complex matrix equation (2.7): YV = I = — (2.7) V* where Y is the nodal admittance matrix, V is the unknown complex node voltage vector, I is the nodal current injection vector and S=P+jQ is the apparent power nodal injection vector representing specified load and generation at nodes. In addition, equations for area interchange control, generator reactive power limits, bus voltage control, tap changer control, HVDC links, static var compensators, etc. are solved. The main methods for solving the power flow problem are the Newton-Raphson method and fast decoupled methods discussed elsewhere [9] ,[23]. Long-term forms of voltage instability phenomena can be analyzed as a steady-state problem using power flow methods. Simulations of "snapshots" in time are performed for the power system after the disturbance or during the load buildup. Apart from these so called "post-disturbance" power flows, there are two other power flow methods widely used: P-V and V-Q curves. These methods are used to obtain a steady state loadability limits which are closely related to voltage stability. P-V curves are useful for conceptual analysis of voltage stability. P can be a total load in a studied area where V is the voltage at some representative bus. P can also be a power transfer across the transmission ties between the 15 systems. P-V curves are mostly useful to analyze the load characteristics as a function of voltage. For example, an active power of a resistive load can be represented by a V 2 / R curve, while the constant power load is a vertical line in a P-V coordinate system. As area load or power transfer are increased, P-V curves for some representative bus may be obtained from the power flow solutions for each operating state. The power will reach its maximum point and the voltage at that particular point is called the critical voltage. Critical voltage is an important aspect for voltage security assessment. At points along the P-V curve, voltage stability can be assessed by linearized system sensitivity analysis or steady-state modal analysis. The disadvantage of this power flow method is that the solution will diverge near this maximum power point on a P-V curve. Another way of analyzing voltage instability is by V-Q curves. These are similar to the P-V curves and are used to test the system robustness by adding reactive load, since reactive load and voltage performance of the system are closely related. V-Q curves describe the behaviour of voltage versus reactive power at the test bus. The bus is connected to a fictitious synchronous condenser, and power flow solutions are performed to obtain the reactive power of a condenser for a set of voltage schedules. The advantage of V-Q curve analysis is that it gives the reactive power margin at the test bus, which is an important aspect used in voltage security analysis. The artificial P-V bus also minimizes power flow divergence problems. The problem with V-Q analysis is that it is generally not known at which buses the curves should be generated. By focusing on a small number of 2 also called nose, maximum power point 16 buses, system wide problems may not be recognized. However, many utilities still use V-Q curves as the base method for voltage stability analysis. In the Western Systems Coordinating Council (WSCC) each member is required to conduct P-V as well as Q-V analysis to ensure that the minimum required margins are met [3]. Each analysis is needed to confirm the results of the other (for example, P-V analysis is needed to confirm the results of Q-V analysis and vice versa). Recent research in the area of voltage stability static analysis is directed to improving the convergence of the power flow computations close to the maximum power point when generating P-V curves. The problem can be overcome by adjusting the power step size. Other methods of importance include the continuation method [4], optimal power flow techniques and combinations of power flow methods with artificial intelligence. 2.4.2.2 Modal analysis for voltage stability evaluation Modal analysis or eigenvalue analysis for voltage stability application is another approach greatly discussed in the literature, i.e. [15] and [12], and used by some utilities. Modal analysis is very similar to sensitivity analysis. Sensitivity analysis in voltage stability studies looks at the sensitivity of complex voltages to small changes in real and reactive power using the linearized model of the network around the solution point. Sensitivity analysis is normally utilized for design of voltage controls and reactive power compensation. Different sensitivity formulations can be used to study the effects of other controls. In modal analysis for voltage stability evaluation eigenvalues and the associated eigenvectors of reduced Jacobian matrices [12] are computed. The 17 eigenvalues of the Jacobian are used to identify different modes through which the system could become unstable. If all eigenvalues are positive the system is considered to be voltage stable. The sign of the eigenvalues criteria is opposite from the one used in small signal rotor angle stability and is related to the fact that V-Q sensitivity on each bus has to be positive for the system to be stable (higher reactive power injection, higher voltage). The magnitude of the eigenvalues provides a relative measure of how close the system is to instability and eigenvectors provide information related to the mechanisms of voltage instability. Because of the non-linearity of the problem, the absolute stability margin cannot be determined. Instead the system is stressed incrementally until it becomes unstable and modal analysis is applied to each operating point to determine critical areas and participating factors to voltage instability following a major disturbance in the system. 2.4.2.3 Time-domain simulation in voltage stability studies The most recent research trends in voltage stability analyses are moving away from the pure static traditional models. Highly loaded power network conditions require the consideration of voltage instability as a dynamic problem for a more adequate and accurate assessment. Dynamic analyses include time-domain simulations and linearized eigenvalue analysis. While power flow analyses may be adequate for simulation of slower forms of voltage instability, dynamic simulation is essential for example when the system operation is close to stability limits or when the fast dynamic process associated with voltage collapse phenomena is analyzed. Time-domain methods also offer higher modelling accuracy, and therefore the ability to study instability mechanisms not captured by static analysis. Results obtained by time-domain simulations are of higher value, for 18 example they contain information on the sequence of events leading to voltage instability or collapse. These methods are often referred to as multi-time-scale simulations. Time-domain simulation of voltage phenomena requires numerical integration of a large set of differential-algebraic continuous-discrete time equations using different integration methods. This is in general very demanding in computation time and gathering of data. Recent work in this area includes computational solutions such as parallel processing. Short-term voltage stability simulations use the same models developed for transient rotor angle stability studies with proper account for load behaviour. Long-term voltage stability simulations are more ambiguous as the problem of model stiffness comes into account. The dynamics of some parts of the instability process are very fast when compared to the interval over which the whole process is studied. The common approach to this problem is variation of the time step size. 19 Chapter 3 POWER SYSTEM SECURITY ASSESSMENT Security assessment is a more general term that covers all forms of stability analyses discussed in the previous chapter. The main aspect of security evaluation of a power system is to compute the security limits and determine the transmission capacity available in the network above already committed uses. This is found to be a challenging task in already highly loaded transmission networks. Dynamic tools for stability evaluation discussed earlier may offer the solution to high requirements for accuracy of those computations. This chapter provides definitions and criteria for power system security assessment. Issues related to the available transfer capabilities (ATC) determination are presented and the example of A T C determination and voltage security assessment of a test system using static power flow analysis is given at the end of the chapter. 3.1 Definitions and Criteria Security analysis is concerned with the system ability to withstand the selected number of "credible" contingencies. The system operating states for the purpose of security evaluation are classified in at least three operating states: normal operating state, emergency (alert) operating state and a restorative state. If load demand is met and all the operating constraints are satisfied the system is in a normal state. After a disturbance, the system may 20 settle into a new normal state, or it may enter an emergency operating state where some of the constraints are not met. Finally, in the restorative state operating constraints are satisfied but not all of the load demand is met due to a blackout or a load shedding action. It is convenient to distinguish between static and dynamic emergency states, or static and dynamic security. In a static emergency state the system reaches a new steady state where some of the operating constraints are violated, so the operation of the system is not permissible for an extended time period. The dynamic emergency state occurs when the disturbance causes the system to become unstable and operating and load constraints are not satisfied. Correspondingly, static security deals with the ability of the system to withstand disturbances without entering the static emergency state, while for dynamic security the system should not enter the dynamic emergency state. Stability is one aspect of dynamic security. The phenomenon is normally too fast for the operators to take actions and bring the system from emergency state back to normal state, therefore automatic controls are implemented. The objective of these security controls is to keep the system in normal state. Three levels of security assessment can be distinguished: security monitoring, security analysis and security margin determination. Security monitoring is checking for violation of operating constraints. Security analysis is checking the system ability to withstand disturbances. In practice, system security is checked for a set of contingencies with a reasonable probability of occurrence. Security margin determination is concerned with how far the system can move from the operating point and still remain secure. 21 One of the key factors in power systems operation, which determination is based on limits obtained by security analysis, is the available transfer capability (ATC) of a transmission network. The following section describes issues related to A T C determination. 3.2 Available Transfer Capability determination Available Transfer Capability (ATC), as defined by the North American Reliability Council (NERC), is the measure of the transfer capability remaining in the physical transmission network for further commercial activity over and above already committed uses. The units of transfer capability are in terms of electric power expressed in megawatts. The maximum power transfer is limited by the physical and electrical characteristics of the system such as thermal limits and voltage and rotor angle stability limits. Thermal overload limits are determined by the maximum amount of electrical current that a transmission line or electrical facility can conduct over a specified time period before it sustains permanent damage from overheating or before it violates public safety requirements. Violation of low voltage limits is encountered for the power transfer that causes at least one bus to reach the low voltage level defined by the operators. As explained in the previous chapter, transient stability limits are the measure of the transmission network capability to survive disturbances through transient and dynamic time periods (from milliseconds to several minutes) following a disturbance. Mathematically, A T C is defined as the Total Transfer Capability (TTC) less the Transmission Reliability Margin (TRM) less the sum of 22 Existing Transmission Commitments (ETC) and the Capacity Benefit Margin (CBM): A T C = T T C - T R M - E T C - C B M (3.1) Total Transfer Capability (TTC) is defined as the amount of electric power that can be transferred over the interconnected transmission network in a reliable manner while meeting all items in a specific set of defined pre- and post-contingency system conditions. Transmission Reliability Margin (TRM) is defined as the amount of transmission transfer capability necessary to ensure that the interconnected transmission network is secure under a reasonable range of uncertainties in system conditions. Capacity Benefit Margin (CBM) is defined as the amount of transmission transfer capability reserved by load serving entities to ensure access to generation from interconnected systems to meet generation reliability requirements. Existing Transmission Commitments (ETC) term essentially includes all normal (pre-transfer) transmission flows included in the given case. The ETC and C B M quantities can be further described in terms of their contractual firmness using terminology such as "recallable", "non-recallable" (or "firm" and "non-firm"), "scheduled" and "reserved" transfers. In an Energy Management System (EMS) the A T C program interfaces with the following modules: State Estimator (SE), Security Analysis (SA), Current Operating Plan (COP) and the OASIS, as shown iri Figure 3. 23 Figure 3. Framework for on-line A T C calculation The current system state is obtained from the State Estimator, the contingency list is obtained from the Security Analysis, and the Current Operating Plan is the source for load predictions, generation schedules, and outage equipment information. Calculated A T C values are transferred for posting at the internet site of the Open Access Same-Time Information System (OASIS). Typical values transferred to the OASIS include the interface identifier, the date and time of the run, the list of constraining facilities, TTC and A T C . The contractual nature of the transactions can influence the level of assurance needed for the transfer to take place. This means that more contingencies may have to be studied for more firm transactions. As the impact of ETC and C B M terms on the A T C quantity is unique for each case, methodologies for A T C determination are mainly focused on determination of T R M and TTC. Some of the major points in the concept for determining TTC 24 are: system conditions, critical contingencies, system limits, parallel path flows and non-simultaneous and simultaneous transfers. Base system conditions are identified and modelled for the period being analyzed, including projected customer demands, generation dispatch, system configuration and base scheduled transfers. The base system conditions may need to be modified as the system conditions change. Generation and transmission system contingencies throughout the network are evaluated to determine the most critical contingencies for the transfer being analyzed. The evaluation process should include a variety of system operating conditions. Concerning the system limits, as discussed earlier, the transfer capability of the transmission network may be limited by the physical and electrical characteristics of the system including thermal, voltage and stability considerations. Once the critical contingencies are identified, their impact on the network must be evaluated to determine the most restrictive of those limitations: TTC = Minimum of {Thermal Limit, Voltage Limit, Stability Limits} Parallel path flows can affect all systems of an interconnected network, especially those near the transacting systems. This may result in transmission limitations in systems other than the transacting systems, which in turn can limit the transfer capability between the two contracting areas. The term "simultaneous" refers to more than one transaction and "non-simultaneous" refers to a single transaction. No simple relationship exists between simultaneous and non-simultaneous transfer capabilities. The 25 simultaneous transfer capability may be lower than the sum of the individual non-simultaneous transfer capabilities. Total Transfer Capability is defined as the amount of electric power transferred in a reliable manner under the following conditions: • For the existing or planned system configuration and in normal (pre-contingency) operating conditions, all facility loadings are within normal ratings and all voltages are within normal limits. • The electric systems are capable of absorbing the dynamic swings and remaining stable following a disturbance that results in the loss of any single component such as transmission line, transformer, or generating unit. • After dynamic swings damp out following a disturbance and after the operation of any automatic system but before any post-contingency operator-initiated adjustments are implemented, all transmission facility loadings are within emergency limits and all voltages are within emergency limits. The procedure of computing TTC may proceed as follows: • Definition of a base case: It can be a current or forecasted condition, existing or planned configuration. It is important to specify areas that may include one or more generators. If there are more generators in the area, appropriate units dispatch must be specified for increase and decrease of outputs. Base cases for future hours are obtained from the Current Operating Plan (COP) which includes: (i) hourly load predictions, (ii) 26 commitment of generating units in the system and their incremental costs, (iii) amount and duration of all contractual agreements on an hourly basis and (iv) schedule of transmission and generation outages. This information is used to build power flow base cases to compute the transmission flows for facilities involved in the A T C computation. Specification of contingencies: The list of contingencies could vary from single outages such as a loss of line or a generator to complex switching scenarios. Determination of network response: A computer simulation is done to determine how the specified generation changes impact transmission lines flows, system voltages, and stability margins. This has to be done for the base case with normal limits enforced plus all specified contingencies with emergency limits enforced. Limitations checked for the base case are: normal branch flow limits, normal corridor transfer limits, bus voltage limits and voltage collapse condition. For contingency cases, the limitations checked are the emergency branch flow limits, emergency corridor transfer limits, bus voltage limits, bus voltage change limits and voltage collapse condition. Finding the maximum transfer: If the base case condition and configuration does not satisfy normal constraints on line flows, voltage and stability, some modification of parameters or conditions should be obtained to make the base case secure as a starting point. From that point, a systematic procedure of 27 increasing the specified transfer can be used to determine the maximum transfer that satisfies above criteria. • Repeat for alternate cases: Previous steps must be repeated for all possible cases that may be in effect at the time TTC is used. Since T R M is designed to account for uncertainty in the model configuration and operating conditions, the alternate cases could be used to compute TRM. The available capability for a given transfer may be smaller for the alternate cases considered. If they are weighted by the likelihood of occurrence, T R M could be determined for the base case. The T R M associated with the base case is determined from the most limiting TTC of the alternate cases. The challenges of A T C computation lie in the need to consider all likely base cases, all likely contingencies, and systematically compute the maximum transfer capability. In the operating environment A T C numbers are posted on a short-term basis, so the number of base cases should be smaller than for a planning environment. A T C numbers must be updated after every transaction, and full A T C calculation of large systems (e.g. 15000 buses) considering up to 7000 contingencies for 500 different transaction directions could take up to 24 hours even when linear methods are used. 3.2.1 Example: A T C computation and static voltage security analysis of the 30 bus test system To demonstrate issues related to the A T C determination a case of static security analysis using commercial power flow program for two bilateral transactions is presented. Analyzed power system model is the 30-bus test 28 system shown in Figure 4 consisted of three generators (200 MW, 150 MW and 80 M W units at buses 1, 2 and 5, respectively), three synchronous condensers (40 MVar, 24 MVar and 24 MVar at buses 8, 11 and 13, respectively), thirty five transmission lines, seven transformers and total load of 278 M W and 125 MVar. The original IEEE test system is modified by replacing the compensator at bus 5 with a generator and reducing the load at bus 30 by 50% to achieve better voltage profile as a starting point. Generators at buses 1, 2 and 5 are assumed to be operated by economic dispatch [27]. Figure 4. Modified IEEE 30-bus test system Two bilateral contracts between independent power producers (IPP) and their customers (LOAD) in amounts of 10 M W are considered. The power factor for both transactions is fixed at 0.95 lagging. In the base case all 29 voltages are within the chosen normal operational limits (0.95 and 1.05 p.u), and transmission lines are loaded less than 60% of their rated power as shown in Figure 5. 70.0 60.0 50.0 40.0 —\ 30.0 20.0 10.0 0.0 1. CN CN C N C M C M e M C M C N C N C M O i CN CN Line Figure 5. Transmission line congestion in the base case In the base case, the compensators at buses 11 and 13 are heavily loaded, working close to their reactive power generation limits. Reactive power generation of the generators is shown in Figure 6. 100 — 90 — 80 — 70 — 60 — 50 — 40 — r 30 — 20 — 10 — 0 ,— • 5 8 Bus Number Figure 6. Reactive power loading of the pool generation units in the base case 30 The TTC calculation for the bilateral transaction between buses 24 and 14 was performed by calculating the optimal power flow solution for increasing transfer with a constant power factor until the first violation of operational limits occurs. The algorithm is presented in Figure 7. Operational limits are chosen arbitrarily for bus voltages 0.95 and 1.05 p.u. and for transmission lines congestion at 100% of their rated apparent power. ( S T A R T ) /Transaction^ Y E S V not feasible Pre-contingency optimal power flow| for the next bilateral transaction Network congested or voltage limits violated? NO Optimal power flow for the next contingency with enforced YES z Security limits violated? NO All contingencies tested?. YES All bilateral cases; tested? NO NO YES pEND|§ Figure 7. Algorithm for T T C calculation For this bilateral case, the first violation is encountered in the voltage at bus 14, for the transfer of 350% of its contracted value. Transfer was then decreased until the violation is eliminated at 320% (32 MW). This amount of power transfer indicates the base TTC. Further on, the power flow is solved for n-1 transmission line outages to check for the emergency limits violations. Emergency limits of 0.9 and 1.1 31 p.u. are chosen for bus voltages and 150% of rated apparent power for transmission lines overload. The most limiting contingency is at the line between buses 12 and 14, and voltage violation is eliminated by decreasing the transfer to 201% (20.1 MW). This amount indicates the transmission security TTC, and because it is lower than the base TTC, it indicates the total transfer capability of the system for the given source-sink bus pair. Changes in voltage profile from the base case due to this transfer are depicted in Figure 8. 1.09 0.97 0.95 , , ,—, , , , , , , ,_, , , , ,_, * 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Bus Number Base Case 24-14 Figure 8. Changes in the voltage profile due to the bilateral transaction 24-14 The bilateral transaction is now scheduled between buses 17 and 27 in amount of 10 MW, coscp=0.95. Base TTC is determined by increasing the transferred power between buses 17 and 27 until the violation of operational limits is encountered at 250% of the contracted transfer value. Transfer is then 32 decreased to eliminate voltage violation at bus 30. Determined base TTC value is 213% (21.3 MW). Power flow is solved for n-1 transmission line outages to check for the emergency limit violation. As there were no limit violations, this transfer is limited by the normal operation voltage limits, not voltage security as in the first case. Therefore, the total transfer capability of the system for the given source and sink is 21.3 MW. The changes in the voltage profile are shown in Figure 9. 1.09 0.95 , — , , _ , r—, , , , , — , , _ , , — , . . , , r—, , — , ; — , 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Bus Number Base Case 17-27 Figure 9. Changes in the voltage profile due to the bilateral transaction 17-27 The bilateral transaction between buses 17 and 27 from the previous case is now scheduled while the transfer 24-14 is taking place. By similar procedure, TTC is determined to be limited by normal operation voltage limits and its amount is 23.4 MW. Changes in voltage profile for this case are depicted in Figure 10. 33 0.97 0.95 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Bus Number _ Base Case 24-14 & 17-27 Figure 10. Changes in the voltage profile due to bilateral transactions 24-14 and 17-27 From the comparison of TTC in case 2 and 3 we can notice that the value has increased for 2.1 M W in case when the transfer 24-14 is already taking place. IPP placed at bus 24 improves the base voltage profile of the system, therefore positively contributing to the voltage security. 34 Chapter 4 POWER SYSTEM MODELLING FOR STABILITY ANALYSIS USING T H E EMTP The Electromagnetic Transients Program EMTP is well known and has been used for the analysis of transients in electric networks caused by switching operations, lightning phenomena or different types of faults commonly encountered in the power system operation. Over the years detailed models of power systems components have been developed for that purpose in the EMTP. Use of integration techniques for solving nonlinear systems has been proven as a valuable tool for power systems analysis. EMTP methodology has tremendous potential in any type of studies involving electric systems. Stability studies of large power systems are not ordinarily done using a full time-domain solution. Traditional approaches rely mostly on 60 Hz frequency domain analyses, as it was noted in Chapter 2. In this chapter the classical model of a synchronous machine and static load model often used in stability studies are presented and their implementation in the EMTP is explained. The last section deals with the models for transmission lines and transformers that can be used for this type of studies carried out in the EMTP. 4.1 Synchronous machine modelling Synchronous generator characteristics and their limitations are crucial factors in all forms of instabilities. In rotor angle stability, a disturbance causes 35 an imbalance between the electromagnetic and mechanical torque which in turn causes the rotor angle to deviate from its original position thus leading to possible instability. In voltage stability, generators are primary sources of reactive power and they are to a great extent responsible for maintaining a good voltage profile throughout the system. As stated in [24], almost all voltage incidents have involved at least one generator operating with limited reactive power capabilities. In the following subsections the modelling of a synchronous generator is described and its implementation through a user defined subroutine in the EMTP is explained. 4.1.1 Classical model of a synchronous machine The most typical approach to synchronous machine modeling is to transform the stator equations using Park's transformation (dqO transformation). This method is well described in the literature and it will not be explained here. Since stability studies are mostly concerned with balanced system conditions, only direct (d) and quadrature (q) windings are considered, while the 0 winding, which plays the role in unbalanced conditions, is neglected. If currents and voltages in a three-phase system are assumed to be balanced, its equivalent can be represented with a single-phase network as it is mostly done in stability studies. This may be justified in the following way: instability phenomenon is mostly triggered by actions of protection devices, like opening of a line due to a fault. The fault itself may be unbalanced, such as a single line to ground fault, but as the breaker will open all three phases3 the system conditions can be considered balanced and single-phase 3 In some cases breakers are designed to open only a faulted phase, in order to minimize the supply disruption. 36 representation adequate. Conditions during the unbalanced fault are accounted for using appropriate connection of sequence networks. It is worth mentioning here (as pointed out in [10]) that this approach is questionable as negative and zero sequence currents flowing in the generator are not included in analysis. The simplified model of a generator used in stability studies and often referred to as a classical generator model is derived from a detailed machine model by making some assumptions. By neglecting the changes in the flux linkage in the field circuit, assuming that damper winding currents have died out and neglecting the saliency, a synchronous machine can be represented as a voltage source E' behind the transient reactance X/ (Figure 11). The classical model of a synchronous machine was implemented in the C++ programming language using the CONNEC subroutine that will be described in the following chapter. At each time step the MicroTran version of EMTP calls CONNEC and passes the values of the Thevenin equivalent of the network R,hev and Vopen to the subroutine. The equivalent circuit can be depicted as follows: E, -o Figure 11. Simplified transient model of a synchronous machine 37 'geA) t^hev^k) -A/W m (^J) v o p e n(t k) Figure 12. Thevenin equivalent of the network as seen from the C O N N E C subroutine Where the generator terminal current 7, is (4.1) The current returned from the subroutine to the main program is calculated at each discrete time point 4 from the internal generator voltage E'(t0 = Emax-cos(co-tk + 8) as follows: E\tk)-VopM (4.2) The initial values of internal generator voltage magnitude Emax and rotor angle 8 can be obtained from the steady state power flow solution. Emax is kept constant during the transient simulation, and the rotor angle is calculated using the swing equation. 38 4.1.2 The swing equation The fundamental equation in power system stability analysis that describes the effect of imbalance between the electromagnetic torque Te and the mechanical torque Tm of a synchronous machine was previously presented in (2.3) as: J— = T -T dt "' ' and dS_ dt = 6)-CD d2S dt2 da> dt (4.3) where 5 is the rotor angle (rad), co is the angular rotor speed (rad/s) and cos is the synchronous angular speed of the turbine-generator on the electrical side (2-7C-60 rad/s) dictated by the system frequency. The equations stated above are valid both for the mechanical and electrical sides. The relationship that converts mechanical to electrical quantities is where Se is the rotor angle in electrical radians, 8m is the rotor angle in mechanical radians and p is the number of poles. The rotor angle is measured with respect to the synchronously rotating reference frame. 39 In terms of power, the swing equation (3) can be rewritten in the following form: jda = K_K ( 4 5 ) dt co co In stability theory, the swing equation is often used in terms of power and in per unit form. Torque and power are used interchangeably (.Po^T^const., Pe=Te in per unit) as stated in [15]. This approximation is based on the assumption that in a perturbed system co is equal to cos when looking at the right side of equation (4.5), yet on the left side it is assumed to produce significant rotor angle deviation as — = co - cos. This assumption dt may not be sufficiently correct for certain types of analyses such as a full time-domain solution in the EMTP. In connection with reference [10], it is reasonable to take a closer look at the approximation stated above. Starting from the assumption that the mechanical power in a turbine will not change during the transient simulation4 and also with the assumption that the angular speed co=cosis used on the right side of the equation (4.5), it is implied that the term PJcos is a constant, and therefore Tm is a constant. That assumption of constant torque in a turbine, however, is not realistic. More correctly, torque/angular speed characteristics can be represented by the linear relationship: — ^mioperaling poinl) ~ ^Hmpie-speed ' (® ~ ^ S ' ) (4-6) 4 Constant as long as gate or valve posi t ion is not changed by a governor. 40 where Tm(operating p0mt) is the turbine torque at synchronous speed (Nm), and Dtorquespeed is the slope of the torque/angular speed curve at the operating point (J-s/rad2). If we replace the expression PJco in (4.5) with (4.6) and recall that Tm(0peratmgpoint) = PJOJS, the swing equation becomes: J ^ + D u , „ d i p - a ^ i B L . L (4.7) dt a>s co To be exact, we also have to analyze the electrical term Te=Pe/co. If we linearize the function of Te(co) around the operating point using Taylor's series expansion: ( rlT \ ' at cos with Te=Pe/o) Tt^-^(s) (4.9) cos cos The second term on the right-hand side of equation (4.9) introduces negative damping in the swing equation dt cos cos cos Approximation of neglecting that negative damping in the swing equation is possible if loads are not modelled as frequency dependant. Further analysis is continued with this assumption and the swing equation is written as 41 The damping constant D t o r q u e . s p e e d ranges from 0.5 to 2.0 p.u of the turbine rating for hydro turbines. Good agreement with the field tests has been obtained with the value of 1.5 p.u. [10]. As noted before, the swing equation is often used in per unit values. Although such approach is prone to misinterpretation of data, and physical units would perhaps be preferable, the present practice of most utilities is the analysis using the per unit system. Equation (4.11) can be written in per unit of active power by multiplying it with cos/Sbase'-J_coJLdaL+ Dl0„ea.-cos {(O_0)s) = InL__^_ (412) Instead of moment of inertia J, a more commonly used quantity is the inertia constant H (W-s/VA) defined as the kinetic energy at rated speed divided by the base power Sbase-H = -Jcol-— (4.13) Z ^hase Introducing (4.13) into (4.12) we get the familiar form of the swing equation used in stability studies: W - ^ B i f L ^ . ^ . p P ( 4. 1 4 ) cos at cos 42 Constants H and D(pMj are often given normalized to the machine rated power (Sralecj), which differs from the system base Sbase used in the study. Conversion can be done using base(H or D) = rated(H or D) • (4.15) Equation (4.14) in general form contains other damping effects than due to torque-speed characteristics. As a common practice is to represent all sources of damping in terms of damping torques proportional to speed deviations, the value of coefficient D can be chosen appropriately to account for these damping effects. 4.1.3 Numerical integration of the swing equation The next step in the modelling of a synchronous machine is applying a numerical integration technique to solve the nonlinear differential equation (4.14). In this work the implicit trapezoidal rule will be used because of its numerous advantages [17]. The application of the trapezoidal rule in solving electrical networks was first introduced in the EMTP [8] and over the years has proven to give reliable solutions. To demonstrate the use of the trapezoidal rule, the basic relationship of voltage and current for an inductance is considered: v,Xt) = L^P- (4.16) at Rewriting the equation and integrating both sides we get 43 (4.17) where ^ is a discrete time point, and At is a time step. The integral on the right hand side is trivial, but to calculate the integral on the left side we need to know the function vL(t) between points vL(t0 and vL(tk-At). The trapezoidal rule assumes that the function is a straight line connecting these points; therefore the integral is equal to = area = VA(O + V / . ( A -a 0 2 (4.18) v L(t) \ ( t k ) t k-At t Figure 13. Illustration of integration with the trapezoidal rule Finally for the inductance we have: 44 At e„ ('*) = - v t (tk - At) - ^ • iL (tk - At) At (4.19) A similar approach can be applied for the capacitance to obtain discrete time model of that circuit element, and the resistance does not need to be discretized. The swing equation can be rewritten in the form where Aco(t)=co(t)-a>s as: cos dt cos — = Aco (4.21) dt where Pa is the acceleration power (p.u.). The subscript (p.u.) in (4.20) will be omitted from now on, and all such quantities will be regarded as in per unit values on the base Sbase-Differential equation (4.20) resembles the well known differential equation describing an RL circuit: L ^ + R-i(t) = v(t) (4.22) dt where L corresponds to 2H/a>s, R to D/cos, i(t) to Aco(t) and v(t) to Pa(t). The equivalent circuits are shown in Figure 14. 45 R t i(t) v(t) O D/COg t Aco(t) (a) 2H/coc t KU R -VvV 2L At ^ e h ( t k ) D/cos - W r fAoo(t k) P A> 0 4H COgAt e h ( t k ) (b) Figure 14. Equivalent circuits for differential equations in (a) continuous and (b) discrete time representation By solving the difference equation represented by (b), the following expression for a change in angular velocity of a rotor is obtained: where Aoo(tk) = - - Aco(tk - At) + ±Pa(tk) + -Pa (tk - At) (4.23) a a a D AH D AH a = — H and b = cos cosAt 0)s cosAt In a similar manner, equation (4.21) integrated with the trapezoidal rule gives 46 (4.24) Equations (4.23) and (4.24) depend on values of electrical and mechanical variables from the previous time step, except for the acceleration power Pa(t0. The accelerating power is the difference between the mechanical power of a turbine, assumed to be constant during the transient simulation, and the electrical power drawn from the generator. To determine the electrical power output from the generator in the present time step some sort of a prediction scheme must be employed. In the next section, the calculation of active power in full time-domain simulations will be explained. 4.1.4 Calculation of active power in single-phase system When considering the power delivered by a generator to the system, it is sufficient to look at the instantaneous values of voltages and currents. Instantaneous power is defined as If v(t) and i(t) are sinusoidal functions of the same frequency, instantaneous power will be a sinusoidal function of double frequency with average value equal or greater than zero. This average value or dc component ofp(t) represents active power. When the three-phase electric power system is represented by its three phase model, calculation of electrical active power output from the generator is simple. By adding instantaneous power for all three phases in each time step, and as a result of the 120° displacement between them, the total representation using full time-domain solution p(t) = v(t)-i(t) (4.25) 47 instantaneous power of a three-phase system will be a constant equal to three times the average power per phase for balanced conditions. That electrical power is equal to the mechanical power supplied by a turbine minus associated losses, and mathematically it can be expressed as follows Pe=Pa(?) + Pb{t) + pc{t) (4.26) When the analysis of a three-phase system is done using single-phase steady state phasor representation active power in physical units is obtained by P«=l-V„m •/„,,. -cos6 (A21) where 6 is an angle between phase voltage Vrms and phase current Irms. When working in per unit system the above equation simply becomes ^ ( P , , ) = 4 , V ( , . , , ) - ^ , V ( / ) , , ) - C O S 6 ? (4.28) If a single-phase system representation is used with the analysis done in time-domain, the task of calculating active power is more challenging. The instantaneous power p(t) in one phase pulsates with twice the line frequency, therefore, its average value can be calculated using the well known expression valid for periodic functions Pe=]-)p(t)dt (4.29) o where T is a 120 Hz period. Pe at each instant 4 can be calculated as the average value of the function p(t) over the period [tk-T,tk]. When discretized with the trapezoidal rule, the above equation becomes the average of a sum of areas over the past period T or 48 I n = J V - l H = 0 (4.30) a r g a W = ^ - " ) + ^ - " " A ° - A r 2 where N=T/At is the number of samples within the interval T. T Figure 15. Calculating electrical power output Pe(tk) 4.1.5 Prediction scheme As observed in section 4.1.3, equation (4.23) requires prediction of electrical power Pe in the present time step 4. Pe(t0 can be obtained by the procedure explained in the previous section, and the unknown variable is the instantaneous powerp(t0=It(t^*E'(t^. That basically leads to the prediction of the rotor angle. Reference [8] suggests predicting angular speed using linear extrapolation, in the following form 49 A*V„ (tk) = 2-Aco(tk- At) -A(D(tk- 2 At) (4.31) Rotor angle is calculated with this predicted value using equation (4.24). With Spred, the internal voltage E'(t0 can be calculated, as well as the current Igen(tiJ, using equation (4.2). From there the electrical power Pe(t0 can be obtained and, therefore, the acceleration power Pa(t^=Pm-Pe(t^. The change in angular speed and rotor angle are then obtained from the swing equation (4.23) and (4.24). At this point predicted and calculated rotor angles can be compared and the procedure repeated if the error is larger that some threshold value. If the time step is sufficiently small, the iterative procedure may not be necessary. 4.1.6 Block diagram of a classical generator model implemented in C O N N E C At the end of the section on synchronous machine modelling, the procedure described in the previous subsections is summarised in the form of a block diagram. Figure 16 shows the block diagram of a classical generator model for stability studies implemented using the CONNEC subroutine in MicroTran. The case studies for testing the classical synchronous generator model using the examples from the literature are presented in Chapter 6, Sections 6.1 and 6.2. 50 yes Main program call to subroutine CONNEC at each time step \ Data passed from the main program R t h e v V operT A t For generator i read H, D, E m a x and S 0 from input file A(oi(tk)=AcDpred 8i(tk)-8pred Calculate generator power output P A ) P A ) = P m f P A ) Calculate AcO|(tk) and 8j(tk) from the swing equation Figure 16. Classical generator model block diagram 51 4.2 Load modelling Another important aspect of power system stability analysis is load modelling. Load characteristics have a significant effect on system performance. Results of power flow, voltage and rotor angle stability simulations are recognized to be highly dependant on the type of load characteristics assumed. Accurate modelling of loads is a difficult task due to several factors, mainly the large number of diverse load components, the lack of precise information on the composition of loads, and the lack of accurate procedures for load models identification. In stability studies frequency dependence of loads is seldom modelled. More frequently only voltage dependence of loads is addressed using so-called exponential load characteristics. In the following subsections the model of static exponential load is described and its implementation using the user defined subroutine CONNEC in MicroTran is explained. 4.2.1 Exponential load Since power system loads are composed of millions of different individual components, the common practice is to represent them as composite load characteristics, as seen from the bulk power delivery points. Load characteristic is an expression relating the load active and reactive power components as a function of voltage and an independent variable z that corresponds to the amount of connected equipment [24]. Therefore the (voltage) load characteristic can be formulated in general as: P = P(V,z) (4.32) Q = Q(V,z) 52 The widely used exponential load characteristic has the form: P = zPn Q = zQ0 (4.33) where VQ is the reference voltage, zPo and zQo are active and reactive power consumed at the reference voltage, and a and B are the exponents representing different types of load. There are three distinctive values for these exponents: • a=B=2 This is a constant impedance load model (often noted Z ) where the power varies directly with the square of voltage magnitude. At conditions below nominal voltage, this type of load will draw less current and vice versa, so its effect on the system stability is favourable. Residential loads for example have mainly constant impedance characteristics. • a=B=\ This is often referred to as a constant current load model (T). Here the load power varies directly with the voltage magnitude. In the absence of data this is the most acceptable load representation. • a=B=0 In the constant power load model (P) load power does not vary with changes in voltage magnitude. This type draws higher current under 53 low voltage conditions, and therefore it is very unfavourable for system stability. Induction motors have a tendency towards a constant power characteristic. Referring to [7], it is justified to use static load models in time-domain simulations for large system stability studies. The comparison of static and dynamic load model voltage responses to a three-phase fault has shown an agreement within 0.5% in the post disturbance voltages. Emphasis should be made here that the exponential load model is not valid for large voltage deviations. Reference [23] notes that this static load models are valid for voltage changes of ±10%. For composite loads exponents a may range from 0.5-1.8 and B between 1.5 and 6. In system wide studies the recommendation is to use a constant current (a=l) for active power load models and constant impedance (3=2) for reactive power load model. 4.2.2 Extracting time-varying RMS voltage response from the instantaneous values Equation (4.33) describes behavior of power systems loads in terms of their consumed active and reactive power dependance on RMS (effective) values of voltage. RMS values of power systems voltages and currents are 60-Hz-based steady-state quantities and when performing full time-domain solutions, they have to be somehow extracted from the instantaneous values. Since power system modelling in the EMTP accounts for transient conditions, voltages and currents during the simulation will in general contain higher harmonics. Those harmonics need to be filtered out, and only the 60-Hz component of the voltage is considered when accounting for the behavior of static exponential loads. 54 It should be noted here that the most accurate system representation would include all three phases with mutual couplings, plus dynamic load characteristics that would differentiate between fast and slow acting loads. However, for large system stability studies the approach explained here can be considered appropriate. A simple RLC parallel circuit in resonance may be used to filter out the 60 Hz component of voltage. Figure 17 depicts the circuit configuration for parallel resonance. Figure 17. Circuit configuration for filtering 60 H z component of VTN(t) The admittance of the circuit can be expressed as At resonant condition susceptance is equal to zero which leads to the well known expression for resonant frequency (4.34) 1 (4.35) 55 It can be noticed that the resonant frequency is entirely specified by the choice of LF and C>. Filtered voltage magnitude will vary with frequency of the input voltage as sketched in Figure 18. Figure 18. Variation of filtered voltage magnitude as a function of frequency Bandwidth of the resonant curve is noted as COBW- By forming the ratio of the resonant frequency to the bandwidth we obtain the factor which is a measure of the selectivity or sharpness of the tuning of the parallel RLC circuit. For a parallel resonant circuit the quality factor is ft = 3 _ = a)0RFCF (4.36) With the resonant frequency set to 60 Hz (377 rad/s) and with the desired quality factor, the values of RF, LF and CF can be determined from the above equations by freely choosing either one of them. Although it may seem 56 appropriate to choose a very high quality factor for high selectivity and narrow bandwidth, this filter would have a very slow response to sudden changes in the input voltage. Therefore a compromise has to be made when deciding about the appropriate values of filter circuit elements. For the purpose of filtering the 60 Hz component in power systems, where other frequencies are introduced mostly by switching manipulations and are of significantly higher value than rated frequency (few hundred Hz to kHz range), the bandwidth of 30 Hz was found to be an adequate choice. The circuit in Figure 17 can be discretized using the trapezoidal rule of integration in a similar manner as explained in the section on synchronous machine modelling. The equivalent circuit is depicted in the following figure. At At 2C C vF(tk) Figure 19. Equivalent circuit of the filter discretized with the trapezoidal rule where 4C, A^ vF(tk-At)-hc(tk-At) (4.37) 57 By solving the circuit, the filtered voltage is W ~ 'l At 2C„ ( 4 3 8 ) + + '-R,, 2LF At To calculate the RMS value of the filtered voltage a procedure similar to the one for calculating average power was applied. The RMS or effective value of a periodic function of voltage can be expressed as K,,=x^U0dt (4.39) x 0 When discretized with the trapezoidal rule the above integral becomes the sum of areas underneath the squared voltage. In a similar manner, we can assume that the RMS value of the voltage at each time step tk is the square root of the average of vF2(t) during the interval [tk-T, 7], or K,,M = *- Z area[n] (4.40) where N=T/At is the number of samples within the period T. Implemented as a MicroTran subroutine, the equivalent circuit for filtering the voltage at the point in the network where the load is connected is depicted in Figure 20. 58 VFOk) R t h e v ( t k ) V0pen(tk) Figure 20. Equivalent circuit for filtering the voltage as seen from the C O N N E C subroutine By similar analogy as in the case of calculation of active power, the value of v/yv has to be predicted in the present time step tk. Linear extrapolation of the following form is used: v w ^ ( ^ ) = 2 - v / A , ( ^ - A / ) - v w ( ^ - 2 A 0 (4.41) This value is then filtered using (4.38) to obtain VFpred(t0 and the RMS value of the voltage at tk can then be calculated using (4.40). 4.2.3 Representation of P Q loads in the E M T P Implementation of PQ loads in the EMTP can be realized using passive shunt admittances whose values are changing in accordance with changes in the load's active and reactive power consumption. 59 v r m s(t k) L O A D Figure 21. Load representation by passive admittances The equivalent admittance is equal to YL=G,+jBL= — + - 1 R, jcoL, (4.42) where R KM and L-YiM P Qco (4.43) and P and Q are calculated according to exponential load characteristic (4.33). If the load equivalent admittance is discretized with the trapezoidal rule the following equivalent circuit is obtained as seen from the CONNEC subroutine where K. ifk) = -j- v m (h - A0+hi ft* - A 0 (4.44) 60 loaA) R t h e v ( t k ) (j) V open ( t k ) Figure 22. Equivalent circuit as seen from the C O N N E C subroutine Finally, the load current Iioad(tk) returned to the main program at each time step is calculated as *load\h)~ D — K,hev (4.45) Vopen(tk) hu(tk) Rlh H—-——-r— ^ + + l J_ + AL RL 2LI. R, 2L, 4.2.4 Block diagram of an exponential load model implemented in C O N N E C The exponential load model described in the previous subsections was written in a C++ subroutine for the MicroTran version of the EMTP. The block diagram is depicted in Figure 23. The case study for testing the 61 exponential load model using the example from the literature is presented in Chapter 6, Section 6.3. Main program call to subroutine CONNEC at each time step tk Data passed from the main program Rthev Voper> A t For load i read Po. Qo. vO. « and p from input file 1 Filter out 60 Hz component of the load voltage 1 Calculate rms value of the filtered voltage Calculate P and Q at V n n A ) Find equivalent R L and L, Calculate l l o a d ( i ) 1 R e t u m 'load(i) t o t h e main program Figure 23. Block diagram of an exponential load model 62 4.3 Transmission lines and transformers In the last section of this chapter on power system models for stability studies, some issues related to transmission lines and transformers are addressed. The EMTP has well developed models for transmission lines that can be used for this purpose. Traditional stability studies use single-phase 71-circuit representations. In this work a constant-parameter (CP) line model is used to overcome the limitations on the number of user defined elements connected through the CONNEC subroutine. In the future, some modifications of the main program are expected to simplify this procedure. Both models will be briefly described in the following subsections. Another important issue that has to be considered is three-phase switching in a single-phase time-domain system representation. At the end of the chapter, the modelling of transformers with off-nominal turns ratio is presented. 4.3.1 Transmission lines in stability studies Transmission lines are characterized by four parameters: series resistance due to conductor resistivity, shunt conductance due to leakage currents between the phases and ground, series inductance due to magnetic field surrounding the conductors and shunt capacitance due to electric field between conductors. In general these parameters cannot be treated as lumped elements, but rather as distributed parameters per length of the line. Derivation of transmission lines models starts from the wave propagation equations [17]. With certain assumptions valid for 60 Hz analysis and for short line lengths, transmission lines can be approximated with nominal 7t-circuits, as depicted in Figure 24. 63 z o- •o Y 2 Y 2 o- -o Figure 24. Nominal 71 equivalent circuit of a transmission line With Z = Z'-I = R'-l + jcoV-l (4.46) Y = Y'-l = G'-l + jcoC'-l where Z and F are the total series impedance and shunt admittance of the line. Parameters R \ L',G' and C' are in per unit length, and / is the length of the line. As mentioned in the introduction, transmission lines are normally modelled as 7t-circuit equivalents in stability studies. A more accurate model of a transmission line is a constant parameter line model which is based on an ideal line representation plus lumped losses and with parameters that are constant and independent of frequency. The ideal line model is described with quantities such as the characteristic impedance of the line Zc (Q) and the travelling time z (s). The equivalent circuit of an ideal line is shown in Figure 25. 64 W Z c ko — • YvV vk(tk) z c W 6 6 vm(t k) Figure 25. Dommel's equivalent circuit of a constant-parameter ideal line model as seen from the line ends where ekh(tk) = vm(tk-T) + Zc-illl(tk-T) ifk) = vk ('* - T) + Zc • ik (tk - T) (4.47) Z C = A | — , and r = l-yjL'-C The series resistance of the line is modelled as lumped and inserted into three places, as depicted in Figure 26. lumped distributed lumped distributed lumped R/4 L ' .C R/2 L ' .C R/4 Figure 26. CP line model The CP line model is adequate for analysis at a single frequency, which is the case in stability studies. A particular advantage of using the CP line model over the equivalent n circuit for the models developed in this work is the decoupling of the network equations. That means that conditions at one 65 end of the line are not seen at the other end before the travelling time x. The CONNEC subroutine used for development of the synchronous machine and load models presented in the previous sections was originally designated for user-defined nonlinear elements. By decoupling, each nonlinear element can be solved independently of each other, as every one of them has its own decoupled Thevenin equivalent. In time, this requirement of separation of nonlinear elements with travel time was overcome by introducing M-phase nonlinear models [8], with the maximum number of phases equal to six. For the models presented here this would be the limitation on the maximum number of loads and generators together if not separated by the travelling time. The CP line model, conveniently, introduces separation of these models by the travelling time, but as a disadvantage it also limits the time step size to values smaller than normally used in stability studies (since the time step cannot be larger than the travelling time). This can significantly prolong the needed simulation time which is an important issue especially for on-line security assessment. Modifications of the MicroTran main program could resolve this issue and are expected in the future. 4.3.2 Three-phase switching in single-phase time-domain system simulations The single-phase system representation in the EMTP presented here has some drawbacks that have to be carefully examined. In stability studies many approximations are made to keep the simulated system responses pure 60 Hz sinusoidal functions. That means all the switching transients are ignored, as well as possible dc offsets. To give an example of such a case, when looking at a three-phase short circuit at the terminals of a synchronous machine, the fault currents of each phase will consist of decaying 60 Hz 66 component as well as a dc component that decays in several cycles [15],[9]. In stability studies the practice is to neglect this decaying dc offset by neglecting the transformer voltage terms in the stator voltage equations. Decaying components of the fundamental 60 Hz are accounted for through the division of the study time in three distinct time periods: subtransient, transient and steady state. When dealing with the single-phase system equivalent in full time-domain simulations, similar approximations have to be made. The dc offset in the system currents and voltages due to a fault has to be avoided. To accomplish this, the fault can be conveniently applied at the moment when no dc-offset will appear. In a largely inductive network this will be the moment when the voltage of the network at the place where the fault is applied is at its maximum. The dc offset is a consequence of the nature of the network elements and it is such that the current cannot change instantaneously through an inductance. The issue of avoiding the dc offset implies a limitation of fault analysis using single-phase equivalents in the EMTP as compared to what is traditionally done in stability studies: stability of the system can be studied only for clearing times of a fault equal to a multiple of half of the 60 Hz period. 4.3.3 Representation of transformers in the E M T P for stability studies In stability studies three-phase transformers are represented with their single-phase equivalents. An equivalent circuit with the secondary side referred to the primary side is shown in Figure 27. 67 jX 2 ' R 2 -ideal transformer Figure 27. Equivalent circuit of a single-phase transformer The magnetizing branch is normally neglected in power flow and stability studies. If the ratio of voltages V//V2 is equal to the ratio of base voltages for p.u. based values, the ideal transformer disappears from the above model. That is often not the case and the off-nominal turns ratio has to be included. Data for transformers provided by electric utilities for power flow and stability studies are usually given in terms of per unit short circuit impedance ZpM. and off-nominal turns ratio t. Z= R,+R2 ) + j\X,+X2 ^ p.u. ~ Z ' 2 ~ Rpu. + p.u. V, l base (4.48) V, t = l p.u. 2 p.u. (4.49) The single-phase transformer can be represented by its branch matrices [Rp.u.] and [Xp.u.]'' [10]. 68 R 1 p.u. 0 R 2 p.u. (4.50) M - • p.u. 1 X X p.U. 1 .e p.u. X p.u. where Rip.u-R2p.u.=Rp.u/2 is often assumed. The equivalent Tt-circuit is depicted in Figure 28 where Y = — (4.51) p.u. tY p.u. V 1base V. 2base Figure 28. Equivalent 7i-circuit of a transformer with off-nominal turns ratio t This transformer model can easily be implemented in the EMTP using the option for multiphase 7t-circuits and inductively coupled branches by entering matrices [Rp.u.] and [Xp.u.]"1 in "Linear and nonlinear elements" card with option "INVERSE". 69 The CP line model described in this section was used in the case studies for testing the classical synchronous machine model and the exponential load model presented in Chapter 6, sections 6.1, 6.2, and 6.3. Single-phase equivalent 7t-circuits for modelling transmission lines and transformer models with off-nominal turns ratio were used for a large power system simulation in the EMTP presented in section 6.4. This concludes the chapter on power system modelling for stability studies using the EMTP. 70 Chapter 5 ADAPTATION OF THE EMTP FOR ON-LINE STABILITY SIMULATIONS The "on-line simulations" term implies the possibility for the user to interact with running simulations. For power system analysis that primarily means manipulation with switches (or breakers). Stability analysis involves studying a large number of contingencies such as opening of transmission lines or loss of generators. This interactive type of analysis is not feasible with the standard off-line EMTP program. The need for an interactive simulation interface has not been essential for the type of studies traditionally carried out with the EMTP, but the concept can be very useful for many other purposes, including stability studies. In this chapter two subroutines available in MicroTran that can be used to implement this interactivity are described, and their use for adapting the EMTP for interactive on-line stability studies is explained. The subsequent sections present on-line switches manipulation through the Dynamic Data Exchange in Microsoft Windows® and user interface for simulation monitoring and switches manipulation. 71 5.1 User-supplied subroutines in the EMTP version MicroTran There are two user-supplied subroutines in MicroTran which are extensively used in this work: the ALPHA subroutine and the CONNEC subroutine. Both of them are compiled and linked with the main program through a dynamic link library (DLL), and executed at each time step of the simulation. The following subsections describe the main features of these subroutines and the adjustments made for the purpose of on-line simulation. 5.1.1 A L P H A subroutine The ALPHA subroutine interacts with the main program through the data card used to represent thyristors in the MicroTran input file and it communicates with the main program through its argument list. Its original application is for user modification of firing angles. A thyristor is in essence a controllable diode that allows current to flow only in one direction when a pulse is applied to the gate. A simple schematic is depicted in Figure 29. For a thyristor to conduct, voltage v k has to be higher than v m when a pulse is present at the gate. p switch v k ~2~ o VvV anode Figure 29. Representation of a thyristor in MicroTran gate | Kswitc • > f — - " A 2. swit h 'km catode 72 The thyristor model was adapted to serve as a controllable switch. Setting the argument KSPEC=-5 for a particular thyristor which ID number is different than zero in the input file, the current direction during conduction is ignored as well as the voltage polarity and the pulse width. In other words the thyristor will start conducting at instant TFIRE defined by the user regardless of the circuit condition. Another parameter of importance is CUROFF which forces the thyristor to open at any time if the current margin is set to a very high number. With these parameters in effect, the thyristor becomes a switch and it is completely controllable from within the subroutine. The main program can also pass to the ALPHA subroutine the variables and parameters that the user requests in the input file. In the request lines the user can ask for node and branch voltages, branch currents or instantaneous power, and synchronous machine parameters. These parameters are passed to the subroutine ALPHA in the array TACS at each time step, in the order requested by the user. Request lines may be impractical in the case of a large system with many variables of interest, and it is suggested that the EMTP code should be modified to pass to the subroutine all the variables calculated internally by the program and make them accessible to the user. The ALPHA subroutine and the MicroTran code are written in the FORTRAN programming language. Even though FORTRAN is still used, and new versions of this compiler are available on the market (e.g., Digital FORTRAN), C and C++ are by far currently the most recognized languages for software development. The C++ language was therefore used to write the subroutine, and to interact with other modules developed here. The issue that had to be resolved was the interfacing of two programs written in different languages. The way this was done was to modify the 73 existing subroutine code in FORTRAN to call A subroutine in C/C++. Special considerations had to be made for calling conventions that determine how the FORTRAN program makes a call to a subroutine, how the arguments are passed and how the routines are named. In a mixed-language program, different languages cannot share the same header files. As a result, if FORTRAN and C/C++ routines that use different calling conventions are linked, the error will not be apparent until the bad call is made at run-time. During execution the bad call causes indeterminate results and/or a fatal error, often somewhere in the program in a place that has no evident relation to the actual cause: memory/stack corruption due to calling errors. The caller routine uses a calling convention to determine the order in which to pass arguments to another routine. The called routine uses a calling convention to determine the order in which to receive arguments passed to it. These conventions in FORTRAN can be specified using the INTERFACE statement. The calling procedure used in this work is presented in Figure 30. This code declares a subroutine named ALPHAC with the C property and the external name _ A L P H A C set with the ALIAS property. The ALIAS directive is useful to establish naming convention as C/C++ are case sensitive and FORTRAN is not. The FORTRAN program uses the C attribute to call the C++ subroutine with the correct calling convention. This calling convention is used to pass all data by value except arrays. Arrays can only be passed by reference (as pointers). By specifying argument options as V A L U E or REFERENCE, arguments will be passed by value or by reference regardless of the calling convention. After adjusting the calling conventions, subroutine 74 A L P H A C in the C++ language is called from FORTRAN code by a C A L L statement. INTERFACE TO SUBROUTINE ALPHAC [C,ALIAS:'_ALPHAC] (TACS.KTACS, 1 INDEX,LINDEX,TFIRE,T,DELT AT,THYPRD,TON,KPOS.CUROFF, 2 ACLOSE,AOPEN,IDNO,ISTART,IWHERE,KSPEC,IHALF,IREAD, 3 LU5.LU6) REAL*8 TACS[REFERENCE],TFIRE[REFERENCE],T[VALUE],DELTAT[VALUE], 1 THYPRD[REFERENCE],TON[REFERENCE],CUROFF[REFERENCE], 2 ACLOSE[REFERENCE],AOPEN[REFERENCE] INTEGERS KTACS[VALUE],INDEX[REFERENCE],LINDEX[VALUE], 1 KPOS[REFERENCE],IDNO[REFERENCE],ISTART[REFERENCE], 2 IWHERE[REFERENCE],KSPEC[REFERENCE],IHALF[VALUE], 3 IREAD[VALUE],LU5[VALUE],LU6[VALUE] DIMENSION TACS(*),INDEX(*),TFIRE(*),THYPRD(*),TON(*),KPOS(*), 1 CUROFFn,KSPECn,ACLOSEn,AOPENn,IWHEREn,IDNO(*) , ISTART(*) END CALL ALPHAC (TACS.KTACS, 1 INDEX,LI NDEX,TFIRE,T,DELT AT,THYPRD,TON,KPOS.CUROFF, 2 ACLOSE,AOPEN,IDNO,ISTART,IWHERE,KSPEC,IHALF,IREAD, 3 LU5.LU6) END Figure 30. Calling C / C + + subroutine (ALPHAC) from the F O R T R A N code The diagram of the interaction between different modules is depicted in Figure 31. MicroTran passes a set of arguments to the ALPHA subroutine at each time step and receives back the position of the switches (open or closed). The FORTRAN subroutine ALPHA serves only to call ALPHAC. With modifications of the MicroTran code the call could be made directly to the A L P H A C subroutine. The programs connec.for and alpha.c are compiled together as a dynamic link library connec.dll that is executed at each time step of the simulation. 75 cbnnec.dll mt32 exe (main program) connector ^ alphac.c w (fortran subroutine (C++ subroutine ALPHA) Ijli tMffjpl ALPHAC) Figure 31. Block diagram of the interaction of subroutine modules with the main program through a dynamic link library 5.1.2 C O N N E C subroutine The CONNEC subroutine interacts with the main program through the input card used to represent nonlinear elements with the compensation method. In compensation-based methods [8], nonlinear elements are simulated as current injections super-imposed on the linear network after a solution without nonlinear elements is found. In other words, the network without nonlinear elements is solved and reduced to a Thevenin equivalent circuit between nodes k and m (Figure 32). ^thev y i \ Linear part of openl I the network 'km m (i) 'km Figure 32. Nonlinear element connected to a linear network In general for the multi-phase Thevenin equivalent circuit we can write 76 (5.1) The subroutine CONNEC communicates with the main program through its argument list. Arguments RTHEV and VOPEN are passed to the subroutine at each time step. RTHEV changes with switches position, and VOPEN changes from step to step. CONNEC returns the parameter CURR=ikm calculated internally within our subroutine. Our subroutine was implemented in C++ in a similar manner as the ALPHA subroutine. It has been used to develop the synchronous machine model and the exponential load model explained in the previous chapter. CONNEC and ALPHA subroutines are both compiled with the FORTRAN subroutine connec.for and linked to the main program through a dynamic link library (DLL). 5.2 On-line switches manipulation through Dynamic Data Exchange in Microsoft Windows This section will give an overview of the Dynamic Data Exchange in Microsoft Windows and its application to on-line switches manipulation in the EMTP. 5.2.1 Dynamic Data Exchange Dynamic Data Exchange (DDE) is a feature of Windows that allows two programs to share data or send commands directly to each other. DDE can be thought of as a direct conversation between two application programs. In most cases, one application is providing some form of data to another 77 application. The application that is the source of the data is called the "server" and the application that is receiving the data is called the "client", although the exchange of data in the opposite direction is also possible. Each data item that a server application can provide has a unique identifier consisting of three parts, a DDE Application Name, a DDE Topic, and a DDE Item Name. The DDE Application Name is almost always the executable filename for the server application. The DDE Topic typically identifies a group or category of data in the server application and each data item that a server can provide has a unique DDE Item Name. Thus, the Application Name, Topic, and Item Name identify the exact source of the data in a server application that is to be linked. DDE links are always initiated in the client application. The client initiates a DDE link by broadcasting a message containing a DDE Application Name, a DDE Topic, and optionally a DDE Item to all other applications currently running. If a server application that can provide the data is running, it responds to the "DDE initiate" and the Windows operating system opens a link between the two applications. DDE also allows a client application to send commands to a server. Either application involved in a DDE conversation can terminate the link. Closing either one of the linked applications causes all links between the two programs to be terminated. 5.2.2 Implementation of the D D E for on-line switches manipulation in the E M T P The DDE programs communication link was implemented using LabWindows/CVI. The ALPHA subroutine that controls the switches is assigned the role of a client, and the user interface program created in 78 LabWindows/CVI represents a server. LabWindows/CVI has an extensive library for DDE that can also be used in Microsoft Visual C++. The block diagram of the communication links between the programs is shown in Figure 33. USER INTERFACE (server) lositions croTran monitored variables lositions run Mi monitored variables switch F i argument list ALPHA MICROTRAN SUBROUTINE switch positions (client) Figure 33. Block diagram of the communication between the programs By launching MicroTran from within the interface, a DDE server is established and ready to receive the connection from a client. After the simulation has started MicroTran calls the subroutines at each time step. At the first call, the ALPHA subroutine registers as a client (only once during the simulation). The link between the server and the client is then established by a unique item name and a topic. DDE callback functions (server and client) provide mechanisms for sending and receiving data or commands from each other. A callback function in a client application can respond to only two types of DDE messages: D D E D I S CONNECT and D D E D A T A R E A D Y . With the hot link 79 established, this callback function is called with the D D E D A T A R E A D Y message whenever the data values change in the server application. A DDE callback function used in the server application can be triggered in several ways. For example, whenever a client application attempts to connect to the server or requests information from the program, the server callback function is executed to process the request. The user can provide opening and closing times of all controllable switches from within the interface. Data are sent through the hot link5 to the ALPHA subroutine and from there to MicroTran. A block diagram of the server-client communication is shown in Figure 34. The term "hot link" implies that client is advised immediately whenever the data item changes at the server site. 80 SERVER register as DDE server launch MT server I callback if DDE_CONNECT connect client no I else if DDE_DISCONNECT eliminate connection else if DDE_DATAREADY receive monitored data and plot on the screen • DDE DATA READY DDE C O N N E C T DDE DATAREADY W MT call to ALPHA connect to DDE server CLIENT--.-til set up hot link for switches position pass monitored data from MT to server , i •' • i yes '^•i.client.Tcallback . I if DDE_DISCONNECT eliminate connection . i else if DDE_DATAREADY pass switches position from server to MT Figure 34. Block diagram of a server-client communication 81 5.3 User interface for on-line simulation monitoring and switches manipulation Once a DDE link is established between the EMTP subroutine and another application, all the variables calculated internally and requested by the user in MicroTran become available for monitoring. The advantage of this is that the user can run the simulation for an extended period of time, apply unlimited number of contingencies within the same simulation by manipulation of the switches, and possibly record interesting results for later analysis. This would be especially advantageous when working with large power system models. While building the interface for large power systems is a challenging task, all the features can be examined on a small test system such as a circuit consisted of a voltage source, controllable switch and a resistor (Figure 35). e(t)=1.41sin(cot)^^ R=2 Q Figure 35. Circuit for testing the switches manipulation and on-line monitoring A simple user interface was developed using LabWindows/CVI. Monitored variables are node and branch voltages and a branch current. These 82 variables are passed from the ALPHA subroutine to the user interface through a DDE link at each time step and plotted on the screen using a strip chart. In this particular example, time waves are plotted on the screen. The interface can be easily adapted to plot RMS values using the method described in section 4.2.2. The user defines opening and closing times of a switch by entering the desired values in the provided numeric fields. A hot link established between the interface and the ALPHA subroutine will immediately process the information and pass it to MicroTran. In case of a large system it may be impractical to enter the data manually and the list of opening and closing times for all the switches could be read for example from a separate data file. The test user interface is depicted in Figure 36. This section concludes the chapter on the adaptation of the EMTP for on-line stability studies. i'^ " HI'. 98 83294 " n , e Zoom Out I """Fast 148:99940 Slow Zoom In; Figure 36. Interface for switches manipulation and on-line simulation monitoring of a test system 83 Chapter 6 CASE STUDIES In this chapter the models developed in Chapter 4 are tested on three examples. Examples of a generator connected to a large power system (an infinite bus case) [15] and the WSCC three-generator nine-bus system [5] are used to check the validity of the synchronous machine model. A load connected to a large power system (infinite bus) [15] is used to test the exponential load model comparing different values of exponents a and p\ Since these examples are taken from the literature, it is possible to compare the results obtained here with those from traditional stability analysis, and verify the validity of the proposed methodology. The feasibility of a large power system simulation in the EMTP for power flow and stability studies is examined in the last section using the reduced WSCC system model. 6.1 Generator-infinite bus case A transient stability of a generator supplying power to an infinite bus through two transmission lines is studied. The system equivalent circuit is shown in Figure 37. 84 Figure 37. Equivalent circuit of a generator-infinite bus case The initial system operating conditions and generator parameters are given in per unit with Sbase= 2220 M V A and Vbase = 24 kV: P = Q.9 Q = -0.436 E, =1.0| 28.34° EB = 0 .90081[£ X'd=0.3 H = 3.5 MWslMVA D = 0 From initial conditions, voltage £ ' can be calculated using the following equations in phasor quantities: E' = E,+jXdll j_P+jQ Expressing the circuit variables in form of sinusoids we have: e \t) = yf2E 'cos(cot + 5) = \ .6442 cos(W + 41.77°) els(t) = yf2EH cos(a>t) = 1.2739cos(ft>r) 85 The generator-infinite bus case was built in MicroTran. The infinite bus is represented by an ideal voltage source eR.(t) and the generator is represented by the classical generator model developed in Chapter 4. A three-phase fault on a transmission line between buses 3 and 4 was applied, and the fault is cleared by opening breakers of the faulted line. Transient stability of a generator rotor angle is studied for different clearing times and the results are presented next. For a clearing time tc = 0.067 s and with a damping coefficient D=0, the system is first-swing stable, and becomes unstable in the third swing as depicted in Figure 38. This is in disagreement with the results obtained by step-by-step stability analysis presented in reference [15]. By introducing damping £>=1.5 p.u. due to the torque-speed characteristics, the response of the rotor angle becomes stable. The response of active power for both cases is depicted in Figure 39. Fault voltage and current are shown in Figure 40. It can be noticed that the fault current has a small dc-offset. This is because the switching at the exact moment when voltage is at its maximum is not possible due to a discretization. Figures 41-44 present responses for fault clearing times of7c=0.083 and/c-0.092. 86 Time (s) Figure 38. Rotor angle response to a disturbance with clearing time tc=0.067 s 1.5 1 0.5 & 0 CD 0--0.5 -1 -1.5 D-1.5 ^ \ ____Q=P_ \ _ V | 2 3 Time (s) Figure 39. Generator active power after the disturbance with clearing time tc=0.067 s 87 8.0 6.0 4.0 2.0 3 0.0 CX, v[3',gnd](1) i[3',gnd](1) -2.0 -4.0 -6.0 -8.0 -V 800.0 880.0 960.0 1,040.0 1,120.0 Time (ms) 1,200.0 Figure 40. Fault current and voltage for tc=0.067 s Figure 41. Rotor angle response to a disturbance with clearing time tc= 0.083 s 88 1.5 1 -0.5 s 0 I (D 0_ -0.5 -1 -1.5 -I 0 1 2 3 4 5 Time (s) Figure 42. Generator active power after the disturbance with clearing time tc=0.083 s 20 0 -I 0 1 2 3 4 5 Time (s) Figure 43. Rotor angle response to a disturbance with clearing time tc=0.092 s 89 1.5 Time (s) Figure 44. Generator active power after the disturbance with clearing time tc—0.092 s We can see in Figure 41 for rc=0.083 s and D=0 that the rotor angle is again stable for the first swing and becomes unstable in the second swing. Introduction of some damping gives stable responses. With classical stability programs this clearing time gives a stable response for D=0. Finally for the last clearing time tested, fc=0.092 s, the rotor angle is transient unstable in the first swing with or without added damping. This matches the results provided by the stability program. The question rises about the cause for the discrepancy in the results. While for the first swing all the presented results match results provided by traditional stability analysis, it is not apparent where the later instability comes from. A full time-domain analysis provides more accurate results involving all the transients. In classical stability theory transients are neglected in the 90 network as well as in the synchronous machine model by neglecting the transformer voltage terms in the stator voltage equations. The effect of angular speed variations due to a disturbance is neglected as well. Approximated models that work reasonably well in stability type of analysis may need to include different assumptions when a full time-domain solution is performed. Inclusion of a damping term in the swing equation may be one of them. Anderson and Fouad in their book [5] note that the classical generator model can give reasonably accurate results only during the first second after a disturbance. It is unclear why other references such as [15] give the results of analyses with the classical generator model for several seconds. If the above statement is accepted as valid, the classical generator model can only assess the instability of the system if this instability were to happen during the first rotor angle swing. The classical model of a synchronous machine tested in our approach can be considered to be valid even when the damping term in the swing equation is neglected. 6.2 Three-generator case The same model of a synchronous machine was tested on the bigger system depicted in Figure 45. Generator data for the three machines are given in Table 1. Load and transmission lines data are given in Table 2. 91 18 kV 230 kV - • L o a d C 2 3 0 k V 1 3 8 k V Load A 230 kV Load B - 4 m 16.5 kV 1 & Figure 45. Three-generator nine-bus system Generator 1 2 3 rated MVA 247.5 192 128 rated kV 16.5 18 13.8 power factor 1 0.85 0.85 type hydro steam steam speed r/min 180 3600 3600 Xd 0.146 0.8958 1.3125 Xd' 0.0608 0.1195 0.1813 Xq 0.0969 0.8645 1.2578 Xq' 0.0969 0.1969 0.25 XI (leakage) 0.0336 0.0521 0.0742 TdO' 8.96 6 5.98 TqO' 0 0.535 0.6 Stored energy at rated 2364 640 301 speed MWs rated H (inertia 9.55 2.586 2.35 constant) MWs/MVA Table 1. Generator data 92 Transformer X 1-4 jO.0576 2-7 jO.0625 3-9 J0.0586 Line Z B/2 4-5 0.010+J0.085 jO.088 5-7 0.032+J0.161 J0.153 7-8 0.0085+J0.072 J0.0745 8-9 0.0119+J0.1008 jO.1045 9-6 0.039+J0.170 J0.179 6-4 0.017+jO.092 J0.079 Load Y A 1.2610-j0.5044 B 0.8777-J0.2926 C 0.9690-J0.3391 Table 2. Transmission lines and load data Transmission lines are modelled as constant parameter lines for the reasons explained in section 4.3.1. Loads are modelled as equivalent shunt admittances, and for generators the classical model presented in 4.1.1 is used. Initial system conditions are obtained from the steady state solution in MicroTran. Internal voltages of generators 1, 2 and 3, and their initial angles are calculated in a similar manner as in the previous case. The disturbance initiating the transient is a three-phase fault occurring near bus 7 at the end of line 5-7. The fault is cleared in five cycles (0.083 s) by opening line 5-7. Figure 46 and Figure 48 depict relative rotor angles for three machines for the damping coefficient set to 0 and 1.5, respectively. It can be noticed that the rotor angles are not stable for D=0, although the response is considerably different than in the case of a generator and infinite bus (more of an oscillatory type). Since the classical model gives only the information on the first-swing stability of the system, it can be concluded that the system is first-swing stable. This agrees very well with the results presented in [5]. Rotor angle responses to a disturbance were also tested for damping Z>=1.5 p.u. 93 (based on generator ratings). Similar to the previous case, the responses become stable (Figure 48). Active power outputs of the three machines are shown in Figure 47 and Figure 49. Generator 2 is supplying the most power to the system prior to the disturbance (163 MW). As the fault location is at the bus 7, and the short circuit current from generator 2 to the fault flows through pure reactances, electrical power output of the generator is zero during the fault. A l l generators are accelerating during the fault, but generator 1 for its high inertia constant will accelerate the least. After the fault is cleared by opening the line 5-7, generators 2 and 3 will start to decelerate to transfer the kinetic energy gained during the fault back to the system. It can be noticed that generator 1 starts acting as a motor, absorbing active power from the network. 120 -20 J I 0 1 2 3 4 5 Time (s) Figure 46. Relative rotor angle differences for three machines (a) S2-St, (b) S3-S, for D=0 94 -1.5 2 3 4 T i m e (s) Figure 47. Electrical power outputs of generators 1, 2 and 3 for D—0 90 80 70 60 x> c 50 - 40 -ro i _ 30 -o o a: 20 10 -0 -/ \ / \ (b) ) 1 2 3 4 T i m e (s) Figure 48. Relative rotor angle differences for three machines (a) S2-dn (b) Sj-S, for D= 1.5 95 I -1.5 J _____ 1 0 1 2 3 4 5 Time (s) Figure 49. Electrical power outputs of generators 1,2 and 3 for D=1.5 6.3 Load connected to an infinite bus case The exponential load model developed in Chapter 4 was tested on the simple system shown in Figure 50. The initial operating state of the system is described in per unit values by P 0=0.8 e 0 =0.4 V0 = 0.8554 a0= 13.52° 96 1 /o° A B , _TYYV\_ Z=j0.25 t P+jQ Figure 50. Load connected to an infinite bus The load is modelled as a constant power exponential load {cc=B=Q). For the initial operating state, time varying sinusoidal voltages at nodes A and B are depicted in Figure 51. The RMS value of the load voltage over a longer simulation time is depicted in Figure 52. 1.0 ? 0.0 -1.0 \ *A V[B] A A A A / 1 V | 7 T J I T V \ l l V II 0.0 50.0 T ime (ms) 100.0 Figure 51. Bus voltages for the initial operating state with constant power load (a=6=0) 97 1 -0.9 -0.8 -0.7 -0.6 -ri. 0.5 -CO 0.4 -> 0.3 -0.2 -0.1 -0 -0.00 0.15 0.30 0.45 0.60 Time (s) 0.75 0.90 1.05 Figure 52. RMS value of the load voltage at bus B over a longer period of time for («= /?= 0) If power consumed by the load is changed to P=1.8430 and Q=0.\5, voltage at bus B collapses after about 0.6 s, as depicted in Figure 53. This is in agreement with the results provided in reference [15]. 98 5v[B](1) 1.0 200.0 400.0 600.0 800.0 1,000.0 Time (ms) Figure 53. Voltage collapse due to a high consumption of a constant power load (a=B= 0) at bus B The simulation was repeated for the same load level based on lp.u. voltage, but with the load represented as a constant current (a=3=\) and constant impedance (a=B=2) exponential load model. The RMS value of the voltage at bus B is shown in Figure 54. As expected, the constant power load is unfavourable for the voltage stability of a system. This model gives an approximate behaviour of induction motors that tend to draw higher currents to maintain their consumption level if the voltage is decreased. 99 1 -, 0.00 0.15 0.30 0.45 0.60 0.75 0.90 1.05 Time (s) Figure 54. RMS voltage at bus B for cases of (a) constant power (a=B= 0), (b) constant current («= /?= 1) and (c) constant impedance («= /?= 2) load model for power consumption of P= 1.8430 and Q = 0.15 at the voltage level of 1 p.u. The results of testing the exponential load model in full time-domain simulations has shown a good agreement with traditional stability analysis. This model is fully applicable for stability analysis in the EMTP. 6.4 Full time-domain simulation of a large power system To examine the feasibility of simulating large power systems for power flow and stability studies with the EMPT, a case of a reduced Western Systems Coordinating Council (WSCC) system was built for the MicroTran version of the EMTP. The power flow data for this system were obtained from BCHydro in the PTI (Power Technologies Inc.) raw format commonly used by utilities. A summary of this case is provided below: 100 Transmission lines 3521 Transformers 1167 Buses 3615 Generators 362 Line shunts 244 Switched shunts 132 Loads 1734 The schematic of the 500 kV transmission network is depicted in Figure 55. The system model includes voltage levels from 2 kV to 500 kV and covers the area of British Columbia, Washington, Oregon, Idaho and Montana. The rest of the WSCC system is represented by equivalent loads. 101 Figure 55. 500 kV transmission network of the reduced W S C C system The transmission lines are modelled as equivalent n circuits and transformers with their off-nominal turns ratio are modelled in the way explained in chapter 4.3.3. Generators are represented as ideal voltage sources, and loads as constant impedances. The comparison of the results from the EMTP full time-domain simulation after 0.3 s and 1.2 s of simulation time, and the steady state solution obtained with a power flow program [21] is shown in Figure 56 and Figure 57. It can be noticed that most of the transients have died out after 1.2 s, and the agreement with the power flow solution is within 0.05 p.u. 102 Comparison of bus voltages between full time-domain solution in the EMTP after 0.3 s of simulation time and a power flow steady-state solution 0.2 0.15 (p.u.) 1 ll O L L o C L •in — o-i to" * o r- -»1 uk WW} ^ — c i C O i O * o »- p <£ o < i cu i 3 — <\ 3 O C y w o cb IOJ o in o ( O I O C O ^ O ' - ' - O J iir o in ui r-