Analys is of Electromagnetic and Electromechanical Power System Transients w i t h D y n a m i c Phasors by S E B A S T I A N H E N S C H E L Abitur mit Latinum, Lily-Braun-Oberschule, Berlin, Germany, 1988 Diplom-Ingenieur, Technische Universitat Berlin, Germany, 1992 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Doctor of Philosophy in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Electrical & Computer Engineering) We accept this thesis as conforming to the required standard The University of British Columbia February 1999 © Sebastian Henschel, 1999 in presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract Over the last 50 years, digital simulation of electric power systems has become an integral part for planning, design and operation in the power industry. The number of possibilities with respect to the purpose of a study, investigated frequency ranges, etc. with, in the past, limited computer resources has resulted in a spectrum of simulation tools, designed to handle very specific tasks. Simplifying assumptions were often needed to facilitate such a simulation. Recent system failures and power outages, partly due to increasingly sensitive operating conditions, have created a demand for more comprehensive studies and more general sim-ulation tools that overcome former limitations. With regards to time-domain simulation, this demand has led to combining the areas of transient, mid-term and long-term stability. Confronted with concerns about black start and system restoration due to a global trend to a deregulated power market, several power utilities suggested to also include the area of electromagnetic transients. However, previously made assumptions as well as technical limitations complicate the implementation of this idea: Stability programs are based on the assumption that power transfer takes place at system frequency and are therefore unable to represent rapid elec-tromagnetic transients. Electromagnetic transients programs, on the other hand, are very accurate but use too small simulation step sizes for an efficient simulation of electromechan-ical transients. A new method for simulating both types of transient phenomena with complex signals and dynamic phasors is presented in this thesis. Whereas in previous work three-phase transformations had been used to accomplish this task, this new method is applied directly in the phase-domain and not restricted to balanced three-phase systems. Several numerical aspects such as an appropriate variable representation, integration method and a control mechanism for variable simulation step sizes have been addressed. Important modeling ii A B S T R A C T in constraints for transmission lines and synchronous machines have been thoroughly analysed and are considered in this work. Finally, several simulation examples are presented that compare the results of the new simulation program with the results obtained from the E M T P . These case studies demon-strate that the methods presented in this thesis are suited to combine the areas of electro-magnetic and electromechanical transients simulation. Contents Abst rac t i i Lis t of Tables v i i i L i s t of Figures ix Acknowledgement x i i 1 In t roduct ion 1 1.1 General Background .1 1.1.1 Electromagnetic Transients 1 1.1.2 Electromechanical Transients 3 1.2 Review of Unified Time-Domain Simulation 4 1.3 Thesis Motivation and Objective 8 1 Numerical Considerations 11 2 C o m p l e x Var iab le Representat ion 12 2.1 Evaluation of Existing Variable Representations : 12 2.1.1 Park's Transformation 13 2.1.2 Clarke's Transformation 15 2.1.3 Summary 16 2.2 Complex Signals and Dynamic Phasors 17 2.2.1 Real Signals and Phasors 18 2.2.2 Introduction to Complex Signals and Dynamic Phasors 19 2.3 Spectral Properties of Complex Signals 27 2.3.1 Spectral Averages 27 2.3.2 Instantaneous Frequency 29 iv C O N T E N T S I 2.3.3 Peculiarities Regarding Complex Signals 31 3 Numerical Integration 33 3.1 Integration with Complex Signals 33 3.1.1 Complex Backward Euler Algorithm 34 3.1.2 Complex Trapezoidal Algorithm 35 3.1.3 Relaxed Trapezoidal Algorithm 35 3.1.4 Relaxed Convolution Algorithm 37 3.2 Integration Errors 39 3.2.1 Accuracy 39 3.2.2 Local Errors 45 3.2.3 Global Errors 49 3.3 Numerical Stability 51 3.3.1 Backward Euler Algorithm 52 3.3.2 Trapezoidal Algorithm 53 3.3.3 Relaxed Trapezoidal Algorithm 54 3.3.4 Relaxed Convolution Algorithm 54 4 Variable Integration Step Sizes 55 4.1 Step Size Prediction 56 4.1.1 Maximum and Minimum Frequency 56 4.1.2 Magnitude Deviation 57 4.1.3 Simulation Mode Selection 58 4.1.4 Step Size Prediction 60 4.2 Prediction Reliability 61 4.2.1 Sampling Quality 62 4.2.2 Frequency Fluctuation 63 4.2.3 Simulation Mode 63 4.3 Step Change Efficiency 64 4.3.1 Step Ratio 65 4.3.2 Time Delay 66 4.4 Decision Making Process 67 5 Structure of the Simulation Process 69 5.1 Overall System Equations 69 C O N T E N T S XI 5.2 Program Structure 72 5.3 Time-Domain Simulation 73 5.3.1 Prediction of Voltages and State Variables 73 5.3.2 Forcing Functions 75 5.3.3 Discontinuity Checks 76 5.3.4 Construction of the Nodal Admittance Matrix 77 I I Component Modeling 8 0 6 T r a n s m i s s i o n L i n e 8 1 6.1 Derivation of the Model Equations 82 6.2 Lossless Line Model 84 6.2.1 E M T P Model for Lossless Lines 85 6.2.2 Interpolated Model for At > r 86 6.3 Lossy Line Model 88 7 S y n c h r o n o u s M a c h i n e 92 7.1 Model Equations with Real Signals 93 7.2 Transition to Complex Variables 96 7.2.1 Rotor Structure Quantities 98 7.2.2 Three-Phase Armature Quantities 99 7.3 Model Equations with Complex Signals 100 7.3.1 Flux Linkage Equations 101 7.3.2 Voltage Equations 103 7.3.3 Mechanical Equations 104 7.4 Time-Domain Simulation 104 7.5 Steady-State Initialization 105 8 S y n c h r o n o u s M a c h i n e P a r a m e t e r s 1 0 8 8.1 Frequency Response Tests 109 8.2 Transfer Functions I l l 8.2.1 The Q-Axis Admittance I l l 8.2.2 The D-Axis Admittance 112 8.2.3 Parameter Extraction 114 8.3 Practical Considerations 117 C O N T E N T S vii 8.3.1 Pre-Processing Filter 117 8.3.2 Adjustment of Lfkd 117 8.3.3 Relating Cd(s) and sG{s) 118 8.4 Results 120 III Simulation Examples 123 9 S imula t ion Case Studies 124 9.1 Transmission Line Energization 125 9.2 Subsynchronous Resonance Case 129 9.2.1 Steady-State Initialization 130 9.2.2 Electromagnetic Transients Simulation 134 9.2.3 Stability Analysis Simulation 140 9.3 Evaluation of Program Performance 143 10 Conclus ion 145 10.1 Recommendations for Future Research 145 10.1.1 Improved Simulation Speed 145 10.1.2 Component Modeling 146 10.1.3 Online Network Aggregation 147 10.1.4 System Studies 147 10.2 Summary 148 Lis t of Symbols 151 Bib l iography 155 Lis t of Tables 2.1 Comparison of different variable representations 25 3.1 Transfer admittance functions for an ideal integrator with the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algorithms 41 3.2 Truncation errors for the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algorithms 46 8.1 Comparison of identified electrical synchronous machine parameters 122 9.1 C I G R E line energization case data 126 9.2 IEEE subsynchronous resonance benchmark case data 130 9.3 Program running times for several test cases 143 viii List of Figures 1.1 Typical frequency ranges for electromagnetic and electromechanical transient phenomena 2 2.1 RL-circuit 18 2.2 Continuous amplitude spectrum of a bandpass signal 20 2.3 Discrete amplitude spectrum of a bandpass signal 21 2.4 Amplitude spectrum of in-phase and quadrature components 22 2.5 Peculiarities regarding complex signals 31 3.1 Numerical integration with the trapezoidal and relaxed trapezoidal algorithms 37 3.2 Numerical accuracy of the backward Euler (BE), trapezoidal (TR), and re-laxed trapezoidal (RTR) rules of integration 43 3.3 Numerical accuracy of the backward Euler (BE), trapezoidal (TR), and re-laxed convolution (RCV) rules of integration 44 3.4 Local integration errors for u> > UJ0 48 3.5 Local integration errors for to < LUQ 48 3.6 Global integration errors for u> > UJQ 50 3.7 Global integration errors for ui < UIQ 50 3.8 Stability region for the backward Euler integration rule (shaded area) 53 3.9 Stability region for the trapezoidal, relaxed trapezoidal and relaxed convolu-tion rules of integration (shaded area) 53 4.1 Membership functions for step size prediction 59 4.2 Monotonic chaining method to assess the prediction reliability. 62 4.3 Two membership functions for prediction reliability 64 4.4 Monotonic chaining method to assess the step change efficiency 65 4.5 Membership functions for efficiency evaluation 67 4.6 Denazification 68 5.1 Discretization of power system components 71 ix LIST OF FIGURES x 5.2 Flowchart: Program structure 73 5.3 Flowchart: Time step simulation loop 74 5.4 Open-circuit admittance matrix 78 6.1 Lossless line model using Bergeron's method for A t < r m ; „ 85 6.2 Lossless transmission line model for At > rmin 88 6.3 Schematic of a lossy transmission line model 89 7.1 Synchronous generator d and g-axis models 92 8.1 Admittance Aq(s) of the Lambton generator and its final approximation by a third-order polynomial 112 8.2 Admittance Aq(s) of the Lambton generator approximated by a first-order function of type Ti(s) 114 8.3 First residue of Aq(s) approximated by a first-order function of type ^(s). . 116 8.4 Parameters Lfkd and lfd as functions of frequency. The cross indicates the best fit, yielding Lfkd = -0.0282 p.u. and lfd = 0.831 p.u 118 8.5 Best concurrence of functions RIGHT (solid line) and LEFT (dashed line) indicates the optimum value of Ljkd 119 8.6 Measured operational d-axis inductance Cd(s) with its second-order analytical (solid line) and iterated approximations (dashed line) 120 8.7 Measured armature-to-field transfer function sG{s) with its second-order an-alytical (solid line) and iterated approximations (dashed line) 121 8.8 Measured operational (/-axis inductance Cq(s) with its third-order analytical (solid line) and iterated approximations (dashed line) 121 9.1 Network diagram for C I G R E line energization test case 125 9.2 Simulation results with a C I G R E line energization test case 127 9.3 Step size adjustment (0-800 ms) 128 9.4 Network diagram for IEEE subsynchronous resonance benchmark case. . . . 129 9.5 Steady-state mechanical torque 132 9.6 Steady-state rotor speed 133 9.7 Transients simulation (0-1 s): Mechanical torque 135 9.8 Transients simulation (0-1 s): Rotor speed 136 9.9 Transients simulation (0-10 s): Mechanical torque 137 9.10 Transients simulation (0-10 s): Rotor speed 138 9.11 Transients simulation: Automatic step size adjustment . 139 9.12 Stability simulation (0-1 s): Mechanical torque and rotor speed 141 LIST OF FIGURES xi 9.13 Stability simulation (0-10 s): Mechanical torque and rotor speed 142 Acknowledgement For their constant help and guidance throughout this work, I would like to express my deepest gratitude to my thesis supervisors, Dr. Hermann Dommel and Dr. Prabha Kundur, whom I respect and admire tremendously. Both their background and expertise in the areas of electromagnetic transients simulation and power system stability analysis not only made this work possible but also enriched my personal and professional life. Thank you both very much for your encouragement and support. Furthermore, I am very grateful for the help I received from Dr. Jose Marti, whose course and seminars provided valuable information on numerical integration and error anal-ysis, which have certainly influenced much of Chapter 3 of this work. Many thanks also to Dr. Takahide Niimura whose expertise in artificial intelligence and fuzzy logic helped me enormeously during the development of the automatic simulation step size control, in be-coming aquainted with the theoretical background necessary to understand and design fuzzy reasoning techniques. To Dr. A l i Moshref for his encouragement and technical advice to questions concern-ing black start and power system restoration. I further appreciate the cooperation and joint work with my dear colleagues and friends, Daniel Lindenmeyer, whose suggestions for the improvements of the intelligent decision making process had been gladly accepted; and Awad I. Ibrahim, whom I thank for his patience and perseverence in setting up many E M T P test cases needed to verify the linearly interpolated transmission line model for large simula-tion step sizes. At this point, I would also like to thank W. Scott Meyer, whose handwritten notes of 1973 served as an inspiration and provided valuable information for the development of this line model. I am much indebted to Guido Hempel who took on the intense and at times difficult task of setting up the subsynchronous resonance test case for simulations with the E M T P and Powerfactory simulation software. To Dr. E. Schmieg and his company for continu-ously keeping us up-to-date with the newest version of his software as well as for providing immediate advice in regards to technical and software related issues. The financial assistance of the Natural Science and Engineering Research Council of Canada, and of B.C. Hydro & Power Authority, through funding provided for the NSERC-B.C. Hydro Industrial Chair in Advanced Techniques for Electric Power Systems Analysis, Simulation and Control, are gratefully acknowledged. Finally, I would like to point out the friendly and cooperative atmosphere in the power group, which helped me in the accomplishment of this task. To all the people and friends who directly or indirectly contributed to this thesis project: Thank you so much; but fore-most I must thank my beautiful and magnificent wife Corina as well as her and my parents for their never-ending patience, constant encouragement and sensitive understanding. I am grateful for their love and support through the years. Vancouver, B.C. , Canada Sebastian Henschel February 9, 1999 xn Chapter 1 Introduction 1.1 General Background Major disturbances on electric power systems are a continuing source of concern to utilities. Programs for the computation of transients have therefore become an integral part for system design and control, in order to maximize the ability of a system to withstand the impact of a disturbance. Over time, these tools have become very efficient and reliable so that, today, many of them enable us to perform a wide range of power system analysis. Most commercial power system simulation software packages provide a number of different algorithms for specific tasks. In regard to time-domain simulation programs, the different tasks are characterized by their specific time frames, which are directly related to the frequency ranges that are encoun-tered with certain types of power system disturbances. These disturbances can be broadly classified as electromagnetic transients and electromechanical transients. Typical frequency ranges for various transient phenomena are shown in Fig. 1.1. 1.1.1 Electromagnetic Transients Electromagnetic transient phenomena are usually triggered by changes in the network con-figuration, which may be caused by the closing or opening action of circuit breakers or 1 1.1. General Background 2 Frequency f in Hz cu o> c <5 Cc o c CT fi CM <D O) c <5 cc o c CD Cr fi - 10 6 - 10 5 10 4 ) \- 10 3 10 2 I- 10 1 10° 10-1 I- IO"2 "IO"3 '— 10"4 SF 6 transients Wave propagation, lightning Switching overvoltages Transformer saturation Steady-state power flow Subsynchronous resonance Transient stability: machine ro-tor dynamics Interarea oscillations Mid-term and long-term stability: automatic generation control n c CD E o C CU s. Q. CD c o> re E o u 0) c u E o c 0) (0 o c (0 u 0) E o o 0) Figure 1.1: Typical frequency ranges for electromagnetic and electromechanical transient phenomena power electronic equipment, or by equipment failure or faults, such as a lightning stroke on a transmission line. Because these transients typically decay within a few cycles, their area of propagation is usually confined to "electric network equipment" such as transmission lines, cables, transformers and associated protective devices. Hence, for an accurate simulation of electromagnetic transients, these components must be modeled in great detail. On the other hand, due to the long time constants associated with the dynamics of "power plant equipment", such as generators and turbines, simplified models of such devices are often suf-ficient. In most cases, an ideal voltage source behind the transient or subtransient machine impedance is a reasonable representation of a power station for the time frame of typical electromagnetic transients studies. 1.1. General Background 3 The study of electromagnetic transient phenomena includes switching surges, transient recovery voltage, and ferroresonance. Electromagnetic transients programs ( E M T P 1 ) are common digital computer tools for the analysis of this class of transients. Typically, the simulation step size is of the order of tens of microseconds, although it may vary considerably depending on the type of electromagnetic phenomenon that is being studied. 1.1.2 Electromechanical Transients Electromechanical transients are slower transients that are caused by a mismatch between power production and consumption, and therefore involve the oscillation of machine rotors because of an unbalance between turbine and generator torques. The analysis of this class of transients is known as stability simulation. With respect to time-domain simulation, one can distinguish three different classes of stability analysis that have evolved in the past: 1. Transient stability. This type of stability analysis investigates the electrodynamic behaviour during the "first rotor swing", i.e. during the first few seconds after a disturbance. The disturbance is usually assumed to be a fault in the transmission system that may result in a temporary shut-down of equipment as a protective measure, thereby altering the flow of power which may impact the system stability and cause one or more generators to lose synchronism and fall out of step. 2. Long-term stability. This analysis is concerned with the long-term dynamic response of a system to severe system upsets. These may result in large excursion of voltage, fre-quency and power flows, which thereby invoke slow actions of controls and protections not modeled in conventional transient stability studies. The focus of this analysis is thus directed toward boiler dynamics of thermal units, penstock and conduit dynamics of hydro units, etc., assuming that inter-machine synchronizing power oscillations have subsided and uniform system-frequency has returned [55]. The time frame of typical studies of this type is in the range of several tens of minutes. 1 I n this text, this acronym is used as a general term for any digital program that can simulate electro-magnetic transients. Sometimes, we may write the E M T P in order to refer to one version of a simulation program that was first developed by H . W . Dommel in cooperation with Bonneville Power Administration, Portland, Oregon, U.S .A. , and which was originally named Electromagnetic Transients Program [6, 23]. 1.2. Review of Unified Time-Domain Simulation 4 3. Mid-term stability. Studies of this type are concerned with the transition of a system from a transient dynamic response to the long-term response. Hence, they focus on synchronizing power oscillations between machines, including the effects of some of the slower phenomena. Typical time frames for these studies lie in the minute range. An important assumption for all these types of analysis is that the exchange of energy between generators and other dynamic equipment takes place with the electric network re-maining at power frequency. It is immediately clear that, with this assumption, electromag-netic transients cannot be properly represented and must therefore be "filtered out" from the solution process. This is achieved using the quasi-stationary or quasi-steady state approach, where the network is modeled using the conventional phasor technique [79] while, contrary to this technique, the phasors are allowed to change in time, thus accounting for the dynamic response associated with the rotary and other mechanical equipment. This quasi-steady representation allows time-domain stability simulation programs to use relatively large step sizes in the range of 10 ms up to 100 s, depending on which type of phenomenon is studied. 1.2 Rev iew of Uni f ied T i m e - D o m a i n S imulat ion It is generally difficult to make a clear distinction between the above transient phenomena as they evolve with time and change their characteristic. For instance, what may have started out as an electromagnetic transient phenomenon in a transmission system may cause a nearby generator to fall out of step, creating a severe unbalance between power generation and consumption. This may cause a system to fall apart into several islands which, in turn, may impact associated boiler dynamics before normal operating conditions can be restored. It is evident that such a system would experience the entire spectrum of transients until it returns to steady state2. Having the common focus of investigating a system's power balance, stability simulation programs have emerged during the last couple of decades that aim at a synthesis of tran-2In this document, the term "steady state" signifies the balanced stationary operating condition of a power system at nominal frequency. Notice, however, that according to System Theory there are three additional forms of steady state, namely the periodic solution with superimposed harmonics, the quasi-periodic steady state, and chaos [70]. 1.2. Review of Unified Time-Domain Simulation 5 sient, mid-term and long-term stability analysis. In terms of modeling requirements, these simulation tools make use of conventional transient stability models as well as of an array of protection and control device models that are required for the mid-term and long-term stability simulation. In general, these programs provide modes for either fast or slow transient phenomena and switch between these modes using controlled step sizes, e.g. as described in [32]. In conjunction with an implicit integration method such as the trapezoidal rule, high-frequency components are filtered out during large integration steps, while the method is numerically stable and accurate for low-frequency modes. Simulation programs that gradually change their step size according to the system dynamics, e.g. in [31, 80], prefer Gear's integration method. Here, the step size is controlled by observing the truncation error of the integration polynomial. During high-frequency phenomena, the step size is reduced to a value below that corresponding to the largest eigenvalue. As high-frequency oscillations decay, the step size is gradually increased. Historically, electromagnetic transients programs have not been included in this unifying trend. Several reasons have contributed to this exclusion: 1. Study objective. Traditionally, studies of electromagnetic transients and of stability had different objectives. One was to investigate very fast transient phenomena that rapidly decay and usually have no impact on slower system dynamics. The other was to study the impact of a severe disturbance on the power balance of a system. 2. Simulation time. In the past, computers were relatively slow and the time spent for a simulation was therefore expensive. Even if it had been technically feasible to simulate a stability case study with E M T P , the simulation time requiring step sizes of approximately 50 LLS would have been far too long as compared to a simulation with, say, A t = 10 ms. In addition, the simulation time spans for a stability simulation are significantly longer than those for an electromagnetic transients simulation. 3. System size. As mentioned earlier, stability programs approximate the network side of an electric power system by employing phasor techniques. Assuming balanced system conditions, the network can be further reduced by using a single-phase representation. 1.2. Review of Unified Time-Domain Simulation 6 Traditionally, this reduction was necessary to accommodate large power system models on a digital computer. The system models for E M T P , on the other hand, could be kept comparatively small because electromagnetic transients decay rapidly and only have a local impact. The small size of the E M T P system model thus permits a detailed representation with all three phases. 4. Technical difficulties. Similar to unified stability programs, there are questions related to the representation of variables, integration methods and component modeling to fa-cilitate an efficient simulation of both electromagnetic and electromechanical transient phenomena. The fact that power systems have become increasingly sensitive to ill-conditioning and abnormalties, which have already resulted in several system failures and power outages in recent years, has stimulated the desire to improve E M T P regarding the capability to perform transient stability analysis. This was addressed in [12, 13] with respect to an evaluation of available analysis tools for power system restoration. It strongly emphasizes the notion for also considering electromagnetic transients in the context of a unified time-domain simulation program. Probably the earliest contribution in this direction-was made in 1979 by B. Kulicke who developed the simulation program N E T O M A C in cooperation with Siemens A G [50, 51]. Only recently was this idea adopted in some other power system simulation packages, e.g. in [22]. Essentially, the software provides a possibility for the program user to specify a time when the simulation should switch from an electromagnetic transients solution algorithm to a stability analysis. This option alone is a simplistic solution to a much more complicated problem, because the technical difficulties that were mentioned above are not resolved by this approach. This becomes clearly visible whenever the program switches from the EMTP-mode to the stability mode: Resulting plots and graphs show noticeable "jumps" where these changes occur [54]. In addition to being misleading in the final print-out, these artificial "jumps" may also become dangerous to the solution process itself, e.g. a protective relay model may interpret it as a sudden overvoltage and trip, thereby altering the system configuration incorrectly for 1.2. Review of Unified Time-Domain Simulation 7 the following simulation. The underlying causes for these "jumps" are the differences in the solution algorithms for E M T P and stability programs (cf. quasi-steady approximation). These fundamental differences, and assumptions made for both types of transients, have so far precluded a combined analysis of both electromagnetic and electromechanical transient phenomena. A possible way to avoid the "jumps" in the solution was presented by J. V . Sarlashkar et al. [76]. Whereas this work maintains the quasi-steady approximation for the simulation of slow transients, the transition between a transient-stability analysis and an E M T P simulation is carried out using signal modulation and demodulation techniques. The decision for this transition is based on an inspection of the spectral content of the signals. Other approaches combine the areas of electromagnetic and electromechanical simulations by representing the full network and machine differential equations without making the quasi-steady approximation. The difficulty then is to find an efficient solution algorithm that is accurate and numerically stable for a wide frequency range. Attempts to overcome this difficulty were proposed in [39], and in the recent work of V . Venkatasubramanian [84], and of U. Linnert [58, 59]. Common to all these approaches is the use of a complex variable representation. Since voltages, currents, etc. are real signals, the task of finding a meaningful imaginary part is not trivial. Linnert's method uses the space-vector transformation [21] to obtain a complex signal, plus a real-valued zero sequence signal, from three-phase quantities. The problem with this approach is that it requires balanced networks in order to decouple and isolate the z-component for the simulation [48]. Also, the fault modeling in this transformed rep-resentation is rather complicated as it depends on certain configurations of inductances and capacitances. A more general concept was presented by Venkatasubramanian who demonstrated an example of a three-phase power system model consisting of ideal voltage sources and linear RLC-elements. He suggested the use of ejult as a source function instead of the usual real-valued sinusoidal functions. However, as a future obstacle, he envisioned the application of this method to elements with a nonlinear characteristic. 1.3. Thesis Motivation and Objective 8 A l l the authors could not resolve the difficulties with respect to modeling: One of the major differences between electromagnetic transients simulation and the analysis of stability lies in the representation of the transmission network. While RLC-equivalents may suffice in stability studies, higher frequency components in electromagnetic transients simulation may be not adequately modeled by linear RLC-equivalent circuits. The research work described in this thesis has focussed on the development of a unified solution method for the simulation of both types of transients. Several aspects in regard to modeling, and difficulties associated with this method, have been investigated and are presented here. 1.3 Thesis M o t i v a t i o n and Objective With the availability of powerful computer hardware as well as improved analytical tech-niques, a distinction between fast and slow dynamics becomes less significant, and already many of today's time-domain simulation programs offer a unified approach for assessing the transient, mid-term and long-term stability of a system. In this light, it seems unjusti-fied that the trend to a unified simulation program stops short of including electromagnetic transients in such a solution algorithm. Considering recent trends towards economic competition and deregulation of the power industry worldwide, power systems have become increasingly vulnerable. Major power sys-tem blackouts in recent years and their subsequent restoration have been characterised by a multitude of different phenomena and abnormal conditions. While power system planning, design and operation all depend almost exclusively on digital computer simulation, it may be beneficial to revise some of the fundamental assumptions that simulation programs are based upon and adjust them to this shift of paradigm faced by the power industry. The need to expand electromagnetic transients programs to facilitate transient stability considerations was pointed out by Electricite de France and the Hydro-Electric Commission of Tasmania in a survey conducted by C I G R E in 1992 [12]. The intention of this questionnaire was to collect information about problems to be taken into account in restoration and black start analysis, and to gain insight into the present situation of analytical tools for this type 1.3. Thesis Motivation and Objective 9 of analysis. During system restoration, the systems are initially very weak so that electromagnetic phenomena, which may be considered insignificant during normal system operation, may have a substancial impact on the system, and can easily provoke a response of the slow dynamics as well. Because of high frequency phenomena, the quasi-steady approximation of the transmission system in conventional stability programs may not be realistic. Cases of subsynchronous resonance [47] have clearly shown that, it is not always possible to make this quasi-steady approximation. Although these cases can be computed by elec-tromagnetic transients programs, the aforementioned limitations indicate that they do not present an efficient tool for a substantial analysis of this kind. The demand for a tool to efficiently combine the fields of electromagnetic and electrome-chanical transients simulation into one program was the incentive for the research work described here. Hence, the objective of this thesis work is to develop a concept for the combined simulation of electromagnetic and electromechanical transient phenomena. This concept must not only resolve difficulties associated with the numerical integration but also consider modeling aspects in order to find a consistent method to reflect both types of transient phenomena throughout the simulated power system. At first, existing difficulties for a combined simulation must be carefully evaluated before a suitable variable representation can be found. This is described in Chapter 2. The follow-ing chapter is concerned with the numerical integration of this new variable representation. A new integration method is introduced, and its numerical errors and stability are analyzed. Chapter 4 presents an efficient step size control that employs fuzzy logic and artificial intel-ligence decision making processes to adjust the time step according to the dynamics in the system. How all the different parts are being interfaced to make up a simuation algorithm is the focus of Chapter 5. Chapter 6 presents a new transmission line model that allows a smooth transition from electromagnetic transients phenomena to very slow dynamics and vice versa. Synchronous machine models are the subject of Chapters 7 and 8, where the former chapter discusses the modeling aspects while the latter introduces a new method to obtain the machine parameters 1.3. T h e s i s M o t i v a t i o n a n d O b j e c t i v e 10 f r o m s t a n d s t i l l f r e q u e n c y r e s p o n s e m e a s u r e m e n t s . I n C h a p t e r 9 t h e c o n c e p t t h a t h a s b e e n d e v e l o p e d i n t h i s r e s e a r c h w o r k i s t r i e d o u t o n p r a c t i c a l c a s e s t u d i e s . T h e r e s u l t s a r e e v a l u a t e d a n d c o m p a r e d t o a r e f e r e n c e s o l u t i o n o b t a i n e d w i t h t h e E M T P . T h e l a s t c h a p t e r c o n t a i n s c o n c l u s i o n s d r a w n f r o m t h i s w o r k , as w e l l as r e c o m m e n d a t i o n s f o r f u t u r e r e s e a r c h w o r k . Part I Numerical Considerations 11 Chapter 2 Complex Variable Representation 2.1 Evaluat ion of E x i s t i n g Variable Representations The variable representation commonly used in electromagnetic transients programs is by means of real signals (instantaneous values is used as a synonym). The advantage of this representation is that real voltages, currents, fluxes, etc. are measurable quantities and can be easily compared with field test measurements. Since the time steps used in an electro-magnetic transients simulation are very small, i.e. in the order of 50 fxs, instantaneous values can be integrated with standard integration methods, for example with the Adams-Moulton, Runge-Kutta or Gear family of algorithms. Numerical integration, however, is bound by maximum step sizes. While the maximum possible step size for bandlimited signals is dictated by the Nyquist criterion or sampling theorem for the highest appearing frequency, 1 Atmax — ——- , (2-1) " I max a direct application of this theorem to power system analysis is not straightforward, as the highest appearing frequency is not always known prior to the simulation, nor can we generally assume the presence of bandlimited signals. Therefore, a safety margin is introduced by the following recommendation: 1 L^*max — cj , n > ^ & Jmax where k = 5 [24]. The maximum frequency depends, of course, on the phenomenon that is being studied. If there is no such phenomenon present during a simulation, or the solution 12 2.1. Evaluation of Existing Variable Representations 13 has returned to the steady state, the maximum possible step size is determined by the frequency of the forcing functions, i.e. the sources or generators. In a 60 Hz A C system the step size can therefore not exceed A t m a x = 1.67 ms in conjunction with the real variable representation. It must be emphasized that variable step size integration rules, such as the Runge-Kutta or Gear family of algorithms, are also bound by this maximum step size when they are applied to the real variable representation, although this limit may be marginally smaller or larger for the specific rule in use. Due to this limitation, it is clear that this representation is not well suited for stability simulations. We also see that the use of large step sizes not only depends on the integration method used but also on the variable representation. Since conventional stability programs use time steps ten times as large, let us now investigate how this is achieved: Stability programs use a complex variable representation. The power frequency, i.e. the 60 Hz component, is filtered out and only the frequency deviations from this frequency are observed. The filtering process is performed differently for the network and electric machines, and is described below. 2.1.1 Park's Transformation Synchronous machines are usually best described with dqO coordinates. To obtain these coordinates, a transformation is used that maps three-phase abc quantities, such as the instantaneous values of voltages or currents, onto a reference frame that rotates with the electrical rotor speed tvr(t) = -j- com{t) (2.3) where np is the number of machine poles, and ujm is the mechanical rotor speed. The angle of this rotating reference frame (d-axis) with respect to phase a of the stationary abc coordinate system is given by The transformation from the stationary armature three-phase abc reference system to the rotating dqO-system was first proposed by Blondel, and further developed by Doherty, Nickle 2.1. Evaluation of Existing Variable Representations 14 and Park [24]. It is well known and can be found in almost every textbook on electric machines. Its main purpose is to achieve constant inductance values in the machine equations in the dqO-system, which used to be a crucial requirement for the use of transient network analyzers (TNA). In computer programs this transformation is no longer required, but there are still advantages in using the dqO coordinates [55]: • They are easy to be physically interpreted. • For balanced conditions, zero sequence quantities disappear. • The parameters associated with the direct axis and quadrature axis can be directly measured from terminal tests. This last point is illustrated in more detail in Chapter 8. In this thesis, we will utilize the power-variant transformation which is widely used by the electrical utility industry and by generator manufacturers for reasons explained in [38]. The transformation from abc to dqO coordinates can be written as sd{t) sq{t) | = s0(t) COS0 cos ( 0 - ^ ) cos ( f l+^ ) s i n ( 0 + f ) - s in# - s i n ( 0 - f ) (2.5) Here, sa(t), Sb(t), sc(t) can be taken as the abc components of any quantity, e.g. voltage, current, flux, etc. The inverse transformation is given by - s i n 0 1 ~ cos 9 cos (6 — cos {9+\) • s i n ( 0 - f ) 1 3 , s i n ( 0 + f ) 1 2TT> 3 / (2-6) The angle 9 is also called the transformation angle. Whereas, in electromagnetic transients studies, this angle is taken as 9 — j3(t) [24], in stability simulations it is mandatory to introduce 9 = (3(t) — to0t in order to filter out the power frequency component U)Q. This implies that the quantities sa(t), Sb(t), and sc(t) take on constant values at steady state, i.e. when During rotor oscillations, of course, these values vary according to the frequency deviation Au; 2.1. Evaluation of Existing Variable Representations 15 2.1.2 Clarke's Transformation In order to interface the frequency-filtered variables with the network equations, one more transformation is applied in stability programs: The a(30 transformation, due to Edith Clarke [14]. Again, we use the power-variant form of this transformation, which is defined as sp(t) I = s0(t) 1 - \ 0 f 1 l 2 2 _ 1 2 v/3 2 1 2 (2.7) Applying this transformation to the abc components obtained with (2.6), yields the following direct transformation from dqO quantities to af30 components: cos 9 — sin 9 0 sin 6 cos 9 0 0 0 1 (2-Assuming balanced system conditions, the zero sequence components are always zero and can therefore be discarded from the above equation. Since the a and (3 components as well as the d and q components are orthogonal, they can be written in complex form so that the above transformation simplifies to: sa(t) +jsp{t) [sd(t)+jsq(t)}. (2-9) It can be verified that, with this variable representation, the above complex quantities become identical to phasors when the simulated system is at steady state. Yet, during electrome-chanical oscillations, these "phasors" deviate from their steady-state value. As we will see later in this chapter, these time-varying "phasors" in stability programs are the single phase equivalent of dynamic phasors'. Because it was traditionally assumed that electromagnetic transients have decayed before the response of the slower machine dynamics takes place, stability programs model the transmission network with its steady-state, single-phase equivalent in the frequency domain. This is identical to the ordinary phasor representation used in power flow and other steady-state computations. However, due to the fact that the "phasors" are time-varying in stability simulations, the network representation is often called a quasi-steady equivalent. In any case, 2.1. Evaluation of Existing Variable Representations 16 it must be observed that the network model is an approximation which is only correct if the time deviations of these dynamic phasors are small. Since the dynamic phasors vary according to the frequency deviation ALO = u>r—uj0, which only depends on the slowly changing rotor speed, the simulation time steps for a stability simulation can be very large. However, the magnitude of dynamic phasors may also change during the transient period which permits practical step sizes of approximately 5-10 ms for a transient stability simulation. 2.1.3 Summary Let us summarize the observations made from the two presented variable representations and draw some conclusions for the work to follow. Aside from the approximations made regarding the transmission network, it appears that the variable representation made in stability simulation is superior to the use of instantaneous values, when large time steps are desired for the solution process. For a combined simulation of both electromagnetic and electromechanical transients, the use of complex variables seems more efficient, provided that a more realistic solution for the quasi-steady state network approximation can be found. The reason for this efficiency can be characterized by two steps: First, a complex number (with the synonyms vector or phasor, depending on the method) is constructed from two values which are derived through similarity transformations from given physical three-phase quantities (e.g. currents, voltages, fluxes, etc.). Second, for an efficient simulation of quasi-steady states, the oscillating motion of this vector in the complex plane is slowed down by system frequency through multiplication by a factor e~julot. The process of slowing the oscillations is also referred to as relaxation in this thesis. Interestingly, the method by U . Linnert to combine electromagnetic and electromechani-cal transients simulation [58, 59] is based on the same approach, but utilizes the space-vector transformation [21] instead of Park's transformation to realize the filtering of the system fre-quency component. 2 . 2 . Complex Signals and Dynamic Phasors 17 Based on the above observations, research on a simulation w i th increased step sizes falls into the following two categories: • Representation of system variables: Variables such as voltages and currents are measurable quantities and therefore real. In order to slow down the oscillatory motion of such variables, one must be able to, at least approximately, extract information about the frequency or spectral content f rom this variable. Because of this, the variable must be translated into an adequate representation that also allows the subtraction of a part of its spectral content. The latter could be easily achieved by a complex representation through mult ipl icat ion by e~^Uot, yet the construction of a complex variable f rom a real voltage or a real current is rather challenging. • Integration of differential equations: The question is whether the variables in their new representation can st i l l be integrated. Unfortunately, this is not always the case. Secondly, i t is of great interest to find a numerical method that is most suited for integration because the gain in simulation step size achieved through the new variable representation may be lost i f the integration rule cannot accommodate large step sizes and becomes unstable. Both areas are closely related. In the following section, the author introduces a feasible and flexible representation for system describing variables. An analysis concerning their numerical integration is presented in Chapter 3. 2.2 C o m p l e x Signals and D y n a m i c Phasors This section expands on the idea of using complex numbers for system describing variables. Whi le such complex representations have already been used in stabi l i ty simulations to simu-late AC power systems w i th large step sizes (cf. previous sections), the following derivation of complex signals and dynamic phasors differs f rom these representations in one important aspect: Instead of applying three-phase similari ty transformations to obtain one complex number and a real zero sequence value (which is usually neglected in stabi l i ty simulations), the dynamic phasor approach expresses each phase quanti ty as a complex signal. This way, 2.2. Complex Signals and Dynamic Phasors 18 the method is not restricted to three-phase systems but can be applied to an arbitrary number of phases. This greater flexibility may be significant in cases where systems are unbalanced and/or electromagnetic coupling involves more than three phases, for instance, in two parallel three-phase transmission lines. The idea of using complex signals to describe power system dynamics is based upon techniques for simultaneous time-frequency analysis in signal processing [15, 20, 57]. In this section, we show how this technique can be applied to power system simulation and explain how complex signals are constructed from real quantities. Dynamic phasors are obtained by multiplying the complex signals by e--7"0*, thus relaxing their oscillating motion. Before beginning to define complex signals and dynamic phasors, let us first look closely at the more familiar relationship between real signals, i.e. instantaneous values, and phasors. It will then be easy to establish a similar relationship between complex signals and dynamic phasors. 2.2.1 Real Signals and Phasors i(t) R e(t)W Figure 2.1: RL-circuit Let us consider the steady-state solution of the RL-circuit shown in Fig. 2.1, as often needed for the ini-tial conditions of a transient simulation. For the pur-pose of this example, a single phase representation is assumed, albeit a notation in matrix form easily ex-pands the derivations to the poly-phase case. The transient behavior of the circuit can be described by the equation v(t) = Ri(t) + L^i(t), (2.10) where v(t) is equal to the user-specified A C source function e(t), which could be defined, say as e(t) = va cos [tJot\. The system is thus described by v(t), i(t), and e(t) which are all real and measurable quantities. Generically, these quantities are referred to as real signals or instantaneous values and denoted with the symbol s(t). At the beginning of most transient studies, initial voltage and current quantities are 2.2. Complex Signals and Dynamic Phasors 19 obtained from either a power flow [3, 36] or some other type of steady-state solution [24]. In either case, the solution is presented in phasor form [79] with a real and imaginary part for each quantity. Let us briefly review the relationship between real signals and their associated phasors. Phasors are the frequency representation of real signals under the assumption that the latter oscillate in sinusoidal motion of constant frequency U>Q and are of constant amplitude. In the frequency domain, (2.10) becomes where V ( t J o ) and I(u0) are the phasors. The equation establishes the position of the phasors relative to each other but it does not give information about the relationship between phasors and their corresponding real signals. This relationship is established by the projection of the phasors onto the real axis1. In general, denoting phasors as S(cv0), we have where t = t0 is the time of the beginning of the transient simulation. Conventionally, this time is set to to = 0 in most power system simulation programs. Different definitions for phasors have been reported in the past: Whereas some describe them as constant quantities with respect to time [37], others attribute a uniform rotation of synchronous speed u>0 to them [56, 78, 79]. In this monograph, we assume the former definition and characterize phasors as the time-invariant frequency transformation of a si-nusoidally alternating steady-state signal. Independent of which definition is preferred, it is important to realize that phasors are a representation of steady-state signals. Next, we are going to introduce two different ways to describe transient signals, namely as complex signals or dynamic phasors. 2.2.2 I n t r o d u c t i o n t o C o m p l e x S i g n a l s a n d D y n a m i c P h a s o r s The concept of representing signals in complex form had been developed early this century by various researchers. Since that time, the theory has been much refined and led to several 1 Often the projection is accompanied by a scaling factor of \J2 to account for the difference between the amplitude's peak value and its effective or R M S value, or \J2/Z to make phasors power-invariant (in three-phase systems). Such factors are omitted in this analysis, thus letting phasors correspond to the sinusoidal peak value of their associated real signal of one single phase. V(LV0) = [ R + j cv0 L] I(u0) (2.H) s{t0) = M[S(LU0)}, (2.12) 2.2. Complex Signals and Dynamic Phasors 20 discoveries in engineering fields such as signal processing and radio transmission. However, with the exception of steady-state phasors, this theory had long been ignored in power engineering. The advantage and usefulness of complex signals and dynamic phasors for transient power system analysis was first studied and pointed out in 1994 in a publication by Venkatasubramanian [84]. In that paper, time-varying phasors (the term used there for complex signals) were used instead of steady-state phasors to decompose three-phase quan-tities into symmetrical components. With this representation he could achieve an efficient transient simulation of unbalanced faults and/or network conditions. The work presented here is not in any way tied to three-phase transformations. It rather attempts to generalize existing concepts in order to overcome the three-phase limitation. This is achieved by combining the complex variable representation with a powerful integration technique in the phase domain. Let us now introduce the mathematical background for this complex representation. Power system transient phenomena are typically characterized by spectra containing very high frequency components. Yet, as the system returns to a steady state, the spectrum of the signals becomes increasingly centered around system frequency co0, with its half-bandwidths being generally less than ui0, as shown in Fig. 2.2. Such signals are called bandpass signals [20]. Using Fourier analysis, a bandpass signal can be visualized as the sum of an infinitely large number of very closely spaced sine waves, as demonstrated in Fig. 2.3. Thus t |s(co)| - c o n 0 C0n CO Figure 2.2: Continuous amplitude spectrum of a bandpass signal. 2.2. Complex Signals and Dynamic Phasors 21 Aco |s(co)| -cof C O Figure 2.3: Discrete amplitude spectrum of a bandpass signal. s(t) = lim a,i cos \(co0 + iAco) t + Au ; ->0 — ' i=—oo Applying the trigonometric identity cos A + B = cos A cos B — sin A sin B the above signal can be rewritten as s(t) lim £ cii cos [iAco t + (f)i\ COSCO0t Thus, we have s(t) — S[(t) cosLu0t — SQ(t) sinw0^, where sj(t) = lim > di cos [iAco t + 4>i A w - > 0 ^ — ' i=—oo oo S<?C0 = ^ m a, sin [z Au; £ + </>; Au—»0 ' (2.13) (2.14) (2.15) (2.16) (2.17) (2.18) Equation (2.16) represents the bandpass signal in terms of the two lowpass signals s/(i) and SQ(£) and two sinusoidal carrier signals coswot and s'mco0t. Both carrier signals are of the same frequency but are 90° out of phase, i.e. in phase quadrature. Hence, the two lowpass signals si(t) and SQ(£) are called the in-phase and quadrature components of the bandpass 2.2. Complex Signals and Dynamic Phasors 22 signal. Their two-sided amplitude spectra are identical, and they are also the same as the positive amplitude spectrum of s(t) but shifted down in frequency LOQ. This is depicted in Fig. 2.4. s ( c o ) | I - c o r CO Figure 2.4: Amplitude spectrum of in-phase and quadrature components Since si(t) and SQ(£) contain all the information about the bandpass signal s(t), one possible way of representing the real signal s(t) is by the means of its dynamic phasor: @[S(t)] = Sl(t)+jsQ(t). (2.19) In signal processing, @[S(t) ] is also called the complex envelope of the signal s(t). The key advantage of introducing dynamic phasors as a lowpass representation of a bandpass signal lies in the requirement of far less samples than might be expected on the basis of applying the Nyquist criterion to the highest signal frequency of s(t). Another representation of the real signal s(t) is the analytic signal, which is also termed complex signal in the context of this thesis. It is defined as S(t) = @[S{t)]ejulot. Substituting (2.19) into this equation yields = [si(t) + j sQ(t)) [cos uj0t + j sin co0t] = si(t) cosu>ot — SQ(t) sinu0t + j [si(t) sin u)0t + s<g(£) cosco0t] = s(t) + j H[s(t)] where (2.20) H[s(t)] = S(T) TX .1 t — T oo dr (2.21) (2.22) 2.2. Complex Signals and Dynamic Phasors 23 denotes the Hilbert transformation. A formal proof for (2.21) is presented by L. Cohen in [15]. The important implication of this equation is that the real signal s(t) is retained in the real part of the complex signal S(t), while the transition from a dynamic phasor to the complex signal can be easily obtained by phase modulation with the carrier frequency UJQ. Let us next answer how these concepts can be applied to power system analysis. While most of the theory on complex signals has been developed in the area of signal processing, the techniques for their application to power system simulation differ widely in the following aspect: In signal processing, one is primarily concerned with finding the imaginary counter-part from the given real signal in order to construct a complex signal, e.g. by evaluating the Hilbert transform. For power system simulation, the signal—real or complex—is initially unknown and must be generated in the course of the simulation. This can be done with prior knowledge of the system differential equations, ofthe forcing functions and of the initial conditions, i.e. steady-state conditions. Let us begin in the reverse order: For a time-domain simulation to commence we seek the steady-state initial conditions using complex signals and dynamic phasors. At steady state, all signals are sinusoidal and have only one frequency component, namely system frequency cu0. In this case, equations (2.17) and (2.18) for the in-phase and quadrature components reduce to Sj(t = t 0) = a0 cos</>o (2.23) sQ(t = t0) = a0sin</>0 (2.24) Consequently, the corresponding dynamic phasor is ®[S(t = to)] s/(t0) +jsQ(t0) = a0e- (2.25) which can be readily determined from the steady-state phasor: a0 = | S(u)0) (po = arctan S[S{u0)] &[SM] (2.26) The initial complex signal is then obtained with S(t = t0) @[S(t0)]e J ^ 0 to (2.27) 2.2. Complex Signals and Dynamic Phasors 24 It should be noted that usually t 0 = 0, in which case the initialization of dynamic phasors and complex signals is identical. Previously, it was assumed that the forcing functions, such as the source voltage of the example on page 18, were user-defined, real functions. However, with the introduction of complex signals, these functions now also need to be complex. The imaginary part can thus be constructed using the Hilbert transformation. For the sinusoidally oscillating source func-tion of the aforementioned example, we have e(t) — va cos [cu0t]. The Hilbert transformation gives: Similarly one could proceed for any arbitrary forcing function. It should be mentioned, at this point, that the evaluation of the Hilbert integral is not always easily achieved. Yet, for the simulation of A C systems, the use of sinusoidal forcing functions should be sufficient in most cases. In later chapters of this thesis, we also introduce DC source functions, e.g. to represent a constant generator field voltage or a constant turbine torque at the machine shaft. However, since these already are D C lowpass signals with their spectra centered at OJ = 0, they cannot be demodulated to a lower carrier frequency. Hence, in these cases, it is sufficient to only specify the real, instantaneous forcing function. Finally the system differential equations have to be modified to accommodate complex signals and dynamic phasors instead of real signals. Since the real and imaginary components of complex signals are orthogonal, they can be solved independently and consecutively with (2.10). Thus, again refering to the RL-circuit in Fig. 2.1, we have U[e(tj\ = va sin [oo0t] (2.28) so that the new, complex forcing function is given by E{t) = e(t)+jn[e(t)]=vae jul0 t (2.29) (2.30) (2.31) V(t) = RI(t) + L-I(t) (2.32) 2.2. Complex Signals and Dynamic Phasors 25 The orthogonality also allows the same mathematical operations, e.g. summation, multipli-cation, integration, etc. to be performed on complex signals as well as instantaneous values. Several methods for numerical integration of the above differential equation are presented in Section 3.1. Applying the phase demodulation (2.20) for dynamic phasors to the above equation, we obtain 9\ V(t)] = + j OJ0L^J ®[ lit) } + Ljf$)[ I(t)] (2.33) It can be seen that the resulting equation resembles very much the phasor equation (2.11) in the frequency domain but retains the time-derivative from (2.32). A comparison and summary of the variable representations discussed so far can be found in Table 2.1. Table 2.1: Comparison of different variable representations I n s t a n t a n e o u s V a l u e s v(t) = Ri(t) + Ljti(t) C o m p l e x P h a s o r s V(cu0) = (R + jco0L) I(co0) C o m p l e x S i g n a l s V(t) = RI(t) + LJtI(t) D y n a m i c P h a s o r s &[V{t)] = (R + JLu0L)®{I(t)}+Ljt@[I(t)} steady-state part transient part Let us now conclude this section with a brief definition of complex signals and dynamic phasors from the point of view of power system analysis: D e f i n i t i o n 2.1 A complex signal is a means of representing system describing quantities such as voltage, current, flux, etc. The real part of the signal is identical to the real signal (or instantaneous value), which generally corresponds to a measurable quantity. There does not exist such comprehensive definition for the imaginary part. The only requirement is that, at steady state, it must complement the real part in such a fashion as to construct a 2.2. Complex Signals and Dynamic Phasors 26 synchonously rotating complex signal of constant magnitude. This initial condition together with the system equations and forcing functions fully determine the time-function of the imaginary part during a transient simulation. Definition 2.2 A dynamic phasor is the frequency modulation of a complex signal onto the negative system carrier frequency —UQ. Its spectrum is therefore the spectrum of the complex signal shifted by system speed —OJQ to suppress parts of its oscillatory dynamics. Consequently, during steady-state initialization, dynamic phasors are identical with their corresponding phasor. With the above concepts and definitions, we have introduced useful tools for the simu-lation of steady states and quasi-steady states of A C power systems. This shall be briefly demonstrated with the following example: Example 2.1 After loss of generation or after load pick-up during a system restoration process, the system frequency may experience deviations which are typically within 2% (58.8-61.2 Hz in 60 Hz systems). During severe outages sustained variations of as much as 5-10% were reported [65]. Since most integration techniques for power system simulation assume a straight line over one integration step, the simulation accuracy of such a system depends on how well the present frequencies are linearized or, in other words, on how small the integration step must be to approximate all the dynamics during this step with a straight line. Assuming that 20 steps per cycle give a good resolution, the step size for this case is restricted to about 0.8 ms. Dynamic phasors suppress the 60 Hz component. In order to accurately simulate the dynamics of the above case, the allowable simulation step size would be approximately 8.0 ms using the integration method for dynamic phasors (see Chapter 3). Use of this larger step size thus increases the speed ofthe simulation by a factor of 10. Notice that transient stability programs use about the same integration step size, but do not account for frequency deviations in the network. This example shows that dynamic phasors can be a powerful tool compared to the tech-niques currently implemented in power system simulation programs, particularly if the sim-2.3. Spectral Properties of Complex Signals 27 ulation times involved are long. For short simulation times, however, there will not be a significant speed advantage by using dynamic phasors. In fact, due to complex arithmetic, the solution is expected to be slightly slower than conventional E M T P programs with real arithmetic. This is especially true since for most simulations there is a fault at the beginning that causes high frequency oscillations for which small time steps are required, no matter whether dynamic phasors are used or not. Yet, once these oscillations begin to decay, the techniques described here become beneficial. 2.3 Spectral Propert ies of C o m p l e x Signals Knowledge of the transient state of a system is important if variable step sizes are used to simulate the network. The step size must then be adjusted according to the transients of the system to guarantee accuracy of the solution. Therefore, knowledge of the frequencies that partake in the simulated phenomenon are of great interest. Let us, therefore, investigate some of the spectral properties of complex signals and evaluate their usefulness for the simulation process. 2.3.1 Spectral Averages Let us first consider real signals. The spectrum of a real signal is given by its Fourier transform: oo s{ui) = — L f s(t) e~ju}tdt (2.34) V27T J —oo The average frequency (u) of a signal s(t) over the entire duration of the phenomenon it describes is the sum over all frequencies multiplied by their particular spectral density |s(w)|2. It can be expressed either in the frequency domain or in the time domain. The following two relationships define average frequency [15]: oo oo (u)= J u\s(u)\2duj = -j J S*(t)^S(r) dt. (2.35) T=t Notice that the above equations are true for real signals s(t) as well as for complex signals S(t). The spectrum of a real signal satisfies s(—u>) — s*(u>) and therefore the energy 2.3. Spectral Properties of Complex Signals 28 density spectrum |s(o>)|2 is always symmetric about the origin. Because of this symmetry, the average frequency will always come out to be zero. Albeit correct in a purely mathematical sense, such a result is physically not meaningful as it does not indicate which frequencies are involved in the phenomenon. In general, there are two approaches to remedy this dilemma: 1. To calculate the average frequency, the negative part of the spectrum can be ignored so that the average frequency indicates, roughly speaking, the maximum spectral density of positive frequencies alone. Thus, oo (co) = J co\s(uo)\2dco (2.36) 0 The integral from zero to infinity takes the place of the one from minus to plus infinity. The method of setting the negative spectrum to zero was pursued by Laplace, Carson, and Heaviside. 2. A new signal S(t) can be defined whose spectrum is only positive. If this signal is the complex signal, we can decompose it into its magnitude and phase, and write S(t) = \S(t)\eJ*W. Calculating the average frequency of the complex signal by integrating from — oo to oo, we arrive at (co) oo / { m- jlw§\^ im2dt (2' 37) Since the average frequency is always real, the imaginary second term must be zero. It can be seen that this term is a perfect differential that, indeed, integrates to zero. Hence, 00 (co) = J <p(t) \S(t)\2dt (2.38) — 00 This integral represents the sum over all quantities ip(t) multiplied by the density IS i^)!2 of the signal, i.e. the square of its magnitude. The quantity ip(t) must thus be the instantaneous value of what we are building the average of, hence it can be rightly called an "instantaneous frequency" or the frequency at each time [15]. 2.3. Spectral Properties of Complex Signals 29 The average frequency of the phenomenon described by the real signal s(t) is contained in the phase modulation of the complex signal S(i). The concept of instantaneous frequencies is explained in more detail in the next paragraph. 2.3.2 Instantaneous Frequency The spectrum of a signal can be obtained by utilizing Fourier's transformation. Numerically, it has become relatively fast and reliable [73]. Yet, Fourier transform pairs, such as a signal being a function of time and its spectrum being a function of frequency, are bound to Heisen-berg's time-bandwidth product theorem, or uncertainty principle. Briefly summarized it says that the bandwidths of both frequency and time cannot be made infinitely small at the same time, that is, if one bandwidth is small the other one is big. Therefore, this transformation cannot provide the information of which frequencies are present at a particular instant. However, "frequencies are related to time and they change with time, as can be noticed every day in terms of e.g. light of changing color or sounds of varying pitch" [15]. For the purpose of knowing the current transient state of a system or network, we are more interested in a close relationship between time and dominant frequencies rather than knowing all frequencies that appear in the course of the simulation without being able to pinpoint exactly when they were present. Previously it was shown that the average over all instantaneous frequencies ip(t) is iden-tical to the average frequency of the positive spectrum. Yet, this does not automatically imply that the instantaneous frequency is the average of the positive spectrum at a par-ticular time! It does demonstrate though that, at a particular time, there exists only one instantaneous frequency per signal while there can be a multitude of frequencies present in the spectrum. Hence, the instantaneous frequency is not identical with a signal's spectrum, although both may coincide. While Heisenberg's principle applies to Fourier transform pairs such as the signal and its spectrum, it does not in any way restrict the concept of instanta-neous frequency. A detailed analysis of properties of time, frequency and their interactions, as well as a reference to other research work done in this area (signal processing) is given by Cohen [15]. 2.3. Spectral Properties of Complex Signals 30 Definition 2.3 Instantaneous frequency is the time-derivative of the phase angle of a com-plex signal. Hence, it can be understood as the angular speed of the signal. In this document, instantaneous frequency is denoted as u)i(t) = (p(t). Notice that the above is a rather mathematical definition, based on whatever method is used to create a complex signal [33, 69, 83]. Since signals are real by nature, the choice of a corresponding imaginary part is not obvious. The concept of instantaneous frequency thus requires careful interpretation and an intuitive insight into the physical behavior of the studied phenomena. The work presented in this thesis utilizes instantaneous frequency as an indicator as to whether the transient solution has arrived at steady state or is shortly before arriving there, the so-called quasi-steady state. The steady state or quasi-steady states can be de-tected by observing the magnitude and instantaneous frequency of a complex signal. After discretization, if \S(t)\ « const, and Ui(t) « uQ (2.39) then the solution is near the steady state. From the paradoxes below it becomes evident that it is possible to create scenarios in which the spectral content of the magnitude \S(t)\ equals the system frequency o>0 while the instantaneous frequency is zero. Physically, this could represent a 60 Hz component on a DC system and be falsely interpreted as steady state. In A C systems, however, complex signals are analytic in the sense that they satisfy the Cauchy-Riemann conditions for differentiability, and in [15] it is shown that the spectral content of ejv^ is always higher than that of \S(t)\ for analytic signals. Hence, the above approximations are a good indicator for steady or quasi-steady states in A C systems. This section has demonstrated that, for the detection of the transient state of a network, complex signals are very helpful. They provide this information easily with their magnitudes and phase angles. If this information were to be extracted from real signals, a Short-Time Fourier analysis would have to be conducted to obtain the spectrum of the signal. This process, to be performed at every simulation time step or so, is quite a computational burden, not to mention that this has to be done not only for one signal but for every variable (voltage, current, flux, etc.) of the system. 2.3. Spectral Properties of Complex Signals 31 2.3.3 Peculiarities Regarding Complex Signals The application of complex signals for solving a system may produce results that at first seem contradictory [15]. Since these phenomena are intrinsic properties of complex signals, they are briefly demonstrated with the following example. y 1 1. 1 1 1 1 1 1 1 1 -600' ' 1 1 1 1 1 1 1 1 ' 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 time t [s] time t [s] Figure 2.5: Peculiarities regarding complex signals E x a m p l e 2 .2 Let us consider a real signal that is composed of two different frequencies s(t) = S\ cos [LOI t] + s2 cos [OJ2 t] so that the corresponding complex signal is: S(t) = S1(t) + S2(t) = S i e i W l 4 + S2ej^1 = \S(t)\e^{t) where the two amplitudes S\ and s2 are constants and u>\ = 377 rad/s and u2 = 1131 rad/s are positive. The signal could represent a 60 Hz signal with a superimposed third harmonic. Its spectrum is thus composed of four peaks (Dirac pulses) at uj\ = ± 3 7 7 rad/s and uo2 — ± 1131 rad/s. Solving for the phase angle and amplitude, S\ sin [UJI t] + s2 sin [<x>21] ip(t) = arctan . . L S i cos [LOI t\ + s 2 cos [u>2 t\ 1 \S(t)\ = \Js\ + s2 + 2 S i s 2 cos [0J2 t — cuit], 2.3. Spectral Properties of Complex Signals 32 the instantaneous frequency is obtained by building the derivative of the phase angle: With different values for the amplitudes, we observe the following four characteristics: 1. Instantaneous frequency may not be one ofthe frequencies in the spectrum. 2. If the spectrum is composed of distinct lines at only a few frequencies, the instantaneous frequency may be continuous and range over an infinite number of values. 3. Although the frequencies of both superimposed signals can be positive, the instanta-neous frequency may be negative. 4. If the signal is bandlimited, the instantaneous frequency may go outside the band. These phenomena are illustrated in the two graphs in Fig. 2.5. In (a) Si = 0.2 p.u. and s 2 = 1.0 p.u. The instantaneous frequency is continuous and ranges outside the bandwidth. In (b) Si = 1.0 p.u. and s 2 = 0.5 p.u. Although the negative spectrum is suppressed, the instantaneous frequency may become negative. Chapter 3 Numerical Integration 3.1 Integration w i t h C o m p l e x Signals The previous chapter introduced two new ways of representing system variables: complex signals and dynamic phasors. Their properties were analyzed, and advantages and disadvan-tages were pointed out and compared to other variable representations. Their application to network simulation is described in this section. The integration of differential equations is a central part of every simulation program. Special care has to be exercised as neither integration nor differentiation come natural to a digital computer. Therefore, certain rules must be established to translate an integral or differential expression to its discretized equivalent that can be interpreted by the computer. Over the years, many sophisticated integration rules have been developed. An extensive list of first to sixth-order algorithms including the algorithm "families" by Runge-Kutta, Adams-Bashforth, Adams-Moulton and Gear is presented in [70]. Considering this tremendous amount of possibilities, it appears rather ironical that it is the simplest rules that have proven most powerful for applications in power system simulation. Among these few rules are the backward Euler algorithm, trapezoidal algorithm and, in some applications, the fourth-order Runge-Kutta-Fehlberg algorithm. Due to their simplicity, accuracy and numerical stability, only the first two rules are used in the context of this work, even though the author believes that the techniques presented here could be applied to any other integration rule as well. To begin with, it should be noted that the only modifications necessary to accommodate 33 3.1. Integration with Complex Signals 34 a simulation with complex signals is to allow the algorithms that are currently implemented in power system simulation programs to perform complex arithmetic. For the presentation of the backward Euler and trapezoidal algorithms, we again refer to the circuit in Fig. 2.1. Using complex signals, the differential equation is given by (2.32). Sometimes it is more convenient to view this equation in the state-space domain, where it would be formulated as jl(t) = -L~L RI(t) + L-1 V(t). (3.1) Rewriting it as an integral equation, we seek a discretization for t I{t) = j - L - 1 R I ^ + L-1 V{T) dr (3.2) o Applying any of the aforementioned integration rules, this otherwise continuous equation is discretized into a form that can then be implemented on a digital computer. 3.1.1 Complex Backward Euler Algorithm Let us first introduce the function F(t) for the integrand. From (3.2) we thus have F(t) =-L--1 RI(t) + L'1 V{t). (3.3) With this function, the backward Euler algorithm gives the following discretized formula, as can be verified in almost any textbook on numerical integration: I(t)=I(t-At) + AtF{t). (3.4) Backsubstitution of the above function gives the complex current: I(t) B=E (R+ 1 V(t) + H(t-At), (3.5) where H(t-At)=^R+-^j j^I(t-At). (3.6) The vector H(t — At) represents the so-called history values from the previous time step. In this document as well as in other literature on numerical power system analysis, the discretized integration rule is broken up into two equations: One expresses the current 3.1. Integration with Complex Signals 35 injections into the circuit as a function of voltage and a history vector. The other gives a formulation for the step-by-step updating of the history vector. The first equcition (3.5) must be solved simultaneously for the entire system, while the history values are calculated independently for each model. Details of the general structure of a simulation procedure are explained in Chapter 5. 3.1.2 Complex Trapezoidal Algorithm The trapezoidal algorithm is often used due to its very good numerical stability. It has the ability to produce results that stay "near" the true solution even if too large time steps are used. The discretization of (3.2) with this rule gives J ( t ) T = R J ( t - A t ) + ^ ( F ( i ) + F ( * - A t ) ) , (3.7) where F(t) is the same function as before. Backsubstitution results in the following equations for the current injections and history sources: I(t) = ( R+^) 1 V(t) + H(t-At), (3.8) 2L\ and H(t-M) = ( R+ — j or V(t-At)-[R- — )I(t-At) (3.9) Except for using complex arithmetic, (3.5), (3.6), (3.8) and (3.9) are identical with the discretization for real signals. Therefore, the properties of the backward Euler and trape-zoidal algorithms in their complex form are not different from the properties of their real formulation. The properties of these rules regarding their numerical error and stability are investigated in greater detail in Section 3.2. 3.1.3 Relaxed Trapezoidal Algorithm Due to the orthogonality of their real and imaginary parts, complex signals S(t) can be integrated with the backward Euler algorithm and trapezoidal rule. However, we have not gained any advantage so far over what we were able to do with only real signals: Instead, by introducing complex arithmetic, we now occupy twice as much memory as before and require 3.1. Integration with Complex Signals 36 about double the number of operations per time step. The advantage of using complex signals comes from "relaxing" them in such a way that the integration is then carried out on a slowed signal. As was shown before, dynamic phasors are identical to phasors at steady state. Therefore, integration at steady state could be carried out with At —)• oo and still yield the correct answer because all dynamic phasors remain constant. Rewriting (2.33) in the state-space form, we have jt®[I(t)} = {-L-lR-juJol) ^{I^ + L-1 &{V{t)} (3.10) Here, A = — L~xR — j cu0l is the system matrix or, in this case, the RL-series branch matrix. Applying the trapezoidal rule to this differential equation yields after backsubstitution with complex signals I(t) ^ R ^ ( l - ^ ) ~ l L - 1 V(t) + H(t-At), (3.11) and -1 r i f ( t - A t ) = e ^ ° A t - -Tj-AJ ~L~l V(t-At) + (l + ^A^j I(t-At) (3.12) Let us give a physical description for this integration formula. For real variables, the trapezoidal algorithm can be visualized in a two-dimensional plane where the true integral would represent the area that is enclosed by the signal and the time-axis on two sides and the projection ofthe signal onto this axis at t^—At and ti on the two other sides. The trapezoidal approximation is based on the assumption that the signal is a straight line between these two projected points. An example for this approximation is shown in Fig. 3.1(a) for the integration of the signal s(t) = 1 p.u. sin [u01] between the times ti — At and t\. The relaxed trapezoidal integration rule derived in this section can be visualized in a three dimensional space spanned by the real axis, imaginary axis and time axis. A correct integral in these three dimensions would represent the surface limited by the trace of the signal on one side, the time-axis on the opposite side, and the projections of the signal onto the time-axis at ty-At and t\ on the remaining two sides. The assumption made for the above integration rule is that the dynamic phasor (see Definition 2.2) can be approximated by a 3.1. Integration with Complex Signals 37 straight line between these points. The complex signal corresponding to the above sinusoid is depicted in Fig. 3.1(b). The integration is indicated by the striped surface. Since the signal in this example is of constant amplitude and oscillates with co0 the relaxed integration yields the correct result independent of how large At is chosen. 3.1.4 Relaxed Convolution Algorithm In [58], Linnert investigates some properties of the relaxed trapezoidal rule and argues that the optimum integration rule would converge to I(t) = (R + JUJ0L)-1 V{t) (3.13) for At —> co, assuming that the step size can be made infinitely large at steady state. This requirement follows from the definition of dynamic phasors in (2.20) and (2.26), and the phasor equation (2.11). Linnert shows that (3.11) does not have this property but retains a history term. Therefore, he suggests to factorize the eigenvalues and applies the convolution theorem to solve the resulting differential equation, which then produces the desired property. Again starting with (3.10) and subtituting J{t) =exp[-At] ^[/(£)], (3.14) 3.1. Integration with Complex Signals 38 we obtain the equation yJ{t) = exp[-At]L-1@[V(t)]. (3.15) After writing the integral form, backsubstitution yields t @[I(t)} = exp[At] y exp [-AT] L ' 1 ®[V{T)) dr. (3.16) o Assuming that this integral has been evaluated until the time t—At, the current at this time is known as @[I(t — At) ] and can be used to solve for the current injection at time t. Then, t @[I(t)]=exp[AAt]@[I(t-At)] + exp[At] j exp[-A r] L 1 @[ V(r)} dr. (3.17) t - A t Applying the product rule, the integral can be expanded to: t j exp[-A T] L 1 9[ V(T) ] dr t-At = -A^expf-At] (L-l9[V{t)}-exp[AAt] L ~ L ®{V{t-At)^ + ... t ... +A-1 I exp[-Ar] L~X ®[V{T)} dr (3.18) t - A t where this last integral can be approximated as t J exp[-Ar] L-x9\V{r)\ dr / e x p [ - A r ] L - 1 ^ — — — i ii dr (3.19) t - A t t-At « - A - exp[-A.] (7 - exp[A At]) L - * [ " ( « ) ] - S W * - * ) I. ( 3 . 2 0 ) The exponential matrix is a square matrix that commutes with the branch matrix so that A_1exp[—At] = exp[—At] A - 1 . It is obtained with the Pade approximation also imple-mented as part of the M A T L A B routines. The terms in the above equations are now collected to express the current as a function of voltages and history values. Transformed back to give the complex current, this is: I ( £ ) R = v - ( l + ^ (l-exp[AA<])) A"1 L ~ L V(t) + H(t-At), (3.21) 3.2. Integration Errors 39 where the history values are determined by: H(t-At) = exp[AAt] ej u°AT I {t - At) + ... . . . + (exp[AA*] + ^ - ( l - e x p [ A A t ] ) ) eju}°Ai A~l L~Y V(t-At) (3.22) This equation was first developed in [39] for lossless elements. For balanced three-phase systems, a similar expression was found by Linnert [58], who called this integration rule SDLV3 and used it in conjunction with the space-vector transformation [21]. 3.2 Integration Errors The natural response of a continuous-time system is a time trajectory. Integration algo-rithms, however, produce a discrete series of points. The important question is: Do these points bear any resemblance to the trajectory of the real system? The answer to this ques-tion is provided by an analysis of the numerical errors, which can be divided into the study of accuracy, local and global integration errors as well as numerical stability. In this section, the backward Euler, trapezoidal, relaxed trapezoidal and relaxed convolu-tion rules of integration are analyzed with respect to their numerical error. While the back-ward Euler and trapezoidal rules of integration have already been thouroughly investigated in other literature, e.g. [63, 70], all references on the relaxed integration rules [39, 58, 59, 84] fail to give a rigid analysis of their numerical errors. Besides finding out if the relaxed trapezoidal and relaxed convolution rules are at all suitable integration methods for this research project, the purpose of such analysis is to determine the maximum allowable simulation step size for each of the rules, as the integration error must always stay below a certain tolerance margin. A l l these issues are addressed here. 3.2.1 Accuracy The accuracy of an integration algorithm is determined by comparing the frequency response of a simulated system to the response of the natural system. To assess the accuracy of the aforementioned integration rules, let us consider a pure integrator or, in electrical terms, a 3.2. Integration Errors 40 single-phase inductive branch. For a combined series resistance and reactance, the current-voltage relationship with the relaxed trapezoidal rule is given by (3.11). Setting R = 0, this equation becomes H T R A t L - i + H{t-At), (3.23) with the history values being i f ( t - A t ) = f * e . L - 1 V ( t - A t ) + A e ^ ° A t / ( t - A t ) . (3.24) 2 4- j w 0 A t 2 + j LOQAI Applying the Z-transformation's displacement property Z[f(t — kAt)] = z~kZ[f(t)] for discrete time series /(£) to the above equations, yields the relaxed trapezoidal integration algorithm in the Z-domain: JW = R 2 T 7 ^ L " I n z ) + ffW- ( 3" 2 5 ) w 2 + jcu0AtzL y ; 2 + joj0At z y j y ' In the Z-domain, we can now combine these two equations and construct the admittance transfer function Y(z) = I(z)/V(z) for this integrator: y , x R T R z -t- e . . l Z J ~ L (2 + j u ; o A t ) ^ - ( 2 - i w 0 A t ) e ^ o A t ^ - J We proceed similarly for the other rules of integration. The admittance transfer functions of these rules are listed and compared in Table 3.1. A detailed analysis of the numerical accuracy of the backward Euler and trapezoidal rules in comparison to Simpson's and Gear's second order rules was presented in 1989 by J . R. Marti [63] in the context of an improved integration method for the E M T P . Assuming a sinusoidal excitation V(t) — e^wt, the response in the time-domain is given by /(t) = Y(z) ejujt , with z = e J w A t . (3.28) The accuracy of the integration rules can be evaluated from a comparison with the exact frequency response of an ideal integrator in the continous-time system: Y(s) = l/(sL) with 3.2. Integration Errors 41 Table 3.1: Transfer admittance functions for an ideal integrator with the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algo-rithms. B E T R R T R At z + pjuo&t Y(~) - ^ K J L (2 + j LO0At) z - (2 - j UJ0At) ei[u">At R C V ] LOQL \ j U!QAt (z — 1) / s = jto. Thus, building the ratio Y(z)/Y(s), an exact integration rule should produce the value 1.0 for all time steps A l Notice that the reciprocal value of this ratio produces the frequency response of the corresponding ideal differentiator. Since the integration rules are applied only once to produce the transfer function Y(z), the ratio Y(z)/Y(s) is the accuracy of one integration step. Whereas the accuracy of the backward Euler and trapezoidal integration rules only depends on the step-size-to-period ratio At/T — Atu/{2ix), the accuracy of both relaxed algorithms also depends on the ratio of excitation frequency to nominal frequency LO/LOQ. Example 3.1 Suppose the frequency ofthe sinusoidal excitation is identical with the nom-inal frequency of the circuit, i.e. to = LOQ and z = exp[j u)QAt\. With this value of z, the admittance transfer functions (see Table 3.1) for both the relaxed trapezoidal rule and the relaxed convolution algorithm become: This value is identical with the exact solution for an ideal integrator, and the ratio Y(z)/Y(s) is therefore one. This means that, at steady state, both relaxed integration rules are 100% accurate. 3 . 2 . I n t e g r a t i o n E r r o r s 4 2 T h e r e s u l t s o f Y(z)/Y(s) f o r t h e f o u r i n t e g r a t i o n a l g o r i t h m s a r e d e p i c t e d i n F i g . 3 . 2 a n d F i g . 3 . 3 . T h e p l o t s s h o w t h e m a g n i t u d e s a n d p h a s e s o f t h e y ( z ) / Y ( s ) - r a t i o s as a f u n c t i o n o f t h e s t e p - s i z e - t o - p e r i o d r a t i o At/T w h i c h i s b o u n d b e t w e e n Atmin = 0 a n d t w o s t e p s p e r c y c l e , i . e . AtNyquist = T / 2 , t h e m a x i m u m p o s s i b l e s t e p s i z e a c c o r d i n g t o t h e N y q u i s t c r i t e r i o n ( 2 . 1 ) . N o t i c e t h a t t h e g r a p h s h a v e n o t b e e n d i s p l a y e d b e y o n d t h e l o c a t i o n o f t h e i r first p o l e . T h e p o l e s a r e i n d i c a t e d b y a c h a n g e o f s i g n i n t h e p h a s e p l o t ( i . e . v e r t i c a l l i n e s i n F i g . 3 . 2 ) . F i g u r e 3 . 2 s h o w s t h e n u m e r i c a l a c c u r a c y o f t h e r e l a x e d t r a p e z o i d a l a l g o r i t h m i n c o m p a r -i s o n w i t h t h e a c c u r a c y o f t h e b a c k w a r d E u l e r a n d t r a p e z o i d a l r u l e s o f i n t e g r a t i o n . S e v e r a l c u r v e s a r e p r e s e n t e d for d i f f e r e n t v a l u e s o f co. S i m i l a r l y , F i g . 3 . 3 c o m p a r e s t h e r e l a x e d c o n -v o l u t i o n r u l e w i t h b o t h t h e b a c k w a r d E u l e r a n d t r a p e z o i d a l a l g o r i t h m s . N e x t , l e t u s c a r e f u l l y d i s c u s s a n d i n t e r p r e t t h e d e p i c t e d g r a p h s . I n g e n e r a l , t h e f o l l o w i n g q u a l i t a t i v e o b s e r v a t i o n s c a n b e m a d e f r o m t h e c u r v e s : • I n t e r m s o f m a g n i t u d e , t h e b a c k w a r d E u l e r i n t e g r a t i o n a l g o r i t h m s h o w s a h i g h e r a c c u -r a c y t h a n t h e t r a p e z o i d a l a l g o r i t h m . H o w e v e r , d u e t o t h e i n c r e a s i n g m i s m a t c h i n t h e p h a s e a n g l e , t h e r u l e i s s u i t a b l e o n l y f o r r e l a t i v e l y s m a l l s t e p s i z e s . • T h e t r a p e z o i d a l r u l e a l s o s h o w s d e v i a t i o n s f r o m t h e e x a c t s o l u t i o n as t h e s t e p s i z e i n c r e a s e s . T h e s e , h o w e v e r , a r e c o n f i n e d t o t h e m a g n i t u d e . O v e r a l l , t h i s r u l e p e r m i t s l a r g e r s t e p s i z e s t h a n t h e b a c k w a r d E u l e r a l g o r i t h m . • T h e r e s u l t s o b t a i n e d w i t h t h e r e l a x e d t r a p e z o i d a l i n t e g r a t i o n r u l e r e q u i r e a m o r e d e -t a i l e d d i s c u s s i o n : A t s t e a d y s t a t e , co = COQ, t h e m i s m a t c h i s z e r o f o r a l l s t e p s i z e s . F o r co > u; 0 , t h e m i s m a t c h i s less t h a n t h e o n e o f t h e t r a p e z o i d a l r u l e , y e t i t g r o w s w i t h h i g h e r f r e q u e n c i e s u n t i l , a t co—>oo, i t b e c o m e s e q u a l t o t h e d e v i a t i o n s s e e n w i t h t h e t r a p e z o i d a l a l g o r i t h m . F o r co < co0, t h e d e v i a t i o n s b e c o m e m o r e s i g n i f i c a n t as t h e f r e q u e n c y d e c r e a s e s . A t co = | w 0 , t h e m i s m a t c h i s a l r e a d y l a r g e r t h a n w i t h t h e t r a p e z o i d a l o r t h e b a c k w a r d E u l e r a l g o r i t h m s . F o r a l l f r e q u e n c i e s , t h e m i s m a t c h o n l y o c c u r s i n t h e m a g n i t u d e , n o t i n t h e p h a s e a n g l e . • T h e r e l a x e d c o n v o l u t i o n a l g o r i t h m a l s o p r o d u c e s n o m i s m a t c h a t s y s t e m f r e q u e n c y . F o r co > COQ, t h e a c c u r a c y i s s i m i l a r t o t h e c h a r a c t e r i s t i c o f t h e r e l a x e d t r a p e z o i d a l r u l e , 3.2. Integration Errors 43 Figure 3.2: Numerical accuracy of the backward Euler (BE), trapezoidal (TR), and relaxed trapezoidal (RTR) rules of integration. 3.2. Integration Errors 44 Figure 3.3: Numerical accuracy of the backward Euler (BE), trapezoidal (TR), and relaxed convolution (RCV) rules of integration. 3.2. Integration Errors 45 even though it is slightly lower. For co < co0, the accuracy is generally better than with the relaxed trapezoidal rule, and for | co0 < to < co0 it is even better than all the other integration methods. However, below this frequency, the deviations are generally larger. Again, the mismatch only occurs in the magnitude. Although the low frequency behavior of both relaxed integration methods appear very inaccurate compared to the two other integration methods, the phenomenon is not at all surprising: While both the backward Euler and the trapezoidal rules become very accurate as co —>• 0, the relaxed rules were specifically designed to obtain zero error with to = co0. To achieve this property, a frequency shift from to = COQ to to — 0 was necessary. This frequency transformation must be taken into account for the calculation of the maximum step size on the right-hand side of the abscissa, according to the Nyquist criterion. While, in Figures 3.2 and 3.3, this limit occurs at At/T = 0.5, the above frequency shift would produce a maximum step size of A t a; . -T = V\ 1 ' 3 ' 2 9 ) 1 1 \to — COQ\ depending on the frequency. Wherever this limit is lower than the right-hand side of the abscissa, it was indicated with a star (*) in Fig. 3.3. Notice that the accuracy at the indicated points is generally higher than the accuracy produced by the backward Euler and trapezoidal rules at the abscissa limit. In practice, though, the low-frequency behavior of the relaxed integration algorithms is not very important because the instantaneous frequency in A C systems is generally larger or around system forcing frequency co0. Observing the previous figures, a quantitative statement about the largest possible in-tegration step sizes for the investigated integration methods is difficult, considering that mismatches occur not only in the magnitudes but also in the phase angles. The following analysis will provide more information about this topic. 3.2.2 Local Errors In order to determine the maximum permissible step size for an integration algorithm more accurately, it is necessary to consider the numerical error that is produced by applying the 3.2. Integration Errors 46 algorithm for one integration step. This local error is composed of a round-off error and a truncation error. Round-off errors are generated by a digital computer when representing real numbers in terms of binary data with finite precision. Since this precision varies with each numerical processor, round-off errors are disregarded in the following discussion. The truncation error derives its name from the integration rules that are based on the Taylor series (e.g., Runge-Kutta). The error is introduced when the infinite series is truncated after a number of terms. More generally, the term is used for the absolute mismatch between the result of one step with the integration rule and the exact solution. Thus, using the ratio Y(z)/Y(s) from the previous section, \Y(z) e = Y(s) (3.30) Again, assuming a sinusoidal excitation z = e J t j A i , the truncation error for the relaxed trapezoidal rule becomes w A t w 0 At + 2 tan [^(w — w 0 1 (3.31) Table 3.2: Truncation errors for the backward Euler (BE), trapezoidal (TR), relaxed trape-zoidal (RTR) and relaxed convolution (RCV) algorithms. B E e (uAt\ e x p [ j ( ^ ) ] V 2 J sin[*f*] T R < = (uAt\ cos[^ ] V 2 J sin[^ ] 1 R T R e = to At o;0 A t + 2 tan [^ (w-w0)] R C V e = w - w 0 w sin sin [f (w - w0)] uJ0 OJ0 ("ft) sin[wf ] The error functions of all studied integration methods are summarized in Table 3.2. Notice that all functions except the relaxed trapezoidal rule have a pole at At/T = 1. This implies a very important constraint for the use of large integration step sizes in conjunction with the relaxed convolution rule: The step size A t cannot exceed the duration of one cycle at frequency w. This result contradicts Linnert, who states in [58] that 3.2. Integration Errors 47 "[...] symmetrical systems that contain no zero sequence values can thus be simulated without error using the relaxed convolution rule with arbitrarily large step sizes." This statement is true only in the steady-state case where LO = LOQ because then the truncation error is indeed zero for all step sizes. Otherwise, one can easily construct cases where the relaxed convolution rule fails to give the correct response. With regard to Linnert's concern about the convergence of the relaxed trapezoidal rule to I(t) = (R + JLU0L)-1 V(t) (3.32) for A t —> oo (see Section 3.1.4), it can be verified that the rule does converge to this steady-state equation if the solution at t — A t is also at steady state. Because of the step size limitation of the relaxed convolution algorithm we do not further consider it in the following analysis of practical step sizes. Let us now look at the interval { A t / T | 0 < At/T < 1.5} more closely, in order to find some guidelines to determine a suitable simulation step size A t for a given frequency to. The graphs of the truncation errors are depicted in Fig. 3.4 for frequencies above system frequency, and in Fig. 3.5 for frequencies below system frequency. Since the error of the relaxed trapezoidal rule varies with different frequencies, a formula must be found to determine the optimum integration step size. As indicated in Section 2.1, the errors of the backward Euler and trapezoidal rules of integration are equal for all fre-quencies so that the maximum step width is determined according to (2.2): 7T At m a x (u; j) = - — (3.33) With a tolerated integration error of e = 3%, the resolution factor k for the trapezoidal rule becomes approximately kTR = 5, which is the recommended value for ordinary E M T P simulations [63]. The resolution factor for the backward Euler rule of integration must obviously be much larger and, as will be demonstrated in the next section, this method is therefore not well suited for regular integration. For the relaxed trapezoidal method, the maximum permissible step size for a given in-stantaneous frequency can be computed for the location where the graphs in Figures 3.4 3.2. Integration Errors 48 I I 1 I 0 0.5 1 1.5 A i / T [p.u.] Figure 3.4: Local integration errors for cu > uo0. o 0.5 1 1.5 A i / T [p.u.] Figure 3.5: Local integration errors for u> < COQ. 3.2. Integration Errors 49 and 3.5 cross the tolerated error margin €RTR = £• This location is given by w A t " = £ l (3.34) cu 0At + 2 tan [ ^ ( w - w 0 ) ] which must be solved for At . Assuming £ < 1, we can make a cubic approximation of the tangent function for small arguments, i.e. 1 3 tan x ~ x + - x . With this approximation, the formula for the step size calculation becomes: A+ I \ / 1 2 £ U J i i c { l - £ for W i > W 0 fo oc\ A t - M ) = \Ja^--of ^ = I 1 + E for U i < c 0 ( 3 ' 3 5 ) The crossover points of the local integration errors of the relaxed trapezoidal algorithm with the 3% error margin are indicated in the above figures. It can be clearly noticed that, within frequency deviations of ±2%, much larger step sizes can be used with the relaxed trapezoidal rule than with any other of the presented integration methods. 3.2.3 Global Errors In this section we consider the global effects of the local errors. As before, the impact of the round-off error is ignored and only the truncation error is studied. Whereas local errors are produced by integrating the system with one step from t 0 to t 0 4- At , the global error represents an accumulation of local errors by integrating the system for a number of steps. Using the step size At , we now integrate the system for one period T = 27r/u of the frequency co. The total number of integration steps is N = T/At. The accumulated global error is equal to the error per step, i.e. the local error, times the number of steps: e = Ne=^- (3.36) uAt K ' The resulting graphs are shown in Fig. 3.6 and Fig. 3.7. Looking at the global errors provides helpful additional information to our understanding of the investigated integration rules: 1. Contrary to the other rules of integration, the backward Euler method has an accumu-lated integration error even when used with extremely small step sizes. This clearly 3.2. Integration Errors 50 3.3. Numerical Stability 51 shows that this rule is not well suited for regular integration. However, due to its excellent numerical stability (see next section) and damping property [24, 51, 63], the method is used to avoid numerical oscillations when resuming the integration after a discontinuity. 2. For low frequencies, the relaxed trapezoidal rule requires much smaller step widths than the regular trapezoidal rule. Therefore, the latter should be used whenever u> < \OJQ. For higher frequencies, the relaxed trapezoidal rule generally produces less errors so that it allows the use of larger step sizes than the trapezoidal method. 3.3 Numerical Stability When analyzing the local error of an integration algorithm, the inaccuracy of only one step is considered, assuming that the results of all previous steps were located exactly on the analytical trajectory. This, however, is not true if the previous points were calculated with the same integration routine. Instead, initial inaccuracies caused by the local errors will accumulate throughout the rest of the integration. Although the global error gives a good indication as to how the local errors accumulate, it remains yet unanswered whether these accumulations stay bounded or unbounded to the original trajectory produced by the system. The stability criterion provides the answer to this question. If an integration algorithm is unstable, initially small errors propagate through the inte-gration and may become unbounded. If, on the other hand, it is stable, the accumulation of errors is clamped through partial cancellation such that the computed results tend to follow the low-frequency shape of the trajectory, while they may be inaccurate in following high-frequency "ripples". The standard method for testing numerical stability is to apply the integration algorithm to the first-order linear test equation ^-S(t) = XS(t), S(t0) = S0 (3.37) 3.3. Numerical Stability 52 with the solution cf>(S0,t) = S0ext (3.38) where S, So, and A may be real valued or complex. Complex values are allowed since higher-order differential equations can have complex eigenvalues. An integration algorithm applied to the system (3.37) produces a discrete-time system with a fixed point at the origin and a characteristic multiplier m [70]. The stability of the fixed point determines the stability of the integration algorithm. The integration algorithm is numerically stable if the fixed point is stable, i.e. \m\ < 1. The integration algorithm is numerically unstable if \m\ > 1. In the following, the previously discussed integration methods are investigated with respect to their numerical stability properties. 3.3.1 Backward Euler Algorithm Since the above test equation represents an unforced system, the forcing function in the backward Euler algorithm from (3.5) and (3.6) can be set to zero such that only the following term remains: S(t) = S(t - At) (3.39) The characteristic multiplier of the fixed point is the factor mBE = (1 — A A t ) - 1 , and thus the backward Euler method is stable whenever 1 > |(1 - A A t ) _ 1 | < 1 ( 1 - A At) | < |(»[AAt]-l)+j(3[AA*])| ^ < m\At] - l ) 2 + (3[AAt]) 2. The right-hand side of this last relationship describes a circle of radius one with its center being at A A t = 1. It means that, for any choice of A A t lying outside of this circle, the backward Euler rule of integration is numerically stable. The stability region is depicted in Fig. 3.8. 3.3. Numerical Stability 53 IS.:'/ ' i l i i i i i i i i E 1 i 1 12 Re[XAt] lillilSSl!PI18lllliSlHllHII!| 1 2 R e [ X A t ] Figure 3.8: Stability region for the backward Figure 3.9: Stability region for the trape-Euler integration rule (shaded area). zoidal, relaxed trapezoidal and relaxed con-volution rules of integration (shaded area). 3.3.2 Trapezoidal Algorithm The discretization with the trapezoidal rule was formally derived in (3.8) and (3.9). To apply it to the linear test equation, the forcing function is again set to zero, and it remains The characteristic multiplier is mTR = (2 + AAt)/(2 — A At) , and its fixed point is stable if |2 + AAt | < | 2 - A A t | |(2 + 5R[AAt]) + j(3f[AAt])| < | (2 - t t [AAt] )+ j (3 [AAt] ) | (2 + K[AAt]) 2 + (3[AAt]) 2 < (2 - K[AAt]) 2 4- (S[AAt]) 2 (3.42) (2 + $R[AAt]) < (2 - f t [AAi j ) 3?[AAt] < -5?[AAt]. This last condition is true for all AAt whose real part is negative. Hence, the stability region for the trapezoidal integration rule includes the entire negative half-plane of AAt, and is thus smaller than the stability region of the backward Euler algorithm. For the integration of damped linear systems, any A t > 0 can be chosen to produce bounded simulation results. The stability region is rule is shown in Fig. 3.9. 3.3. Numerical Stability 54 3.3.3 Relaxed Trapezoidal Algorithm With zero forcing functions, the expression for the relaxed trapezoidal algorithm in (3.11) and (3.12) simplifies to the following equation: 2 + XAt • S(t) = 2 - XAt (3.43) Hence, the method is stable if 2 + XAt J w At < 1. (3.44) 2 - XAt Because of | ejtjAt | = 1, the analysis of the numerical stability for this integration rule is identical to the derivations made in (3.42). Consequently, the stability regions of both trapezoidal rules are equal. 3.3.4 Relaxed Convolution Algorithm Setting the forcing function to zero, the relaxed convolution rule of integration (3.21) reduces to the following expression: S{t) = eXAtS(t- At). (3.45) Hence, the characteristic multiplier is mRCv = eXAt, and the integration method is stable when 1 > |e A A t | 0 > K[AAt] Similar to the ordinary and relaxed trapezoidal rules , the relaxed convolution integration method is numerically stable for all choices of At > 0, as long as the real part of XAt is negative. This, however, is always the case in linear systems with sufficient damping. The stability region is identical to those of the two trapezoidal integration rules (see Fig. 3.9). Chapter 4 Variable Integration Step Sizes This chapter describes a method to automatically adjust the step size A t of the integration algorithm according to the system dynamics. In previous chapters, we already mentioned variable step size algorithms that are commonly in use, e.g. fourth order Runge-Kutta [73] and Gear's second order integration algorithms [34, 31, 80]. These algorithms are based on a straight line approximation or polynomial expansion of the underlying oscillations and adjust their step size according to the produced truncation error. An important drawback of these methods is that the step size is continuously varied. This may be very practical if the system's differential equations were expressed in the state-space form, however, for equations in the nodal form, which is common to most E M T P -type applications, continuously changing step sizes is very inefficient. The reason for this inefficiency is that the nodal admittance matrix used for the network representation depends on the step size A t and must hence be recomputed each time the step size is changed. Since the time for rebuilding the matrix is much larger than the time required to advance the solution from time t — A t to t, this may considerably slow down the performance of the program. A practical step size control is therefore subject to several conflicting notions: While the step size must be sufficiently small to maintain accuracy, it should also be as large as possible in order to provide high simulation efficiency. In addition, the simulation step size cannot be altered too frequently because each change requires additional time for rebuilding the nodal admittance matrix. The vague and partially contradictory nature of these requirements has 55 4.1. Step Size Prediction 56 led to the development of an intelligent step size control, which is described in [42] and which uses fuzzy logic [81] and risk assessment tools [17], in order to provide an accurate and efficient solution algorithm. The various parts of this algorithm, which are concerned with predicting a new simulation step size, evaluating the reliability level of the prediction, and rating the efficiency of a step size change for the solution process, all contribute to the overall decision making process that determines if the present step size is maintained or changed to the predicted size. The underlying concepts for this decision making are explained next. 4.1 Step Size P r e d i c t i o n For an accurate solution, the simulation step size must be inversely proportional to the maximum frequency present in the system. This section describes how we can obtain the highest system frequency at every time step of a simulation, and how it is forecast for future time steps, so that the step size can then be reduced or increased accordingly. 4.1.1 Maximum and Minimum Frequency In Chapter 2 the concept of complex signals was introduced. One advantage of this repre-sentation of system variables is that the instantaneous frequency of a quantity S(t) can be approximated by using the derivative of the phase angle with respect to time. For practical purposes the discretization is sufficiently accurate. Comparing the instantaneous frequencies of each variable, the ap-proximate maximum frequency u>max(t) > to0 can be obtained at each instant. In order to optimize the step size, it is not advantageous to utilize the highest occurring frequency deviations of each step, but rather use an average of these maximum frequencies over a certain period of time, say 3/4 of a cycle at nominal frequency. The reason is that individual instantaneous frequencies may fluctuate considerably during the transient phase after a disturbance, where the step size is not supposed to change at every computed time. Vi(S, t) angle[S(t)] - angle[S(t - At)] At (4.1) 4.1. Step Size Prediction 57 Hence, the occurring frequency deviations must be sampled and evaluated to predict a new step size. This is done using linear regression to fit a straight line Ulinear(t). = Ua + OCa t (4.2) to the frequency samples. Since the samples are equally spaced at the A t currently in use, the regression formula [72] becomes NoF to, a = ^(^ + K{t + f)-KNoFf\iomax(tl) (4.3) NoF 4=1 where _ 6(NoF-2i + l) A t NoF (NoF - 1) (NoF + 1) ^ = t - (NoF-i) A t and where NoF is the number of samples taken. While the regression line represents a mean, linear distribution function through all samples, it must be shifted until it cuts through the highest sampled frequency. This procedure is necessary to ensure that we are tracing the maximum frequencies. With the resulting tracing line we are able to make a linear prediction about what maximum frequency can be expected within the next window of NoF samples. This frequency value, topredictedtmax, is then used to obtain a prediction for the optimum simulation step size. In addition, we also sample the lowest frequency deviation from power frequency, i.e. umin(t) < OJQ. The use of these values is important because the maximum step size after (3.35) is different for frequencies above or below the power frequency toQ. Applying the above procedure to the additional frequency samples, we can thus forecast to predicted,min-4.1.2 Magnitude Deviation With the above sampling of instantaneous frequencies, we have obtained a criterion to make a step size prediction. As an additional criterion, we will now also evaluate the deviations that occur in the magnitudes of complex signals. This way, a step size prediction is possible where 4.1. Step Size Prediction 58 the frequency information is not availabe (e.g. in D C systems). The magnitude deviation is approximated by: \S(t)\-\S(t-At)\ AAi(S, t) = arctan A t S„ (4.5) where Sn is the nominal value of the corresponding variable. The arcus tangent function is used to obtain a normalized magnitude deviation for all variables, which ranges from zero to 7 r / 2 . Similarly as with the above frequency prediction, the highest magnitude deviations A-Amax of each time step are sampled over 3/4 of a period. Instead of applying the linear regression method to forecast the highest deviation, the maximum sample is taken as this value. Hence, A A predicted = max [ AAmax } (4.6) The fuzzy rule for a step size control with this magnitude criterion is: "If the magnitude deviations are close to TT/2, the used simulation step size should be small. If the magnitude deviations are close to zero, the step size can be large." Mathematically, this is reflected by the fuzzy membership function for "The magnitude deviation is small", which is defined as -< 2 AApredicted , . _\ ^•magnitude V*'' I TX The function becomes ^magnitude — 1-0 if there are no deviations in magnitudes of all variables, and it is zero for AApredicted -> 7r/2. It is shown in Fig. 4.1. 4.1.3 Simulation Mode Selection The simulation program that resulted from this thesis research uses the backward Euler, trapezoidal and relaxed trapezoidal rules of integration. The backward Euler rule is only used immediately after discontinuities, with a fixed, small step size to avoid numerical os-cillations [51, 63]. It can therefore be left out from the discussion of self-adjusting step sizes. 4 . 1 . Step Size Prediction 5 9 Magnitude deviation is low Frequency gradient is low Frequency is around UJ0 • • • J l_ 1 0 TT / 4 TT /2 0 50 100 0 1 2 ^ p r e d i c t e d [rad] \aa| [% of W 0 / S E C ] Upredicted/Uo [p.U.] Figure 4 . 1 : Membership functions for step size prediction. In previous chapters, it was emphasized that the research work described here aims at overcoming the discrepancy between electromagnetic and electromechanical transients simulation. However, for the appropriate adjustment of integration time steps, the knowledge about the type of simulated transient phenomenon can be helpful so that two modes are introduced to distinguish between both types of transients: "If the magnitude deviation and frequency gradient are both low and the predicted frequency is around nominal frequency, then the program is in the transient sta-bility mode. Otherwise, the program should use the electromagnetic transients mode." The rule consists of three parameters: the magnitude deviation AApredicted, the gradient of the frequency samples, i.e. absolute slope \aa\ of the regression line, and the ratio of predicted to nominal frequency. To build this ratio, let us first define ^predicted — max \_CO predicted,max j 2 C J Q LOpredicted,min] • ( 4 - 8 ) The ratio is then given by copredicted (t)/wo- The membership functions of the three parameters are shown in Fig. 4 . 1 . In order to determine the membership value for the stability mode, the geometric average h'stab = v ^ l h-2 ^magnitude ( 4 - 9 ) is applied, where fii and [i2 are the two membership values for the "low gradient" and "frequency close to nominal" fuzzy sets, respectively. The membership value for the elec-tromagnetic transients mode is obtained by negating the value for the stability mode, as 4 . 1 . Step Size Prediction 60 explained in [81]: fJ'emtp = fistab — 1 ~~ t^stab- (4-10) These two membership degrees allow us to classify the behavior of a transient as an electro-magnetic or electromechanical phenomenon. 4.1.4 Step Size Prediction The step size control was designed to use the complex trapezoidal rule (see Section 3 .1 .2 ) during the presence of electromagnetic transients, whereas, for simulating quasi-steady and steady states, the integration is carried out with the relaxed trapezoidal method. Depending on which of the integration rules is used, the corresponding maximum step size can be determined with ( 3 . 3 3 ) or ( 3 . 3 5 ) . However, to achieve an optimum between accuracy and simulation speed, only half of this step size is used in the simulation, i.e. gopt = 2. Hence, AtTR = — (4 .11 ) Qopt Cpredicted, max A , 0 .61 I" {-"Jpredicted,max / . ., 0 \ & t R T R , i = - — . / - -g (4.12) Qopt V ^predicted,max ~~ ^0 | A i 0 . 5 9 / LOpredicted,min /. 1 0 \ A-tRTR.,2 = . ' — — T3- (4.16) Qopt V \tOpredicted ,min - to0\ where, during the stability mode, the lesser ofthe two relaxed trapezoidal step sizes is used. Thus, Atstab = min [AtRTRtl ; AtRTR>2]. ( 4 . 14 ) For the electromagnetic transients mode, the magnitude deviation is applied as an addi-tional constraint for the trapezoidal step size. Therefore, the electromagnetic transient step size becomes: Atemtp = min [AtTR J (1 - ^magnitude)}- ( 4 . 15 ) U>o By weighing the electromagnetics time step Atemtp and the stability time step Atstab with their respective membership degrees ^emtp and (J,stab, we obtain the recommended time step: Atpredicted — f^emtp Atemtp "f f-^stab Atstab ( 4 . 16 ) 4.2. Prediction Reliability 61 At each calculated time step, this formula provides the algorithm with a new step size forecast. Notice that the final step size is always limited by the user-defined boundaries AtTnin and /\tmax. 4 .2 P r e d i c t i o n Rel iab i l i ty The judgement about the reliability of the above step size prediction is based on three key questions: 1. Is the sampling quality sufficient? 2. Do the frequency samples deviate much from the regression line? 3 . Is the simulated transient phenomenon clearly of an electromagnetic or electromechan-ical nature? As a fuzzy rule, a high degree of prediction reliability is defined as follows: "If the sampling is good, if the frequency fluctuation is small, and if a transient phenomenon can be well classified, then the reliability of the step size prediction is high. Otherwise it is rather low." Introducing the membership degrees ^sample, Hfluct and /J,mode for the three criteria, respec-tively, the membership for high reliability can be obtained using scalable monotonic chain-ing [17]. This approach to fuzzy reasoning has the advantages of being faster than conven-tional fuzzy reasoning techniques, and of avoiding the problem of saturating the solution fuzzy set. The application of scalable monotonic chaining to determine the reliability of the step size prediction is schematically shown in Fig. 4.2. The composition of the three individual membership functions is described below. In order to prevent the algorithm from over-reacting due to sudden changes in this reliability measure, the membership values of every time step are accumulated and averaged over 3/4 of a period. The result from the averaging is the overall prediction reliability rating fi reliability 4.2. Prediction Reliability 62 frequency samples N reliability level Figure 4.2: Monotonic chaining method to assess the prediction reliability. 4.2.1 Sampling Quality The frequency samples are read into a vector with a length of at most 50 values. In order to improve the accuracy, a larger vector would be possible but would also cause longer processing times for computing the regression line, standard deviation and other related quantities. If the vector contains less than 50 sampled frequency values, the performance of the sampling is lower than with a full vector. This is reflected in a linear relationship for a high number of samples NoF where, as before, NoF is the actual number of samples. Using this as the only criterion for good sampling would create a difficulty: Even with a full vector of sampled values (NoF = 50), the sampling window W = NoF At could be very small, depending on the current step size At . Whenever such a situation occurs, the algorithm discards every second sample and remains with NoF = 25, now sampling the maximum frequency of only every second time step. This process is repeated until the sampling window contains three quarters of one cycle. Another linear membership function was created to reflect the influence of the 4.2. Prediction Reliability 63 sampling window: H4 = NoF At ; 0 < ^ 4 < 1 (4-18) 8 7T The two membership functions are combined utilizing the geometric average to give the total membership function iisampie = yjv-z f'4 f ° r the fuzzy set: "The sampling is good." 4.2.2 Frequency Fluctuation The reliability of the step size prediction is increased if the frequency samples do not deviate much from the linear regression line. For this assessment, it is sufficient to compare the samples of the maximum instantaneous frequency tomax(ti). The fluctuation is measured by computing the standard deviation of the samples, thus a = X N o F 2 \ J^^(umax{U) -Ulinearitij} • (4.19) \| ° i=l The associated membership function for "Frequency fluctuation is low" was chosen to be a straight line that goes through 1.0 if a = 0 and is zero if a > 20% of the predicted frequency Wpredicted- Thus, a fJ-fluct = 1 - f-z ; 0 < Hftuct < 1 (4.20) Kj.L predicted The function for low frequency fluctuation is depicted in Fig. 4.3. 4.2.3 Simulation Mode The step size prediction is also more reliable if the simulated phenomena can be clearly classified as either electromagnetic or electromechanical transients. This classification can be quantified by the two simulation modes \iemtp a n d Ustab, considering that one of them becomes 1.0 whenever a phenomenon can be clearly categorized, e.g. if pstab is either close to zero or to 1.0, a mode is well defined. Between these two values, the mode can be pictured as being in a transitional state. Mathematically, this relationship can be expressed by the following formula: Hmode = 4 (flstab ~ 0.5)2 ; 0 < \lmode < 1 (4-21) 4.3. Step Change Efficiency 64 The membership function for the fuzzy set "The program mode is well defined" is shown in Fig. 4.3. Frequency fluctuation is low Program mode is well defined 0 50 100 0 0.5 1 Std. variation a [% of fpredicted] Membership degree f i s t a b [p-u.] Figure 4.3: Two membership functions for prediction reliability. 4.3 Step Change Efficiency In this section, we describe the process that evaluates how efficient it is for the program to change its simulation step size. Two aspects are important for this analysis: 1. The ratio of the proposed step size to the currently used step size /\tprediced/At. This ratio should show whether the difference in step sizes is big enough to encourage a change of step size. If the difference is too small, such a change would not be justified. 2. The ratio of time required to rebuild the nodal admittance matrix and the time required for the solution of one time step TrebuM/Tsoiution. The overall efficiency rating is composed of two individual efficiency levels for the above two aspects. They can be represented by the two fuzzy sets \xsiev and p delay as explained below. Similarly, as for the prediction reliability, monotonic chaining is used in order to determine the efficiency coefficient. A sketch of this process is shown in Fig. 4.4. As was done with the prediction reliability, this efficiency rating is calculated at every time step 4.3. Step Change Efficiency 65 efficiency level Figure 4.4: Monotonic chaining method to assess the step change efficiency. and accumulated. The average over 3/4 of a period is taken as the overall efficiency coeffi-cient (^efficiency 4.3.1 Step Ratio Generally, the efficiency for applying a new step size is increased if the ratio between the proposed and currently used step size rs = Atprediced/At is large. More specifically, for ratios rs greater than 1.0, the efficiency was defined to grow exponentially. Thus, 1.2 for rs > 1, (4.22) until the membership degree is 1.0 at rs > 10. For step ratios rs slightly below 1.0, the stepsize need not be immediately changed as the solution may still be reasonably accurate. As a result, the efficiency estimate for changing to the new time step must be very low. However, if the ratio drops below a certain margin Qmin/Qopt, the efficiency membership degree for a step size change should be 1.0. This is reflected in 2 for r s < l , (4.23) where the margin is defined through gmin = 10. The membership function for fuzzy set "The step ratio is high" is depicted in Fig. 4.5. Pstep — 10 instep Qmin Qopt 4.3. Step Change Efficiency 66 4.3.2 Time Delay Next, we introduce a time delay that will prevent the program from rebuilding the matrix too frequently, a process that requires much computation time and may considerably slow down the program's performance. The evaluation of this efficiency is based on two time values: The time required to rebuild the nodal admittance matrix rre(,Ujw, and the time needed to complete the solution of one time step rsoiution. The program automatically measures both durations when it first sets up the network matrix and then computes the solution of the first time step. But how many time steps n are needed so that TTebund is negligible compared to nTsoiution? If we choose a small tolerance £ = 0.005, then we can assess n = T r e b u M (4.24) £ ^solution This value represents the number of steps that must be computed with the old step size At until the time for rebuilding the matrix appears negligible. In addition, we must assume that the time delay nAt should always be constant, independent of what step size is used. Therefore, nAt = nAtprediCted = constant, and n = n A * . (4.25) ^''predicted In order to rate the efficiency, the number of steps K since the last step change is counted and compared with the number of required steps n. The resulting membership degree p, delay < 1 is obtained with H delay = ^ for Ts > 1 (4.26) For rs < 1, the algorithm automatically sets p, delay equal to 1.0 because, if the predicted optimum step size is significantly smaller than the one presently in use, the program should react immediately (with no time delay) and reset the step width to a smaller size. The membership function for the fuzzy set "The time since the last step change is high" is shown in Fig. 4.5 for n = 12,000. 4.4. Decision Making Process 67 Step ratio is high Time delay is high IO"1 10° 101 0 6,000 12,000 Step ratio rs Number of steps K Figure 4.5: Membership functions for efficiency evaluation. 4.4 Decis ion M a k i n g Process So far, we have explained the techniques used to forecast an optimum step size, to evalu-ate the reliability of this prediction and to obtain an efficiency level for using the predicted step instead of the present one. Defuzzification is now required to arrive at a final—crisp— decision: "use the predicted step size" or "maintain the current step size". This is accom-plished by the application of the inference technique schematically depicted in Fig. 4.6. The decision is based on yet another fuzzy set, which determines the likelihood for using a new step size. The fuzzy rules for this set are: "If the prediction reliability is high, the likelihood level for a new step size is increased," and "If the the step change efficiency is high, this likelihood level is also increased." Both rules are evaluated, again using scalable monotonic chaining, after the reliability and efficiency levels are added up as shown in Fig. 4.6. The resulting membership value ^change depends on the function that is chosen for the fuzzy set. While the fuzzy sets for the reliability and efficiency ratings have been purely linear, the membership function for the combined level is now modified by a hedge, which we call aggressiveness A. With this hedge, the program user can influence the tendency of the program to accept step size changes faster or slower: If the user is particularly interested in a detailed represen-4.4. Decision Making Process 68 prediction reliability step change efficiency hedge' Defuzzification step change certainty [ < 0.6V yes | V y \ no old At new At Figure 4.6: Defuzzification tation of electromagnetic transients, he would set A to a lower value, so that the program maintains smaller step sizes longer after a disturbance. If he instead wishes to speed up the simulation process after a disturbance, because his interest are the dynamics associated with machine rotor swings, the factor would be close to A = l. The membership degree for the fuzzy set "Likelihood of a step change is high" is then determined by ( (A reliability ~t~ (^efficiency where j3 = 1.2 — 0.4 A is the exponent to represent the hedge. If the likelihood of a new step size exceeds the value (iChange — 0.6, the program accepts Atpredicted as its new step size and rebuilds the nodal admittance matrix before continuing with the simulation. Otherwise, the simulation continues with the step size in use. The performance of this step size control algorithm has been demonstrated with the practical case studies in Chapter 9. A detailed discussion of these results is also presented in [42]. (4.27) Chapter 5 Structure of the Simulation Process While the previous chapters have already described some details implemented in the solution process of the simulation program, the overall structure has not yet been presented. This is therefore the subject of this chapter, which also covers some additional details about the way in which the equations are being interfaced in the algorithm. 5.1 Overall System Equations In general, an electric power system can be described by a set of differential equations and another set of algebraic equations. Since higher-order differential equations can be reduced to ordinary first-order differential equations, they can be expressed in their continuous form where g(-) is a function of the state variables X(t) and the forcing functions U(t). This equation is then discretized as part of the solution algorithm. Most stability programs use this representation to build the system equations [28, 55]. Since the transmission network is usually assumed to be at power frequency, it can be represented as an algebraic equation of the form which can be easily interfaced with the discretized form of the above equation. The use of (5.1) is advantageous for an accurate representation of continuous nonlinearities, which is d di X{t)=g(X,U,t) , (5.1) I(t) = YV(t), (5.2) 69 5.1. Overall System Equations 70 the reason why this part usually contains the machine and plant model equations in stability studies. The form is not well suited to represent discontinuous nonlinearities such as a changing topology through circuit breaker operation, because this would require a complete restructuring of the state-space equation. In stability studies, it is therefore assumed that all switching action occurs in the transmission network, thus affecting only the algebraic equation (5.2). Most electromagnetic transients programs discretize the differential equations at the branch stage, so that they become algebraic equations with a history term H(t — At), thus I(t) = YV(t) + H(t - At) (5.3) This representation is comparatively flexible to incorporate changes in the topology. How-ever, continous nonlinearities must be represented as piecewise linear approximations. This is admissible because electromagnetic transients programs use very small time steps so that the changes caused by nonlinear characteristics are typically small. For the new simulation program, it was therefore decided to use both representations (5.1) and (5.3), thus allowing for more flexibility during the modeling of individual network components.-While the algebraic form is probably better for the modeling of network com-ponents, the state-space form may be of advantage for controller models, because they can then be directly imported from stability programs. Figure 5.1 demonstrates how the system equations are separated into nodal and state-space equations, and with what integration rule they are usually discretized. With these equations, we can express the current injections and forcing functions as functions of the state variables and network voltages so that I(X,t) = Y V(t) + H(t - At) (5.4) X(t) = g(X,V,t) (5.5) In order to discretize this last equation, we apply the trapezoidal rule (3.7) on dynamic phasors to obtain the functional form of the relaxed trapezoidal rule: 9[X(t)) R= ®[ X(t-At)} + y (@[X{t)} + ®[X(t-At) ]\. (5.6) 5.1. Overall System Equations 71 Electrical synchronous machine equations: State-space equations, discretized with the Trapezoidal Rule K - o ) . Mechanical synchronous machine equations: State-space equations, discretized with the Trapezoidal Rule K = o). A C network equations: Nodal admittance matrix, discretized with the Relaxed Trapezoidal Rule (con = 377 rad/s). Figure 5.1: Discretization of power system components Backsubstitution to complex signals yields: (l + ± Q ) X ( t ) R=R (1 - | n ) e x p [ f i ] X ( t - A t ) . . . A t ^X(t) + X{t-At)e^p[n} ), (5.7) ' 2 where Q = j At diag[cjn] is a diagonal matrix whose elements consist of the product of j At times the nominal frequency of each corresponding state variable. The values for X(t) are obtained with (5.5) while X(t — At) is still available from the previous time step and need not be re-ccomputed. Let us introduce F(V,X,t) = I(X,t) -YV(t) - H(t- At) (5.8) G(v,x,t) = ( i + | n ) x ( t ) - ( i - | n ) e x P [ n ] x ( t - A t ) . . . . . . - ~ (s(X,V,t)+s(X,V,i-At)exp[ft]) . (5.9) At the solution point F(V,X,t) = 0 (5.10) G(V,X,t) = 0. (5.11) Notice that these equations are, in general, nonlinear. The E M T P does not require an iterative solution to nonlinear equations because of the small time steps used. A good prediction of the respective quantities is usually sufficient [24]. With larger simulation time 5.2. Program Structure 72 steps, a simultaneous predictor-corrector solution can be found with Newton's method, where the state and voltage vectors are iteratively improved according to the following equation for the (k + l)st iteration [55]: \ i r T ^ / J A I r A T T V J A I (5.12) " v(t) -• v(t) - 1 ' AV(t) " X(t) k+1 X(t) k AX(t) The improvements are obtained from F(V,X,t) G(V,X,t) dF dV dG dV dF aX aG dX ' AV{t) ' - k AX(t) (5.13) Convergence is quickly achieved with this method, provided the starting values are suffi-ciently close to the final solution. These are obtained by linearly extrapolating the dynamic phasors of past values as explained below. When nonlinear elements are present, some parts of the Jacobian matrix need to be continuously updated, although there are certain tech-niques to skip the updating from time to time to obtain a faster solution, e.g. with the "dishonest Newton method" [28]. 5.2 P r o g r a m Structure The structure of the new simulation program is shown in Fig. 5.2. It can be divided into three main parts, which are: reading of input data, steady-state initialization, and transient simulation. Optionally, the post processing of output data could be added to the schematic. Input data must be provided to accurately describe each network component model for a power flow, electromagnetic transients and stability simulation. In addition, the program user must provide minimum and maximum step sizes Atmin and Atmax, respectively, as well as the end time of the simulation tend and the aggressiveness factor A (see Section 4.4). Power flow and steady-state in i t ia l iza t ion . Before the beginning of a time-domain simulation, the program performs a balanced three-phase power flow calculation. The results of this calculation are then used to obtain the initial voltage and state variable vectors at time t 0 = 0, according to (2.27). At the same time, the durations TTebuM and T s o i u t i o n for 5.3. Time-Domain Simulation 73 < ^ Start ^ > Read input data / Setup \ frequency-domain^ ^ matrices / Power flow or steady-state initialization Output data ( Stop ) Figure 5.2: Flowchart: Program structure constructing the nodal admittance matrix—or Jacobian matrix—and for solving the Newton iteration, respectively, are measured to initialize the automatic step size control as explained in Chapter 4. 5.3 T i m e - D o m a i n S imulat ion A simplified flowchart of the time simulation is depicted in Fig. 5.3. It consists of a loop that is being executed until the simulated time span exceeds the user-specified end time tend. In the following paragraphs, some of the details of this time step loop are explained. 5.3.1 Prediction of Voltages and State Variables The use of Newton's method as an iterative correction to an initial guess requires the pre-diction of the voltage and state variable vectors. To allow fast convergence, great care must 5.3. Time-Domain Simulation 74 < Start here > Predict node voltages and state variables for /. Compute forcing functions at t. Compute model history values and collect to nodal history vector. Solve the system at time ffor unknown node voltages and current injections. For each discontinuous model, check if limiting values have been violated. Determine a new step size A t . Rebuild the nodal admittance matrix with new topology. Interpolate for the time tx of this violation, and set a discontinuity flag for this instant. For each model store computed voltages and currents for calculation of new history values. Rebuild the nodal admittance matrix with Af ... . Set tn = t, Af = Af O ' new and t=t0 + At Stop Figure 5.3: Flowchart: Time step simulation loop 5.3. Time-Domain Simulation 75 be exercised with the prediction of variables, and several prediction formulae have been pre-sented in the past [24, 54, 88]. For the applications presented here, linear extrapolation was sufficient. Similar to the prediction in stability studies [55], it is carried out on dynamic phasors, thus 9\SP(t + Atp)) = 9\S(t)]-z*E9[S{t- At)] , (5.14) where @[Sp(t + Atp)] is the predicted variable for the predicted new time step Atp. Using (2.20) to transform to complex signals, we have Sp(t + Atp) = A t ^ A t p e ^ A t " S(t) - ^ e W A t + A t f ) ) 5 ( t _ Aty ( 5 1 5 ) Here, u>n is the nominal frequency of the variable that is being predicted. For instance, for voltages in an A C transmission system un = uo, while state variables associated with the machine models are DC quantities so that con = 0 there. At moments of discontinuity, the voltages of the predominantly inductive transmission network may change instantaneously whereas the fluxes cannot change at an instant. There-fore, the vector containing the nodal current injections prior to the disturbance I(X,t-) is used to predict the voltages after the disturbance by solving I{X,t-)=YV{t+) + H(t+-At), (5.16) where the changed topology has been realized in the nodal admittance matrix Y as well as in the history vector. 5.3.2 Forcing Functions Forcing functions are the source functions that continuously "drive" the power system. So far, three different types of sources have been considered: 1. Ideal voltage sources 2. Ideal current sources 3. Dependent current sources 5.3. Time-Domain Simulation 76 While the functions ofthe ideal source types are determined by the program user, the current injection of a dependent source depends on the actual values of the state variables and node voltages. The term "source" for the last type may be misleading, as it could be, for instance, a model of an inductance represented as a flux linkage / ( ^ t ) = 1/L \I>(£) with a voltage equation 4f(t) = V(t) in state-space form. This is not an active source although, technically, the program understands it as such. The benefit of this way of modeling is that it allows the adequate representation of non-linear model characteristics in the form I = f(X,t) where / can be any nonlinear function. This is realized in the model of a synchonous generator (see Chapter 7) where the terminal current injections into the network are nonlinear functions of the stator and rotor fluxes, and of the position of the rotor with respect to the d and q axes. While the source functions of the ideal sources can be directly computed at the new time step, the dependent sources use predicted voltage and state variable vectors to produce an initial approximation of their current injection. The source functions are then written into the vector I(X, t) at the corresponding nodes. At the nodes without any of the above source types, the vector is zero. In addition to calculating the source functions, the computation of the forcing functions also includes an update of the state variable derivatives X(t). 5.3.3 Discontinuity Checks There are two situations in which the routine investigates if there has been any sudden disturbances or discontinuities during the last time step: First, it must be determined if there have been any user-scheduled discontinuities, such as an opening or closing of a circuit breaker or a fault somewhere in the network. This, of course, must be considered before the network solution part is carried out. After a solution has been obtained, the algorithm checks for any limit violations that may have occured, e.g. the breakdown voltage across a surge arrester. Notice that these limits can be nonlinear functions. Next, we shall explain how discontinuities are treated with this algorithm: Suppose a discontinuity occurs at tx, where (ti — At) < tx < t\. In case of a limit violation, t x > i is 5.3. Time-Domain Simulation 77 determined by linear interpolation. This may not be very accurate, particularly if large time steps are used. Suppose tX ) 1 is a not so accurate first guess for the moment of discontinuity, and also suppose that (ti — At) < tX}i < tx < tx. In this case, the program computes the solution at t = tx>\ without introducing the disturbance. Afterward, linear interpolation is again used to locate tx. Since the interval [ t x , i ,^ i ] is smaller than the original interval \{t\ — At) , ty], the next interpolation at tX ] 2 is likely to be more accurate. This way, the algorithm breaks up the above intervals as often as needed to be sufficiently close to the actual moment t = tx of the disturbance. Then, the simulation algorithm uses a backward Euler step with the user-defined min-imum step size Atmin. Since there are always two checks for discontinuities per computed time step, they may intercept each other so that the algorithm proceeds with backward Euler steps until the solution has become "smooth"—another fuzzy set that has not been explained in Chapter 4. This way, the program can overcome the difficulty of subsequent flashovers during the switching instant, pointed out in [29] as an unresolved issue. After the transition from the backward Euler to the trapezoidal rule of integration, the automatic step size control described in Chapter 4 comes into effect, until the solution process is again interrupted by a discontinuity. 5.3.4 C o n s t r u c t i o n o f t h e N o d a l A d m i t t a n c e M a t r i x The nodal admittance matrix is obtained with the following three steps: 1. Construction of the open-circuit admittance matrix, 2. Node collapsing for null-impedance connections, 3. Node reduction and sanity checks. The idea for this algorithm was formally introduced in 1982 by F. Alvarado [1], while some additional features had been incorporated by the author during the implementation of this algorithm. The main ideas are explained below. 5.3. Time-Domain Simulation 78 Construction of the open-circuit admittance matrix. All component models con-tribute their respective branch admittace matrices to the open-circuit admittance matrix. At this stage it is assumed that all switch models and null-impedance connections are open or not connected. Notice that the ground node is still part of the matrix and has not been eliminated up to this point. A systematic way of adding a model's branch matrix Y to the open-circuit admittance matrix of the whole system is schematically presented in Fig. 5.4 [24]. K K K rn, m. rn a D c a b c K m b rn [Y] -[Y] -fYl fYl Figure 5.4: Open-circuit admittance matrix Node collapsing for null-impedance connections. In this new simulation program, there exists a model subtype "null-impedance wire", which is a impedance-less connection between two nodes. Ideal switches that are closed adopt this subtype, as do purely inductive branches in DC networks during the steady-state initialization. This null-impedance con-nection is realized by collapsing the two connected nodes into one single network node, after the open-circuit admittance matrix has been obtained. The row and column of the duplicate node are eliminated only after a copy of the row has been retained for the calculation of the branch current through this impedance-less connection. All other node collapses are then also performed on this copy. This procedure ensures that the branch current is correct, 5.3. T i m e - D o m a i n S i m u l a t i o n 79 e v e n i f s e v e r a l n u l l - i m p e d a n c e w i r e s a r e c o n n e c t e d i n p a r a l l e l o r a n y o t h e r c o n f i g u r a t i o n . I n p r a c t i c e , h o w e v e r , s u c h c o m b i n a t i o n s a r e r a t h e r u n u s u a l . Node reduction and sanity checks. T h e n o d a l a d m i t t a n c e m a t r i x s t i l l c o n t a i n s t h e g r o u n d n o d e ( s ) a n d i s , t h e r e f o r e , s i n g u l a r . A p r o c e d u r e w a s i m p l e m e n t e d t h a t s u b s e q u e n t l y d e t e r m i n e s i f a s u b n e t w o r k i n t h e m a t r i x is s i n g u l a r a n d i n v e s t i g a t e s p o s s i b l e r e a s o n s f o r t h i s c o n d i t i o n . I f a f l o a t i n g n e t w o r k is d e t e c t e d , t h e n t h e p r o g r a m a u t o m a t i c a l l y s e l e c t s a g r o u n d n o d e w i t h z e r o p o t e n t i a l as a d e f a u l t v a l u e . T h e c o r r e s p o n d i n g r o w a n d c o l u m n a r e t h e n e l i m i n a t e d a n d v o l t a g e s a r e e x p r e s s e d w i t h r e s p e c t t o t h e c h o s e n g r o u n d . A t t h e s a m e t i m e , t h e p r o g r a m p e r f o r m s s a n i t y c h e c k s f o r a n y p a r a l l e l v o l t a g e s o u r c e s o r s e r i e s c o n n e c t i o n s o f c u r r e n t s o u r c e s . Part II Component Model ing 80 Chapter 6 Transmission Line In earlier chapters, we presented ways to model an RL-series branch for different integration rules and for different frequency ranges. The methods for linear capacitive elements are analogous. While most power system equipment can, at least in parts, be represented with combinations of these linear branches, there are other power system components that require different modeling techniques. Cables and transmission lines, for instance, introduce time delays for the current and voltage waves which travel from one end to the other. This time delay cannot be modeled with a combination of RLC-branches alone. In this chapter, we therefore discuss modeling aspects for transmission lines. Models for the electromagnetic transients simulation of cables or transmission lines have been well developed over the past three decades, but even today, they present an area for more research. The Department of Electrical and Computer Engineering at the University of British Columbia alone has come forth with five line and cable models for electromagnetic transients studies in the past 20 years [9, 60, 61, 62, 64, 68]. However, while research on line modeling has mainly focussed on very detailed models for the study of electromagnetic transients, the quasi-steady state, nominal 7r-circuit used in stability analysis has hardly changed at all. While switching between electromagnetic transients line models and their equivalents for stability analysis is theoretically possible, it is not always practical, as both model types may have a different internal topology. Whenever models are exchanged during a simulation, this difference in topology may cause rapid transients which, in turn, can create unrealistic 81 6.1. Derivation ofthe Model Equations 82 behavior ofthe simulated network (e.g. tripping a protective device) [54]. A consistent model for all time frames is therefore preferable. However, a symbiosis of both types of line models is not as straightforward as with other models: The quasi-steady state, nominal 7r-circuit used in stability studies is not sufficiently accurate to represent wave propagation in electromagnetic transients studies. On the other hand, line models preferred in electromagnetic transients programs, such as frequency-dependent or constant parameter line models, require a simulation step size of less or equal to the lowest wave travel time. With this constraint, EMTP line models are not well suited for stability analysis. In this chapter, a new modeling approach is proposed to overcome the step size constraint for EMTP transmission line models. In order to keep the derivations as simple as possible, the constant parameter line model is used as a basis for the new development [41]. It is believed, however, that the method can be applied to the more accurate frequency-dependent line models as well. 6.1 Der ivat ion of the M o d e l Equations Before describing further modeling aspects, let us briefly derive the line model equations and describe their transformation into the modal domain, which permits an easier understanding and modeling of the polyphase case. The equations of a lossless polyphase transmission line can be decomposed into a number of single phase line equations, which can then be solved analytically with d'Alembert's method for one dimensional wave propagation. The equations describing wave propagation on lossless transmission lines are given by d2 d2 v(t,x)-L'C ^v(t,x) = 0 dx2 v ' ' dt2 (6.1) d2 d2 — i(t,x)-C'L' — i(t,x) = 0 The transmission line is characterized by its inductance matrix L' and capacitance matrix C per unit length, where the size of both matrices is determined by the number of phases. For a 6.1. Derivation of the Model Equations 83 solution of the above differential equations, they must first be diagonalized. Using eigenvalue decomposition, an orthogonal transformation T~l = can be found such that A = T-lL'C'Tv = T7xC'L'Ti is a diagonal matrix, and v = T~l v and i = T71 i are the modal values of line voltages and currents [87]. Hence, the n-th order wave equation (6.1) can be decomposed into n scalar equations, where each of them has the form d2 d2 — v(t,x)-\ — v(t,x) 0 d2 - d2 -(6.2) The French mathematician Jean Le Rond d'Alembert (1717-1783) first described a solu-tion [7] to the equation d2 d2 — u{t,x) - c 2 — u(t,x) = 0, where c = +y/X is the propagation constant, at time * i as tl+CX . . lf(ti + cx) + ip(ti - cx) 1 . , . . , u{tux) = ^ 2 2c~ I ^ d a -(6.3) t\ —cx The functions ip and ip are given through the boundary conditions d u(t, x) (p(t) = u(t,x) m = 2 = 0 dx x=0 Let us next apply this method to the first equation in (6.2). For voltages and currents at node (K) of the line, i.e. at x = 0, v(t,0) = vk(t) and i(t,0) = ik{t), while the quantities at the other end © at x = £ are v(t,£) = vm(t) and i(t,£) = -1m(t). Then, <p(t) m vk{t) -L s ! t ( 0 where L is the mode inductance found from the diagonal matrix [86] T~lL'Ti. Similarly, T7lC'Tv is also diagonal with its elements C being the mode capacitance. Substituting the above functions into (6.3) yields _ t l + C X _,. \ ^k(ir + cx) +vk(ti -cx) L f -. . , w(£i,x) = — '— - 2c I ^6 ) tl—cx 6.2. Lossless Line Model 84 Letting x — I, we obtain the voltage at terminal @ at time t\ vm(ti) = Vkih+r) +vk(tl-r) L 2v /A « k ( * i + r ) - i k ( * i - T ) (6.5) where r = £ A / A represents the wave travel time for the mode. Repeating this process for the modal currents in the second equation in (6.2) results in *m(*l) = - g + c 2x/A Vk(h + T) - v k ( t i - r ) (6.6) Substituting the factors with the modal characteristic wave conductance the terms Vk{ti+T) and ik(ti + r) are eliminated to yield the two equations Gc 't)m(ti) - i m ( t i ) = G c w k ( t ! - r ) + i k ( i i - T ) G c 'Wk(*i) - *k ( * i ) = G c w m ( t i - r ) + i m ( t i - r ) from which (6.8) and (6.9) are later constructed with vector notation. Notice that these equations do not contain approximations, but are an exact analytical solution to the trans-misision line wave equations. (6.7) 6.2 Lossless L ine M o d e l The modal decomposition presented in the previous section enables us to reduce a polyphase transmission line into several single phase lines, which can then be separately modeled. In this section, we will hence only consider a single phase line, i.e. one mode of a polyphase line. A l l voltages and currents in the modal domain are indicated with superbar notation. In the modal domain, each mode is characterized by its wave travel time r and characteristic impedance Rc (purely resistive for lossless lines). For easier notation we will also make use of its reciprocal value, the characteristic conductance Gc. More details on modal transfor-mation, transmission line equations and line parameters can be found in the E M T P Theory Book [24]. Before introducing losses, let us now discuss the ideal case of lossless transmission lines. Albeit not very accurate, this case is particularly interesting since the line equations for a 6.2. Lossless Line Model 85 distributed line inductance and capacitance can be solved analytically. Early in this century, d'Alembert's equations were applied to wave propagation on transmission lines. The resulting modeling technique is commonly known as Bergeron's method [5]. In the following, we will describe a model of a single phase (mode), lossless transmission line for electromagnetic transients studies. It is shown that, for this model, the simulation step size cannot exceed the line's travel time. Therefore, two extensions are suggested that overcome this constraint and permit the line to be simulated with larger step sizes, provided the transient state of the network tolerates such an increase. 6.2.1 E M T P Model for Lossless Lines The lossless model implemented in the E M T P is based on Bergeron's method. Figure 6.1 shows the equivalent circuit from which the model's most important property is immediately ® L(t) o-QI ijt) ® a ik(t) Ut) 7 7 ) \ v (m) Figure 6.1: Lossless line model using Bergeron's method for At < r„ clear: Both terminals ® and © are galvanically separated. The currents into the two 6.2. Lossless Line Model 86 terminals (K) and © are given by Mt) im(t) Gr Gr vk{t) vm(t) + hk(t) hm(t) (6.8) " hk(t) ' Gc' vk(t-r) 1 " Mt~r) _ hm(t) _ _ vm{t-r) _ 1 . im{t-r) _ where the history vector is known and can be calculated from previous current and voltage values as follows: i r n i r - u _ \ i r i i r ~ u _\ i (6.9) This last equation shows that the conditions at one end depend on what happened at the other end at one travel time r earlier. As long as the simulation time step is less than r, those conditions at the far end at time t — r will appear at the near end at instant t. This leads to the galvanic separation which is exploited in the E M T P , as shown in Fig. 6.1. To maintain the galvanic separation, it is necessary that A t < T„ and mm for i = 1... n (6.10) where n is the number of modes, i.e. number of phases. This condition does not impose a heavy constraint for simulating electromagnetic transients because the time steps are usually much smaller than rmin, yet it precludes this model from being useful for stability studies. Considering the storage requirements for this model, the values of hk(t) and hm(t) have to be calculated and stored in a history memory at each time step, until these values are used r times later to make up the history vector. The history memory for each mode requires d entries, where d = ceil T At and ceil is a rounding function to give the nearest integer value that is greater than or equal to the argument. 6.2.2 Interpolated Model for At > r The integration method used in the E M T P is the trapezoidal rule of integration. With this, the program assumes that a variable changes linearly as a straight line between two adjacent points. Therefore, if two solutions s(t) and s(t — At) were known, they could be linearly 6.2. Lossless Line Model 87 interpolated to yield s(t —r) provided that A t > r. Here, the variable s stands for any of the line voltages or currents. The concept of this linearly interpolated line model was based on hand-written notes by W. S. Meyer. In a combined work with A . I. Ibrahim, the model equations were derived and implemented for the UBC-version of EMTP, MicroTran®. Let us now take the linearly interpolated model one step further, by applying the linear interpolation formula to the dynamic phasors of the currents and voltages. This should allow us to accurately simulate the transmission line with large time steps during low frequency oscillations. In Section 5.3.1 we had already derived an expression for the extrapolation of dynamic phasors. Similarly, we obtain for the interpolation At - T $>[S(t) 9[S{t- At)) , A t " L - W J A F Substitution with (2.20) yields the interpolation formula for complex signals (6.11) Si(t - r) A t - r S(t) T ,ju>n(At-r) S(t - At) , (6.12) A t " ~ w ' A t where ujn is the nominal frequency ofthe transmission line. Thus, introducing the two factors A t - i V -3 unr ju!n{At-r) (6.13) A t ' •* A t the past values in (6.9) can be replaced by application of this formula to give the new model for the lossless line. As before, the terms are grouped to give the current injections Gc . im(t) _ 1-p2 1+p2 -2p -2p 1+p2 vk(t) vm(t) + hm(t) (6.14) and a history vector whose values are determined from the past time step: hk{t) hm(t) qGc 1 - p 2 V - i - 1 p vk{t-At) vm(t-At) + 1-p2 p -1 -1 p 4 ( t - A t ) v ( t - A t ) (6.15) The equivalent circuit of this new model is presented in Fig. 6.2. We see that the conductance matrix in (6.14) is now complex and also possesses off-diagonal elements, so that the line equations can no longer be solved independently for both ends. Notice that the line's history vector requires only values of the previous time step, so that the history memory for this model is extremely short. 6.3. Lossy Line Model 88 R c2p rmv (m) HK(V Figure 6.2: Lossless transmission line model for A t > r„ With this new model, the simulation time step can be increased to large and very large simulation step sizes, according to the dynamics in the system. However, it should be em-phasized that the interpolation formula (6.12) fails for A t < r. Should the chosen simulation time step fall below this margin, the program must automatically replace the above equations with equations (6.8) and (6.9) by setting p = 0 and q — 1. Different representations of transmission line models for fast and for slow transients has so far prevented a straightforward symbiosis of electromagnetic transients studies and stabil-ity analysis. While the nominal 7r-circuit dominates in stability studies, it does not properly reflect frequency variations, and is also too inaccurate for electromagnetic transient investiga-tions. Constant or frequency dependent parameter line models do provide very good results in electromagnetic transients simulations, but their maximum step size constraint (6.10) has precluded their usage in stability programs until now. The above extension to the EMTP constant parameter model therefore represents an important development for the simulation of transmission lines over a wide frequency range. 6.3 Lossy Line Model In the previous section it was shown how a transmission line can be modeled when losses are negligible. In this section we will demonstrate how transmission losses can be included in the line model. Due to the limited applicability of the linearly interpolated model, this and 6.3. Lossy Line Model 89 subsequent sections will focus on the model equations given by (6.8), (6.9), (6.14) and (6.15). Because of its simplicity, the constant parameter model with lumped resistances [23] is chosen for the loss representation, but the method can likely be extended to the more sophisticated frequency dependent line models. This model is depicted in Fig. 6.3. It consists of two lossless H R / 4 H Ml R/2 Idl Uty R/4 >Jt\ ® ® ® ® ® ® Figure 6.3: Schematic of a lossy transmission line model line sections, one between nodes Q and (2) and the other between (3) and (§). Depending on the simulation time step At and on the travel time for each mode r, either (6.8) and (6.9) or (6.14) and (6.15) are used. Notice that for this purpose the indices in the previous equations for nodes (2) and @ must be replaced with Q - @ in the appropriate places, and the wave travel times now are _ 1 7* 2 ^ i°ssy ' as each of the two lossless sections now has only half of the total line length. Transmission losses are accounted for by including lumped resistances between the two lossless sections and between each lossless section and the respective line terminal. In the more general case of polyphase transmission lines, the resistance per unit length is given as a matrix R' whose size is determined by the number of phases. For easier notation we introduce the conductance matrix to represent the lumped resistance in the following derivation. The goal is to treat the four nodes Q , (2), (3) and @ as internal nodes of the transmission line model, while (£) and © are its external nodes that interface the line model with the rest of the network. Hence, we seek expressions for the terminal currents ik and im which depend only on the terminal voltages i>k and vm, and on a history vector that can be calculated from past values. From Fig. 6.3 it follows immediately that i k (t) = i^t) = T j i i ( t ) . . im(t) = i4(t) = TiU(t) 6.3. Loss}' Line Model 90 The conductance or admittance matrices of the lossless line segments in (6.8) and (6.14) depend on the simulation step size and the wave travel time of the respective mode. Both matrices can be characterized by a self admittance YS for the two diagonal elements and a mutual admittance YM for the two off-diagonal elements for one mode. For n modes, we obtain two diagonal matrices Ys and Ym. The transformation back to phase quantities yields the matrix Y Y Y Y Ti Ys T " 1 TiYnT-1 T l Y m T ~ l Ti Ys T~l (6.17) for the representation of the two lossless sections. Thus, we can write the complete nodal equation for the lossy transmission line model utilizing (6.17), with the lumped resistances connected at the correct places: k(t) -hx(t) -h2(t) -hA(t) where the submatrices are given as A B BT D ' vk(t) ' vm(t) V1(t) Mt) v4(t) (6.18) " G O ' - G 0 0 0 0 G , ti — 0 0 0 - G G + Ys Y 0 0 Y \G + YS 0 0 -W l2G + Ys Y 0 0 Y G + Ys D In order to simplify the notation, let the terminal voltages and currents be vm(t) vt(t) = and the internal voltages and history values ; it{t) k{t) im(t) Vi(t) = [vl(t),v2(t),v3(t),v4(t)}T hi(t) = [Mi),Mt),Mi),Mt)] T 6.3. Lossy Line Model 91 Then, the current injections from the line model into the rest of the network through termi-nals ® and @ can be formulated as it{t) = [A - BDlBT] vt(t) - BD^hiit) (6.19) v v ' ' v ' Y,ine(At) hltne(t) For a simulation with the lossy transmission line model, its matrix Ynne(At) is inserted into the nodal admittance matrix of the entire system, and the history values hnne(t) are placed into the full system's history vector at the appropriate places for line terminals @ and @, as explained in Chapter 5. After the system equations are solved, all bus voltages including vt(t) are known. With known terminal voltages, the terminal currents can be calculated with (6.19). In order to update the internal history vector, (6.9) or (6.15) can be used depending on the selected time step and the wave travel time for each mode. For this, the internal modal voltages as well as the modal currents through the lossless line sections are needed. The internal phase voltages can be obtained using Vi{t) = —D~1BTvt(t) - D^hiit), (6.20) which was derived from (6.18). Modal transformation yields the corresponding voltages vl in the modal domain. The phase currents can be calculated by means of relationship (6.16) and the following equation for the remaining currents: 1 G -G ' v2(t) " . i 3(t) . ~ 2 -G G . v3{t) _ (6.21) Again, the modal equivalents ii are obtained through modal transformation. After the internal history values are updated and transformed back to phase quantities, the model is prepared for the next time step. Chapter 7 Synchronous Machine Since early this century, the modeling of synchronous machines has been a main research topic in power system engineering so that, today, there exist a large number of different models. The choice of the "right" model therefore greatly depends on the anticipated fre-quency range for the study at hand. For the purpose of this work, we will restrict ourselves to the model whose d and g-axis rotor structures are shown in Fig. 7.1, and which seems to be widely accepted in electromagnetic transients studies as well as in stability analy-sis [2, 16, 19, 22, 24, 49, 50, 55, 77]. Figure 7.1: Synchronous generator d and q-axis models 92 7.1. Model Equations with Real Signals 93 While the issue of correct modeling has already been addressed in many of these pub-lications, this chapter is more concerned with adapting the model to the complex variable representation used for the solution process. This task is challenging from two points of view: 1. A generator "transforms" mechanical rotation in space into electrical oscillations in time. Since the motion of complex signals corresponds to the electrical oscillation in time, we must find out if and how it can also represent physical rotations. 2. A generator also "transforms" DC quantities such as the field voltage and current to A C quantities. The concept of complex signals had been developed with respect to A C systems. The question whether this concept is also meaningful for D C signals has not yet been addressed. In the definition of complex signals on page 25, the real part was defined as the real or instantaneous value of a physical quantity. Therefore, a description of the synchronous machine with complex signals must, in its real part, be identical to existing machine model equations. Hence, the main difficulty lies in finding a physically meaningful imaginary part to these equations. 7.1 Model Equations with Real Signals In this section, we will present the electrical and mechanical machine equations for use with real signals. This will enable us to expand on these relationships later, using complex signals. In the following, we assume that the characteristic machine parameters can be obtained through measurements and are known. A new approach for identifying these parameters from standstill frequency response measurements is presented in Chapter 8. With a few modifications, the notation and formulae for describing the machine model are adopted from P. Kundur [55]. The following assumptions are made: • The direct (d) axis is centered magnetically in the center of the north pole, 7.1. Model Equations with Real Signals 94 • The quadrature (q) axis leads the direct axis by 90 electrical degrees in the turning-direction of the rotor. • The armature currents are assumed to be positive (real parts) when they flow out of the machine. For the derivations in this thesis, we assume the field current and electromagnetic torque to be positive when they flow into the machine. • Nonlinear effects generated by low-flux iron saturation are neglected for this analysis in order to take advantage of linear circuit analysis and superposition. Suggestions on how to account for these effects are discussed in [24, 46, 55, 67]. While, in practice, most machine models are limited to two windings on the rotor in each axis, we will assume Nd = v + 1 d-axis windings (v amortisseur windings and one field winding) and Nq = p g-axis windings on the rotor structure. This way, the models can easily be extended from second order to higher orders should this become necessary. The synchronous machine model equations are typically divided into a set of flux link-age equations and another set of voltage equations. The voltage equations for the real or instantaneous values of the three phases and the rotor circuits are (7.1) (7.2) (7.3) (7.4) for K=1...U (7.5) for K = 1... fi. (7.6) The flux linkage for each generator winding is composed of a linkage due to the stator windings and a linkage due to the windings on the rotor structure. Hence, the instantaneous ea{t) = jti)a{t) - Raia(t) eb(t) = ^Mt) ~ Raib{t) ec{t) = jtMt) - Raic{t) efd{t) = ^fd{t) + Rfdifd{t) 0 = ^ ^ ( t ) + Rnd^Kdit) 0 = ^TpKa{t) + RKqiKq{t) 7.1. Model Equations with Real Signals 95 flux linkage e.g. in the phase a winding can be divided into 1pa{t) = 1pa,stator(t) + ^a,rotor(t) (7.7) where 1pa,stator{t) = -laaia{t) ~ Lb ib(t) ~ LcUit) v H Ipa.,rotor (t) = lfd lfd{t) + YJ LKCI iKd{t) + YJ I A K C I kqit) K=l K = l The flux linkages of the other windings are analogous. Notice that the self and mutual inductances are not constant as they depend on the electric rotor position 8(t) with respect to the winding of phase a. Substituting the expressions for rotor position dependence, the stator flux linkages for the three armature windings become 4^a,slator it) LaaO LabO LabO ' ia(t) ' ^b,stator(t) — LabO LaaO LabO ib(t) + 'Pcslatori't) LabO LabO Laa0 ic{t) Laa2 x cos [20] cos [29 + TT/3] cos [28 - vr/3] cos [20 +TT/3] cos [2(0 - 2vr/3)] cos [29 - T T ] cos [29 - TT/3] cos [29 - TT] cos [2(0 + 2TT/3)] " ia(t) ' ib(t) (7-* and the corresponding linkages with the rotor windings are Lafd ifd(t) + y~] LaKd iKd(t) Lpa,rol,or it) "Pb,rotor it) ^c,rotor(t) K=l X cos [9] cos [9 - 2TT/3] C O S [0 + 2TT/3] ] LaKq kqit) j x K = l sin [0] sin [0 - 2TT/3] sin [0 + 2TT/3] (7.9) For the windings on the rotor structure, the flux linkage with the armature windings is determined by / 2ir 2ir \ fpfdelator (t) = - L a f d ( ia{t) C O S [9] + ib(t) C O S [0 - — ] + ic(t) C O S [0 + — ] J (7.10) / 2ir 2ir \ 1pKd,rtator(t) = ~LaKd I ia(t) C O S [ 0 ] + lb(t) C O S [0 - — ] + lc(t) C O S [0 + y ] j (7.11) / 2ir 2TT \ 'ip*q,stator it) = +LaKq lia(t) sin [0] + ib(t) sin [6-—] + ic(t) sin [0 + —], J (7.12) 7.2. Transition to Complex Variables 9 6 where the index is re = 1... v for the d-axis amortisseur windings and re = 1... / i for the g-axis amortisseur windings. For the rotor flux linkages, it is assumed that there is no electromagnetic coupling between windings on the d-axis and those on the g-axis. This is reflected in the following two matrix equations for the d and g-axis flux linkages: 4>fd,rotor{t) Lffd Lfid • 1pld,rotor{t) Lfu Liu • i\d{t) '4)i'd,rotor (t) Lfud L\vd • _ iud{t) _ lplq,rotor (t) L\lq L\2q • • • LiM faq,rotor {t) — Li22q • ' • L 2 M iftnq, rotor (t) _ LiM L2M • • (7.13) (7.14) The instantaneous electrical behavior of a synchronous machine can thus be described by equations (7.1)-(7.14). They have been presented in the literature, and are the starting-point for simplifications made in stability analysis to reduce the number of model equations, and also to interface them with the network equations [55]. The equation for the electromagnetic machine torque is given by Tm(t) = (Mt) *,(*) - Mt) »-(*)) , (7-15) where ipd(t) and ipq(t) as well as id(t) and iq(t) are obtained applying Park's transforma-tion (2.5) to the armature flux linkage and currents, respectively. The electromagnetic torque can be converted to its electrical equivalent, Te(t), through division by the number of pole pairs np/2. 7.2 Transi t ion to C o m p l e x Variables Let us now derive the machine equations for use with complex signals. Notice that, by definition, the real part must always be identical to the equations presented in the previous section (see definition of complex signals on page 25). In order to gain insight into how complex signals can be used to describe physical rotation, let us recall the use of complex 7.2. Transition to Complex Variables 97 signals to describe a steady-state signal: In previous chapters, we wrote S(t) = S(to0) eJwo{i-to) (7.16) to describe a uniform sinusoidal oscillation in time. The starting time of a simulation is usually set to io = 0 so that this term can be dropped in the following analysis. Originally, S{LOQ) was taken from the steady-state results in the frequency-domain. However, in the time-domain, this quantity essentially represents a complex constant. Similarly, we are now going to introduce the signal S(t) = SN eje = SneW°+UTt) (7.17) to represent the time-function of a constant quantity SN, undergoing a uniform rotation in space with constant speed ojr. Here, SN could be a steady-state D C quantity such as the current in the field winding of a machine. In general, the quantity is complex so that it can be decomposed into a real and imaginary part. For quantities associated with the machine rotor structure, it is common to identify orthogonal d and g-axis components. Hence S(t) = sdeie + j s g e j e . (7.18) The above equation still describes a uniform sinusoidal oscillation. More generally, sd and sq can both be independent functions of time, thus S(t) = sd(t)e^0 + jsq(t)e^ (7.19) Notice that Sn(t) = sd(t) + j sg(t) is the time-function of a complex rotor quantity with a contribution in the d-axis as well as in the g-axis. However, in practice, the rotor windings are hypothetically assumed to be either on the d-axis or on the g-axis, therefore having no mutual flux linkage at all. Depending on the considered axis, either sd(t) or sq(t) is always zero. For "complex" signals representing d-axis quantities, it follows that Snd(t) = sd(t), while "complex" signals associated with g-axis quantities are Snq(t) = sq(t). With this definition, complex signals associated with a rotor quantity are identical to the instantaneous values of this quantity, and the imaginary part is always zero. Hence, we can continue to use real values for rotor structure quantities. 7.2. Transition to Complex Variables 98 7.2.1 Rotor Structure Quantities Observing this identity, we can now apply (7.17) to express the stator flux linkage with the physically revolving rotor windings in complex notation: ^ a,rotor (t) ^b,rotor (t) ^c,rotor (t) ,j(9-2*/3) ,j(9+2ir/3) + LaKa iKq(t) j X K=r 3 j 0 eJ(0-2*/3) pj{0+2-n/3) (7.20) A comparison shows that the real parts of this equation are equal to the instantaneous flux linkages in (7.9). Applying Park's transformation (2.5), we obtain the complex dqO coordinates of the above flux linkage: ^'d,rotor (t) ^ q,rotor (t) ^0,rvtor(t) A-'afd ifd{t) + J2L™di*d^ x K=r -J 0 + ^ ] LgKg JKq(t) J X K=1 / 1 0 (7.21) Because of the identity of complex and instantaneous rotor structure quantities, the equations describing the flux linkage among the rotor windings, (7.13) and (7.14), remain unchanged. Similarly, the rotor flux linkage with the armature windings also remains the same. However, we do assume complex signals for the armature currents, and must therefore extract the real part before substituting them in (7.10)-(7.12). For shorter notation, we temporarily omit the time-dependence of the currents. Tpf delator (t) = -Lafd (K[/a] cos [9] + %[Ib] COS [8 - ~] + U[IC] COS [8 + — \ o 3 A.d,Stator(t) = ~LaKd (Wa}cOs[9} + ^ Ib}COs[9-~} + U[Ic]cOs[9+~-\ 3 3 iPwtatorit) = +LaKq U[Ia] sin [8] + 5R[/6] sin [9 - ~] + &[/ c] sin [9 + — V o 3 (7.22) (7.23) (7.24) With these equations we have achieved a translation of D C signals on the rotor structure to A C signals in the armature windings, using the complex variable representation. As a 7.2. Transition to Complex Variables 99 result of this derivation, it was established that a D C signal has no imaginary part and is therefore identical to its corresponding real or instantaneous signal. This observation is consistent with our previous understanding of complex signals: Let us consider (7.16) and set LO0 — 0. Consequently, the complex DC signal becomes SDC(t) = S(co0 = 0) , (7.25) which is the D C steady-state solution. At zero frequency, the steady-state network equations in the frequency-domain are purely resistive and introduce no phase shift. Since DC source functions are real functions, it follows that the steady-state solution and thereby the DC signal are real functions as well. 7.2.2 Three-Phase Armature Quantities Let us now consider the implications of the complex variable representation for the flux linkages among the three armature windings. In Chapter 2 we had demonstrated that the nodal equation for an RL-series branch is the same for a representation with real as well as with complex signals (see Table 2.1). This is acceptable if the branch impedance is constant. For a rotating machine, this is not the case, because the inductance seen by the arma-ture currents changes with the rotor position. To illustrate this fact, let us apply Park's transformation (2.5) to the armature flux linkage equation (7.8). This gives (7.26) ,lPd,stator{t) ' Ld ' id(t) ' '4>q,stator(t) = — '$0,stator(t) L0 io(t) where Ld Lq L0 + Laoo + - Laa2 LaaQ + LabO ~ ^ Laa2 = L aaO 2L abO It clearly shows that the d-axis and g-axis components of the armature flux linkage are influenced by different inductive paths. Notice that both components are orthogonal. 7.3. Model Equations with Complex Signals 100 Introducing complex signals for the dqO components, we have *o(*) = V>o(<)+j '3[*o(*)] (7.27) (7.28) (7.29) Because of the orthogonality of d and g-axis components as well as of the above real and imaginary parts, ipd(t) and are both alined along the d-axis, whereas ipq(t) and j $s[$d(t)] are located on the g-axis. Hence, ipd(t) and jS[\I>9(£)] share the same magnetic path through Ld, while ipq(t) and j 3 [ f ^i)] both flow through Lq. The inductive path for the zero sequence remains unchanged. With this observation, we can now easily express the complex armature flux linkage in dqO coordinates: ^d,stator{t) ' Ld ' Wd{t)\' ^q,stator (^) = — Lq -j ^0,stator{t) Lo Wo(t)] Ln La Wait)} 9[/o(t)] (7.30) This equation completes the transition from the real flux linkage equations to their com-plex equivalent. We have thus created a physically meaningful, imaginary counterpart to the electrical machine equations. 7.3 M o d e l Equations w i t h C o m p l e x Signals In the previous paragraphs, we derived a complex representation of the machine equations. In this section, these equations are properly arranged in order to provide a consistent model that can be interfaced with the network equations. Keeping in mind that the variables are all functions of time, we can omit the time-dependence in the notation for the following analysis. 7.3. Model Equations with Complex Signals 101 7.3.1 Flux Linkage Equations Considering the four equations for the flux linkage, (7.13) (7.14), (7.21) and (7.30), we can write three real matrix equations: Ld Lafd Laid La^d ' Wd\ -Ld —Lafd — Laid ' —Laud 3 2 Lafd Lffd Lfid Lfvd kd 3 2 Laid Lfid L\\d L\ud kd i>vd -3 2 Laud Lfvd Llvd Lund kd ~Lq Lalq La2q " " Laftq ' Wo) ' ~Lq Lalq La2q " " La/iq 3 T 2 ^al1 Lllq Ll2q ' • " LiM kq 3 T ~2 Ll2q L22q - • • L 2 m kq -1 T 2 ^ajiq Ll/iq L2liq • • Ljxiiq (7.31) (7.32) and ' -L0 ' K[ / 0 . 3 [ * o ] . 3[ / 0 (7.33) Notice that the d and g-axis equivalent circuits shown in Fig. 7.1 do not correctly repre-sent the flux linkage relationships obtained with these equations. This would require more peripheral leakage inductances, similar to Lfkd-, between the other rotor RL-branches. The purpose of these inductances is to properly scale the currents in the rotor branches. Since the exact current-splits in the amortisseur windings are of no concern for most transient studies, these inductances can be generally omitted. It is important however to retain Ljkd in the of-axis model to ensure an accurate representation of the field winding [8, 18]. With the equivalent circuit shown, we must now scale the equations associated with flux linkage in the damper windings so that the peripheral leakage terms disappear. Substituting Lana • , , 2 La£)q iDqK = j and ipDqK = — ipKq L/aDq <-> ^aKq 7.3. Model Equations with Complex Signals 102 in (7.32) for the variables of the Ac-th g-axis amortisseur winding, results in the following expression: LaDq LaDq LaDq ~Lq LaDq LaDq LaDq '<pDql Lo\\q LD12q • • Loijiq i-Dql i>Dq2 — LaDq Lmiq Lm2q • • LD2HQ i-Dq2 i>Dqv. — LaDq Lo\jLq LD2p,q " LDnnq i-Dq/j. (7.34) where LDxyq sponding to 2 Lxyq LaDq) / (3 Laxq Layq) • The resistance and inductance values corre-L DqK L DKKq La Dq and RDqK = 2L aDq o 7-2 K<? 6 nanq for K — 1 . . . fj, are depicted in the g-axis circuit diagram in Fig. 7.1. Their values can be ob-tained from field test measurements as explained in Chapter 8. The magnitizing inductance can be calculated from LaDq — Lq — where the leakage inductance Li is available from manufacturer's data. Similarly, we now proceed with the scaling of the d-axis equations. The difference there is that the field quantities must not be scaled if they ought to maintain their correct values at the exciter terminals. Substituting ^D(1K — r ^kq OK in (7.31), where results in the following matrix equation: and ipDdK = - 5K 'ijjKd X — ^ LfDd Lafd °K — \ l TTT T ) Z J^fKd, J^aKd n^d] -Ld Lafd Lafd • • Lafd ' Wd\' — Ld —Lafd — Lafd - - • -Lafd I'^Jd -Lafd 2 T LfDd • • • LfDd -Lafd LfDd L\\d LfDd i-Ddi '4>Ddv _ -Lafd LfDd LfDd • • Lpvd ^Ddv (7.35) 7.3. Model Equations with Complex Signals 103 This scaling always works for d-axis models up to the third order. For higher orders, we must assume that the mutual inductances among the damper windings, i.e. Lxyci in (7.31) are all equal. Since higher-order representations of these windings are hypothetical, we can make this assumption without loss of generality. The resistance and inductance values in Fig. 7.1 can be related to the above quantities according to the following equations: 2 2 LDOIK, = LDKKCI — LjDd and RDCLK — ^ <5K RKd where K = 1... v and where Lafd = Ld — L( is the magnitizing inductance. The peripheral leakage inductance Ljkd — LfDd—Lafd is related to Canay's characteristic inductance through the relationship Lc = Ld — LfDd. Due to the lack of proper measurement methods, this inductance is usually set to Lc = Lu in traditional machine representations. However, newer measurement techniques, such as standstill frequency response test, provide the means to unambigously define this inductance (see Chapter 8). 7.3.2 Voltage Equations While the voltage equations of the rotor windings can still be expressed using real signals, the equations of the three machine armature windings must be reformulated for complex signals. For convenience, the three-phase quantities are transformed to dqO coordinates using Park's transformation. Notice that we must now use the scaled resistance values in order to be consistent with the above flux linkages. Thus, Ed = IVi-Ralt-Ur*, (7-36) at E q = ± y n - R a I q + u;r*d (7.37) at efd = ^ipfd + Rfd'^/d (7-39) 0 = ^TlpDdK + RDdnlDdK forK-=l . . . i> (7.40) dt 0 = ^1pDgK + RDq^Dq* for K = 1 . . . U . (7.41) 7.4. Time-Domain Simulation 104 7.3.3 Mechanical Equations Similar to the electrical rotor quantities, the electromagnetic torque also represents a D C signal. Consequently, it has no imaginary component and equals the real definition in (7.15). However, let us explicitly denote the real parts of the d and g-axis currents and fluxes. Hence, 3 77 Tm = - ^ ( B [ * d ] H [ J , ] - » [ * , ] » [ J r f ] ) • (7.42) The deviations in the mechanical rotor position are given by 1 A0 m = — (ur - UJ0 ) (7.43) dt np 2 = oom wo , (7.44) where ojm is the mechanical rotor speed and where np/2 is the number of pole pairs. The moment of inertia of a turbine-generator is not directly considered in the model equations. Rather, it is part of a rotor shaft model, which can be represented with linear RLC-branches for most types of simulation. This technique was also realized with the uni-versal machine model in the EMTP [24]. It allows the translation of the mechanical equations into their electrical equivalents. For instance, for the torque: Tm(t) = Jmjtujm{t) i(t) = Cjtv{t) (7.45) where Jrn = np Ekinetic/(2u)0) is the moment of inertia. The inertia constant H is defined as the kinetic enerergy Ekinetic divided by the power rating of the machine. While most power system analysis tools take the machine's apparent power rating Srating to build this ratio, a normalization with the real power rating Prating was chosen in [22]. 7.4 T i m e - D o m a i n S imulat ion In Chapter 5 we presented the solution process of a transient simulation. This section describes how the synchronous generator model can be integrated into this solution process. 7.5. Steady-State Initialization 105 Since a generator model is nonlinear, the dependent current source subtype will be used for it, which was introduced in Section 5.3.2. The source functions are the three armature currents Ia{t), hit) and Ic(t). the field current ifd(t), and the electromagnetic torque Tm(t). As we can see, these quantities are functions of the flux linkages and the rotor position 0(t) = 90 + A9{t) + u0t . (7.46) The currents can be obtained in dqO coordinates through (7.34), (7.35), and (7.33). With this, we have the field current ifd(t) as well as 5R[/rf] and 5R[/g], needed to calculate the electromagnetic torque according to (7.42). Applying Park's inverse transformation (2.6) gives the complex armature currents Ia(t), I0(t) and Ic(t). Besides these injections, the generator model also has a set of state-space equations, which consists ofthe voltage equations, (7.36)-(7.41), in addition to the rotor angle equation (7.44). These state-space equations are used to update the derivatives at each time step. Notice, that um in this last equation is interpreted as the voltage across the mechanical mass terminal, where the electromagnetic torque is treated like a current injection into this terminal. Both quantities are therefore part of the voltage and current vectors, respectively. In addition to updating the current source functions and state variable derivatives, the algorithm for the generator model must also provide a function to compute its Jacobian matrix. 7.5 Steady-State Ini t ia l izat ion During steady state, the dqO transformation of balanced steady-state armature quantities produces complex, constant dqO quantities. The quantities associated with the rotor wind-ings are real and constant. Therefore all time derivative terms disappear. Under balanced conditions, the zero sequence quantities are all zero and can be dropped from this analysis. 7.5. Steady-State Initialization 106 With the electrical rotor speed the voltage equations (7.36)-(7.41) become Ed = -Rald-OJo^a (7.47) Eq = -Ralq+UoVd (7.48) efd = Rfdifd (7.49) 0 = RDCLK iDdn for K = 1...U (7.50) 0 = RDqniDqK for K = \ . . . fl . (7-51) From these equations it follows immediately that the amortisseur currents of the d and g-axis damper windings are zero. The flux linkage equations yield jipfd = -Lafdft[Id] + lLffdifd (7.52) 4>DdK = -Lafd$l[Id} + LfDdifd ior K - 1 . . . U (7.53) •tpDqn = -LaDqU[Iq] f o r « = l . . . / z (7.54) and for the stator flux linkages: Vd = -Ld$t[Id}-jLq%[Id] + Lafdifd (7.55) Vq = -LqM[Iq]-jLd%{Iq}-jLafdifd (7.56) Because of \&d = j tyq, both equations are linearly dependent. In order to initialize the above quantities at steady state, we must first determine the rotor position with respect to one of the phases. Phase a is typically selected as a reference. If we define the internal rotor angle 5l as the angle by which the g-axis leads the phase a reference axis, the following relationships can be established: Eq = En e*a-W ; / , = /„ e*(°+*-*) (7.57) where E(u>o) = En eja and I(u>o) = In e^a+^ are the steady-state phasors obtained with e.g. a power flow solution. Substituting (7.55) in (7.48), we have Eq + RaIq+ jUJ0 Lq 3[Id] = -OJQ Ld U[Id]+ uj0 Lafdifd . (7.58) By expanding both sides by the term LJQ LqU[Id] and observing the relationship Id — j Iq, the equation becomes Eq + (Ra + joj0Lq) Iq = cu0 [Lq- Ld) $l[Id]+ UJ0Lafdifd . (7.59) 7.5. Steady-State Initialization 107 The right-hand side of this equation consists only of real values whereas the left-hand side is, in general, complex. This allows us to determine the internal rotor angle. Let us first introduce the following notation for the right-hand side: Kl = UJ0 ( Lq - Ld ) In sin [ a + </> - Sl ] + co0 Lafd ifd . (7.60) Then, we can determine Ki and the angle Si with Ki ejSi = E(u0) + {Ra+juj0Lq) I(co0) (7.61) Once these two quantities are obtained from a phasor solution, the other armature and rotor quantities of the synchronous machine model can be easily determined using the equa-tions at the beginning of this section. In this chapter, we have demonstrated that the standard synchronous machine model for electromagnetic and electromagnetic transients studies can be modeled using the complex variable representation. It was shown that this representation allows a consistent transition from complex A C signals to real DC signals, and that it is suitable to describe oscillations in time as well as in space. The next chapter presents a new identification method for the electrical synchronous machine parameters from standstill frequency response measurements. Since these measure-ments only involve real signals, the analysis with complex signals and dynamic phasors is concluded at this point. Chapter 8 Synchronous Machine Parameters In this chapter a new method to identify the electrical parameters of a synchronous machine from standstill frequency response tests is presented. In the past two decades, frequency response tests, carried out at both standstill and synchronous speed, have become widely accepted as an alternative to sudden short-circuit tests for the determination of synchronous machine parameters. While short-circuit tests have been well defined and standardized [44], the measurement and calculation' procedures for frequency response tests currently undergo constant improvement and revision. Besides the avoidance of mechanical stress on the machine shafts during testing, a major benefit is that test results allow a more accurate modeling of the electrical machine parts [45]. Yet, there are some difficulties involved with extracting the necessary parameters from the measured responses. As long as the order of the model does not exceed two, the parameters can be obtained analytically [16, 55, 77]. For higher orders, the I E E E Standards [46] suggest the use of iterative techniques. In practical applications it is typical to first extract the time constants by curve fitting, and then to identify the parameters by solving a set of nonlinear equations [19, 43, 82]. A recent paper [85] pointed out the weakness of this approach: The order of the model must be known before the data is analyzed, and an initial estimate of the parameters and/or time constants is needed to lead the iterations to the desired solution. Due to the nonlinear nature of the problem, there is no unique solution that produces the measured responses. With improper starting values, the resulting parameters may nicely fit the measured graphs 108 8.1. Frequency Response Tests 109 but may not necessarily be physically meaningful, e.g. negative time constants, resistances, or inductances. To avoid such a 'blindfold' approach of iterating on perhaps an unfortunate initial guess, a method is thus needed to provide a good starting point that will converge to a physically meaningful solution. This reference [85] also presented the underlying physical relationships between the mea-sured frequency responses and the characteristics of common machine equivalent circuits, which encouraged this work. The technique presented in this chapter identifies the machine parameters without iterations, and also determines the order—i.e. the number of rotor cir-cuits, Nd and Na, on the right-hand side of Node A in the diagrams of Fig. 7.1—directly from the data without prior knowledge [40]. The generator model used in this analysis is the one from the previous chapter. As before, the nonlinear effects due to low-flux iron saturation and thermal variations, while taking the measurements, have been neglected. The presented techniques have been validated with the standstill frequency response data obtained for Ontario Hydro's 555 M V A , 24 kV, 0.9 P F Lambton generator [19, 55]. 8.1 Frequency Response Tests The outcome of standstill frequency response tests are four transfer functions that describe the machine behavior over a certain frequency range. The procedures necessary to mea-sure these functions are well defined in the IEEE Standards [46]. This section is intended to establish the mathematical relationship between the measured signals and some of the parameters involved. The operational inductance is measured with the field winding being shorted so that Ae / d = 0. In the Laplace domain A is used to denote a small change, and s = ju is the Laplace operator. The operational inductance in <i-axis direction is: 8.1. Frequency Response Tests 110 Similarly, we obtain for the operational inductance in g-axis direction: The armature-to-field transfer function establishes the ratio of field current to armature current over the observed frequency range. More precisely, it is defined as ffW - (8.3) s Aid with Q(s) being referred to the stator side. What is actually measured is: G(s) = hG(s), (8.4) where h = §N a/Nyy is the turns ratio between the armature and the field winding. Notice that this ratio is only needed if the per unit values of model parameters are desired (i.e. L a d-base reciprocal system [55, 75]). If the correct unit values of the previous chapter are used then set h = 1. The fourth function is the armature-to-field transfer impedance. Since measured data for this function are usually not available, it will not be used in the parameter identification process. Then, how can the electrical machine parameters be determined from the remaining three transfer functions? To begin with, some parameters can be directly obtained by extrapolating the transfer functions towards s = 0. Thus we have: Ld = l i m £ d ( s ) Lq = l im£ , ( s ) s—>0 G0 = h ^ = nmg<s) Rfd In the absence of a test for the leakage inductance Z^, it is common practice to use the manufacturer's calculated values. Knowledge of Li leads to the mutual inductances: Lafd — Ld — Le LaDq = Lq — Ll-The armature resistance Ra is obtained as part of the process of constructing the operational inductances Cd(s) and Cq(s). The remaining unknown parameters are those associated with the rotor windings, more specifically, the values of Rodn, RDqK, LDdK, and LDqK that constitute the rotor RL-branches in addition to Lfkd. A possible way to identify these parameters is the subject of the following paragraphs. 8.2. Transfer Functions 111 8.2 Transfer Functions Suppose that a transfer function can be represented as a polynomial with a certain number of poles and zeros. Applying partial fraction expansion, these polynomials can be expressed as a sum of first-order derivative terms: M n («+ZK) N Hs) = F0^ = E 7 J ^ T (8-5) K = l Each first-order derivative term can be represented as a series connection of a resistance and inductance: The differential equation of the re-th RL-branch is diK vK = -RKiK - L K ~ ^ (8.6) with the current flowing out of the branch. After the Laplace Transform, we have - A z K l/LK (8.7) AvK (s + RK/LK)' The sum over N such terms thus corresponds to paralleling TV RL-branches to assemble a circuit resembling the right-hand side of Node A of the synchronous machine equivalent circuits shown in Fig. 7.1. Hence we can write for both g-axis and d-axis models: A ^ i = y W ( 8 . 8 ) A - JVd-1 _ J _ J -1X1 Ad _ LDdK Lfd /g g\ AvAd ^ (S+ f ^ ) (S+P-) In these equations, the voltage and current changes correspond to voltages and currents at Node A of the models, always with the currents directed toward the machine terminals. The objective is now to express both admittance polynomials through the known parameters and transfer functions on the left side of this node. 8.2.1 The Q-Axis Admittance In the g-axis model, the values of all inductances to the left of Node A are already known and can be used to express the transfer function (8.8) as a function of the ratio of stator 8.2. Transfer Functions 112 Lambton GS: Magnitude and Phase of Aq(s) CO f in Frequency [Hz] Figure 8.1: Admittance Aq(s) of the Lambton generator and its final approximation by a third-order polynomial. voltage Avq and terminal current ( — Aiq). Finally replacing this ratio by sCq(s), yields the g-axis transfer admittance Aq(s) = -AiAq 1 1 1 (8.10) AvAq s \Cq(s) - Le LaDq which is needed to extract poles and zeros, in order to evaluate the remaining parameters of the g-axis model. Figure 8.1 depicts the g-axis transfer admittance of the Lambton generator. 8.2.2 The D-Axis Admittance If the series inductance Ljkd were omitted in the d-axis model, its characteristic equation could be written in the very same fashion as (8.10) for the g-axis model, since both Nodes A and B (see Fig. 7.1 on page 92) would collapse into a single node. In this case the armature-to-field transfer function Q(s) would not be needed. With the inductance in place, however, the transfer function only satisfies the conditions at Node B so that we can define -AiBd 1 / 1 1 Bd(s) AvBd s \£d(s) - Le La :fd (8.11) 8.2. Transfer Functions 113 where the voltage and current relationship between the two Nodes A and B is given by •%r- = im - si'<" <812> -AiAd Bd(s) Notice that AiAd and AiBd are identical. If L^d had a known value, we could easily write the characteristic equation to extract the poles and zeros from, and find the remaining parameters. Yet the value of this inductance first needs to be found. For its determination, we make use of the armature-to-field transfer function. Using (8.3), (8.4), and (8.11) to substitute the following expression Aid = AiAd - ^ (8.13) SJ^afd we arrive at a transfer function describing the ratio of the current through Lfkd, namely AiAd to the field current Aifd: C(s) = ^ = J - ( L d - C d { s ) ) . (8.14) h Aifd sg(s) V Lafd J Expanding both numerator and denominator by the term AvAd allows us to define yet another transfer function: f ^ _ 1 AvAd _ 1 AiAd AvAd L f d [ S ) ~ sh (-Aifd) - sh Aifd (-AiAd) ( 8 - 1 5 ) The first part of this equation characterizes the voltage-to-current relationship accross the field winding. Applying (8.7) to the field winding results in exchanging the index K by fd. Building the inverse and relating the quantities to the field winding yields 1 AvAd s + rjd/lSd h{-Aifd) l/lfd • where rfd = \ Rfd = 4^ -. h Q This equation together with (8.12) and (8.14) inserted into (8.15) and solved for lfd gives: ' ^ U T T - 1 ' " ^ (817) v v ' b a Since both the unknown parameters lfd and Lfkd ought to be real and constant with respect to s, we can obtain them by dividing real and imaginary parts of above equation: 3(a) 8.2. Transfer Functions 114 Lambton GS: Magnitude and Phase of the 0-th Residue m -a P 'S a3 cS Frequency [Hz] Figure 8.2: Admittance Aq(s) of the Lambton generator approximated by a first-order func-tion of type .7-1 (s). Theoretically, there is a chance that S(6) could become zero causing Lfkd to be undefined; however in practice, this inductance tends to be close to zero such that normally \$s{b)\ is rather large compared to |S(a)|. Practical aspects of the parameter extraction process are discussed in more detail in Section 8.3. Having determined Lfkd, the characteristic equation of the d-axis transfer function becomes: -AiAd Bd(s) Ad(s) = (8.20) AvAd l-sLSkdBd(s) From this equation, poles and zeros are extracted as shown next in order to obtain the remaining d-axis parameters through partial fraction expansion. 8.2.3 Parameter Extraction Both the d-axis and g-axis models represented in Fig. 7.1 always have as many open-circuit time constants as they have short-circuit time constants, even though the number of d-axis model time constants may differ from that of the g-axis model. This observation follows directly from eigenvalue calculations. Consequently, the order of the numerator and denom-inator polynomials constituting the transfer functions Cd(s) and Cq(s) with respect to s are 8.2. Transfer Functions 115 equal. This implies that the order of Ad(s) and Ag(s) is higher by one in the denominator than in the numerator. Dropping the index, we can write for both functions: A(s) = ^ ^Y[p^4, (8.21) s+PN i=i (s+pKy where N is the number of the rotor circuits of the ci-axis or g-axis model, i.e. either or Nq. Hence, there are two different types of fractions appearing in the polynomial: s +p T2(s) = kS-±^ (8.23) s + p The pole p and gain k of the first function can be obtained through the following two relationships: p— c—— and k=p\\m.T\(s) tan^i *->o Principially, the phase angle <pi could be taken anywhere together with the associated fre-quency w i . Due to measurement inaccuracies, small phase angles are difficult to interpret so that we select u\ to be the frequency where dcft/dco has its first negative peak value, and 4>\ is the corresponding phase shift at this point. The approximation of Aq(s) by a first-order function T\(s) is shown in Fig. 8.2. In order to extract its pole p from the transfer function, the latter is divided by T\(s), which takes care ofthe term in front of the product in (8.21). The remainder is called the first residue, from which we can now extract a function T2(s) of the second type as is depicted in Fig. 8.3. The gain k as well as pole p and zero z of a first-order term can be determined from the first occuring maximum phase shift (f)2 and its associated frequency co2: , 1 + sin </>2 1 - sin </>2, k = min ( , j 1 — sin 02 1 + sin <f>2 r~ ^2 p = uj2vk and z =-j= Having found the pole-zero-pair, the new residue is formed by dividing the present one by the most recently determined function T2(s). This extraction procedure is repeated until there are no more noticeable phase shifts present in the residuals. The transfer function can 8.2. Transfer Functions 116 Lambton GS: Magnitude and Phase of the 1st Residue OH 10 10 10 10 10 Frequency [Hz] Figure 8.3: First residue of Aq(s) approximated by a first-order function of type ^ ( s ) . then be written in the form (8.21). The count of poles determines N, the number of rotor circuits of the model. To convert the polynomial form into a sum of first-order terms, partial fraction expansion is used: N-l F] (s + zK) N A{s) = A0 K = l N K=l kK (8.24) where kK = A0 N N-l n (zn - PK) II (P/x - PK) The RL-branch parameters can then be easily determined from the first-order summands as indicated in (8.8) and (8.9). The complete approximation of Aq(s) by a third-order polynomial of the order Nq = 3 was already shown in Fig. 8.1. 8.3. Practical Considerations 117 8.3 P r a c t i c a l Considerations Unfortunately, there are limits to the theory discussed in Section 8.2. Since the transfer functions Cd(s), Cg(s) and sQ(s) are all directly obtained from measurements, there will be some measurement inaccuracies superimposed on the true data. This may cause some calcu-lated functions to misbehave and conceal their poles and zeros in amplified noise. Therefore, it is necessary to consider some practical aspects which are discussed below. 8.3.1 Pre-Processing Filter It is the nature of both functions Ad(s) and Aq(s) that even small inaccuracies in the measured data may cause substancial oscillations of these two functions, primarily on the low-frequency side of the phase angles. Yet the phase angles are the most important tool in extracting the parameters. Hence, the measurement data needs to be pre-processed before it can be analyzed for its poles and zeros. The objective is to smoothen the data, which can be done with a digital filter as simple as the 5s formula [35], i.e. a five-point averaging formula. For a higher level of accuracy and sophistication, one could apply the pole-zero extraction method described earlier directly on the measured data, obtaining a smooth polynomial for each of the transfer functions as in [85]. Since the objective of our algorithm is to provide a good starting point for iterative refinement, we prefer the simpler method for smoothing the data. 8.3.2 Adjustment of Earlier, we described how to determine the inductance Lfkd. Yet, only in theory does the right-hand side of (8.18) give a constant. Due to measurement and/or filtering inaccuracies, this equation describes the inductance as a function of s. Where in the frequency range lies the correct value of Ljkd? Going through the entire frequency range of measured and filtered data points, Lfkd and lfd are calculated following (8.18) and (8.19). For the data of the Lambton generator, the results are shown in Fig. 8.4. With each value of Lfkd we can evaluate the right-hand side of (8.15), while ljd and rfd are used to determine the left-hand side, using (8.16). In theory, both sides should give the same results. In practice, however, .3. Practical Considerations 118 100 10 10 10 10' 1000 -1000 -2000 10" 10" 0 10- 10u Frequency [Hz] 10 10' Figure 8.4: Parameters Ljkd and lfd as functions of frequency. The cross indicates the best fit, yielding Lfkd = -0.0282 p.u. and lfd = 0.831 p.u. there is a mismatch between both. Calling the right-hand side RIGHT(Lfkd,s) = C(s) sBd(s) and the left-hand side LEFT(Lfkd,s) = lfd + LfkdC(s) rfd (8.25) (8.26) an error function can be established as NoF (8.27) , . / RIGHT'(Lfkd, St) , ( J W ~ 2-. ( LEFT{Lfkd, Si) ~ where NoF is the number of discrete frequency points at which the data were measured. The pair of Lfkd and lfd that gives the least error is then selected for further calculations. The cross in Fig. 8.4 marks the occurrence of the optimum parameters which produce the mismatch shown in Fig. 8.5. 8.3.3 Relating Cd(s) and sQ(s) As can be seen from Fig. 8.5, in practical applications ri(Lfkd) tends to be greater than zero even with perfectly selected inductances Lfkd and lfd. This may have two reasons: First, 8.3. Practical Considerations 119 Lambton GS: Magnitude and Phase of Lfd(s) bD CO JS PH Frequency [Hz] Figure 8.5: Best concurrence of functions RIGHT (solid line) and LEFT (dashed line) indicates the optimum value of Lfkd. the separate measurement and independent filtering of Cd(s) and sQ(s) may cause both transfer functions to be slightly uncorrelated, as a result of which RIGHT(Lfkd, s) is not a pure first-order term as demanded by (8.26). Second, the machine might be built in a way that the model of the field winding requires a more complex structure than a first-order RL-branch. If the latter is the reason for the inconsistency, then the model we are trying to fit is not suitable for the task at hand. Before giving up on it, we assume measurement and filtering inaccuracies to cause the mismatch, implying that either Cd(s) or sQ(s) or both are slightly inaccurate to give the desired first-order term. Let us define the mismatch between RIGHT(Lfkd, s) and LEFT(Lfkd, s) as a residual function g(s) such that 1 LEFT(Lfkd, s) RIGHT(Lfkd, s) g(s) aCd(s) - b g{s)sg(s) (8.28) (8.29) with a = 1 L fkd b T - L L f k d T L>1 + Ld. La}d Lafd Unfortunately, there is no way of determining which part of g(s) is due to a mismatch in Cd(s) and which is caused by inaccuracies in sQ(s). Because of the difference of scale—Ld is 8.4. Results 120 much smaller than G0—the operational inductance Cd(s) is far less sensitive to errors than the armature-to-field transfer function sQ(s) so that we assume that the mismatch g(s) is caused entirely by Cd(s). Hence, we can write LEFT(Lfkd, s) = - J — (aCd(s) - b) , (8.30) sy(s) v ' where Cd(s) denotes the corrected operational inductance. Its use ensures a pure first-order RL-branch to represent the field winding. Setting (8.29) and (8.30) equal, yields £ W = ^ W - * ( JL (8.31) g(s) a \g(s) ) Using this inductance instead of Cd(s), we can now calculate the transfer admittance Ad{s) according to (8.20) and obtain the parameters as described earlier. 8.4 Results For the Lambton generator case study used in this chapter, Figures 8.6 to 8.8 depict the three transfer functions that served as input as a direct result from standstill frequency response tests (asterisks) beside the functions calculated from the identified parameters (solid line). Lambton GS: Magnitude and Phase of Ld(s) 10' 3 10 Z 10'' 10° 101 10 2 Frequency [Hz] Figure 8.6: Measured operational d-axis inductance Cd(s) with its second-order analytical (solid line) and iterated approximations (dashed line). 8.4. Results 121 Lambton GS: Magnitude and Phase of sG(s) Oi : : : • • • • •. : : • : : : : 1 : : : :: 10'3 10'2 10'' 10° 101 102 100i OH 10"3 10'2 10'1 10° io 1 i o 2 Frequency [Hz] Figure 8.7: Measured armature-to-field transfer function sQ(s) with its second-order analyt-ical (solid line) and iterated approximations (dashed line). Notice that there has been no additional refinement on the parameters presented here. For comparison, iterative results obtained by Dandeno and Poray [19] are also shown (dashed line). Lambton GS: Magnitude and Phase of Lq(s) 10'3 10'2 10"' 10° 101 102 Frequency [Hz] Figure 8.8: Measured operational g-axis inductance Cq(s) with its third-order analytical (solid line) and iterated approximations (dashed line). 8.4. Results 122 The corresponding parameters are listed in Table 8.1. The right-hand column presents the values obtained iteratively by the SSFR 2.2 and SSFR 3.3 models in [19] for the d-axis and g-axis parameters, respectively. The value of the leakage inductance Lt was specified manually. Table 8.1: Comparison of identified electrical synchronous machine parameters. Parameters Values Ref. [19] p.u. 0.155 0.155 Lafd p.u. 1.7351 1.8580 LaDq p.u. 1.6450 1.7620 No- 2 li 0.2 n/a Rfd p.u. 0.000856 0.00108 R-Ddl p.u. 0.4482 0.01065 Lfkd p.u. -0.0282 0.1328 Lfd p.u. 0.1840 0.0105 Lndl p.u. 0.1924 0.0114 Nq 3 RDql p.u. 0.00505 0.01060 Roq2 p.u. 0.01343 0.0301 Roq3 p.u. 0.10445 0.2610 Luql p.u. 3.8647 2.1010 Ll)q2 p.u. 0.4842 0.4504 L®q3 p.u. 0.1239 0.1131 Notice that Lafd and LaDq in the right column are already adjusted to air-gap values to account for saturation. Although the values listed here bear little resemblance, some interesting observations can be made: By comparison, it was found that the parameters of the d-axis model yield a better fit to the given transfer functions Cd(s) and sQ(s) than the parameters obtained iteratively. This is especially true for the high-frequency end of the phase curve of Cd{s). Larger deviations can be noticed in the curves for the operational g-axis inductance Cq(s). Overall it is not clear if the iterated solution or the analytical results provide a better match with the measurements. However, a close resemblence of the values for the g-axis inductances can be found in Table 8.1. This indicates that the analytical parameters provide a good initial estimate which would be likely converge to the optimum solution during iterative refinement. Part III Simulation Exampl 123 Chapter 9 Simulation Case Studies In this chapter, two practical case studies are presented. The simulations are carried out with the new simulation program that resulted from this thesis research. The simulations are compared with reference solutions, obtained by the following two commercial power system simulation programs: M i C R O T R A f # The University of British Columbia's version of the E M T P . This version is mainly distributed among universities, but it is not significantly different from other E M T P versions. The program computes a steady-state solution prior to the time-domain simulation with instantaneous values or real signals. Being a reference for our simulations, the program is being used with very small simulation step sizes of At = 10 us. POWERFACTORY® Power system analysis package by DIGSILENT GmbH, Gomaringen, Germany. The program is widely distributed in Eastern Europe and South America. It calculates a power flow solution prior to time-domain simulation. The user can determine if the program should operate as an electromagnetic transients algorithm (EMT) using instantaneous values or as a stability program using RMS values. Besides a graphical user interface, the program offers many additional functions. More details can be found in [22] and on the internet under "http://www.digsilent.de". The PowERFACTORYprogram is only used for the subsynchronous resonance case study. There, it serves mainly as a representative of a transient stability program (RMS-option). 124 9.1. Transmission Line Energization 125 However, to ensure the accuracy of the results, we also make use of the EMT-option of this program and compare it with the results of the EMT-solution as well as with our own simulations. 9.1 Transmission Line Energizat ion A line energization for a C I G R E Working Group 13.05 test case [11] was simulated with the variable step size simulation program that resulted from this thesis research [41]. This case study was chosen in order to test and verify the new transmission line model of Chapter 6. The simulated system is shown in Fig. 9.1, and the parameters are listed in Table 9.1. A C 202.8 km @ V \ A / r^r^ (^ lo ea,b,c(V ^ T h e v ' "Thev ^ l ine' ^ l i n e Figure 9.1: Network diagram for C I G R E line energization test case. The feeding network is modeled with three ideal voltage sources behind a Thevenin impedance. The new simulation program expresses the source functions with complex sig-nals, hence Ea(t) = 1.0 e j { l O 0 t - 7 T / 2 ) p.u. (9.1) Eb(t) = 1 .0e j ( w o t + 5 7 r / 6 ) p.u. (9.2) Ec(t) = 1 . 0 e , ( ^ / 6 ) p.u.. (9.3) where UJ0 = 314.16 rad/s. The Thevenin impedance is Zxhev — 6.75 + j 127 Q in all three phases. The voltage across this impedance is given as V(t) = RThev I(t) + LThev jt I(t), where I(t) is the current through it. This differential equation is discretized by the program with the relaxed trapezoidal rule of integration for large time steps, (3.11) and (3.12). Thus: I(t) = YThevV(t) + H(t) (9.4) H{t) = FThevV(t-At) + GThevI{t-At) (9.5) 9.1. Transmission Line Energization 126 Table 9.1: C I G R E line energization case data. Equipment / Events Parameters Values Voltage sources frequency 50 Hz peak magnitude 1.0 p.u. phase shift in a 270° phase shift in b 150° phase shift in c 30° Thevenin impedance RThev 6.75 Q X-Thev 127.0 Q for all sequences Transmission line length I 202.8 km RLe,Pos 0.04 Q/km X\m^os 0.318 Q/km C[ine,Pos 11-86 nF /km Kne,zero 0.26 Q/km XLe,zero ^ tyl™ C\ine,zero 7.66 nF /km transposed, balanced Circuit breakers initial status open at £ = 3.05 ms in a closed at £ = 5.55 ms in c closed at £ = 8.05 ms in b closed where A t / At \ - i Thev — y-L'Thev ^ ^ Thev A J * Thev — 6 -^-[L'Thev 7^ T h e v J a n m = e - ^ f { l - f A y { l + * A ) , and where A = -L~lhevRThev - juQl. The circuit breakers in phases a, b, c close at 3.05 ms, 8.05 ms, and 5.55 ms, respectively. The voltage at the receiving end in phase b is shown in Fig. 9.2 (dotted line). Each dot represents a computed time step. The results were compared with a reference solution (solid line) obtained with M I C R O T R A N , with a constant step size of A£ = 10 /us, and A t = 5 ps during the switching instances. After the first switching instant, the simulation routine for the line model automatically chooses (6.8) and (6.9) for small time steps of 10 ps (At < r m i n ) . Except for the voltage 9.1. Transmission Line Energization 127 2.5 1.5 p 0.5 CU n hO u co Simulation results (0-40 ms) > -0.5 -1.5 -2 -2.5 ~ i r solid line: M I C R O T R A N dotted line: New program 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Time [s] 2.5 - 2 h -2.5 Simulation results (0-800 ms) -j i_ solid line: M I C R O T R A N dotted line: New program _ i i 0 0.1 0.2 0.3 0.4 0.5 Time [s] 0.6 0.7 0.8 Figure 9.2: Simulation results with a C I G R E line energization test case. 9.1. Transmission Line Energization 128 and current representation as complex numbers, the model is then identical to the constant parameter model of the EMTP [23, 24], with the same degree of accuracy. The graphs show the solution points as instantaneous values. For stability analysis it is common to use RMS voltages and currents instead. Since we chose complex signals for our variable representation, RMS values are contained in the magnitude of these signals and can be easily obtained at any instant. In Fig. 9.3, the simulation step size is monitored. It can be seen that the program decreases the step size immediately after the switching instant to Atmin = 5 fis for the backward Euler steps and then continues with 2Atmin = 10 /is. However, once the high frequency components of the electromagnetic transients disappear, the simulation step size is gradually increased until it reaches the user-defined maximum value Atmax = 50 ms at steady-state (circled dots). The lossless transmission line sections are then represented by (6.14) and (6.15) for all three modes. From a comparison with the reference solution in Fig. 9.2, it can be seen that the new transmission line model is also very accurate with large time steps. , C I G R E Line Energization Test Case 1 1 )• • I I I I-: New simulation program -. I : : : : : : : : : : : : : : : : : ; : : : : : : : : : : : ; : : ::::::::::::::::: ::: i : : : : : : : : : :-— :::::: : : : r: : ::::::::::: :::: : ::::: :::::::::: -: : : : : : MlCROTRAN 1 1 l l l l l 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time [s] Figure 9.3: Step size adjustment (0-800 ms). 9.2. Subsynchronous Resonance Case 129 9.2 Subsynchronous Resonance Case In this section, we investigate electromagnetic and electromechanical transient phenomena with the new simulation program. The IEEE subsynchronous resonance simulation bench-mark case [26, 47] was chosen for the testing of the new synchronous generator model of Chapter 7 and of the simultaneous integration of nodal, algebraic network equations (5.4) together with the state-space equations (5.5) for the machine variables. The subsynchronous resonance case study is also a useful example to illustrate differences in the results by elec-tromagnetic transients programs and stability analysis software. Figure 9.4: Network diagram for IEEE subsynchronous resonance benchmark case. A circuit diagram of this case is shown in Fig. 9.4. The data is presented in Table 9.2. For the representation of the step-up transformer, the excitation current and short-circuit losses are ignored. The transformer branch matrix then becomes L 1 u0 M V A rating 3X trafo V 2 VHV VLV -V-1 V~l vHV v LV VLV (9.6) for the high-voltage (HV) and low-voltage (LV) terminals of each phase. Considering voltage ratings as well as the connection type Y N d 01, the values of the above voltages are given as VHV = 500.0/v^ kV and VLV = 26.0 kV. For simplicity and in order to compare the results with a stability analysis of this case, 9.2. Subsynchronous Resonance Case 130 Table 9.2: IEEE subsynchronous resonance benchmark case data. Equipment / Events Parameters Values Voltage sources (infinite bus) frequency 60 Hz peak magnitude 361.9 kV phase shift in a 266.187° phase shift in b 146.187° phase shift in c 26.187° Thevenin reactance XThev 16.809 fi in all three phases Series compensation Cseries 25.5178 (J.F in all three phases Transmission line Rline,pos 5.603 fi X u „ s 140.07 fi Rline,zero 140.07 fi Xline,zero 437.02 fi transposed, balanced Step-up transformer M V A rating 892.4 M V A voltage rating H V 500.0 kV voltage rating LV 26.0 kV connection type Y N d 01 Xtrafo 0.14 p.U. three single-phase transformers Fault reactance Xfault 11.206 fi in all three phases Fault sequence initial status no fault at £ = 0.0 ms three-phase fault after £ = 75.0 ms clear fault at next current-zero continued on page 181 the transmission line model was reduced to a three-phase RL-series branch. An investigation of the effect of using different transmission line models can be found in [25]. 9.2.1 Steady-State Initialization A peculiarity of most time-step simulation programs is a very small oscillation in the machine variables at steady state. The reason for this oscillation is a numerical error caused by a different network representation for the steady-state initialization and for the time-domain simulation: For the initialization, the network is represented as an algebraic equation in the 9.2. Subsynchronous Resonance Case 131 Continuation of Table 9.2 Equipment / Events Parameters Values Synchronous generator M V A rating 892.4 M V A voltage rating 26.0 kV power factor 0.9 (lag.) connection type Y N , grounded 0.000 p.u. 0.130 p.u. xd 1.790 p.u. xq 1.710 p.u. x0 0.130 p.u. K 0.169 p.u. x\ 0.228 p.u. K 0.135 p.u. K 0.200 p.u. rp, 1do 4.300 s V qo 0.850 s rpll 1do 0.032 s rpll 0° 0.050 s H 2.894 s pole pairs 1 Excitation field rating if0 600 A constant voltage Turbine Power rating 803.16 M W constant torque, one mass frequency domain, e.g. (2.11). Under balanced conditions, this representation is absolutely exact and produces no error. In contrast, during a time-domain simulation, the differential equations associated with the transmission network are discretized with an integration algorithm such as the trape-zoidal rule. With finite step sizes, this discretization introduces a numerical error, as was investigated in Section 3.2. During a steady-state, time-domain simulation, this error prop-agates through the solution, thereby causing small oscillations. The shape and magnitude of these oscillations depends on the integration algorithm, as well as on the studied system. In general, the oscillations disappear in predominantly resistive networks and become more significant in inductive cases. Figures 9.5 and 9.6 compare the steady-state oscillations obtained with MICROTRAN 9.2. Subsynchronous Resonance Case 132 2.145 2.14 2.135 x 10 Comparison with M I C R O T R A N — 2.13F 2.125 2.12 2.115 2.11 2.105 l\ i >\ " l"1 l\ 1 J: l \ 1 | ,l: 1 1 1 . . 1 1 1 I ' • i» • 1 1 1 1 ' I I I 1 , / 1 \ : A X . V U L 1 | 1 i' i i 1 i 1 i' i ! ' , 1 ll II 1 1 1 1 ' w 1 1 w ll' 1 1 1 1 1 1 1 ' U i " )> l' \ > • ' i t V \ . . . v . . . . . . . . II • • • / • ! • v' 1 ' i \ 1 1 V i n ./ ' . . / h ^ !\; " i i 1 i • i i ' : i i l 1 ' i ' i i i 1 [ ' l , i H i ''' il " : i ' : i1 das soli hed lint d line: : : M I C R O T R A N ^ew program 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time [s] 2.135 2.134 2.133 2.132 o 2.131 2.13 2.129 2.128 2.127 2.126 2.125 x 10 Comparison with PoWERFACTORY (EMT) I I I I ! - -' ' ' ] ' • • I , I 1 L , , ' V 1 -' 1 ' : I r ' ' I I ; ' . I i ll 1 i Ii ll -u . M y i ' V • J 1 J -_l 1 J : ' — : ' J dashed line:: P O W E R F A C T O R Y solid line: New program i I I I I 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time [si Figure 9.5: Steady-state electromagnetic torque 9.2. Subsynchronous Resonance Case 133 (dashed line) and the new simulation program (solid line). The graphs show the mechanical torque and rotor speed of the synchronous generator. Notice that the new test program does not introduce any oscillations but produces the correct solution. The reason is that a network discretization with the relaxed trapezoidal algorithm or the relaxed convolution algorithm is free of numerical errors at steady state. The same oscillating phenomenon can be observed in a comparison with the POWER-FACTORYsoftware. The electromagnetic torque is shown in Fig. 9 .5 (dashed line). Using the electromagnetic transients setting (EMT) of this program, and the same step size as with M I C R O T R A N , that is, A t = 1 0 /JS, the program produces oscillations of much smaller ampli-tude. The reason may be the use of a different integration algorithm. The information what algorithm is used by the program has not been disclosed by DIgSILENT. The low resolution of this graph is caused by the truncation after the fourth decimal digit in this program's output file. Comparison with M I C R O T R A N 3601 | I I 1 1 1 1 1 1 1 3600.8 - -3600.6 - -3600.4 - -g. 3600.2 - ^ • • - i . - j . - ^ - : -i _ 3600 j : ; ; \ • CO \ . -g 3599.8 - -3599.6 - -3599.4 - -3599.2 h dashed line : M I C R O T R A N : solid line: I s'e\v program 35gg I 1 I I I I I . I I I I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [s] Figure 9 .6 : Steady-state rotor speed 9.2. Subsynchronous Resonance Case 134 9.2.2 Electromagnetic Transients Simulation After the steady-state initialization, we now apply a three-phase fault at t = 0 s. The fault is cleared after 75 ms when the instantaneous value of the fault current in each phase becomes zero. Notice that the PowERFACTORY program clears the fault at t = 75 ms simultaneously in all three phases. To compare the results of this simulation, we selected the electromagnetic transients option of this program. For the new test program of this thesis research, a minimum step size of Atmin = 5 /is and a maximum step size of Atmax = 16.67 ms, i.e. one cycle at power frequency, were specified. The aggressiveness (see Section 4.4) was kept at A = 0.5 for all simulations. The simulated time span was one second. The sudden short-circuit induces a subsynchronous resonance between the series com-pensator and the remaining network. The resulting torsional vibrations in the generator torque and rotor speed are shown in Fig. 9.7 and Fig. 9.8, respectively. We see that the peak generator torque of 7 x 10 6 Nm/rad occurs shortly after the fault clearing. This value corresponds to approximately 3.25 p.u. of the steady-state torque. The peak rotor speed occurs at t = 100 ms with a value of 3626 rpm. A comparison of the simulation results shows a very good match of the three programs. Next, the simulation time was increased to 10 s in order to study the long-term effect of this subsynchronous resonance phenomenon. The plots from M I C R O T R A N and our new simulation routine are shown on pages 137 and 138. It can be seen that, after approximately 5 s, the system in Fig. 9.4 has returned to steady state. The simulation results from the new program are practically identical with those obtained with M I C R O T R A N . However, while M I C R O T R A N used a constant simulation step size of At = 10 fis during this simulation time span, the step size of the new simulation program was continuously adjusted according to the dynamics in the system. This is shown in Fig. 9.11. These graphs illustrate how the simulation step size is increased after the fault, until the program uses the user-defined maximum step size of one full cycle, starting at t = 2 s. At this time, the program correctly identifies only stability type of transients in the solution and proceeds with Atmax for a simulation of the electromechanical transients. 9.2. Subsynchronous Resonance Case 135 Figure 9.7: Transients simulation (0-1 s): Electromagnetic torque. 9.2. Subsynchronous Resonance Case 136 Comparison with P O W E R F A C T O R Y ( E M T ) 3620 g_ 3600 3 o 3580 3570 1 i - / u i \ ?, / o \ 1 dash solic ed line: line: IS POWERFACTC ew program ):RY i i i I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [s] Figure 9.8: Transients simulation (0-1 s): Rotor speed. 9.2. Subsynchronous Resonance Case 137 Figure 9.9: Transients simulation (0-10 s): Electromagnetic torque. 9.2. Subsynchronous Resonance Case 138 Simulation with the new program 36301 1 1 1 1 1 1 1 ! r 3620 3570 I 1 1 1 1 1 1 1 1 1 1 0 1 2 3 4 5 6 7 8 9 10 Time [s] Simulation with M I C R O T R A N 3630 | i | 1 1 1 1 1 ! r 3620 3570 1 1 1 L — 1 1 1 1 1 1 0 1 2 3 4 5 6 7 8 9 10 Time [s] Figure 9.10: Transients simulation (0-10 s): Rotor speed. 9.2. Subsynchronous Resonance Case 139 10 10" IEEE Subsynchronous Resonance Case 10 10 10 • 1 1 • New simulation program M I C R O T R A N _l I I I I I 1 I L_ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time [s] 10 10 ' In 10" 10" 10 IEEE Subsynchronous Resonance Case 10 : New simulation program -j M I C R O T R A N _ i i i _ 0 1 2 3 4 5 6 7 8 9 10 Time [s] Figure 9.11: Transients simulation: Automatic step size adjustment. 9.2. Subsynchronous Resonance Case 140 From t = 0-300 ms, we can observe a relatively periodical switching between the two step sizes A t = 10 /xs and A t « 500 /xs. The step size increases whenever there is a local min imum or maximum in the rotor speed curve, because, at these values, the change of rotor speed and hence the change of instantaneous frequency is near zero. Yet, dur ing large acceleration or deceleration of the rotor, the instantaneous frequency also changes significantly and causes the simulation t ime step to drop to A t = 2Atmin. In comparison to Fig. 9.3 we also notice that step size changes more often than during the line energization test case. The reason is that the ratio of durations rsoiution for the solution of one t ime step and Treouud for the construction of the nodal admittance matr ix is much larger for the subsynchronous resonance case than for the line energization study. This is mainly caused by a number of Newton iterations for the solution of each t ime step. These iterations were not necessary for the line energization test case because no state-space equations were required there. 9.2.3 Stability Analysis Simulation Let us now compare the previous results w i th the results f rom a transient stabil i ty analysis. For this purpose, the setting of the POWERFACTORY program was changed to RMS for the transient stabi l i ty simulation. In addit ion, the simulation step size was increased to A t = 1.0 ms. The graphs for the generator torque and rotor speed are shown in Fig. 9.12 (dashed line) and in Fig. 9.13. Owing to the network representation at constant power frequency uo, its dynamic effects and contr ibut ion to the electromagnetic torque and rotor speed are missing. This can be found from the absence of the 2nd subharmonic component in the solution, f rom the instan-taneous drop of generator torque during the fault and from the lacking "backswing" effect in the rotor speed. Consequently, the peak electromagnetic torque is only about 1.3 p.u., which is 2.5 times less than the actual peak torque. In the remainder of the simulation, the mechanical torque and rotor speed by the stabil i ty program approximately bui ld the mean value to the sub-harmonic oscillations. After these oscillations have disappeared, the solution obtained from 9.2. Subsynchronous Resonance Case 141 Figure 9.12: Stability simulation (0-1 s): Electromagnetic torque and rotor speed. 9.2. Subsynchronous Resonance Case 142 3580 3570 Time [s 10 Figure 9.13: Stability simulation (0-10 s): Electromagnetic torque and rotor speed. 9.3. Evaluation of Program Performance 143 the stability program converges to the same steady-state solution as obtained with the other simulation programs. 9.3 Eva luat ion of P r o g r a m Performance The previous simulation examples have demonstrated that the new test program gave very accurate results. It was also shown that the step size was automatically adjusted in ac-cordance with the dynamics in the system. In this section we will compare the simulation performance of the three aforementioned simulation programs with respect to their running time. The running times required by the programs to simulate the test cases of this chap-ter are listed in Table 9.3. A l l simulations were carried out on a 200 MHz Pentium Pro computer. Table 9.3: Program running times for several test cases. Simulated case M I C R O T R A N P O W E R ] EMT-mode F A C T O R Y RMS-mode New program Line energization 0:23 min 0.3 ms/step n/a n/a 0:08 min 7.2 ms/step Steady-state simulation 1:29 min 0.9 ms/step 11:14 min 6.7 ms/step 2:13 min 133 ms/step 0:17 min 21.9 ms/step Subsynchronous resonance (1 s duration) 1:30 min 0.9 ms/step 11:14 min 6.7 ms/step 0:08 min 8 ms/step 16:23 min 208.2 ms/step Subsynchronous resonance (10 s duration) 14:31 min 0.9 ms/step n/a 0:28 min 28 ms/step 18:05 min 122.4 ms/step From this table it can be seen that the running times of the new program for the line energization test case as well as the steady-state simulation are superior to the times required by the other programs. However, the times needed for the completion of the subsynchronous resonance test cases are rather high. The reason is that the new program performs several Newton iterations on the solution of each simulation time step. Generally, only about 2-3 iterations are needed for the con-vergence of the pre-fault solution. However, during the fault and particularly after the fault clearing, the number of iterations may vary from 5 to 20 for each time step. This effect is 9 . 3 . Evaluation of Program Performance 144 by itself a heavy burden on the program's efficiency, yet it is aggravated by the use of very small time steps during the post-fault simulation, as was shown in Fig. 9.11. However, the performance of the new program is very much increased as soon as the simulation step size becomes larger again: For the simulation from 1-10 s the program requires less than two minutes. This timing is also very good compared to the 20 seconds needed by the stability program, considering the detailed representation of all three phases by the new simulation program. Chapter 10 Conclusion 10.1 Recommendations for Future Research In this thesis, a concept for the combined simulation of electromagnetic and electromechan-ical transient phenomena has been developed. The results shown in the previous chapter validate the applicability of this concept to time-domain power system simulation. This opens up the opportunity for several areas of future research work, of which the most pre-vailing ones are being addressed next. 10.1.1 Improved Simulation Speed In Section 9.3, the program's slow performance with small simulation step sizes was pointed out as a disadvantage. Albeit this does not in any way compromise the goals for this thesis research, it offers room for improvement. Since the code implemented for the test program is not optimized, it can probably be improved in several places. However, most of its efficiency is lost during the Newton iterations in combination with small step sizes. In the scope of power flow and stability analysis, there have been several contributions to improve Newton iterations [27, 28, 30, 71]. Whether or not these ideas can be applied to the simulation of electromagnetic transient phenomena needs to be investigated. Perhaps the Newton iterations can be replaced by other iteration schemes that do not rely on a Jacobian derivative matrix. The current iteration or Gauss-Seidel method in power flow computations would be an example for such schemes. Since no Jacobian matrix is required 145 10.1. Recommendations for Future Research 146 there, these methods are faster per iteration loop. However, generally they demand a higher number of iterations to arrive at the final solution point. It is therefore difficult to say if such methods could provide a faster solution. From the results in the previous chapter it can be seen that the E M T P yields very good and fast results by the variable prediction and correction method described in the E M T P Theory Book [24]. Similar ideas were presented in [52, 58]. No Newton iterations are required in these cases. Perhaps such solution could be applied to variable step sizes as well. 10.1.2 Component Modeling So far the program possesses models of linear RLC-branches, ideal voltage and current sources, circuit breakers, a constant parameter transmission line, and a synchronous machine. It is clear that many more models are needed in the future, particularly transformer models, load models, excitation and generation controls, as well as models to represent saturation effects. In Chapter 6, an extension of the new transmission line model to include frequency dependent parameters was suggested for future research. With the simultaneous solution of nodal equations and state-space equations, it is hoped that most excitation and generation control models can be imported from stability programs without difficulty. In addition, most existing E M T P models can be easily modified to allow for the new variable representation with complex signals and dynamic phasors. Intensive future research is envisioned for the inclusion of models for fast acting controls and power electronic equipment concerning H V D C transmission. Such devices continuously perform several switching actions per cycle, thereby inhibiting an efficient simulation with large simulation time steps. Despite this continuous switching activity, these electronic devices have a well defined steady-state behavior that does not depend on the simulation step size used. It thus seems feasible to develop a model that may properly reflect the steady-state and quasi-steady state behavior even if very large step sizes are used. A modeling technique that may be useful for this purpose is described in [66]. 10.1. Recommendations for Future Research 147 10.1.3 Online Network Aggregation The networks studied in a stability analysis are considerably larger than the networks of an electromagnetic transients study. With the automatic step sizes adjustment of Chapter 4, the gain through increased simulation steps during electromechanical transients can be used to investigate larger power systems. However, during the simulation of electromagnetic transient phenomena, very small sim-ulation step sizes are used so that large system sizes become a heavy burden to the solution process. Since typical electromagnetic transients only have a local impact, not all parts of the network need to be represented in detail and can therefore be aggregated. During the presence of slower dynamic oscillations, these aggregated networks must be inflated in order to accurately represent the system in a particular subarea. Network aggregation methods have been developed mainly for use in stability programs [4, 10, 74]. The applicability of these methods to electromagnetic transients types of analysis needs further investigation. To complicate matters, subarea network aggregation and in-flation would have to be performed online, i.e. during a simulation, since the transients propagate through parts of the system and generally change with time. 10.1.4 System Studies Time-domain system studies have so far been carried out with electromagnetic transients programs or with stability analysis programs. In some cases, however, a combination of both types of software is desirable. This lack of appropriate simulation tools was recognized in the context of power system restoration analysis [12, 13]. The main goal of this thesis was to develop the concepts for a combined simulation of electromagnetic and electromechanical power system transients. Even though the case studies of the previous chapter have demonstrated the applicability of these concepts, many more system studies need to be carried out in order to improve the new test program and develop it into a useful tool for a combined transients analysis, e.g. for system restoration studies. 10.2. Summary 148 10.2 S u m m a r y With regards to time-domain power system simulation, the areas of electromagnetic tran-sients studies and stability analysis have traditionally been isolated for several reasons. How-ever, recent system failures and power outages have created a demand for more comprehen-sive studies and more general simulation tools. In stability analysis, this has led to the combination of transient, mid-term and long-term stability simulation tools. Yet, difficulties concerning the numerical integration as well as modeling constraints have so far prevented an inclusion of electromagnetic transients simulation in this unifying trend. In this thesis, a concept has been developed to overcome these difficulties. The main achievements towards this goal are briefly summarized. A variable representation with complex signals and dynamic phasors was chosen and combined with the relaxed trapezoidal rule of integration in order to be able to simulate power system dynamics including the dynamics of the transmission network with large step sizes. The advantage of this variable representation compared to common three-phase simi-larity transformations is that it does not impose restrictions on the number of phases or on the symmetry of the network. Furthermore, this complex variable representation provides knowledge of the system's instantaneous frequency, which is important in determining the transient state of the system. Two integration algorithms based on the dynamic phasor representation have been de-rived. Their properties with respect to numerical errors as well as numerical stability have been thoroughly analyzed and compared with the properties of the backward Euler and trapezoidal algorithms, which are used in the E M T P [24, 63]. While both new algorithms prove to be numerically stable, only the relaxed trapezoidal algorithm permits step sizes larger than a cycle of power frequency. In order to adjust the simulation time step according to the system dynamics, a variable step size control was developed [42], which uses fuzzy logic and artificial intelligence for its decision making process. The advantage of this method over traditional techniques, which monitor the truncation error of an integration algorithm, is that the duration for building the nodal admittance matrix can be taken into account: While traditional methods are likely 10.2. Summary 149 to adjust the step size at almost each time step [31], the new approach may prevent a step size increase if its contribution to a faster solution is small compared to the time needed to rebuild the nodal admittance matrix. In the past, there have been discussions whether or not the state-space form of differential equations is preferable to a nodal formulation [24]. Both forms have certain advantages and disadvantages so that it was decided to use both nodal and state-space equations for the representation of the power system. This way, it is possible to implement component models with whatever formulation is best suited in the specific case. With the introduction of a complex variable representation, component models that used to rely on a description with instantaneous values had to be adapted accordingly. Existing transmission line models created a major challenge in this regard because models in electromagnetic transients studies do not permit step sizes larger than the smallest wave travel time constant. Contrarily, the 7r-model in stability analysis cannot properly reflect wave propagation phenomena. A new transmission line model [41] was therefore developed to overcome these difficulties. A practical line energization case study has demonstrated the accuracy of the new model in combination with large simulation step sizes. In Chapter 7, a synchronous generator model was derived for complex variables in all three phases. While the real part of the generator equations corresponds to the detailed generator models in stability analysis and electromagnetic transients studies, the difficulty was to define a corresponding imaginary part. Furthermore, since the concept of complex signals and dynamic phasors had been developed mainly for A C systems, a meaningful expansion of this concept for the DC rotor circuits was sought. This was accomplished with the new synchronous machine model for complex signals. In addition, an analytic method for the identification of the electrical machine parameters from standstill frequency response measurements has been proposed [40]. While existing methods are commonly based on a least square fitting of the model equations to the measured frequency characteristics, it is shown that the parameters can be analytically approximated without prior knowledge of the number of rotor circuits. This method was verified with the Lambton generator by Ontario Hydro [19, 55]. 10.2. Summary 150 Finally, the concepts developed in this thesis have been implemented in a simulation program, which was tested with a C I G R E line energization case study [11] and the IEEE subsynchronous resonance simulation benchmark case [47]. The results of these test cases demonstrate that the new program achieves the same accuracy as electromagnetic transients programs while the simulation step size is automatically adjusted according to the transient state of the system. Concluding it can be stated that the concepts and theory presented throughout this work are useful for detailed study cases involving long simulation times, where predominantly slow power system transients are present during most of the simulated time. Typical applications include off-line time-domain studies of medium-size systems, such as SSR studies and power system restoration. By not limiting the scope of this analysis to only the technical issues, such as the variable representation or integration method, the author has also improved and expanded existing modeling techniques for a wider frequency range. Aside from the difficulties regarding the program's run-time requirement for small step sizes, the simulations have validated the concepts presented in this thesis. The author hence considers this work a useful contribution to the combined study of electromagnetic and electromechanical transient phenomena. List of Symbols Alphanumeric Symbols Ms) Direct-axis transfer admittance Ag{s) Quadrature-axis transfer admittance A Aggressiveness factor for simulation Bd(s) Armature transfer admittance C(s) Rotor current to field tranfer function & [ . . . ] Dynamic phasor ea(t), eb(t) ec(t) Generator terminal voltages in phases abc. Ed(t), Eq(t), E0(t) Complex generator voltages in dqO coordinates efd(t) Generator field voltage fmax Maximum frequency Q(s) Armature-to-field transfer function Gc Characteristic wave conductance H(t-At) Vector of complex history values h Turns ratio between armature and field winding H Inertia constant @{I(t)} Current vector (dynamic phasor) Imaginary part IM Current vector (phasor) i(t) Current vector (real signal) I(t) Current vector (complex signal) ia(t), ib(t), ic(t) Generator terminal currents in phases abc. Id(t), Iq(t), I0(t) Complex generator currents in dqO coordinates yd{t) Generator field current j Imaginary unit Jin Moment of inertia V Gain of a transfer function (Chapter 8) t Length of a transmission line L Inductive matrix Ms) Direct-axis operational inductance £jd{s) Field winding operational inductance Quadrature-axis operational inductance Li Armature leakage inductance LaDq Quadrature-axis armature to rotor mutual inductance Lafd Direct-axis armature to rotor mutual inductance LDdK Direct-axis amortisseur inductance 151 L i s t o f S y m b o l s 152 Q u a d r a t u r e - a x i s a m o r t i s s e u r i n d u c t a n c e LfDd D i r e c t - a x i s field t o a m o r t i s s e u r m u t u a l i n d u c t a n c e lfd F i e l d i n d u c t a n c e r e f e r r e d t o t h e field Lfd F i e l d i n d u c t a n c e Lfkd P e r i p h e r a l l e a k a g e i n d u c t a n c e rn. C h a r a c t e r i s t i c m u l t i p l i e r ( C h a p t e r 3) Nd N u m b e r o f d - a x i s r o t o r w i n d i n g s np N u m b e r o f p o l e s . T w i c e t h e n u m b e r o f p o l e p a i r s Nq N u m b e r o f g - a x i s r o t o r w i n d i n g s NoF N u m b e r o f d i s c r e t e f r e q u e n c y p o i n t s V I n t e r p o l a t i o n f a c t o r ( C h a p t e r 6) V P o l e o f a t r a n s f e r f u n c t i o n ( C h a p t e r 8) q I n t e r p o l a t i o n f a c t o r ( C h a p t e r 6) R R e s i s t i v e m a t r i x » [ . . . ] R e a l p a r t Ra A r m a t u r e r e s i s t a n c e rfd F i e l d r e s i s t a n c e r e f e r r e d t o t h e field Rfd F i e l d r e s i s t a n c e S(u0) P h a s o r s L a p l a c e o p e r a t o r S(UJ) S p e c t r u m o f a s i g n a l s(t) R e a l s i g n a l S*(t) C o n j u g a t e o f a c o m p l e x s i g n a l t T i m e T O n e p e r i o d , c y c l e to S t a r t i n g t i m e , e .g . o f a s i m u l a t i o n tend S i m u l a t i o n e n d t i m e Tm{t) E l e c t r o m a g n e t i c g e n e r a t o r t o r q u e 9[V(t)] V o l t a g e v e c t o r ( d y n a m i c p h a s o r ) V(cu0) V o l t a g e v e c t o r ( p h a s o r ) v(t) V o l t a g e v e c t o r ( r e a l s i g n a l ) V(t) V o l t a g e v e c t o r ( c o m p l e x s i g n a l ) X(t) C o m p l e x s t a t e v a r i a b l e v e c t o r X D i s t a n c e o n a t r a n s m i s s i o n l i n e Z[...] Z - t r a n s f o r m a t i o n o f a f u n c t i o n z Z e r o o f a t r a n s f e r f u n c t i o n ( C h a p t e r 8) Jreek Symbols P A n g l e f r o m p h a s e a r e f e r e n c e a x i s t o d - a x i s A S m a l l c h a n g e o f a q u a n t i t y AAprediciCd P r e d i c t e d m a x i m u m m a g n i t u d e d e v i a t i o n At S i m u l a t i o n s t e p s i z e , t i m e s t e p A t m i n M i n i m u m s i m u l a t i o n t i m e s t e p Atpredicted P r e d i c t e d s i m u l a t i o n t i m e s t e p A t r n a x M a x i m u m s i m u l a t i o n s t e p s i z e t L o c a l n u m e r i c a l e r r o r List of Symbols 153 e Global numerical error e Tolerance margin for local numerical error r/(...) Error function 8(t) Angle from phase a reference axis to ci-axis 8i(t) Internal machine angle from phase a-axis to g-axis K, Index for generator rotor windings A Eigenvalue A Eigenvalue matrix p, Number of q-axis amortisseur windings, see Nq Pi Membership degree "Frequency gradient is low" p2 Membership degree "Frequency is close to nominal" p3 Membership degree "Number of samples is high" fj,4 Membership degree "Sampling window is large" pchange Membership degree "Likelihood of a step change is high" pdelay Membership degree "Duration since last step change is long Pefficiency Membership degree "A step change is efficient" pemtp Membership degree "Simulation is in EMTP mode" Pfiuct Membership degree "Frequency fluctuation is low" Pmagnitude Membership degree "The magnitude deviation is small" pmode Membership degree "Simulation mode is well defined" Preiiabie Membership degree "Forecast reliability is high" psample Membership degree "Sampling quality is high" pstab Membership degree "Simulation is in stability mode" pstep Membership degree "Step ratio is high" v Number of d-axis amortisseur windings Wj(t) Instantaneous frequency g(s) Residual function a Standard deviation of frequency samples T Modal wave propagation time Trebxdid Time required to rebuild the admittance matrix Tsoiution Time required to advance the solution by one time step TMVN Minimum modal wave propagation time 4> Phase angle, phase shift </>(...) Analytical solution ip(t) Partial solution to wave propagation equation (Chapter 6) ip(t) Partial solution to wave propagation equation (Chapter 6) '0a(t), 'tpb(t), '^cif) Generator fluxes in phases abc ^d{t), ^q{t), ^o{t) Complex generator fluxes in dqO coordinates (u) Average frequency LO Angular frequency wo System angular frequency, power frequency ^unear(t) Regression line through frequency samples upredicted,max Predicted maximum frequency u^dieted,mm Predicted minimum frequency upredicted Absolute predicted frequency u>m(t) Mechanical rotor speed List of Symbols 154 0Jr(t) Acronyms A C B E C I G R E I E E E E M T P N E T O M A C R C V R T R T N A T R U B C U B C Nominal frequency of a quantity Electrical rotor speed Alternating current Backward Euler algorithm International Conference on Large High Voltage Electric Systems, Conference Internationale des Grands Reseaux Electriques a Haute Tension Institute of Electrical and Electronics Engineers Electromagnetic Transients Program [6, 24] Networks, Torsion and Machines [50, 53] Relaxed convolution algorithm Relaxed trapezoidal algorithm Transient network analyzer Trapezoidal algorithm University of British Columbia University of British Columbia Bibliography [1] F. L. Alvarado, "Formation of Y-Node Using the Primitive Y-Node Concept," IEEE Transactions on Power Apparatus and Systems, vol. PAS-101, no. 12, pp. 4563-4571, December 1982. [2] P. M . Anderson and A . A . Fouad, Power System Control and Stability, Power System Engineering Series. IEEE Press, Piscataway, NJ , U.S.A., 1994, ISBN 0-7803-1029-2. [3] J. Arillaga and C. P. Arnold, Computer Analysis of Power Systems, John Wiley & Sons, Ltd., New York, N Y , U.S.A., 1990. [4] D. T . Bargiotas and J. S. Lawler, "Effects of Aggregation Methods on Individual Modes in Reduced Order Power System Models," in South East Conference, 1988, IEEE Conference Proceedings, pp. 579-586. [5] L. Bergeron, Water Rammer in Hydraulics and Wave Surges in Electricity, John Wiley & Sons, Ltd., New York, N Y , U.S.A., 1949, (1961). [6] Bonneville Power Administration, BPA, Portland, Oregon, Electromagnetic Transients Program Rule Book, Apri l 1982. [7] I. N . Bronstein and K . A . Semendjajew, Taschenbuch der Mathematik, Verlag Harri Deutsch, Frankfurt/M., Germany, 6th edition, May 1966, Edition Leipzig, Germany. [8] I. M . Canay, "Causes of Discrepancies on Calculation of Rotor Quantities and Exact Equivalent Diagram of Synchronous Machines," IEEE Transactions on Power Appara-tus and Systems, vol. PAS-88, pp. 1114-1120, July 1969. [9] F. Castellanos, Full Frequency-Dependent Phase-Domain Modeling of Transmission Lines and Corona Phenomena, Ph.D. thesis, The University of British Columbia, Department of Electrical Engineering, Faculty of Applied Science, Vancouver, B.C. , Canada, Apri l 1997. [10] J. H. Chow, Ed., Time-Scale Modeling of Dynamic Networks with Applications to Power Systems, Lecture Notes in Control and Information Sciences. Springer-Verlag Berlin-Heidelberg-New York, 1982, ISBN 3-540-12106-4. [11] C I G R E Working Group 13.05, "The Calculation of Switching Surges (I). A Comparison of Transient Network Analyzer Results," Electra, vol. 19, pp. 67-78, 1971. 155 B I B L I O G R A P H Y 156 [12] C I G R E Working Group 38.02, T F 02, "Modelling and Simulation of Black Start and Restoration of Electric Power Systems: Results of a Questionnaire," Electra, vol. 131, pp. 156-169, 1992. [13] C I G R E Working Group 38.02, T F 02, "Modelling and Simulation of Black Start and Restoration of Electric Power Systems," Electra, vol. 147, pp. 21-42, 1993. [14] E. Clarke, Circuit Analysis of AC Power Systems, vol. 1, John Wiley & Sons, Inc., New York, N Y , U.S.A., 1943. [15] L. Cohen, Time-Frequency Analysis, Prentice Hall Signal Processing Series. Prentice Hall P T R , Englewood Cliffs, NJ 07632, U.S.A., 1995, ISBN 0-13-594532-1. [16] M . E. Coultes and W. Watson, "Synchronous Machine Models by Standstill Frequency Response Tests," IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 4, pp. 1480-1489, Apri l 1981. [17] E. Cox, The Fuzzy Systems Handbook, Academic Press, Inc., Cambridge, M A , U.S.A., 1994. [18] P. L. Dandeno, P. Kundur, A . T. Poray, and H. M . Zein El-Din, "Adaptation and Validation of Turbogenerator Model Parameters through On-Line Frequency Response Measurements," IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 4, pp. 1656-1664, Apri l 1981. [19] P. L. Dandeno and A. T. Poray, "Development of Detailed Turbogenerator Equivalent Circuits from Standstill Frequency Response Measurements," IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 4, pp. 1646-1655, Apri l 1981. [20] P. Denbigh, System Analysis & Signal Processing, Addison Wesley Longman Ltd., Edinburgh Gate, Harlow, Essex CM20 2JE, U .K. , 1st edition, 1998, ISBN 0-201-17860-5. [21] Deutsches Institut fur Normung e. V . , Berlin, Germany, DIN 13321 Komponenten in Drehstrornnetzen, September 1980. [22] DIgSILENT GmbH, Gomaringen, Germany, DIgSILENT Power System Analysis, Ver-sion 10.31, 1997. [23] H. W. Dommel, "Digital Computer Solution of Electromagnetic Transients in Single and Multiphase Networks," IEEE Transactions on Power Apparatus and Systems, vol. PAS-88, pp. 388-399, Apri l 1969. [24] H. W. Dommel, EMTP Theory Book, Microtran Power System Analysis Corporation, 4689 W. 12th Avenue, Vancouver, B .C . V6R 2R7, Canada, 2nd edition, May 1992. [25] H. W. Dommel, "Introduction to the Use of MicroTran and Other E M T P Versions," Course Notes (ELEC 553), Dept. of Electrical & Computer Engineering, U B C , Vancou-ver, B.C. , Canada, 1998, available at http://www.ece.ubc.ca/power/e553.pdf. B I B L I O G R A P H Y 157 H. W. Dommel, J. R. Marti, and L. Marti, "Subsynchronous Resonance: Fact Sheet No. 1," Tech. Rep., MicroTran Power System Analysis Corporation, 4427 W. 6th Av-enue, Vancouver, B.C. , Canada V6R 1V2, February 1991. H. W. Dommel and N . Sato, "Fast Transient Stability Solutions," IEEE Transactions on Power Apparatus and Systems, vol. PAS-137, pp. 1643-1650, July/August 1972. Electric Power Research Institute, "Power System Dynamic Analysis - Phase I," Final Report EPRI EL-484 (Research Project 670-1), Palo Alto, C A 94304, U.S.A., July 1977. Electric Power Research Institute, E M T P Development Coordination Group, Palo Alto, C A , U.S.A., Electromagnetic Transients Program (EMTP): Rule Book, Apri l 1986, EPRI EL-4541-CCMP. V. V . Ermakov and N . N . Kalitkin, "The Optimal Step and Regularization for Newton's Method," U.S.S.R. Comput. Maths. Math. Phys., vol. 21, no. 2, pp. 235-242, 1981. 0. B. Fosso, "Implementation of an Algorithm for Numerical Integration (Gear's Method)," Project Report, Power System Planning Division, Ontario Hydro, 700 Uni-versity Avenue, Toronto, Canada, June 1993. R. J. Frowd, J. C. Giri , and R. Podmore, "Transient Stability and Long-Term Dynamics Unified," IEEE Transactions on Power Apparatus and Systems, vol. PAS-101, pp. 3841-3850, October 1982. D. Gabor, "Theory of Communication," Journal ofthe IEE, vol. 93, pp. 429-457, 1946. W. C. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Inc., Englewood Cliffs, NJ 07632, U.S.A., 1971. R. W. Hamming, Digital Filters, Prentice-Hall Signal Processing Series. Prentice-Hall, Inc., Englewood Cliffs, NJ 07632, U.S.A., 2nd edition, 1983, ISBN 0-13-212506-4. E. Handschin, Elektrische Energieiibertragungssysteme, Teil 1: Stationarer Betriebszu-stand, Huthig-Verlag, Heidelberg, Germany, 1983. H. Happoldt and D. Oeding, Elektrische Kraftwerke und Netze, Springer-Verlag Berlin-Heidelberg-New York, 5th edition, 1978, ISBN 3-540-08305-7. M . L. Harris, P. J. Lawrenson, and J. M . Stephenson, Per-Unit Systems with Special Reference to Electric Machines, IEE Monograph Series 4. Cambridge University Press, London, U .K. , 1970. S. Henschel, "Long-Term Evaluation for Power Systems under Consideration of Elec-tromagnetic and Electromechanical Transient Phenomena," Ph.D. thesis proposal, Uni-versity of British Columbia, Department of Electrical Engineering, Vancouver, B.C. , Canada, November 1994. S. Henschel and H. W. Dommel, "Noniterative Synchronous Machine Parameter Identi-fication from Frequency Response Tests," IEEE Transactions on Power Systems, May 1998, accepted for publication, ref. no. PE-003-PWRS-0-05-1998. B I B L I O G R A P H Y 158 [41] S. Henschel, A . I. Ibrahim, and H. W. Dommel, "Transmission Line Model for Variable Step Size Algorithms," International Journal of Electrical Power & Energy Systems, vol. 21, no. 3, pp. 191-198, January 1999. [42] S. Henschel, D. Lindenmeyer, T. Niimura, and H. W. Dommel, "Intelligent Step Size Control for Fast Time-Domain Power System Simulation," in Power Systems Compu-tation Conference, Trondheim, Norway, 1999, accepted for publication, ref. no. 22. [43] J. D. Hurley and H. R. Schwenk, "Standstill Frequency Response Modeling and Evalu-ation by Field Tests on a 645 M V A Turbine Generator," IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 2, pp. 828-836, February 1981. [44] IEEE, "Standard 115-1983," in IEEE Guide: Test Procedures for Synchronous Ma-chines. Institute of Electrical and Electronics Engineers, Inc., September 1983. [45] IEEE Joint Group, Chairman: P. L. Dandeno, "Supplementary Definitions and Associ-ated Test Methods for Obtaining Parameters for Synchronous Machine Stability Study Simulations," IEEE Transactions on Power Apparatus and Systems, vol. PAS-99, no. 4, pp. 1625-1633, July/August 1980. [46] IEEE Standard 115A-1987, "IEEE Standard Procedures for Obtaining Synchronous Machine Parameters by Standstill Frequency Response Testing," in IEEE Standards Collection 'Electric Machinery". Institute of Electrical and Electronics Engineers, Inc., August 1995, ISBN 1-55937-533-7. [47] IEEE Subsynchronous Resonance Task Force of the Dynamic System Performance Working Group, Power System Engineering Committee, "First Benchmark Model for Computer Simulation of Subsynchronous Resonance," IEEE Transactions on Power Apparatus and Systems, vol. PAS-96, no. 5, pp. 1565-1572, September/October 1977. [48] A. Inan, "Integration von Raumzeigergrofien mit faktorisierten Eigenwerten," Etz-Archiv 3, vol. 8, pp. 279-281, May 1.981. [49] E. W. Kimbark, Power System Stability, vol. 3. Synchronous Machines, John Wiley & Sons, Ltd., New York, NY, U.S.A., 1956. [50] B. Kulicke, "Digitalprogramm N E T O M A C zur Simulation elektromechanischer und -magnetischer Ausgleichsvorgange in Drehstromnetzen," Elektrizitatswirtschaft 78, vol. 1, pp. 18-23, 1979. [51] B. Kulicke, "Simulationsprogramm N E T O M A C : Differenzenleitwertverfahren bei kon-tinuierlichen und diskontinuierlichen Systemen," Siemens Forschungs- und Entwick-lungsberichte, vol. 10, no. 5, pp. 299-302, 1981. [52] B. Kulicke, "Simulationsprogramm N E T O M A C : Nachbildung von Synchron- und Asyn-chronmaschinen," Siemens Forschungs- und Entwicklungsberichte, vol. 11, no. 3, pp. 156-161, 1982. [53] B. Kulicke, "Programmsystem zur Berechnung elektromechanischer und -magnetischer Ausgleichsvorgange," Etz-Archiv 6, vol. 12, pp. 439-443, 1984. B I B L I O G R A P H Y 159 [54] B. Kulicke, "Digitale Simulation elektromagnetischer und -mechanischer Ausgleichs-vorgange," Course Notes, Institut f. elektrische Energietechnik, T .U . Berlin, Berlin, Germany, 1991. [55] P. Kundur, Power System Stabiliy and Control, EPRI Power System Engineering Series. McGraw-Hill, Inc., New York, N Y 10020, U.S.A., 1994, ISBN 0-07-035958-X. [56] K . Kupfmuller, Einfiihrung in die theoretische Elektrotechnik, Springer-Verlag Berlin-Heidelberg-New York, 10th edition, 1973, ISBN 3-540-06021-9. [57] N . Levanon, Radar Principles, John Wiley & Sons, Inc., 1st edition, 1988, ISBN 0-471-85881-1. [58] U . Linnert, Berechnung von Ausgleichsvorgangen in Elektroenergiesystemen unter Ver-wendung grofitmoglicher Schrittweiten, Ph.D. thesis, Universitat Erlangen-Niirnberg, Technische Fakultat, Erlangen, Germany, 1995. [59] U . Linnert, "Correct and Fast Network Simulation Using the New Calculation Method SDLV3 with Advancing Implied Euler Step," in IPST '95 International Conference on Power System Transients, Lisbon, Portugal, September 1995, pp. 215-220. [60] F. J. Marcano, "Modelling of Transmission Lines Using Idempotent Decomposition," M.S. thesis, The University of British Columbia, Department of Electrical Engineering, Faculty of Applied Science, Vancouver, B.C. , Candada, August 1996. [61] J. R. Marti, The Problem of Frequency Dependence in Transmission Line Modelling, Ph.D. thesis, The University of British Columbia, Department of Electrical Engineering, Faculty of Applied Science, Vancouver, B.C. , Canada, April-1981. [62] J. R. Marti, "Accurate Modelling of Frequency-Dependent Transmission Lines in Elec-tromagnetic Transients Simulations," IEEE Transactions on Power Apparatus and Sys-tems, vol. PAS-101, pp. 147-157, January 1982. [63] J. R. Marti and J. Lin, "Suppression of Numerical Oscillations in the E M T P , " IEEE Transactions on Power Systems, vol. PWRS-4, no. 2, pp. 739-747, May 1989. [64] L. Marti, Simulation of Electromagnetic Transients in Underground Cables with Frequency-Dependent Modal Transformation Matrices, Ph.D. thesis, The University of British Columbia, Department of Electrical Engineering, Faculty of Applied Science, Vancouver, B.C. , Canada, November 1986. [65] R. M . Mathur and R. K . Varma, "Modeling Effects of System Frequency V ariation in Transient Stability Studies," Final Report P.O. 957595726-11 (contract), Power System Planning Division, Ontario Hydro, 700 University Avenue, Toronto, Canada, September 1993. [66] P. Mattavelli, G. C. Verghese, and A. M . Stankovic, "Phasor Dynamics of Thyristor-Controlled Series Capacitor Systems," IEEE Transactions of Power Systems, vol. 12, no. 3, pp. 1259-1267, August 1997. B I B L I O G R A P H Y 160 S. H . Minnich, "Small Signals, Large Signals, and Saturation in Generator Modeling," IEEE Transactions on Energy Conversion, vol. EC-1, no. 1, pp. 94-102, March 1986. H. V . Nguyen, Simulation of Lightning Surges on Transmission Lines, Ph.D. thesis, The University of British Columbia, Department of Electrical Engineering, Faculty of Applied Science, Vancouver, B.C. , Canada, February 1996. A. H. Nuttall, "On the Quadrature Approximation to the Hilbert Transform of Modu-lated Signals," Proceedings ofthe IEEE (Letters), vol. 54, pp. 1458-1459, 1966. T. S. Parker and L. O. Chua, Practical Numerical Algorithms for Chaotic Systems, Springer-Verlag New York-Heidelberg-Berlin, 1989, ISBN 0-387-96689-7. C. D. Pawloski, M . H. Khammash, V . Vittal , and K . Zarley, "Reduction in the Compu-tation Time for Robust Stability Analysis of Power Systems," in 28th North American Power Symposium, M.I.T., Cambridge, M A , U.S.A., November 1996, pp. 467-474. R. C. Pfaffenberger and J. H. Patterson, Statistical Methods for Business and Eco-nomics, Richard D. Irwin, Inc., Homewood, IL, U.S.A., 1977, ISBN 0-256-01797-2. W. H . Press, B. P. Flannery, S. A . Teukolsky, and W. T. Vetterling, Numerical Recipes - The Art of Scientific Computing, Cambridge University Press, Cambridge CB2 1RP, Great Britain, 1986, ISBN 0-521-30811-9. A. D. Ranjit and J. H. Chow, "Aggregation Properties of Linearized Two-Time-Scale Power Networks," IEEE Transactions on Circuits and Systems, vol. 38, no. 7, pp. 720-730, July 1991. A. W. Rankin, "Per Unit Impedance of Synchronous Machines," AIEE Transactions, vol. 64, August 1945. J. V . Sarlashkar, F. L. Alvarado, and J. Mahseredjian, "Interfacing Time-Domain and Phasor-Domain Simulations," in 31st Universities Power Engineering Conference, Irakleion, Greece, September 1996. R. P. Schulz, W. D. Jones, and D. N . Ewart, "Dynamic Models of Turbine Generators Derived from Solid Rotor Equivalent Circuits," IEEE Transactions on Power Apparatus and Systems, vol. PAS-92, pp. 926-933, May/June 1973. H. H. Skilling, Electric Networks, John Wiley k Sons, Inc., New York, NY, U.S.A., 1974, ISBN 0-471-79420-1. C. P. Steinmetz, Theoretical Elements of Electrical Engineering, McGraw-Hill, Inc., 1915. M . Stubbe, A . Bihain, J. Deuse, and J. C. Badder, "STAG: A New Unified Software Program for the Study of the Dynamic Behavior of Electrical Power Systems," IEEE Transactions on Power Systems, vol. PWRS-4, no. 1, pp. 129-138, February 1989. K . Tanaka, An Introduction to Fuzzy Logic for Practical Applications, Springer-Verlag New York-Heidelberg-Berlin, 1997, ISBN 0-387-94807-4. B I B L I O G R A P H Y 161 [82] S. D. Umans, J. A. Mallick, and G. L. Wilson, "Modelling of Solid Rotor Turbogenera-tors," IEEE Transactions on Power Apparatus and Systems, vol. PAS-97, pp. 269-291, January/February 1978. [83] D. E. Vakman and L. A . Vainshtein, "Amplitude, Phase, Frequency—Fundamental Concepts of Oscillation Theory," Soviet Phys. Uspekhi, vol. 20, pp. 1002-1016, 1978. [84] V . Venkatasubramanian, "Tools for Dynamic Analysis of the General Large Power System Using Time-Varying Phasors," International Journal of Electrical Power & Energy Systems, vol. 16, no. 6, pp. 365-376, 1994. [85] A. Walton, "Characteristics of Equivalent Circuits of Synchronous Machines," IEE Procdeedings on Electrical Power Applications, vol. 143, no. 1, pp. 31-40, January 1996. [86] L. M . Wedepohl, "Application of Matrix Methods to the Solution of Traveling Wave Phenomena in Polyphase Systems," IEE Proceedings, vol. 110, pp. 2200-2212, December 1963. [87] J. H . Wilkinson and C. Reinsch, Linear Algebra, vol. 2 of Handbook for Automatic Computation Series, Springer-Verlag New York-Heidelberg-Berlin, 1971. [88] R. Zurmiihl, Praktische Mathematik fur Ingenieure und Physiker, Springer Verlag, Berlin-Heidelberg, Berlin, Germany, 1965.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of electromagnetic and electromechanical power...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of electromagnetic and electromechanical power system transients with dynamic phasors Henschel, Sebastian 1999
pdf
Page Metadata
Item Metadata
Title | Analysis of electromagnetic and electromechanical power system transients with dynamic phasors |
Creator |
Henschel, Sebastian |
Date Issued | 1999 |
Description | Over the last 50 years, digital simulation of electric power systems has become an integral part for planning, design and operation in the power industry. The number of possibilities with respect to the purpose of a study, investigated frequency ranges, etc. with, in the past, limited computer resources has resulted in a spectrum of simulation tools, designed to handle very specific tasks. Simplifying assumptions were often needed to facilitate such a simulation. Recent system failures and power outages, partly due to increasingly sensitive operating conditions, have created a demand for more comprehensive studies and more general simulation tools that overcome former limitations. With regards to time-domain simulation, this demand has led to combining the areas of transient, mid-term and long-term stability. Confronted with concerns about black start and system restoration due to a global trend to a deregulated power market, several power utilities suggested to also include the area of electromagnetic transients. However, previously made assumptions as well as technical limitations complicate the implementation of this idea: Stability programs are based on the assumption that power transfer takes place at system frequency and are therefore unable to represent rapid electromagnetic transients. Electromagnetic transients programs, on the other hand, are very accurate but use too small simulation step sizes for an efficient simulation of electromechanical transients. A new method for simulating both types of transient phenomena with complex signals and dynamic phasors is presented in this thesis. Whereas in previous work three-phase transformations had been used to accomplish this task, this new method is applied directly in the phase-domain and not restricted to balanced three-phase systems. Several numerical aspects such as an appropriate variable representation, integration method and a control mechanism for variable simulation step sizes have been addressed. Important modeling constraints for transmission lines and synchronous machines have been thoroughly analysed and are considered in this work. Finally, several simulation examples are presented that compare the results of the new simulation program with the results obtained from the EMTP. These case studies demonstrate that the methods presented in this thesis are suited to combine the areas of electromagnetic and electromechanical transients simulation. |
Extent | 8095772 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065161 |
URI | http://hdl.handle.net/2429/9908 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1999-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1999-388948.pdf [ 7.72MB ]
- Metadata
- JSON: 831-1.0065161.json
- JSON-LD: 831-1.0065161-ld.json
- RDF/XML (Pretty): 831-1.0065161-rdf.xml
- RDF/JSON: 831-1.0065161-rdf.json
- Turtle: 831-1.0065161-turtle.txt
- N-Triples: 831-1.0065161-rdf-ntriples.txt
- Original Record: 831-1.0065161-source.json
- Full Text
- 831-1.0065161-fulltext.txt
- Citation
- 831-1.0065161.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065161/manifest