Multilevel M A T E Algorithm for the Simulation of Power System Transients with the OVNI Simulator by Mazana Lukic Armstrong B.Sc , Fakultet Elektrotehnike i Racunarstva, Zagreb, Croatia, 1995 M.A.Sc., The University of British Columbia, Vancouver, Canada, 2001 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 M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Doctor of Philosophy in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University Of British Columbia December, 2006 © Mazana Lukic Armstrong 2006 Abstract The existing M A T E (" M u l t i Area Thevenin Equivalent") concept implemented in the real-time simulator O V N I provides a means of partit ioning large systems of equations into subsystems connected through links. The subsystems are solved independently, with the overall solution integrated at the level of the links. M A T E partitioning enables fast simulation of power system networks in two distinct ways: first, it allows parallel processing in a multi-machine environment and second, it allows integration of different solution techniques for individual subsystems. In this thesis, we generalize the concept of M A T E wi th the new Mult i level M A T E concept, in which each subsystem becomes the basis for a new level of M A T E par-titioning. The new concept allows subsystems to be further partitioned into sub-subsystems that, for example, are of a constant and changing nature. Parti t ioning at the subsystem level leads to higher overall solution efficiency. Furthermore, power system components can be described using Mult i level M A T E by their nodal and/or branch equations on the subsystem level in O V N I . A three-phase induction machine model is an example of a power system component that is naturally described with branch equations. Mult i level M A T E also provides a convenient framework for the incorporation of controllers as well as nonlinear elements. The software implementation of O V N I written in this work is based on the M u l t i -level M A T E algorithm. Models of power system components have been implemented and tested using the newly developed concept. A test case describing a doubly-fed induction generator wind turbine system has been modelled and studied as a practical example of the new solution scheme's capabilities. i i T a b l e o f C o n t e n t s Abstract i i Table of Contents i i i List of Tables v i i List of Figures viii Acknowledgements x i i i 1 Introduction 1 1.1 Background 1 1.2 Motivat ion and Objectives 2 1.3 Thesis Organization • • 3 1.4 Thesis Contributions 4 2 A n Overview of Power System Simulation Tools and Techniques 6 2.1 Introduction 6 2.2 Techniques for Simulation of Large Power System Networks 7 2.2.1 Sparsity . 8 2.2.2 Diakoptics 10 2.2.3 Modif ied Noda l Analysis 10 2.3 Summary 12 3 Multi level M A T E 13 3.1 Introduction 13 3.2 M A T E System Part i t ioning 14 i i i Table of Contents 3.3 Mult i level M A T E Concept 16 3.3.1 Example: A n Ideal Switch Implementation wi th Mult i level M A T E 20 3.4 Computat ional Efficiency of Mult i level M A T E Part i t ioning 22 3.5 Commentary on Computational Efficiency of the Mult i level M A T E Part i t ioning 25 3.6 Solution of Nonlinearities within the Mult i level M A T E Concept . . . 25 3.6.1 Example: A Nonlinear Control Function Implementation with Mult i level M A T E 28 3.7 Mult i level M A T E Block Diagram 31 3.8 Summary 32 4 Modell ing with Multi level M A T E 34 4.1 Introduction ; 34 4.2 Model l ing of Ideal Voltage Sources 35 4.2.1 Commentary on Modell ing of Ideal Voltage Sources 36 4.3 Model l ing of Ungrounded Circuits and Subsystems 36 4.3.1 Ungrounded Circuits 37 4.3.2 Ungrounded Subsystems 42 4.4 Model l ing of Electric Machines , 46 4.4.1 Phase-Domain Induction Machine Model 47 4.4.2 Phase-Domain Synchronous Machine Mode l 64 4.4.3 Commentary on Modell ing of Electric Machines 68 4.5 Model l ing of Controllers 68 4.5.1 Example: Model l ing of a Single-Machine Infinite-Bus System in dq Coordinates wi th Mult i level M A T E 70 4.5.2 Commentary on the Simultaneous Solution of Controller Equa-tions 74 4.6 Summary 77 5 Modell ing and Simulation of a D F I G W i n d Turbine System in O V N I 78 5.1 Introduction 78 5.2 W i n d Turbine Generators 79 iv Table of Contents 5.3 Doubly Fed Induction Generator W i n d Turbines 80 5.4 D F I G Model Structure 80 5.5 Doubly Fed Induction Generator Model 81 5.5.1 Two-mass Representation of a D F I G W T G Shaft 81 5.6 Voltage Converter Mode l and Control 83 5.6.1 Stator-side Converter Model 85 5.6.2 Stator-side Converter Control 88 5.6.3 Rotor-side Converter Model 89 5.6.4 Rotor-side Converter Control 92 5.7 W i n d Turbine Model 94 5.8 M a x i m u m Power Tracking 95 5.9 Implementation of a D F I G W i n d Turbine System in O V N I 97 5.10 D F I G W i n d Turbine Test System 98 5.11 Simulation Results 99 5.11.1 Transient Response of the D F I G W i n d Turbine to a Step De-crease in W i n d Velocity 99 5.11.2 Transient Response of the D F I G W i n d Turbine to a Three-phase Short Circui t 107 5.11.3 Commentary on the Comparison of Transient Simulation Re-sults Obtained wi th E M T P - t y p e and Stability-type Tools . . . 108 5.12 Summary 108 6 Conclusion 119 6.1 Summary of Contributions 119 6.2 Recommendation for Future Research 120 Bibliography 122 Appendices 128 A Trapezoidal Integration Rule 129 B Reference Frame Transformations 131 Table of Contents B . l Clarke Transformation 131 B.2 dq Transformation 132 C D i g i t a l P a s s - B a n d F i l t e r for F i l t e r i n g S t a t o r - F l u x 134 D Test S y s t e m P a r a m e t e r s 137 D . l Parameters of the 2250 hp Induction Motor from Chapter 4.4.1 . . . 137 D.2 Parameters of the 835 M V A Synchronous Generator from Chapter 4.4.2138 D.3 Parameters of the 7.5 k W Doubly Fed Induction Generator from Chap-ter 5.10 139 v i L i s t of T a b l e s 3.1 Part i t ioning of subsystems 24 3.2 Iterations for the nonlinear controller equation 31 D . l 2250 hp induction motor parameters 137 D.2 835 M V A steam turbine generator parameters 138 D.3 7.5 k W doubly fed induction generator wind turbine system parameters 139 v i i L i s t o f F i g u r e s 2.1 Independent voltage source 11 3.1 A n example of a system that demonstrates the M A T E concept . . . . 14 3.2 General description of a link 15 3.3 A n example of a system that demonstrates the Mult i level M A T E concept 16 3.4 A n example system to demonstrate ideal switching using the Mult i level M A T E concept 20 3.5 Mul t i level M A T E vs. single-level M A T E computation speedup . . . . 24 3.6 A n example system wi th a nonlinear control function 29 3.7 Flow chart of iteration process for nonlinear controllers 32 3.8 Block diagram of Mult i level M A T E implementation 33 4.1 A n example system to demonstrate modelling of ideal voltage sources 35 4.2 A n example system to demonstrate modelling of ungrounded circuits 37 4.3 Model l ing of an ungrounded transformer circuit wi th the Mult i level M A T E approach 38 4.4 Introducing a common reference point to the ungrounded secondary side of an ideal transformer 39 4.5 Introducing a common reference point to the ungrounded network with two ungrounded ideal voltage sources 40 4.6 Introducing a common reference point to the ungrounded network with two grounded ideal voltage sources . . 42 4.7 Grounded system with an ungrounded subsystem 43 4.8 Rotor and stator circuits of an induction machine 48 4.9 Coupled inductors 49 4.10 Discrete phase domain induction machine electrical model 51 v i i i List of Figures 4.11 Electrical circuit configuration for testing the phase-domain induction machine model implemented wi th Mult i level M A T E 54 4.12 Stator phase current during free acceleration of a 2250 hp induction motor 59 4.13 Rotor phase current during free acceleration of a 2250 hp induction motor 59 4.14 Electromagnetic torque during startup of a 2250 hp induction motor . 60 4.15 Rotor speed during startup of a 2250 hp induction motor 60 4.16 Torque-speed characteristics during startup of a 2250 hp induction motor 61 4.17 Torque-speed characteristics during step changes in load torque of a 2250 hp induction motor 61 4.18 Stator phase current during step changes in load torque of a 2250 hp induction motor 62 4.19 Rotor phase current during step changes in load torque of a 2250 hp induction motor 62 4.20 Electromagnetic torque during step changes in load torque of a 2250 hp induction motor 63 4.21 Rotor speed during step changes in load torque of a 2250 hp induction motor 63 4.22 Rotor and stator circuits of a synchronous machine 64 4.23 Stator phase current during a step increase in input torque of a 835 M V A steam turbine generator . . 65 4.24 F ie ld current during a step increase in input torque of a 835 M V A steam turbine generator 65 4.25 Electromagnetic torque during a step increase in input torque of a 835 M V A steam turbine generator 66 4.26 Rotor speed during a step increase in input torque of a 835 M V A steam turbine generator . 66 4.27 Rotor angle during a step increase in input torque of a 835 M V A steam turbine generator . 67 4.28 Torque-rotor angle characteristics during a step increase in input torque of a 835 M V A steam turbine generator 67 4.29 Single-machine infinite-bus equivalent network 71 4.30 Thyristor excitation system with A V R and PSS . 71 ix List of Figures 4.31 Mult i level M A T E partitioning of the single machine infinite bus system 73 4.32 Transient response of the single-machine infinite-bus system with man-ual control: (a) terminal voltage, (b) active power, (c) rotor angle . . 75 4.33 Transient response of the single-machine infinite-bus system wi th A V R / P S S : (a) terminal voltage, (b) active power, (c) rotor angle 76 5.1 Schematic of a doubly fed induction generator wind turbine system . 81 5.2 Two-mass representation of a D F I G W T G shaft for accurate simulation of torsional shaft oscillations 82 5.3 Rotor and stator-side P W M converters of a D F I G 84 5.4 Geometric representation of the relationship between stationary abc and a(3 reference frame stator voltages, and definition of stator voltage oriented dq reference frame 86 5.5 Stator-side converter controller 88 5.6 Geometric representation of the relationship between stationary abc and a/3 reference-frame flux linkages, and definition of stator-flux-oriented dq reference frame 90 5.7 Rotor-side converter controller 93 5.8 W i n d turbine characteristics with maximum energy capture curve . . 96 5.9 D F I G wind turbine test case 99 5.10 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) rotor angular velocity, (b) electromagnetic torque, (c) rotor speed, (d) stator flux 100 5.11 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) D F I G active power, (b) D F I G reactive power, (c) machine stator active power, (d) machine stator reactive power . . . . 101 5.12 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator-side converter active power, (b) stator-side converter reactive power, (c) rotor-side converter (machine rotor) ac-tive power, (d) rotor-side converter (machine rotor) reactive power . . 102 5.13 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) D F I G stator-side R M S voltage, (b) D F I G stator-side R M S current, (c) D F I G rotor-side R M S voltage, ( d ) D F I G rotor-side R M S current 103 List of Figures 5.14 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator phase to neutral voltage, (b) stator phase current, (c) rotor phase to neutral voltage, (d) rotor phase current, (e) stator-side converter phase voltage, (f) stator-side converter phase current . 104 5.15 Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator voltages dq components, (b) stator currents dq components, (c) rotor voltages dq components, (d) rotor currents dq components, (e) stator-side converter voltages dq components, (f) stator-side converter currents dq components 105 5.16 Comparison of a transient response of the D F I G wind turbine to a step decrease in wind velocity against the experimental results from the literature: (a) rotor speed, (b) d-axis rotor current, (c) q-axis rotor current 106 5.17 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) rotor angular velocity, (b) electromagnetic torque, (c) rotor speed, (d) stator flux 109 5.18 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) D F I G active power, (b) D F I G reactive power, (c) machine stator active power, (d) machine stator reactive power 110 5.19 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) stator-side converter active power, (b) stator-side converter reactive power, (c) rotor-side converter (machine rotor) active power, (d) rotor-side converter (machine rotor) reactive power . . I l l 5.20 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) D F I G stator-side R M S voltage, (b) D F I G stator-side R M S current, (c) D F I G rotor-side R M S voltage, (d) D F I G rotor-side R M S current 112 5.21 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) stator phase to neutral voltage, (b) stator phase current, (c) rotor phase to neutral voltage, (d) rotor phase current, (e) stator-side converter phase voltage, (f) stator-side converter phase current . . . . 113 5.22 Transient response of the D F I G wind turbine D C - l i n k voltage to a three-phase short circuit 114 x i List of Figures 5.23 Transient response of the D F I G wind turbine to a three-phase short circuit: (a) stator voltages' dq components, (b) stator currents' dq components, (c) rotor voltages' dq components, (d) rotor currents' dq components, (e) stator-side converter voltages' dq components, (f) stator-side converter currents' dq components 115 5.24 Transient response of the D F I G wind turbine to a three-phase short circuit: rotor phase voltages and currents 116 5.25 Transient response of the D F I G wind turbine to a three-phase short circuit obtained wi th the stability tool T S A T : (a) D F I G stator terminal voltage, (b) D F I G rotor angular velocity, (c) D F I G active power, (d) D F I G reactive power 117 5.26 Transient response of the D F I G wind turbine to a three-phase short circuit obtained wi th the stability tool T S A T : (a) D F I G rotor dq volt-ages and currents (from top to bottom - Vdr, idr, vqr and iqr), (b) D F I G stator current 118 A . l Illustration of integration with the trapezoidal rule . . 130 C . l Circui t configuration for filtering 60 Hz component o f i/» a(t) 134 C.2 Variat ion of filtered flux-linkage magnitude as a function of frequency 135 C.3 Equivalent circuit of the parallel R L C filter discretized wi th the trape-zoidal rule 136 x i i Acknowledgements I have been extremely lucky to meet and to get to know so many wonderful people since my arrival in Canada almost eight years ago. I would sincerely like to thank every one of you who has crossed my path on this long journey. M y sincere gratitude goes to my supervisor, Prof. Jose R. Mar t i , , whose extensive knowledge and brilliant ideas have inspired this work. I thank h im for his continuous ; support and understanding, especially during the difficult times. I would also like to express my sincere gratitude to Prof. Hermann W . Dommel , whose course on electromagnetic transients opened my eyes to the world of power systems simulation. His willingness to help and his invaluable experience is very much appreciated. I would like to acknowledge Prof. Luis R. Linares for introducing me to O V N I and helping me get started wi th the tremendous task of programming the simulator. His kindness and good heart have helped me endure through many difficult times. I am extremely grateful to Dr . Prabha Kundur , who provided me wi th many opportunities at Powertech Labs, and helped me see the bigger picture of many problems we have discussed. His guidance and understanding have been invaluable. Many thanks to Dr . A l i Moshref for al l the time he spent helping me learn about stability simulations. A special place in my heart belongs to my dear friends who shared the ups and downs of graduate student life wi th me. M y special thanks goes to my dear friend Richard Rivas, for al l the discussions about work and life, as well as his honesty and kindness. I am grateful to my friend Khosro K a b i r i for always being there for me. I would like to thank Stephane Bib ian , Benedito Bonatto, Jesus Calviho, Maninder Dhaliwal , Jorge Hol lman, A w a d Ibrahim, Hee-Sang K o , Daniel Lindemeyer, Fernando Moreira, D a i Nan, K e n Wicks and many others in the old group for their friendship and care. I would also like to thank Mel iha Selak for al l the help,she provided to me and my family along the years. In the new Power Group, I would like to give special thanks to Tom de Rybe l for always being so helpful and for organizing fun events. Thanks to Hafiz, Lucy, Quanhong, Marcelo, Michael and everyone else; your thoughts and comments have been appreciated. I would like to thank my parents, M i r a and Ivan, my sister Danijela and my brother Drazen for their love and their incredible support through al l of my endeavors. They have inspired me throughout my life to excel in everything I do. I would like to thank my wonderful husband Gordon and my son David for always x i i i Acknowledgements being there for me, for their love and patience. Y o u fill my life wi th joy and happiness. Finally, I would like to acknowledge the following people and institutions whose invaluable financial support over the years have enabled me to complete this thesis work: Dr . M a r t i through a research assistantship; U B C through a University Gradu-ate Fellowship and a L i Tze Fong Memorial Fellowship; and Powertech Labs through a research assistantship and several part-time contracts. , Mazana Armstrong, October 2006 xiv C h a p t e r 1 I n t r o d u c t i o n 1.1 Background Object V i r t u a l Network Integrator (OVNI) is a simulation tool aimed at obtaining very fast, real-time solutions of large power system networks [1], [2], [3]. O V N I is based on the E M T P program (Electromagnetic Transients Program), used to simulate electromagnetic transients [4]. Quoting from [5], electromagnetic transients "occur on a scale of microseconds (e.g., ini t ia l transient recovery voltage), milliseconds (e.g., switching surges), or cycles (e.g., ferroresonance). B y nature, these phenomena are a combination of traveling-wave effects on overhead lines and cables, and of oscillations in lumped-element circuits of generators, transformers, and other devices." The term electromagnetic transients is used to distinguish them from the electromechanical oscillations of generators in transient stability studies. In the E M T P , basic lumped circuit components, such as inductors and capac-itors, are described wi th ordinary differential equations. Differential equations are discretized using the trapezoidal rule of integration (Appendix A ) to obtain alge-braic difference equations. Discrete time step sizes depend on the type of phenomena studied and the highest frequency expected in the results. They range from a few nanoseconds for very fast transients to an upper l imit of about one millisecond dic-tated by the steady-state base system frequency of 50 or 60 Hz . The E M T P methodol-ogy normally requires a three-phase system representation as opposed to the positive sequence single-phase representation used in stability programs which assume that voltages and currents can be represented at the base system frequency as phasors. The system of equations solved in the E M T P for a network wi th n nodes is formed as: [G] • [v(t)} = [h(t)} (1.1) where • [G] is a nodal conductance matrix, • [v(t)] is a vector of nodal voltages, • [h(t]] is a vector of accumulated currents representing independent current sources and history terms. 1 Chapter 1. Introduction The solution for nodal voltages is obtained by first triangularizing the [G] matrix through a Gaussian reduction procedure and then backsubstituting for the updated values of [h(t)]. Sparsity techniques 1 are used to reduce the number of operations. It is. apparent even from this very rudimentary description of the E M T P and its applications that a major drawback of E M T P simulations is their speed. The first steps in creating the new O V N I simulator were aimed at "speeding-up" the E M T P simulation [2]. The solution in [2] offered a software approach to topological decou-pling of power system networks, using the transmission line models in the E M T P , and a hardware approach to parallelize computation and achieve real-time performance. Soon after this solution, the new concept of M A T E emerged, which eliminated the need for decoupling using the transmission line model [7]. Object-oriented software implementation of O V N I was proposed in [8]. Finally, cluster-based hardware solu-tions for real-time simulations in O V N I were proposed in [9], [10]. 1.2 Motivation and Objectives As stated in the previous section, the original idea behind O V N I was to speed up the E M T P simulation, for the following two reasons: (i) to allow larger or more computationally demanding systems to be simulated wi th the E M T P at acceptable speeds in an off-line environment, and (ii) to allow real-time testing of power system equipment in a closed-loop environment. B y combining the concept of system parti-tioning (Diakoptics 2 ) wi th the concept of system reduction (Thevenin equivalents 3), the M A T E (Mul t i -Area Thevenin Equivalent) methodology has opened the door to new horizons of investigation of computational procedures and modelling techniques used in simulating large power systems. Prior to the beginning of this work O V N I existed as a special purpose simulator built for specific tasks to test for real-time performance, as described in [3]. The ultimate goal, however, was and is to develop O V N I as a general purpose power system simulator capable of handling any transient time-domain power system simulation. W i t h that in mind, an object-oriented approach to software programming of O V N I was proposed in [8]. The initiative for the work presented in this thesis started wi th the notion of building a general purpose simulator. The simulator would, along wi th the hardware solutions for speed, provide engineers with an all-in-one tool capable of performing 1 Sparsity applications in power system analysis were first introduced by N. Sato and W. F. Tinney in 1963 [6]. 2The concept of Diakoptics was first introduced by Gabriel Kron in 1952 [11]. 3 The concept of a Thevenin equivalent was first introduced by Hermann von Helmholtz in 1853 [12]. In 1883, Leon Charles Thevenin, unaware of Helmholtz's work, published the same result [13], [14]. 2 Chapter 1. Introduction the simulation of electromagnetic (EMTP- type ) transients, as well as electromechan-ical (stability) transients, wi th a wide range of available models of power system components. The model complexity would range from detailed models of, for ex-ample, frequency-dependent transmission lines and phase-domain electric machine models, to the approximate, simplified classical machine and pi-circuit transmission line representations. Justification for an integrated approach to power system analyses comes from the trends that have been emerging in the last decades: new technologies involving power electronics equipment ( H V D C transmission, F A C T S devices) and controls that make standard stability analyses much more complex. A n example of a small-scale system that embeds the complexity and challenges of today's power systems is, for example, a wind turbine generation system consisting of an induction generator, two back-to-back pulse-width-modulated ( P W M ) voltage source inverters, and complex controls to regulate the turbine's responses. Ideally, engineers studying integration of large wind farms into the power system grid would be able to choose the complexity of the models for their simulation, and perform their analyses at different levels of detail for the most accurate picture of an overall state of the system. Ultimately, a combination of detailed and approximate modelling wi th an overall goal of fast and accurate simulation results would lead to a more reliable design and operation of the system. The objectives of this work are closely related to the development of the O V N I software, including the computational methods for system partit ioning and modelling of power system components. The solutions presented in this thesis, however, are general and can be applied to many existing programs that simulate power system networks. 1.3 Thesis Organization This thesis is divided into six chapters. The first chapter briefly describes the sequence of events that led to the realization of this work, and states the objectives of this research. The second chapter provides background information on power system simulation and analysis tools in use today, then concentrates on a literature overview of the methods used for solving large systems of equations associated wi th the simulation of power system networks. We provide an overview of Sparsity, Diakoptics and Modified Nodal Analysis, and describe their applications. In Chapter 3, we first review the concept of M A T E and introduce the new concept of Mult i level M A T E . Mult i level M A T E allows multilevel parti t ioning of the system of equations, enabling the system to be divided according to its nature into changing, 3 Chapter 1. Introduction constant and nonlinear subsystems. Multilevel M A T E introduces the concept of a functional sublink that allows convenient modelling of power system components that are better described by branch equations. We discuss the computational efficiency of the multilevel approach, as well as the solution of nonlinearities. We provide examples of modelling an ideal switch and a nonlinear control function to demonstrate the concept. Using a block diagram of Multilevel M A T E , we depict the computation requirements and the flow of information from the lowest level to the highest level of system partitioning. In Chapter 4, we describe modelling with Multilevel M A T E , and present a gen-eral approach to modelling ideal voltage sources connected node-to-ground or node-to-node as sublinks. We provide a solution for matrix computations with Multilevel M A T E in the case of ungrounded circuits and ungrounded subsystems, and demon-strate it through examples. We describe and test an extensive and complete derivation of a phase-domain induction machine model implementation with Multilevel M A T E , and show the results obtained with a similarly derived phase-domain synchronous machine model. Finally, we describe and show the implementation of a controller's equations with Multilevel M A T E on a standard transients stability test case of a synchronous machine connected to an infinite bus modelled in dq coordinates. In Chapter 5 , we test our modelling methodology and computational solutions on a complex doubly fed induction generator (DFIG) wind turbine system reproduced from an experimental test case described in the literature. Following a brief introduction to wind turbine technology, we provide a detailed description of the test case and the modelling of its components, then perform the tests to validate the modelling. We then compare our results against the experimental results in the literature. In Chapter 6, the final chapter of the thesis, we state the main contributions of our research and provide recommendations for future work. 1.4 Thesis Contributions This thesis describes novel power system modelling and implementation solutions suit-able for integrated analyses of electromagnetic and stability transients in the newly developed, EMTP based, general purpose simulator OVNI. The main contributions of this thesis can be summarized as follows: • The new concept of Multilevel M A T E is proposed. It improves computational speed of the existing single-level M A T E partitioning approach, and allows mod-elling of power system components with branch equations. • An efficient framework for iterating nonlinear equations using the Multilevel M A T E concept is proposed. 4 Chapter 1. Introduction • A new technique for conditioning of ungrounded or weakly grounded circuits and subsystems is presented and demonstrated on examples. • Branch modelling of power system components using Multilevel MATE is pro-posed. Induction and synchronous machines, ideal voltage sources, ideal switches and controllers are modelled and then validated through an extensive example of a doubly fed induction generator wind turbine system. • The new concepts described in this thesis, as well as numerous new power system components models, have been implemented into a new version of the general purpose simulator, OVNI. 5 Chapter 2 An Overview of Power System Simulation Tools and Techniques 2.1 Introduction Computer simulation of power systems is a necessary tool for both planning and everyday operation of electric power networks. In the early days of power system analysis, the need for computational aids in analyzing power system networks led to the design of the first analog computer, as early as 1929 [15]. The A C network analyzer was able to determine power flows and voltages during normal and emergency conditions and switching operations. The earliest application of digital computers to power system problems dates back to the 1940s; however, it was rather l imited in scope due to the l imited capacity of the punch card computers used then. Large-scale digital computers became available in the mid-1950s, and the success of load-flow programs led to the development of programs for short-circuit and stability analyses. Phenomena analyzed wi th different simulation methods and tools range from slower widespread network events, such as voltage instability, to fast electromag-netic transients affecting networks locally due to, for example, surges from lightning strikes. Methods for power system simulation range from load-flow [16] programs that use approximate, simpler network representations, to highly sophisticated and complex modelling by the Electromagnetic Transients Program ( E M T P ) [4]. Power-flow calculations are used to determine the steady-state of the power gen-eration/transmission system for the given loading condition. The power and voltage constraints make the load-flow problem nonlinear. The numerical solution there-fore must be iterative, leading to convergency and other solution problems. Possibly thousands of load-flow methods have been developed over the years, streaming from continuous efforts to improve the original load-flow concept and address the needs of, for example, power system reliability and security. The stability of power systems was associated early on wi th rotor-angle stability of synchronous machines. More recently, voltage stability has emerged as one of the major concerns in modern power systems characterized by highly loaded transmission networks operating close to their capacity limits. Dynamic simulations aimed at 6 Chapter. 2. An Overview of Power System Simulation Tools and Techniques studying the stability of power system networks require more complex modelling than do power-flow studies. Dynamic simulations are currently an integral part of planning and operation studies of power system networks. Electromagnetic transients simulations are the most common methods for simu-lation of power system surges and fast-circuit transients. Electromagnetic transients are rarely an important factor in system planning or operating decisions, but are essential factors in the design process. These types of simulations require extensive modelling detail and small simulation time steps due to the small time constants associated wi th the studied phenomena. W i t h the ever-increasing computer power available today, the E M T P has become accepted for the solution of a wide range of transient problems. In particular, the E M T P - t y p e programs are used extensively to analyze system disturbances in power system protection design. Regardless of the simulation method used and the complexity of the associated modelling, most of the tools employed today describe power system transmission networks wi th nodal equations that are particulary suitable for digital simulation. Moreover, they a l l use similar programming techniques to efficiently handle large ma-trices. The following section discusses the historical development of these techniques and their applications. 2.2 Techniques for Simulation of Large Power System Networks Nodal analysis implements Kirchhoff's current law 4 , which states that the sum of the currents entering the node in a circuit is equal to the sum of the currents leaving the same node. A s its name suggests, nodal analysis produces nodal (node) voltages as a solution. The nodal approach has an advantage over the mesh approach derived from Kirchhoff's voltage law in that it introduces fewer variables and equations in the case of power system networks, therefore leading to smaller system matrices. The system of equations associated wi th nodal analysis has the general form: [Y] • [V] = [/] (2.1) where • [Y] is a nodal admittance matrix, • [V] is a vector of nodal voltages, • [I] is a vector of nodal currents. 4Kirchhoff's laws were first announced by Gustav Robert Kirchhoff in 1845. 7 Chapter 2. An Overview of Power System Simulation Tools and Techniques To solve this system of equations for nodal voltages, one has to determine the inverse of [Y]. For a large number of nodes, this task is extremely demanding com-putationally, especially in the past when computer power was l imited and expensive compared to today. Practical ly from the beginning of the digital-computer simulation era, power systems engineers have been developing techniques to enable faster solu-tions of large systems of equations. Two main very distinctive streams established themselves in the 1960s, Diakoptics and Sparsity Techniques. Sparsity techniques, highly favoured over Diakoptics, have developed roots in the majority of the simula-tion tools in existence today. To overcome the limitations of nodal analysis in the modelling of some system components, modified nodal analysis ( M N A ) was formulated in the mid-1970s. Be-cause of its relevance to this work, we wi l l discuss M N A in the last part of this section. 2.2.1 Sparsity In 1978, W . F . Tinney, referring to his previous papers [6] and [17] wrote in a com-mentary: "The aim of sparsity methods is to exploit the property that the coefficients of many large systems of simultaneous equations are mostly zero. Standard matrix methods take no advantage of the zeros, storing and processing them the same as non-zero numbers, and causing most, or al l , of them to become non-zero in the solu-tion process... Today, awareness of sparse matrix methods is practically inescapable. Papers, books, and symposia on sparsity abound and applications exist in almost every field. Bu t only a decade ago, when our paper was published, the concepts of sparsity were vir tual ly unknown. The advantages are so great (10 to 100 or more) and the idea is so simple, it is puzzling that it took so long to be discovered." [18] The first paper on sparsity appeared in 1957 in Management Science [19]. The first sparsity method applied to the solution of power system networks was developed by Sato and Tinney in 1961, and published in 1963 [6]. The work in these papers was not recognized or accepted after published, and remained unknown until the mid 1960s. In 1967, when the second paper was published by Tinney and Walker [17], sparsity was finally acknowledged, and has since become widely accepted and applied in power system simulation tools. The original method proposed in the 1960s remains valid today, wi th most improvements consisting of marginally more efficient algorithms. B y 1977 when a survey paper on sparse matr ix research was published [20], there were already over six hundred published references on the topic of sparsity. In nodal analysis the power system network is described by a large sparse matrix of coefficients referred to as its "admittance matrix". The task of solving the system for nodal voltages involves obtaining the inverse of the matr ix of coefficients, which is a very computationally demanding task for large power systems including 20,000 8 • Chapter 2. An Overview of Power System Simulation Tools and Techniques to 50,000 nodes. B y appropriately ordered triangular decomposition, the inverse of a sparse matrix can be expressed as a product of sparse matr ix factors, which brings an advantage in computational speed, storage requirements and accuracy. The method of sparsity follows these essential steps: (1) store and process only non-zero terms; (2) factor out the coefficient matrix into upper and lower triangular matrices; (3) order the factorization to approximately minimize new non-zero terms; and (4) use the factors to obtain the direct solutions by forward and back substitution. The most evident advantage that can be taken of in sparse systems is information storage and retrieval. This aspect alone allowed simulation of large systems with more accurate models. The computation procedures are organized to avoid operations with zeros (or ones), to minimize overhead such as accessing of arrays and to reduce the arithmetic operations involved. Sparsity algorithms are designed to preserve the system's ini t ia l sparsity as much as possible. The main drawbacks of sparsity methods are related to implementation issues. In order to achieve the advantages of sparsity, the programming scheme must store and process only non-zero elements. In a matrix, the address of each element in the computer's memory can be related to its row and column indices. B y using sparsity, in addition to the memory allocation of non-zero matrix elements, tables of index-ing information are required to identify the elements and facilitate their addressing. This procedure, if not properly planned, may diminish the time-saving advantages of sparsity methods. The most effective way of ordering the operations is performed as i f it were being done by visual inspection and manual computation. It is very important to realize that, although the word "optimal" is often used to describe ordering algorithms, no known algorithms are optimal in a global sense for general sparse matrices. A t the completion of the ordering algorithm, the exact form of the table of factors is established, and this information is recorded to guide the elimination process. During the elimination, no operation is performed that would lead to a predictable zero result, and no memory is allocated for a predictable zero element. The original matrix, the table of factors, and al l indexing tables contain only non-zero elements. Rows from the original matrix are transferred to a compact working row in which elements to the left of the diagonal are eliminated by appropriate linear combination with previously processed rows from the partially completed table of factors. When work on the row has been completed, it is added to the table of factors. For nonsym-metric matrices, the derived elements to the left of the diagonal must be stored as well. Finally, a direct solution of the system of equations is obtained by mult iplying the matrices derived from the table of factors. 9 Chapter 2. An Overview of Power System Simulation Tools and Techniques 2.2.2 Diakoptics In Diakoptics, the networks are partitioned into subnetworks that can be analyzed independently, then the subnetwork solutions are combined to obtain the solution of the entire network. In the literature, Diakoptics is also referred to as a solution of networks by tearing. The concept of Diakoptics was first proposed by Gabriel K r o n in 1953 [11], in the Journal of Applied Physics. His sixteen-page paper describes the advantages of system tearing and gives hints about a multi-machine solution of large-system matrix equations. It also points out a reduction in the computation time of matrix inversion on a single machine simply due to partitioning. His brilliant ideas, which are in fact one of the essences of the O V N I simulator, were well ahead of his time and were not understood by his contemporaries. Kron 's Diakoptics, nevertheless, was not completely forgotten. A number of pub-lications expanded on his ideas for applications in analysis of power systems, electrical networks and structural and other large-scale systems. In Kron ' s original form, D i -akoptics required the explicit computation of inverses of submatrices representing the individual subnetworks. As a result, sparsity of subnetwork equations could not be exploited. A number of works published in the 1970s, just after sparsity established it-self, investigated network analysis by tearing using sparse matr ix solution techniques. A good overview of these attempts is given in [21]. To this day, however, Diakoptics, however, has been questioned as to whether more efficient computations are possible with the use of sparsity and solving the network as a whole. The only situations when system tearing was deemed necessary involved very large networks, the equations for which could not be stored on a single computer, even though sparse matrix solution techniques were being used. Tearing was also recognized as necessary in cases of repetitive identical networks when the equations of only one such network needed to be solved and stored. In addition, Diakoptics was identified by only a few researchers as a solution to allow parallel processing and the use of latency concepts in finding the network solution. Possibly the main reason why Diakoptics never developed roots in the simulation of large power system networks was due to its ambiguous description in the many papers that could not deliver clear application examples and justify its use. Moreover, a number of papers, as stated in [20], claimed that Diakoptics tearing could never be more advantageous than sparsity techniques. 2.2.3 Modified Nodal Analysis As stated in [22], nodal analysis "treats voltage sources inefficiently and is incapable of including current-dependent elements, linear or nonlinear... Another disadvantage of the nodal method is that branch currents are not accurately or conveniently obtained 10 Chapter 2. An Overview of Power System Simulation Tools and Techniques as part of the output of the program." In this 1975 paper, the new concept of modified nodal analysis (MNA) is proposed, which allows branch equations to be incorporated into the system of nodal equations, eliminating these disadvantages. Certain circuit components, such as voltage sources, transformers and dependent sources, cannot be incorporated into nodal analysis in a natural way. In addition, certain variables, such as source or transformer currents, have to be obtained by postprocessing the results of nodal analysis. In the modified nodal analysis, all non-natural elements are removed in order to write a conventional set of nodal equations, and are then reintroduced with branch equations. • a Z« b > - - -—vw—(+3—- ~ T -: • i ^ -Figure 2.1: Independent voltage source The system of equations associated with MNA can be demonstrated through an example of an ideal voltage (Vs) and series impedance (Zs) connected between the nodes a and b in Fig. 2.1. The original set of nodal equations (2.1) is affected in two ways. First, the new constraint equation is introduced relating the two nodal voltages, Va and Vj,. Second, the new unknown current, Is, is permitted to flow. The new equation will increase the number of rows and columns in the system matrix by one. Taking into account the flow of the current and the polarity of the voltage source, as depicted in-Fig. 2.1, the following MNA system of equations is obtained: YR +1 ' va' ' la ' -1 vb = h (2.2) +1 -1 - Z s \ Is where Ia and lb are the nodal injection currents of the system. In (2.2), matrix [YR] is a reduced form of the nodal matrix [Y], excluding the contribution due to voltage sources, current controlling elements, etc. For any given circuit, the MNA matrix dimension is the sum of the number of nodes and the number of currents as outputs. In situations where the series impedance Zs is zero, conventional nodal analysis would not produce a solution because of an undefined l/Zs entry into the admittance matrix, unless one of the nodes, a or b, is reduced from the system of equations. With the modified nodal analysis, a solution can be found by employing a suitable pivoting strategy to maintain the strong diagonal of the MNA matrix, which is desirable from 11 Chapter 2. An Overview of Power System Simulation Tools and Techniques an accuracy and efficiency point of view. The M N A equations are ordered in such a way as to reduce both execution time and storage requirements. The row interchange process is applied to the node equation of each node con-nected to one or more voltage sources. The node equation is always interchanged with one of the branch equations having the smallest number of non-zeros to mini-mize fill-ins. For example, in (2.2) the bottom equation would be interchanged either wi th the node equation for Va or H , depending on their number of non-zero elements dictated by the composition of the adjacent network. Modified nodal analysis combined with sparsity techniques is commonly used in to-day's simulation tools. It removes the limitations imposed by classical nodal analysis and is appropriate for the symbolic and numeric analysis of linear electrical circuits. Adopt ing M N A results in an increased order of the system of equations. These branch equations are generally highly sparse and can be efficiently dealt wi th using modern sparsity-based computing packages. 2.3 Summary This chapter provides a brief overview of the existing methods and tools for power system simulation. Further discussion on existing techniques for efficient computation of large and sparse systems of equations, namely Sparsity and Diakoptics, is presented. Modified nodal analysis and its advantages is also presented. 12 Chapter 3 Multilevel M A T E 3.1 Introduction In this chapter, we introduce our novel concept of Mult i level M A T E . Mult i level M A T E expands the original concept of M A T E by allowing subsystems to be partitioned with sublinks. A subsystem, for example, can be partitioned into constant and changing nature sub-subsystems, where only the changing nature subsystem's matrix needs to be recalculated at every simulation time step. A n example of a power system element of changing nature is a phase-domain synchronous machine model. Mult i level M A T E also allows for convenient modelling of power system compo-nents whose function is better described by branch equations, as opposed to the nodal equations typically used in the E M T P . A n example of such an element is a switch. If the change of a system's topology due to switching is confined to branch equations, the subsystem's nodal equations are unaffected. Subsystem partit ioning and branch (functional) modelling bring a significant improvement to the simulation efficiency of the O V N I code. Mult i level M A T E provides an efficient framework for implementing nonlinearities that require an iterative procedure. Power system components that can be conve-niently modelled wi th in the Mult i level M A T E framework include switches, controllers, transformers, synchronous and induction machines, voltage-source inverters etc. The research contributions reported in this chapter include the following: • Formulation of the new Multilevel MATE concept. • Functional modelling of power system components with Multilevel MATE. • Description of the new and efficient framework for iterating nonlinear equations with Multilevel MATE. 13 Chapter 3. Multilevel MATE 3.2 M A T E System Partitioning The Mul t i -Area Thevenin Equivalent ( M A T E ) concept was first introduced by M a r t i [23] in an internal report to his research group. Reference [1] presents a detailed explanation of the algorithm M A T E that extends the main ideas of Diakoptics by recognizing that the subsystems split by the branch links can be represented by Thevenin Equivalents. To demonstrate the M A T E concept, let us consider the system in F i g . 3.1. A n y system can be partitioned into subsystems by introducing links. For Figure 3.1: A n example of a system that demonstrates the M A T E concept example, the system depicted in F i g . 3.1 is partitioned into three subsystems A, B and C, by introducing six links. The hybrid system of modified nodal equations has the following form: 'A 0 0 p' ' vA ' hA 0 B 0 Q VB hB 0 0 c r vc hc Pl 9* rt —z (3-1) where: • [A], [B] and [C] are the subsystems' admittance matrices •• \p]i [Q] a n d [r] are the subsystems' current injection (or connectivity) arrays • b*L a n d [ r t] are transpose arrays of p, q and r • [z] is a matrix of the links' Thevenin impedances 14 Chapter 3. Multilevel MATE • IAB]> [^cl a r e vectors of the subsystems' accumulated currents • [Va] is a vector of the links' Thevenin voltages • [VA], [VB], [VC] are the subsystems' nodal voltages • [ia] is a vector of the links currents. Connectivity arrays for the six links connecting the subsystems are constructed so that a particular link current is assigned a direction of flow where the "from" node is assigned a "1" in the connectivity array of the "from" subsystem, and the "to" node is assigned a in the connectivity .array of the "to" subsystem. A n y component (or branch) in the system can become a link described by its Thevenin impedance ([z]) and its Thevenin voltage ([V^]). The most general representation of a link is depicted in F i g . 3.2. Note that [z] and/or [Va] can be zero. "from" z "to" • Figure 3.2: General description of a link Parti t ioned matrices in (3.1) are manipulated to obtain the subsystems' Thevenin equivalents as seen from the l inking nodes [1]. The result of the manipulation is the M A T E system of equations, in the following form: I 0 0 a 0 / 0 6 0 0 I c 0 0 0 za where a = A b = B c = C - l P VA eA VB v c ia ea &A --= A~ lhA eB z = B~ lhB e c --= c- lhc (3 .2 ) a ea pla + qtb + rtc + z •T?eA + qteB + Ttec + Vtt and [/] is the identity matrix. The individuali ty of each subsystem is preserved in (3.2). Vector [e^] and array [a] represent, respectively, the Thevenin equivalent voltage vector and the Thevenin 15 Chapter 3. Multilevel MATE impedance array of subsystem A as seen from the l inking nodes. The solution for link currents is now independent from the solution for the nodal subsystems' voltages. The interaction of the subsystems' Thevenin equivalents is solved at the links level and returned to each subsystem in the form of injected l ink currents at the l inking nodes: ia = za ea The subsystems' nodal equations are then solved independently from each other at the subsystem level as: v A - e A - a-ia vB = eB - b • ia vc - ec - c-ia 3.3 Multilevel M A T E Concept The concept of Mult i level M A T E introduced in this thesis can be demonstrated on the same system depicted in F ig . 3.1. In this case, the system is further divided into ten sub-subsystems wi th a total of six links and twenty sublinks, as shown in F i g . 3.3. The complete system of equations for the depicted system is shown in (3.3). Matrices [A], [B] and [C] represent the conductance matrices of the corresponding subsystems, while \pA], [<7B], [rc] are the connectivity arrays of the sub-subsystems ( A i , • • • Bi, B2, • • • Ci, C2) wi thin the corresponding subsystems. Subsystem B Subsystem C Figure 3.3: A n example of a system that demonstrates the Mult i level M A T E concept 16 Chapter 3. Multilevel MATE The hybrid system of equations describing the system shown in F i g . 3 .3 has the following form: A pA ZA 0 PA 0 where [A] A, 0 0 and 0 B QB -ZB 0 [q* 0] C c 0 0 rc -zc P VA hA rj 1Asublink VAsublink q vB hB 0 1B sublink = VBsublink r vc he 0 }Csublink_ _ Vc_sublink —z i a ( 3 . 3 ) 0 0 A2 0 0 A3 [B]=-0 0 0 0 0 B2 0 0 0 0 0 Bz 0 0 0 0 0 BA 0 0 0 0 0 Bb [C] Ci 0 0 c2 PAX ' ' Pi ' ' VA\. ' ' hAX ' P A = PA2 P = P2 VA = VA2 hA = PAZ Pz VA3 hA3 for subsystem A, wi th similar expressions for subsystems B and C. Matrices [ZA], [ZB] and [zc] correspond to the sublink's Thevenin impedances, while vectors [V/Lsu&imk]) [VBsubUnk] and \Yc_subUnk] correspond to the sublink's Thevenin voltages. The dimen-sion of vectors [iAsubUnk] and [V]4_SUMtTik] i s that of the number of sublinks in subsystem A. To solve the Mult i level M A T E equations, first the M A T E concept is applied to 17 Chapter 3. Multilevel MATE each subsystem to obtain the system of equations (3.4): 17 A~l" PA 0 pAA-lpA + zA 0 0 [p* o] vA ^Asublink VB ^Bsublink vc }C sublink 0 0 q'gB-iqB + Zg 0 [q* 0] ptAA~XhA + VA_sUbiink B~lhB QBB~lhB + VejuftJinA: rcC"XhC + Vc_sublink -vn 0 0 [r* 0] P\A~XP ' B-Xq ' qBB-lq ' c-v ' rcC-xr Note, that the sublink currents for subsystems A , B and C , respectively, are: (ptAA~1pA + ZA) • iA_sublink + PAA~lp • ia •• {qBB-lqB +.zB) • iBsubiink + QBB~lq • ia •(rtcC~lrc +.zc) • icsubiink + rcC~lr • ia p\A lhA + VA_sublink -- qBB hB + VB sublink rcC-lhc + Vc sublink (3.4) (3-5) Before applying the M A T E concept to the entire system, we note that the sublink currents at the subsystem level in (3.3) do not need to be known at the system level.. The computational efficiency of the M A T E solution is best maintained if the sublink branch currents are "somehow" made invisible to the top level of the system's partitioning. Going back to Equation (3.3), we can remove the equations for the sublink currents from the system of equations (3.4) and transfer these currents to the right-hand side, as shown in (3.6). A 0 0 p vA 0 B 0 q vB 0 0 c r vc pt q* rt —z hA pA " iAsublink hB ~ qB • iB_sublink he — rC • ic.sublink (3-6) Equations (3.5) and (3.6) fully describe the system. Equat ion (3.6) closely resem-bles the system of original equations for the single-level M A T E system partitioning (3.1). Next, the M A T E equations (3.2) are.applied to (3.6) to obtain the following 18 Chapter 3. Multilevel MATE equations: I 0 0 a VA e A e Asublink 0 I 0 b VB eB &Bsublink 0 0 I c Vc ec &C sublink 0 0 0 za la ^oisublink (3.7) where e Asublink — A pA ' iAsublink &B sublink = B~1qg • IB sublink CC sublink = C~lrc • %C sublink ^•asublink ' P ^-Asublink "t" Q ^Bsublink "T" f &C'sublink Vectors [eA_subHnk], [eBsublink] and [eCsubiink] represent the contributions of the sub-link currents to the Thevenin equivalent voltages of each subsystem as seen from the linking nodes. B y substituting (3.5) into (3.8), we can eliminate the sublink currents from the overall system of equations and obtain modified Thevenin impedances and voltages. For subsystem A, for example, [a] and [e^] represent a Thevenin equivalent as seen from the l inking nodes, as if the sublinks were not present in the system. The Modified Thevenin Equivalent that includes the sublinks' contributions is represented by [CLMTE] and [eA_MTE]. The system of equations for the Modified Thevenin Equivalents is: where 7 0 0 a M T E 0 1 0 bMTE 0 0 1 CMTE 0 0 0 Za_MTE O-MTE = a — A a bMTE — b - Ab CMTE = C - A C VA &A-MTE VB Z-B-MTE v c eC-MTE eA_MTE (3.9) e-A.MTE = eA — AeA e-B-MTE = £B — AeB ec-MTE = ec - A e c (3.10) Za.MTE = PlaMTE + Q^MTE + ^CMTE + * Za-MTE = P^A-MTE + Qt£B.MTE + rteC-MTE + Va and A a = A~lpA {p\A~lpA + zA) -1 t p\-a AeA = A xpA {p\A lpA + zA) (pA • eA + VA_subHTlk) for subsystem A, wi th similar expressions for subsystems B and C. (3-11) A t every time step, each subsystem wi l l pass its modified Thevenin impedances and voltages (i.e., [CIMTE] and [eA_MTE] for subsystem A) to the links. The system of links is solved first, and the link currents are returned to their corresponding subsys-tems to obtain the solution for the subsystem voltages. In the case where subsystem 19 Chapter 3. Multilevel MATE links are introduced for the purpose of subsystem parti t ioning (partitioning sublinks), the sublink currents do not need to be solved. However, i f some of the sublinks are functional sublinks (e.g., nonlinear elements) and it is important to monitor their currents (as in the case of switches), these sublinks need to be solved at every time step. In the following example, an ideal switch implementation wi th the Multi level M A T E approach is presented. 3.3.1 Example: An Ideal Switch Implementation with Multilevel MATE Multi level M A T E allows switches to be modelled as sublinks wi th in each subsystem. Switching operations involve a change of the subsystem's topology. However, this change can be confined to terms in (3.11) that, in effect, modify the subsystem's Thevenin equivalent as seen from the connecting link nodes. F i g . 3.4 depicts a simple example system to demonstrate the principles of ideal switching on a subsystem level. Sub_subsystem A| gA12=2 S \ 1 - A A A M ^ — * / . t Sub-subsystem A 2 ^ ,(t)=l A (f) , gA34=2 S 5 A / v V ^ 4 gA40=l S / . Rlink=l CI . , gB56=4 s 6 r v V v v l — - — W \ A ilink(t) 1 s> rf) i6(t)=i A Subsystem A i—r / \ Subsystem B Figure 3.4: A n example system to demonstrate ideal switching using the Multi level M A T E concept The system is composed of two subsystems, A and B, connected through a single link. A n ideal switch is placed within subsystem A and its branch equation is treated as a sublink. The hybrid system of equations is written from (3.3) in general terms as: AX 0 PAX Px VAX tlAX 0 A2 PA2 0 P2 VA2 h<A2 PA -ZA 0 ^Asublink ^Asublink 0 B q vB hB [P* 1 o] ql —z ia -va (3.12) 20 Chapter 3. Multilevel MATE For the example system it is,written as: 3 - 2 0 0 0 - 2 2 0 0 1 0 0 2 - 2 - 1 0 0 - 2 3 0 0 1 - 1 0 0 0 [ 0 0 0 1 | o] 4 - 4 - 4 5 [ - 1 0 ] 0 Vl 1 0 V2 0 0 V3 0 1 V4 0 0 "^switch 0 - 1 : V5 ' 0 ' 0 . V* . 1 - 1 ^link 0 W i t h the switch taken out of the system, we apply M A T E as in (3.2), and obtain the subsystems' Thevenin equivalents as seen from the link nodes: za = pla + qlb + z = 1 + 1.25 + 1 = 3.25 [fi] ea = pleA + qleB + Va = 0 - 1 + 0 = - 1 [V] (3.13) where a = A7lPl 21P2 A'1 A7lh A^h Al A2 b = B-lq eB - B l h B The modified Thevenin equivalent is obtained for the case where the switch is closed from (3.10) and (3.11). B y referring to (3.12), the case where a switch in subsystem A is closed implies that the corresponding diagonal element of the [ZA] matrix becomes zero (Rswitch = 0)- For an ideal switch opening, the corresponding column of the subsystem's connectivity array [p^] and the row of \pA] are set to zero, and a diagonal of [zA] is set to one (iswitch = 0). The change of topology due to switching w i l l modify the A terms in (3.11) or, in a more general case with nonlinearities, the A terms in (3.17). For example, when the switch is closed we obtain: z*MTE = (pla - p*Aa) + 9*6 + z.= (l- 0.33) + 1.25 + 1 = 2.92 [fi] ea.MTE = (p*eA - P*Ae A ) + qleB + Va = (0 + 0.33) - 1 + 0 = -0 .67 [V] When the switch is open, the modified Thevenin equivalent is the same as (3.13). The switch current can be calculated at any time during the simulation by solving the equation for sublink currents. Model l ing a switch using the Multi level M A T E concept has several advantages over the previous approach used in O V N I [8]. The subsystem admittance matrix (subsystem topology) is not affected by switching, therefore it does not need to be 21 Chapter 3. Multilevel MATE retriangularized and/or prestored. The size of the subsystem matrix is not increased by adding branch equations that are solved simultaneously wi th nodal equations, as is the case wi th the modified nodal analysis. M N A not only increases the size of the subsystem matrix, but also mixes nodal and branch equations. Finally, wi th the new approach the switch current can be calculated at every time step, or only when needed. The case of an ideal switch implementation with Mult i level M A T E can be extended to any element that contributes to a subsystem with branch equations. Such elements include, for example, ideal transformers, controllers, average models of voltage source inverters, and so on. These branch equations are treated as sublinks in O V N I and contribute to the subsystems' Thevenin equivalents obtained from nodal equations to form Modified Thevenin equivalents ( M T E s ) . 3.4 Computational Efficiency of Multilevel M A T E Partitioning The M A T E system partit ioning is the core feature of O V N I that maps the partitioned network onto a hardware configuration of P C clusters to allow for very fast simulation speeds [9]. Each subsystem and the system of l ink (branch) equations are solved on separate P C clusters that communicate with each other at a predefined rate. High M A T E partit ioning of the network produces smaller-sized subsystems that are solved more efficiently (with a higher hardware price), but which also results in an increased number of links that have the opposite effect on the simulation speed. The size of the system of links, because of its communication requirements, is a l imit ing factor in achieving real-time simulation speeds with O V N I [10]. The a im of M A T E partitioning is to attain an optimal compromise between maximizing the computation speed of the subsystems and minimizing the communication wi th the links system by keeping its size reasonably small. Opt imizing system partit ioning is beyond the scope of this, thesis. However, papers dealing wi th the subject have been published since the early days of Diakoptics, wi th the most recent one published in conjunction with the development of O V N I [24]. Mult i level M A T E partitioning recognizes the fact that the partitioning (and its corresponding links) that may be required to separate parts of a system due to their different natures (i.e., different simulation step requirements, constant vs. changing component natures) may not in general correspond to optimal partitioning for the most efficient computational time. Moreover, the links' system size, being a constraint for real-time simulation, should not be affected by the presence of, for example, nonlinearities, controllers and switches in subsystems. The Mult i level M A T E concept offers an efficient solution for dealing with branch equations at the subsystem level 22 Chapter 3. Multilevel MATE without affecting the size of the links or the size of the subsystems' nodal equations. Mult i level M A T E introduces a new type of link at the subsystem level, referred to as sublinks. Sublinks are divided into two categories, partitioning sublinks, which are introduced to topologically parti t ion a subsystem, and functional sublinks, which come from the modelling requirements of certain power system components. Parti t ioning sublinks are introduced to increase the efficiency of the computation, for example, by dividing a subsystem into sub-subsystems of constant and changing natures. A good example of a power system component with a changing nature that needs to be updated at every simulation time step is a phase-domain induction or synchronous machine model [25]. Mult i level partitioning also allows for the use of latency [26]. For example, a power electronics device model that requires high frequency switching and therefore a very small time step, ps, can be interfaced through sublinks to the rest of the subsystem, which can be accurately simulated at a slower rate (ms). Functional sublinks are branch equations introduced by the modelling requirements of power system components. For example, a controller or a switch wi l l introduce branch equations. Other power system components that can introduce branch equations include synchronous and induction machines and transformers. Examples are given in subsequent chapters. The benefits of second-level M A T E wi l l be briefly demonstrated on the example shown in F i g . 3.1 and F ig . 3.3. Figure 3.1 depicts a system separated into three subsystems and six links based, for example, on some optimal partit ioning algorithm for real-time simulation. It is assumed that each subsystem (A, B, C) is solved on a separate P C cluster processor, and that the links' system of six links communicates with the subsystems at a predefined rate. In F i g . 3.3, each subsystem is partitioned into sub-subsystems. Because the links' system is unchanged, the difference in sim-ulation speed between the single level and multilevel approach wi l l occur inside the subsystems. We have made several assumptions for the purposes of comparing the computa-tional efficiencies of the solution of non-partitioned and partitioned subsystems. The first assumption is that subsystems A, B and C have no functional sublinks. The sublinks that are introduced are only for partitioning purposes. The second assump-tion is that a l l of the subsystems are of a changing nature and that their admittance matrices need to be updated at every time step. This situation is the most computa-tionally demanding scenario. We compare the three subsystems assuming the same number of nodes (90, 900, 1800 or 3600), and a different number of partitions and links. Based on F i g . 3.3, Table 3.1 is derived. The result of our comparison of the computational efficiency of the single-level and multilevel approaches in floating point operations (flops) are shown in F i g . 3.5. Higher partit ioning of subsystems leads to faster computational time (comparison between subsystems A, B and C in F i g . 3.5 along y-axis), as long as the sublinks are 23 Chapter 3. Multilevel MATE Subsystem Number of Sub-subsystems Number of Sublinks A 3 6 B 5 13 C 2 1 Table 3.1: Part i t ioning of subsystems - • - A - s - B - A - C 25 , : 0 ! - - r — , -90 900 1800 3600 Number of Nodes Figure 3.5: Mult i level M A T E vs. single-level M A T E computation speedup not too numerous relative to the size of the sub-subsystems. For example, in the case where subsystems A, B and C have ninety nodes each, the optimal partitioning among the three subsystems is the partitioning of subsystem A wi th three sub-subsystems of thirty nodes and six sublinks. The speedup over the single-level M A T E with no subsystem partit ioning is 2.88 times, as opposed to the speedup for partitioning of subsystems B and C of 2.35 and 2.79 times, respectively. Improvement in the computational speed over single-level M A T E is greater for large networks wi th many nodes (comparison between different number of nodes for individual subsystems in F ig . 3.5 along x-axis). Mult i level M A T E reduces the number of computations over twenty times in the case of the largest number of nodes for partitioning of subsystem B . If some of the subsystem matrices were constant by nature, their inverses would be pre-calculated prior to the simulation run time, and the Mult i level M A T E partit ioning would show an even greater increase in calculation speed. If sub-subsystems Bi, B2, B3 and 234 were constant and B5 was the only subsystem of changing nature, Mult i level vs. single-level M A T E would result in a computational speedup of over a hundred times for the case of 3600 nodes. 24 Chapter 3. Multilevel MATE In general, for subsystems comprised of both changing and constant nature ele-ments and for subsystems that are sparse, the application of the Mult i level M A T E concept wi l l always lead to a significant computational speedup over the single-level M A T E concept. If conventional sparsity techniques were applied for subsystems with changing topology, each topological change (e.g. switching) wi thin a particular region of the subsystem would require the retriangularization of the entire subsystem ma-trix. W i t h the Mult i level M A T E concept, the change of topology is confined to the sublink equations, therefore the subsystem's matrix need not be recalculated when the change of topology takes place. 3.5 Commentary on Computational Efficiency of the Multilevel M A T E Partitioning The M A T E partit ioning and its computational advantages are demonstrated in several papers related to the real-time performance of the O V N I simulator [3], [2]. In this work, the advantages of the original M A T E concept are extended to the level of a subsystem, bringing along the associated efficiency and flexibility of representation. The application of the newly developed Mult i level M A T E concept, therefore, results in further computational speedup through subsystem partit ioning and/or more efficient and convenient component modelling. The analytical bases of the proposed approach are derived here, but the extent of its applications reaches far beyond the scope of this thesis work. For example, the approach could allow components of different nature requiring different integration step sizes and/or different integration methods to be combined in the same subsystem. It could also allow the integration of phasor and phase-domain modelling. Exploi t ing these aspects would result in even more efficient algorithms for the solution of large power system networks. 3.6 Solution of Nonlinearities within the Multilevel M A T E Concept The power system network is nonlinear by nature, and modern simulation tools must allow for accurate modelling of its nonlinear behavior. Nonlinearities are associated with certain power system components; for example, an under-load tap changing ( U L T C ) transformer wi l l introduce nonlinear behavior wi th its changing turns ra-tio. The iron core an electric machine is subject to saturation, which makes the input/output characteristic of the machine nonlinear. In controllers, the limiters, the function of multiplication, etc., wi l l also result in nonlinear system behavior. Nonlinearities in general can either affect the coefficients of the system equations 25 Chapter 3. Multilevel MATE (e.g., machine saturation) or the state variables (e.g., nonlinear controllers). Nonl in-ear elements that cause the coefficients of difference equations to change with time are simpler to understand and implement. A n U L T C transformer changes its turns ratio depending on the measured voltage on its low-voltage terminals, which is reflected in the change of coefficients in the system's matrix. The change of a system's coefficients with time is referred to as the system's "changing nature". Another type of nonlinearity is described by nonlinear functions applied to state variables of the system, for example, the controller's sine function, or the multiplier. The nonlinear torque equation in the modelling of a synchronous or induction machine also belongs to this type of nonlinearity. The simultaneous solution of the system's equations and the nonlinearities that operate on state variables requires an iterative procedure. The basic principle for solving nonlinearities requiring iterations within the M u l -tilevel M A T E framework is explained in general terms in the small example of (3.14). A n example of a physical system that reflect this configuration is given at the end of this chapter. In the context of Mult i level M A T E , nonlinearities iterate the solution against the link equations alone. For the general derivation of the equations it is sufficient to consider only one nonlinear subsystem, as in (3.14). •A PA PA-nonlin ' P ' VA h A PA ~ZA 0 0 0 1 Asublink ^Asublink jP Ajnonlin ^ ~ZAjnonlin _ _ 0 _ 1> Ajnonlin = ^Ajnonlin o. B q h B [pt 0 0 ] ql —z ia -va l _ u j — — I t - — I I— _ i (3.14) The nonlinear subsystem A equations are composed of: • nodal equations described by the conductance matr ix [A], the nodal voltages [VA], and the vector of accumulated currents [/ia] • branch sublink equations described by the connectivity arrays \pA] and \pA], the sublink Thevenin impedance matrix [ZA], the sublink currents [lAsubiink] and the sublink Thevenin voltages vector [VA.5U6HTI*;] • branch nonlinear equations described by the connectivity arrays \pAjnoniin] and [PAjnonlin], matr ix [zAjnonlin], the vector of subsystem's nonlinear currents [ ^ Ajnonlin] and the nonlinear voltage vector [ V ^ . n o n ; i n ] The meaning of the nonlinear equations' matrices and vectors wi l l be clarified in the example at the end of this section. Here we apply Mult i level M A T E to subsystem A and extract a sublink equation similar to (3.5), now containing a term for the 26 Chapter 3. Multilevel MATE nonlinear currents [i Ajnonlin]'-(pAaA + zA) ^Asublink ~ i " {pAaA .nonlin) ' ^ Ajnonlin ~ l~ VA®1 ' ^ ot — PA&A "I - VA_sublink (3.15) In a procedure equivalent to.that described in Section 3.3, the sublink currents are reduced from the subsystem's equations in (3.14) and the following system of equations is obtained: I Q-MTEnonlin 0 &MTE vA 0 VAjnonlirflMTEnrnilin "1" Z Ajnonlin PA-nonlinaMTE 1 Ajnonlin 0 I b VB o 0 ptaMTE + qlb + Z la &A.MTE P A-nonline A.MTE + VA_nonnn eB Ple-A.MTE + qteB + Va where aA = A~lpA ^Ajnonlin = A PA_nonlin a = A~lp eA = A~lhA (3.16) A(lA_nonlin — CLA (pAQ,A + ZA) pA • &A_nonlin Aa = aA {pAaA + zA)~l p\ • a AeA = aA (pAaA + zA)~l • (p\ • eA + VA_subunk) (^MTEnonlin — ^ -Ajnonlin AaA_nonnn &MTE — a, — Aa &A.MTE = eA — AeA (3.17) The nonlinear currents equation is extracted from (3.16) in the following general form: nonlin ' ^ MTE-nonlin ~\~ ZAjnonlin) ' ^ Ajnonlin H~ PAjnonlin ' &MTE ' = PAjnonlin ' eAJATE + VAjnonlin (3.18) The system of equations is then reduced by moving the unknown nonlinear currents to the right-hand side to obtain the following: "1 0 &MTE 'VA' 0 1 b • VB = 0 0 P^MTE + qtb + z ZA-MTE — e Ajnonlin eB Pl {eAJATE — eAjnonlin) + q*eB + Va (3.19) 27 Chapter 3. Multilevel MATE where 6Ajnonlin — QAjnonlin ' 1 Ajnonlin From (3.19) we notice that the nonlinearities contribute only to the subsystem's Thevenin voltage.vector. A t each solution time step subsystem A assumes the last known value of vector [iAjnoniin} and calculates the subsystem's Thevenin voltage [e-A-MTE — eAjnonlin]- The link currents are calculated from: ia — Za.MTE ieoi.MTE €-ajnonlin) (3.20) where ^ajnonlin — P ' ^-Ajnonlin and [za_MTE] and [ea_MTE] are as in (3.10). The link currents [ia] are returned to subsystem A and the nonlinear equations are recalculated to produce a new vector [iAjnoniin]- Iterations are repeated unti l the nonlinear equations satisfy their solution tolerance, then the simulation proceeds to the next time step. The nonlinear equations iterate against the links system alone, and the overall solution for nonlinearities is integrated at the links level. A n example of a nonlinear function implementation is shown in the following subsection. 3.6.1 Example: A Nonlinear Control Function Implementation with Multilevel MATE The simultaneous solution of a nonlinear control function and the network requires an iterative procedure. The Mult i level M A T E approach offers a computationally efficient framework for the iterative solution of nonlinear equations and is applicable to real-time simulations. We now examine the system in F i g . 3.6 to demonstrate the basic principles of solving a nonlinear control function equation simultaneously with the network solution. Mul t ip l ica t ion of two state variables, for example, makes the control scheme in this example nonlinear. The system in question is partitioned into two subsystems, A and B, wi th subsystem A containing a nonlinear control equation that controls the voltage of node 3. The partitioned system equations can be written in general terms from (3.14) as: A Petri 0 P VA hA Petri ~zctrl 0 Ictrl — Vctrl 0 B q hB [ p* 0 0 ] ql —z (3.21) where Vctri = V\ • v2 is a nonlinear control function that corresponds to the vector VAjnonlin in (3.14). Similarly, [pctri], [p*td a n d Wtri] correspond, using more general 28 Chapter 3. Multilevel MATE li(t)=l A Subsystem A \ / Subsystem B v- , ( t)--j0^( t) <pv,(t) gA23=2 S A / \ A A — / / Rlin^lQ ' , gB45=2S J_Uy\A/N | 7 /\/V\A- 5 gA20=l S ilink(t) gA30=l S« / -t I gB5o=i s: \ / T 1 A \ / \ Figure 3.6: A n example system with a nonlinear control function notation, to \pA_noniin], \PA.noniin\ a n d [zA.noniin\, respectively. For the example in F i g . 3.6, the hybrid system of nodal and branch equations can be written as: 3 - 2 0 0 - 2 5 - 2 0 0 - 2 3 - 1 0 0 - 1 0 0 0 1 [ 0 0 1 I 0 ] 2 - 2 - 2 3 [ - 1 0 ] -1 v3 ''Ctrl V4 V5 _ llink 1 0 0 Vi • v2 0 1 0 The asterisk (*) sign denotes variables that are recalculated during the iteration pro-cess. The nonlinear equation relating the link currents wi th the nonlinear controller function is extracted as in (3.18): (Petri • actrl + Zctrl) • ^ctrl + Petri ' ° ' *a = Petri ' eA + Vc ctrl (3.22) where actri = A • Petri- For the example system, the nonlinear equation is: 0.5238 • ietri - 0.5238 • iunk = -0.1905 + vi • v2 The remaining linear part of the system is represented by the M A T E matrix in 29 Chapter 3. Multilevel MATE general terms as: " 1 0 a VA . &A &ctrl 0 1 b vB = (3.23) _ 0 0 pta + qtb + z IQ, _ p* {eA - ectrl) + o*eB + Va _ where ectTi = actri • ictri is the contribution of the nonlinear control function to the Thevenin equivalent voltage of subsystem A. For our example, the system of linear equations is: 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0.1905 0.2857 0.5238 - 1 . 5 ' - 1 3.0238 Vl V2 vA I'link 0.5238 + 0.1905 • w ' 0.2857 + 0.2857 • ictri 0.1905 + 0.5238 • ictri ' 1 " 1 -0.8095 + 0.5238 • ictri Equations (3.22) and (3.23) completely describe the example system of (3.21) with the clear separation between nonlinear and linear parts, both represented by their Thevenin equivalents as seen from the system's l inking nodes 3 and 4. To recapitulate, the two equations that iterate against each other at each time step of the simulation are in general terms: \Pctrl ' actrl + Zctrl) • ictrl = p\trl ' eA + Vctrl ~ p\trl ' a ' *o (3.24) za • ia — e a &a-ctri (3.25) where, for this case, ea_ctri — pl • actri • ictri, ox evaluated: 0.5238 • ictri = -0.1905 • v2 + 0.5238 • iHnk 3.0238 • iHnk = -0.8095 + 0.5238 • ictri A t each iteration step, subsystem A assumes the last known value of the nonlinear controller current ictr[ and calculates its contribution to ea_ctri, which modifies the Thevenin equivalent voltage. The modified Thevenin equivalent of subsystem A is then passed to the links to obtain the link current from (3.25). The calculated link current is returned to subsystem A. The nonlinear controller requests the value of the link current (ia) from the subsystem as well as the voltages v\ and v2 needed to recalculate its nonlinear equation (3.24). The new value of ictri is passed back to subsystem A in order to recalculate its modified Thevenin equivalent voltage. The modified Thevenin equivalent voltage of subsystem A is again passed to the links, and the new value of the link current is obtained. Iterations are repeated unti l the 30 Chapter 3. Multilevel MATE nonlinear controller satisfies its solution tolerance and raises a flag for the simulation to proceed to the next time step. The iteration procedure is depicted as a flow chart in F i g . 3.7. For our example, the nonlinear iterations.produce the following output: i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 llM -0.2677 -0.3082 -0.3295 -0.3401 -0.3453 -0.3478 -0.3489 -0.3495 -0.3498 -0.3499 -0.3500 -0.3500 -0.3500 -0.3500 Vi 0.5748 0.5380 0.5186 0.5090 0.5043 0.5020 0,5009 0.5004 0.5002 0.5001 0.5000 0.5000 0.5000 0.5000 Va 0.3622 0.3069 0.2779 0.2634 0.2564 0.2530 0.2514 0.2506 0.2503 0.2501 0.2500 0.2500 0.2500 0.2500 -0.2340 -0.3567 -0.4180 -0.4478 -0.4621 -0.4690 -0.4722 -0.4737 -0.4745 -0.4748 -0.4750 -0.4750 -0.4751 -0.4751 Table 3.2: Iterations for the nonlinear controller equation The nonlinear controller iterates the solution for its current against the link equa-tions alone. If there were more nonlinear controllers in the system, each one would perform its iterating procedure in the same manner, wi th the overall solution for nonlinearities being integrated at the links level. For slow-changing control variables in power system networks (e.g., R M S values of system voltages and currents), the algorithm automatically recognizes that the iteration procedure is not needed. In other words, the links solution is based on the present time step of the subsystems' Thevenin equivalents and the previous time step of the nonlinear currents. The values of the present time step nonlinear controller currents are calculated from the nonlinear equations and are used to obtain, the present time step nodal subsystem voltages. The iterative procedure that we used for solving nonlinearities corresponds to the fixed point iteration method. The main reason we chose this method is its simplicity and generality. The fixed point iteration method requires less computation time at each iteration step than the more complex Newton type iteration methods. Also, because it does not require calculation of the Jacobian matrix, the method is more general for the implementation of user-built nonlinear controller equations. It should be noted that even though the convergence of the Newton type methods could be faster, the applicability of these methods to solving nonlinear systems within the Mult i level M A T E concept needs to be investigated in the future. 3.7 Multilevel M A T E Block Diagram Fig . 3.8 shows a block diagram of the Mult i level M A T E implementation. The diagram describes the computational requirements and the flow of information from the lowest level (single element) to the highest level (subsystem's Thevenin equivalent and the links). 31 Chapter 3. Multilevel MATE t= 0, At, 2 At... T,„. i-i-I. Calculate link'currents la. ca y-a ca'_arl J Calculate controller's voltage(s) ;,(••) • • Solve nonlinear controller equation yes • Solve for subsystem's noiLiI voltiiucs Figure 3.7: F low chart of iteration process for nonlinear controllers 3.8 Summary This chapter introduces the new concept of Mult i level M A T E , which is suitable for implementation in the real-time simulator O V N I . The advantages of subsystem par-ti t ioning are presented in terms of computational speedup, and the advantages of modelling wi th branch equations are demonstrated through a switch implementation and a nonlinear control function implementation. 32 Chapter 3. Multilevel MATE ELEMENT IN SUB-SUBSYSTEM A, -5 -S i nut SUB-SUBSYSTEM /4, Knows: Calculates: A"' = A' • a, = A;1 • />, = A' • K SUBSYSTEM A + SUBLINKS Knows: PA> <M' ^A^wblink^ eA = Calculates: «A1 <v noniin , a a2 . a ^ = AA2 . „ » „ . a 3 . ."AK non/in _ P'A O - aA_aontiR eA+VAMM) Aa = a A A , , A^ rt = aA ' *A_\ubImk ^^A_nonlin ^A ^A„i\onl'm A ^ MTE ^A^iionlin f MTE yes - A a , "ME = " - A a eA_WE =eA ~^ eA A_ntiulin MTE_lumlin I M^non/i 'rt_nnn/rn ( P Al_iumlin aAfTE ^ iiontm ^4_iwirf«i J (rt rt itwrf/n c rt _ MTE. iT^rt ^ «onto >: i rt i mwttfn r*AfT£ . fl Calculate nonlinear node?~(nn) voltages and sublink branch currents: vA (nn) = ( e x _ * m : ( n n ) - ^ _ n o n „ n ( n n ) ) - f l w n r ( w i ) • / „ ^ rt ^ mMii«i;rt £ atm/iit' P aMTE> P {eA_MTE eA_nonlin) ZA_ NONLINEAR ELEMENTS ' Nonlinear • equations SYSTEM + LINKS Knows: Z, Va ,p aMJE, q bUJE, r C M T E , p {eAmE ~eA_nonlin)' 1 [eB_MTE ~ eB_nonlin)' r [eC_MTE ~ eC _nonlin) Calculates: Z<z_MTE = P'aMTE +<l'bMTE + r'' CMTE + * ea.UIE = P' {eA_MTE ~ eA_nonlin) +V (eB_ MTE CB_ nonlin J M r e eC _ nonlin ) + V a 'a - <-a_MTE ca_MTE Figure 3.8: Block diagram of Mult i level M A T E implementation 33 Chapter 4 Modelling with Multilevel M A T E 4.1 Introduction The Mult i level M A T E concept offers many advantages in the modelling of power system networks and their components, as wi l l be shown in this chapter. Some of the advantages are associated wi th higher computational efficiency of the models compared wi th other current modelling techniques. A h example demonstrating this efficiency is the model of a switch described in the preceding chapter. Multi level M A T E also provides a convenient and very general approach to the modelling of ideal voltage sources, whether grounded or ungrounded, by the means of sublinks. Other advantages of Mult i level M A T E are associated wi th modelling accuracy. This is especially evident in cases of modelling of ungrounded or weakly grounded circuits and subsystems. The nodal analysis typically used in the E M T P type of solution may suffer ill-conditioning problems when weak connections or no connec-tions to the common reference point (ground) exist. This is reflected in the system's admittance matr ix becoming singular. A similar problem arises in the modelling of networks that are grounded only through voltage sources in the context of the M A T E solution framework in the O V N I simulator. In this chapter, we present a new ap-proach to achieving good conditioning even under these l imi t ing situations. It solves the problem of ill-conditioning for many power system components. Finally, we present the Mult i level M A T E - b a s e d modelling of power system compo-nents used in O V N I , which offers a simple, more general and more intuitive modelling approach. The research contributions reported in this chapter include the following: • Description of a new technique for conditioning of ungrounded or weakly grounded circuits and subsystems. • Implementation of power system components such as ideal voltage sources, phase-domain induction and synchronous machines, and controllers with functional sublinks. 34 Chapter 4. Modelling with Multilevel MATE 4.2 Modelling of Ideal Voltage Sources W i t h i n the concept of Mult i level M A T E , ideal voltage sources are naturally modelled as functional sublinks. This approach has several advantages, the most important one being generality. Whether an ideal voltage source is connected to the ground or between the nodes in the circuit, its equation is constructed and solved in exactly the same manner. Also, the choice of calculating the voltage source current (and therefore the power supplied to the network) is inherently available and remains with the user. A simple example to demonstrate ideal voltage source modelling is shown in F ig . 4.1. v,=2 v + g!2=l S - A A A A -v2=i v "g2o=2S g34=l S - A / V V ^ • g4o=2 S Figure 4.1: A n example system to demonstrate modelling of ideal voltage sources The circuit contains two voltage sources, grounded V\ connected between node "1" and ground and ungrounded V2 connected between nodes "2" and "3". The system of equations wi th voltage sources modelled as sublinks can be written as 012 -012 0 0 - 1 0 0 -012 012 + 020 0 0 0 - 1 v2 0 0 0 034 -034 0 +1 0 0. 0 -034 034 + 040 0 0 0 - 1 0 0 0 0 0 ii -Vi 0 - 1 + 1 0 0 0 i-2 -v2 (4.1) where (4.1) corresponds to the M A T E subsystem matrices and vectors described as A PA VA hA PA - Z A . 1 Asublink ^Asublink (4.2) Referring to the procedure of solving the system of hybrid equations using the Mult i level M A T E concept described in Section 3.3, the solution for the system's 35 Chapter 4. Modelling with Multilevel MATE nodal voltages and ideal voltage source currents can be found. Voltage sources wi l l contribute to the Thevenin equivalent voltage vector CA.MTE of the subsystem to obtain: iAsubiink = (PAA~XPA + ZA)'1 • (PAA-^A + VA_Subiink) = [ 1-2727 0.1818 f vA = eA.MTE = A-lpA; iAsubiink ='[ 2.0000 0.7273 -0 .2727 -0.0909 ] 4.2.1 Commentary on Modell ing of Ideal Voltage Sources B y examining the hybrid system of equations (4.1) one can observe that in the tra-ditional E M T P formulation, the system of nodal equations would be smaller. The reason for this is that in the E M T P ideal voltage sources are internally converted into Norton equivalents. To be more specific, an ideal voltage source would "borrow" the series impedance from the network, converting it in essence into its internal source impedance. In O V N I , we rely on object-oriented design, where the identity of each component must be preserved. Therefore, the internal impedance of a voltage source is its property and does not depend on the connected network. Regardless of whether a voltage source has an internal impedance or not, it w i l l contribute to the hybrid system of equations wi th a branch equation (functional sublink). If a source is non-ideal, the internal nodal voltage (e.g., node "1" in F ig . 4.1) would not be calculated, which corresponds to the situation in the E M T P (there is no increase in the number of nodal equations over the E M T P approach). If a voltage source is ideal, its internal impedance wi l l be set to zero and the system of equations wi l l have the form as in (4.1). In both cases, the calculation of the source current is optional. To recapitulate, the user has the flexibility of specifying a voltage source as ideal, in which case its nodal voltage would be included in the solution, or non-ideal, in which case its internal node would be non-existent to the surrounding network, resulting in a smaller system of nodal equations. This approach is more general than that implemented in the E M T P and is especially useful in O V N I when, for example, a subsystem (at any level of partitioning) consists of voltage sources only (ideal or non ideal), that simply represent the Thevenin equivalent of the corresponding subsystem. 4.3 Modelling of Ungrounded Circuits and Subsystems Modell ing of ungrounded circuits and subsystems using nodal equations results in singularity of their conductance matrices in O V N I . To overcome this problem, dif-ferent solution alternatives are presented for ungrounded circuits and ungrounded subsystems. In this thesis the term ungrounded circuit refers to a circuit configura-36 Chapter 4. Modelling with Multilevel MATE tion where the reference point connection (ground) does not exist. A n example of such a configuration is an ungrounded secondary side of a power system transformer. The term ungrounded subsystem refers to a subsystem network that is connected to the ground through links only. In other words, the reference point is unknown when solving for a subsystem's network, but becomes known once the solution for the links is found. A n example of this configuration would be, for example, an ungrounded three phase induction machine connected through links to a grounded network. Sub-systems and circuits are also considered to be ungrounded i f grounded through ideal voltage sources only. 4.3.1 Ungrounded Circuits To demonstrate the proposed solution for modelling ungrounded circuits, a simple network configuration containing a single-phase ideal transformer is shown in F ig . 4.2. The system has four nodes, two of which ("3" and "4") belong to the ungrounded circuit. 1 gi2=2 S 2 — V W — — 1:10 I,=2-cos(cot) (Q^ ) ^ g l 0 = l S g34=l S Figure 4.2: A n example system to demonstrate modelling of ungrounded circuits The ideal transformer model introduces a branch equation (functional sublink) describing either one of the currents ip or is that are related to each other with the transformer's turns ratio n = ip/is. The hybrid system of equations for this particular case can be written in general terms as: 010 + 012 -012 0 0 0 -012 012 0 0 n 0 0 034 -034 - 1 0 0 -034 034 +1 0 n - 1 + 1 0 V3 h 0 0 _0_ 0 (4.3) The system's conductance matrix in (4.3) is singular, and therefore the system cannot be solved using the M A T E concept. Let us first comment on the cause of 37 Chapter 4. Modelling with Multilevel MATE the matrix singularity. Nodal voltages v3 and v 4 can only be calculated relative to a common reference point. A s the common reference point does not exist on the secondary side of the transformer, there are an infinite number of solutions for these two voltages that satisfy the ideal transformer relationship v3 — v± = n • v2- Unless the common reference point is introduced artificially, the solution for nodal voltages cannot be obtained. Now when the problem and its solution have been identified, the common refer-ence point needs to be introduced without affecting the solution of the system. One approach is to ground one of the nodes (e.g., node "4") on the transformer's secondary side. In this way the voltage V4 becomes known and the system of nodal equations has a unique solution. This approach, however, affects the system's solution in the case of unbalanced network conditions (e.g., line-to-ground fault) when simulating the operation of a three-phase transformer, by allowing the flow of a zero sequence current through the ground. To maintain the approach general enough for any num-ber of transformer windings and winding connections, the approach shown in F ig . 4.3 is proposed. Figure 4.3: Model l ing of an ungrounded transformer circuit wi th the Mult i level M A T E approach Shunt impedances zcsn = 10, can be inserted from nodes "3" and "4" to the common reference point, in this case chosen to be the ground potential of 0 V . To compensate for the extra current that the branch voltage vs = u 3 — v± has to supply to the shunt impedances, voltage-dependent current sources of vs/2 are added to the circuit at nodes "3" and "4". In this work, the shunt impedances, whose purpose is to introduce a common reference point in ungrounded circuits, are referred to as 38 Chapter 4. Modelling with Multilevel MATE the compensating shunt impedances. The circuit of F i g . 4.3 can also be depicted in a more compact representation as shown in F ig . 4.4. vS © v s /2 (J) g34=l S Figure 4.4: Introducing a common reference point to the ungrounded secondary side of an ideal' transformer The system of equations describing the secondary side network configuration shown in F i g . 4.3 or F i g . 4.4 can be written as 010 + 012 -012 0 0. 0 " -012 012 0 0 n 0 -n/2 034 + Vcsh -034 - 1 0 n/2 -034 034 + Vcsh +1 0 n - 1 + 1 0 Vl Vi h 0 0 _0_ 0 (4.4) A t t = 0 of the simulation, the circuit of F i g . 4.2 wi th the artificially introduced reference point on the secondary side of the transformer w i l l produce the following solution for its currents and voltages: is = iAsublink = (PA • aA + ZA) 1 - (PA • eA + VA.sublink) = 0.13245 A Vl V2 V3 Vi — VA — e A aA ' 1 Asublink — 0.67550 ' 0.013245 0.066225 -0.066255 V where (4.4) corresponds to the M A T E matrices and vectors described in (4.2) and i-i„ aA = A vpA IA eA = A~lhj 39 Chapter 4. Modelling with Multilevel MATE Note that with this approach, the nodal solution can be obtained even when the secondary side of the transformer is unloaded, by setting the conductance between nodes "3" and "4" to zero (# 3 4 = 0) and changing the sublink equation to is = 0 in (4.4). The correctness of the solution can be checked against that obtained using the branch voltage equation vs-g3i~ is = 0, instead of the two equations for nodal voltages ^3 and Vi, in the secondary transformer circuit. Equation (4.3) then becomes: " #io + 9u -9\2 0 0 1 T vx 1 [" J i " - # i 2 #12 0 n v2 _ 0 0 0 #34 - 1 vs 0 . °~ " - 1 | (T~J [ is J [ 0 with the solution is = 0.1325 A [Vl v2 v s ] T = [ 0.6755 0.0132 0.1325 ]T V 4.3.1.1 Ungrounded Ideal Voltage Sources in Ungrounded Circuits In a more general case, an ungrounded circuit that contains ungrounded (branch) ideal voltage sources can be artificially grounded in the same manner as shown in the case of the ungrounded transformer circuit in F ig . 4.4. However, there can be only one artificially introduced reference point in the ungrounded circuit, meaning that the compensating shunt impedances can be added to only one ungrounded voltage source. To demonstrate this, a simple network wi th two ungrounded ideal voltage sources is considered in F i g . 4.5. _ | i j g i2=l S 2 g 2 3=2 S 3 I i l • V \ A A - ^ - A A A A ^ J~s i S i n V , = 5 V Q ) V l / 2 0 L.,,, 4 g4 5=3 S ^ g 5 6=4 S 6 Figure 4.5: Introducing a common reference point to the ungrounded network with two ungrounded ideal voltage sources B y adding no compensating shunt impedances (shown dashed in F ig . 4.5) to the ungrounded voltage source V i , the system's conductance matr ix is non-invertible and 40 Chapter 4. Modelling with Multilevel MATE the solution for nodal voltages cannot be found. B y introducing the common reference point, the system's equations become 012 + Vcsh -012 0 0 0 G - 1 0 " ~Vl ' r & 1 2 -012 012 + 023 -023 0 0 0 0 0 0 0 -023 023 0 0 0 0 - 1 vz 0 0 0 0 045 + Vcsh -045 0 +1 0 V4, _ * i 2 0 0 0 -045 045 + 056 -056 0 0 v5 0 0 0 0 0 -056 056 0 + 1 ' 0 - 1 0 0 +1 0 0 0 0 h -vx 0 0 - 1 0 0 + 1 . 0 0 iz - V 3 and can now be solved, wi th the following solution: [ii iz f = [ - 2 . 4 2.4 ]T A [ vx v2 v3 V i v5 ve]T=[2.5 4.9 6 . 1 - 2 . 5 - 3 . 3 - 3 . 9 ]T V The current flowing through the circuit ix = — i 3 can also be calculated from Ki rch -hoff's Voltage Law to verify the solution of the ungrounded network VI-V3 923 934 fl45 556 -2.4 A If a second reference point were introduced through the compensating shunt impedances at voltage source V^, the obtained system solution would be incorrect. For the exam-ple system it would give the solution ix — —2.5704 A. 4.3.1.2 Grounded Ideal Voltage Sources in Ungrounded Circuits Let us now consider the system in F ig . 4.6, which has two grounded ideal voltage sources. The nodal conductance matrix of the system is singular, unless a shunt impedance is introduced. W i t h a similar approach as for the ungrounded voltage sources previously discussed, the shunt impedance is inserted at the voltage source node, and the extra current supplied to the newly added shunt element is compensated for by a voltage-dependent current source. Note that in the case wi th grounded voltage sources, the circuit is anchored to a fixed reference point. Therefore, inserting compensating shunt elements at more than one grounded voltage source node wi l l not modify the overall network solution (but is, of course, unnecessary). 41 Chapter 4. Modelling with Multilevel MATE 1 g 1 2=l S 2 g23=2 S 3 i 3 f V,=5 V ( t ) v, ri) 11 n in>(£ v3 rt) v 3 = i o v Figure 4.6: Introducing a common reference point to the ungrounded network with two grounded ideal voltage sources For the system in F i g . 4.6, the hybrid system of equations can be written as 012 + Vcsh -012 0 - 1 0 -012 012 + 023 -023 0 0 0 -023 023 + Vcsh 0 - 1 - 1 0 0 0 0 0 0 - 1 0 0 V i v2 0 V3 = V3 ii -V h -V3 with the solution [ i ! i3 ]T = [ -3.3333 3.3333 ]T A [ vx v2 v3 ]T = [ 5.0 8.3 10.0 ]T V 4.3.2 Ungrounded Subsystems In the conventional M A T E approach, i f a subsystem is not connected to a common reference point (ground), its conductance matrix cannot be inverted, and therefore the Thevenin equivalent of the subsystem cannot be found. Unlike the case of ungrounded networks, the connection to the reference point does exist outside of the ungrounded subsystem, and there is a unique solution for the network's nodal voltages. Ungrounded subsystems can also be thought of as subsystems wi th a reference point that is different from the ground. The nodal solution for the subsystem volt-ages can only be obtained i f the reference voltage can somehow be calculated prior to obtaining the solution for its nodal voltages. W i t h this in mind, we derive the follow-ing example to demonstrate how ungrounded subsystems can be handled in O V N I . The ungrounded wye-connected three-phase resistive load in F i g . 4.7 represents a good example of an ungrounded subsystem connected to a grounded network. The voltage sources on the grounded side of the network are chosen so that the system is perfectly balanced. The simple configuration of Subsystem B can be considered as the Thevenin equivalent of a more complex subsystem, as seen from the l inking 42 Chapter 4. Modelling with Multilevel MATE Ungrounded Subsystem A gi4=2S A A A A g24=2S A A A A g34=2S A A A A ra_iink—l & - M A / ^ la_link\ fb_link=l & rc link=l - A A A A link Grounded Subsystem B g58=2S 8 Va=2V A A A / u - - 0 - n g6P=21S' 9 Vb=-1 V A A A / ^ — O -A A A / ^ - O ^ Figure 4.7: Grounded system with an ungrounded subsystem nodes. Let us now find the Thevenin equivalent of the ungrounded subsystem A . The Thevenin equivalent of any ungrounded subsystem as seen from the l inking nodes (in this case nodes 1, 2 and 3) can be found by choosing any one of the nodes to be a subsystem's reference point. In the example, we arbitrarily choose node 3 to be the reference point for the ungrounded subsystem. Subsystem A ' s equations are derived with respect to the new reference point u 3 r e / as follows: 014 0 - 0 1 4 Vl "-1 0 0" ^aJink "if 0 0 024 - 0 2 4 V2 + 0 - 1 0 lbJink 0 — 0 . - 0 1 4 - 0 2 4 014 + 024 + 034. 0 0 0 J'cJink_ 0 . - 0 3 4 . ' Vzref (4.5) This equation can be written in more general terms as: A • vA + p • ia = h A - hA_ref using the familiar notation for M A T E partitioning, and introducing the new term hA_ref that corresponds to the contribution of the reference voltage to the subsystem's vector of accumulated currents: hAjref = PA-Vrl,f • Vref B y mult iplying the system of equations with A'1, we obtain the following expression: (4.6) vA + a • ia = eA -- 1 - 1 - 1 Jref 43 Chapter 4. Modelling with Multilevel MATE From (4.6) we see that the reference voltage increases al l nodal voltages in the un-grounded subsystem by the same amount. The nodal voltages in subsystem A are therefore referenced to the voltage of node 3. Therefore, the Thevenin equivalent as seen from the l inking nodes 1, 2 and 3 is: 1 0.5 0 0.5 0 0 1 0 0 2.5 0.5 0 za = pta + qtb + z = 0.5 1 0 + 0 0.5 0 + 0 1 0 = 0.5 2.5 0 0 0 0 0 0 0.5 0 0 1 0 0 1.5 ea = P^A + QteB + VA = 0 2 0 2 0 + - 1 + 0 = - 1 0 - 1 0 - 1 (4.7) The reference voltage of subsystem A is calculated at the links level and is uniquely defined by the relationship for link currents that is valid for ungrounded subsystems only: the sum of link currents coming into the subsystem has to be equal to the sum of the link currents coming out of the subsystem. If we combine this relationship with the links system of equations representing the system's Thevenin equivalent as seen from the l inking nodes, we wi l l obtain, for the example system, the following hybrid system of equations: " 2.5 0.5 0 +1 1 I" iaMnk 1 I" 2 0.5 2.5 0 +1 _ ibjink = - 1 0 0 1-5 +1 JcJjnk - 1 _ + i + i + i | o J [ v3 j [ o which can be written, using the more general M A T E notation, as: Pvref %a CQ, Pvref ZVref _ . Vref . . Kef . (4.8) B y applying the M A T E inversion concept to these modified link equations, ref-erence voltage equation can be extracted and solved independently from the link currents. The connectivity array pVref that resembles the p and q connectivity arrays on the subsystem level can be calculated from PVref = P^^PA-Vref (4-9) and P l r e f = \Pvref]T (4-10) 44 Chapter 4. Modelling with Multilevel MATE For completeness, the equation expressions to calculate zVref and hVref are given as (4.11) zvre f — PA.vref A PA.vref zAvref Kef =PA.vrefA~lhA-VAvref (4.12) where the terms zAvreS and VAVref are defined by a general expression for subsystem A reference voltages as PAvTefV.A + zAVrefvref = VAVref (4.13) The reference voltage equation for the ungrounded subsystem can be extracted and solved as follows v3r - 1 r 1 - 1 Za PVref ' ZVre >) •\P Vref za "2.5 0.5 0 " - 1 ~+l + 1 +1] 0.5 2.5 0 +1 _ 0 0 1.5 _+l_ "2.5 0.5 0 " - 1 ' 2 ' + 1 +1] 0.5 2.5 0 - 1 0 0 1.5 - 1 - 1 -[o] (4.14) -0.25 V The link currents can then be calculated from (4.8), taking into account the contri-bution of the reference voltage to the Thevenin equivalent as: taJink IbJink i(x Za • (e Q Pvref ' ^Ve/) 1 -0.5. -0.5 (4-15) The links system passes the values of the link currents back to subsystems A and B and, in addition, the value of the reference voltage to the ungrounded subsystem A . Finally, the solution for nodal voltages can be calculated for subsystems A and B as: VA = eA ~ A lpAvref • Vref vB = eB - b • ia = 0.5 -0.25 0 V 1.5 -0.75 -0.75 (4.16) V In a more general approach, the reference node does not need to be reduced from the ungrounded subsystem equations but can be duplicated in order to keep the size 45 Chapter 4. Modelling with Multilevel MATE and structure of the system of nodal equations untouched. The effect of duplicating the node appears in the equations as adding a conductance of one Siemens between the original node 3 and its reference duplicate 3 re / . However, since vs — v3ref there wi l l be no current flowing between nodes 3 and 3ref, meaning that the newly introduced reference point does not affect the solution of the system. Such a case (4.5) wi l l have the following form instead: " #u 0 0 0 924 0 0 0 ^ 4 + 1 —#14 ~924 — 934 Even though this more general approach is the one adopted for O V N I , the approach with the reduced reference node provides a more intuitive understanding of what, for example, pvref represents. The bottom line is that both approaches produce identical systems of equations (4.8). Note that the choice of reference point for the ungrounded subsystems is arbitrary. In fact, in the particular case examined here, the choice of the neutral point of a three-phase resistor's wye connection (node 4) would have been the more natural choice. Under balanced conditions, as was the case here, the voltage of the ungrounded neutral point is equal to zero. As wi l l be shown in the next section, a model of a three-phase electric machine wi th ungrounded stator windings falls into the category of ungrounded subsystems. 4.4 Modelling of Electric Machines In this section, we describe the modelling of three-phase induction and synchronous machines suitable for transient analysis with the Mult i level M A T E concept. Electric machines are the most complex electric power components and therefore the most difficult to accurately model for detailed simulation analyses. Over the years, certain models have been accepted as a good machine representation for transient stability analysis, as well as for the studies using the electromagnetic transients programs ( E M T P ) [27],[28]. Bo th approaches adopt machine modelling in dqO coordinates, which leads to less demanding computational requirements. The dqO modelling is especially convenient for the implementation into transient stability programs, where power system net--#14 Vl -#24 -#34 v3 #14 + #24 + #34. + " - 1 0 0 " r • "i "0" " 0 " + 0. - 1 o. iadink 0 0 0 0 - 1 ibdink — 0 - 1 0 0 0 J'cdink_ 0 0 Vsref (4.17) 46 Chapter 4. Modelling with Multilevel MATE works are represented wi th their single-phase equivalents. The situation when simu-lating a dqO machine model wi th the E M T P - t y p e of programs is less favourable, as the machine model has to be interfaced with the three-phase detailed network represen-tation. Also, since the E M T P - t y p e of programs are generally used for a very accurate simulation of locally spread network phenomena, the dqO model is l imited in its repre-sentation of the internal machine structure (internal faults, saturation effects). W i t h this l imitat ion in mind, phase-domain synchronous and induction machine models were proposed in [25], [29]. The phase-domain models of three-phase induction and synchronous machines of [25], [29] have been implemented in O V N I exploiting the Mult i level M A T E con-cept. The detailed description of the derivation of the phase-domain induction ma-chine model is presented and its implementation wi th Mult i level M A T E is explained. Because of the similarities in modelling of the induction and synchronous machine phase-domain models that wi l l be outlined in the corresponding section, only a brief description of the implementation of the synchronous machine phase-domain model with the Mult i level M A T E is given here. The implementations of the induction and synchronous machine phase-domain models were tested and successfully compared against results from available literature. 4.4.1 Phase-Domain Induction Machine Model 4.4.1.1 Electrical Equations of an Induction Machine Stator and rotor electrical circuits of the induction machine are, in essence, mutually coupled three-phase inductances, where the stator circuit is stationary and the rotor circuit rotates at an angular velocity ujr wi th respect to the stator, as shown in F ig . 4.8. 5 From fundamental circuit theory, the voltages across the windings (phases) of an induction machine can be expressed for the stator as: ( , = IMh. 4. /j> j a — dt lhal,a eb = ^ + Rbib (4.18) • g _ 4 , D A c Bt Lc°c and for the rotor as: eB ec 5 Figure courtesy of P. Kundur: Power System Stability and Control = ^ + RAiA = ^ + RBiB (4.19) = d-^ + Rdc 47 Chapter 4. Modelling with Multilevel MATE Figure 4.8: Rotor and stator circuits of an induction machine where the flux linkages in the stator phase a and .the rotor phase A are given by: ipa = Laaia + Labib 4- Lacic + LaAiA + LaBiB + LaCic fpA = LaAia + LbAib + LcAic + LAAiA + L A B i B + LAcic (4.20) This can be written in a matrix form for all stator (a, b, c) and rotor (A, B, C) phases as: ~1pa i>A lf>B where Laa Lab Lac LaA^r) -^afi(^r) Lac(9r) Lab Lbb Lbc LbA(6r) LbB(@r) Lbci^r) Lbc LCc LcA(@r) LcB(9r) Lcc{6r) LaA{0r) LbA(6r) LcA(6r) LAA LAB LAc LO.B{0T) LbB{9r) LcB{Qr) LAB LBB LBC Lac(6r) Lbc(6r) Lcc(0r) LAc LBC Lcc ia k ic iA iB Laa = Lbb = Lcc = Lis + Lms Lab = Lac — Lbc = ^ Lms LaA = LbB = LcC = Lsr cos(6r) LAB - LbC - LcA - Lsr cos(8r + 120°) Lac = LbA = LcB — Lsr cos(6>r - 120°) LAA = LBB = Lcc — Lir + Lmr LAB — LAc = LBc = —\Lmr (4.21) and • Lis and Lir are leakage inductances of stator and rotor windings respectively • Lms is a magnetizing inductance of the stator winding 48 Chapter 4. Modelling with Multilevel MATE • Lmr is a magnetizing inductance of the rotor winding • Lsr is the amplitude of the mutual inductance between stator and rotor windings From (4.21) it can be noted that the mutual inductances between the rotor and stator windings depend on the rotor angular position 6r(t) and are, therefore, time-dependent. i, + e> —>— h — » > — + e2 Figure 4.9: Coupled inductors In the general case of two coupled inductors, as depicted in F i g . 4.9, the voltage equations C l(t) = ^ e 2 (i) - *m with the flux linkages (4.22) dt ^(t) = Ln(t) • h(t) + L12(t) • i2(t) M*) = Li2(t) • H{t) + L22(t) • i2{t) can be discretized by applying the trapezoidal implici t integration rule, described in Appendix A , to obtain the following system of equations: 2 e2(t)_ ~ A t Ln(t) Ll2(t) Ll2(t) L 2 2 ( t ) »i(*y «2(*) + ehi{t) e/,2(t) (4.23). where the history terms are calculated as ' L n ( t - A t ) L 1 2 ( t - A t ) ' L 1 2 ( t - A t ) L 2 2 ( t - A t ) ehiit) 2 eft2(t)_ ~ ~ A t \(t- A t ) ' "ei(<- At)" > 2 (* - At)_ e 2 ( t - A t ) (4.24) To shorten the notation, in the equations to follow (t) is omitted and (t — A t ) is replaced wi th ' to denote the values of variables from the preceding time step. From the system of equations describing coupled stator (4.18) and rotor (4.19) windings of an induction machine, we can obtain the discretized equations for the 49 Chapter 4. Modelling with Multilevel MATE branch voltages in a similar manner: ea Ra 0 0 0 0 0 " 0 Rb 0 0 0 0 H e c 0 0 Rc 0 0 0 ic eA 0 0 0 RA 0 0 iA e B 0 0 0 0 RB 0 is ec 0 0 0 0 0 Rc ic Laa Lab Lac LaA^r) LaB{6r) Lac{@r) 'ia' eha Tab Lbb Lbc LbA^r) LbB^r) Lbc{@r) ib ehb ^ At Lac LaA(6r) Lbc LbA(Qr) LCc LCA{QT) LCA(9T) LAA LCB{QT) LBA Lcc(9r) LCA ic iA + ehc ehA LaBifir) LbB(9r) LcB{9r) LAB LBB LBC is ens Lacier) Lbc(6r) Lcc{Qr) LAC LBC Lcc }c_ euc Writ ten in compact form we have: &abc &ABC Rs 0 0 R r where the history terms are: + _2_ Xt labc lABC + &h abc &hABC (4.25) eha Ra 0 0 0 0 0 ehb 0 Rb 0 0 0 0 ehc 0 0 Rc 0 0 0 ehA 0 0 0 RA 0 0 ehB 0 0 0 0 RB 0 ehc 0 0 0 0 0 Rc r ~i I a j h h .i %A .1 %B J bc\ 2_ At Laa Lab Lac LaA{9r) LAB(Or) Lac{Qr) 'i' 'e'a Lab Lbb Lbc LbA(9'r) LCA(QT) LbB(e'r) Lbc(0'r) Lcc(0'r) i'b t % Lac Lbc LCc LcB(e'r) V 1 eF aA{Qr) LbA(e'r) LcA(e'r) LAA LBA LCA V eA aB(Qr) ac(8r) LbB(e'r) Lbc{Qr) LCB{e'T) LAB LBB LBC i'p 6P Lcc{8r) LAC LBC Lcc . i'c. .e'c. which written in compact form are: ®/i abc R . 0" 2 ' L S T ^abc f eafac £hABC_ 0 R r ~ At ml }ABC. .eABC. (4.26) 50 Chapter 4. Modelling with Multilevel MATE The above system of equations represents the branch voltage equations for cou-pled stator and rotor windings of an induction machine. This system of equa-tions does not assume any particular connection of the stator and rotor windings (a grounded/ungrounded wye or a delta). The discretized phase-domain induction machine model can be visualized as depicted in F i g . 4.10. ib 24, At stator 24, At eha A A A c - 0 - ^ — \ \ ill, 2L„, -At At At 2A,c\ At X 24, At 'ehc 2iW At <-— 2LA At 24« ^Ah ehA r A + rotor 2LAA 24, Ar ^r-G- /WV 24 A? At ehB 2L„ y ™ \ " At '2L„ r ——11 . At =^£, At \ At * 2LA, At lA IB < . B IC C Figure 4.10: Discrete phase domain induction machine electrical model 4.4.1.2 Mechanical Equations of an Induction Machine The system of equations (4.25) describing the electrical part of an induction machine contains terms L s r , which depend on the value of the rotor angular position 6r at a time instant t. Since the rotor is rotating, the rotor angular position is constantly increasing and is related to the rotor angular velocity ojr as where tor is obtained from the well-known relationship between the machine's accel-erating torque (Ta) and the rotor angular velocity (uir) Ta = Te-Tm = j ^ + DUm = j ( ^ ^ + D ( ^ U r (4.28) where • T e is the electromagnetic torque 51 Chapter 4. Modelling with Multilevel MATE • Tm is the mechanical torque • J is the polar moment of inertia of the rotor and the connected load • pf is the number of poles of the machine • ojm is the angular velocity of the rotor in mechanical radians per second • D is the damping coefficient The differential equations (4.27),(4.28), discretized wi th the implici t trapezoidal rule, become difference equations of the following form: 9r(t) - ^ur(t) = 6r(t- At) + ^-ujr(t - At) (4.29) / 4 7 2D —— + — ) ur(t)-Te(t)+Tm(t) \PfAt pf = (-4 J 2D \PfAt pf ur(t-At)+Te(t-At)-Tm(t-At) (4.30) The electromagnetic torque Te is calculated from the energy stored in the machine coupling field Wf, which is in fact the energy stored in al l mutual and self inductances of rotor and stator windings, excluding leakage inductances. The energy stored in the coupling field can be written in the following form: Wf = statorjenergy + mutual stator-rotor + rotor .energy 1 2 [ia ib ic] 'Lu 0 0 " 'ia' 0 Lls 0 ib 0 0 Lls + [ia ib i c ] - L ^ Aref) A bB Aref) %C + 1 2 Aref) Aref) Aref) •fjref) _ >(re/) Llr 0 0 0 (re, Ir ' 0 Aref) 0 0 - (re/) Jlr Aref) %A Aref) lB Aref) lC (4.31) with rotor quantities referred to the stator side using the effective turns ratio of the stator and rotor windings, and indicated with the superscript (ref\ The electromagnetic torque Te is obtained by solving the following relationship 'Pf\ 8Wf-0dr (4.32) 52 Chapter 4. Modelling with Multilevel MATE which results in the following expression valid in the continuous and discrete time domains: T -L e ia ib i i 1 . JLL(re/) . Cj ddr ST 'ArefY lA Aref) lB Aref) }c ' ("7 )^ ' Lms ' [ia ib ic] " sin 9, sui\9r — sin(0 r + in(0 r - 120°) 120°) sin(0 r sin(0 r + 120°) sin(0 r - 120°)' sin6>r sin(0 r + 120°) Aref) lA Aref) lB A™f) (4.33) 4.4.1.3 Implementation of the Phase-Domain Induction Machine Mode l A system of equations that completely describes the induction machine phase-domain model includes the electrical part, as described by branch equations (4.25), (4.26) and the mechanical part (4.29), (4.30), together with the nonlinear torque equation (4.33). These equations can be combined in a matrix representation of hybrid modified nodal equations of the induction machine in the following general form: A PA P Ajnonlin VA V PA.vref hA / A -ZA 0 1 Asublink + 0 •ia + 0 • Vref = ^Asublink P Ajnonlin 0 ZAjnonlin _ _ 1 Ajnonlin _ _o_ 0 _ ^Ajnonlin _ (4.34) where hA_ref = PA-vref ' Vref Using these equations, the circuit depicted in F i g . 4.11 is modelled using Multi level M A T E to test the phase-domain induction machine model. The circuit consists of a three-phase induction motor (subsystem A) connected to a three-phase sinusoidal voltage source (subsystem B) through links (ias, if,s, ics)- Stator and rotor windings are delta- and wye-connected, respectively. For generality, the rotor side is modelled to allow access to the rotor terminals through links iAr, iBr, icr- In this way, the rotor terminals can either be shorted or connected to an external network, as in the case of a doubly fed induction generator, without any modifications to the induction machine model. The system of branch equations (4.25) describing the electrical part of the induc-t ion machine contains nonlinear terms which depend on the value of the rotor angular position 9r at time instant t. W i t h a sufficiently small integration step, the prediction 53 Chapter 4. Modelling with Multilevel MATE Subsystem A lbs bs Subsystem B vbs as v ' as CS Figure 4.11: Electr ical circuit configuration for testing the phase-domain induction machine model implemented with Mult i level M A T E of the electrical angular velocity u>r can be made using linear extrapolation from the past values ur.pred(t) = 2ur(t - At) -cur(t- 2At) (4.35) and used to calculate the value of the rotor angular position according to (4.29). This procedure is justified, as the mechanical system that includes the rotating masses has slower time responses (the order of seconds) than the electrical circuit, time responses of which are much faster (the order of milliseconds). The induction machine implementation wi th the Mult i level M A T E concept uses the prediction of 6r, but the accuracy of this prediction can be checked and corrected by iterations i f required. W i t h this in mind, the electrical and mechanical systems of equations are coupled only through the nonlinear torque equation relating the machine's electromagnetic torque T e to the rotor and stator phase currents. Therefore, subsystem A, representing the induction machine, is divided according to its nature into subsubsystems Ai and A 2 , where the former represents the electrical part, and the latter represents the mechanical part of the machine. The induction machine coupled branch equations are modelled as sublinks of subsystem A, and the nonlinear torque equation- is modelled as a nonlinear sublink equation. Following the process of building the system of nodal equations A • vA = hA we notice that the machine rotor nodes "A", "B", " C " and " i V " , and stator nodes "a", "6" and "c", form ungrounded subsubsystems due to their connection to the links and/or sublinks. A s a result, the conductance matrix of subsystem A is zero. In such a case, one node of each ungrounded subsubsystem has to be declared as its reference 54 Chapter 4. Modelling with Multilevel MATE point and solved at the links level. Since the ungrounded subsubsystems in this case consist of one node each, al l nodes of the induction machine model are duplicated as reference points vref for their respective subsubsystems. For the circuit in F i g . 4.11, nodal equations for the stator node "a" and the rotor node "A" (taking into account the delta connection of the stator windings and the wye connection of the rotor windings) can be written as Va la ^b ^as Varef ~ VA + lA- iAr ~ VAref = 0 0 (4.36) Since va = varef and VA = VAref, these equations simply state that the sum of the link and sublink currents coming into nodes "a" and "A" is zero. B y expanding this relationship to the remaining induction machine nodes, we obtain the following arrays and vectors associated wi th the nodal equations of the electrical part of the induction machine Ai • VAI + PAX • iAXsubiink + V • ia = hAx - hA_ref (4.37) expanded as 1 0 0 0 0 0 0" ~VA~ 0 1 0 0 0 0 0 Vb 0 0 1 0 0 0 0 VC 0 0 0 1 0 0 0 VA + 0 0 0 0 1 0 0 VB 0 0 0 0 0 1 0 vc 0 0 0 0 0 0 T vN + 0 0 0 0 0 0 -1 0 0. 0 0 0 0 +1 - 1 0 0 0 0 0 0 0 - 1 0 0 0 0 0 +1 - 1 0 0 0 0 0 0 0 -1 0 0 0 - 1 0 +1 0 0 0 0 0 0 0 0 - 1 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 +1 0 0 0 +1 - 1 - 1 - 1 0 0 0 0 0 - 1 0 ^as ibs %Ar iBr lCr ia ic iA iB ic + 0" Varef 0 —Vbref 0 Vcref 0 — —VAref 0 —VBref 0 — VCref 0 —VNref_ Pseudo nodal equations for the mechanical system can be written from (4.29), (4.30), to obtain the following arrays and vectors associated wi th the mechanical part of the induction machine: A-2 • VA2 + PA2jnonlin " iA2.nonlin — h A2 (4.38) 55 Chapter 4. Modelling with Multilevel MATE and _ At 47 ? 2D pfAt pf Ulr + 0„ + ^ w , 2 *~r The mechanical torque Tm(t) acts as a forcing function, and is therefore included on the right-hand side of the system of equations, which corresponds to sources. The sublink equations of subsystem A represent the coupled branch equations of the machine written in the following form, which takes into account the stator and rotor winding connections: PAI VAl — ZA1 - ^Al sublink In the expanded form: = -V, Alsublink (4.39) +1 0 - 1 - 1 +1 0 0 - 1 +1 0 0 0 0 0 0 0 0 0 R 0 0 0 +1 0 0 0 0 R r 0 0 0 0 +1 0 + At 0 0 0 0 0 T T 0 1 Va 0 Vb 0 Vc - 1 VA - 1 VB - 1 _ vc vN 'ia ib L s r ic iA ic _4_ A t 14 "C i'b .1 V i'p i'c. e ? + e'f B'A eP e'c. Finally, the equation representing the nonlinear torque equation can be expressed in the following matr ix form: " 2<A2jnonlin ' 1>A2jnonlin -V, A2.nonlin (4.40) and [1] [Te] = PfEms r . . . -i —Hi— | > lb *c| s i n 0 r sin(0 r + 120°) sin(0 r - 120°)' sin(0 r - 120°) s i n 0 r sin(0 r + 120°) sin(0 r + 120°) s i n ( 0 r - 1 2 O ° ) s i n 0 r Nr 'iA iB where Nr/Ns is the rotor winding to stator winding effective turns ratio. To summarize, the induction machine phase-domain model represented as a six-56 Chapter 4. Modelling with Multilevel MATE terminal device wi l l have the following hybrid form: 0 PAX 0 Px 0 A2 0 P Aljnonlin 0 PAX 0 -zAx 0 0 0 0 0 ZA2jnonlin 0 VAX VA2 ^AXsublink 1 A2jnonlin 0 PAX.vref hA2 0 ^AXsublink 0 _ VA2-nonlin _ 0 'Jref (4.41) B y applying the Mult i level M A T E approach, we obtain the Mult i level M A T E system of equations as follows: I 0 0-AX 0 ax 0 I 0 OjA2-nonlin 0 0 0 PAXAAX + ZAX 0 0 p 0 0 ZA2jnonlin 0 VAX VA2 1> AXsublink i A2jnonlin 0 . " aAXjvref eA2 0 AX sublink 0 _ VA2jnonlin _ 0 •VAXref (4.42) A t every time step of the simulation, the induction machine model w i l l perform the following procedure: • predict the value of 9r using linear extrapolation (4.35) • recalculate [ZAX] based on the predicted value of 9r: ZAX R s 0 0 R r Lsr(9r) 2 + At [Ljr(9r) L • recalculate the contributions of the sublinks to the Thevenin equivalent of sub-subsystem Ai\ Aai = ai - aA1 • {PAX^AX + zAx) 1 PAiai AeAX = eAX ~ aAX- - (PAXaAX + ZAX)'1 {p\xeAX + VAXsublink) (4.43) • calculate the modified Thevenin equivalent of subsystem A as seen from the l inking nodes 6 and submit it to the links \ptaMTE]^ ^^MMTE], where &MTE — dx — Acii &A.MTE = ZAX — AeAx (4.44) receive the links currents [ia] and reference voltages [vref] from the links and 6 The Multilevel M A T E algorithm will recognize that the mechanical system is decoupled from the electrical, therefore CLMTE = &\.MTE and eAMTE = SAI.MTE 57 Chapter 4. Modelling with Multilevel MATE recalculate its nodal voltages: VA — e-AMTE — 0,MTEvre{Vref — CLMTE^a 0>MTEvref — O-Alvref — &&lvref = Q>Alvref ~ aA {PAaA + ZA) PAaA\vref • update the history terms [/i^] a n d [VAi_s«6hnfc] for the next simulation step A t every time step of the simulation the links system w i l l perform the following procedure wi th respect to subsystem A : • receive subsystem's A modified Thevenin equivalent • calculate [ Z Q . M T E ] and [ea_MTE] from Za-MTE — PtO-MTE + (f^MTE + Z Ea.MTE — PTCA-MTE + Q^CB-MTE + Va • calculate the reference voltages of subsystem A from Vref = {Pvref ' Za).MTE ' Pvref ~ zvTef^ ' (pl^f ' Za.MTE ' ea.MTE ~ Kef) • calculate the l ink currents from ia = za^MTE ' {ea.MTE ~ Pvref ' Vref) • return [vref] and [ia] to subsystem A . The Mult i level M A T E algorithm wi l l automatically recognize that the mechanical system is decoupled from the electrical, therefore no iterations are required to solve the nonlinear torque equation. Furthermore, the induction machine model w i l l appear as a 6x6 Thevenin equivalent as seen from the machine's terminals (3x3 i f the rotor terminals were shorted and not exposed to the network). The phase-domain induction machine model presented here has been tested and successfully compared wi th results from the literature [30]. The machine parameters are given in Appendix D . l . In the following graphs, results are presented first for the free acceleration of a 2250 hp induction motor from stall, and then for step changes in load torque from zero to nominal (8900 Nm) at 0.5 s and from nominal back to zero at 2.5 s, wi th the machine init ial ly operating at synchronous speed. 58 Chapter 4. Modelling with Multilevel MATE Figure 4.13: Rotor phase current during free acceleration of a 2250 hp induction motor 59 Chapter 4. Modelling with Multilevel MATE 60 Chapter 4. Modelling with Multilevel MATE 61 Chapter 4. Modelling with Multilevel MATE -10001 1 1 1 1 1 1 0 0.5 1 1.5 2 2.5 3 t(s) Figure 4.18: Stator phase current during step changes in load torque of a 2250 hp induction motor 62 Chapter 4. Modelling with Multilevel MATE 2000 0 -2000 1.5 t(s) Figure 4.20: Electromagnetic torque during step changes in load torque of a 2250 hp induction motor 1830 1820 1810 1780 1760 ! ! V I 11 I I I I 1 1 0 0.5 1 1.5 2 2.5 3 t(s) Figure 4.21: Rotor speed during step changes in load torque of a 2250 hp induction motor 63 Chapter 4. Modelling with Multilevel MATE 4.4.2 Phase-Domain Synchronous Machine Model A phase-domain model of a synchronous machine [25] has been implemented in O V N I in an analogous manner to the induction machine model, using sublink equations and the Mult i level M A T E concept. The synchronous machine model is composed of the stator phase windings (a, b, c), equivalent to those of an induction machine, and the field rotor winding (fd) supplied by a D C excitation voltage. In addition, the synchronous machine model has equivalent damper windings in the d and q axes to model the paths for induced rotor currents. The number of damper windings varies depending on the type of synchronous machine to be represented. A generic diagram of the electrical part of a wye-connected synchronous machine is depicted in F i g . 4.22. The modelled machine is a steam turbine synchronous generator with non-salient poles. The equations representing the electrical part are derived in a similar manner as described for the induction machine model in the preceding section. The mechanical equations of the synchronous machine are equivalent to those of the induction machine. Figure 4.22: Rotor and stator circuits of a synchronous machine The phase-domain synchronous machine model implemented wi th Mult i level M A T E has been tested and successfully compared against results from the literature [30]. The machine parameters are given in Appendix D.2. Results of dynamic performance of a 835 M V A steam turbine generator during a step increase in input torque from zero to fifty percent of the rated value are shown in F i g . 4.23 to F i g . 4.28. 64 Chapter 4. Modelling with Multilevel MATE — 0 Figure 4.23: Stator phase current during a step increase in input torque of a 835 M V A steam turbine generator Figure 4.24: Field current during a step increase in input torque of a 835 M V A steam turbine generator 65 Chapter 4. Modelling with Multilevel MATE A ! ! 1 0 1 2 3 4 5 6 7 8 9 10 Time Figure 4.25: Electromagnetic torque during a step increase in input torque of a 835 M V A steam turbine generator 378.5 377.5 377 0 1 2 3 4 5 t(s) 6 7 8 Figure 4.26: Rotor speed during a step increase in input torque of a 835 M V A steam turbine generator 66 Chapter 4. Modelling with Multilevel MATE 601 1 1 ! 1 1 1 1 1 r _ 1 0 L _ 1 1 1 1 1 1 1 1 1 1 0 1 2 3 4 5 6 7 8 9 10 1(s) Figure 4.27: Rotor angle during a step increase in input torque of a 835 MVA steam turbine generator _2I i i i I I I I -10 0 10 20 30 40 50 60 S (el. deg) Figure 4.28: Torque-rotor angle characteristics during a step increase in input torque of a 835 M V A steam turbine generator 67 Chapter 4. Modelling with Multilevel MATE 4.4.3 Commentary on Modelling of Electric Machines Our mathematical description of the implementation of a phase-domain machine model wi th Mult i level M A T E may give rise to questions about how the efficiency of computation for these models can be evaluated. The main disadvantage of a phase-domain machine model is that its impedance matr ix is time-dependent due to rotation of the rotor. Furthermore, i f nodal analysis is used, time-dependent coupled branch equations have to be transformed into nodal equations at every time step. Even if modified nodal analysis is used, the machine's branch equations wi l l increase the size of the system of equations to be solved, and cause the matr ix of coefficients to change at every time step of the simulation.. In the approach proposed in this the-sis, the changing impedance matrix of the machine modelled wi th branch equations is recalculated at every time step separately from the subsystem's nodal equations, therefore the subsystem's nodal admittance matrix does not need to be inverted at ev-ery time step. The contribution of the machine's branch equations to the subsystem's nodal equations is integrated at the level of the subsystem's Thevenin equivalent (see equations (4.43), (4.44)) which makes this approach more efficient than the M N A or nodal analysis implementation approaches. Implementation of the phase-domain machine models in the E M T P - t y p e of sim-ulators provides significant advancement over the commonly used dqO model with respect to the simulation time step necessary to interface the machine model and the network equations. In this work the simulation step of 1 ms was used, whereas a typ-ical time step used in the E M T P simulations with a dqO model is in the order of 50 fj,s or less. Small time step wi th the dqO model is necessary to avoid instabilities due to predictions used in interfacing the model with a three-phase network representation. 4.5 Modelling of Controllers A traditional approach to simulating controllers in power system networks uses an artificial time step delay between the network solution and the solution of the control systems to decouple the two systems. In recent years, different approaches have been proposed to achieve a simultaneous solution of the network and controller equations [31],[32],[33],[34]. W i t h i n O V N I , controller equations of the form A^i • x = hctri (4.45) . where • [Actri] is a matrix of controller coefficients • [x] is a vector of controller variables (inputs + outputs + internal variables) 68 Chapter 4. Modelling with Multilevel MATE • [hctri] is a vector of controller history terms can be solved simultaneously using a hybrid system of nodal and branch equations describing the electric network. Controllers are generally represented by block diagrams showing interconnections among different control blocks (elements), such as transfer functions, limiters, etc. Discretized equations of controller elements can be obtained in OVNI by mapping the Laplace operator s to the z-domain using the bilinear transformation (trapezoidal discretization), with a discrete time step At: 2 1 — z " 1 s = 4-- T—^ (4-46) At 1 + v ; For example, discretization of the first-order transfer function of the general form H(s) = Kb-°±^ = • (4.47) v ' a0 + a i 5 Xi(s) v ' where X\{s) and X2(s) are the input and output of a transfer function in the Laplace domain, is obtained by substituting (4.46) into (4.47). The following first-order dif-ference equation is obtained with and -B0x1(t) + A0x2(t) = h{t) (4.48) h(t) = BlXl{t - At) - Axx2{t - At) (4.49) 2 2 A0 = ao + ai— Ai - a0 - fli — The term h(t) in (4.49) depends on the values of the input and output from the previous time step, therefore it is referred to as a history term of a first-order transfer function. Controller equations at the subsystem level (local control) can be incorporated into the MATE system of equations as functional sublinks. In the case of nonlinear controllers containing, for example, limiters, their equations can be treated as non-linear sublinks as described, in Section 3.6. In the most general representation, the 69 Chapter 4. Modelling with Multilevel MATE controller equations at the subsystem level wi l l have the following form: A PA Petri VA hA P \ - Z A 0 0 P lA - V A Petri — ZA_ctrl Zctrl • ictri = — Vctrl (4.50) 0 B Q VB hB pt —z ia ~V°-where • [Petri] is a connectivity array of controllers in subsystem A , and Petri ' VA - ZA.ctrl ' U ~ Zctrl ' ictri = ~Vctrl (4-51) represents the system of controller equations where [p*t7.(], generally speaking, is not the transpose of the connectivity array \pctri}- As can be noticed from (4.50), the controller equations introduce asymmetry into the sublink equations. If a controller operates at the system level (global control), its equations can be incorporated into the M A T E system of equations as links. Such a controller uses global measurements of power system variables from multiple subsystems as inputs and/or controls variables from multiple, subsystems as outputs. A global controller would provide, for example, damping of inter-area oscillations. The third type of controller is implemented at the level of a power system com-ponent (component control). Its action is not transparent to the subsystem, but it is embedded in the information passed by the power system component to the subsystem. The implementation of controllers in O V N I using Mult i level M A T E is demon-strated in the following example and described in [33]. 4.5.1 Example: Modelling of a Single-Machine Infinite-Bus System in dq Coordinates with Multilevel M A T E A single-machine infinite-bus power system from [27] is modelled using the M u l t i -level M A T E approach. The system consists of an equivalent synchronous generator representing four generating units of a power plant, a step-up transformer and a double-circuit transmission line connecting the generator to a large power system represented by an infinite bus. The synchronous generator is modelled in the dq ref-erence frame, including the rotor dynamics. The model includes a bus-fed thyristor 70 Chapter 4. Modelling with Multilevel MATE excitation system wi th an automatic voltage regulator ( A V R ) and power system sta-bilizer (PSS). In the example, network transients and stator transients are neglected as is traditionally done in power system transient stability analyses. The electrical discrete-time circuit of the example system is shown in F i g . 4.29. Excitation system with A VR and PSS kd SM Electrical System kq e/d_hist CO ekd hist (Z) Zhjjtist SM Mechanical Svstem Synchronous Machine (SM)fNetwork Interface Figure 4.29: Single-machine infinite-bus equivalent network Aa)r \ + sTw 1 l + sTR I: > l + sTx 1 + sT F min Figure 4.30: Thyristor excitation system wi th A V R and PSS The bus-fed thyristor excitation system with the A V R and the PSS controller are shown in F i g . 4.30. Nonlinear equations associated wi th the limiters are described as: -x3(t) + xs(t) = 0 for x S m i n <xs < x S m a x max for xs max < XS — Xs min for xs < xs min 71 Chapter 4. Modelling with Multilevel MATE Efd(t) + KA (xi(t) - xs(t)) = KAVref for E F m i n < Efd < E F m a x Efd — EFm&x for E F m a x < Efd Efd = E F m m for Efd < Ep min A complete system of equations of.the excitation system discretized using the trapezoidal rule (bilinear transformation) is shown in (4.52). 0 0 1 KA(0,0) 0 0 -KA(O,O)" - 1 0 0 0 0 -KSTAB iTw£i) 0 0 0 0 0 0 0 0 0 0 0 0 - 1 ( 0 , 0 ) 0 0 1 Et(t) Aur(t) Efdjt) xi(t) x2(t) xs{t) xS{t) KAVref {EFmax, EFmin) Et(t-At)-(1 TRAt -KSTAB (TW£) A w r ( i - At) t - L 0 {xsmax; 3-Smin) ( l - T W At (l - T,£) x2(t At) - (l - T2±) x3(t At) 2)x2(t At) At) (4.52) where • Et is the synchronous machine stator terminal voltage • Auir is the per unit relative rotor angular velocity • Efd = ^jrefd is the exciter output voltage Mult i level M A T E system partitioning is next applied to separate the main com-ponents of the system: the network and the synchronous machine wi th the exciter. The system partit ioning structure is shown in F i g . 4.31, where • [VA] represents the network's nodal voltages in the dq reference frame • [iA] represents the network's branch currents for ideal switches (i swli isw2 ' ' ' i n Fig . 4.29) in the dq reference frame • [VB] represents the synchronous machine stator and rotor voltages, linear me-chanical variables (5, Ator) and linear exciter variables (xi,x2 and x3 in F ig . 4.30) 72 Chapter 4. Modelling with Multilevel MATE •A network PA P P ' A -ZA BSM IB q q'B -zB P' 4 -z IA IB -VA Figure 4.31: Mult i level M A T E partitioning of the single machine infinite bus system • [is] represents the synchronous machine stator and rotor currents, nonlinear mechanical variables (Te) and nonlinear l imited exciter variables (xs, Efj and Et in F i g . 4.30) • [ia] represents the terminal dq l ink currents and the rotor angle 5 The synchronous machine model in the dq reference frame is better described by branch equations, as wi th the phase-domain model. The machine's branch equations are included as functional sublinks of subsystem B representing the synchronous machine. The stator-side nodal equations are derived in a similar way as for the phase-domain induction machine model (4.36). The rotor-side nodal equations state that the nodal voltage of the field winding is equal to the exciter voltage (with appropriate scaling) and that the nodal voltages of the damper windings are zero. The nonlinear equation for calculating the terminal voltage Et: as well as the nonlinear torque equation and the equations of the exciter's limiters are included as nonlinear sublinks of subsystem B and are solved by iterations at every time step. Subsystem A represents the connecting network wi th the infinite bus, al l vari-ables decomposed into dq coordinates. Since the subsystem does not have a direct connection to ground (only through the sublinks and the ideal voltage source repre-senting the infinite bus), the procedure described earlier for ungrounded subsystems is applied to reference one nodal voltage for each ungrounded subsubsystem. The 73 Chapter 4. Modelling with Multilevel MATE nonlinear equations associated wi th the network relate the infinite bus voltage to the rotor angle 5 as eBd - EB sin 5 eBq = EB cos 5 The transient response of the system variables to a three-phase short circuit ap-plied on Circui t 2 of the transmission line is simulated using the Mult i level M A T E approach. The results for the system's currents and voltages, as well as for the syn-chronous machine's mechanical variables, are depicted in F i g . 4.32, 4.33. The results agree well wi th the time responses presented in [27]. 4.5.2 Commentary on the Simultaneous Solution of Controller Equations A time delay between the solution of control and network equations can cause instabil-ity, inaccuracy and numerical oscillations [31]. These problems are well documented in the literature, wi th different solutions suggested to either reduce or overcome them. In the case of the E M T P - b a s e d programs, implementation of the T A C S [35] module is an example of a non-simultaneous approach to solving control equations. The approach to the simultaneous solution of control and network equations de-scribed in this thesis takes advantage of the proposed Mult i level M A T E concept and modelling wi th functional sublinks. As nonlinearities iterate individually only against the system of links, which has to be optimally small due to requirements for the multimachine implementation of O V N I [10], the framework for iterating nonlinear equations is computationally efficient. For linear controller equations, no iterations are necessary. For nonlinear controller equations, iterations w i l l be performed within the framework described in Section 3.6. W i t h the fixed point type of iteration, an average value of one to two iterations were needed in test cases described in this thesis, with the number increasing to about five to ten iterations for a few time steps following a sudden dramatic change of system conditions (for example, a short circuit, or during init ialization from zero ini t ial conditions, Section 3.6.1). The results published in [34] using Newton iteration show iteration counts of the same order as those obtained in this work. Newton iterations, however, have a higher solution overhead due to re-calculation of a Jacobian matrix at each simulation time step. 74 Chapter 4. Modelling with Multilevel MATE (a) . Generator terminal voltage for the exciter with manual control 5(PU) \.c • 1 0.8 0.6 0.4 0.2 A \ y ~ — — ^ u 3 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t(s) (b) Active power response for the exciter with manual control Pe (MW) 2500 2000 1500 1000 500 0 D 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t(s) (c) Rotor angle response for the exciter with manual control 8 (deg) 120 100 80 60 40 20 u 3 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t(s) Figure 4.32: Transient response of the single-machine infinite-bus system with manual control: (a) terminal voltage, (b) active power, (c) rotor angle 75 Chapter 4. Modelling with Multilevel MATE (a) Generator terminal voltage for the exciter with AVR and PSS 5(pu) A rx 1 0.8 0.6 0.4 0.2 n u -D 0.5 1 1.5 2 2.5 3 3.5 4 4.5 t(s) (b) Active power response for the exciter with AVR and PSS Pe(MW) 3000 1000 - - - -500 : o I : LI : 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 UsL (c) Rotor angle response for the exciter with AVR and PSS 8 (deg) 140 120 100 80 60 40 20 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t(s) Figure 4 . 3 3 : Transient response of the single-machine infinite-bus system with A V R / P S S : (a) terminal voltage, (b) active power, (c) rotor angle 76 Chapter 4. Modelling with Multilevel MATE 4.6 Summary The concept of Mult i level M A T E is further explored in this chapter. New solutions are presented for modelling of ideal voltage sources and ungrounded circuits and subsystems. The modelling of electric machines with branch equations by means of sublinks is then described. Finally, modelling of controllers wi th Mult i level M A T E is explained and demonstrated through an example. 77 Chapter 5 Modelling and Simulation of a DFIG Wind Turbine System in OVNI 5.1 Introduction Doubly fed induction generator (DFIG) wind turbines are increasingly used in new wind turbine installations al l over the world. Growing concerns about the impact of a large number of wind turbine generators ( W T G ) on transient and voltage stability of power system networks has led engineers to revisit modelling and simulation practices traditionally used for system stability analyses. Because O V N I is based on the E M T P methodology for accurate detailed modelling, and the Mult i level M A T E concept, which, combined wi th hardware solutions, allows for fast simulation of large power system networks, it represents an ideal tool for testing and developing benchmark models of different wind turbine installations. The advantages of the detailed modelling and simulation of D F I G wind turbine systems wi th O V N I include: (1) three-phase representation allowing for accurate sim-ulation of balanced and unbalanced network conditions; (2) phase-domain full-detail machine modelling, including the modelling of torsional shaft oscillations; (3) mod-elling of the W T G control systems; (4) and the modelling of power electronic voltage source inverters. The purpose of this chapter is to demonstrate how an experimental D F I G wind turbine system from [36] is modelled in O V N I and to describe solutions for modelling challenges associated with the full time-domain ( E M T P - t y p e ) simula-tion of power systems. Finally, the results of the simulations are compared against traditional stability analysis simulation results obtained wi th the Transient Stability Assessment Tool ( T S A T ) . The research contributions reported in this chapter include the following: • Description and implementation of detailed modelling of a doubly fed induction generator wind turbine system with Multilevel MATE. • Comparison of the OVNI simulation results against the experimental results and transient stability simulation. 7 8 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 5.2 Wind Turbine Generators Variable-speed turbines with doubly fed induction generators have become a new standard for wind turbines of installed capacity above 2 M W . As the ratings of such wind farms connected to the power system grid become closer to the ratings of tra-ditional generating units, and their share in the total installed generating capacity of the power system becomes considerable, it becomes necessary to perform studies of the impact of large wind farm connections to the power system network. The main concern is to study the responses of wind farms to power system faults and their impact on overall system stability. Variable-speed operation of the turbines is achieved through the use of power electronics converters that can also be used to improve the grid integration aspects. It is anticipated that in the future it may be possible to request specific responses of wind farms to network disturbances from the manufacturers to help system recovery. Detailed model parameters of converter-controlled wind turbines can only be provided by manufacturers, and these control details are usually confidential and not readily available. However, efforts are being made to create reasonably accurate general models of doubly fed induction generators that can produce realistic results of wind farm responses to system disturbances and the influence of the associated controls. Fixed-speed wind turbines are easier to model, as their responses mainly depend on the nature of their components, such as the response of induction generators, shaft stiffness and inertia of the rotor parameters, the details of which are normally available. There are four types of wind turbines in use today: • Fixed-speed wind turbine wi th a conventional induction generator connected directly to the grid. This turbine was originally only stall-regulated, but in recent years the larger MW-s ize turbines have had an active p i t ch 7 mechanism as well. • Pitch-regulated variable rotor resistance induction generator has a variable rotor resistance that can vary the slip in the range of 2-10% and allow a sufficiently variable speed operation to reduce output power fluctuations. • Doubly fed induction generator has become a new standard wind turbine for larger units. A slip-ringed induction, machine is used as a generator, with a converter connected between the rotor slip rings and the network. The stator of the machine is directly connected to the network. 7 Wind turbine blades have a controllable pitch angle to turn the blades out of the wind or into the wind. 79 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI • Multi-pole synchronous generator connected to the power network v ia a con-verter. The synchronous generator is a direct-drive type wi th low speed and a higher number of poles. The wind turbine and the generator rotate at the same mechanical speed v ia the same shaft. 5.3 Doubly Fed Induction Generator Wind Turbines Physically, the machine of a doubly fed induction generator is a conventional wound rotor induction machine. The key distinction is that this machine is equipped with a solid-state A C excitation system. The A C excitation is supplied through an A C - D C -A C voltage converter. Doubly fed induction machines have a significantly different dynamic behaviour than conventional synchronous or induction machines. The fundamental frequency electrical dynamic performance of a doubly fed in-duction generator is completely dominated by the field converter. Conventional as-pects of the generator's performance related to the rotor angle, excitation voltage and synchronism are largely irrelevant. The electrical behaviour of the generator and converter is that of a current-regulated voltage-source inverter. The converter makes the wind-turbine behave like a voltage behind a reactance that produces the desired active and reactive current delivered to the device terminals. A schematic of a doubly fed wind turbine system wi th two voltage-source inverters and accompanying controls is depicted in F i g . 5.1. 5.4 DFIG Model Structure To construct a D F I G wind turbine system model, the following model structure is normally considered: • Doubly fed induction generator model based on a phase-domain induction ma-chine model • Voltage converter model based on an average model of a back-to-back P W M voltage-source inverter wi th stator and rotor-side converter control • W i n d model that maps the wind speed to the shaft mechanical power for the turbine. Emulat ion of wind disturbances such as gusts and ramps by varying input wind speed to the wind-power module • Crowbar protection 8 0 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI rotor A- A | voltage converter ADC-link D C D C A C J i i i PWM VDC \ r PWM rotor-side converter controller stator-side converter controller • \ i stator flux angle calculation stator voltage angle calculation VDc iqi va,vb,vc to gri Abie, Figure 5.1: Schematic of a doubly fed induction generator wind turbine system • Mechanical controls, including blade-pitch control and over/under speed trips 5.5 Doubly Fed Induction Generator Model A doubly fed induction generator can be modelled as an induction machine in phase coordinates. The difference between a conventional induction machine and a doubly fed induction machine is that the rotor of a doubly fed induction machine is connected to the grid via a voltage converter. The electrical equations of the two machines are identical, as discussed in Section 4.4.1. W i t h respect to the mechanical equations, the detailed modelling of D F I G wind turbines requires a two-mass shaft representation, which we describe in the following subsection. 5.5.1 Two-mass Representation of a DFIG W T G Shaft Torsional shaft oscillations of variable-speed wind generators are not transferred to the grid under normal operating conditions due to the fast active power control of the converters. In this single-mass shaft model representation is usually sufficient. However, when the system response to heavy disturbances is analyzed (such as short circuits in the network), the generator and turbine acceleration can be simulated wi th sufficient accuracy only if shaft oscillations are included in the model. In this case, the 81 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI shaft has to be approximated by at least a two-mass model, one mass representing the turbine inertia and the other mass representing the generator inertia. The inertia of the gear-box is usually not modelled separately but rather included in the generator inertia. The mechanical coupling of the turbine and generator through the gear-box can be represented by two spring-connected rotating masses described by the rotational form of Newton's second law [28]: [ J ] d%[d] + [ D ] i t [ 6 ] + [ K ] [0] = [ T t ] _ [ T a ] d9_ dt (5.1) - OJ where • [J] is a diagonal matr ix of the turbine and generator moments of inertia ( J t , J 9 ) • [9] is a vector of angular positions • [to] is a vector of angular velocities • [D] is a matr ix of damping coefficients • [K] is a matr ix of stiffness coefficients (spring constants) • [Tt] is a vector representing the turbine mechanical torque • [T9] is a vector representing the generator electromagnetic torque turbine A ftgear box+generator Ktg h-AAA-D,0 Jo Figure 5.2: Two-mass representation of a D F I G W T G shaft for accurate simulation of torsional shaft oscillations 82 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI The above equation for the two mass turbine-generator coupling depicted in F i g . 5.2 can be writ ten as: 'Jt 0" d2 ~et~ + 0 J9. dt2 69. Dt + Dtg -Dtg -Dt„ Dg + Dtg d_ ~dt Ktg -Ktg which when substituted wi th [co] = [%'g] (5.2) becomes T t = J t l t + D t U J t + D t g { U t ~ ^ + K t g { - 9 t ~ 0 g ) 0t 0 0 T (5.2) dt - Jg^ + DgUJg - Dtg(Ut - Ug) ~ Ktg(9t - 6g) These two equations can be written in the following form: - duig T m - T g = Jg-jj- + DgLOg Tt-fm = Jt^r + Dtut dt (5.3) where Tm = DtgiUt ~ Ug) + Ktg(9t - Qg) (5.4a) (5.4b) (5.5) Equation (5.4a) represents the generator inertia and is equivalent to (4.28). The shaft model represents the turbine inertia (5.4b) and the coupling between the tur-bine and the generator (5.5). The difference equations describing the generator and turbine's inertia are produced analogously to (4.30) presented in Section 4.4.1.2. Note that the turbine torque Tt is calculated from the mechanical power extracted from the wind using the following equation: Tt = wind (5.6) 5.6 Voltage Converter Model and Control The converter and its controls highly influence the dynamic response of a D F I G . In this section, rotor- and stator-side converter modelling w i l l be described, as well as the modelling of their controls. The models represent typical equipment and control structures derived from [36]. The voltage converter that supplies the rotor of a doubly fed induction generator consists of two voltage-source inverters linked v ia a D C - l i n k capacitor, as shown in F i g . 5.1. The voltage-source inverter is connected to the network and to the D F I G stator terminals. It is, therefore, referred to as a stator-side P W M converter. The 83 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI voltage-source inverter connected to the D F I G rotor circuit is referred to as a rotor-side P W M converter. The rotor and stator-side P W M converters are self commutated converters and are set up by six pulse bridges, as depicted in F i g . 5.3. VB O -rotor-side PWM converter UL1 DC-link <DCr 'DCs IDC C i VDC stator-side PWM converter F T - 1 / ,VaI \Vbl Vcl Figure 5.3: Rotor and stator-side P W M converters of a D F I G Assuming an ideal D C voltage and an ideal P W M modulation (infinite modulation frequency), the fundamental frequency line-to-line stator-side converter voltage R M S value (VsUine) and the D C voltage (vDc) can be related to each other as V 3 VglJine — ~~7^Prn,VDC 2^/2 (5.7) with a similar expression for the rotor-side converter voltage-transfer characteristic. The pulse-width modulation factor Pm is the control variable of the stator-side P W M converter. The converter model is completed by the power-conservation equation (assuming a lossless converter) expressed for the stator side as VDciDC - v / 3 R e ( V v s i j i „ e / * 1 J i n e ) = 0 (5.8) The rotor side has a similar expression, where Vsuine and Isuine are phasor quanti-ties of stator-side converter A C line voltage and current, and * denotes a complex-conjugate value. The switching frequency of P W M converters is usually several hundreds Hz, and the average switching losses are proportional to the square of VDC- The switching losses can be taken into account by including a resistance between the two D C ter-minals shown in F i g . 5.3. The approximate fundamental frequency modelling approach of P W M convert-ers (average model) neglects al l switching operations occurring within the voltage-source inverters and represents the converter as an ideal, lossless DC-to-fundamental-frequency A C converter complying with expression (5.7). According to the results presented in several publications, the high-frequency ripple due to switching harmon-ics caused by the P W M operation of the voltage converter is of no significance for 84 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI studying the performance of D F I G W T G in response to network events [37], and in practice is small and further l imited by the inclusion of supply-side inductors [36]. Hence, in this example the average model of P W M converters is considered sufficient. Because the induction generator and the network are both modelled in phase coordi-nates (abc reference frame), the converter equations are modelled in phase coordinates as well. 5.6.1 Stator-side Converter Model When deriving the equations of the stator-side converter, it is assumed that the converter is connected to the grid v ia a three-phase line wi th per phase inductance L and resistance R (see F i g . 5.1). The voltage balance across the line is Va Vb Vr R 1>al ibl id + L d_ dt hi id + Val Vbl Vd (5.9) where va, v0, vc are the D F I G stator terminal phase voltages, vai, ff,i, wci are the stator-side converter A C terminal phase voltages, and iai, ibi, ici are the stator-side converter A C terminal phase currents. The power conservation equation between the D C link and the stator-side con-verter, wi th losses neglected, can be written in phase coordinates as P, Valial + Vblibl + Vdicl = VDCiDCs (5.10) where al l quantities are instantaneous values and PCOnv is the three-phase active power of the voltage converter. The calculation of active power from instantaneous values of phase currents and voltages is valid under all operating conditions, transient and steady state operation, balanced and unbalanced conditions and sinusoidal or non-sinusoidal waveforms [38]. The stator-side converter is vector controlled in the stator voltage-oriented dq reference frame. The angular position of the stator voltage vector V 5 can be calcu-lated from the stationary a/3 stator voltage components. Transformation from the stationary abc reference frame to the stationary a/3 reference frame (5.11) is done using Clarke's transformation (Appendix B . l ) . The vectors of corresponding stator voltages are shown in F i g . 5.4. The variables in the reference frame are scaled to have the same magnitude as the original phase-domain variables. A scaling factor of 2/3 is used in the equations. Vsa _ 2 1 cos yS0_ ~ 3 0 sin cos Va Vb Vr. (5.11) 85 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI b-axis \ \ q-axis ws \ \ 3/2 v: fi-axis d-axis 3/2 \ s ^ \ / Ois / / / Va/ 3/2.vsa a-axis, a-axis c-axis Figure 5.4: Geometric representation of the relationship between stationary abc and a/3 reference frame stator voltages, and definition of stator voltage oriented dq refer-ence frame The angular position 6S of the stator voltage vector V s can be calculated in the most general case as 9S = y ^ d i ^ t a n - 1 ^ (5.12) where us is the angular frequency of the stator voltage. This equation can be written in terms of the stator phase voltages using the inverse Clarke transformation as 9S = t a n - 1 [ 2 \ 2 i c ] (5.13) In the particular situation considered here, the stator voltages' magnitude and fre-quency are dictated by the infinite bus. Therefore, the stator voltage' angular fre-quency is constant, and 9S = ust. Transformation of the phase variables in the stationary abc reference frame to a 86 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI rotating stator voltage oriented dq reference frame, where dq variables are scaled to have the same amplitude as the R M S value of phase variables, can be written for the stator voltages as Vd 1 2 V2 ' 3 cos 9S — sin #, cos (9S - f ) cos (9, + f ) - s i n ( 0 s - f ) -sm(9s + f) Vb Vr (5.14) The d-axis of the rotating dq reference frame is aligned along the position of the stator voltage vector V s , as shown in F ig . 5.4. The equations relating stator voltages in the abc phase-domain frame and the dq rotating reference frame can be written using the inverse dq transformation as va = %/2V s cos (cost) = V2vd cos 9S — V2vq sin 9S 2TT\ („ 2TT Vb V2VS cos f cod — = y/2vd cos ^9S vc = V2VS cos ( co, ( 2n 9S —— t + y ^ = V2vd cos ^9S + y ) - V2vq sin (^9S + y (5.15) where the dq stator voltages (vd,vq) are scaled to have the same magnitude as the R M S values of the stator phase voltages (Vs). From (5.15), we notice that by aligning the d-axis of the rotating reference frame with the stator voltage vector position, the (/-component of the stator voltage is equal to zero. The amplitude of the d-component is equal to the R M S value of the stator phase voltage and, under balanced network conditions, is constant. Reference frame transformations are described in more detail in Appendix B . If the losses of the three-phase line connecting the converter wi th the grid are neglected (R = 0), the R M S phase values of the stator voltages (Vs) and the converter voltages (Vsi) are equal. B y taking into account the stator-side converter voltage transfer characteristic (5.7), we can write Va = vd-2\/2 • vDC (5.16) The power-balance equation for the stator-side converter can be written in dq coordinates from (5.10), wi th the appropriate scaling of the dq variables, as VDC^DCS = 3vdidi (5.17) B y combining (5.16) and (5.17), we obtain a direct relationship between the D C 87 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI stator-side converter current IDCS and the d-axis stator-side converter current idi, as 3 iDCs — 2V2 Pmidl (5.18) From Kirchoff's 1 s t law for the D C link currents of the voltage converter, we obtain the following expression, wi th the resistive losses neglected, „dvDC . 3 C— - lDCs ~ IVCr - ^= PmUl — iDCr (5.19) where C is the value of a D C - l i n k capacitor as depicted in F i g . 5.3. Hence the DC- l ink voltage VDC can be controlled v ia idi-5.6.2 Stator-side Converter Control The objective of the stator-side converter control is to keep the D C - l i n k voltage constant, regardless of the magnitude and direction of the rotor power. A vector-control approach is used wi th the reference frame oriented along the stator-voltage vector position, enabling independent control of active and reactive power flowing between the stator and the stator-side converter. The stator-side P W M converter is current regulated, wi th the direct-axis current component idi used to regulate the D C - l i n k voltage and the quadrature-axis current component iqi used to regulate the reactive power. vDC Stator-side DC voltage controller (Stage 2) idi K, V I z - 1 idi + • PI iql + o hi Stator-side current controller (Stage 1) f z - 0-i \ \ z - 1 ) PI Ki z - a t \ z - l l PI — • "9 Figure 5.5: Stator-side converter controller The first stage of the stator-side converter controller is a current, controller that 88 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI contains two current control loops, for the d- and g-axis current components. The reference of the direct-axis current component is set by a slower D C voltage-control loop (Stage 2), as depicted in F i g . 5.5. The reference of the quadrature-axis cur-rent component corresponding to the reactive power can be used for optimum,power sharing between the generator and the grid-side converter or kept to a constant value. The outputs of current controllers v'd and v'q represent the regulated voltage drops across the stator-side connecting line. (L,R), which is used to calculate the reference values vdl and vql of the stator-side converter. Equation (5.9) can be written in the dq reference frame as vd = Ridl + L^- - cosLiql + vdl dt r-,. T diQi Vq = Rlqi + LjT~ + 0JsLldi + V, (5.20) 91 Taking into account that vq — 0 in the stator-voltage-oriented dq reference frame, reference voltages of the stator-side converter are calculated from (5.20) as v*di = -v'd + {usLiql + vd) v*qi = ~v'q - (WsLidl) (5.21) Finally, the reference voltages are converted back to phase domain using the in-verse dq transformation, as follows: Ki Ki COS0, cos (8 2?r 2l — sin 9S sin (9* — cos(0 S + 4? - s i n (9S + ^ 2TT 3 JJ Jdl (5.22) where the multiplier y/2 accounts for the scaling of the dq variables to R M S phase values. 5.6.3 Rotor-side Converter Mode l The rotor-side converter is connected directly to the rotor terminals of a D F I G . There-fore, the rotor-side converter A C terminal voltages are the voltages across the rotor windings. The power conservation equation between the D C link and the rotor-side converter, wi th the losses neglected, can be written similarly to (5.10) as Pconv = V A i A + VBil3 + VCic = VDCiDCr (5.23) where al l voltages and currents are instantaneous phase values. The rotor-side converter is vector controlled in the stator flux-oriented dq reference 89 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI frame. The angular position of the stator-flux vector ^ s can be calculated from the stationary afl stator-flux components in an equivalent way as presented for the stator-voltage vector angular position in Section 5.6.1. However, i n this case the instantaneous values of the stator-flux linkages ipa, ipb and ipc (calculated from (4.21)) have to be filtered out using a digital pass-band filter to eliminate D C offsets. Filtered values are denoted as ipa_F, ipb.F and iftc_F in F ig . 5.6, depicting relationships between stator-flux linkage vectors in different reference frames. Note that by aligning the d-axis wi th the stator-flux vector position, the ^-component of the stator-flux vector is equal to zero. Dig i ta l filtering of instantaneous stator-flux linkages is described in Appendix C . ft-axis I d-axis a-axis i i / / / / / / / ' / i i i i / i i i / c-axis / Figure 5.6: Geometric representation of the relationship between stationary abc and a/3 reference-frame flux linkages, and definition of stator-flux-oriented dq reference frame The angular position of the stator-flux vector 9sf is obtained from the instanta-neous values of the filtered stator-flux linkages as \1pa-F ~ \i>b.F ~ J 90 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI The transformation of the rotor-phase variables into a rotating dq reference frame, where the dq variables are scaled to have the same amplitude as R M S values of phase variables, can be written for the rotor voltages as Vdr 1 2 COS (6siip) Vqr - sin (dsup) VA VB vc (5.25) where 9siip = 9sf — 9r is the slip's angular velocity and 9r is the rotor's angular position, defined in Section 4.4.1. The transformation is valid for the rotor voltages and currents as well. Equations relating rotor voltages in the abc phase domain and the dq rotating reference frame can be written using the inverse dq transformation as vA = V2Vr cos (usiipt + ip) = V2vdr cos (0sKp) - V sin (0sMp) V B = y/2Vr cos (^jsHpt + <p- = ^ v d r cos ^ 9sHp - y - V2vqr sin ^ V c = y/2Vr COS (bJslipt + (p + = ^ V d r C 0 S (0»liP + Y ) _ sin I 9sHp - — • 1/1 2 7 r sin I 8sHp + y (5.26) where UJSUp = cosf — cor is the relative angular velocity between the rotor and the reference dq axes. The magnitude and phase angle of the rotor-voltage vector with respect to the rotating dq reference frame are defined as Vr <p = tan' - i vqr Vdr (5.27) Similar to its presentation for the stator-side converter model, the rotor-side con-verter voltage transfer characteristic can be written as Vr 2V2 VDC (5.28) where Pmr the modulation depth of the rotor-side converter. Since the rotor-side converter is controlled to the desired value in the stator-flux-oriented dq reference frame, it is convenient to decompose the modulation depth into dq coordinates, to 91 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI obtain the following expressions: Pmrd Vdr = —7=VDC 7 (5-29) rmrq v = V f " D C The power-balance equation for the rotor-side converter (5.23) can be written in dq variables wi th appropriate scaling as Pconv = VDClDCr = 3 {vdridr + V V ) (5.30) B y combining (5.29) and (5.30), we obtain the following expression: 3 ^DCr = 2y^2 (^~>mr^dr "F -Prnrqiqr) (5.31) which, in combination wi th (5.19), gives a differential equation describing the D C link of the voltage converter in dq coordinates as C — — 2V^2^~>M^1 2 V*2~ ^^>mr^^r ^mrqiqr) (5.32) 5.6.4 Rotor-side Converter Control The P W M converter inserted in the rotor circuit allows for fast and flexible control of the D F I G by modifying the magnitude and phase angle of the rotor voltage. A s men-tioned in the previous section, the induction machine is controlled in a synchronously rotating dq reference frame wi th the d-axis aligned along the stator-flux vector po-sition. In this way, decoupled control between the electrical torque and the reactive power is obtained. The rotor-side converter is controlled by a two-stage controller. Stage 1 consists of very fast current controllers regulating the rotor d-axis and g-axis currents to their reference values. The reference of the quadrature-axis current component is specified by a slower speed controller (Stage 2) as shown in F i g . 5.7. The q component of the rotor current directly influences the torque, so iqr can be used for torque or active power control. The d-axis. component is a reactive current component, so idr can be used for reactive power or voltage control. Assuming that al l reactive power is supplied by the stator, i*dr may be set to zero. The outputs of the current controllers v'dr and v'qr, in a similar way to the stator-side converter controller, are used to calculate the reference values of v*dr and vqr of 92 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Rotor-side speed controller (Stage 2) + O-Rotor-side current controller (Stage 1) PI — • PI PI Figure 5.7: Rotor-side converter controller the rotor-side converter. The voltage equations of the rotor circuit in dq reference frame can be writ ten as _ n • , (r Ll \ d l d r (T LI \ . LM)J dt " " " " ^ L{ref)J^ J s ' " K " s ' (5.33) V = + ( L r - - J ^ % + uslip (Sntif + f^r , L l '•>s L ^ J dt S U P \L(ref)"ms <\~r ^ ( r e / ) where Rr, L r , L 0 and L s r e ^ are rotor-referred quantities related to per-phase resis-tances and reactances of an induction machine (Section 4.4.1.1) as — Ra L r = L\T 3 La 2 jiref) 3 L 2 ^ mr (5.34) 'Nr\2 and ims^ ^ is the rotor-referred R M S value of the stator magnetizing current, iml^ can be calculated from the following expression: ^=M%i ^ 93 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Equation (5.33) takes into account that in the stator-flux-oriented reference frame, the dq components of the stator flux are ipsd = I/Js = L0ims and tpsq = 0. The reference voltages of the rotor-side converter are calculated from (5.33) as Vdr = V-dr CO, Vqr = Vqr + "si -Aref) 0qr + L jjref) Idi (5.36) The reference voltages are converted back to phase domain using the inverse dq transformation wi th appropriate scaling as follows: yh. cos 9, COS l Usup 3 27T slip e 2*cos {9sup + — sin 9S sin (9shp sin 2TT 3 9slip 3 J J Jdr qr. (5.37) 5.7 Wind Turbine Model The kinetic energy of a mass of air m of velocity vwinci is given by a well-known equation: • . ^ = y ; vlind (5-38) The power of the moving air mass is the derivative of the kinetic energy with respect to time, and can be writ ten as dEk _ 1 dm 2 _ 1 2 / 5 3 9 ) •M) — — g Vwind ~~ 2 W i n d {O.OV) where q represents the mass flow given by the following expression: q = P -Vyjind • A (5.40) where • A is a cross section of the air flow • p is air density . Only a fraction of the wind kinetic energy can be converted into rotational power at the wind turbine shaft. This fraction of power (Pwind) depends on the wind velocity, rotor speed and blade position (if the turbine is pitch or active-stall controlled) and on the turbine design. It is usually quantified by a constant C p , representing aerodynamic 94 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI efficiency v = p (5.41) For a specific turbine design, Cp is normally calculated as a function of the pitch-angle P and the tip-speed ratio A. The tip-speed ratio is given as A = ^ (5.42) where R is the radius of the turbine blades and cot is the turbine speed. CP(X, P) is normally defined in the form of a two-dimensional lookup characteristic for different values of f3 and A. Alternatively, analytical approaches for approximating the Cp characteristic can be used, such as the following expression [39]: Cp = (0.44 - 0.01670) sin * £ X _ ~ ^ - 0.00184(A - 2)p (5.43) The mechanical power extracted from the wind (with A = irR2) is calculated from Pwind =\-p-n-R2- CP(X, P) • vlind . (5.44) with the mechanical torque driving the induction generator calculated from (5.6). 5.8 Maximum Power Tracking The family of PWind — w r curves can be derived from (5.44) for various values of wind velocity, with ur representing the shaft speed referred to the generator side of the gearbox. For the 7.5 k W wind turbine of [36], the curves are derived for fixed P, as shown in F i g . 5.8.. Curve Popt defines the maximum energy captured from the wind, and is characterized by the factor Kopt in the following expression: Popt = KoptL0Zr (5.45) The objective of the tracking control is to keep the turbine speed fixed on the Popt curve as the wind velocity varies. For wind velocity higher than the turbine's rating, the turbine energy captured has to be limited by applying pitch control or driving the machine to the stall point. One method of maximum power tracking is realized by observing the mechanical torque on the generator shaft, and driving the D F I G to the opt imum power curve by 95 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI setting the reference rotor speed of the rotor-side converter controller (Fig. 5.7) to T K, ( 5 . 4 6 ) opt This method of maximum power tracking is referred to as a speed-mode control. For extracted wind power higher than rated, the D F I G is driven to a stall point for the particular wind velocity, where C0„ In F ig . 5.8 the rated power of P m Q X = 7 . 5 k W defines the stall point for wind velocities higher than the rated (10 m/s) . 14000i 12000H 10000H g" 8000H 6000-4000-2000 500 1000 1500 2000 2500 3000 n (rpm) Figure 5.8: W i n d turbine characteristics wi th maximum energy capture curve Another method of maximum power tracking is a current-mode control where for the given shaft speed an electromagnetic torque is imposed on the D F I G (after compensating for friction losses) given from: T* = Koptio2 — Bcor by regulating the q-axis rotor current i*r to (5.48) 2 L ( r e / ) r , qr 3pfL2J (ref) ms (5.49) 9 6 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI The choice of a speed-mode or current-mode control is made by choosing the desired position of a switch in the rotor-side converter controller (Fig. 5.7). 5.9 Implementation of a DFIG Wind Turbine System in OVNI The doubly fed induction generator wind turbine system model consists of the phase-domain induction generator discrete model described in Section 4.4.1.3, the discrete phase-domain model of the voltage converter and the voltage-converter controllers operating in the dq reference frame. The phase-domain voltage converter, model is composed of the voltage transfer characteristics of the stator and rotor-side converters plus the differential equation of the D C link. The required values of the converter's terminal voltages are calculated as follows: p* p* * mdl n m<7l • rt vai = -^VDC COS 9S f~VDC s in 6S vbi = -^VDC COS \d3 - — \ ^-vDCsm 19, - — vcl = —[T~VDC cos I Os + y I 2 V D C S m V T p* p* (5.50) * mdr n mqr • n vA = —^-vDC cos esup 2 V d c S m s l i p * Pmdr fn 2T{\ Pmqr • f n 27T vB = —g— vDc cos I 6 slip ^-vDCsm I0slip - — * Pmdr (n M Pmqr . / „ 2TT vc = —^vDCcos I 9snp + - I ^-vDCsm I 9sHp + — where Pmdi, Pmqi, Pmdr and Pmqr are control variables of the voltage converter defined by (5.16), (5.28). 97 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI The D C link equation in the abc reference frame has the following form: c^- =ial (% cos es - % sin es dt V 2 2 M^-('.-S)-¥-(*-y)) +id ( ^ c o s I es + ^ ) - - ^ s i n (es + ^ •p! V p V 2 , V 3 / 7 (5.51) -ZB ^ COS ^ - ~ ^ Sin - f -<c (^f1 cos (eslip + *f) - ^ sin ^ + ^ and is discretized using the trapezoidal integration rule. These equations are com-bined with the equations describing the phase-domain induction machine presented in Section 4.4.1. 5.10 D F I G Wind Turbine Test System In this thesis, we have modelled a doubly fed induction generator wind turbine system from the literature using the Mult i level M A T E concept, and implemented it in O V N I . Reference [36] describes a doubly fed induction generator application to variable-speed wind-energy generation using back-to-back P W M converters in an experimental setup. The paper provides details on the control setup of the stator- and rotor-side converters and provides the data necessary for modelling of the wind turbine system in the phase domain. Model l ing details have been presented in the preceding sections of this chapter. We extended the D F I G wind turbine system case to include a step-up transformer and a 10 k V double-circuit transmission line connecting the wind turbine to a strong network modelled as an infinite bus. The single-line diagram of the test network is depicted in F i g . 5.9. For the purpose of comparison, the wind turbine generator defined in Appendix D is connected to a 250 V bus to replicate the situation in [36]. In addition, the inclusion of the connecting network allows us to study the response of the wind turbine to a three-phase fault applied in the middle of Circui t 2 of the 10 k V , 10 k m transmission line. We obtained the data for modelling the step-up transformer and the transmission line from typical values [40]. In the following section, simulation results are presented showing the response of 98 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 10 km line 0.25/10 kV / 7 W CCT1 CCT2 eR=1.6% ex=1.6% i Figure 5.9: D F I G wind turbine test case the D F I G wind turbine system to a change in wind velocity following the maximum speed tracking system described in Section 5.8. Subsequently, we show the responses of a D F I G system to a three-phase fault, and compare them against the results calculated wi th the Transient Stability Assessment Tool T S A T [41]. 5.11 Simulation Results 5.11.1 Transient Response of the DFIG Wind Turbine to a Step Decrease in Wind Velocity Initially, the induction machine is operating in steady state wi th a wind velocity of 9 m/s, corresponding to the optimal shaft speed of 1350 rpm and mechanical torque of 38 Nm. The maximum power captured at this wind speed is 5400 W. The maximum power-tracking system works in current-mode control. A t t = 2 s, the wind velocity decreases instantaneously to 5 m / s , which corre-sponds to the opt imal shaft speed of 750 rpm and mechanical torque of 12 Nm. The maximum power captured from the wind at this speed is 932 W. Reduced wind power, and therefore mechanical torque, cause the D F I G to decelerate, wi th the deceleration torque being the difference between the turbine's mechanical torque and the torque given by the opt imum curve (5.48). A similar situation wi th opposite effects would occur in the case of an increase in wind velocity. This scenario corresponds to the scenario described in [36] for the experimental setup. B y simulating the response of the D F I G model built wi th O V N I to the same type of disturbance, a comparison of the results can be obtained and conclusions can be drawn about the accuracy of the modelling approach. 99 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 0 .84 Figure 5.10: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) rotor angular velocity, (b) electromagnetic torque, (c) rotor speed, (d) stator flux In current-mode control, the D F I G reaches the new steady state after about sixty seconds. The D F I G is ini t ial ly running above the synchronous speed of 1000 rpm (rated frequency of the system is 50 Hz) for vWind — 9 m/s, and in the new steady state it runs below the synchronous speed for vWind = 5 m/s. The electromagnetic torque corresponding to the mechanical torque on the shaft after compensating for the friction losses decreases from the ini t ial 30 Nm to around 7.5 Nm, corresponding to the new steady state wi th a wind velocity of 5 m/s. A slight decrease of the stator flux reflects the decrease in voltage at the stator terminals. 100 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 4000 3000 <22000 1000 -2900 -2950 < > O-3000 -3050 -3100 4000 3000 -^ooo 1000 -2900 Figure 5.11: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) D F I G active power, (b) D F I G reactive power, (c) machine stator active power, (d) machine stator reactive power Initially, for a wind velocity of 9 m/s, the D F I G generates around 3800 W of active power and consumes around 3100 VAr of reactive power. Most of the active power of the D F I G is being delivered by the stator (3100 W), and the rest is delivered by the rotor v ia the converter. For the new steady state wi th wind speed of 5 m/s, the D F I G generates 500 W of active power and consumes about 200 VAr less of reactive power. In this latter case, the stator side delivers 720 W of active power and the rotor side consumes 220 W. We demonstrate here the capability of the D F I G to deliver active power both v ia the stator and the rotor of the machine. 101 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Figure 5.12: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator-side converter active power, (b) stator-side converter reactive power, (c) rotor-side converter (machine rotor) active power, (d) rotor-side converter (machine rotor) reactive power In F i g . 5.12(a), we notice the change of the stator-side converter active power from positive (generating) to negative (consuming) as the rotor speed moves from supersynchronous to subsynchronous due to the decrease in wind velocity. The stator-side converter active power depicted in F i g . 5.12(a) is equal to the rotor-side converter active power in F i g . 5.12(c) (with the converter losses neglected). The point where the active power of the converter is zero corresponds to the point where the rotor speed is equal to the synchronous speed of the induction generator. The reactive power of the stator-side converter is regulated to zero by the stator-side converter control (iqi = 0). 102 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) (b) 0 20 40 60 0 20 40 60 (c) (d) t (s) t (s) Figure 5.13: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) D F I G stator-side R M S voltage, (b) D F I G stator-side R M S current, (c) D F I G rotor-side R M S voltage, (d) D F I G rotor-side R M S current A decrease in the D F I G power output causes a decrease in the R M S values of the electrical variables at the D F I G terminals. For the R M S rotor voltage value depicted in F ig . 5.13(c), the minimum value corresponds to the time instant when the rotor speed is equal to the synchronous speed. 103 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 0 20 40 60 0 20 40 60 t (S) t (s) Figure 5.14: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator phase to neutral voltage, (b) stator phase current, (c) rotor phase to neutral voltage, (d) rotor phase current, (e) stator-side converter phase voltage, (f) stator-side converter phase current Figure 5.14 depicts the D F I G currents and voltages in the phase domain. F ig -ures 5.14(c) and F i g . 5.14(d) show the smooth operation of the D F I G through syn-chronous speed. B y examining F ig . 5.14(f) we note the change of "direction" of the stator-side converter current, indicating the change in converter power from genera-tion to consumption. 104 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI c - 200 CO CO C D « O > cr 1— o s 55 CO 0) C O eg "o > cr T 3 o CO CO 100 -100 300 200 100 -100 50 (a) V d s U S V qs (b) CO cu C D £ o > cr " D i _ O O rr -50 20 40 60 (c) V V %1 20 40 60 (e) / v J V 20 40 60 < ° CO 1-2 o cr TD o CO -4 -6 < •£ CD » o cr •o O CO CO -1 -2 _10 < cu 5 C J ^ 0 o rr d1 I qs — — — L _ _ 20 40 (d) 60 20 40 (f) 60 o. t(s) 20 40 t(s) 60 Figure 5.15: Transient response of the D F I G wind turbine to a step decrease in wind velocity: (a) stator voltages dq components, (b) stator currents dq components, (c) rotor voltages dq components, (d) rotor currents dq components, (e) stator-side converter voltages dq components, (f) stator-side converter currents dq components Figure 5.15 shows the dq components of the D F I G electrical quantities. Stator quantities are shown in the stator voltage-oriented reference frame, and rotor quan-tities are shown in the stator flux-oriented reference frame. The control of the q-axis stator-side converter current iq\ determines the displacement factor at the D F I G ter-minals, and it is kept at zero in this case (no reactive power is supplied to the network by the converter). For supersynchronous operation of the D F I G , this implies a phase displacement of 180° between the stator-side converter phase voltage and current, and for subsynchronous operation it gives a unity displacement factor. The control of the d-axis stator-side converter current i^i determines the active power exchanged with the network for the purpose of maintaining the converter D C - l i n k voltage at a 105 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI constant value. The induction machine is controlled through rotor-side control in the dq refer-ence frame, wi th the d-axis current component controlling the machine's reactive power, and the q-axis current component being proportional to the machine's elec-trical torque. Assuming that al l reactive power is supplied to the machine by the stator, the d-axis rotor current is set to zero. B y controlling the q-axis rotor current a reference torque or a reference speed can be imposed on the D F I G . simulation e^xperiment Figure 5.16: Comparison of a transient response of the D F I G wind turbine to a step decrease in wind velocity against the experimental results from the literature: (a) rotor speed, (b) d-axis rotor current, (c) q-axis rotor current Our simulation results compared against the experimental results published in [36] for a decrease of wind velocity are shown in F i g . 5.16. It can be seen from the figure that the results showed good agreement. Differences can be attributed to a number of factors, such as using an average model to model the P W M converter, neglecting saturation, neglecting details of the machine design contributing to the generator winding/slot harmonics, etc. 106 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 5.11.2 Transient Response of the DFIG Wind Turbine to a Three-phase Short Circuit Initially, the induction machine operates in steady state wi th a wind velocity of 8 m/s, generating 2700 W of active power that is transmitted to a "strong" power system modelled as an infinite bus. The D F I G stator-side converter is regulated to maintain its q-axis current at zero, meaning that the reactive power of the induction generator is entirely supplied from the network. A t t = 0.2 s, a three-phase fault is applied in the middle of circuit C C T 2 , as shown in F i g . 5.9. The fault is removed after 0 .1s without t r ipping the circuit. The response of the wind turbine system to the fault is depicted in F i g . 5.17 - F i g . 5.24. Large disturbances cause large ini t ial fault currents both in the stator and ro-tor. Moreover the D C link wi l l experience overvoltages due to reduced voltage at the stator-side converter terminals. High currents flowing through the converter may ac-tivate the crowbar protection in order to preserve the converter which may then lead to the wind turbine being disconnected from the network. It has been recognized in the literature that stability models of D F I G wind turbines may not adequately esti-mate the transients in the rotor currents [42], therefore they may not provide correct responses to system disturbances. To verify this statement, a transient response of a D F I G to a three-phase fault has been simulated wi th the Transient Security Assess-ment Tool ( T S A T ) [41], which is a good representative of commonly used stability simulation tools. The exact same case as tested with O V N I could not be replicated due to the design of T S A T [41], which allows accurate computations in the M V A range rather than in the k V A range. Therefore the machine and the network parameters were scaled to reflect the power flow at the M V A level. In particular, the D F I G machine generating 2.7 k W was scaled to generate 2.7 M W of active power and consume 3 M V A r of reactive power. Default controls of a D F I G in T S A T were modified to reflect similar conditions as in the O V N I model. The D F I G model used in T S A T is a user-defined model of type 3 that represents a rotor voltage-controlled D F I G wi th flux transients. The shape and general trends of al l the responses computed wi th T S A T correspond to those calculated with O V N I . The most significant difference in the results is found in the oscillatory behavior caused by the transient in the converter's D C link. As is commonly done in transient stability modelling, the stator-side converter and the D C link are neglected and only the rotor-side converter control is modelled. B y examining the results shown in Figs. 5.25 and 5.26, one can conclude that this is probably an oversimplification imposed on the model. This type of behavior, which is not well represented wi th stability models, is crucial for determining the correct responses of protective devices, such as the aforementioned crowbar protection. A more interesting comparison could be made between the stability model that includes the stator-side converter model and the 107 Chapter 5 . Modelling and Simulation of a DFIG Wind Turbine System in OVNI detailed model described in this work, both implemented in O V N I . This task, however, is left for future work. 5.11.3 Commentary on the Comparison of Transient Simulation Results Obtained with EMTP-type and Stability-type Tools The preceding section compares transient simulation results obtained with O V N I , which is an E M T P - t y p e of tool, and T S A T , which represents a typical stability sim-ulation tool. A similar comparison between any other E M T P - and stability-type tool wi l l have a similar outcome. It is very well-known that E M T P - t y p e of modelling is more accurate than stability modelling. Smaller scale systems, such as the D F I G W T G tested i n this work, can be suc-cessfully replicated in a few commercially available E M T P - t y p e simulators. However, . these simulators do not implement partitioning techniques. Our O V N I simulator implements the concept of M A T E (or, in the more general case, Mult i level M A T E described in Chapter 3 ) , which allows parallel processing and efficient computation at the subsystem level. Mul t i level M A T E can be applied to simulations of very large power system networks. Stabil i ty tools use approximate modelling to cope with the size of the network, while O V N I uses multilevel system parti t ioning and efficient modelling solutions to achieve this objective. Therefore, for the same objective of coping wi th system size, Mult i level M A T E allows more detailed modelling than do the stability tools. The solutions proposed in this thesis are intended to make O V N I a general purpose simulator capable of performing E M T P and stability simulations of large power systems. Other E M T P - t y p e tools are not oriented toward this purpose. 5.12 Summary In this chapter, a detailed example of modelling of a doubly fed induction generator wind turbine system is presented. Using our simulation tool we have replicated and simulated an experimental setup from the literature. Our comparison of responses to a step change in wind speed has shown good agreement wi th the experimental results. In addition, we have the same system simulated for a three-phase fault near the D F I G terminals. We have also made a qualitative comparison to a similar case setup using a commercially available transient stability simulation tool. 108 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Figure 5.17: Transient response of the D F I G wind turbine to a three-phase short circuit: (a) rotor angular velocity, (b) electromagnetic torque, (c) rotor speed, (d) stator flux 109 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI g-2500 (d) a- 2000 -1000 -1500 _-2000 < > -2500 O -3000 -3500 -4000 0.2 0.4 t(s) 0.6 0.8 Figure 5.18: Transient response of the D F I G wind turbine to a three-phase short circuit: (a) D F I G active power, (b) D F I G reactive power, (c) machine stator active power, (d) machine stator reactive power 110 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Figure 5.19: Transient response of the D F I G wind turbine to a three-phase short cir-cuit: (a) stator-side converter active power, (b) stator-side converter reactive power, (c) rotor-side converter (machine rotor) active power, (d) rotor-side converter (ma-chine rotor) reactive power 111 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) (b) 0 0.2 . 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 t (s) t (s) Figure 5.20: Transient response of the D F I G wind turbine to a three-phase short circuit: (a) D F I G stator-side R M S voltage, (b) D F I G stator-side R M S current, (c) D F I G rotor-side R M S voltage, (d) D F I G rotor-side R M S current 112 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Figure 5.21: Transient response of the D F I G wind turbine to a three-phase short circuit: (a) stator phase to neutral voltage, (b) stator phase current, (c) rotor phase to neutral voltage, (d) rotor phase current, (e) stator-side converter phase voltage, (f) stator-side converter phase current 113 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 700 650 600 o550 o 500 450 400 0.1 Iw I ll A A A / ^ ~ ~ 0.2 0.3 0.4 t(s) 0.5 0.6. 0.7 0.8 Figure 5.22: Transient response of the D F I G wind turbine D C - l i n k voltage to a three-phase short circuit 114 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) 200 r > 100 co > ° " •S > 0 •100 0 300 > 200 CT > 100 > 0 -100 C 50 > 0 CT > > -50 -100 C ds V qs 0.2 0.2 0.2 0.4 (C) 0.4 (e) 0.4 t(s) 0.6 0.6 0.6 0.8 Vd1 I AA. \s V . ' — q1 0.8 1 1 IP qr 0.8 0 < ~ 2 CO CT *CO . _ ° - 4 (b) 1 qs I A A / W - — r AAA,- • < 5 CT -P 0 0.2 0.4 t(s) 0.6 • 1 R • IAAA/NIAA/V 'dr 0.8 Figure 5.23: Transient response of the D F I G wind turbine to a three-phase short circuit: (a) stator voltages' dq components, (b) stator currents' dq components, (c) rotor voltages' dq components, (d) rotor currents' dq components, (e) stator-side converter voltages' dq components, (f) stator-side converter currents' dq components 115 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) Figure 5.24: Transient response of the D F I G wind turbine to a three-phase short circuit: rotor phase voltages and currents 116 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) Generator terminal voltage (pu) 1.000 • 44.800 44.960 Generator active power (MW) 2.800 4^.120 45.280 Time (sec) (C) 44.800 44.960 45.120 45.280 Time (sec) (b) 5.440 45.600 Generator speed (Hz) 53.990 • 44.800 44.9 \: : : V ^ — Generator reactive power (MVAR) -2.000 (d) 45.440 45.600 44.800 44.960 45.120 45.280 Time (sec) 45.120 45'.280 45".440 45.600 Time (sec) V : i T 1 1 45.440 45.600 Figure 5.25: Transient response of the D F I G wind turbine to a three-phase short circuit obtained wi th the stability tool T S A T : (a) D F I G stator terminal voltage, (b) D F I G rotor angular velocity, (c) D F I G active power, (d) D F I G reactive power 117 Chapter 5. Modelling and Simulation, of a DFIG Wind Turbine System in OVNI Output of WTGUDM block 0.300 • (a) 44.800 44.960 45.120 TUt Time (sec) Line current magnitude (pu) 0.042 (b) 45.440 45.600 / Time (sec) Figure 5.26: Transient response of the D F I G wind turbine to a three-phase short circuit obtained wi th the stability tool T S A T : (a) D F I G rotor dq voltages and currents (from top to bot tom - Vdr, idr, vqr and iqr), (b) D F I G stator current 118 Chapter 6 Conclusion This chapter summarizes the main contributions of the work presented in the preced-ing chapters and provides suggestions for the continuation of this research. 6.1 Summary of Contributions This thesis makes the following contributions to the field of power system transient simulations: • Development of the new Multilevel M A T E Concept. We have presented a new concept of multilevel system parti t ioning and func-tional modelling for the simulation of power system networks. The proposed concept by its nature allows a more efficient solution of the subsystems' equa-tions, and provides a means for general modelling of power system components with nodal and branch equations. The advantages of subsystem partitioning wi th the Mult i level M A T E concept introduced in this work are presented in terms of the computational speedup over the existing single-level M A T E . The greatest improvements are realized when simulating subsystems containing components of a changing nature. The advantages of modelling with branch equations are demonstrated through the implementations of an ideal switch and a nonlinear control function. W i t h the Mult i level M A T E concept, the subsystem's admittance matr ix is not affected by a changing topology due to switching, and its order is not increased by adding branch equations. W i t h this approach, branch currents can be calculated at every time step or only when needed. Nonlinear functions of state variables of the system require solutions with itera-tions at every time step of the simulation.. W i t h our new approach, the nonlinear equations iterate against the linear systems' Thevenin equivalents alone, which, in the M A T E approach, are calculated at the level of the links. • Modell ing of Power System Components with Multi level M A T E . M u l -tilevel M A T E provides a very general approach to modelling of voltage sources, grounded or ungrounded, by means of sublinks. A solution for conditioning 119 Chapter 6. Conclusion of near-singular or singular nodal equations is provided by the introduction of. compensating shunt impedances. Also , the solution for calculating ungrounded subsystem reference voltages is given. We describe branch (sublink) implementation of phase-domain induction and synchronous machine models wi th Mult i level M A T E . The advantage of branch implementation is evident from the fact that the changing nature of the ma-chine's coupled branch equations does not affect the subsystem's admittance matrix. Also , the transformation of coupled branch equations into nodal equa-tions is no longer required, which greatly simplifies the modelling approach. Mult i level M A T E enables the controller's equations (linear and nonlinear) to be conveniently incorporated. Depending on the scope of the controller, global and/or local control is easily implemented by means of links and sublinks. • Validation of Model l ing with Multi level M A T E through an Example of a Doubly-Fed Induction Generator W i n d Turbine System. We test our proposed approach to modelling and simulation on an example of a doubly fed induction generator wind turbine system. We describe modelling of the wind turbine system components and control in detail. We have replicated an experimental setup from the literature [36] in our simulation tool and tested it for two types of disturbances: decrease in wind velocity and a three-phase fault in the connecting double-circuit transmission line. The results were successfully compared against the results in the literature and against a traditional stability simulation tool. The comparison shows the advantages of using more detailed modelling, especially when control and protection devices play a major role in the system's response. 6.2 Recommendation for Future Research This thesis addresses a broad range of issues associated wi th the modelling and simula-tion of power system networks. Many aspects of this work can be further investigated and developed, including the following: 1. Diakoptics versus Sparsity? This question has been posed since the 1960s, and some of the reviewers of our recent paper have reinstated it . Sparsity techniques have not been used in this thesis, but we noted that the efficiency of the multilevel partit ioning approach depends on the sparse nature of the network. The actual comparison of the two approaches and the possibility of employing sparsity wi thin the M A T E concept is s t i l l to be investigated. 2. Opt imal network partitioning algorithm. Opt imal partit ioning that takes into account the multilevel approach needs to be developed. The algorithm has to take into account the particular system's topology and characteristics. 120 Chapter 6. Conclusion 3. Fixed point versus Newton-type iteration method? The iterative proce-dure for nonlinearities implemented in this work corresponds to the fixed point iteration method. Newton's method requires more processing time due to the calculation of the Jacobian matrix, but its convergence is faster. The possibil-ity of implementing Newton's method for the iteration of nonlinear equations wi thin the Mult i level M A T E concept needs to be further investigated. 4. Expansion of the D F I G wind turbine system model. The wind turbine system implemented and tested in O V N I does not include models of the pro-tection devices, for example crowbar protection. Blade pitch control and over and under speed trips should also be included. 5. Stability studies of integration of a D F I G wind turbine farm into the power system network. The wind turbine model described in this thesis could be used for stability analysis of large wind farms interconnected with the traditional power system network. It could also be used for benchmarking the approximate wind-farm models used in transient stability tools. 6. Implementation of latency and eigenvalue analysis into O V N I . Meth-ods have already been developed in our Power System Group to take advantage of the slow and fast nature of the subsystems as reflected in the size of the time step used for their integration, eventually resulting in even more efficient simulation speeds. Eigenvalue analysis provides a valuable insight for stability considerations of power system networks. These approaches need to be imple-mented into the current version of O V N I . 7. Development and implementation of power electronic models into O V N I . This thesis uses an average model of a voltage-source inverter, which is considered appropriate for stability analysis considerations. Detailed models, however, need to be developed and implemented in O V N I that take advantage of the multilevel partit ioning approach. The availability of detailed models is necessary for specific types of studies requiring highly detailed modelling or to validate approximate models and their application limitations. 8. Coupling of O V N I - N E T hardware solution to the new version of O V N I and benchmarking of a large power system network model. Fast simulation of large power system networks is attainable through a combi-nation of software and hardware solutions. The software solution for O V N I ' s implementation wi th Mult i level M A T E needs to be coupled wi th the hardware solution of O V N I - N E T and tested on an extensive example of a large power system network. 121 Bibliography [1] J . R. M a r t i , ' L . R . Linares, J . A . Hollman, and F . A . Moreira . O V N I : Integrated software/hardware solution for real-time simulation of large power systems. In Conference Proceedings of the 14th Power Systems Computation Conference, PSCC02, Seville, Spain, June 2002. [2] J . R . M a r t i and L . R . Linares. Real-time E M T P - b a s e d transients simulation. IEEE Transactions on Power Systems, 9(3):1309- 1317, August 1994. [3] S. Acevedo, L . R . Linares, J . R. M a r t i , and Y . Fujimoto. Efficient H V D C con-verter model for real time transients simulation. IEEE Transactions on Power Systems, 14(1):166-171, 1999. [4] H . W . Dommel . Digi ta l computer solution of electromagnetic transients in single and multiphase networks. IEEE Transactions on Power Apparatus and Systems, PAS-88(4):388- 399, A p r i l 1969. [5] H . W . Dommel and W . S. Meyer. Computat ion of electromagnetic transients. Proceedings of the IEEE, 62(7):983-993, 1974. . [6] N . Sato and W . F . Tinney. Techniques for exploiting the sparsity of the network admittance matrix. Proceedings of the IEEE, PAS-82:944-950, 1963. ' [7] J . R . M a r t i , L . R. Linares, J . Calvino, H . W . Dommel , and J . L i n . O V N I : A n object approach to real-time power system simulators. In Proceedings of the 1998 International Conference on Power System Technology, Powercon'98, pages 977-981, Beijing, China , August 1998. [8] L . R . Linares. OVNI: (Object Virtual Network Integrator) A Fast Algorithm for Simulation of Large Electric Networks in Real Time. P h D thesis, The University , of Br i t i sh Columbia , 2000. [9] J . A . Hol lman and J . R . M a r t i . Real time network simulation wi th PC-cluster. IEEE Transactions on Power Systems, 18(2):563-569, 2003. [10] T . De Rybe l , J . A . Hol lman, and J . R. M a r t i . O V N I - N E T : A flexible cluster interconnect for the new O V N I real-time simulator. In Conference Proceedings of the 15th Power Systems Computation Conference, PSCC05, Liege, Belgium, August 2005. 122 Bibliography 11] Gabriel K r o n . A set of principles to interconnect the solutions of physical systems. Journal of Applied Physics, 24(8):965-980, 1953. 12] H . Helmholtz. Uber einige Gesetze der Vertheilung elektrischer Strome in korperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche [Some laws concerning the distribution of electrical currents in conductors with applications to experiments on animal electricity]. Annalen der Physik und Chemie, 89(6):211-233, 1853. 13] L . Thevenin. Extension de la loi d 'Ohm aux circuits electromoteurs com-plexes [Extension of Ohm's law to complex electromotive circuits]. Annales Telegraphiques, 10:222-224, 1883. 14] L . Thevenin. Sur un nouveau theoreme d'electricite dynamique [On a new theo-remof dynamic electricity]. C. R. des Seances de lAcademie des Sciences, 97:159-161, 1883. [15] I E E E recommended practice for industrial and commercial power systems anal-ysis. IEEE Std 399-1990, 1990. [16] J . B . Ward and H . W . Hale. Digi ta l computer solution of power-flow problems. Transactions of AIEE, PAS-75:398-404, June 1956. [17] W . F . Tinney and J . W . Walker. Direct solutions of sparse network equations by optimally ordered triangular factorization. Proceedings of the IEEE, 55(11):1801-, 1809, 1967. [18] W . F . Tinney. Ci ta t ion classic - direct solutions of sparse network equations by optimally ordered triangular factorization. Current Contents - Engineering Technology and Applied Sciences, 18:12-12, 1979. [19] H . M . Markowitz . The elimination form of the inverse and its application to linear programming. Management Science, 3:255-269, 1957. [20] I. S. Duff. A survey of sparse matrix research. Proceedings of the IEEE, 65(4) :500-535, 1977. [21] I. Hajj . Sparsity considerations in network solution by tearing. IEEE Transac-tions on Circuits and Systems, 27(5):357-366, 1980. [22] Chung-Wen Ho, A . Ruehl i , and P. Brennan. The modified nodal approach to network analysis. IEEE Transactions on Circuits and Systems, 22(6):504-509, 1975. [23] J . R . M a r t i . Mult i -area Thevenin Equivalent ( M A T E ) . A circuit concept. Re-search note to H . W . Dommel and L . Linares, August 1993. 123 Bibliography [24] P. Zhang, J . R . M a r t i , and H . W . Dommel. Network parti t ioning for real-time power system simulation. In Conference Proceedings of the 6th International. Conference on Power Systems Transients (IPST'05), Montreal , Canada, June 2005. [25] J . R. M a r t i and K . W . Louie. A phase-domain synchronous generator model including saturation effects. IEEE Transactions on Power Systems, 12(1):222-229, 1997. [26] F . A . Moreira and J . R . M a r t i . Latency techniques for time-domain power system transients simulation. IEEE Transactions on Power Systems, 20(l):246-253, February 2005. [27] P. Kundur . Power Sysem Stability and Control. The E P R I Power System Engi-neering Series. M c G r a w - H i l l , Inc., 1994. [28] H . W . Dommel. EMTP Theory Book. Microtran Power System Analysis Corpo-ration, 2nd edition, 1996. [29] J . R. M a r t i and T . O. Myers. Phase-domain induction motor model for power system simulators. In Proceedings of the IEEE Communications, Power, and Computing Conference WESCANEX 95, 1995. [30] P. C . Krause, O. Wasynczuk, and S. D . Sudhoff. Analysis of Electric Machinery and Drive Systems. W i l e y - I E E E Press, 1st edition, 2002. [31] A . E . A . Araujo, H . W . Dommel, and J . R . M a r t i . Simultaneous solution of power and control systems equations. IEEE Transactions on Power Systems, 8(4):1483-1489, 1993. [32] B . D . Bonatto and H . W . Dommel. A circuit approach for the computer mod-elling of control transfer functions. In Conference Proceedings of the 14th Power Systems Computation Conference, PSCC02, Seville, Spain, June 2002. [33] L . Linares, M . Armstrong, and J . R. Mar t i . Software implementation of controller representation in the O V N I simulator. In Conference Proceedings of the 15th Power Systems Computation Conference, PSCC05, Liege, Belgium, August 2005. [34] J . Mahseredjian, L . Dube, M i n g Zou, S. Dennetiere, and G . Joos. Simultaneous solution of control system equations in E M T P . IEEE Transactions on Power Systems, 21(1):117-124, 2006. [35] L . Dube and H . W . Dommel. Simulation of control systems in an electromagnetic transients program wi th T A C S . In Proceedings of IEEE PICA Conference, pages 266-271, 1977. 124 Bibliography [36] R . Pena, J . C . Clare, and G . M . Asher. Doubly fed induction generator using back-to-back P W M converters and its application to variable-speed wind-energy generation. IEE Proceedings-Electric Power Applications, 143(3):231-241, May 1996. [37] J . Morren, S . W . H . de Haan, P. Bauer, J . T . G . Pierik, and J . Bozelie. Comparison of complete and reduced models of a wind turbine wi th doubly-fed induction generator. In Tenth European conference on power electronics and applications, EPE 2003, Toulouse, France, September 2003. [38] J . M . Aller , A . Bueno, and T . Paga. Power system analysis using space-vector transformation. IEEE Transactions on Power Systems, 17(4):957-965, 2002. [39] Seu l -Ki K i m , Eung-Sang K i m , Jae-Young Yoon , and Ho-Yong K i m . P S C A D / E M T D C based dynamic modeling and analysis of a variable speed wind turbine. In IEEE Power Engineering Society General Meeting, 2004. [40] D . Beeman. Industrial Power Systems Handbook. M c G r a w - H i l l Book Company Inc., 1st edition, 1955. [41] Powertech Labs Inc. TSAT - Transient Security Assessment Tool, Version 5.1, 2005. [42] V . Akhmatov. Analysis and dynamic behaviour of electric power systems with large amount of wind power. P h D thesis, Technical University of Denmark, 2003. [43] C . W . Taylor. Power System Voltage Stability. The E P R I Power System Engi-neering Series. M c G r a w - H i l l , Inc., 1994. [44] T . V . Cutsem and C . Vournas. Voltage Stability of Electric Power Systems. Power Electronics and Power Systems Series. Kluwer Academic Publishers, 1998. [45] P. M . Anderson and A . A . Fouad. Power System Control and Stability. I E E E Power Systems Engineering Series. Institute of Electr ical and Electronics Engi -neers, 1993. [46] S. B . L ippman. C++ Primer. Addison Wesley Longman, Inc., 2nd edition, 1996. [47] J . Arr i l laga and N . R . Watson. Computer Modelling of Electrical Power Systems. John Wi ley and Sons, L t d . , 2nd edition, 2001. [48] E . Clarke. Circuit Analysis of A-C Power Systems. Wiley, 1950. [49] M . Lukic . Feasibility of security assessment in power systems using full time-domain solutions in the E M T P . Master's thesis, The University of Bri t i sh Columbia, December 2000. 125 Bibliography [50] B. D. Bonatto. EMTP Modelling of Control and Power Electronic Devices for Electric Power Quality Assessment. PhD thesis, The University of British Columbia, 2001. [51] A. Petersson. Analysis, Modelling and Control of Doubly-Fed Induction Gener-ators for Wind Turbines. PhD thesis, Chalmers University of Technology, 2005. [52] Powertech Labs Inc. VSAT - Voltage Security Assessment Tool, Version 1.2, 2000. [53] I. A. Hiskens. Analysis tools for power systems - contending with nonlinearities. Proceedings of the IEEE, 83(11):1573- 1587, November 1995. [54] H. W. Dommel. Techniques for analyzing electromagnetic transients. Computer Applications in Power, IEEE, 10(3): 18-21, 1997. [55] R. Pena, J. C. Clare, and G. M . Asher. A doubly fed induction generator us-ing back-to-back P W M converters supplying an isolated load from a variable speed wind turbine. IEE Proceedings-Electric Power Applications, 143(5):380-387, September 1996. [56] M . Armstrong, J. R. Marti, L. R. Linares, and P. Kundur. Multilevel M A T E for efficient simultaneous solution of control systems and nonlinearities in the OVNI simulator. IEEE Transactions on Power Systems, 21(3):1250-1259, 2006. . [57] J. V. Mitsche. Stretching the limits of power system analysis. IEEE Computer Applications in Power, 6(1):16-21, 1993. [58] F. Wu. Solution of large-scale networks by tearing. IEEE Transactions on Circuits and Systems-, 23(12):706-713, 1976. [59] K. W. Chan, R. W. Dunn, and A. R. Daniels. Efficient heuristic partitioning algorithm for parallel processing of large power systems network equations. IEE Proceedings-Generation, Transmission and Distribution, 142(6):625-630, 1995. [60] L. Chua and Li-Kuan Chen. Diakoptic and generalized hybrid analysis. IEEE Transactions on Circuits and Systems, 23(12):694-705, 1976. [61] H. H. Ffapp. Diakoptics - the solution of system problems by tearing. Proceedings of the IEEE, 62(7):930-940, 1974. [62] N. Rabbat and Hsueh Ffsieh. A latent macromodular approach to large-scale sparse networks. IEEE Transactions on Circuits and Systems, 23(12):745-752, 1976. 126 Bibliography [63] L . M . Wedepohl and L . Jackson. Modified nodal analysis: an essential addition to electrical circuit theory and analysis. Engineering Science and Education Journal, l l (3) :84-92 , 2002. [64] Anonymous. Modified nodal analysis and the incorporation of multiphase, cou-pled networks. Engineering Science and Education Journal, l l (6):254-256, 2002. [65] R . J . Koessler, S. P i l lu t l a , L . H . Trinh, and D . L . Dickmander. Integration of large wind farms into uti l i ty grids (part 1 - modeling of D F I G ) . In IEEE Power Engineering Society General Meeting, 2003. [66] P. Pourbeik, R . J . Koessler, D . L . Dickmander, and W . Wong. Integration of large wind farms into uti l i ty grids (part 2 - performance issues). In IEEE Power Engineering Society General Meeting, 2003. [67] H . W . Dommel. State of the art of transient stability simulations for electric power systems. Unpublished report not publically availabe, 1976. [68] J . R . M a r t i . Synchronous generator equivalent circuit. Unpublished internal report, January 1996. [69] H . W . Dommel . Notes on power system analysis. Course notes for E E C E 459, 1975. [70] J .R . M a r t i . Network analysis and simulation. Course notes for E E C E 560, 1999. 127 Appendices 128 Appendix A Trapezoidal Integration Rule The application of the implicit trapezoidal integration rule in solving electrical net-works was first introduced by Dommel in [4]. To demonstrate the use of the trape-zoidal rule, the basic relationship of voltage and current for an inductance is consid-ered: Mt) = (A.1) Rewriting the equation and integrating both sides, we obtain J vL{t)dt = L j diL(t) (A.2) t-At ' tk-At where tk is a discrete time point and A t is a discretization time step. The integral on the right-hand side is trivial, but to calculate the integral on the left side we need to know the function vL(t) between the points vL(tk) and vr,(tk — t). The trapezoidal rule assumes that the function is linear between two consecutive discretization steps (Fig. A . l ) , therefore the integral can be calculated as tk f vL(tk) + vL(tk - At) . / vL(t)dt = area — —— • At (A.3) tk-At A similar approach is applied for the capacitance to obtain its discrete-time equiv-alent. 129 Appendix A. Trapezoidal Integration Rule Appendix B Reference Frame Transformations The introduction of reference frame theory in the analysis of electrical machines has proved useful not only for analysis but has also provided a powerful tool for the implementation of sophisticated control techniques. We present an overview of the most commonly used reference frame transformations. For control purposes, only symmetrical operation is considered. B . l Clarke Transformation Three-phase A C machines are conventionally modelled using phase-variable notation. However, for a three-phase, star-connected machine, the phase quantities are not independent variables but, for symmetrical operation, ia(t) + ib(t) + ic(t) = 0 va(t) + vb{t) + vc{t) = 0 (B.l) where the variables denote stator phase currents, voltages and flux linkages, respec-tively. As a result of this redundancy, it is possible to transform the phase-variable rep-resentation into an equivalent two-phase representation. The transformation from three-phase to two-phase quantities proposed by Clarke [48], is written in matrix form as 2 ' 1 _ i,p(t) _ ~ 3 0 c o s ( ^ '2j 3 COS ( f ) sin ( f ) ia{t) h(t) ic(t) (B.2) The transformation is equally valid for the voltages and flux linkages. The stator current space vector is defined as the complex quantity h(t) = isa(t) +jiafi(t) (B.3) 131 Appendix B. Reference Frame Transformations Equation (B.2) can be written in a more compact form by introducing a vector rotation factor, a — e JT", as follows: I s(i) = H [ia(t) + aib{t) + a\(t)] (B.4) The choice of the constant in the transformations (B.2), (B.4) is somewhat arbi-trary. The value of 2/3 has the advantage that magnitudes are preserved across the transformation. A sinusoidal phase current with a peak magnitude of Im w i l l produce a current space vector wi th a peak magnitude of Im. Since any digital current control scheme is designed to control current amplitude, this form has proved to be more suitable for control purposes. For the Clarke transformation, the inverse relationship can be written as ' ia(t) ' 2 ib(t) Q ic(t) O COS cos sin s i n ( f isp(t) (B.5) B.2 dq Transformation The stator current, voltage and flux-linkage space vectors are complex quantities defined in a reference frame whose real axis is fixed to the magnetic axis of stator winding a. The corresponding quantities defined for the rotor circuit of a three-phase A C machine are similarly defined in a reference frame fixed to the rotor. In the analysis of electrical machines, it is generally necessary to adopt a common reference frame for both rotor and stator. For this reason, a second transformation, known as dq transformation, is formulated that rotates space-vector quantities by a known angle 9. Whi l e in the Clarke transformation the axes of the space-vector plane are stationary, here a new reference frame is defined where the axes are made to rotate at the same rate as the angular frequency of the phase quantities. This results in stationary current, voltage and flux-linkage space vectors. The current space vector in the dq rotating reference frame can be written as hq = isd + jiSQ = he"36 (B.6) which can be writ ten in matr ix form as cos(0) sin(0) " ^sct Isq - sin(0) cos(0) _ (B.7) The real component of the current space vector in this new reference frame is 132 Appendix B. Reference Frame Transformations called the direct (d-) axis component, while the imaginary component is called the quadrature (q-) axis component. The main advantage of the dq transformation is the elimination of position dependency from the machine electrical variables. The inverse dq transformation can be written in matr ix form as ^sct " cos(0) - sin(0) ' isd sin(0) cos(0) 1sq (B.8) The dq transformation performed directly on the phase quantities has the following well-known form: isd 2 Isq ~ 3 cos (0) - sin (0) cos (0 - f ) - s i n ( 0 - f ) cos (0 - f ) - s i n ( 0 - f ) ib ir (B.9) 133 Appendix C Digital Pass-Band Filter for Filtering Stator-Flux To accurately calculate of the stator-flux vector position 9sf in Section 5.6.3, it is necessary to filter out D C offsets in the stator flux. Instantaneous values of stator flux linkages ipa, ipi, and ipc are filtered using the digital pass-band filter described here. A simple R L C parallel circuit in resonance may be used to eliminate the D C component of the stator-flux linkages. F i g . C . l depicts the circuit configuration for parallel resonance. Figure C . l : Circui t configuration for filtering 60 Hz component of ^ ( t ) The admittance of this circuit can be expressed as In the resonant condition, the susceptance is equal to zero, which leads to the following well-known expression for the resonant frequency: The resonant frequency is entirely specified by the choice of Lp and CF- The magnitude of the filtered flux linkage ^ a _ F w i l l vary as a function of the frequency of the input flux linkage \& a , as depicted in F i g . C.2. (C.l) 1 (C.2) CO = UQ = 134 Appendix C. Digital Pass-Band Filter for Filtering Stator-Flux l 0.707 (Oj G>o 0)2 • G)(rad/s) Figure C.2: Variat ion of filtered flux-linkage magnitude as a function of frequency The bandwidth of the resonant curve is noted as OJBW- B y forming the ratio of the resonant frequency to the bandwidth, we obtain a factor that is a measure of the selectivity or sharpness of tuning of the parallel R L C circuit. For a parallel resonant circuit, the quality factor is W i t h the resonant frequency set to 50 Hz (314 rad/s) for the particular example in Chapter 5 and wi th the desired quality factor, the values of RF, LF and CF can be determined from the above equations by freely choosing either one of them. The circuit in F i g . C . l can be discretized using the trapezoidal rule of integration. The equivalent discrete circuit is depicted in F i g . C.3 , where Qo = — LOQRFC'F (C.3) LOBW hLF{t) = - ^ ^ ( t - A t ) + hLF(t - At) hCF{t) = ^ aAt ~ A t ) - hCF(t - At) (CA) B y solving the circuit, the filtered stator-flux linkage of phase a is + hLF(t) + hCF{t) 1pa.F(t) (C.5) 135 Appendix C. Digital Pass-Band Filter for Filtering Stator^Flux Figure C.3: Equivalent circuit of the parallel R L C filter discretized wi th the trape-zoidal rule 136 Appendix D Test System Parameters D . l Parameters of the 2250 hp Induction Motor from Chapter 4.4.1 Rating: 2250 hp Line-to-line voltage: 2300 V Poles: 4 Speed: 1786 rpm Base torque: 8.9 • 10 3 Nm Rated current: 421.2 A Inertia: J = 63.87 kg m2 Rs = 0.029 ft R'r = 0.022 ft Xu = 0.226 ft X'lT = 0.226 ft Xm - 13.04 ft Table D . l : 2250 hp induction motor parameters 137 Appendix D. Test System Parameters D.2 Parameters of the 835 MVA Synchronous Generator from Chapter 4.4.2 Rating: .835 MVA Line-to-line voltage: 26 kV Poles: 2 Speed: 3600 rpm Combined inertia of generator and turbine: J = 0.0658 • 10 6 J s2 Rs = 0.00243 fi Xis = 0.1538 fi Xq = 1.457 fi Xd = 1.457 fi Rkql = 0.00144 fi R'fd = 0.00075 fi X'lkql = 0.6578 fi X'lfd = 0.1145 fi R'kq2 = 0.00681 fi R'kd = 0.01080 fi X'lko2 = 0.07602 fi X'lkd = 0.06577 fi Table D.2: 835 M V A steam turbine generator parameters 138 Appendix D. Test. System Parameters D.3 Parameters of the 7.5 kW Doubly Fed Induction Generator from Chapter 5.10 Rated power = 7.5 k W Stator voltage = 415 V Rotor voltage = 440 V Rated stator current = 19 A Rated rotor current = 11 A Pole pairs =3 Rated speed = 970 rmp Base frequency = 50 Hz Ns/Nr= 1.7 J = 7.5kgm2 stator connection = delta rotor connection = wye Machine resistances and inductances per phase: R s = 1.06 ft Rr = 0.80 ft L s = 0.0664 H L 0 = 0.0810 H L r = 0.0320 H Control parameters: Stator-side converter Rotor-side converter Kv=0.12 tfw=0.49 a„=0 .9248 a w=0.988 #,-=4.72 # i r = 2 0 Oi=0.96 a i r =0.985 Table D.3: 7.5 k W doubly fed induction generator wind turbine system parameters 139
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multilevel MATE algorithm for simulation of power system...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multilevel MATE algorithm for simulation of power system transients with the OVNI simulator Armstrong, Maz̆ana Lukić 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Multilevel MATE algorithm for simulation of power system transients with the OVNI simulator |
Creator |
Armstrong, Maz̆ana Lukić |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | The existing MATE (" Multi Area Thevenin Equivalent") concept implemented in the real-time simulator OVNI provides a means of partitioning large systems of equations into subsystems connected through links. The subsystems are solved independently, with the overall solution integrated at the level of the links. MATE partitioning enables fast simulation of power system networks in two distinct ways: first, it allows parallel processing in a multi-machine environment and second, it allows integration of different solution techniques for individual subsystems. In this thesis, we generalize the concept of MATE with the new Multilevel MATE concept, in which each subsystem becomes the basis for a new level of MATE partitioning. The new concept allows subsystems to be further partitioned into subsubsystems that, for example, are of a constant and changing nature. Partitioning at the subsystem level leads to higher overall solution efficiency. Furthermore, power system components can be described using Multilevel MATE by their nodal and/or branch equations on the subsystem level in OVNI. A three-phase induction machine model is an example of a power system component that is naturally described with branch equations. Multilevel MATE also provides a convenient framework for the incorporation of controllers as well as nonlinear elements. The software implementation of OVNI written in this work is based on the Multilevel MATE algorithm. Models of power system components have been implemented and tested using the newly developed concept. A test case describing a doubly-fed induction generator wind turbine system has been modelled and studied as a practical example of the new solution scheme's capabilities. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0100317 |
URI | http://hdl.handle.net/2429/30708 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-266724.pdf [ 9.08MB ]
- Metadata
- JSON: 831-1.0100317.json
- JSON-LD: 831-1.0100317-ld.json
- RDF/XML (Pretty): 831-1.0100317-rdf.xml
- RDF/JSON: 831-1.0100317-rdf.json
- Turtle: 831-1.0100317-turtle.txt
- N-Triples: 831-1.0100317-rdf-ntriples.txt
- Original Record: 831-1.0100317-source.json
- Full Text
- 831-1.0100317-fulltext.txt
- Citation
- 831-1.0100317.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0100317/manifest