Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analysis of electromagnetic and electromechanical power system transients with dynamic phasors Henschel, Sebastian 1999

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

A n a l y s i s of Electromagnetic a n d E l e c t r o m e c h a n i c a l P o w e r S y s t e m Transients w i t h D y n a m i c Phasors by SEBASTIAN HENSCHEL 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 O F THE REQUIREMENTS FOR T H E D E G R E E OF  Doctor of Philosophy in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Electrical & Computer Engineering)  We accept this thesis as conforming to the required standard  The University of British Columbia February 1999 © Sebastian Henschel, 1999  in  presenting  degree freely  at  this  the  thesis  in  partial  fulfilment  University  of  British  Columbia,  available for  copying  of  department publication  this or of  reference  thesis by  this  for  his thesis  and study. scholarly  or for  her  The University of British Vancouver, Canada  DE-6 (2/88)  Columbia  I further  purposes  the  requirements  I agree  that  agree  may  representatives.  financial  permission.  of  It  gain shall not  be is  that  the  Library  permission  granted  by  understood be  for  an  advanced  shall for  the that  allowed without  make  it  extensive  head  of  my  copying  or  my  written  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 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. ii  Important modeling  ABSTRACT  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  EMTP.  These case studies demon-  strate that the methods presented in this thesis are suited to combine the areas of electromagnetic and electromechanical transients simulation.  Contents Abstract  ii  L i s t of Tables  viii  L i s t of F i g u r e s  ix  Acknowledgement 1  xii  Introduction  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 Variable Representation  12  2.1  2.2  2.3  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  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  Spectral Properties of Complex Signals  27  2.3.1  Spectral Averages  27  2.3.2  Instantaneous Frequency  29 iv  I  CONTENTS 2.3.3  Peculiarities Regarding Complex Signals  3 Numerical Integration 3.1  3.2  3.3  4.2  4.3  4.4  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  Integration Errors  39  3.2.1  Accuracy  39  3.2.2  Local Errors  45  3.2.3  Global Errors  49  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  55  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  Prediction Reliability  61  4.2.1  Sampling Quality  62  4.2.2  Frequency Fluctuation  63  4.2.3  Simulation Mode  63  Step Change Efficiency  64  4.3.1  Step Ratio  65  4.3.2  Time Delay  66  Decision Making Process  5 Structure of the Simulation Process 5.1  33  Integration with Complex Signals  4 Variable Integration Step Sizes 4.1  31  Overall System Equations  67  69 69  CONTENTS  II 6  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  Component Modeling  81  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  Lossy Line Model  88  Synchronous Machine  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  8  80  Transmission Line  6.3 7  XI  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  Synchronous Machine Parameters  108  8.1  Frequency Response Tests  109  8.2  Transfer Functions  Ill  8.2.1  The Q-Axis Admittance  Ill  8.2.2  The D-Axis Admittance  112  8.2.3  Parameter Extraction  114  8.3  Practical Considerations  117  CONTENTS  8.4  III 9  vii  8.3.1  Pre-Processing Filter  117  8.3.2  Adjustment of L  117  8.3.3  Relating C (s) and sG{s)  118  fkd  d  Results  120  Simulation Examples  123  S i m u l a t i o n 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  10 C o n c l u s i o n 10.1 Recommendations for Future Research  143 145 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  L i s t of S y m b o l s  151  Bibliography  155  L i s t o f 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  Truncation errors for the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algorithms  46  3.2  8.1  Comparison of identified electrical synchronous machine parameters  122  9.1  C I G R E line energization case data  126  9.2  I E E E 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 relaxed trapezoidal (RTR) rules of integration  43  Numerical accuracy of the backward Euler (BE), trapezoidal (TR), and relaxed convolution (RCV) rules of integration  44  3.4  Local integration errors for u> > UJ  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 convolution 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  3.3  0  ix  x  LIST O F F I G U R E S 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 ; „  85  6.2  Lossless transmission line model for At > r  88  6.3  Schematic of a lossy transmission line model  89  7.1  Synchronous generator d and g-axis models  92  8.1  Admittance A (s) of the Lambton generator and its final approximation by a third-order polynomial  112  Admittance A (s) of the Lambton generator approximated by a first-order function of type Ti(s)  114  8.2  m  min  q  q  8.3  First residue of A (s)  .  116  8.4  Parameters Lfkd and lf as functions of frequency. The cross indicates the best fit, yielding L = -0.0282 p.u. and l = 0.831 p.u  118  Best concurrence of functions RIGHT indicates the optimum value of Lj  119  q  approximated by a first-order function of type ^(s). d  fkd  8.5  fd  (solid line) and LEFT  (dashed line)  kd  8.6  Measured operational d-axis inductance Cd(s) with its second-order analytical (solid line) and iterated approximations (dashed line)  120  Measured armature-to-field transfer function sG{s) with its second-order analytical (solid line) and iterated approximations (dashed line)  121  Measured operational (/-axis inductance C (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 I E E E subsynchronous resonance benchmark case.  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  8.7 8.8  q  9.10 Transients simulation (0-10 s): Rotor speed 9.11 Transients simulation: Automatic step size adjustment 9.12 Stability simulation (0-1 s): Mechanical torque and rotor speed  . . . 129  138 . 139 141  LIST O F F I G U R E S 9.13 Stability simulation (0-10 s): Mechanical torque and rotor speed  xi  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 analysis, 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 becoming 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 concerning 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 simulation 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 continuously 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 N S E R C 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 foremost 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 February 9, 1999  Sebastian Henschel  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 encountered 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 configuration, which may be caused by the closing or opening action of circuit breakers or  1  1.1. General Background  2  Frequency f in Hz  co> u c  <5 Cc  o c  SF - 10  6  - 10  5  6  transients  CD  Wave propagation, lightning  10  4  \- 10  3  E o C  CU  s. Q.  CT fi  n c  )  Switching overvoltages Transformer saturation  CD  c o> re E o u 0)  10  2  Steady-state power flow I- 10 CM <D O)  Subsynchronous resonance  1  Transient stability: machine rotor dynamics  10°  c  <5 cc  Cr fi  0) (0  10-  Interarea oscillations  1  o c  (0  o c  CD  c u E o c  I- IO"  2  "IO" 3  '— 10"  Mid-term and long-term stability: automatic generation control  u 0) E o o  0)  4  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 sufficient. 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.  3  1.1. General Background  The study of electromagnetic transient phenomena includes switching surges, transient recovery voltage, and ferroresonance.  Electromagnetic transients programs ( E M T P ) are 1  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. disturbance.  during the first few seconds after a  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, frequency 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. I n this text, this acronym is used as a general term for any digital program that can simulate electromagnetic 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  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 remaining at power frequency. It is immediately clear that, with this assumption, electromagnetic 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  R e v i e w of U n i f i e d T i m e - D o m a i n S i m u l a t i o n  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 state . 2  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 tranIn 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]. 2  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.  6  1.2. Review of Unified Time-Domain Simulation  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 facilitate 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  NETOMAC  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  7  1.2. Review of Unified Time-Domain Simulation 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 quasisteady 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 representation 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 e  jult  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 a n d O b j e c t i v e  With the availability of powerful computer hardware as well as improved analytical techniques, 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 unjustified 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 system 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 electromagnetic 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 electromechanical 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 following 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 intelligence 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.  Thesis Motivation and  from standstill frequency response  In  Chapter  9 the  10  Objective  concept  that  measurements.  has been  developed  in this research  w o r k is t r i e d  o n p r a c t i c a l case studies.  T h e results are e v a l u a t e d a n d c o m p a r e d to a reference  obtained w i t h the E M T P .  T h e last chapter contains conclusions d r a w n f r o m this work,  w e l l as r e c o m m e n d a t i o n s for f u t u r e r e s e a r c h w o r k .  out  solution as  Part I N u m e r i c a l Considerations  11  Chapter 2 Complex Variable Representation 2.1  E v a l u a t i o n 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, At  max  1 — ——-  ,  (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 tv (t) = -j- co {t) r  m  (2.3)  where n is the number of machine poles, and uj is the mechanical rotor speed. The angle of p  m  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  14  2.1. Evaluation of Existing Variable Representations  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 s {t) s {t)  cos ( 0 - ^ )  COS0  d  q  -sin#  | =  cos(fl+^)  -sin(0-f)  sin(0+f)  (2.5)  s (t) 0  Here, s (t), Sb(t), s (t) can be taken as the abc components of any quantity, e.g. voltage, a  c  current, flux, etc. The inverse transformation is given by cos 9  -sin0  cos (6 — 3 cos  1 ~  /  •sin(0-f)  1  {9+\) ,  sin(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)  —  to t in order to filter out the power frequency component U)Q. This 0  implies that the quantities s (t), a  i.e. when  Sb(t), and s (t) take on constant values at steady state, c  During rotor oscillations, of course, these values vary according to the  frequency deviation Au;  15  2.1. Evaluation of Existing Variable Representations  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  s (t) p  I =  1  - \  0  f  1 2  s (t) 0  _ 1 2  v/3 2 1 2  l 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  (2-  1  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: s (t) a  +js {t) p  [sd(t)+js (t)}. q  (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 electromechanical 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 steadystate 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,  16  2.1. Evaluation of Existing Variable Representations  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> —uj , r  0  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 quasisteady states, the oscillating motion of this vector in the complex plane is slowed down by system frequency through multiplication by a factor e~ . The process of slowing the julot  oscillations is also referred to as relaxation in this thesis. Interestingly, the method by U . Linnert to combine electromagnetic and electromechanical 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 frequency component.  17  2 . 2 . Complex Signals and Dynamic Phasors  Based on the above observations, research on a simulation w i t h 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. I n order to slow down the oscillatory motion of such variables, one must be able t o , at least approximately, extract information about the frequency or spectral content from this variable. Because of this, the variable must be translated into an adequate representation t h a t also allows the subtraction of a part of its spectral content.  T h e latter could be easily achieved by a complex  representation through multiplication by e~^ , yet the construction of a complex Uot  variable from a real voltage or a real current is rather challenging.  • Integration of differential equations:  T h e question is whether the variables i n  their new representation can still be integrated. Unfortunately, this is not always the case. Secondly, i t is of great interest to find a numerical method t h a t is most suited for integration because the gain i n 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.  B o t h areas are closely related.  I n the following section, the author introduces a feasible  and flexible representation for system describing variables.  A n analysis concerning their  numerical integration is presented i n Chapter 3.  2.2  C o m p l e x Signals a n d D y n a m i c Phasors  This section expands on the idea of using complex numbers for system describing variables. While such complex representations have already been used i n stability simulations to simulate A C power systems w i t h large step sizes (cf. previous sections), the following derivation of complex signals and dynamic phasors differs from these representations i n one important aspect: Instead of applying three-phase similarity transformations to obtain one complex number and a real zero sequence value (which is usually neglected i n stability simulations), the dynamic phasor approach expresses each phase quantity as a complex signal. This way,  18  2.2. Complex Signals and Dynamic Phasors  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 - " *, thus relaxing their oscillating motion. -  7  0  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  Let us consider the steady-state solution of the R L circuit shown in Fig. 2.1, as often needed for the initial conditions of a transient simulation. For the purpose of this example, a single phase representation is  e(t)W  assumed, albeit a notation in matrix form easily expands the derivations to the poly-phase case.  The  transient behavior of the circuit can be described by Figure 2.1: RL-circuit  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) = v cos [tJot\. The system is thus described by v(t), i(t), and e(t) which are all a  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 V(LV ) = [ R + j cv L] 0  where V ( t J o ) and I(u ) 0  0  (2.H)  I(u ) 0  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 axis . In general, denoting phasors as S(cv ), we have 1  0  s{t ) = M[S(LU )}, 0  0  (2.12)  where t = t is the time of the beginning of the transient simulation. Conventionally, this 0  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> to them [56, 78, 79]. In this monograph, we assume the former 0  definition and characterize phasors as the time-invariant frequency transformation of a sinusoidally 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  Introduction  to C o m p l e x Signals and D y n a m i c  Phasors  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 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 threephase 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. 1  20  2.2. Complex Signals and Dynamic Phasors 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 quantities 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 co , with its half0  bandwidths being generally less than ui , as shown in Fig. 2.2.  Such signals are called  0  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  -co  n  0  |s(co)|  C0  n  Figure 2.2: Continuous amplitude spectrum of a bandpass signal.  CO  21  2.2. Complex Signals and Dynamic Phasors  Aco  -co  |s(co)|  CO  f  Figure 2.3: Discrete amplitude spectrum of a bandpass signal.  s(t) = lim  Au;->0  a,i cos \(co + iAco) t +  (2.13)  0  — ' i=—oo  Applying the trigonometric identity cos A + B = cos A cos B — sin A sin B  (2.14)  the above signal can be rewritten as lim  s(t)  £  cii cos [iAco t + (f)i\ COSCO t 0  (2.15) Thus, we have s(t) — S[(t) cosLu t — SQ(t) sinw ^, 0  0  (2.16)  where sj(t)  =  lim  Aw->0  >  ^ — ' i=—oo  di cos [iAco t + 4>i  (2.17)  a, sin [z Au; £ + </>;  (2.18)  oo  <?C0  S  =  ^  m  Au—»0  '  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'mco t. Both carrier signals are of the 0  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  22  2.2. Complex Signals and Dynamic Phasors  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(co)|  I -co  CO  r  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)] = (t)+js (t). Sl  (2.19)  Q  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)]e .  (2.20)  julot  Substituting (2.19) into this equation yields =  [si(t) + j s (t)) [cos uj t + j sin co t]  =  si(t) cosu>ot — SQ(t) sinu t + j [si(t) sin u) t + s<g(£) cosco t]  =  s(t) + j H[s(t)]  Q  0  0  0  0  0  (2.21)  where H[s(t)] =  S(T)  dr  TX .1 t — T oo  (2.22)  23  2.2. Complex Signals and Dynamic Phasors  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 counterpart 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. steadystate 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. A t steady state, all signals are sinusoidal and have only one frequency component, namely system frequency cu . In this case, equations 0  (2.17) and (2.18) for the in-phase and quadrature components reduce to j(t = t )  S  =  0  s (t = t ) Q  a cos</>o  (2.23)  0  (2.24)  = a sin</>  0  0  0  Consequently, the corresponding dynamic phasor is s/(t ) +js (t )  ®[S(t = to)]  0  Q  0  =  a e-  (2.25)  0  which can be readily determined from the steady-state phasor: S[S{u )] (po = arctan  a = | S(u) ) 0  0  0  &[SM]  (2.26)  The initial complex signal is then obtained with S(t = t ) 0  @[S(t )]e J ^ 0 0  to  (2.27)  24  2.2. Complex Signals and Dynamic Phasors  It should be noted that usually t = 0, in which case the initialization of dynamic phasors 0  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 function of the aforementioned example, we have e(t) — v  cos [cu t]. The Hilbert transformation  a  0  gives: U[e(tj\ = v  a  sin [oo t]  (2.28)  0  so that the new, complex forcing function is given by E{t) =  e(t)+jn[e(t)]=v e jul t a  0  (2.29)  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 D C 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 (2.30) (2.31)  V(t) = RI(t)  +  L-I(t)  (2.32)  25  2.2. Complex Signals and Dynamic Phasors  The orthogonality also allows the same mathematical operations, e.g. summation, multiplication, 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 OJ L^J ®[ lit) } + Lj $)[ 0  I(t)]  f  (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  Instantaneous  Values  v(t)  =  Ri(t)  +  Lj i(t) t  Complex  Phasors  V(cu ) 0  =  (R + jco L) I(co )  Complex  Signals  V(t)  =  RI(t)  Dynamic  Phasors  &[V{t)]  =  (R +  0  +  0  LJ I(t) t  JLu L)®{I(t)}+Lj @[I(t)}  steady-state part  0  t  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  Consequently, during steady-state initialization,  to suppress parts of its oscillatory dynamics. dynamic phasors are identical with their  corresponding phasor. With the above concepts and definitions, we have introduced useful tools for the simulation 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.861.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 techniques currently implemented in power system simulation programs, particularly if the sim-  27  2.3. Spectral Properties of Complex Signals  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  S p e c t r a l P r o p e r t i e s 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  f s(t) e~ dt  s{ui) = — L V27T  (2.34)  ju}t  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)| . It can be expressed either in the frequency domain or in the time domain. The 2  following two relationships define average frequency [15]: oo  (u)=  oo  J u\s(u)\ duj = -j 2  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  28  2.3. Spectral Properties of Complex Signals  density spectrum |s(o>)| is always symmetric about the origin. Because of this symmetry, the 2  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)\ dco 2  (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 oo  (co)  / { - w§\^ m  jl  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)\ dt 2  (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 Heisenberg'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 identical 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 particular 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 instantaneous 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].  30  2.3. Spectral Properties of Complex Signals  Definition 2.3 Instantaneous frequency is the time-derivative of the phase angle of a complex 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 detected by observing the magnitude and instantaneous frequency of a complex signal. After discretization, if \S(t)\ « const,  and Ui(t) « u  Q  (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> while the instantaneous frequency is zero. Physically, this 0  could represent a 60 Hz component on a D C 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 e ^ jv  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.  31  2.3. Spectral Properties of Complex Signals  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 0  1  1.  1  1  1  1  1  1  1  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  0.045  1 0.05  -600' 0  ' 0.005  1  0.01  1  0.015  1  1  0.02  time t [s]  0.025  1  0.03  1  1  0.035  '  1  0.04  0.045  0.05  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 s(t)  =  S\  frequencies  cos [LOI t] + s cos [OJ t] 2  2  so that the corresponding complex signal is: S(t)  =  S (t) +  =  Sie  = where the two amplitudes are positive.  S\ and s  2  S (t)  1  i W l 4  2  S2e ^ j  +  1  \S(t)\e^  {t)  are constants and u>\ = 377 rad/s and u  2  The signal could represent a 60 Hz signal with a superimposed  = 1131  third  rad/s  harmonic.  Its spectrum is thus composed of four peaks (Dirac pulses) at uj\ = ± 3 7 7 rad/s and uo — 2  ± 1131 rad/s.  Solving for the phase angle and  ip(t) =  amplitude,  S\ sin [UJI t] + s sin [<x>1]  arctan . . L S i cos [LOI t\ + s cos [u>2 t\ 1 2  2  2  \S(t)\  =  \Js\ + s + 2 S i s cos [0J2 t — cuit], 2  2  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 instantaneous 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 = 1.0 p.u. The instantaneous frequency is continuous and ranges outside the bandwidth. 2  In (b) Si = 1.0 p.u. and s = 0.5 p.u. Although the negative spectrum is suppressed, the 2  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 disadvantages 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. A n extensive list of first to sixth-order algorithms including the algorithm "families" by Runge-Kutta, AdamsBashforth, 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  34  3.1. Integration with Complex Signals  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)  RI(t)  -L~  L  + 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:  + AtF{t).  I(t)=I(t-At)  (3.4)  Backsubstitution of the above function gives the complex current:  I(t) = (R+ B  where  E  H(t-At)=^R+-^j  1  V(t) + H(t-At),  (3.5)  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  35  3.1. Integration with Complex Signals  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) = J(t-At) + ^ T  (F(i) + F(*-At)),  R  (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  and  H(t-M)  2L\ = ( R+ — j  1  V(t) + H(t-At),  (3.8)  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 trapezoidal 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  36  3.1. Integration with Complex Signals  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 j ®[I(t)}  = {-L- R-ju l)  ^{I^  l  t  Jo  + L- &{V{t)} 1  (3.10)  Here, A = — L~ R — j cu l is the system matrix or, in this case, the RL-series branch matrix. x  0  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  if (t-At) = e ^ °  A t  - -Tj-AJ  r  ~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. A n example for this approximation is shown in Fig. 3.1(a) for the integration of the signal s(t) = 1 p.u. sin [u 1] between the times ti — At and t\. 0  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  37  3.1. Integration with Complex Signals  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 co the relaxed integration 0  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 + JUJ L)-  1  0  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)  38  3.1. Integration with Complex Signals we obtain the equation yJ{t)  (3.15)  = exp[-At]L- @[V(t)]. 1  After writing the integral form, backsubstitution yields  y t  @[I(t)} = exp[At]  exp [-AT]  L '  1  ®[V{T))  (3.16)  dr.  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-At  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- 9[V{t)}-exp[AAt]  L ~ ®{V{t-At)^  l  I  L  + ...  t  ...  +A-  1  exp[-Ar]  (3.18)  dr  L~ ®[V{T)} X  t-At  where this last integral can be approximated as t  J  t-At  exp[-Ar] L- 9\V{r)\  dr  x  / exp[-Ar]L  - 1  ^———i  i i dr  (3.19)  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 exp[—At] = exp[—At] A . It is obtained with the Pade approximation also imple_1  - 1  mented as part of the  MATLAB  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 ( £ ) = - ( l + ^(l-exp[AA<])) A" L ~ V(t) + H(t-At), R  v  1  L  (3.21)  39  3.2. Integration Errors where the history values are determined by: H(t-At) ...  = exp[AAt] e °  j u AT  I {t - At) +  ...  + (exp[AA*] + ^ - ( l - e x p [ A A t ] ) ) e ° ju}  Ai  This equation was first developed in [39] for lossless elements.  A~ L~ l  Y  V(t-At)  (3.22)  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 E r r o r s  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 question 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 convolution rules of integration are analyzed with respect to their numerical error. While the backward 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  40  3.2. Integration Errors  single-phase inductive branch. For a combined series resistance and reactance, the currentvoltage relationship with the relaxed trapezoidal rule is given by (3.11). Setting R = 0, this equation becomes HTR  At  L  -i  +  {t-At),  (3.23)  H  with the history values being if (t-At) = f * . L - V(t-At) + 2 4- j w A t 2+j e  A ^°  1  0  LOQAI  A t  e  /(t-At).  (3.24)  Applying the Z-transformation's displacement property Z[f(t — kAt)] = z~ Z[f(t)]  for  k  discrete time series /(£) to the above equations, yields the relaxed trapezoidal integration algorithm in the Z-domain:  J  W = 2T7^ " R  L  2 + jcu AtzL  w  y  0  ;  -  I n z ) +  ffW  2 + joj At  z  0  (3  y  j  y  "  25)  '  In the Z-domain, we can now combine these two equations and construct the admittance transfer function Y(z) = I(z)/V(z) y,  x l Z J  RTR  ~  for this integrator:  z -t- e L (2 + j u ; o A t ) ^ - ( 2 - i w A t ) e ^ o A t 0  . . ^- 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  EMTP.  Assuming a sinusoidal excitation V(t) — e^ , the response in the time-domain is given wt  by /(t) = Y(z) e  jujt  ,  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  41  3.2. Integration Errors  Table 3.1: Transfer admittance functions for an ideal integrator with the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algorithms. BE  TR  RTR  Y(~) K  At  z + pjuo&t  L  (2 + j LO At) z - (2 - j UJ At) ei ">  -  ^  J  RCV ]  LOQL  0  \  0  j U! At (z — 1) Q  [u  At  /  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  Example  3.1  LO/LOQ.  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 10 accurate.  Integration Errors  3.2.  Y(z)/Y(s)  T h e results o f Fig. 3.3.  4 2  T h e plots show the magnitudes  of the step-size-to-period c y c l e , i.e. ( 2 . 1 ) .  for the four i n t e g r a t i o n a l g o r i t h m s are d e p i c t e d i n F i g .  At  =  Nyquist  ratio  T / 2 ,  At/T  a n d phases of t h e y ( z ) / Y ( s ) - r a t i o s  At i  w h i c h is b o u n d between  m n  =  3 . 2  and  as a f u n c t i o n  0 a n d t w o steps p e r  the m a x i m u m possible step size a c c o r d i n g t o the N y q u i s t c r i t e r i o n  Notice that t h e graphs have n o t been displayed beyond t h e location of their  first  pole.  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 ) . Figure 3 . 2 shows t h e numerical accuracy of the relaxed trapezoidal a l g o r i t h m i n comparison w i t h the accuracy  of t h e b a c k w a r d E u l e r a n dtrapezoidal rules o f integration.  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 volution rule w i t h b o t h the backward E u l e r a n d trapezoidal  •  the relaxed  con-  algorithms.  N e x t , let us carefully discuss a n d interpret t h e depicted graphs. qualitative observations  Several  I n general, t h e following  c a n b e m a d e f r o m t h e curves:  I n terms of magnitude, t h e backward E u l e r integration a l g o r i t h m shows a higher  accu-  racy than the trapezoidal algorithm. However, due to the increasing mismatch i n the p h a s e angle, t h e r u l e is 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 sizes.  •  T h e t r a p e z o i d a l rule also shows increases.  These,  however,  deviations  f r o m t h e exact s o l u t i o n as t h e step  are confined to the magnitude.  Overall, this rule  size  permits  larger step sizes 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 results obtained w i t h t h e relaxed trapezoidal integration rule require a more detailed discussion:  A t steady  state,  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  sizes.  F o r co > u ; , t h e m i s m a t c h i s l e s s 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 0  w i t h higher frequencies the  trapezoidal  the frequency  u n t i l , a t co—>oo,  algorithm.  decreases.  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 seen  F o r co < co , t h e d e v i a t i o n s  A t co = | w  0  0  , t h e m i s m a t c h is already  trapezoidal or the backward Euler algorithms. occurs i n t h e magnitude, n o t i n t h e phase  •  T h e relaxed  become  more  with  s i g n i f i c a n t as  larger t h a n w i t h the  F o r a l l frequencies,  the mismatch only  angle.  c o n v o l u t i o n a l g o r i t h m also produces  no mismatch at system  frequency.  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 < co , the accuracy is generally better than 0  with the relaxed trapezoidal rule, and for | co < to < co it is even better than all the 0  0  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 = co . To 0  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 At -T  1  a; =  V\ 1 \to  . 1'  —  3  '  2  9  )  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 co . 0  Observing the previous figures, a quantitative statement about the largest possible integration 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  46  3.2. Integration Errors  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, e =  \Y(z) Y(s)  Again, assuming a sinusoidal excitation z = e  (3.30) J t j A i  , the truncation error for the relaxed  trapezoidal rule becomes w At w At + 2 tan [^(w — w 0  (3.31)  1 0  Table 3.2: Truncation errors for the backward Euler (BE), trapezoidal (TR), relaxed trapezoidal (RTR) and relaxed convolution (RCV) algorithms. BE  e  TR  < =  RTR  e =  (uAt\  exp[j(^)]  V 2 J sin[*f*] (uAt\ cos[^] V 2 J sin[^] to  At  o; A t + 2 tan [^(w-w)] 0  RCV  e =  1  0  w- w  uJ  0  sin sin [f (w OJ ("ft) sin[ f ] w  w )] 0  w  0  0  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  47  3.2. Integration Errors "[...] 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 + JLU L)-  V(t)  1  0  (3.32)  for A t —> oo (see Section 3.1.4), it can be verified that the rule does converge to this steadystate 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 frequencies so that the maximum step width is determined according to (2.2): At  7T  max  (u;j) = - —  (3.33)  With a tolerated integration error of e = 3%, the resolution factor k for the trapezoidal rule becomes approximately k R = 5, which is the recommended value for ordinary E M T P T  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 instantaneous frequency can be computed for the location where the graphs in Figures 3.4  48  3.2. Integration Errors  I  I  1  I  0  0.5  1  1.5  A i / T [p.u.]  Figure 3.4: Local integration errors for cu > uo . 0  o  0.5  1  A i / T [p.u.]  Figure 3.5: Local integration errors for u> < COQ.  1.5  3.2. Integration Errors  49  and 3.5 cross the tolerated error margin w  €RTR  A  = £• This location is given by "=  t  cu At + 2 tan [ ^ ( w - w ) ] 0  £  (3.34)  l  0  which must be solved for A t . Assuming £ < 1, we can make a cubic approximation of the tangent function for small arguments, i.e. 3  1  tan x ~ x + - x . With this approximation, the formula for the step size calculation becomes: I  A+ A  t  -  M  \  )  /  1  2  £  U  J  i  i  = \Ja^--of  {  c  ^  =  l - £  for  I 1 + E for  Wi>W  0  < c  0  U i  fo  ( 3  oc\  '  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 to t 4- A t , the global 0  0  error represents an accumulation of local errors by integrating the system for a number of steps. Using the step size A t , 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=^-  uAt  K  (3.36) '  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 accumulated integration error even when used with extremely small step sizes. This clearly  3.2. Integration Errors  50  51  3.3. Numerical Stability  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 integration 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(t ) 0  = S  0  (3.37)  52  3.3. Numerical Stability with the solution cf>(S ,t) = S e  (3.38)  xt  0  0  where S, So, and A may be real valued or complex. Complex values are allowed since higher-order differential equations can have complex eigenvalues. A n 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 m E = (1 — A A t ) , and thus - 1  B  the backward Euler method is stable whenever 1  >  <  |(1 - A A t ) | 1 ( 1 - A At) | |(»[AAt]-l)+j(3[AA*])|  <  m\At]  <  _ 1  ^  - l ) + (3[AAt]) . 2  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  E  IS.:'/ '  i  liiiiiiii 1  i  1  lillilSSl!PI18lllliSlHllHII!|  12 Re[XAt]  Figure 3.8: Stability region for the backward Euler integration rule (shaded area).  3.3.2  Trapezoidal  1  2  Re[XAt]  Figure 3.9: Stability region for the trapezoidal, relaxed trapezoidal and relaxed convolution rules of integration (shaded area).  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 m  TR  = (2 + AAt)/(2 — A At), and its fixed point is stable if  |2 + AAt| |(2 + 5R[AAt]) + j(3f[AAt])| (2 + K[AAt]) + (3[AAt]) (2 + $R[AAt]) 3?[AAt] 2  2  < |2-AAt| < |(2-tt[AAt])+j(3[AAt])| < (2 - K[AAt]) 4- (S[AAt]) < (2-ft[AAij) < -5?[AAt]. 2  2  (3.42)  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.  54  3.3. Numerical Stability  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: S(t) =  2 + XAt • 2 - XAt  (3.43)  Hence, the method is stable if 2 + XAt  J w At  < 1.  2 - XAt Because of | e  jtjAt  (3.44)  | = 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) = e S(t-  At).  XAt  Hence, the characteristic multiplier is m v RC  (3.45)  = e , and the integration method is stable XAt  when 1 > 0  |e  AAt  |  > 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. A n 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 statespace 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  56  4.1. Step Size Prediction  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 representation 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 Vi(S, t)  angle[S(t)] - angle[S(t - At)] At  (4.1)  is sufficiently accurate. Comparing the instantaneous frequencies of each variable, the approximate maximum frequency u> (t) > to can be obtained at each instant. max  0  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.  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). = U + OC t a  (4.2)  a  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\io (t ) max  (4.3)  l  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, to  ,  predictedtmax  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. u in(t) m  < 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 to . Applying the above Q  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  58  4.1. Step Size Prediction  the frequency information is not availabe (e.g. in D C systems). The magnitude deviation is approximated by: AAi(S, t) = arctan  \S(t)\-\S(t-At)\ A t S„  (4.5)  where S is the nominal value of the corresponding variable. The arcus tangent function is n  used to obtain a normalized magnitude deviation for all variables, which ranges from zero to 7r/2. Similarly as with the above frequency prediction, the highest magnitude deviations A-A  of each time step are sampled over 3/4 of a period. Instead of applying the linear  max  regression method to forecast the highest deviation, the maximum sample is taken as this value. Hence, A A predicted = max [ AA  }  max  (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 < ^•magnitude  2 AAp di ted re  c  TX  , . _\ V*'' I  The function becomes ^magnitude — 1-0 if there are no deviations in magnitudes of all variables, and it is zero for AA  predicted  4.1.3  - > 7r/2. It is shown in Fig. 4.1.  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 oscillations [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 UJ  0  •  J 0  TT/4 ^ p r e d i c t e d  TT/2  [rad]  0  50 SEC \aa| [% of W 0 / ]  100  0  1 Upredicted/Uo  •  •  l_  1 2  [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 stability mode. Otherwise, the program should use the electromagnetic transients mode." The rule consists of three parameters: the magnitude deviation AAp  ,  redicted  the gradient of  the frequency samples, i.e. absolute slope \a \ of the regression line, and the ratio of predicted a  to nominal frequency. To build this ratio, let us first define ^predicted — max \_CO predicted,max j 2 CJQ The ratio is then given by copredicted (t)/wo-  LOpredicted,min] •  ( 4 - 8  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 [i are the two membership values for the "low gradient" and 2  "frequency close to nominal" fuzzy sets, respectively. The membership value for the electromagnetic 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 electromagnetic 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.33)  or  (3.35).  However, to achieve an optimum between accuracy and  simulation speed, only half of this step size is used in the simulation, i.e. g  opt  At  TR  =  — Qopt Cpredicted, max  R , i  =  -  A , & t  R  0.61 T  I"  —  A i A-tRTR.,2  (4.11)  {-"Jpredicted,max  / . ., \ (4.12) 0  . /-  Qopt  = 2. Hence,  -g  V ^predicted,max  ~~ ^ 0 |  0.59  / LOpredicted,min . ' — — T3Qopt V \tOpredicted ,min to \  =  /. \ (4.16) 1  0  0  where, during the stability mode, the lesser ofthe two relaxed trapezoidal step sizes is used. Thus,  At  stab  ; At  = min [At  ].  (4.14)  RTR>2  RTRtl  For the electromagnetic transients mode, the magnitude deviation is applied as an additional constraint for the trapezoidal step size. Therefore, the electromagnetic transient step size becomes: At  emtp  = min [At  TR  J  (1 -  ^magnitude)}-  (4.15)  U>o  By weighing the electromagnetics time step At  emtp  their respective membership degrees  ^ p emt  and the stability time step At  stab  with  and (J, t b, we obtain the recommended time step: s  a  Atpredicted — f^emtp Atemtp "f f-^stab Atstab  (4.16)  61  4.2. Prediction Reliability  At each calculated time step, this formula provides the algorithm with a new step size forecast. At i Tn  n  4.2  Notice that the final step size is always limited by the user-defined boundaries  and /\t  .  max  Prediction Reliability  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 electromechanical 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,  de  mo  for the three criteria, respec-  tively, the membership for high reliability can be obtained using scalable monotonic chaining [17]. This approach to fuzzy reasoning has the advantages of being faster than conventional 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 overreacting 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  62  4.2. Prediction Reliability  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 A t . 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: H =  NoF At  4  8  ; 0< ^4 < 1  (4-18)  7T  The two membership functions are combined utilizing the geometric average to give the total membership function  4.2.2  ii ie samp  = yjv-z  f'4  f ° the fuzzy set: "The sampling is good." r  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  to (ti). max  The fluctuation is measured by  computing the standard deviation of the samples, thus  a= \  X  N  o  2  F  J^^( max{U)  \|  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, fJ-fluct = 1 - f-z Kj.L  a  ; 0 < Hftuct < 1  (4.20)  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 \i mtp  a  e  n  d Ustab, considering that one of them  becomes 1.0 whenever a phenomenon can be clearly categorized, e.g. if p b is either close to sta  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 < \l ode < 1 m  (4-21)  64  4.3. Step Change Efficiency  The membership function for the fuzzy set "The program mode is well defined" is shown in Fig. 4.3. Frequency  0  fluctuation  50 Std. variation a [% of  is low  Program  100  0  fpredicted]  mode  is well  0.5 Membership degree  defined  1 fi  s t a b  [p-u.]  Figure 4.3: Two membership functions for prediction reliability.  4.3  Step C h a n g e 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 /\t diced/At. pre  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  T /T i . rebuM  so ution  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  \x  siev  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  65  4.3. Step Change Efficiency  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 coefficient (^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 r = At s  /At  prediced  is large. More specifically, for ratios  r greater than 1.0, the efficiency was defined to grow exponentially. Thus, s  1.2  Pstep —  for  10  r > 1, s  (4.22)  until the membership degree is 1.0 at r > 10. s  For step ratios r slightly below 1.0, the stepsize need not be immediately changed as the s  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  instep  Qmin  for  r <l, s  (4.23)  Qopt  where the margin is defined through g  min  step ratio is high" is depicted in Fig. 4.5.  = 10. The membership function for fuzzy set "The  66  4.3. Step Change Efficiency  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 r (, jw, and the time needed to re  complete the solution of one time step r i so  .  U  The program automatically measures both  ution  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 T n Tebu  d  is negligible compared to  nT i ? so ution  If we choose a small tolerance £ = 0.005, then we can assess n =  (4.24)  T r e b u M  £ ^solution  This value represents the number of steps that must be computed with the old step size A t until the time for rebuilding the matrix appears negligible. In addition, we must assume that the time delay n A t should always be constant, independent of what step size is used. Therefore, n A t = nAt i  pred Cted  = 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  For r  s  delay  = ^  for  T  s  > 1  (4.26)  < 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  IO"  10°  1  Time delay is high  10  1  0  Step ratio r  s  6,000  12,000  Number of steps K  Figure 4.5: Membership functions for efficiency evaluation.  4.4  D e c i s i o n M a k i n g Process  So far, we have explained the techniques used to forecast an optimum step size, to evaluate 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 accomplished 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-  68  4.4. Decision Making Process  Defuzzification hedge' step change certainty  prediction reliability  step change efficiency  [ < 0.6V V y  yes |  old At  \  no  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  (4.27)  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  (i hange C  —  0.6, the program accepts  At dicted pre  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].  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  d di  X{t)=g(X,U,t)  ,  (5.1)  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 I(t) =  YV(t),  (5.2)  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 69  70  5.1. Overall System Equations  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. However, 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 components, 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) X(t)  =  Y V(t) + H(t - At)  (5.4)  = 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)  71  5.1. Overall System Equations  Electrical synchronous machine equations: A C network equations:  State-space equations, discretized with the Trapezoidal Rule  Nodal admittance matrix, discretized with the Relaxed Trapezoidal Rule  K-o).  (co = 377 n  rad/s).  Mechanical synchronous machine equations: State-space equations, discretized with the Trapezoidal Rule  K  =  o).  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  '  2  ^X(t) + X{t-At)e^p[n}  ),  (5.7)  where Q = j At diag[cj ] is a diagonal matrix whose elements consist of the product of j At n  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 [ n ] x ( t - A t ) ... P  ... - ~  (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  EMTP  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]. W i t h larger simulation time  72  5.2. Program Structure  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]: " v(t) X(t)  \  i  r  -  k+1  T ^ / J A I • v(t) X(t) k  r  1  A TTVJA  ' AV(t) AX(t)  I  "  (5.12)  The improvements are obtained from F(V,X,t)  dF dV  dF aX  G(V,X,t)  dG dV  aG dX  ' AV{t)  '  (5.13) -k  AX(t)  Convergence is quickly achieved with this method, provided the starting values are sufficiently 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 techniques to skip the updating from time to time to obtain a faster solution, e.g. with the "dishonest Newton method" [28].  5.2  Program 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 d a t a  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 At  min  and At ,  respectively, as well  max  as the end time of the simulation t d and the aggressiveness factor A (see Section 4.4). en  P o w e r flow a n d steady-state i n i t i a l i z a t i o n .  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  T buM Te  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  (  Stop  Output data  )  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  Time-Domain Simulation  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 t  . In  end  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 prediction of the voltage and state variable vectors. To allow fast convergence, great care must  74  5.3. Time-Domain Simulation  <  Start here  >  Predict node voltages and state variables for /.  Compute forcing functions at t.  Compute model history values and collect to nodal history vector.  Rebuild the nodal admittance matrix with new topology.  Solve the system at time ffor unknown node voltages and current injections.  Interpolate for the time t of this violation, and set a discontinuity flag for this instant. x  For each discontinuous model, check if limiting values have been violated.  Determine a new step size A t .  For each model store computed voltages and currents for calculation of new history values.  Rebuild the nodal admittance matrix with Af ... .  Set t = t, Af = Af n  O  '  and t=t + At  new  Stop  0  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 presented 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\S (t + At )) = P  9\S(t)]-z*E9 {t-  p  At)] ,  [S  (5.14)  where @[S (t + At )] is the predicted variable for the predicted new time step At . Using p  p  p  (2.20) to transform to complex signals, we have S (t + At ) = p  p  A  t  ^  A  t  p  e^  A  t  " S(t) - ^  e  WAt A +  t f ) )5  (  t  _ y At  ( 5  1  5  )  Here, u> is the nominal frequency of the variable that is being predicted. For instance, for n  voltages in an A C transmission system u = uo, while state variables associated with the n  machine models are D C quantities so that co = 0 there. n  At moments of discontinuity, the voltages of the predominantly inductive transmission network may change instantaneously whereas the fluxes cannot change at an instant. Therefore, 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 nonlinear 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  5.3.3  X(t).  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 t , where (ti — At) < t x  x  < t\.  In case of a limit violation, t i is x>  5.3. Time-Domain Simulation  77  determined by linear interpolation. This may not be very accurate, particularly if large time steps are used. Suppose t  X)1  is a not so accurate first guess for the moment of discontinuity,  and also suppose that (ti — At) < t i < t X}  x  < t. x  In this case, the program computes the  solution at t = t \ without introducing the disturbance. x>  is again used to locate t . x  Since the interval  \{t\ — At), ty], the next interpolation at t  X]2  [t ,i,^i] x  Afterward, linear interpolation  is smaller than the original interval  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 = t of the disturbance. x  Then, the simulation algorithm uses a backward Euler step with the user-defined minimum step size At . min  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 of the N o d a l A d m i t t a n c e  Matrix  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 contribute 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 a  D  rn, m. rn  c  a  b  c  K  m  b  [Y]  -[Y]  -fYl  fYl  rn  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 connection 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.  79  Time-Domain Simulation  even if several n u l l - i m p e d a n c e wires are c o n n e c t e d i n p a r a l l e l or a n y o t h e r c o n f i g u r a t i o n . practice, however,  s u c h c o m b i n a t i o n s are r a t h e r u n u s u a l .  Node reduction and sanity checks. ground node(s)  In  The  nodal admittance  m a t r i x still contains  a n d is, therefore, singular. A p r o c e d u r e was i m p l e m e n t e d t h a t  the  subsequently  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 . If 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  value.  The  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 are  e l i m i n a t e d a n d voltages are expressed w i t h respect to the chosen g r o u n d .  A t the same  then time,  t h e p r o g r a m p e r f o r m s s a n i t y checks for a n y p a r a l l e l v o l t a g e sources or series c o n n e c t i o n s current  sources.  of  Part II Component M o d e l i n g  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  82  6.1. Derivation ofthe Model Equations  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  D e r i v a t i o n of the M o d e l E q u a t i o n s  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 d dx  2  d ^v(t,x) dt 2  v(t,x)-L'C ' '  2  d —  v  2  2  d — i(t,x)  = 0  (6.1)  2  i(t,x)-C'L'  = 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  83  6.1. Derivation of the Model Equations  solution of the above differential equations, they must first be diagonalized. Using eigenvalue decomposition, an orthogonal transformation T~ =  can be found such that  l  A = T- L'C'T  =  l  v  T7 C'L'Ti x  is a diagonal matrix, and v = T~ v and i = T7 i are the modal values of line voltages and l  1  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 d —  d —  2  2  v(t,x)-\  v(t,x)  0 (6.2)  d  2  -  d  -  2  The French mathematician Jean Le Rond d'Alembert (1717-1783) first described a solution [7] to the equation d —  2  d — u(t,x) = 0, 2  u{t,x) - c  2  where c = +y/X is the propagation constant, at time * i as .  .  u{t x)  lf(ti + cx) + ip(ti - cx) =^  u  tl+CX  1  .  ,. . ,  2c~ I ^  2  t\  d  a  —cx  (6.3)  -  The functions ip and ip are given through the boundary conditions  m=  (p(t) = u(t,x) 2  d u(t, x) dx  =0  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) = v (t) and i(t,0) = i {t), while the quantities at the k  k  other end © at x = £ are v(t,£) = v (t) and i(t,£) = -1 (t). m  Then,  m  <p(t)  v {t)  m  -L  k  s  ! t  (0  where L is the mode inductance found from the diagonal matrix [86] T~ L'Ti. l  T7 C'T l  v  Similarly,  is also diagonal with its elements C being the mode capacitance. Substituting the  above functions into (6.3) yields _  _,.  \  ^k(ir + cx) +v (ti  w(£i,x) = —  '—  k  -  -cx)  L  tl+CX  2c  f  I  tl—cx  -.  .  ,  ^  6  )  84  6.2. Lossless Line Model Letting x — I, we obtain the voltage at terminal @ at time t\ v (ti) m  Vkih+r)  =  L  +v (t -r) k  l  «k(*i+r) -  ik(*i-T)  2v/A  (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) =  -  c + 2x/A  g  Vk(h  + T) -  v (ti-r) k  (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  G 't)m(ti) - i ( t i ) c  G  m  c  'Wk(*i) -  *k(*i)  =  G w (t!-r) +  =  G w (ti-r) + i ( t i - r )  c  c  k  m  i  k  (ii-T)  (6.7)  m  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 transmisision line wave equations.  6.2  Lossless L i n e 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 R (purely resistive for lossless lines). For easier notation we will also make use c  of its reciprocal value, the characteristic conductance G . More details on modal transforc  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  85  6.2. Lossless Line Model  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  ®a  i (t) k  ijt)  Ut)  77)\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  86  6.2. Lossless Line Model terminals (K) and © are given by  Mt) i (t)  v {t)  G  k  r  v (t)  Gr  m  m  +  h (t) h (t)  (6.8)  k  m  where the history vector is known and can be calculated from previous current and voltage values as follows: r  " h (t) 'i _ h (t) _  n i rv - (t-r) u _\  G' c  k  i  _ v {t-r) _  m  r  ~ u _\ i 1i i" r Mt~ ) r  k  m  1  (6.9)  . i {t-r) _ m  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 At <  T„  and  for i = 1... n  mm  (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 r , yet it precludes this model from being useful for stability studies. min  Considering the storage requirements for this model, the values of h (t) and h (t) have k  m  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 At  T  $>[S(t) "  L  -  W  J  9[S{t- At)) , A  (6.11)  F  Substitution with (2.20) yields the interpolation formula for complex signals At - r At "  Si(t - r )  S(t)  ~  w  T  ,ju> (At-r) n  ' At  S(t - At) ,  (6.12)  where uj is the nominal frequency ofthe transmission line. Thus, introducing the two factors n  V  At At  i  -3 u r  ju! {At-r)  (6.13)  n  n  ' •* 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  G 1-p  1+p -2p  2  c  . im(t) _  2  -2p 1+p  v (t) v (t) k  2  m  +  hm(t)  (6.14)  and a history vector whose values are determined from the past time step: k  qG  m  1-p  h {t) h (t)  c 2  V - i -1 p  v {t-At) v (t-At) k  m  + 1-p  2  p -1  -1 p  4(t-At) v(t-At)  (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 offdiagonal 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  Rc 2 p  r m v (m)  K(V  H  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 emphasized 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 stability 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 investigations. 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  89  6.3. Lossy Line Model  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  HR/4  ®  H  R/2  Ml  ®  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 A t 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° y ss  '  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 i and i k  i>k and v , m  m  which depend only on the terminal voltages  and on a history vector that can be calculated from past values. From Fig. 6.3  it follows immediately that i (t) i (t) k  m  = =  i^t) i (t) 4  = =  Tjii(t) TiU(t)  .  .  90  6.3. Loss}' Line Model  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 Y for the two diagonal elements and a S  mutual admittance Y for the two off-diagonal elements for one mode. For n modes, we M  obtain two diagonal matrices Y  s  and Y .  The transformation back to phase quantities  m  yields the matrix Y  Y  Ti Y  Y  Y  T Y  s  l  m  T" T~  TiYnT-  1  1  Ti Y  l  (6.17)  T~  l  s  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)  A  ' v (t)  B  k  '  v (t) m  -h (t)  (t)  x  V1  -h (t)  B  D  T  2  (6.18)  Mt)  -h (t)  v (t)  A  4  where the submatrices are given as  "GO 0 G  ' -G ,  0  Y  G + Y  s  Y  D  ti —  \G  0  0  0  0  0  - G  0  0  + Y  0  S  0  -W  0  0  G +  l 2  Y  Y  s  Y  G + Y  s  In order to simplify the notation, let the terminal voltages and currents be v (t) = t  k{t) i (t)  ; it{t)  v (t) m  m  and the internal voltages and history values Vi(t)  =  [v (t),v (t),v (t),v (t)}  T  l  2  3  4  hi(t) = [Mi),Mt),Mi),Mt)]  T  91  6.3. Lossy Line Model  Then, the current injections from the line model into the rest of the network through terminals ® and @ can be formulated as i {t) = [A - BD B ] l  v (t)  T  t  v  - BD^hiit)  t  v  '  Y, (At)  '  ine  v  h (t)  (6.19) '  ltne  For a simulation with the lossy transmission line model, its matrix Yn (At) ne  is inserted into  the nodal admittance matrix of the entire system, and the history values hn (t) are placed ne  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 v (t) 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~ B v (t) 1  T  t  - D^hiit),  (6.20)  which was derived from (6.18). Modal transformation yields the corresponding voltages v  l  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 . i (t) . ~ 2 3  G -G  -G G  ' v (t) " . v {t) _ 2  (6.21)  3  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 frequency 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 analysis [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 publications, 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" D C 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 turningdirection 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 N = v + 1 d-axis windings (v amortisseur windings and one field d  winding) and N = p g-axis windings on the rotor structure. This way, the models can easily q  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 linkage 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 e {t)  =  j i) {t)  -  e (t)  =  ^Mt)  ~  a  b  t  a  e {t) = j Mt) c  t  R ia(t)  (7.1)  aib{t)  (7.2)  R i {t)  (7.3)  a  -  R  a  c  =  ^fd{t)  +  Rfdifd{t)  0  =  ^ ^ ( t )  +  Rnd^Kdit)  0  =  ^TpKa{t) +  R qiKq{t)  fd{t)  e  K  (7.4) for K=1...U for  K  = 1... fi.  (7.5) (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  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) +  (7.7)  ^a,rotor(t)  where 1pa,stator{t)  =  -l ia{t)  ~ Lb ib(t) ~  aa  LcUit)  v  Ipa.,rotor  (t)  =  l  fd  lf {t)  YJ LKCI iKd{t)  +  d  H  +  K=l  YJ  I  A  K  C  kqit)  I  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)  L aO  LabO  ' ia(t) '  L bO  ib(t)  L bO  L 0  ic{t)  a  a  'Pcslatori't)  L a2  L bO  LabO L aO  a  —  ^b,stator(t)  L bO  a  x  a  a  a  aa  +  cos [20] cos [20 + T T / 3 ]  cos [29 + TT/3] cos [2(0 - 2vr/3)]  cos [29 - TT/3]  cos [29 - TT]  cos [28 - vr/3] cos [29 - T T ]  " ia(t) ' i (t)  (7-*  b  cos [2(0 + 2TT/3)]  and the corresponding linkages with the rotor windings are cos [9]  Lpa,rol,or it) "Pb,rotor it)  L fd ifd(t) + y~] L a  i d(t)  aKd  ^c,rotor(t)  cos [9 - 2TT/3]  X  K  C O S [0 + 2TT/3]  K=l  sin [0] j  ] L q kqit) aK  (7.9)  sin [0 - 2TT/3]  x  sin [0 + 2TT/3]  K= l  For the windings on the rotor structure, the flux linkage with the armature windings is determined by / fpfdelator (t)  =  -L  a  f  2ir ( i {t)  d  a  C O S  [9]  i (t)  +  b  C O S  / 1pKd,rtator(t) 'ip* ,stator it) q  = =  ~L  aKd  +L  aKq  I i (t) / a  li (t) a  COS[0]  +  l (t) b  C O S  [0 - — ] + i (t) 2ir  [0 -  c  — ]  2ir  sin [0] + i (t) sin [6-—] b  2ir \  +  l (t) c  + i (t) c  C O S  C O S  [0 + J 2ir \ —  ]  [0 + y ] j  sin [0 +  2TT \ —], J  (7.10) (7.11) (7.12)  7.2. Transition to Complex Variables  9 6  where the index isre= 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)  L  ffd  Lfid  •  1pld,rotor{t)  Lfu  Liu  •  i\d{t)  '4 i'd,rotor (t)  Lfud  L\vd  •  _ iud{t) _  )  L\lq  lplq,rotor (t) faq,rotor {t)  —  iftnq, rotor (t)  _  Li  M  L\2q  • ••  Li22q  • '•  L2  • •  M  (7.13)  Li  M  L  2  (7.14)  M  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 startingpoint 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 T (t) m  where  ipd(t)  and  ip (t) q  =  as well as  (Mt) *,(*) - Mt) »-(*)) , i (t) d  and  i (t) q  (7-15)  are obtained applying Park's transforma-  tion (2.5) to the armature flux linkage and currents, respectively. torque can be converted to its electrical equivalent,  T (t), e  The electromagnetic  through division by the number  of pole pairs n /2. p  7.2  T r a n s i t i o n t o 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  97  7.2. Transition to Complex Variables signals to describe a steady-state signal: In previous chapters, we wrote S(t) = S(to ) e  -  Jwo{i  0  (7.16)  to)  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) = S e  = SneW°  je  (7.17)  +UTt)  N  to represent the time-function of a constant quantity S , undergoing a uniform rotation in N  space with constant speed oj . Here, S r  could be a steady-state D C quantity such as the  N  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) = s e  + js e  ie  d  j e  g  .  (7.18)  The above equation still describes a uniform sinusoidal oscillation. More generally, s  d  and s can both be independent functions of time, thus q  S(t) = s (t)e^  0  d  + js (t)e^  (7.19)  q  Notice that S (t) = s (t) + j s (t) is the time-function of a complex rotor quantity with a n  d  g  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 s (t) or s (t) is always d  q  zero. For "complex" signals representing d-axis quantities, it follows that S (t) nd  "complex" signals associated with g-axis quantities are S (t) nq  = s (t), while d  = s (t). W i t h this definition, q  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.  98  7.2. Transition to Complex Variables  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)  L  j  i (t)  aKa  3  Kq  j0  (7.20)  J(0-2*/3)  X  e  j{0+2-n/3)  K=r  p  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)  +  ^ ] Lg g K  J2 ™ * ^ L  di  d  x  K=r  J  J q(t) K  K=1  X  /  -J 0  +  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)  (K[/ ]  = -L  a  afd  cos  [9] + %[I ] b  \  (Wa}cOs[9} + ^I } Os[9-~}  aKd  b  \ iPwtatorit)  = +L  aKq  [8 - ~]  +  U[I ] C  o  A.d, tator(t) = ~L S  COS  C  COS  (7.22) [8 +—  3  +  3  (7.23) U[I ]cOs[9+~3 c  U[I ] sin [8] + 5R[/ ] sin [9 - ~] + &[/ ] sin [9 + — (7.24) V o 3 a  6  c  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  99  7.2. Transition to Complex Variables  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 D C signal becomes S (t) = S(co = 0) , DC  (7.25)  0  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 D C source functions are real functions, it follows that the steady-state solution and thereby the D C 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 armature 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  Pd,stator{t)' L '4>q,stator(t)  ' i (t) '  ,l  d  =  d  (7.26)  —  io(t)  L  '$0,stator(t)  0  where + L o +- L 2  L  ao  d  L Q + L bO ~ ^ L 2  Lq L  0  aa  aa  =  a  aa  LaaO 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.  100  7.3. Model Equations with Complex Signals Introducing complex signals for the dqO components, we have  (7.27) (7.28) *o(*)  =  (7.29)  V>o(<)+j'3[*o(*)]  Because of the orthogonality of d and g-axis components as well as of the above real and imaginary parts, ip (t) and  are both alined along the d-axis, whereas ip (t) and  d  q  j $s[$ (t)] are located on the g-axis. Hence, ip (t) and jS[\I> (£)] share the same magnetic d  d  9  path through L , while ip (t) and j 3 [ f ^i)] both flow through L . The inductive path for d  q  q  the zero sequence remains unchanged. With this observation, we can now easily express the complex armature flux linkage in dqO coordinates:  ^d,stator{t) ' L ^q,stator (^)  ' Wd{t)\'  d  =  —  Lq  ^0,stator{t)  -j Lo  Wo(t)]  L  n  La  Wait)} 9[/o(t)]  (7.30)  This equation completes the transition from the real flux linkage equations to their complex equivalent. We have thus created a physically meaningful, imaginary counterpart to the electrical machine equations.  7.3  M o d e l E q u a t i o n s 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.  101  7.3. Model Equations with Complex Signals  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: L  Lafd  d  -L  i>vd  -  3 2  a  —L  —Lafd  — Laid  L fd  Lffd  Lfid  Lfvd  kd  Laid  Lfid  L\\d  L\ud  kd  Laud  Lfvd  Llvd  Lund  kd  d  3 2 3 2  ' Wd\  L ^d  Laid  a  ~L  2 3  ~2  T  Laftq  L lq  L 2q  " "  La/iq  Lllq  Ll2q  ' •"  Ll2q  L22q  Ll/iq  L iq  a  a  ^1 al  T  " "  a  ~Lq 3  -1 T  2 ^ajiq  aud  L 2q  L lq  q  '  a  2l  ' Wo) '  Li  M  kq  m  kq  • •  L  • •  Ljxiiq  -  2  (7.31)  (7.32)  and '  .  3[*o]  'K[/  -L  0  3[/  .  (7.33)  0  0  Notice that the d and g-axis equivalent circuits shown in Fig. 7.1 do not correctly represent 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 •  iDqK = j L/aDq  ,  and  ,  2  ip  DqK  L £)  =—  a  q  ip  <-> ^aKq  Kq  102  7.3. Model Equations with Complex Signals  in (7.32) for the variables of the Ac-th g-axis amortisseur winding, results in the following expression: L Dq a  ~L  L Dq  q  a  '<pDql  where  LaDq  L Dq  L Dq  L Dq  a  a  a  Lo\\q  LD12q • •  Loijiq  i-Dql  i>Dq2  — L Dq  Lmiq  Lm2q • •  LD2HQ  i-Dq2  i>Dqv.  — L Dq  Lo\jLq  LD2p,q  LDnnq  i-Dq/j.  L  Dxyq  a  a  2L  xyq  L )  / (3 L  aDq  L )  axq  "  (7.34)  • The resistance and inductance values corre-  ayq  sponding to L DKKq  L DqK  and  L Dq a  R  DqK  =  2L aDq 6o n  7-2  anq  ?  K<  for K — 1 . . . fj, are depicted in the g-axis circuit diagram in Fig. 7.1. Their values can be obtained from field test measurements as explained in Chapter 8. The magnitizing inductance can be calculated from L Dq — L — a  where the leakage inductance Li is available from  q  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  and  ip  = - 5 'ijj  DdK  K  Kd  OK  in (7.31), where X — °K  —  ^ LfDd L f a  \l  d  TTT ) J^fKd, J^aKd T  Z  results in the following matrix equation: n^d]  -L  — Ld  I'^Jd  '4>Ddv  Lafd  Lafd  —L fd T  — L fd  d  -Lafd  a  2  a  LfDd  -Lafd  LfDd  L\\  _ -Lafd  LfDd  LfDd  •• - -  •  • ••  d  ••  Lafd  ' Wd\'  -L fd a  Lf  Dd  LfDd  i-Ddi  Lpvd  ^Ddv  (7.35)  103  7.3. Model Equations with Complex Signals  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. L i in (7.31) xyc  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 =  LDOIK,  LDKKCI  Lj  —  and  Dd  RDCLK  —  2  ^ <5 R K  Kd  where K = 1... v and where L f = L — L( is the magnitizing inductance. The peripheral a  leakage inductance Lj  d  d  — Lf —L f  kd  Dd  the relationship L = L — Lf . c  d  is related to Canay's characteristic inductance through  a d  Due to the lack of proper measurement methods, this  Dd  inductance is usually set to L = Lu in traditional machine representations. However, newer c  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,  E  d  E q  = IVi-Ralt-Ur*, at = ±y -R I n  a  q  (- ) 7  +  u; * r  36  (7.37)  d  at e  fd  =  0 =  ^ipfd  + Rfd'^/d  (- ) 7  + RDdnlDdK  ^TlpDdK  dt  0 = ^1pDg + RDq^Dq* K  forK-=l...i>  for  K=  1.. . U .  39  (7.40) (7.41)  7.4. Time-Domain Simulation  7.3.3  104  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 T  m  = - ^ (B[* ]H[J,]-»[*,]»[J ]) d  •  r f  (7.42)  The deviations in the mechanical rotor position are given by 1 A0  dt  m  =  — (u - UJ )  n  =  (7.43)  0  r  p  2  oo  wo ,  m  (7.44)  where oj is the mechanical rotor speed and where n /2 is the number of pole pairs. m  p  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 universal machine model in the EMTP [24]. It allows the translation of the mechanical equations into their electrical equivalents. For instance, for the torque:  T (t) = J j uj {t) m  where J  rn  =n E p  /(2u) )  kinetic  the kinetic enerergy E  0  kinetic  m t  i(t) = Cj v{t)  m  (7.45)  t  is the moment of inertia. The inertia constant H is defined as divided by the power rating of the machine. While most power  system analysis tools take the machine's apparent power rating S i rat  ng  to build this ratio, a  normalization with the real power rating P ting was chosen in [22]. ra  7.4  Time-Domain Simulation  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.  105  7.5. Steady-State Initialization  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 I {t), hit) a  current ifd(t), and the electromagnetic torque T (t).  and I (t). the field c  As we can see, these quantities are  m  functions of the flux linkages and the rotor position 0(t) = 9 + A9{t) + u t . 0  (7.46)  0  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[/ ] and 5R[/ ], needed to calculate rf  g  the electromagnetic torque according to (7.42). Applying Park's inverse transformation (2.6) gives the complex armature currents I (t), I (t) and I (t). a  0  c  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 u  m  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 Initialization  During steady state, the dqO transformation of balanced steady-state armature quantities produces complex, constant dqO quantities. The quantities associated with the rotor windings 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.  106  7.5. Steady-State Initialization With the electrical rotor speed  the voltage equations (7.36)-(7.41) become  Ed  =  -Rald-OJo^a  (7.47)  E  =  -Ralq+UoVd  (7.48)  e  =  Rfdifd  (7.49)  q  fd  0 =  RDCLK iDdn  0 =  RDqniDqK  for  1...U  K =  for  (7.50) (7-51)  K = \ . . . fl .  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 =  -L ft[I ]  + lLf i  4>DdK =  -L $l[I }  +  •tpDqn =  -L U[I ]  jipfd  afd  d  afd  d  aDq  (7.52)  fd fd  ior  L ifd fDd  (7.53)  K - 1 . . . U  for« = l.../z  q  (7.54)  and for the stator flux linkages: V  d  V  q  =  -L $t[I }-jL %[I ]  =  -L M[I ]-jL %{I }-jL i  d  d  q  q  q  +L i  d  d  q  (7.55)  afd fd  (7.56)  afd fd  Because of \& = j ty , both equations are linearly dependent. In order to initialize the above d  q  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 5  l  as the angle by which the g-axis leads the phase a reference axis, the following relationships can be established: E =E q  where E(u>o) = E e  e* -W  / , = /„ *(°+*-*)  a  n  ;  and I(u>o) = I e^ ^  ja  are the steady-state phasors obtained with e.g.  a+  n  (7.57)  e  n  a power flow solution. Substituting (7.55) in (7.48), we have E + RI+ q  a  q  jUJ L 3[I ] 0  q  = -OJQ L U[I ]+  d  d  By expanding both sides by the term LJQ L U[I ] q  d  d  uj L i 0  afd fd  .  (7.58)  and observing the relationship I — j I , d  q  the equation becomes E + (R q  a  + joj L ) 0  q  I = cu [L q  0  q  L ) $l[I ]+ UJ L i d  d  0  afd fd  .  (7.59)  107  7.5. Steady-State Initialization  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: K = UJ ( L - L ) I sin [ a + <> / - S ] + co L 0  l  q  d  n  l  0  afd  i  fd  .  (7.60)  Then, we can determine Ki and the angle Si with  Ki e  jSi  = E(u ) + {R +juj L ) 0  a  0  q  I(co ) 0  (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 equations 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 D C 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 measurements 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 measured 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 circuits, N and N , on the right-hand side of Node A in the diagrams of Fig. 7.1—directly d  a  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 measure these functions are well defined in the I E E E 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 A e / = 0. In the Laplace domain A is used to denote a small change, and d  s = ju is the Laplace operator. The operational inductance in <i-axis direction is:  110  8.1. Frequency Response Tests 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 Ai with Q(s) being referred to the stator side. What is actually measured is: d  G(s) = hG(s),  (8.4)  where h = §N /Nyy is the turns ratio between the armature and the field winding. Notice a  that this ratio is only needed if the per unit values of model parameters are desired (i.e. L -base reciprocal system [55, 75]). If the correct unit values of the previous chapter are ad  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: L  d  L  q  =  lim£ (s) d  =  lim£,(s) s—>0  G =h ^ 0  =  Rf  n g< ) m  s  d  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: L fd a  L Dq a  — L — Le d  =  L  q  — Ll-  The armature resistance R is obtained as part of the process of constructing the operational a  inductances C (s) d  and C (s). q  The remaining unknown parameters are those associated  with the rotor windings, more specifically, the values of Rodn, R , DqK  constitute the rotor RL-branches in addition to Lf . kd  parameters is the subject of the following paragraphs.  L , DdK  and L  DqK  that  A possible way to identify these  111  8.2. Transfer Functions  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) = F ^  =E7J^  0  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  di K~^  K  v = -R i K  -  K K  L  (8.6)  with the current flowing out of the branch. After the Laplace Transform, we have -Az  K  Av  K  l/L (s + R /L )' K  K  (8.7)  K  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  JVd-1  A -  Ad _  1X1  Av  Ad  ( 8  _ J _  L (+ f ^ ) DdK  ^  S  .  8 )  J-  Lfd  /g g\  (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  112  8.2. Transfer Functions Lambton GS: Magnitude and Phase of Aq(s)  CO  fin  Frequency [Hz]  Figure 8.1: Admittance A (s) third-order polynomial. q  voltage Av  q  of the Lambton generator and its final approximation by a  and terminal current ( — Ai ).  Finally replacing this ratio by sC (s), yields the  q  q  g-axis transfer admittance  A (s) = q  -Ai Av  Aq  Aq  1  1  1  s \C (s) - L q  e  (8.10)  L  aDq  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 Ljk were omitted in the d-axis model, its characteristic equation d  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 armatureto-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  B (s) d  -Ai Av  Bd  Bd  1 / 1 s \£d(s) - L  e  1 L :fd a  (8.11)  8.2. Transfer Functions  113  where the voltage and current relationship between the two Nodes A and B is given by  •%r- = i m - '<"  <>  si  -Ai  B (s)  Ad  Notice that Ai  Ad  and Ai  812  d  are identical. If L^d had a known value, we could easily write  Bd  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 = Ai  - ^  Ad  (8.13)  SJ^afd  we arrive at a transfer function describing the ratio of the current through Lfkd, namely Ai  Ad  to the field current Aif : d  C(s) = ^ h Ai  =J - ( g( ) V L L  fd  s  d  C  s  d  {  s  afd  )  ) . J  (8.14)  Expanding both numerator and denominator by the term Av d A  allows us to define yet  another transfer function: f L f d [ S )  ^ _ 1 Av ~ sh (-Ai ) Ad  fd  _ 1 Ai Av - sh Aifd (-Ai ) Ad  Ad ( 8  Ad  -  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 Av h{-Ai )  s + rj /l l/l •  Ad  d  fd  where  Sd  fd  r  = \ R = 4^-. h Q This equation together with (8.12) and (8.14) inserted into (8.15) and solved for l fd  fd  fd  '^UTT- '"^ 1  v  v  '  a  Since both the unknown parameters lf and Lf d  kd  (817)  b  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)  gives:  114  8.2. Transfer Functions Lambton G S : Magnitude and Phase of the 0-th Residue  m -a P  'S a3  cS  Frequency [Hz]  Figure 8.2: Admittance A (s) of the Lambton generator approximated by a first-order function of type .7-1 (s). q  Theoretically, there is a chance that S(6) could become zero causing Lf  kd  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 Lf , the characteristic equation kd  of the d-axis transfer function becomes: A (s) = d  -Ai  Ad  Av  Ad  B (s) d  (8.20)  l-sL B (s) Skd  d  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 denominator polynomials constituting the transfer functions C (s) and C (s) with respect to s are d  q  115  8.2. Transfer Functions  equal. This implies that the order of Ad(s) and A (s) is higher by one in the denominator g  than in the numerator. Dropping the index, we can write for both functions:  A(s) =  ^^Y[p^4, s+  PN  ii  (8.21)  (s+p y  =  K  where N is the number of the rotor circuits of the ci-axis or g-axis model, i.e. either  or  N . Hence, there are two different types of fractions appearing in the polynomial: q  s +p T (s) 2  = k -±^  (8.23)  S  s+p  The pole p and gain k of the first function can be obtained through the following two relationships:  p—  —— and  k=p\\m.T\(s)  c  tan^i  *->o  Principially, the phase angle <pi could be taken anywhere together with the associated frequency 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 A (s) by a first-order q  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 T (s) 2  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) and its 2  associated frequency co : 2  ,  1 + sin </>  1 - sin </>,  2  k = min (  2  ,  1 — sin 02  j  1 + sin <f>  2  r~ ^2 p = uj vk and z =-j= 2  Having found the pole-zero-pair, the new residue is formed by dividing the present one by the most recently determined function T (s). This extraction procedure is repeated until 2  there are no more noticeable phase shifts present in the residuals. The transfer function can  116  8.2. Transfer Functions Lambton GS: Magnitude and Phase of the 1st Residue  OH  10  10  10  Frequency  Figure 8.3: First residue of A (s)  10  10  [Hz]  approximated by a first-order function of type ^ ( s ) .  q  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  (s + z )  F]  A{s) = A  k  N  K  K = l 0  (8.24)  K  N K=l  N-l  n  where  k =A K  0  (n z  PK)  N  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 A (s)  polynomial of the order N = 3 was already shown in Fig. 8.1. q  q  by a third-order  117  8.3. Practical Considerations  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 C (s), C (s) and sQ(s) are all directly obtained from measurements, there will be d  g  some measurement inaccuracies superimposed on the true data. This may cause some calculated 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 A (s) that even small inaccuracies in the q  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 Lf . Yet, only in theory does the kd  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 Lj ?  Going through the entire frequency range of measured and  kd  filtered data points, L  fkd  and l  fd  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 Lf  kd  evaluate the right-hand side of (8.15), while lj and r d  fd  we can  are used to determine the left-hand  side, using (8.16). In theory, both sides should give the same results. In practice, however,  118  .3. Practical Considerations 100  10  10  10  10'  10  10'  1000  -1000 0  -2000 10"  10"  10-  10  u  F r e q u e n c y [Hz]  Figure 8.4: Parameters Lj and lf as functions of frequency. The cross indicates the best fit, yielding L = -0.0282 p.u. and l = 0.831 p.u. kd  d  fkd  fd  there is a mismatch between both. Calling the right-hand side RIGHT(L ,s)  =  fkd  C(s) sB (s)  LfkdC(s)  (8.25)  d  and the left-hand side = l + fd r  LEFT(L ,s) fkd  (8.26)  fd  an error function can be established as , ,  (  . J  W  NoF  / RIGHT'(Lfkd, St) ~ 2-. ( LEFT{Lfkd, Si) ~  (8.27)  where NoF is the number of discrete frequency points at which the data were measured. The pair of Lf  kd  and lf that gives the least error is then selected for further calculations. d  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(Lf ) tends to be greater than zero kd  even with perfectly selected inductances Lf  kd  and lf . This may have two reasons: First, d  119  8.3. Practical Considerations Lambton GS: Magnitude and Phase of Lfd(s)  bD  CO  JS PH  Frequency [Hz]  Figure 8.5: Best concurrence of functions RIGHT indicates the optimum value of Lf .  (solid line) and LEFT  (dashed line)  kd  the separate measurement and independent filtering of C (s) and sQ(s) may cause both d  transfer functions to be slightly uncorrelated, as a result of which RIGHT(Lf , kd  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 C (s) or sQ(s) or both are d  slightly inaccurate to give the desired first-order term. Let us define the mismatch between RIGHT(Lf , kd  s) and LEFT(Lf ,  s) as a residual function g(s) such that  kd  LEFT(L , fkd  1  s)  g(s)  RIGHT(Lf ,  s)  kd  (8.28)  aC (s) - b  (8.29)  d  g{s)sg(s) with  a= 1  L fkd L }d a  b  T -L + T L>1 Lf L f k d  a  L. d  d  Unfortunately, there is no way of determining which part of g(s) is due to a mismatch in C (s) and which is caused by inaccuracies in sQ(s). Because of the difference of scale—L is d  d  120  8.4. Results  much smaller than G —the operational inductance Cd(s) is far less sensitive to errors than 0  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(Lf d, s) = - J — (aC (s) - b) , k  sy(s)  v  d  (8.30)  '  where C (s) denotes the corrected operational inductance. Its use ensures a pure first-order d  RL-branch to represent the field winding. Setting (8.29) and (8.30) equal, yields  £  W  =  ^ W - * ( JL g(s)  (8.31)  a \g(s)  )  Using this inductance instead of Cd(s), we can now calculate the transfer admittance A {s) d  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°  10  1  10  2  Frequency [Hz]  Figure 8.6: Measured operational d-axis inductance C (s) with its second-order analytical (solid line) and iterated approximations (dashed line). d  121  8.4. Results Lambton GS: Magnitude and Phase of sG(s) Oi  :  : : • • • • •.  :  : •  :  : : :  1 : : : ::  3  10'  10''  10°  10  10"  10'  10'  10°  io  10'  2  1  10  2  100i  OH 3  2  1  1  io  2  Frequency [Hz]  Figure 8.7: Measured armature-to-field transfer function sQ(s) with its second-order analytical (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°  10  1  10  2  Frequency [Hz]  Figure 8.8: Measured operational g-axis inductance C (s) with its third-order analytical (solid line) and iterated approximations (dashed line). q  122  8.4. Results  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 p.u. p.u. L fd p.u. L Dq Noli p.u. Rfd R-Ddl p.u. Lfkd p.u. p.u. Lfd Lndl p.u. N p.u. RDql Roq2 p.u. Roq3 p.u. Luql p.u. Ll)q2 p.u. p.u. L®q3 a  a  q  Notice that L f and L D a  d  a  q  Values 0.155 1.7351 1.6450  Ref. [19] 0.155 1.8580 1.7620 2  0.2 0.000856 0.4482 -0.0282 0.1840 0.1924 3 0.00505 0.01343 0.10445 3.8647 0.4842 0.1239  n/a 0.00108 0.01065 0.1328 0.0105 0.0114 0.01060 0.0301 0.2610 2.1010 0.4504 0.1131  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 C (s) and sQ(s) than the parameters obtained iteratively. d  This is especially true for the high-frequency end of the phase curve of C {s). d  Larger deviations can be noticed in the curves for the operational g-axis inductance C (s). q  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. P O W E R F A C T O R Y ® Power system analysis package by D I G S I L E N T 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  125  9.1. Transmission Line Energization  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 L i n e E n e r g i z a t i o n  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.  AC  202.8 km  @  V\A/  a,b,c(V  e  ^Thev  r^r^  (^  lo  '"Thev  ^ line' ^ 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 signals, hence E (t)  =  1.0  E (t)  =  1.0e  E (t)  =  1.0e  a  b  c  where  UJ  0  j { l O 0 t e  -  p.u.  7 T / 2 )  p.u.  (9.2)  p.u..  (9.3)  j ( w o t + 5 7 r / 6 )  , (  ^  / 6 )  (9.1)  = 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) = R hev I(t) + L T  Thev  j  t  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) H{t)  = Y V(t) Thev  = F V(t-At) Thev  + H(t)  (9.4)  + G I{t-At) Thev  (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 peak magnitude phase shift in a phase shift in b phase shift in c  50 Hz 1.0 p.u.  Thevenin impedance  270° 150° 30° 6.75 Q 127.0 Q  RThev X-Thev  for all sequences Transmission line  length I  202.8 km  RLe, os  0.04 Q/km  P  X\ ^ m  0.318 Q / k m 11-86 n F / k m  os  C[ e, os Kne,zero XLe,zero in  P  0.26 Q/km ^ tyl™  C\ ,  7.66 n F / k m transposed, balanced initial status open at £ = 3.05 ms in a closed at £ = 5.55 ms in c closed at £ = 8.05 ms in b closed ine zero  Circuit breakers  where  At  At /  *  a  Thev  —  Thev  —6  n  and where A = -L~ R  m  l  hev  Thev  y-L'Thev  =  ^ Thev  -^-[L'Thev  e  -  ^  \ - i  - ^ f  {  l  J  A  7^  - f  A  y  J  T h e v  {  l  +  *  A  )  ,  ju l. Q  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  127  9.1. Transmission Line Energization  Simulation results (0-40 ms)  2.5  r  ~i  1.5  p  CU hO co  >  0.5  n u  -0.5  -1.5  solid line: M I C R O T R A N dotted line: New program  -2  -2.5 0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  Time [s] Simulation results (0-800 ms)  2.5  solid line: M I C R O T R A N dotted line: New program  -2h  -2.5 0  0.1  -j 0.2  i_ 0.3  _i  0.4  0.5  i  0.6  0.7  0.8  Time [s] Figure 9.2: Simulation results with a C I G R E line energization test case.  128  9.1. Transmission Line Energization  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 R M S voltages and currents instead. Since we chose complex signals for our variable representation, R M S 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 At  = 5 fis for the  min  backward Euler steps and then continues with 2At  = 10 /is. However, once the high  min  frequency components of the electromagnetic transients disappear, the simulation step size is gradually increased until it reaches the user-defined maximum value At  = 50 ms at  max  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:::::::::::::::::;:::::::::::;::  1  0  0.1  ::: i ::::::::::-  :::::: : : : r: :  —  :  :::::::::::::::::  :  :  :  ::::::::::: :: :  :::: :::::::::: -  :  :  :  :  :  MlCROTRAN  1  l  l  l  l  l  0.2  0.3  0.4  0.5  0.6  0.7  Time [s]  Figure 9.3: Step size adjustment (0-800 ms).  0.8  129  9.2. Subsynchronous Resonance Case  9.2  Subsynchronous Resonance Case  In this section, we investigate electromagnetic and electromechanical transient phenomena with the new simulation program. The I E E E subsynchronous resonance simulation benchmark 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 electromagnetic transients programs and stability analysis software.  Figure 9.4: Network diagram for I E E E 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  u M V A rating  V  0  1  3X trafo  -V-  V~  1  2  l  HV  v  HV  V  LV  V  v  LV  (9.6)  LV  V  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^ k V and V  LV  = 26.0 k V .  For simplicity and in order to compare the results with a stability analysis of this case,  130  9.2. Subsynchronous Resonance Case Table 9.2: I E E E subsynchronous resonance benchmark case data. Equipment / Events  Parameters  Voltage sources (infinite bus)  frequency 60 Hz peak magnitude 361.9 k V phase shift in a 266.187° phase shift in b 146.187° phase shift in c 26.187° X 16.809 fi in all three phases C is 25.5178 (J.F in all three phases Rline,pos 5.603 fi X „ 140.07 fi Rline,zero 140.07 fi Xline,zero 437.02 fi transposed, balanced M V A rating 892.4 M V A voltage rating H V 500.0 k V voltage rating LV 26.0 k V connection type Y N d 01 Xtrafo 0.14 p.U. three single-phase transformers X 11.206 fi in all three phases 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  Thevenin reactance Series compensation Transmission line  Thev  ser e  u  Step-up transformer  Fault reactance Fault sequence  Values  s  fault  the transmission line model was reduced to a three-phase RL-series branch. A n 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  131  9.2. Subsynchronous Resonance Case Continuation  of Table 9.2  Equipment / Events Synchronous generator  Values 892.4 M V A M V A rating voltage rating 26.0 k V power factor 0.9 (lag.) connection type Y N , grounded 0.000 p.u. 0.130 p.u. 1.790 p.u. x 1.710 p.u. x 0.130 p.u. x 0.169 p.u. K x 0.228 p.u. \ 0.135 p.u. K 0.200 p.u. K rp, 4.300 s do V 0.850 s rpll qo 0.032 s do rpll 0.050 s 0° H 2.894 s pole pairs 1 field rating if 600 A constant voltage Power rating 803.16 M W constant torque, one mass Parameters  d  q  0  1  1  Excitation Turbine  0  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 trapezoidal 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 propagates 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  132  9.2. Subsynchronous Resonance Case  2.145  Comparison with  x 10  M I C R O T R A N  i  2.14  J: l \ ,l 1  |  . . 1 V  2.135  L  i' ll  2.13F  1  1  i !  1 1  |  1  ii  i'  —  :  1 U  i  1 1  I  1 w  1  ' I I I  '  1  ' ,  1  1  1  / 1 1 " )> 1 1  1  w  •  1  U i l' v'  1  2.125 1  i : i  •  n  ' i  \ 1 1  V i  \ X . it  '  i ./ '  i i ' i l i i [ ' l, i H i ''' il "  V\  \>  • • • / • ! •  v  II  i  A  1  '  1  :  , 1  1  II. . 1 . . .1 . . . ll'  . . .  l"1  ' • i»1 • 11  1  l\  "  >\  l\  h  i'  ^  !\; "  '../  i  1  2.12  1  2.115  :  i 2.11  ' :  das hed lint: : M I C R O T R A soli d line: ^ew program  i  1  2.105  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  N  0.9  Time [s]  2.135  Comparison with PoWERFACTORY ( E M T )  x 10 I  I  I  I  !  2.134  2.133 -  -  2.132  o  ' ' , I  ,  2.131 2.13  ']'•• 1  L  :  u  '  I  V 1 ll ' 1 ' I ' II ; ' . I i1 i . M y i , '  r  -  ll  Ii  ' V •  J  _l  1 1  J  J  -  '—: '  :  J  2.129  2.128  2.127  dashed line:: P O W E R F A C solid line: New program  2.126  2.125  i  0.1  I  0.2  0.3  0.4  0.5  0.6  I  I  0.7  T O R Y  I  0.8  Time [si  Figure 9.5: Steady-state electromagnetic torque  0.9  133  9.2. Subsynchronous Resonance Case  (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 POWERFACTORYsoftware. The electromagnetic torque is shown in Fig. 9 . 5 (dashed line). Using the electromagnetic transients setting ( E M T ) 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 amplitude. 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 3601 |  I  I  1  1  M I C R O T R A N  1  1  1  1  1  3600.8 -  -  3600.6 -  -  3600.4 g.  3600.2 -  i_  3600  ^••-i.-j.-^-:  j :  -  ;  ;  CO  -g  \  •  \ .  3599.8 -  -  3599.6 -  -  3599.4 -  -  dashed line 3599.2 h 35gg I 0  :  M I C R O T R A N :  solid line: Is'e\v program 1 0.1  I 0.2  I 0.3  I 0.4  I 0.5  I 0.6  .  I 0.7  T i m e [s]  Figure 9 . 6 : Steady-state rotor speed  I 0.8  I 0.9  I 1  9.2. Subsynchronous Resonance Case  9.2.2  134  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 At  min  = 5 /is and a maximum step size of At  max  = 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 compensator 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 Nm/rad occurs shortly after the fault clearing. This value 6  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 At  max  for a simulation of the electromechanical transients.  9.2. Subsynchronous Resonance Case  Figure 9.7: Transients simulation (0-1 s): Electromagnetic torque.  135  136  9.2. Subsynchronous Resonance Case  Comparison with  P O W E R F A C T O R Y ( E M T )  i  1  3620  -  /  u  \ i  g_ 3600 3 o  ?, / o  \1 3580  dash ed line: POWERFACTC ):RY solic line: IS ew program 3570 0  0.1  0.2  0.3  0.4  0.5  i  i  i  I  0.6  0.7  0.8  0.9  Time [s]  Figure 9.8: Transients simulation (0-1 s): Rotor speed.  1  9.2. Subsynchronous Resonance Case  Figure 9.9: Transients simulation (0-10 s): Electromagnetic torque.  137  138  9.2. Subsynchronous Resonance Case Simulation with the new program 36301  1  1  1  1 2  1  1  1  1  1  !  r  3620  3570 I 0  1  1  1 3  4  1 5  1  1 6  7  1  1  8  1 10  9  Time [s]  3630 |  i  Simulation with  |  1  1  M I C R O T R A N  1  1  1  !  r  3620  3570 0  1  1  1  2  L —  1  3  4  1  5  1  6  1  7  1  8  1  9  Time [s]  Figure 9.10: Transients simulation (0-10 s): Rotor speed.  1  10  139  9.2. Subsynchronous Resonance Case IEEE Subsynchronous Resonance Case  10  1  • 1  10"  • New simulation program  10  MICROTRAN  10  10  I 0.2  _l  0  0.1  I 0.3  I 0.5  I  0.4  Time [s]  I 0.6  1 0.7  L_  0.9  1  IEEE Subsynchronous Resonance Case  10  10 '  In  I  0.8  :  New simulation program -j  10"  10"  MICROTRAN 10  10  _i  0  1  2  3  i _  i  4  5  6  Time [s]  7  8  9  10  Figure 9.11: Transients simulation: Automatic step size adjustment.  140  9.2. Subsynchronous Resonance Case  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 m i n i m u m or m a x i m u m 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, during large acceleration or deceleration of the rotor, the instantaneous frequency also changes significantly and causes the simulation time step to drop to A t =  2At . min  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 r i so  solution of one time step and T u reou  d  ution  for the  for the construction of the nodal admittance m a t r i x  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 time 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 t h the results from a transient stability analysis. For this purpose, the setting of the P O W E R F A C T O R Y program was changed to R M S for the transient stability simulation.  I n addition, 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 contribution to the electromagnetic torque and rotor speed are missing. This can be found from the absence of the 2nd subharmonic component in the solution, from the instantaneous 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. I n the remainder of the simulation, the mechanical torque and rotor speed by the stability program approximately build the mean value to the subharmonic oscillations. After these oscillations have disappeared, the solution obtained from  9.2. Subsynchronous Resonance Case  Figure 9.12: Stability simulation (0-1 s): Electromagnetic torque and rotor speed.  141  9.2. Subsynchronous Resonance Case  142  3580  3570  10 Time [s  Figure 9.13: Stability simulation (0-10 s): Electromagnetic torque and rotor speed.  143  9.3. Evaluation of Program Performance  the stability program converges to the same steady-state solution as obtained with the other simulation programs.  9.3  E v a l u a t i o n of P r o g r a m P e r f o r m a n c e  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 accordance 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 chapter are listed in Table 9.3. A l l simulations were carried out on a 200 M H z Pentium Pro computer. Table 9.3: Program running times for several test cases. Simulated case Line energization Steady-state simulation Subsynchronous resonance (1 s duration) Subsynchronous resonance (10 s duration)  MICROTRAN  0:23 min 0.3 ms/step 1:29 min 0.9 ms/step 1:30 min 0.9 ms/step 14:31 min 0.9 ms/step  P O W E R ]F A C T O R Y  EMT-mode n/a  RMS-mode n/a  11:14 min 6.7 ms/step 11:14 min 6.7 ms/step n/a  2:13 min 133 ms/step 0:08 min 8 ms/step 0:28 min 28 ms/step  New program 0:08 min 7.2 ms/step 0:17 min 21.9 ms/step 16:23 min 208.2 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 convergence 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  R e c o m m e n d a t i o n s for F u t u r e Research  In this thesis, a concept for the combined simulation of electromagnetic and electromechanical 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 prevailing 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  146  10.1. Recommendations for Future Research  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  EMTP  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  EMTP  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  10.1.3  147  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 simulation 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 inflation 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.  148  10.2. Summary  10.2  Summary  With regards to time-domain power system simulation, the areas of electromagnetic transients studies and stability analysis have traditionally been isolated for several reasons. However, recent system failures and power outages have created a demand for more comprehensive 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 similarity 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 derived. 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  EMTP  [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 D C 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 I E E E 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) Ag{s) A  B (s) C(s) d  & [ . . . ]  e (t),  e (t)  a  b  E (t), E (t), d  q  e (t) E (t) c  0  fd(t)  e  fmax  Q(s) G  c  H(t-At) h H  @{I(t)} IM i(t) I(t) i (t), ib(t), i (t) I (t), I (t), I (t) a  d  yd{t) j Jin  V t  L  Ms) £jd{s) Li LaDq Lafd LDdK  c  q  0  Direct-axis transfer admittance Quadrature-axis transfer admittance Aggressiveness factor for simulation Armature transfer admittance Rotor current to field tranfer function Dynamic phasor Generator terminal voltages in phases abc. Complex generator voltages in dqO coordinates Generator field voltage Maximum frequency Armature-to-field transfer function Characteristic wave conductance Vector of complex history values Turns ratio between armature and field winding Inertia constant Current vector (dynamic phasor) Imaginary part Current vector (phasor) Current vector (real signal) Current vector (complex signal) Generator terminal currents in phases abc. Complex generator currents in dqO coordinates Generator field current Imaginary unit Moment of inertia Gain of a transfer function (Chapter 8) Length of a transmission line Inductive matrix Direct-axis operational inductance Field winding operational inductance Quadrature-axis operational inductance Armature leakage inductance Quadrature-axis armature to rotor mutual inductance Direct-axis armature to rotor mutual inductance Direct-axis amortisseur inductance 151  152  List of Symbols Quadrature-axis amortisseur field  inductance  LfDd  Direct-axis  to amortisseur mutual inductance  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  Field  Lfkd  Peripheral leakage  rn. N n N NoF V V q R »[...]  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)  Ra  Armature  fd  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  inductance inductance  d  N u m b e r o f d-axis rotor windings  p  N u m b e r o f poles. T w i c e t h e n u m b e r o f pole pairs N u m b e r o f g-axis rotor w i n d i n g s  q  N u m b e r of discrete frequency  points  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) 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) 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) Resistive matrix Real part  r  Rfd  Field  S(u ) s  Phasor  0  resistance  resistance  Laplace  operator  S(UJ)  Spectrum of a signal  s(t)  Real signal  S*(t)  Conjugate  t  Time  T to  O n e period, cycle  tend  Simulation end time  T {t) 9[V(t)] V(cu )  Electromagnetic generator Voltage vector  v(t)  Voltage vector (real signal)  V(t) X(t)  Voltage vector (complex  X  Distance o n a transmission line  of a complex signal  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  m  Voltage vector (dynamic  0  (phasor)  C o m p l e x state variable  Z[...] z Jreek  torque phasor)  signal) vector  Z-transformation of a function 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)  Symbols  P  A n g l e f r o m phase a reference  A  S m a l l change of a quantity  axis t o d-axis  AAp di i d  Predicted m a x i m u m magnitude deviation  At  S i m u l a t i o n step size, t i m e step  re  At i m  c C  n  M i n i m u m s i m u l a t i o n t i m e step  Atp dicted re  P r e d i c t e d s i m u l a t i o n t i m e step  At  r n a x  M a x i m u m s i m u l a t i o n step size  t  L o c a l numerical error  List of Symbols e e r/(...) 8(t) 8i(t) K, A A p, Pi p p fj,4 pchange pdelay Pefficiency 3  p tp em  Pfiuct Pmagnitude p ode Preiiabie psample p tab p tep v Wj(t) g(s) a T m  s  s  Trebxdid  Ti T 4> </>(...) ip(t) ip(t) so ution  MVN  '0a(t), 'tpb(t), '^cif) ^q{t),  ^o{t)  (u) LO wo ^unear(t) upredicted,max u^dieted,mm upredicted u> (t) m  Global numerical error Tolerance margin for local numerical error Error function Angle from phase a reference axis to ci-axis Internal machine angle from phase a-axis to g-axis Index for generator rotor windings Eigenvalue Eigenvalue matrix Number of q-axis amortisseur windings, see N Membership degree "Frequency gradient is low" Membership degree "Frequency is close to nominal" Membership degree "Number of samples is high" Membership degree "Sampling window is large" Membership degree "Likelihood of a step change is high" Membership degree "Duration since last step change is long Membership degree "A step change is efficient" Membership degree "Simulation is in EMTP mode" Membership degree "Frequencyfluctuationis low" Membership degree "The magnitude deviation is small" Membership degree "Simulation mode is well defined" Membership degree "Forecast reliability is high" Membership degree "Sampling quality is high" Membership degree "Simulation is in stability mode" Membership degree "Step ratio is high" Number of d-axis amortisseur windings Instantaneous frequency Residual function Standard deviation of frequency samples Modal wave propagation time Time required to rebuild the admittance matrix Time required to advance the solution by one time step Minimum modal wave propagation time Phase angle, phase shift Analytical solution Partial solution to wave propagation equation (Chapter 6) Partial solution to wave propagation equation (Chapter 6) Generator fluxes in phases abc Complex generator fluxes in dqO coordinates Average frequency Angular frequency System angular frequency, power frequency Regression line through frequency samples Predicted maximum frequency Predicted minimum frequency Absolute predicted frequency Mechanical rotor speed q  2  ^d{t),  153  154  List of Symbols  0J (t) r  Nominal frequency of a quantity Electrical rotor speed  Acronyms AC BE CIGRE  IEEE EMTP NETOMAC RCV RTR TNA TR UBC UBC  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. I E E E Press, Piscataway, N J , 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, I E E E 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, B P A , Portland, Oregon, Electromagnetic Transients Program Rule Book, April 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 Apparatus 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, April 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 BerlinHeidelberg-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  BIBLIOGRAPHY  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, N J 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, April 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, April 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, April 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-178605. [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, Version 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, April 1969. [24] H . W . Dommel, EMTP Theory Book, Microtran Power System Analysis Corporation, 4689 W . 12th Avenue, Vancouver, B . C . V 6 R 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 ( E L E C 553), Dept. of Electrical & Computer Engineering, U B C , Vancouver, B.C., Canada, 1998, available at http://www.ece.ubc.ca/power/e553.pdf.  157  BIBLIOGRAPHY  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 A v enue, Vancouver, B.C., Canada V 6 R 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 E P R I 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, April 1986,  E P R I 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 University 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, N J 07632, U.S.A., 1971. R. W . Hamming, Digital Filters, Prentice-Hall Signal Processing Series. Prentice-Hall, Inc., Englewood Cliffs, N J 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 BerlinHeidelberg-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, I E E Monograph Series 4. Cambridge University Press, London, U . K . , 1970. S. Henschel, "Long-Term Evaluation for Power Systems under Consideration of Electromagnetic and Electromechanical Transient Phenomena," P h . 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 Identification from Frequency Response Tests," IEEE Transactions on Power Systems, May 1998, accepted for publication, ref. no. PE-003-PWRS-0-05-1998.  BIBLIOGRAPHY  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 Computation Conference, Trondheim, Norway, 1999, accepted for publication, ref. no. 22. [43] J. D. Hurley and H . R. Schwenk, "Standstill Frequency Response Modeling and Evaluation 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] I E E E , "Standard 115-1983," in IEEE Guide: Test Procedures for Synchronous Machines. Institute of Electrical and Electronics Engineers, Inc., September 1983. [45] I E E E Joint Group, Chairman: P. L . Dandeno, "Supplementary Definitions and Associated 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] I E E E 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] I E E E 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," Archiv 3, vol. 8, pp. 279-281, May 1.981.  Etz-  [49] E. W. Kimbark, Power System Stability, vol. 3. Synchronous Machines, John Wiley & Sons, Ltd., New York, N Y , 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 kontinuierlichen und diskontinuierlichen Systemen," Siemens Forschungs- und Entwicklungsberichte, vol. 10, no. 5, pp. 299-302, 1981. [52] B. Kulicke, "Simulationsprogramm N E T O M A C : Nachbildung von Synchron- und Asynchronmaschinen," 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.  BIBLIOGRAPHY  159  [54] B . Kulicke, "Digitale Simulation elektromagnetischer und -mechanischer Ausgleichsvorgange," Course Notes, Institut f. elektrische Energietechnik, T . U . Berlin, Berlin, Germany, 1991. [55] P. Kundur, Power System Stabiliy and Control, E P R I 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 BerlinHeidelberg-New York, 10th edition, 1973, ISBN 3-540-06021-9. [57] N . Levanon, Radar Principles, John Wiley & Sons, Inc., 1st edition, 1988, 0-471-85881-1.  ISBN  [58] U . Linnert, Berechnung von Ausgleichsvorgangen in Elektroenergiesystemen unter Verwendung 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 Electromagnetic Transients Simulations," IEEE Transactions on Power Apparatus and Systems, 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 ThyristorControlled Series Capacitor Systems," IEEE Transactions of Power Systems, vol. 12, no. 3, pp. 1259-1267, August 1997.  BIBLIOGRAPHY  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 Modulated 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 Computation 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 Economics, 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 vol. 64, August 1945.  Transactions,  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, N Y , U.S.A., 1974, ISBN 0-471-79420-1. C. P. Steinmetz, 1915.  Theoretical Elements of Electrical Engineering,  McGraw-Hill, Inc.,  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.  BIBLIOGRAPHY  161  [82] S. D. Umans, J. A . Mallick, and G . L . Wilson, "Modelling of Solid Rotor Turbogenerators," 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, Berlin-Heidelberg, Berlin, Germany, 1965.  Springer Verlag,  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065161/manifest

Comment

Related Items