Multilevel M A T E Algorithm for the Simulation of Power System Transients with the OVNI Simulator by Mazana Lukic Armstrong B . S c , 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 O F THE REQUIREMENTS FOR 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 T h e existing M A T E (" M u l t i A r e a Thevenin Equivalent") concept implemented i n the real-time simulator O V N I provides a means of partitioning large systems of equations into subsystems connected through links. The subsystems are solved independently, w i t h the overall solution integrated at the level of the links. M A T E partitioning enables fast simulation of power system networks i n two distinct ways: first, it allows parallel processing i n 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 w i t h the new Multilevel M A T E concept, i n which each subsystem becomes the basis for a new level of M A T E partitioning. T h e 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 M A T E by their nodal and/or branch equations on the subsystem level i n 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. M u l t i l e v e l M A T E also provides a convenient framework for the incorporation of controllers as well as nonlinear elements. T h e software implementation of O V N I written i n 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 w i n d turbine system has been modelled and studied as a practical example of the new solution scheme's capabilities. ii Table of Contents Abstract ii Table of Contents iii List of Tables vii List of Figures viii Acknowledgements xiii 1 2 Introduction 1 1.1 Background 1 1.2 M o t i v a t i o n and Objectives 2 1.3 Thesis Organization 1.4 Thesis Contributions 3 4 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 Modified N o d a l Analysis 10 2.3 3 • • Summary 12 Multilevel M A T E 13 3.1 Introduction 13 3.2 M A T E System Partitioning 14 iii Table of Contents 3.3 Multilevel M A T E Concept 3.3.1 16 E x a m p l e : A n Ideal Switch Implementation w i t h Multilevel M A T E 20 3.4 C o m p u t a t i o n a l Efficiency of Multilevel M A T E P a r t i t i o n i n g 3.5 Commentary on Computational Efficiency of the M u l t i l e v e l Partitioning 3.6 22 Solution of Nonlinearities within the Multilevel M A T E Concept 3.6.1 MATE 25 . . . E x a m p l e : A Nonlinear C o n t r o l Function Implementation w i t h Multilevel M A T E 4 28 3.7 Multilevel M A T E Block Diagram 31 3.8 Summary 32 Modelling with Multilevel M A T E 34 4.1 Introduction 4.2 M o d e l l i n g of Ideal Voltage Sources 35 4.2.1 36 4.3 4.4 4.5 34 Commentary on Modelling of Ideal Voltage Sources 36 4.3.1 Ungrounded Circuits 37 4.3.2 Ungrounded Subsystems 42 M o d e l l i n g of Electric Machines , 46 4.4.1 Phase-Domain Induction Machine M o d e l 47 4.4.2 Phase-Domain Synchronous Machine M o d e l 64 4.4.3 Commentary on Modelling of Electric Machines 68 M o d e l l i n g of Controllers 4.5.2 4.6 ; M o d e l l i n g of Ungrounded Circuits and Subsystems 4.5.1 5 25 68 E x a m p l e : Modelling of a Single-Machine Infinite-Bus System in dq Coordinates w i t h Multilevel M A T E 70 Commentary on the Simultaneous Solution of Controller E q u a tions 74 Summary 77 M o d e l l i n g and Simulation of a D F I G W i n d T u r b i n e 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 D o u b l y Fed Induction Generator W i n d Turbines 80 5.4 D F I G M o d e l Structure 80 5.5 D o u b l y Fed Induction Generator M o d e l 81 5.5.1 81 5.6 6 Two-mass Representation of a D F I G W T G Shaft Voltage Converter M o d e l and C o n t r o l 83 5.6.1 Stator-side Converter M o d e l 85 5.6.2 Stator-side Converter C o n t r o l 88 5.6.3 Rotor-side Converter M o d e l 89 5.6.4 Rotor-side Converter C o n t r o l 92 5.7 W i n d Turbine M o d e l 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 i n 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 Decrease i n W i n d Velocity 99 5.11.2 Transient Response of the D F I G W i n d Turbine to a Threephase Short C i r c u i t 107 5.11.3 Commentary on the Comparison of Transient Simulation Results Obtained w i t h E M T P - t y p e and Stability-type Tools . . . 108 5.12 Summary 108 Conclusion 119 6.1 Summary of Contributions 119 6.2 Recommendation for Future Research 120 Bibliography 122 Appendices 128 A Trapezoidal Integration R u l e 129 B Reference Frame Transformations 131 Table of Contents C D B.l Clarke Transformation 131 B.2 dq Transformation 132 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 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 M o t o r from Chapter 4.4.1 . . . 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 D o u b l y Fed Induction Generator from Chapter 5.10 vi 137 139 L i s t of T a b l e s 3.1 P a r t i t i o n i n g 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 vii List of Figures 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 M u l t i l e v e l M A T E concept 16 3.4 A n example system to demonstrate ideal switching using the Multilevel M A T E concept 20 3.5 M u l t i l e v e l M A T E vs. single-level M A T E computation speedup 24 3.6 A n example system w i t h a nonlinear control function 29 3.7 F l o w chart of iteration process for nonlinear controllers 32 3.8 Block diagram of Multilevel 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 M o d e l l i n g of an ungrounded transformer circuit w i t h the Multilevel M A T E approach 38 Introducing a common reference point to the ungrounded secondary side of an ideal transformer 39 Introducing a common reference point to the ungrounded network w i t h two ungrounded ideal voltage sources 40 4.4 4.5 4.6 . . . . Introducing a common reference point to the ungrounded network w i t h two grounded ideal voltage sources . . 42 4.7 Grounded system w i t h an ungrounded subsystem 43 4.8 R o t o r and stator circuits of an induction machine 48 4.9 Coupled inductors 49 4.10 Discrete phase domain induction machine electrical model viii 51 List of Figures 4.11 Electrical circuit configuration for testing the phase-domain induction machine model implemented w i t h Multilevel 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 i n load torque of a 2250 hp induction motor 61 4.18 Stator phase current during step changes i n load torque of a 2250 hp induction motor 62 4.19 Rotor phase current during step changes i n load torque of a 2250 hp induction motor 62 4.20 Electromagnetic torque during step changes i n load torque of a 2250 hp induction motor 63 4.21 Rotor speed during step changes i n load torque of a 2250 hp induction motor 63 4.22 R o t o r and stator circuits of a synchronous machine 64 4.23 Stator phase current during a step increase i n input torque of a 835 M V A steam turbine generator . . 65 4.24 F i e l d current during a step increase i n input torque of a 835 M V A steam turbine generator 65 4.25 Electromagnetic torque during a step increase i n input torque of a 835 M V A steam turbine generator 66 4.26 R o t o r speed during a step increase i n input torque of a 835 M V A steam turbine generator . 66 4.27 R o t o r angle during a step increase i n input torque of a 835 M V A steam turbine generator . 67 4.28 Torque-rotor angle characteristics during a step increase i n input torque of a 835 M V A steam turbine generator 67 4.29 Single-machine infinite-bus equivalent network 4.30 Thyristor excitation system w i t h A V R and P S S ix 71 . 71 List of Figures 4.31 Multilevel M A T E partitioning of the single machine infinite bus system 73 4.32 Transient response of the single-machine infinite-bus system w i t h manual control: (a) terminal voltage, (b) active power, (c) rotor angle . . 75 4.33 Transient response of the single-machine infinite-bus system w i t h 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 w i n d 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 R o t o r 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-fluxoriented dq reference frame 90 5.7 Rotor-side converter controller 93 5.8 W i n d turbine characteristics w i t h m a x i m u m energy capture curve . . 96 5.9 D F I G w i n d turbine test case 99 5.10 Transient response of the D F I G wind turbine to a step decrease i n w i n d 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 i n w i n d 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 i n w i n d 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 . . 102 5.13 Transient response of the D F I G wind turbine to a step decrease i n 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 i n w i n d 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 w i n d turbine to a step decrease i n w i n d 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 w i n d turbine to a step decrease i n w i n d 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 . . Ill 5.20 Transient response of the D F I G w i n d 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 xi List of Figures 5.23 Transient response of the D F I G wind turbine to circuit: (a) stator voltages' dq components, (b) components, (c) rotor voltages' dq components, dq components, (e) stator-side converter voltages' stator-side converter currents' dq components a three-phase short stator currents' dq (d) rotor currents' dq components, (f) 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 w i t h 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 w i t h the stability tool T S A T : (a) D F I G rotor dq voltages and currents (from top to bottom - Vd , idr, v and i ), (b) D F I G stator current 118 A.l Illustration of integration w i t h the trapezoidal rule 130 C.l C i r c u i t configuration for filtering 60 H z component o f i/» (t) 134 C.2 Variation of filtered flux-linkage magnitude as a function of frequency 135 C.3 Equivalent circuit of the parallel R L C filter discretized w i t h the trapezoidal rule 136 r qr qr . . a xii Acknowledgements I have been extremely lucky to meet and to get to know so many wonderful people since my arrival i n C a n a d a 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 . M a r t i , , whose extensive knowledge and brilliant ideas have inspired this work. I thank h i m for his continuous support and understanding, especially during the difficult times. I would also like to express my sincere gratitude to Prof. Hermann W . D o m m e l , 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 w i t h 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 D r . P r a b h a K u n d u r , who provided me w i t h 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. M a n y thanks to D r . A l i Moshref for a l l the time he spent helping me learn about stability simulations. ; A special place i n my heart belongs to my dear friends who shared the ups and downs of graduate student life w i t h me. M y special thanks goes to my dear friend Richard Rivas, for a l 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 B i b i a n , Benedito Bonatto, Jesus Calviho, Maninder Dhaliwal, Jorge H o l l m a n , A w a d Ibrahim, Hee-Sang K o , Daniel Lindemeyer, Fernando Moreira, D a i N a n , K e n W i c k s and many others i n the old group for their friendship and care. I would also like to thank M e l i h a Selak for a l l the help,she provided to me and my family along the years. In the new Power G r o u p , I would like to give special thanks to T o m de R y b e 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 a l l of my endeavors. They have inspired me throughout my life to excel i n everything I do. I would like to thank my wonderful husband G o r d o n and my son D a v i d for always xiii Acknowledgements being there for me, for their love and patience. Y o u fill my life w i t h 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: D r . M a r t i through a research assistantship; U B C through a University Graduate Fellowship and a L i T z e Fong Memorial Fellowship; and Powertech Labs through a research assistantship and several part-time contracts. , M a z a n a A r m s t r o n g , October 2006 xiv Chapter 1 Introduction 1.1 Background Object V i r t u a l Network Integrator ( O V N I ) 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]. Q u o t i n g from [5], electromagnetic transients "occur on a scale of microseconds (e.g., initial 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 i n transient stability studies. In the E M T P , basic lumped circuit components, such as inductors and capacitors, are described w i t h ordinary differential equations. Differential equations are discretized using the trapezoidal rule of integration ( A p p e n d i x A ) to obtain algebraic difference equations. Discrete time step sizes depend on the type of phenomena studied and the highest frequency expected i n the results. T h e y range from a few nanoseconds for very fast transients to an upper limit of about one millisecond dictated by the steady-state base system frequency of 50 or 60 H z . T h e E M T P methodology normally requires a three-phase system representation as opposed to the positive sequence single-phase representation used i n stability programs which assume that voltages and currents can be represented at the base system frequency as phasors. The system of equations solved i n the E M T P for a network w i t h 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 sources and history terms. 1 current 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 are used to reduce the number of operations. 1 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 i n creating the new O V N I simulator were aimed at "speeding-up" the E M T P simulation [2]. T h e solution in [2] offered a software approach to topological decoupling of power system networks, using the transmission line models i n 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 i n [8]. Finally, cluster-based hardware solutions for real-time simulations i n O V N I were proposed i n [9], [10]. 1.2 Motivation and Objectives A s stated i n 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 w i t h the E M T P at acceptable speeds i n an off-line environment, and (ii) to allow real-time testing of power system equipment i n a closed-loop environment. B y combining the concept of system partitioning (Diakoptics ) w i t h the concept of system reduction (Thevenin equivalents ), the M A T E ( M u l t i - A r e a Thevenin Equivalent) methodology has opened the door to new horizons of investigation of computational procedures and modelling techniques used i n simulating large power systems. 2 3 P r i o r 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 i n [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 i n m i n d , an object-oriented approach to software programming of O V N I was proposed i n [8]. The initiative for the work presented i n this thesis started w i t h the notion of building a general purpose simulator. The simulator would, along w i t h the hardware solutions for speed, provide engineers w i t h an all-in-one tool capable of performing Sparsity applications in power system analysis were first introduced by N. Sato and W. F. Tinney in 1963 [6]. The concept of Diakoptics was first introduced by Gabriel Kron in 1952 [11]. 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]. 1 2 3 2 Chapter 1. Introduction the simulation of electromagnetic ( E M T P - t y p e ) transients, as well as electromechanical (stability) transients, w i t h a wide range of available models of power system components. T h e model complexity would range from detailed models of, for example, 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 i n 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 w i n d 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 w i t h 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 partitioning and modelling of power system components. The solutions presented i n 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 i n use today, then concentrates on a literature overview of the methods used for solving large systems of equations associated w i t h the simulation of power system networks. We provide an overview of Sparsity, Diakoptics and Modified N o d a l Analysis, and describe their applications. In Chapter 3, we first review the concept of M A T E and introduce the new concept of Multilevel M A T E . Multilevel M A T E allows multilevel partitioning 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 general approach to modelling ideal voltage sources connected node-to-ground or nodeto-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 demonstrate 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 suitable 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 modelling of power system components with branch equations. • A n 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 proposed. 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 b o t h planning and everyday operation of electric power networks. In the early days of power system analysis, the need for computational aids i n analyzing power system networks led to the design of the first analog computer, as early as 1929 [15]. T h e 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 limited i n scope due to the limited capacity of the punch card computers used then. Large-scale digital computers became available i n the mid-1950s, and the success of load-flow programs led to the development of programs for short-circuit and stability analyses. Phenomena analyzed w i t h different simulation methods and tools range from slower widespread network events, such as voltage instability, to fast electromagnetic 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 P r o g r a m ( E M T P ) [4]. Power-flow calculations are used to determine the steady-state of the power generation/transmission system for the given loading condition. T h e power and voltage constraints make the load-flow problem nonlinear. T h e numerical solution therefore 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. T h e stability of power of synchronous machines. major concerns i n modern networks operating close systems was associated early on w i t h rotor-angle stability More recently, voltage stability has emerged as one of the power systems characterized by highly loaded transmission to their capacity limits. D y n a m i c 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. D y n a m i c simulations are currently an integral part of planning and operation studies of power system networks. Electromagnetic transients simulations are the most common methods for simulation of power system surges and fast-circuit transients. Electromagnetic transients are rarely an important factor i n system planning or operating decisions, but are essential factors i n the design process. These types of simulations require extensive modelling detail and small simulation time steps due to the small time constants associated w i t h 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 i n 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 w i t h nodal equations that are particulary suitable for digital simulation. Moreover, they a l l use similar programming techniques to efficiently handle large matrices. T h e following section discusses the historical development of these techniques and their applications. 2.2 Techniques for Simulation of Large Power System Networks N o d a l analysis implements Kirchhoff's current l a w , which states that the sum of the currents entering the node i n 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. T h e nodal approach has an advantage over the mesh approach derived from Kirchhoff's voltage law i n that it introduces fewer variables and equations i n the case of power system networks, therefore leading to smaller system matrices. The system of equations associated w i t h nodal analysis has the general form: 4 [Y] • [V] = [/] where • [Y] is a nodal admittance matrix, • [V] is a vector of nodal voltages, • [I] is a vector of nodal currents. Kirchhoff's laws were first announced by Gustav Robert Kirchhoff in 1845. 4 7 (2.1) 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 computationally, especially i n the past when computer power was limited and expensive compared to today. Practically from the beginning of the digital-computer simulation era, power systems engineers have been developing techniques to enable faster solutions of large systems of equations. T w o main very distinctive streams established themselves i n the 1960s, Diakoptics and Sparsity Techniques. Sparsity techniques, highly favoured over Diakoptics, have developed roots i n the majority of the simulation tools i n existence today. To overcome the limitations of nodal analysis i n the modelling of some system components, modified nodal analysis ( M N A ) was formulated i n the mid-1970s. Because of its relevance to this work, we w i l l discuss M N A i n the last part of this section. 2.2.1 Sparsity In 1978, W . F . Tinney, referring to his previous papers [6] and [17] wrote i n a commentary: "The a i m 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 a l l , of them to become non-zero i n the solution process... Today, awareness of sparse matrix methods is practically inescapable. Papers, books, and symposia on sparsity abound and applications exist i n almost every field. B u t only a decade ago, when our paper was published, the concepts of sparsity were virtually 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 i n Management Science [19]. The first sparsity method applied to the solution of power system networks was developed by Sato and T i n n e y i n 1961, and published i n 1963 [6]. T h e work i n these papers was not recognized or accepted after published, and remained unknown until the m i d 1960s. In 1967, when the second paper was published by T i n n e y and Walker [17], sparsity was finally acknowledged, and has since become widely accepted and applied i n power system simulation tools. The original method proposed i n the 1960s remains valid today, w i t h most improvements consisting of marginally more efficient algorithms. B y 1977 when a survey paper on sparse m a t r i x 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 m a t r i x " . T h e task of solving the system for nodal voltages involves obtaining the inverse of the m a t r i x 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 m a t r i x can be expressed as a product of sparse m a t r i x factors, which brings an advantage i n computational speed, storage requirements and accuracy. T h e 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 i n sparse systems is information storage and retrieval. T h i s aspect alone allowed simulation of large systems w i t h more accurate models. T h e computation procedures are organized to avoid operations w i t h 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 initial sparsity as much as possible. The m a i n 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 i n 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 m a t r i x elements, tables of indexing information are required to identify the elements and facilitate their addressing. T h i s procedure, i f not properly planned, may diminish the time-saving advantages of sparsity methods. T h e 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 i n 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. D u r i n g the elimination, no operation is performed that would lead to a predictable zero result, and no memory is allocated for a predictable zero element. T h e original matrix, the table of factors, and a l l indexing tables contain only non-zero elements. Rows from the original m a t r i x are transferred to a compact working row i n which elements to the left of the diagonal are eliminated by appropriate linear combination w i t h previously processed rows from the partially completed table of factors. W h e n work on the row has been completed, it is added to the table of factors. For nonsymmetric 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 multiplying 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. T h e concept of Diakoptics was first proposed by Gabriel K r o n in 1953 [11], i n the Journal of Applied Physics. H i s sixteen-page paper describes the advantages of system tearing and gives hints about a multi-machine solution of largesystem m a t r i x equations. It also points out a reduction i n the computation time of matrix inversion on a single machine simply due to partitioning. H i s brilliant ideas, which are i n 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. K r o n ' s Diakoptics, nevertheless, was not completely forgotten. A number of publications expanded on his ideas for applications i n analysis of power systems, electrical networks and structural and other large-scale systems. In K r o n ' s original form, D i akoptics required the explicit computation of inverses of submatrices representing the individual subnetworks. A s a result, sparsity of subnetwork equations could not be exploited. A number of works published i n the 1970s, just after sparsity established itself, investigated network analysis by tearing using sparse m a t r i x solution techniques. A good overview of these attempts is given i n [21]. To this day, however, Diakoptics, however, has been questioned as to whether more efficient computations are possible w i t h the use of sparsity and solving the network as a whole. T h e 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 i n 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 i n finding the network solution. Possibly the m a i n reason why Diakoptics never developed roots i n the simulation of large power system networks was due to its ambiguous description i n the many papers that could not deliver clear application examples and justify its use. Moreover, a number of papers, as stated i n [20], claimed that Diakoptics tearing could never be more advantageous than sparsity techniques. 2.2.3 Modified Nodal Analysis A s stated i n [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 nonnatural elements are removed in order to write a conventional set of nodal equations, and are then reintroduced with branch equations. • a « b Z > - - -—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 (V ) and series impedance (Z ) 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, V and Vj,. Second, the new unknown current, I , 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: s s a s YR +1 - 1 ' v' v = +1 -1 -Z a h b (2.2) Is \ s ' la ' where I and lb are the nodal injection currents of the system. a 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 Z is zero, conventional nodal analysis would not produce a solution because of an undefined l/Z 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 s s 11 Chapter 2. An Overview of Power System Simulation Tools and Techniques an accuracy and efficiency point of view. T h e M N A equations are ordered i n 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 connected to one or more voltage sources. T h e node equation is always interchanged w i t h one of the branch equations having the smallest number of non-zeros to minimize fill-ins. For example, i n (2.2) the bottom equation would be interchanged either w i t h the node equation for V or H , depending on their number of non-zero elements dictated by the composition of the adjacent network. a Modified nodal analysis combined w i t h sparsity techniques is commonly used i n today'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. A d o p t i n g M N A results i n an increased order of the system of equations. These branch equations are generally highly sparse and can be efficiently dealt w i t h using modern sparsity-based computing packages. 2.3 Summary T h i s 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 Multilevel M A T E . Multilevel 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 m a t r i x 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. Multilevel M A T E also allows for convenient modelling of power system components whose function is better described by branch equations, as opposed to the nodal equations typically used i n 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 partitioning and branch (functional) modelling bring a significant improvement to the simulation efficiency of the O V N I code. Multilevel M A T E provides an efficient framework for implementing nonlinearities that require an iterative procedure. Power system components that can be conveniently modelled w i t h i n the Multilevel 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 T h e M u l t i - A r e a Thevenin Equivalent ( M A T E ) concept was first introduced by M a r t i [23] i n an internal report to his research group. Reference [1] presents a detailed explanation of the algorithm M A T E that extends the m a i n ideas of Diakoptics by recognizing that the subsystems split by the branch links can be represented by Thevenin Equivalents. T o 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 i n F i g . 3.1 is partitioned into three subsystems A, B and C, by introducing six links. T h e hybrid system of modified nodal equations has the following form: 'A 0 0 0 B 0 0 0 c P 9* l t r p' Q r —z ' v A ' VB v c h A h h B (3-1) c where: • [A], [B] and [C] are the subsystems' admittance matrices •• \p]i [Q] • b*L a d [r] are the subsystems' current injection (or connectivity) arrays n a n d [ ] are transpose arrays of p, q and r rt • [z] is a m a t r i x of the links' Thevenin impedances 14 Chapter 3. • IAB]> [^cl a r e Multilevel MATE vectors of the subsystems' accumulated currents • [V ] is a vector of the links' Thevenin voltages • [VA], [VB], [VC] are the subsystems' nodal voltages • [i ] is a vector of the links currents. a a 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 " i n the connectivity array of the "from" subsystem, and the "to" node is assigned a i n the connectivity .array of the "to" subsystem. A n y component (or branch) i n the system can become a link described by its T h e v e n i n impedance ([z]) and its T h e v e n i n voltage ([V^]). The most general representation of a link is depicted i n F i g . 3.2. Note that [z] and/or [V ] can be zero. a "from" "to" z • Figure 3.2: General description of a link Partitioned matrices i n (3.1) are manipulated to obtain the subsystems' Thevenin equivalents as seen from the linking nodes [1]. T h e result of the manipulation is the M A T E system of equations, i n the following form: I 0 0 a where eA VA 0 / 0 6 VB 0 0 I c 0 0 0 z (3.2) vc i e a a a a = A- l P b= B = A~ lhA &A -- c = C = e c -- e B t c- h l c t •T?e + q eB t a = B~ hB l z pa + qb + rc + z l a e A + T ec + V t tt and [/] is the identity matrix. The individuality of each subsystem is preserved i n (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 linking nodes. T h e 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 i n the form of injected link currents at the linking nodes: i = z e a The subsystems' nodal equations at the subsystem level as: v v v A B c 3.3 a a are then solved independently from each other -e = e - e A B c a-i - b• i c-i a a a Multilevel M A T E Concept The concept of M u l t i l e v e l M A T E introduced i n this thesis can be demonstrated on the same system depicted i n F i g . 3.1. In this case, the system is further divided into ten sub-subsystems w i t h a total of six links and twenty sublinks, as shown i n F i g . 3.3. T h e complete system of equations for the depicted system is shown i n (3.3). Matrices [A], [B] and [C] represent the conductance matrices of the corresponding subsystems, while \p ], [<7B], [rc] are the connectivity arrays of the sub-subsystems (Ai, • • • Bi, B2, • • • Ci, C2) w i t h i n the corresponding subsystems. A Subsystem B Subsystem C Figure 3.3: A n example of a system that demonstrates the Multilevel M A T E concept 16 Chapter 3. Multilevel MATE The hybrid system of equations describing the system shown i n F i g . 3 . 3 has the following form: A p PA Z A 0 A 0 B 0 QB B 0 0 A VAsublink h v B B = 1B sublink r rc -zc c [q* 0] rj 0 -Z C VA 1Asublink q 0 h P VBsublink v he c 0 }Csublink_ —z i _ Vc_sublink a (3.3) where 0 [A] A, 0 0 0 A 2 0 0 0 A [B]=3 0 0 0 0 B 0 0 0 0 0 B 0 0 0 0 0 B 0 0 0 0 0 B 2 z A Ci [C] 0 0 c 2 b and PAX ' PA = PA2 P = A P2 VA = Pz PAZ ' hX ' ' V \. ' ' Pi ' VA2 VA3 A h A = h3 A for subsystem A, w i t h 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 T h e v e n i n voltages. T h e dimension of vectors [iAsubUnk] and [V]4_ MtTik] i that of the number of sublinks i n subsystem A. ) s SU To solve the M u l t i l e v e l 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 0 A~ "PA l p A- p 0 + z l A A 0 P\A~ P ' B- q ' q B- q X A X 0 0 0 q'gB-iqB 0 + Zg l B ' c-v ' 0 [p* o] 0] [q* r C- r 0] [r* x c v A p A~ h ^Asublink t + V _s biink X A A A U B~ h l VB B ^Bsublink Q B~ h + VejuftJinA: r C" h + Vc_sublink l vc B }C sublink B X c C -v (3.4) n Note, that the sublink currents for subsystems A , B and C , respectively, are: (p A~ p t 1 A {q B- q •(r C~ r B t A A A • i •• p\A a +.z ) • iBsubiink + Q B~ q • i +.z ) • icsubiink + r C~ r • i B c B B h a + l A -- q B l l c + Z ) • i _sublink + P A~ p l A l B h B V _ blink A B c r C- h a +V l c c (3-5) sublink l c su +V c sublink Before applying the M A T E concept to the entire system, we note that the sublink currents at the subsystem level i n (3.3) do not need to be known at the system level.. T h e 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. G o i n g back to E q u a t i o n (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 i n (3.6). 0 0 0 0 c 0 0 pt A B q* t r p q r —z v v A h B h vc he — C • ic.sublink A B p A ~ q B r " iAsublink • i _sublink B (3-6) Equations (3.5) and (3.6) fully describe the system. E q u a t i o n (3.6) closely resembles 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 0 0 0 I 0 0 I 0 0 a VA e b VB e &Bsublink c Vc ec &C sublink z e Asublink A B l a (3.7) ^oisublink a where e Asublink — A p ' iAsublink A &B sublink = B~ qg • IB sublink 1 CC sublink = C~ rc • %C sublink l ' P ^-Asublink "t" Q ^Bsublink "T" f &C'sublink ^•asublink Vectors [e _ nk], [eBsublink] and [e subiink] represent the contributions of the sublink 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 linking nodes, as if the sublinks were not present i n the system. T h e Modified Thevenin Equivalent that includes the sublinks' contributions is represented by [CLMTE] and [e _ ]. A subH C A MTE The system of equations for the Modified Thevenin Equivalents is: 7 0 0 10 0 0 0 0 a 0 E VA &A-MTE b TE VB Z-B-MTE 1 C TE vc eC-MTE 0 Z_ M T M M a (3.9) e _MTE MTE A where O-MTE = a — A a e-A.MTE = e bMTE — b - Ab e-B-MTE = £B — AeB CMTE = ec-MTE = e C- AC Z .MTE = P a TE Za-MTE = P^A-MTE + Q^MTE l M a A a = A~ p + Q £B.MTE l A Ae A = A p {p\A p x A l A — Ae A -A e (3.10) c + * + eC-MTE rt + V a t + z ) - p\-a 1 {p\A~ p l c + ^CMTE t and A A A +z) A (p • e + V _ A A A ) (3-11) subHTlk for subsystem A, w i t h similar expressions for subsystems B and C. A t every time step, each subsystem w i l l pass its modified Thevenin impedances and voltages (i.e., [CIMTE] and [e _ TE] for subsystem A) to the links. T h e system of links is solved first, and the link currents are returned t o their corresponding subsystems to obtain the solution for the subsystem voltages. In the case where subsystem A M 19 Chapter 3. Multilevel MATE links are introduced for the purpose of subsystem partitioning (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 i n the case of switches), these sublinks need to be solved at every time step. In the following example, an ideal switch implementation w i t h the Multilevel M A T E approach is presented. 3.3.1 Example: An Ideal Switch Implementation with Multilevel MATE Multilevel M A T E allows switches to be modelled as sublinks w i t h i n each subsystem. Switching operations involve a change of the subsystem's topology. However, this change can be confined to terms i n (3.11) that, i n 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 | t gA12=2 S \ , -AAAM^—* . 1 / 5 2 gA34=2 S A/vV^ . 4 R =l CI . link , gB56=4 s r v V v v l — - — W \ A ilink(t) gA40=l S ,(t)=l A (f) / Sub-subsystem A ^ Subsystem A i—r / \ 6 s> 1 rf) i (t)=i A 6 Subsystem B Figure 3.4: A n example system to demonstrate ideal switching using the Multilevel 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 w i t h i n subsystem A and its branch equation is treated as a sublink. T h e h y b r i d system of equations is written from (3.3) i n general terms as: 0 A A X 0 PAX PA2 2 -ZA PA [P* 0 0 B 1 o] q l Px P2 0 VAX tlAX VA2 h<A2 ^Asublink ^Asublink v i q —z B a h B -v a (3.12) 20 Chapter 3. Multilevel MATE For the example system it is,written as: 3 -2 -2 2 0 0 0 0 0 2 -2 0 1 -1 0 0 0 -2 0 1 -1 3 0 0 0 0 0 0 1 0 0 [ 0 1 | o] V2 V3 V 4 0 -1 0 -1 4 -4 -4 5 [ - 1 0 ] 0 1 0 0 0 Vl 0 ' 0' 1 "^switch V : 5 . *. V 0 ^link W i t h the switch taken out of the system, we apply M A T E as i n (3.2), and obtain the subsystems' Thevenin equivalents as seen from the link nodes: z = p a + q b + z = 1 + 1.25 + 1 = 3.25 [fi] e = p e + q e + V = 0 - 1 + 0 = - 1 [V] l l a l a (3.13) l A B a where A7 A'2 P2 l Pl a= 1 b = B- q l 1 A7 h Al A^h A2 l e - B h l B B The modified T h e v e n i n equivalent is obtained for the case where the switch is closed from (3.10) a n d (3.11). B y referring to (3.12), the case where a switch i n subsystem A is closed implies that the corresponding diagonal element of the [ZA] matrix becomes zero (R itch = 0)- For an ideal switch opening, the corresponding column of the subsystem's connectivity array [p^] a n d the row of \p ] are set to zero, a n d a diagonal of [z ] is set to one (iswitch = 0). T h e change of topology due to switching w i l l modify the A terms i n (3.11) or, i n a more general case w i t h nonlinearities, the A terms i n (3.17). For example, when the switch is closed we obtain: sw A A z*MTE e .MTE a p*Aa) + 9*6 + z.= (l- 0.33) + 1.25 + 1 = 2.92 [fi] = (p*e - P * A e ) + q e + V = (0 + 0.33) - 1 + 0 = - 0 . 6 7 [V] = (p a l l A A B a W h e n 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. M o d e l l i n g a switch using the Multilevel M A T E concept has several advantages over the previous approach used i n O V N I [8]. T h e subsystem admittance m a t r i x (subsystem topology) is not affected by switching, therefore i t does not need to be 21 Chapter 3. Multilevel MATE retriangularized a n d / o r prestored. The size of the subsystem m a t r i x is not increased by adding branch equations that are solved simultaneously w i t h nodal equations, as is the case w i t h the modified nodal analysis. M N A not only increases the size of the subsystem matrix, but also mixes nodal and branch equations. Finally, w i t h the new approach the switch current can be calculated at every time step, or only when needed. T h e case of an ideal switch implementation w i t h Multilevel M A T E can be extended to any element that contributes to a subsystem w i t h 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 i n 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 partitioning 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]. E a c h subsystem and the system of link (branch) equations are solved on separate P C clusters that communicate w i t h each other at a predefined rate. H i g h M A T E partitioning of the network produces smaller-sized subsystems that are solved more efficiently (with a higher hardware price), but which also results i n an increased number of links that have the opposite effect on the simulation speed. T h e size of the system of links, because of its communication requirements, is a l i m i t i n g factor i n achieving real-time simulation speeds w i t h O V N I [10]. T h e a i m of M A T E partitioning is to attain an o p t i m a l compromise between m a x i m i z i n g the computation speed of the subsystems and minimizing the communication w i t h the links system by keeping its size reasonably small. O p t i m i z i n g system partitioning is beyond the scope of this, thesis. However, papers dealing w i t h the subject have been published since the early days of Diakoptics, w i t h the most recent one published i n conjunction w i t h the development of O V N I [24]. Multilevel 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 i n general correspond to o p t i m a l 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 i n subsystems. T h e M u l t i l e v e l M A T E concept offers an efficient solution for dealing w i t h 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. Multilevel 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 partition a subsystem, and functional sublinks, which come from the modelling requirements of certain power system components. Partitioning 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 w i t h a changing nature that needs to be updated at every simulation time step is a phase-domain induction or synchronous machine model [25]. M u l t i l e v e l 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 w i 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 w i l l be briefly demonstrated on the example shown i n F i g . 3.1 and F i g . 3.3. Figure 3.1 depicts a system separated into three subsystems and six links based, for example, on some optimal partitioning 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 w i t h 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 i n simulation speed between the single level and multilevel approach w i l l occur inside the subsystems. We have made several assumptions for the purposes of comparing the computational efficiencies of the solution of non-partitioned and partitioned subsystems. T h e first assumption is that subsystems A, B and C have no functional sublinks. T h e sublinks that are introduced are only for partitioning purposes. T h e second assumption 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. T h i s situation is the most computationally 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. T h e result of our comparison of the computational efficiency of the single-level and multilevel approaches i n floating point operations (flops) are shown i n F i g . 3.5. Higher partitioning of subsystems leads to faster computational time (comparison between subsystems A, B and C i n 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 B 3 5 2 6 13 1 C Table 3.1: P a r t i t i o n i n g of subsystems -•-A 25 , 0 -s-B - A - C : ! 90 ,- -r— 900 1800 3600 Number of Nodes Figure 3.5: M u l t i l e v e l 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, i n the case where subsystems A, B and C have ninety nodes each, the o p t i m a l partitioning among the three subsystems is the partitioning of subsystem A w i t h three sub-subsystems of thirty nodes and six sublinks. T h e speedup over the single-level M A T E w i t h no subsystem partitioning 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 i n the computational speed over single-level M A T E is greater for large networks w i t h many nodes (comparison between different number of nodes for individual subsystems i n F i g . 3.5 along x-axis). Multilevel M A T E reduces the number of computations over twenty times i n 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 r u n time, and the Multilevel M A T E partitioning would show an even greater increase i n calculation speed. If sub-subsystems Bi, B , B and 23 were constant and B was the only subsystem of changing nature, Multilevel vs. single-level M A T E would result i n a computational speedup of over a hundred times for the case of 3600 nodes. 2 3 4 24 5 Chapter 3. Multilevel MATE In general, for subsystems comprised of both changing and constant nature elements and for subsystems that are sparse, the application of the Multilevel M A T E concept w i 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 w i t h changing topology, each topological change (e.g. switching) w i t h i n a particular region of the subsystem would require the retriangularization of the entire subsystem matrix. W i t h the M u l t i l e v e l M A T E concept, the change of topology is confined to the sublink equations, therefore the subsystem's m a t r i x 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 partitioning 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 Multilevel M A T E concept, therefore, results i n further computational speedup through subsystem partitioning a n d / o r 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 i n the same subsystem. It could also allow the integration of phasor and phase-domain modelling. E x p l o i t i n g these aspects would result i n 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 w i t h certain power system components; for example, an under-load tap changing ( U L T C ) transformer w i l l introduce nonlinear behavior w i t h its changing turns ratio. T h e iron core an electric machine is subject to saturation, which makes the i n p u t / o u t p u t characteristic of the machine nonlinear. In controllers, the limiters, the function of multiplication, etc., will also result i n nonlinear system behavior. Nonlinearities i n 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). Nonlinear elements that cause the coefficients of difference equations to change w i t h 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 i n the change of coefficients i n the system's matrix. T h e change of a system's coefficients w i t h 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. T h e nonlinear torque equation i n the modelling of a synchronous or induction machine also belongs to this type of nonlinearity. T h e simultaneous solution of the system's equations and the nonlinearities that operate on state variables requires an iterative procedure. T h e basic principle for solving nonlinearities requiring iterations w i t h i n the M u l tilevel M A T E framework is explained i n general terms i n 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 Multilevel 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 i n (3.14). •A A PA 0 ~A Z jP Ajnonlin ^ [p t l_ PA-nonlin P u 0 ~ Ajnonlin _ Z o. 0 B 0] j q l — ' P ' VA 0 _0 _ 1 Asublink h ^Asublink = 1> Ajnonlin ^Ajnonlin h q —z i I t - B -v a — A a —I I— _ i (3.14) The nonlinear subsystem A equations are composed of: • nodal equations described by the conductance m a t r i x [A], the nodal voltages [VA], and the vector of accumulated currents [/ia] • branch sublink equations described by the connectivity arrays \p ] and \p ], the sublink Thevenin impedance matrix [ZA], the sublink currents [lAsubiink] and the sublink Thevenin voltages vector [VA.5U6HTI*;] A A • branch nonlinear equations described by the connectivity arrays \pAjnoniin] and [PAjnonlin], m a t r i x [zAjnonlin], the vector of subsystem's nonlinear currents [ and the nonlinear voltage vector [ V ^ . n o n ; ] ^ Ajnonlin] i n T h e meaning of the nonlinear equations' matrices and vectors w i l l be clarified i n the example at the end of this section. Here we apply M u l t i l e v e l 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]'(p a A A + z) A ^Asublink ~i" {p A .nonlin) ' ^Ajnonlin ~l~ VA® ' ^ot — P &A "I V _ link (3.15) In a procedure equivalent to.that described i n Section 3.3, the sublink currents are reduced from the subsystem's equations i n (3.14) and the following system of equations is obtained: a 1 - A I 0 A &MTE PA-nonlin MTE b Q-MTEnonlin 0 VAjnonlirflMTEnrnilin "1" ZAjnonlin I 0 o 0 v A VB + qb + Z t la l M sub 1 Ajnonlin a p a TE A (3.16) &A.MTE P A-nonlin A.MTE + V _ n e e A non n B P e-A.MTE + q e l +V t B a where a = A~ p ^Ajnonlin A P _nonlin a = A~ p A(l _ li e = Ae l A A = A A — CL (p Q, + Z ) n A A A A p • & _nonlin A A Aa = a {p a + z )~ p\ • a l l A non A A~ h l A A (^MTEnonlin — ^-Ajnonlin &MTE — a, — Aa A A = a (p a A A A + z )~ • (p\ • e + V _ u k) l A A A A sub n Aa _ n A non n &A.MTE = e — Ae A A (3.17) T h e nonlinear currents equation is extracted from (3.16) i n the following general form: nonlin ' ^MTE-nonlin ~\~ ZAjnonlin) ' ^ Ajnonlin H~ PAjnonlin ' &MTE ' = PAjnonlin ' AJATE + VAjnonlin (3.18) e T h e system of equations is then reduced by moving the unknown nonlinear currents to the right-hand side to obtain the following: "1 0 0 1 0 0 &MTE b P^MTE ZA-MTE — Ajnonlin e 'VA' • + qb t + z VB = (3.19) eB P {eAJATE — eAjnonlin) + q*e + V l B 27 a 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]- T h e link currents are calculated from: ia — a.MTE Z i oi.MTE (3.20) €-ajnonlin) e where ^ajnonlin — P ' ^-Ajnonlin and [z _MTE] and [e _MTE] are as i n (3.10). T h e link currents [i ] are returned to subsystem A and the nonlinear equations are recalculated to produce a new vector [iAjnoniin]- Iterations are repeated until the nonlinear equations satisfy their solution tolerance, then the simulation proceeds to the next time step. T h e 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 i n the following subsection. a a 3.6.1 a Example: A Nonlinear Control Function Implementation with Multilevel M A T E T h e simultaneous solution of a nonlinear control function and the network requires an iterative procedure. T h e Multilevel 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 i n F i g . 3.6 to demonstrate the basic principles of solving a nonlinear control function equation simultaneously with the network solution. M u l t i p l i c a t i o n of two state variables, for example, makes the control scheme in this example nonlinear. T h e system i n question is partitioned into two subsystems, A and B, w i t h subsystem A containing a nonlinear control equation that controls the voltage of node 3. T h e partitioned system equations can be written i n general terms from (3.14) as: A Petri Petri ~ ctrl 0 [ p* 0 where V i VAjnonlin ctr 0 z = V\ • v 2 B 0 ] q l P VA 0 Ictrl h A — V trl c h q —z (3.21) B is a nonlinear control function that corresponds to the vector i n (3.14). Similarly, [p i], [p*td a ctr 28 n d Wtri] correspond, using more general Chapter 3. Multilevel MATE Subsystem A \ v - --j0^ ,(t) / (t) <pv,(t) gA23=2 S A/\AA— gA20=l S gA30=l S« li(t)=l A / Subsystem B / Rlin^lQ ' , gB45=2S J_Uy\A/N | 7 /\/V\Ailink(t) I / 5 gB5o=i s: 1A \ -t T / \ / \ Figure 3.6: A n example system with a nonlinear control function notation, to \p _noniin], \PA.noniin\ A a n d respectively. [ A.noniin\, z For the example i n F i g . 3.6, the hybrid system of nodal and branch equations can be written as: -2 5 -2 0 3 -2 0 0 0 -2 3 -1 [ 0 0 1 0 0 -1 0 0 0 1 v 3 Vi • v ''Ctrl 2 -2 -2 3 [ - 1 0 ] I0] 1 0 0 2 0 1 0 V 4 V _ 5 -1 link l The asterisk (*) sign denotes variables that are recalculated during the iteration process. The nonlinear equation relating the link currents w i t h the nonlinear controller function is extracted as i n (3.18): (Petri • ctrl a where a tri = A c + Zctrl) • ^ctrl + Petri ' ° ' *a = Petri ' A ctrl + Vc e (3.22) • Petri- For the example system, the nonlinear equation is: 0.5238 • ietri - 0.5238 • iu nk = - 0 . 1 9 0 5 + vi • v 2 The remaining linear part of the system is represented by the M A T E matrix i n 29 Chapter 3. Multilevel MATE general terms as: " 1 0 a 0 1 b _ 0 0 pta + qtb + z &A VA . &ctrl = v B (3.23) _ p* {e - e ) + o*e + V _ IQ, A ctrl B a where e i = a i • i tri 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: ctT ctr c 1 00 0 1 0 0 01 0 0 0.1905 0.2857 0.5238 -1.5 ' -1 3.0238 0 1 0 0 1 0 0.5238 + 0.1905 • w ' 0.2857 + 0.2857 • i i 0.1905 + 0.5238 • i i ' 1 " 1 Vl V2 ctr ctr v A - 0 . 8 0 9 5 + 0.5238 • i i I'link ctr Equations (3.22) a n d (3.23) completely describe the example system of (3.21) w i t h the clear separation between nonlinear and linear parts, b o t h represented by their Thevenin equivalents as seen from the system's l i n k i n g nodes 3 and 4. T o recapitulate, the two equations that iterate against each other at each time step of the simulation are i n general terms: \Pctrl ' ctrl + Z rl) • i l a ct ctr a e t z •i — e a (3.24) = p\ rl ' A + Vctrl ~ p\trl ' ' *o a a &a-ctri (3.25) where, for this case, e _ i — p • a i • ictri, ox evaluated: l a ctr ctr 0.5238 • i tri = - 0 . 1 9 0 5 c 3.0238 • i Hnk • v + 0.5238 • i 2 Hnk = - 0 . 8 0 9 5 + 0.5238 • i i ctr A t each iteration step, subsystem A assumes the last known value of the nonlinear controller current i [ and calculates its contribution to e _ i, which modifies the Thevenin equivalent voltage. T h e modified Thevenin equivalent of subsystem A is then passed to the links t o obtain the link current from (3.25). T h e calculated link current is returned to subsystem A. T h e nonlinear controller requests the value of the link current (i ) from the subsystem as well as the voltages v\ and v needed to recalculate its nonlinear equation (3.24). T h e new value of i i is passed back to subsystem A i n order to recalculate its modified Thevenin equivalent voltage. T h e 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 u n t i l the ctr a ctr a 2 ctr 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. T h e 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 llM-0.2677 3 4 5 6 7 8 9 10 11 12 13 14 -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 equations alone. If there were more nonlinear controllers i n the system, each one would perform its iterating procedure i n the same manner, w i t h 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. T h e 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. T h e iterative procedure that we used for solving nonlinearities corresponds to the fixed point iteration method. T h e main reason we chose this method is its simplicity and generality. T h e 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 Multilevel M A T E concept needs to be investigated i n the future. 3.7 Multilevel M A T E Block Diagram F i g . 3.8 shows a block diagram of the Multilevel M A T E implementation. T h e 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 a. l a c a'_arl J y-a c Calculate controller's voltage(s) ;,(••) • • Solve nonlinear controller equation yes • Solve for subsystem's noiLiI voltiiucs Figure 3.7: F l o w chart of iteration process for nonlinear controllers 3.8 Summary T h i s chapter introduces the new concept of Multilevel M A T E , which is suitable for implementation i n the real-time simulator O V N I . T h e advantages of subsystem partitioning are presented i n terms of computational speedup, and the advantages of modelling w i t h branch equations are demonstrated through a switch implementation and a nonlinear control function implementation. 32 Chapter 3. Multilevel MATE ELEMENT IN SUBSUBSYSTEM A, nut -5 -Si S U B - S U B S Y S T E M /4, Knows: Calculates: A"' = A' • a, = A; • />, 1 = A' • K S U B S Y S T E M A + SUBLINKS Knows: <v noniin «A1 PA> <M' ^A^wblink^ ,Aa = e a . ^= = a non/in _ ."AK . 3. a Calculates: .„»„ A2 A 2 P'A O A ^MTE ^A^iionlin f MTE A_aontiR -a eA+V ) AMM Aa = a A , , A A^rt A ' *A_\ubImk ^^A_nonlin ^A ^A„i\onl'm = a (P 'rt_nnn/rn Al_iumlin AfTE ^iiontmA_ a (rt ^4_iwirf«i Z J rt itwrf/n rt_ MTE. iT^rt ^ «onto >: i rt i mwttfn r*AfT£ . fl c yes -Aa, "ME ="-Aa A_WE A ~^ A e =e e NONLINEAR ELEMENTS ' Calculate nonlinear node?~(nn) voltages and sublink branch currents: MTE_lumlin A_ntiulin P MTE> P a Nonlinear v (nn) = ( e _ * m : ( n n ) - ^ _ „ ( n n ) ) -f l ( w i ) • / „ A x non I M^non/i { A_MTE n wnr • equations ^rt^ mMii«i;rt £ atm/iit' A_nonlin) e e S Y S T E M + LINKS Knows: Z, V ,p a , q b , r a MJE UJE C M T , E p {e ~ A_nonlin)' e AmE 1 [ B_MTE ~ B_nonlin)' e e r [ C_MTE e ~ C e _nonlin) Calculates: Z<z_MTE = P' MTE a a.UIE e +<l' MTE = P' { A_MTE e 'a - <-a_MTE b + '' MTE + * r C ~ A_nonlin) e V + ( B_ MTE e B_ nonlin J C Mre C _ nonlin ) e + V a a_MTE c Figure 3.8: Block diagram of Multilevel M A T E implementation 33 Chapter 4 Modelling with Multilevel M A T E 4.1 Introduction The Multilevel M A T E concept offers many advantages i n the modelling of power system networks and their components, as w i l l be shown i n this chapter. Some of the advantages are associated w i t h higher computational efficiency of the models compared w i t h other current modelling techniques. A h example demonstrating this efficiency is the model of a switch described i n the preceding chapter. Multilevel 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 Multilevel M A T E are associated w i t h modelling accuracy. T h i s is especially evident i n cases of modelling of ungrounded or weakly grounded circuits and subsystems. T h e nodal analysis typically used i n the E M T P type of solution may suffer ill-conditioning problems when weak connections or no connections to the common reference point (ground) exist. T h i s is reflected i n the system's admittance m a t r i x becoming singular. A similar problem arises i n the modelling of networks that are grounded only through voltage sources i n the context of the M A T E solution framework i n the O V N I simulator. In this chapter, we present a new approach to achieving good conditioning even under these l i m i t i n g situations. It solves the problem of ill-conditioning for many power system components. Finally, we present the Multilevel M A T E - b a s e d modelling of power system components used i n O V N I , which offers a simple, more general and more intuitive modelling approach. The research contributions reported i n 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, phasedomain induction and synchronous machines, and controllers with functional sublinks. 34 Chapter 4. 4.2 Modelling with Multilevel MATE Modelling of Ideal Voltage Sources W i t h i n the concept of Multilevel M A T E , ideal voltage sources are naturally modelled as functional sublinks. T h i s approach has several advantages, the most important one being generality. Whether an ideal voltage source is connected to the ground or between the nodes i n 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 w i t h the user. A simple example to demonstrate ideal voltage source modelling is shown in F i g . 4.1. v =i v 2 v,=2 v g!2=l S g34=l S -AAAA- -A/VV^ + "g2o=2S • 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 V connected between nodes "2" and " 3 " . T h e system of equations w i t h voltage sources modelled as sublinks can be written as 2 012 -012 -012 012 + 020 0 0. -1 0 0 0 0 0 -1 0 0 0 0 0 034 -034 -034 034 + 040 0 0 0 -1 0 0 0 +1 0 -1 v 2 +1 0 0 0 ii i-2 0 0 0 0 (4.1) -Vi -v 2 where (4.1) corresponds to the M A T E subsystem matrices and vectors described as A PA PA -Z A VA . 1 Asublink h A ^Asublink (4.2) Referring to the procedure of solving the system of h y b r i d equations using the Multilevel M A T E concept described i n 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 w i l l contribute to the Thevenin equivalent voltage vector CA.MTE of the subsystem to obtain: + ZA)' iAsubiink = (PA ~ PA A v X 1 = e .MTE = A- p ; l A 4.2.1 A A • (PAA-^A iAsubiink ='[ + V _ ubiink) = 2.0000 A S 0.7273 [ 1-2727 -0.2727 0.1818 f -0.0909 ] Commentary on Modelling of Ideal Voltage Sources B y examining the hybrid system of equations (4.1) one can observe that i n the traditional E M T P formulation, the system of nodal equations would be smaller. The reason for this is that i n the E M T P ideal voltage sources are internally converted into Norton equivalents. T o be more specific, an ideal voltage source would "borrow" the series impedance from the network, converting it i n 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 w i t h a branch equation (functional sublink). If a source is non-ideal, the internal nodal voltage (e.g., node " 1 " i n F i g . 4.1) would not be calculated, which corresponds to the situation i n the E M T P (there is no increase i n the number of nodal equations over the E M T P approach). If a voltage source is ideal, its internal impedance will be set to zero and the system of equations w i l l have the form as i n (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, i n which case its nodal voltage would be included i n the solution, or non-ideal, in which case its internal node would be non-existent to the surrounding network, resulting i n a smaller system of nodal equations. T h i s approach is more general than that implemented i n the E M T P and is especially useful i n 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 M o d e l l i n g of ungrounded circuits and subsystems using nodal equations results i n singularity of their conductance matrices in O V N I . To overcome this problem, different 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. Subsystems 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 i g . 4.2. The system has four nodes, two of which ("3" and "4") belong to the ungrounded circuit. 1 I,=2-cos(cot) (Q^) ^ gi2=2 S 2 1:10 —VW—— g l 0 g34=l S =lS 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 i or i that are related to each other w i t h the transformer's turns ratio n = i /i . The hybrid system of equations for this particular case can be written i n general terms as: p p 010 + 012 -012 -012 012 0 0 0 0 n 0 s s 0 0 0 0 034 -034 -034 034 -1 + 1 0 n -1 +1 0 h V3 0 0 _0_ 0 (4.3) The system's conductance m a t r i x i n (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. N o d a l voltages v and v 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 v — v± = n • v - Unless the common reference point is introduced artificially, the solution for nodal voltages cannot be obtained. 3 4 3 2 Now when the problem and its solution have been identified, the common reference 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. T h i s approach, however, affects the system's solution i n 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 number of transformer windings and winding connections, the approach shown i n F i g . 4.3 is proposed. Figure 4.3: M o d e l l i n g of an ungrounded transformer circuit w i t h the Multilevel M A T E approach Shunt impedances z = 10, can be inserted from nodes "3" and "4" to the common reference point, i n this case chosen to be the ground potential of 0 V . To compensate for the extra current that the branch voltage v = u — v± has to supply to the shunt impedances, voltage-dependent current sources of v /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 i n ungrounded circuits, are referred to as csn s 3 s 38 Chapter 4. Modelling with Multilevel MATE the compensating shunt impedances. T h e circuit of F i g . 4.3 can also be depicted i n a more compact representation as shown i n F i g . 4.4. g34=l S v © S v /2(J) 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 i n F i g . 4.3 or F i g . 4.4 can be written as 010 + 012 -012 -012 012 0 0 0. 0 n 0 0 0 -n/2 034 + Vcsh -034 -1 n/2 -034 n -1 034 + Vcsh +1 0 " +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 w i t h 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 = (P Vl V2 V3 Vi — A — A V e • A + ZA) a A 1 - (PA • A + V .sublink) = 0.13245 e A 0.67550 ' 0.013245 0.066225 -0.066255 A ' 1 Asublink — a A V where (4.4) corresponds to the M A T E matrices and vectors described i n (4.2) and a = Ai-i„p e = A~ hj IA v A A l A 39 Chapter 4. Modelling with Multilevel MATE Note that w i t h 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" a n d "4" to zero (# = 0) and changing the sublink equation to i = 0 i n (4.4). 34 s The correctness of the solution can be checked against that obtained using the branch voltage equation v -g ~ i = 0, instead of the two equations for nodal voltages ^3 and Vi, i n the secondary transformer circuit. E q u a t i o n (4.3) then becomes: s 3i s " #io + 9u -9\2 -#i2 #12 0 °~ . 0 0 0 " 0 n 1 Tv 1 x v _ 2 -1 v - 1 | (T~J [ i #34 [" J i " 0 0 s s J [ 0 with the solution i = 0.1325 A s [ v Vl 4.3.1.1 v ] =[ 0.6755 T 2 s 0.0132 0.1325 ] V T U n g r o u n d e d Ideal Voltage Sources i n U n g r o u n d e d Circuits In a more general case, an ungrounded circuit that contains ungrounded (branch) ideal voltage sources can be artificially grounded i n the same manner as shown i n the case of the ungrounded transformer circuit i n F i g . 4.4. However, there can be only one artificially introduced reference point i n the ungrounded circuit, meaning that the compensating shunt impedances can be added to only one ungrounded voltage source. T o demonstrate this, a simple network w i t h two ungrounded ideal voltage sources is considered i n F i g . 4.5. _ | i gi2=l S 2 g =2 S 23 3 i l • V\AA-^-AAAA^ I i Sin J~s V,=5V Q ) j V l / 2 0 L.,,, 4 g =3 S ^ 45 g =4 S 6 56 Figure 4.5: Introducing a common reference point to the ungrounded network w i t h two ungrounded ideal voltage sources B y adding no compensating shunt impedances (shown dashed i n F i g . 4.5) to the ungrounded voltage source V i , the system's conductance m a t r i x 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 0 -1 0 0 -012 012 + 023 -023 -023 023 0 0 0 0 0 0 0 0 0 0 0 0 -1 G 0 0 0 0 0 0 -045 045 + 056 -056 -056 056 045 + Vcsh -045 0 0 0 +1 0 0 +1 -1 0 0 +1 0 0 0 . 0 0 " 0 -1 0 0 +1' 0 0 ~Vl ' vz r & 1 02 0 _*i 02 V4, v 5 0 -v h iz x -V3 and can now be solved, w i t h the following solution: [ii [v iz f = x v 2 v 3 [ -2.4 V i v 2.4 ] T A v ] =[2.5 4.9 T 5 e 6.1-2.5-3.3 -3.9 ] T V The current flowing through the circuit i = — i can also be calculated from K i r c h hoff's Voltage L a w to verify the solution of the ungrounded network x VI-V 923 934 3 3 fl45 -2.4 A 556 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 example system it would give the solution i — —2.5704 A. x 4.3.1.2 G r o u n d e d Ideal Voltage Sources in U n g r o u n d e d Circuits Let us now consider the system i n F i g . 4.6, which has two grounded ideal voltage sources. T h e 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 i n the case w i t h 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 w i l l not modify the overall network solution (but is, of course, unnecessary). 41 Chapter 4. Modelling with Multilevel MATE g =l S 1 f V,=5 V ( t ) v, 12 g =2 S 2 23 i 3 ri) 11 n 3 >(£ 3 rt) in v v =iov 3 Figure 4.6: Introducing a common reference point to the ungrounded network w i t h two grounded ideal voltage sources For the system i n F i g . 4.6, the hybrid system of equations can be written as 012 + Vcsh -012 0 -1 0 -1 0 0 0 0 0 -012 012 + 023 -023 0 0 -023 023 + Vcsh 0 -1 0 0 -1 v 0 0 ii h 2 V3 = Vi 0 V -V -V 3 3 w i t h the solution [ i! i ] [v x 4.3.2 = [ -3.3333 T 3 v 2 v ] T 3 = [ 5.0 3.3333 ] A T 8.3 10.0 ] T V 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 w i t h a reference point that is different from the ground. T h e nodal solution for the subsystem voltages 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 i n m i n d , we derive the following example to demonstrate how ungrounded subsystems can be handled in O V N I . The ungrounded wye-connected three-phase resistive load i n F i g . 4.7 represents a good example of an ungrounded subsystem connected to a grounded network. T h e voltage sources on the grounded side of the network are chosen so that the system is perfectly balanced. T h e simple configuration of Subsystem B can be considered as the Thevenin equivalent of a more complex subsystem, as seen from the linking 42 Chapter 4. Modelling with Multilevel MATE Ungrounded Subsystem A Grounded Subsystem B gi4=2S r _iink—l A A A A g =2S & a 58 A A A / - M A / ^ la_link\ g24=2S & = A A A A g34=2S r a 9 V =-1 V b AAA/^—O- link l = c A A A A 1 V =2V - - 0 - n u g =2 S' 6P fb_link l 8 -AAAA A A A / ^ - O ^ link Figure 4.7: Grounded system w i t h an ungrounded subsystem nodes. L e t us now find the Thevenin equivalent of the ungrounded subsystem A . T h e Thevenin equivalent of any ungrounded subsystem as seen from the linking nodes (in this case nodes 1, 2 a n d 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 w i t h respect to the new reference point u / as follows: 3 r e 014 0 -014 0 024 -024 .-014 -024 Vl V2 "-1 0 0 - 1 0 + 0 014 + 024 + 034. 0 0" 0 ^aJink "if lbJink 0 J'cJink_ 0 0 — 0 ' Vzref .-034. (4.5) This equation can be written i n more general terms as: A •v + p •i A = h a - A h_ A ref using the familiar notation for M A T E partitioning, and introducing the new term hA_ f that corresponds to the contribution of the reference voltage to the subsystem's vector of accumulated currents: re h f = Ajre PA-V ,f rl • Vf re B y multiplying the system of equations w i t h A' , we obtain the following expression: 1 v A + a• i = e a A 43 - -1 -1 -1 ref J (4.6) Chapter 4. Modelling with Multilevel MATE F r o m (4.6) we see that the reference voltage increases a l l nodal voltages i n the ungrounded subsystem by the same amount. T h e nodal voltages i n subsystem A are therefore referenced to the voltage of node 3. Therefore, the T h e v e n i n equivalent as seen from the l i n k i n g nodes 1, 2 and 3 is: z = pa t a + qb + z = e = P^A + Q e t a 1 0.5 0 t B + V = A 0.5 1 0 0 0 0 0 0 0 + 2 -1 -1 + 0.5 0 0 0 0.5 0 0 0 0.5 = 2 -1 -1 0 0 0 + 1 0 0 0 1 0 0 0 1 + = 2.5 0.5 0 0.5 2.5 0 0 0 1.5 (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 linking nodes, we w i l l obtain, for the example system, the following hybrid system of equations: " 2.5 0.5 0 +1 1 0.5 2.5 0 +1 0 0 1-5 +1 _ +i +i +i | I" i aMnk _ 1 i ji k b n I" 2 -1 = JcJjnk o J [ v 3 -1 j [ o which can be written, using the more general M A T E notation, as: Pv f re CQ, % Pvref a Vf _ Z re . ref . V . Kef . (4.8) B y applying the M A T E inversion concept to these modified link equations, reference voltage equation can be extracted and solved independently from the link currents. T h e connectivity array p that resembles the p and q connectivity arrays on the subsystem level can be calculated from Vref PVref = P^^PA-Vref (4-9) and Pl = \Pv ] T r e f ref 44 (4-10) Chapter 4. Modelling with Multilevel MATE For completeness, the equation expressions to calculate z Vref v z re f — PA.v A ref PA.v reS are given as (4.11) Av ref (4.12) =PA.v ~ A-V v A lh ref where the terms z v and VA A reference voltages as A Vref z ref Kef and h A ref are defined by a general expression for subsystem Vref +z PAv f .A V v AVref Te =V ref (4.13) AVref T h e reference voltage equation for the ungrounded subsystem can be extracted and solved as follows v -1 a Z 3r r Z •\P V f z r +1 +1 +1] +1] -1 1 PVref ' V e >) re -1 a "2.5 0.5 0" 0.5 2.5 0 ~+l +1 _ 0 0 1.5 _+l_ "2.5 0.5 0 " 0.5 2.5 0 -1 0 0 1.5 -1 -1 -1 (4.14) ' 2 ' -[o] -0.25 V T h e link currents can then be calculated from (4.8), t a k i n g into account the contribution of the reference voltage to the Thevenin equivalent as: 1 taJink i(x IbJink Z a • (e Q Pv f re -0.5. ' ^Ve/) (4-15) -0.5 The links system passes the values of the link currents back to subsystems A and B and, i n 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: 0.5 VA = e ~ A pAv l A ref -0.25 • Vref 0 1.5 v = e - b•i = B B a -0.75 V (4.16) V -0.75 In a more general approach, the reference node does not need to be reduced from the ungrounded subsystem equations but can be duplicated i n order to keep the size 45 Chapter 4. Modelling with Multilevel MATE and structure of the system of nodal equations untouched. T h e effect of duplicating the node appears i n the equations as adding a conductance of one Siemens between the original node 3 and its reference duplicate 3 r e / . However, since v — v f there w i 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) will have the following form instead: s " #u 0 0 —#14 0 0 0 0 ^4 + 1 ~924 — 934 4 92 -#14 -#24 -#34 #14 + #24 + #34. + 3re Vl v 3 + "-10 0 " r• "i iadink 0. - 1 o. ibdink 0 0 - 1 J'cdink_ 0 0 0 — "0" 0 0 0 " 0 " 0 -1 0 Vsref (4.17) 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, p f represents. T h e b o t t o m line is that both approaches produce identical systems of equations (4.8). vre Note that the choice of reference point for the ungrounded subsystems is arbitrary. In fact, i n the particular case examined here, the choice of the neutral point of a threephase 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. A s will be shown i n the next section, a model of a three-phase electric machine w i t h 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 w i t h the Multilevel M A T E concept. Electric machines are the most complex electric power components a n d 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]. B o t h approaches adopt machine modelling i n dqO coordinates, which leads to less demanding computational requirements. T h e dqO modelling is especially convenient for the implementation into transient stability programs, where power system net- 46 Chapter 4. Modelling with Multilevel MATE works are represented w i t h their single-phase equivalents. T h e situation when simulating a dqO machine model w i t h the E M T P - t y p e of programs is less favourable, as the machine model has to be interfaced w i t h the three-phase detailed network representation. 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 limited i n its representation of the internal machine structure (internal faults, saturation effects). W i t h this limitation i n m i n d , phase-domain synchronous and induction machine models were proposed i n [25], [29]. The phase-domain models of three-phase induction and synchronous machines of [25], [29] have been implemented i n O V N I exploiting the Multilevel M A T E concept. T h e detailed description of the derivation of the phase-domain induction machine model is presented and its implementation w i t h M u l t i l e v e l M A T E is explained. Because of the similarities i n modelling of the induction and synchronous machine phase-domain models that w i l l be outlined i n the corresponding section, only a brief description of the implementation of the synchronous machine phase-domain model w i t h the Multilevel M A T E is given here. T h e 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 M a c h i n e Stator and rotor electrical circuits of the induction machine are, i n essence, mutually coupled three-phase inductances, where the stator circuit is stationary and the rotor circuit rotates at an angular velocity uj w i t h respect to the stator, as shown i n F i g . 4.8. F r o m fundamental circuit theory, the voltages across the windings (phases) of an induction machine can be expressed for the stator as: r 5 , ( = a IMh. e = ^ g and for the rotor as: (4.18) b b , D A Bt c°c = ^ + Ri L A A e B = ^ + Ri B B ec = -^ + d 5 l, 4 _ c aa lh + Ri b • 4. /j> j — dt Rdc Figure courtesy of P. Kundur: Power System Stability and Control 47 (4.19) Chapter 4. Modelling with Multilevel MATE Figure 4.8: R o t o r and stator circuits of an induction machine where the flux linkages i n the stator phase a and .the rotor phase A are given by: ipa = L i fpA = L i +Li + L i aa a 4- L i + L i + Li +L i ab b aA a ac c bA b aA A cA c AA A + L i + L + aB B A B i B L ic + L cic (4.20) aC A This can be written i n a m a t r i x form for a l l stator (a, b, c) a n d rotor (A, B, C) phases as: L {0 ) i>A La Lab Lac LaA^r) Lb Lbb Lbc L (6 ) lf>B -^afi(^r) LbB(@r) L (9 ) Lbci^r) L c{6r) ~1pa a L c(9r) a bA LO.B{0T) L c(6r) LbA(6r) LbB{9r) Lbc(6r) L (6 ) LcB{Qr) L c(0r) L LB Lc LBC LBB L Lcc LBC Lc aA a Lbc Lc cA C L A(@r) r a r c AA c cB r r c A AB A ia k i c iA (4.21) iB A where Laa = Lab LA LB Lc L L = a A a AA AB L = L bb cc = Lis + L ms Lc —L = ^ L s = L = L = L cos(6 ) - L - L - L cos(8 + 120°) = L = L — L cos(6> - 120°) = L = Lcc — Li + L —L c = L c = —\L a bc m bB cC sr r bC cA sr r bA cB sr BB A r r B mr mr and • Lis and L ir • L ms are leakage inductances of stator and rotor windings respectively is a magnetizing inductance of the stator w i n d i n g 48 Chapter 4. Modelling with Multilevel MATE • L is a magnetizing inductance of the rotor w i n d i n g mr • L sr is the amplitude of the mutual inductance between stator and rotor windings F r o m (4.21) it can be noted that the mutual inductances between the rotor and stator windings depend on the rotor angular position 6 (t) a n d are, therefore, timedependent. r i, + > e —>— h —»>— + e 2 Figure 4.9: Coupled inductors In the general case of two coupled inductors, as depicted i n F i g . 4.9, the voltage equations (t) = ^ (4.22) ( i ) - *m Cl e2 dt w i t h the flux linkages ^(t) M*) = L (t) • h(t) + L (t) • i (t) = Li (t) • H{t) + L (t) • i {t) n 2 12 2 22 2 can be discretized by applying the trapezoidal implicit integration rule, described i n A p p e n d i x A , to obtain the following system of equations: 2 e (t)_ ~ A t 2 L (t) L (t) n l2 L (t) L (t) l2 2 2 »i(*y + «2(*) ehi{t) e/, (t) (4.23). 2 where the history terms are calculated as e iit) 2 'Ln(t-At) e (t)_ ~ ~ A t L ( t - A t ) h ft2 1 2 L L 1 2 2 2 (t-At)' (t-At) At)' \(t> ( * - At)_ 2 "ei(<- At)" e ( t - At) (4.24) 2 To shorten the notation, i n the equations to follow (t) is omitted and (t — A t ) is replaced w i t h ' to denote the values of variables from the preceding time step. F r o m the system of equations describing coupled stator (4.18) a n d 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 i n a similar manner: a e 0 Rb 0 0 R 0 0 0 0 0 0 a c e A e 0 R 0 0 0 0 0 e B ec 0 0 0 c 0 0 0 0 RA 0 0 0" 0 0 0 0 RB 0 H ic iA is ic Rc La Lb Lc LaA^r) Tb Lbb Lbc LbA^r) Lc Lbc Lc L A(9T) LaA(6r) LbA(Qr) L A{QT) LAA LaBifir) LbB(9r) LcB{9r) LAB Lacier) Lbc(6r) L c{Qr) LAC a a L B{6r) L c{@r) LbB^r) Lbc{@r) L B{QT) L c(9r) a a ^ At a C C C c LBA LCA LBB LBC 'ia' ib ic iA is LBC Lcc }c_ a a C c eha ehb ehc ehA ens euc + W r i t t e n i n compact form we have: R 0 &abc 0 R s &ABC _2_ labc lABC + Xt r &h abc + (4.25) &hABC where the history terms are: 0 Ra 0 0 0 0 0 eha ehb ehc ehA ehB ehc Rb 0 0 0 0 La Lb Lc 0 0 0 R c 0 0 0 0 0 a R B 0 Lc Lbc Lc a a 0 0 0 0 RA Lb Lbb Lbc a 2_ At 0 0 a C aA{Q ) LbA(e'r) LcA(e'r) LCB{e'T) aB(Q ) LbB(e'r) ac(8 ) Lbc{Q ) L c{8 ) r r r r c r r ~i 0 0 0 0 0 I ja h h .i A % .1 B J bc\ % Rc LaA{9 ) L B(O ) LbA(9' ) LbB(e'r) r A L c{Q ) Lbc(0' ) Lcc(0' ) r a r r L A(Q ) C T Lc (e' ) B 'i' i'b r V V r r LAA LBA LCA LAB LBB LBC LAC LBC Lcc i'p . i'c. ' 'a e t % 1 F e e A P 6 . 'c. e which written i n compact form are: ®/i abc £hABC_ R. 0 2 0" ' R ~ At f L S T }ABC. r 50 afac ^abc e ml . ABC. e (4.26) Chapter 4. Modelling with Multilevel MATE T h e above system of equations represents the branch voltage equations for coupled stator a n d rotor windings of an induction machine. T h i s system of equations does not assume any particular connection of the stator a n d rotor windings (a grounded/ungrounded wye or a delta). T h e discretized phase-domain induction machine model can be visualized as depicted i n F i g . 4.10. rotor stator 24, eA eha At AAAc-0-^— \ At At ib 2L AA A + A At 24 -At 2A,c\ At At At A? X 24« y ehB \ ™ " 2L„ 2L , At At ^Ah '2L„ hc . B r At \ <-— 'e IB ——11 . At =^£, At At A < 2iW 24, lA ^r-G-/WV 2L \ ill, 2L„, 24, r h IC At * 24, C Ar Figure 4.10: Discrete phase domain induction machine electrical model 4.4.1.2 M e c h a n i c a l Equations of an Induction M a c h i n e The system of equations (4.25) describing the electrical part of an induction machine contains terms L , which depend on the value of the rotor angular position 6 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 oj as s r r r where to is obtained from the well-known relationship between the machine's accelerating torque (T ) a n d the rotor angular velocity (ui ) r a r T = T -T a e m = j ^ +D Um where • T is the electromagnetic torque e 51 = j ( ^ ^ +D ( ^ U r (4.28) Chapter 4. Modelling with Multilevel MATE • T is the mechanical torque m • J is the polar moment of inertia of the rotor and the connected load • pf is the number of poles of the machine • oj is the angular velocity of the rotor i n mechanical radians per second m • D is the damping coefficient T h e differential equations (4.27),(4.28), discretized w i t h the implicit trapezoidal rule, become difference equations of the following form: 9 (t) - ^u (t) r 2D —— + — ) \PfAt p / 47 f = 6 (t- r At) + ^-uj (t - At) r 4J u (t)-T (t)+T (t) = (\PfAt r e (4.29) r 2D m u (t-At)+T (t-At)-T (t-At) r p e m f (4.30) T h e electromagnetic torque T is calculated from the energy stored i n the machine coupling field Wf, which is i n fact the energy stored i n a l l m u t u a l and self inductances of rotor and stator windings, excluding leakage inductances. T h e energy stored i n the coupling field can be written i n the following form: e Wf = statorjenergy + mutual stator-rotor + rotor .energy 1 2 [i a ib ic] 'Lu 0 0 0 L 0 ls 0 " 0 L i i ]-L^ b B b c Aref) ls + C % lr L Aref) Aref) Aref) A a >(re/) 1 Aref) 'ia' ib + [i •fjref) _ 2 0 0 (re, Aref) Ir ' 0 0 Aref) 0 A % Aref) 0 B l - (re/) (4.31) Aref) C lr l J w i t h rotor quantities referred to the stator side using the effective turns ratio of the stator and rotor windings, and indicated w i t h the superscript ( f\ re The electromagnetic torque T is obtained by solving the following relationship e 'Pf\ 8Wf0d r 52 (4.32) Chapter 4. Modelling with Multilevel MATE which results i n the following expression valid i n the continuous and discrete time domains: T -L i a e ii 1 . JL (re/) . L ib Cj dd ST ' Lms ' [ia ib ic] " l B l Aref) }c r ' ("7^) 'ArefY A Aref) sin 9, sin(0 r- — sui\9 120°) sin(0 + 120°) r r Aref) A Aref) sin(0 + 120°) sin(0 - 120°)' sin6> sin(0 + 120°) sin(0 r l r r B r l A™f) r (4.33) 4.4.1.3 Implementation of the P h a s e - D o m a i n Induction M a c h i n e M o d e 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 w i t h the nonlinear torque equation (4.33). These equations can be combined in a matrix representation of h y b r i d modified nodal equations of the induction machine i n the following general form: A PA -Z / A A 0 P Ajnonlin P Ajnonlin VA 0 1 Asublink ZAjnonlin _ + _ 1 Ajnonlin _ V 0 h PA.v A ref •i 0 0 + a _o_ • Vref = ^Asublink _ ^Ajnonlin _ (4.34) where h _ref A = PA-v f ' Vref re Using these equations, the circuit depicted i n F i g . 4.11 is modelled using Multilevel M A T E to test the phase-domain induction machine model. T h e circuit consists of a three-phase induction motor (subsystem A) connected to a three-phase sinusoidal voltage source (subsystem B) through links (i , if, , i )- 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 i n the case of a doubly fed induction generator, without any modifications to the induction machine model. as s cs T h e system of branch equations (4.25) describing the electrical part of the induction machine contains nonlinear terms which depend on the value of the rotor angular position 9 at time instant t. W i t h a sufficiently small integration step, the prediction r 53 Chapter 4. Modelling with Multilevel MATE Subsystem A Subsystem B bs bs l as v bs v ' as CS Figure 4.11: Electrical circuit configuration for testing the phase-domain induction machine model implemented w i t h Multilevel M A T E of the electrical angular velocity u> can be made using linear extrapolation from the past values u .pred(t) = 2u (t - At) -cu (t- 2At) (4.35) r r r r and used to calculate the value of the rotor angular position according to (4.29). T h i s 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). T h e induction machine implementation w i t h the Multilevel M A T E concept uses the prediction of 6 , but the accuracy of this prediction can be checked and corrected by iterations i f required. r W i t h this i n m i n d , the electrical and mechanical systems of equations are coupled only through the nonlinear torque equation relating the machine's electromagnetic torque T 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 , where the former represents the electrical part, a n d the latter represents the mechanical part of the machine. T h e induction machine coupled branch equations are modelled as sublinks of subsystem A, and the nonlinear torque equation- is modelled as a nonlinear sublink equation. e 2 Following the process of building the system of nodal equations A • v = h 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 A 54 A Chapter 4. Modelling with Multilevel MATE point and solved at the links level. Since the ungrounded subsubsystems i n this case consist of one node each, a l l nodes of the induction machine model are duplicated as reference points v f for their respective subsubsystems. re For the circuit i n 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 V ref ~ 0 VA + lA- iAr ~ VAref = 0 (4.36) a Since v = v f a n d 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 w i t h the nodal equations of the electrical part of the induction machine a are Ai • VAI + PAX • iAXsubiink + V • i = hx- h_ a A A (4.37) ref expanded as 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0" 0 ~V ~ 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 V +1 -1 A Vb C VA + V B vc T v N + -1 0 0 0 0 0 0 0 0. 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 -1 +1 -1 0 0 0 0 0 0 0 0 -1 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ia 0 0 0 0 0 ic +1 0 -1 0 +1 -1 +1 0 0 -1 -1 0 0 0 0 0 0 -1 0 0 ^as ibs %Ar iBr lCr iA + iB ic 0" Varef —Vbref 0 0 0 Vcref — 0 —VAref —VBref 0 0 — VCref —VNref_ Pseudo nodal equations for the mechanical system can be written from (4.29), (4.30), to obtain the following arrays and vectors associated w i t h the mechanical part of the induction machine: A-2 • VA2 + PA2jnonlin " iA2.nonlin — A2 h 55 (4.38) Chapter 4. Modelling with Multilevel MATE and _ At ? 2D 47 p At + Ul p f r f 0„ + ^2w*~r , The mechanical torque T (t) acts as a forcing function, and is therefore included on the right-hand side of the system of equations, which corresponds to sources. m T h e sublink equations of subsystem A represent the coupled branch equations of the machine written i n the following form, which takes into account the stator and rotor winding connections: VAl — A1 Z PAI - ^Al sublink (4.39) = -V,Alsublink In the expanded form: +1 -1 0 0 0 0 0 +1 -1 0 0 0 -1 0 0 0 +1 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 Va 0 0 0 -1 -1 -1_ 1 Vb V c VA VB vc v N "C 'ia ib R 0 L 0 R + At r T s r T i'b ic _4_ iA At 14 e .1 V + ? 'f e ' B A ic i'p e i'c. e'c. P Finally, the equation representing the nonlinear torque equation can be expressed in the following m a t r i x form: (4.40) " 2<A2jnonlin ' 1>A2jnonlin-V,A2.nonlin and PfE ms [1] [Te] = . . -i —Hi— | > b *c| l where N /N r r. s sin0 sin(0 + 120°) sin(0 - 120°)' sin(0 - 120°) sin0 sin(0 + 120°) sin(0 + 120°) s i n ( 0 - 1 2 O ° ) sin0 r r r r r r 'iA N r r r r is the rotor winding to stator winding effective turns ratio. To summarize, the induction machine phase-domain model represented as a six- 56 iB Chapter 4. Modelling with Multilevel MATE terminal device w i l l have the following hybrid form: 0 PAX 0 0 PAX A 0 0 0 2 VAX 0 Px P Aljnonlin 0 0 0 ZA2jnonlin 0 -z x 0 A 0 VA2 PAX.v 0 0 0 ref hA2 ^AXsublink ^AXsublink _ VA2-nonlin _ 1 A2jnonlin 'Jref (4.41) B y applying the M u l t i l e v e l M A T E approach, we obtain the M u l t i l e v e l M A T E system of equations as follows: I 0 0-AX I 0 0 P X AX + Z X 0 0 0 0 p A A A VAX 0 ax O A2-nonlin 0 0 0 0 ZA2jnonlin 0 VA2 j ." ref eA2 1> AXsublink AX sublink i A2jnonlin _ VA2jnonlin _ AXjv 0 0 0 a •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 9 using linear extrapolation (4.35) r • recalculate [ZAX] based on the predicted value of 9 : r R ZAX 2 0 s R 0 L (9 ) sr At [Lj (9 ) + r r r L r • recalculate the contributions of the sublinks to the T h e v e n i n equivalent of subsubsystem Ai\ A a i = ai - a • {PAX^AX + z x) A1 P ii 1 a A A Ae X = e X ~ a X- - (PAX AX + ZAX)' a A A {p\x AX + V Xsublink) 1 (4.43) e A A • calculate the modified Thevenin equivalent of subsystem A as seen from the linking nodes 6 a n d submit it to the links \p a TE]^ t M ^^MMTE], &MTE — dx — Acii &A.MTE = ZAX where (4.44) — AeAx receive the links currents [i ] and reference voltages [v f] from the links and a re 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 6 57 Chapter 4. Modelling with Multilevel MATE recalculate its nodal voltages: VA — e-AMTE — 0,MTEv V f re{ 0>MTEv — O-Alvref ref &&lvref — • update the history terms [/i^] a n — CLMTE^a re = Q>Alvref ~ A {PA A a + A) a Z PA A\vref a d [VAi_ «6hnfc] for the next simulation step s A t every time step of the simulation the links system w i l l perform the following procedure w i t h respect to subsystem A : • receive subsystem's A modified Thevenin equivalent • calculate [Z .MTE] Q and [e _ TE] a from M Za-MTE — P O-MTE + (f^MTE Ea.MTE — P CA-MTE + Z t T + Q^CB-MTE + V a • calculate the reference voltages of subsystem A from Vref = {Pv ' a).MTE ' Pvref ~ v f^ Z ref z Te ' (pl^f ' a.MTE Z ' a.MTE e ~ Kef) • calculate the link currents from ia = ^MTE z a ' { a.MTE e ~ Pv ref ' V f) re • return [v f] and [i ] to subsystem A . re a The M u l t i l e v e l M A T E algorithm w i 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 w i t h results from the literature [30]. T h e machine parameters are given i n A p p e n d i x 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 N m ) at 0.5 s and from nominal back to zero at 2.5 s, w i t h the machine initially operating at synchronous speed. 58 Chapter 4. Modelling with Multilevel MATE Figure 4.13: R o t o r phase current during free acceleration of a 2250 hp induction motor 59 Chapter 4. Modelling with Multilevel 60 MATE Chapter 4. Modelling with Multilevel 61 MATE Chapter 4. Modelling -1000 1 0 1 0.5 with Multilevel MATE 1 1 1 1 1 1.5 2 2.5 1 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 i n load torque of a 2250 hp induction motor 1830 ! ! 1820 1810 V 1780 1760 11 0 I I I 0.5 1 1.5 t(s) I 2 1 I 1 2.5 3 Figure 4.21: R o t o r speed during step changes i n 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 i n O V N I in an analogous manner to the induction machine model, using sublink equations and the Multilevel 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 i n the d and q axes to model the paths for induced rotor currents. T h e 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. T h e modelled machine is a steam turbine synchronous generator w i t h non-salient poles. T h e equations representing the electrical part are derived i n a similar manner as described for the induction machine model i n the preceding section. The mechanical equations of the synchronous machine are equivalent to those of the induction machine. Figure 4.22: R o t o r and stator circuits of a synchronous machine The phase-domain synchronous machine model implemented w i t h Multilevel M A T E has been tested and successfully compared against results from the literature [30]. T h e machine parameters are given i n A p p e n d i x D . 2 . Results of dynamic performance of a 835 M V A steam turbine generator during a step increase i n input torque from zero to fifty percent of the rated value are shown i n 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 ! ! 1 A 0 1 2 3 4 5 Time 6 7 8 9 10 Figure 4.25: Electromagnetic torque during a step increase i n 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: R o t o r speed during a step increase i n input torque of a 835 M V A steam turbine generator 66 Chapter 4. Modelling with Multilevel 601 _ L_ 1 0 0 1 ! 1 1 1 1 2 1 1 3 1 1 1 1 4 5 1(s) 1 1 6 MATE r 1 1 1 7 8 1 9 1 10 Figure 4.27: Rotor angle during a step increase in input torque of a 835 M V A steam turbine generator _I i i -10 0 10 2 i 20 I S (el. deg) 30 I 40 I 50 I 60 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 w i t h M u l t i l e v e l M A T E may give rise to questions about how the efficiency of computation for these models can be evaluated. T h e m a i n disadvantage of a phase-domain machine model is that its impedance m a t r i x 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 w i l l increase the size of the system of equations to be solved, and cause the m a t r i x of coefficients to change at every time step of the simulation.. In the approach proposed i n this thesis, the changing impedance m a t r i x of the machine modelled w i t h 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 every time step. T h e 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 i n the E M T P - t y p e of simulators provides significant advancement over the commonly used dqO model w i t h 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 typical time step used i n the E M T P simulations w i t h a dqO model is i n the order of 50 fj,s or less. Small time step w i t h the dqO model is necessary to avoid instabilities due to predictions used i n interfacing the model w i t h a three-phase network representation. 4.5 Modelling of Controllers A traditional approach to simulating controllers i n 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 = h ri ct (4.45) . where • [Actri] is a m a t r i x of controller coefficients • [x] is a vector of controller variables (inputs + outputs + internal variables) 68 Chapter 4. Modelling • [h tri] c with Multilevel MATE 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 (4-46) s = 4-A t T—^ 1 + v ; For example, discretization of the first-order transfer function of the general form H(s) = K -°±^ ' a + • b v ai5 0 = Xi(s) (4.47) ' v where X\{s) and X (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 difference equation is obtained 2 -B x (t) 0 with 1 h(t) = B {t lXl + A x (t) 0 (4.48) = h{t) 2 - At) - A x {t x 2 (4.49) - At) and A 0 2 = ao + ai— 2 Ai - a 0 - 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 nonlinear 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 w i l l have the following form: A PA -Z P\ Petri — Petri Z _ trl A 0 0 A c P Zctrl B 0 h lA -V A = ictri • pt VA Q —z A — Vctrl VB h ia ~ °- (4.50) B V where • [Petri] is a connectivity array of controllers i n subsystem A , and Petri ' A V - Z .ctrl A 'U ~ Z trl c ' ictri = ~V trl c (4-51) represents the system of controller equations where [p* . ], generally speaking, is not the transpose of the connectivity array \p tri}- A s can be noticed from (4.50), the controller equations introduce asymmetry into the sublink equations. t7 ( c 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. T h e t h i r d type of controller is implemented at the level of a power system component (component control). Its action is not transparent to the subsystem, but it is embedded i n the information passed by the power system component to the subsystem. The implementation of controllers in O V N I using M u l t i l e v e l M A T E is demonstrated i n 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. T h e 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 i n the dq reference frame, including the rotor dynamics. T h e model includes a bus-fed thyristor 70 Chapter 4. Modelling with Multilevel MATE excitation system w i t h an automatic voltage regulator ( A V R ) and power system stabilizer (PSS). In the example, network transients a n d stator transients are neglected as is traditionally done i n power system transient stability analyses. T h e electrical discrete-time circuit of the example system is shown i n F i g . 4.29. SM Mechanical Svstem Excitation system with A VR and PSS e/d_hist Synchronous Machine (SM)fNetwork Interface kd CO kd hist e kq SM Electrical System (Z) Zhjjtist Figure 4.29: Single-machine infinite-bus equivalent network I>: 1 l + sT R F min l + sT x Aa) r \ + sT 1 + sT w Figure 4.30: T h y r i s t o r excitation system w i t h A V R and P S S T h e bus-fed thyristor excitation system w i t h the A V R and the P S S controller are shown i n F i g . 4.30. Nonlinear equations associated w i t h the limiters are described as: -x (t) + x (t) = 0 3 s max — Xs min for x S m i n <x s for x max < for x s s s < x X S < x min 71 S m a x Chapter 4. Modelling with Multilevel MATE E (t) + K fd A (xi(t) - x (t)) = K V Efd — E for s A for ref E Fm&x Efd = E for F m m E < Ef F m i n F m a x < E < E fd F m a x d Efd < Ep min A complete system of equations of.the excitation system discretized using the trapezoidal rule (bilinear transformation) is shown i n (4.52). 0 -1 0 0 0 1 0 0 -KSTAB i w£i) T 0 0 Et(t) Au (t) Efdjt) xi(t) x (t) xs{t) x {t) A 0 0 0 0 0 0 0 KVf E (t-At)-(1 A r 2 0 K (0,0) re (TW£) 0 -1(0,0) {E t -KSTAB 0 0 0 , Fmax 0 0 1 Fm At) At) )x (t At x (t At) RAt r A E in) T A w ( i - At) -K (O,O)" (l-T 2 2 W (l - T,£) x (tt - LAt) - (l - T ±) 2 2 3 0 {xsmax; 3-Smin) S (4.52) where • E is the synchronous machine stator terminal voltage t • Aui is the per unit relative rotor angular velocity r • Efd = ^jrefd is the exciter output voltage Multilevel M A T E system partitioning is next applied to separate the m a i n components of the system: the network and the synchronous machine w i t h the exciter. The system partitioning structure is shown i n F i g . 4.31, where • [VA] represents the network's nodal voltages i n the dq reference frame • [i ] represents the network's branch currents for ideal switches (i A swli isw2 ' ' ' i n F i g . 4.29) i n the dq reference frame • [VB] represents the synchronous machine stator a n d rotor voltages, linear mechanical variables (5, Ato ) and linear exciter variables (xi,x and x i n F i g . 4.30) r 2 72 3 Chapter 4. •A network PA ' A -ZA P Modelling with Multilevel MATE P IA BSM IB q' -z B q IB B -z 4 P' -VA Figure 4.31: M u l t i l e v e l M A T E partitioning of the single machine infinite bus system • [is] represents the synchronous machine stator and rotor currents, nonlinear mechanical variables (T ) and nonlinear limited exciter variables (x , Efj and E i n F i g . 4.30) e s t • [i ] represents the terminal dq link currents and the rotor angle 5 a T h e synchronous machine model in the dq reference frame is better described by branch equations, as w i t h the phase-domain model. T h e machine's branch equations are included as functional sublinks of subsystem B representing the synchronous machine. T h e stator-side nodal equations are derived i n a similar way as for the phasedomain induction machine model (4.36). The rotor-side nodal equations state that the nodal voltage of the field w i n d i n g is equal to the exciter voltage (with appropriate scaling) and that the nodal voltages of the damper windings are zero. T h e nonlinear equation for calculating the terminal voltage E : t 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 w i t h the infinite bus, a l l variables decomposed into dq coordinates. Since the subsystem does not have a direct connection to ground (only through the sublinks and the ideal voltage source representing 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 w i t h the network relate the infinite bus voltage to the rotor angle 5 as eBd - E sin 5 e = E cos 5 Bq B B The transient response of the system variables to a three-phase short circuit applied on C i r c u i t 2 of the transmission line is simulated using the Multilevel M A T E approach. T h e results for the system's currents and voltages, as well as for the synchronous machine's mechanical variables, are depicted i n F i g . 4.32, 4.33. T h e results agree well w i t h the time responses presented i n [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 instability, inaccuracy and numerical oscillations [31]. These problems are well documented in the literature, w i t h 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 described i n this thesis takes advantage of the proposed M u l t i l e v e l M A T E concept and modelling w i t h functional sublinks. A s 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 w i t h i n the framework described i n Section 3.6. W i t h the fixed point type of iteration, an average value of one to two iterations were needed i n test cases described i n this thesis, w i t h 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 initialization from zero initial conditions, Section 3.6.1). The results published i n [34] using Newton iteration show iteration counts of the same order as those obtained i n 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 y~ — — ^ 0.6 0.4 0.2 A 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 P (MW) e 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 w i t h 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 P (MW) e 3000 1000 - - - 500 - : o I : 0 0.5 LI 1 : 1.5 2 2.5 3 3.5 4 4.5 5 UsL (c) 8 (deg) Rotor angle response for the exciter with AVR and PSS 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 w i t h A V R / P S S : (a) terminal voltage, (b) active power, (c) rotor angle 76 Chapter 4. 4.6 Modelling with Multilevel MATE Summary The concept of M u l t i l e v e l M A T E is further explored i n this chapter. New solutions are presented for modelling of ideal voltage sources and ungrounded circuits and subsystems. T h e modelling of electric machines w i t h branch equations by means of sublinks is then described. Finally, modelling of controllers w i t h Multilevel 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 ( D F I G ) wind turbines are increasingly used i n new wind turbine installations a l l over the world. Growing concerns about the impact of a large number of w i n d 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 M u l t i l e v e l M A T E concept, which, combined w i t h hardware solutions, allows for fast simulation of large power system networks, it represents an ideal tool for testing and developing benchmark models of different w i n d turbine installations. The advantages of the detailed modelling and simulation of D F I G wind turbine systems w i t h O V N I include: (1) three-phase representation allowing for accurate simulation of balanced and unbalanced network conditions; (2) phase-domain full-detail machine modelling, including the modelling of torsional shaft oscillations; (3) modelling of the W T G control systems; (4) and the modelling of power electronic voltage source inverters. T h e purpose of this chapter is to demonstrate how an experimental D F I G w i n d turbine system from [36] is modelled i n O V N I and to describe solutions for modelling challenges associated w i t h the full time-domain ( E M T P - t y p e ) simulation of power systems. Finally, the results of the simulations are compared against traditional stability analysis simulation results obtained w i t h the Transient Stability Assessment T o o l ( T S A T ) . The research contributions reported i n this chapter include the following: • Description and implementation of detailed modelling of a doubly fed generator wind turbine system with Multilevel MATE. • Comparison of the OVNI simulation transient stability simulation. results against the experimental 78 induction results and Chapter 5. Modelling 5.2 and Simulation of a DFIG Wind Turbine System in OVNI Wind Turbine Generators Variable-speed turbines w i t h doubly fed induction generators have become a new standard for w i n d turbines of installed capacity above 2 M W . A s the ratings of such wind farms connected to the power system grid become closer to the ratings of traditional 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 w i n d farm connections to the power system network. T h e 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 i n 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 w i n d 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 w i n d turbines are easier to model, as their responses m a i n l y 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 w i t h a conventional induction generator connected directly to the grid. T h i s turbine was originally only stall-regulated, but i n recent years the larger M W - s i z e turbines have had an active p i t c h mechanism as well. 7 • 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 w i n d turbine for larger units. A slip-ringed induction, machine is used as a generator, w i t h a converter connected between the rotor slip rings and the network. The stator of the machine is directly connected to the network. Wind turbine blades have a controllable pitch angle to turn the blades out of the wind or into the wind. 7 79 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI • Multi-pole synchronous generator connected to the power network v i a a converter. T h e synchronous generator is a direct-drive type w i t h low speed and a higher number of poles. T h e wind turbine and the generator rotate at the same mechanical speed v i a 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. T h e 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. T h e fundamental frequency electrical dynamic performance of a doubly fed i n duction generator is completely dominated by the field converter. Conventional aspects of the generator's performance related to the rotor angle, excitation voltage and synchronism are largely irrelevant. T h e electrical behaviour of the generator and converter is that of a current-regulated voltage-source inverter. T h e 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 w i t h two voltage-source inverters and accompanying controls is depicted i n F i g . 5.1. 5.4 DFIG Model Structure To construct a D F I G w i n d turbine system model, the following model structure is normally considered: • D o u b l y fed induction generator model based on a phase-domain induction machine model • Voltage converter model based on an average model of a back-to-back P W M voltage-source inverter w i t h stator and rotor-side converter control • W i n d model that maps the wind speed to the shaft mechanical power for the turbine. E m u l a t i o n of w i n d disturbances such as gusts and ramps by varying input wind speed to the wind-power module • Crowbar protection 80 Chapter 5. Modelling and Simulation - A rotor A of a DFIG Wind Turbine System in OVNI to gri | voltage converter DC-link A DC Ji PWM DC V rotor-side converter controller • \ stator flux angle calculation DC AC ii PWM \r stator-side converter controller i stator voltage angle calculation Abie, Vc D ii q v ,v ,v a b c Figure 5.1: Schematic of a doubly fed induction generator w i n d 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 i n phase coordinates. T h e 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 v i a a voltage converter. The electrical equations of the two machines are identical, as discussed i n Section 4.4.1. W i t h respect to the mechanical equations, the detailed modelling of D F I G w i n d turbines requires a two-mass shaft representation, which we describe i n 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 n o r m a l 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 i n the network), the generator and turbine acceleration can be simulated w i t h sufficient accuracy only if shaft oscillations are included i n 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. T h e inertia of the gear-box is usually not modelled separately but rather included i n 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 ] it + [6] [ K ] d9_ dt = [0] [ T t ] _ [ T a ] (5.1) - OJ where • [J] is a diagonal m a t r i x 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 m a t r i x of damping coefficients • [K] is a m a t r i x of stiffness coefficients (spring constants) • [T ] is a vector representing the turbine mechanical torque t • [T ] is a vector representing the generator electromagnetic torque 9 turbine Aft gear box+generator Kg Jo t h-AAAD, 0 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 T h e above equation for the two mass turbine-generator F i g . 5.2 can be w r i t t e n as: 'Jt 0" 0 J. 9 d 2 ~e ~ dt t 2 + D + D t -D „ 9. 6 d_ -Dtg Dg + D tg t Ktg coupling depicted i n -K tg 0t 0 ~dt tg 0 T (5.2) which when substituted w i t h [co] = [%' ] (5.2) becomes g T t = J t l dt t - Jg^ + D t U J t + + D DgUJg t - g { U ~^ t D (U tg t + K t g { ~ 9 t 0 ~ K (9 Ug) - - tg t g ) (5.3) - 6g) These two equations can be written i n the following form: T duig m - T g T -f t m = Jg-jj- + = J ^r t dt + DgLOg (5.4a) Du (5.4b) t t where T m = DtgiUt ~ Ug) + K g(9t - Qg) t (5.5) E q u a t i o n (5.4a) represents the generator inertia and is equivalent to (4.28). T h e shaft model represents the turbine inertia (5.4b) a n d the coupling between the turbine and the generator (5.5). T h e difference equations describing the generator and turbine's inertia are produced analogously to (4.30) presented i n Section 4.4.1.2. Note that the turbine torque T is calculated from the mechanical power extracted from the wind using the following equation: t wind T = t 5.6 (5.6) Voltage Converter Model and Control T h e converter a n d its controls highly influence the dynamic response of a D F I G . In this section, rotor- a n d stator-side converter modelling w i l l be described, as well as the modelling of their controls. T h e models represent typical equipment and control structures derived from [36]. T h e voltage converter that supplies the rotor of a doubly fed induction generator consists of two voltage-source inverters linked v i a a D C - l i n k capacitor, as shown i n F i g . 5.1. T h e 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. T h e 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 rotorside P W M converter. T h e rotor and stator-side P W M converters are self commutated converters and are set up by six pulse bridges, as depicted i n F i g . 5.3. rotor-side PWM converter DC-link <DCr IDC V B O - i stator-side PWM converter 'DCs C FT-1/ VDC ,V I a \Vbl UL1 Vcl Figure 5.3: R o t o r 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 (v c) can be related to each other as D V3 VglJine — ~~7^Prn,VDC (5.7) 2^/2 w i t h a similar expression for the rotor-side converter voltage-transfer characteristic. The pulse-width modulation factor P is the control variable of the stator-side P W M converter. T h e converter model is completed by the power-conservation equation (assuming a lossless converter) expressed for the stator side as m VDciDC - v 3Re(V ij „ /* / v s i e 1 J i n e ) = 0 (5.8) The rotor side has a similar expression, where Vsuine and Isuine are phasor quantities of stator-side converter A C line voltage and current, and * denotes a complexconjugate value. The switching frequency of P W M converters is usually several hundreds H z , and the average switching losses are proportional to the square of VDC- T h e switching losses can be taken into account by including a resistance between the two D C terminals shown i n F i g . 5.3. The approximate fundamental frequency modelling approach of P W M converters (average model) neglects a l l switching operations occurring w i t h i n the voltagesource inverters and represents the converter as an ideal, lossless DC-to-fundamentalfrequency A C converter complying w i t h expression (5.7). A c c o r d i n g to the results presented i n several publications, the high-frequency ripple due to switching harmonics 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 i n response to network events [37], and i n practice is small and further limited by the inclusion of supply-side inductors [36]. Hence, i n this example the average model of P W M converters is considered sufficient. Because the induction generator and the network are b o t h modelled i n phase coordinates (abc reference frame), the converter equations are modelled i n phase coordinates as well. 5.6.1 Stator-side Converter Model W h e n deriving the equations of the stator-side converter, it is assumed that the converter is connected to the grid v i a a three-phase line w i t h per phase inductance L and resistance R (see F i g . 5.1). T h e voltage balance across the line is V 1>al a Vb R V ibl + L id r d_ dt Val + hi id (5.9) Vbl Vd where v , v , v are the D F I G stator terminal phase voltages, v i, ff,i, w i are the stator-side converter A C terminal phase voltages, and i i, ibi, i i are the stator-side converter A C terminal phase currents. a 0 c a a c c T h e power conservation equation between the D C link and the stator-side converter, w i t h losses neglected, can be written i n phase coordinates as P, Valial + Vblibl + Vdicl = (5.10) V iDCs DC where a l l quantities are instantaneous values and P nv is the three-phase active power of the voltage converter. T h e 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 nonsinusoidal waveforms [38]. CO T h e stator-side converter is vector controlled i n the stator voltage-oriented dq reference frame. T h e angular position of the stator voltage vector V can be calculated 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 ) . T h e vectors of corresponding stator voltages are shown i n F i g . 5.4. T h e variables i n the reference frame are scaled to have the same magnitude as the original phase-domain variables. A scaling factor of 2/3 is used i n the equations. 5 Vsa y 0_ S _ 2 1 cos ~ 3 0 sin cos V a Vb V. r 85 (5.11) Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI fi-axis b-axis \ d-axis \ q-axis ^ \Oi s w 3/2 \ / / s / s / \ \ 3/2 v : Va/ a-axis, a-axis 3/2.v sa 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 reference frame T h e angular position 6 of the stator voltage vector V S s can be calculated i n the most general case as 9 = y^di^tan - 1 S (5.12) ^ where u is the angular frequency of the stator voltage. T h i s equation can be written in terms of the stator phase voltages using the inverse Clarke transformation as s 9 = tanS 1 [ 2 \ 2 i c ] (5.13) In the particular situation considered here, the stator voltages' magnitude and frequency are dictated by the infinite bus. Therefore, the stator voltage' angular frequency is constant, and 9 = u t. S s Transformation of the phase variables i n 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 1 Vd cos 9 cos (9 - f ) — sin #, - s i n ( 0 - f ) 2 S V2 ' 3 cos (9, + f ) -sm(9 + f) S s (5.14) Vb s Vr T h e d-axis of the rotating dq reference frame is aligned along the position of the stator voltage vector V , as shown i n F i g . 5.4. T h e 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 s v = %/2V cos (co t) = a Vb v c s V2V S = V2V S s cos f cod — cos ( co, t+ V2v cos 9 d 2TT\ — S V2v sin 9 q S („ = y/2v cos ^9 d ( 2TT S 2n 9 S —— (5.15) y ^ = V2v cos ^9 + y ) - V2v sin (^9 + y d S q S where the dq stator voltages (vd,v ) are scaled to have the same magnitude as the R M S values of the stator phase voltages (V ). F r o m (5.15), we notice that by aligning the d-axis of the rotating reference frame w i t h the stator voltage vector position, the (/-component of the stator voltage is equal to zero. T h e 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 i n more detail in A p p e n d i x B . q s If the losses of the three-phase line connecting the converter w i t h the grid are neglected (R = 0), the R M S phase values of the stator voltages (V ) and the converter voltages (V i) are equal. B y taking into account the stator-side converter voltage transfer characteristic (5.7), we can write s s V = a vd 2\/2 • (5.16) v DC T h e power-balance equation for the stator-side converter can be written i n dq coordinates from (5.10), w i t h the appropriate scaling of the dq variables, as VDC^DCS = 3v i i d d (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 — (5.18) Pmidl 2V2 F r o m Kirchoff's 1 law for the D C link currents of the voltage converter, we obtain the following expression, w i t h the resistive losses neglected, s t „dv . DC C— - l DCs 3 ^= PmUl — ~ IVCr - (5.19) iDCr where C is the value of a D C - l i n k capacitor as depicted i n F i g . 5.3. Hence the D C - l i n k voltage VDC can be controlled v i a 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 vectorcontrol approach is used w i t h the reference frame oriented along the stator-voltage vector position, enabling independent control of active a n d reactive power flowing between the stator and the stator-side converter. T h e stator-side P W M converter is current regulated, w i t h the direct-axis current component idi used to regulate the D C - l i n k voltage and the quadrature-axis current component i i used to regulate the reactive power. q v Stator-side D C voltage controller (Stage 2) idi DC K,V Stator-side current controller (Stage 1) fz - idi + I z- 1 0-i \ \ z- 1 ) • — • PI PI o iql + Ki z-a \ z-ll t "9 PI hi 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. T h e reference of the direct-axis current component is set by a slower D C voltage-control loop (Stage 2), as depicted i n F i g . 5.5. T h e reference of the quadrature-axis current 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' and v' represent the regulated voltage drops across the stator-side connecting line. (L,R), which is used to calculate the reference values v and v of the stator-side converter. E q u a t i o n (5.9) can be written i n the dq reference frame as d dl q ql v = Ri + L^- - co Li + v dt r-,. di i 91 V = Rl i + jT~ + 0J Ll i + V, d dl s T ql dl (5.20) Q L q q s d Taking into account that v — 0 i n the stator-voltage-oriented dq reference frame, reference voltages of the stator-side converter are calculated from (5.20) as q v* i = -v' + {u Li d d * i = ~' v - v q s q ql +v) d (5.21) (WsLidl) Finally, the reference voltages are converted back to phase domain using the i n verse dq transformation, as follows: — sin 9 sin (9* — 2TT - s i n (9 + ^ COS0, 2?r cos (8 2l cos(0 + 4? Ki Ki S S S 3 dl J (5.22) JJ 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 M o d e l The rotor-side converter is connected directly to the rotor terminals of a D F I G . Therefore, the rotor-side converter A C terminal voltages are the voltages across the rotor windings. T h e power conservation equation between the D C link and the rotor-side converter, w i t h the losses neglected, can be written similarly to (5.10) as Pconv = V i A A + V il3 B + V ic = C V iDCr DC (5.23) where all voltages and currents are instantaneous phase values. The rotor-side converter is vector controlled i n the stator flux-oriented dq reference 89 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI frame. T h e angular position of the stator-flux vector ^ can be calculated from the stationary afl stator-flux components i n an equivalent way as presented for the stator-voltage vector angular position i n Section 5.6.1. However, i n this case the instantaneous values of the stator-flux linkages ip , ip and ip (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 ip _F, ipb.F and ift _F i n F i g . 5.6, depicting relationships between stator-flux linkage vectors i n different reference frames. Note that by aligning the d-axis w i t h the stator-flux vector position, the ^-component of the stator-flux vector is equal to zero. D i g i t a l filtering of instantaneous stator-flux linkages is described i n Appendix C . s a a b c c ft-axis I d-axis a-axis / i i / / / / / / c-axis / / i i i / i i i i ' / 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 T h e angular position of the stator-flux vector 9 f is obtained from the instantaneous values of the filtered stator-flux linkages as s \1pa-F ~ \i>b.F ~ 90 J 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 1 Vdr 2 VA COS (6 iip) - sin (dsup) Vqr s (5.25) VB vc where 9 i = 9 f — 9 is the slip's angular velocity and 9 is the rotor's angular position, defined i n Section 4.4.1. T h e transformation is valid for the rotor voltages and currents as well. s ip s r r Equations relating rotor voltages i n the abc phase domain and the dq rotating reference frame can be written using the inverse dq transformation as v = V2V cos (usiipt + ip) = V2v cos (0 A r dr V B = y/2V cos V c = y/2V r r (^j t sHp COS (bJslipt + <p- + (p + = = ^ ^ V v sKp C r V sin (0 cos ^9 d r d ) - sHp 0 S (0» P li + ) - V2v sin sinI ^9 - y Y ) sMp sHp qr _ • - — 1/1 sin I 8 p + sH 2 7 r y (5.26) where UJ U = co f — co is the relative angular velocity between the rotor and the reference dq axes. T h e magnitude and phase angle of the rotor-voltage vector with respect to the rotating dq reference frame are defined as S p s r Vr -i <p = tan' (5.27) v qr Vdr Similar to its presentation for the stator-side converter model, the rotor-side converter voltage transfer characteristic can be written as V r VDC 2V2 (5.28) where P the modulation depth of the rotor-side converter. Since the rotor-side converter is controlled to the desired value i n the stator-flux-oriented dq reference frame, it is convenient to decompose the modulation depth into dq coordinates, to mr 91 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI obtain the following expressions: Vdr = Pmrd — =V 7 DC 7 (5-29) mrq r v = V f " D C T h e power-balance equation for the rotor-side converter dq variables w i t h appropriate scaling as 3 {v idr Pconv = VDClDCr = B y combining dr (5.29) and (5.30), we obtain the + (5.23) can be written i n VV) (5.30) following expression: 3 ^DCr = (5.31) 2y^2 (^~ ^dr "F -Prnrqiqr) >mr which, i n combination w i t h (5.19), gives a differential equation describing the D C link of the voltage converter i n dq coordinates as C 5.6.4 — — 2V^2^~ ^ >M 1 2 V*2~ ^^ ^^ >mr r (5.32) ^mrqiqr) Rotor-side Converter Control The P W M converter inserted i n 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 mentioned i n the previous section, the induction machine is controlled i n a synchronously rotating dq reference frame w i t h the d-axis aligned along the stator-flux vector position. In this way, decoupled control between the electrical torque and the reactive power is obtained. T h e 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. T h e reference of the quadrature-axis current component is specified by a slower speed controller (Stage 2) as shown i n F i g . 5.7. T h e q component of the rotor current directly influences the torque, so i can be used for torque or active power control. T h e d-axis. component is a reactive current component, so i can be used for reactive power or voltage control. A s s u m i n g that a l l reactive power is supplied by the stator, i* may be set to zero. qr dr dr T h e outputs of the current controllers v' and v' , i n a similar way to the statorside converter controller, are used to calculate the reference values of v* and v of dr qr dr 92 qr Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Rotor-side current controller (Stage 1) — Rotor-side speed controller (Stage 2) • PI + O- PI PI Figure 5.7: Rotor-side converter controller the rotor-side converter. T h e voltage equations of the rotor circuit i n dq reference frame can be written as _ V n • , (r Ll \dl M)J ' L = ( d r T \ . " " " " ^ " {ref)J^ ' dt " J s LI K s L + ( L - - J ^ % + u (Sntif + f^r r L ^ J dt S U P slip \ (ref)"ms <\~r L , L l ^ '•>s ( r e / (5.33) ) where Rr, L , L a n d L ^ are rotor-referred quantities related to per-phase resistances and reactances of an induction machine (Section 4.4.1.1) as r r 0 e s L r — Ra = L\ 3 2 ^mr L T 3 La (5.34) 2 'N \ 2 jiref) r 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 i n the stator-flux-oriented reference frame, the dq components of the stator flux are ip d = I/J = L i and tp = 0. T h e reference voltages of the rotor-side converter are calculated from (5.33) as s CO, V r = -dr V d qr = qr V V s 0 ms sq 0q r -Aref) + "si + (5.36) L Idi jjref) T h e reference voltages are converted back to phase domain using the inverse dq transformation w i t h appropriate scaling as follows: cos 9, slip COS l U e up s S 27T dr 2TT J 3 shp cos {9sup + yh. 5.7 — sin 9 sin (9 sin 9slip *3 2 3JJ (5.37) qr. Wind Turbine Model The kinetic energy of a mass of air m of velocity v i equation: is given by a well-known winc • . ^ = y ; lind (5-38) v The power of the moving air mass is the derivative of the kinetic energy w i t h respect to time, and can be written as dE k •M) — _ 1 dm —g 2 _ 1 / 2 wind ~~ 2 V W i n d 53 9 ) {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 . O n l y a fraction of the w i n d kinetic energy can be converted into rotational power at the w i n d turbine shaft. T h i s fraction of power (P i d) 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 , representing aerodynamic w n p 94 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI efficiency v = p (5.41) For a specific turbine design, C is normally calculated as a function of the pitchangle P and the tip-speed ratio A. T h e tip-speed ratio is given as p A= ^ (5.42) where R is the radius of the turbine blades and co is the turbine speed. C (X, P) is normally defined i n the form of a two-dimensional lookup characteristic for different values of f3 and A. Alternatively, analytical approaches for approximating the C characteristic can be used, such as the following expression [39]: t P p C = (0.44 - 0.01670) sin * £ _ ~ ^ X p - 0.00184(A - 2)p (5.43) The mechanical power extracted from the wind (with A = irR ) is calculated from 2 P nd wi =\-p-n-R - C (X, 2 P P) • vl ind . (5.44) w i t h the mechanical torque driving the induction generator calculated from (5.6). 5.8 Maximum Power Tracking The family of P i d — w curves can be derived from (5.44) for various values of wind velocity, w i t h u 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 i n F i g . 5.8.. Curve P defines the m a x i m u m energy captured from the wind, and is characterized by the factor K i n the following expression: W n r r opt opt Popt = K L0 Z opt r (5.45) The objective of the tracking control is to keep the turbine speed fixed on the P curve as the w i n d 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. opt One method of m a x i m u m power tracking is realized by observing the mechanical torque on the generator shaft, and driving the D F I G to the o p t i m u m 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 ( F i g . 5.7) to T (5.46) K,opt T h i s method of m a x i m u m power tracking is referred to as a speed-mode control. For extracted w i n d power higher than rated, the D F I G is driven to a stall point for the particular w i n d velocity, where C0„ In F i g . 5.8 the rated power of P = 7 . 5 k W defines the stall point for wind velocities higher than the rated (10 m / s ) . m Q X 14000i 12000H 10000H g" 8000H 6000- 4000- 2000 500 1000 1500 n (rpm) 2000 2500 3000 Figure 5.8: W i n d turbine characteristics w i t h m a x i m u m energy capture curve Another method of m a x i m u m 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* = K io — Bco 2 opt r (5.48) by regulating the q-axis rotor current i* to r (re/) , 2 L qr r 3p L J 2 f 96 (ref) ms (5.49) 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 i n the rotor-side converter controller ( F i g . 5.7). 5.9 Implementation of a D F I G Wind Turbine System in O V N I The doubly fed induction generator wind turbine system model consists of the phasedomain induction generator discrete model described i n Section 4.4.1.3, the discrete phase-domain model of the voltage converter and the voltage-converter controllers operating i n 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. T h e required values of the converter's terminal voltages are calculated as follows: p* mdl * ai = -^VDC v • rt <7l m COS 9 f~ DC sin 6 V S S bi = -^VDC COS \d - — \ ^-v sm v cos I Os + y I 2 v 3 = — T~ DC V cl [ * v A * v B * v p* mdr = —^-v DC cos e up s mdr fn P mdr where P , P i, P by (5.16), (5.28). mdr 2 2T{\ (n C S V m c S • n m T s l i p mqr • P DC M s mqr D ^-v sm I 9 np + - I and P V (5.50) d D DC mq mqrV = —g— v c cos I 6 slip P 19, - — DC p* n = —^v cos c mdi p* n mqr . P ^-v sm DC f n 2 I0 slip /„ 7T - — 2TT I 9 p+ — sH are control variables of the voltage converter defined 97 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI The D C link equation i n the abc reference frame has the following form: c^dt =i al (% V 2 cos e - % 2 s sin e s M^-('.-S)-¥-(*-y)) +id (^ c o s •p! I e + ^ ) - - ^ s i n (e + ^ s s p V V , 2 V ^ Sin -ZB ^ COS ^ -<c (^f cos (e - 1 slip ~ + *f) - ^ sin ^ 3 / 7 (5.51) -f + ^ and is discretized using the trapezoidal integration rule. These equations are combined w i t h 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 w i n d turbine system from the literature using the Multilevel M A T E concept, and implemented it i n O V N I . Reference [36] describes a doubly fed induction generator application to variablespeed wind-energy generation using back-to-back P W M converters i n an experimental setup. T h e paper provides details on the control setup of the stator- and rotor-side converters and provides the data necessary for modelling of the w i n d turbine system in the phase domain. M o d e l l i n g details have been presented i n 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 w i n d turbine to a strong network modelled as a n infinite bus. T h e single-line diagram of the test network is depicted i n F i g . 5.9. For the purpose of comparison, the w i n d turbine generator defined i n A p p e n d i x D is connected to a 250 V bus to replicate the situation i n [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 i n the middle of C i r c u i 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 0.25/10 k V CCT1 /7W CCT2 e =1.6% e =1.6% R x 10 km line i Figure 5.9: D F I G w i n d turbine test case the D F I G w i n d turbine system to a change i n w i n d velocity following the m a x i m u m speed tracking system described i n 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 w i t h the Transient Stability Assessment T o o l 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 i n steady state w i t h a w i n d velocity of 9 m/s, corresponding to the o p t i m a l shaft speed of 1350 rpm and mechanical torque of 38 Nm. T h e m a x i m u m power captured at this w i n d speed is 5400 W. T h e m a x i m u m power-tracking system works i n current-mode control. A t t = 2 s, the w i n d velocity decreases instantaneously to 5 m / s , which corresponds to the o p t i m a l shaft speed of 750 rpm and mechanical torque of 12 Nm. T h e m a x i m u m 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, w i t h the deceleration torque being the difference between the turbine's mechanical torque a n d the torque given by the o p t i m u m curve (5.48). A similar situation w i t h opposite effects would occur i n the case of an increase i n wind velocity. T h i s scenario corresponds to the scenario described i n [36] for the experimental setup. B y simulating the response of the D F I G model built w i t h 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 i n 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. T h e D F I G is initially running above the synchronous speed of 1000 rpm (rated frequency of the system is 50 Hz) for v i d — 9 m/s, and i n the new steady state it runs below the synchronous speed for v i = 5 m/s. T h e electromagnetic torque corresponding to the mechanical torque on the shaft after compensating for the friction losses decreases from the initial 30 Nm to around 7.5 Nm, corresponding to the new steady state w i t h a w i n d velocity of 5 m/s. A slight decrease of the stator flux reflects the decrease i n voltage at the stator terminals. W n W 100 nd Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 4000 -2900 3000 -2950 < > <22000 O-3000 1000 -3050 -3100 -2900 4000 3000 -^ooo 1000 Figure 5.11: Transient response of the D F I G wind turbine to a step decrease i n 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 w i n d 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 i a the converter. For the new steady state w i t h w i n d 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 b o t h v i a 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 w i n d turbine to a step decrease i n 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 i n w i n d velocity. T h e statorside converter active power depicted i n F i g . 5.12(a) is equal to the rotor-side converter active power i n F i g . 5.12(c) (with the converter losses neglected). T h e 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. T h e reactive power of the stator-side converter is regulated to zero by the stator-side converter control (i i = 0). q 102 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) 0 20 (b) 40 60 (c) t 0 20 40 (d) (s) t (s) Figure 5.13: Transient response of the D F I G wind turbine to a step decrease i n 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 i n the D F I G power output causes a decrease i n 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 i g . 5.13(c), the m i n i m u m value corresponds to the time instant when the rotor speed is equal to the synchronous speed. 103 60 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI 0 20 40 60 t (S) 0 20 40 t (s) Figure 5.14: Transient response of the D F I G w i n d turbine to a step decrease i n 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 i n the phase domain. F i g ures 5.14(c) and F i g . 5.14(d) show the smooth operation of the D F I G through synchronous speed. B y examining F i g . 5.14(f) we note the change of "direction" of the stator-side converter current, indicating the change i n converter power from generation to consumption. 104 60 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (b) (a) 200 c- < ° CO CO CO 100 C D « V 1-2 dUsS O > o cr cr TD o qs s -100 55 CO 20 300 0) 200 V eg "o > cr 100 CO 40 60 ———L__ -4 -6 20 (c) (d) 40 60 40 60 40 60 < V •£ C O »CD o cr •o T 3 o O CO CO %1 CO CO qs o V 1— I -100 20 40 60 (e) 50 -1 d1 -2 20 (f) _10 < CO cu cu 5 £ C D o > cr CJ O O rr -50 J / v V o 0 o. ^ "D i_ rr 40 20 60 20 t(s) t(s) Figure 5.15: Transient response of the D F I G wind turbine to a step decrease i n 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 i n the stator voltage-oriented reference frame, and rotor quantities are shown i n the stator flux-oriented reference frame. T h e control of the q-axis stator-side converter current i \ determines the displacement factor at the D F I G terminals, and it is kept at zero i n 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. T h e control of the d-axis stator-side converter current i^i determines the active power exchanged w i t h the network for the purpose of maintaining the converter D C - l i n k voltage at a q 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 reference frame, w i t h the d-axis current component controlling the machine's reactive power, and the q-axis current component being proportional to the machine's electrical torque. A s s u m i n g that a l 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 ^experiment Figure 5.16: Comparison of a transient response of the D F I G w i n d turbine to a step decrease i n w i n d 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 i n [36] for a decrease of w i n d velocity are shown i n 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 w i n d i n g / s l o t 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 i n steady state w i t h a w i n d velocity of 8 m/s, generating 2700 W of active power that is transmitted to a "strong" power system modelled as an infinite bus. T h e 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 i n the middle of circuit C C T 2 , as shown i n F i g . 5.9. T h e fault is removed after 0 . 1 s without t r i p p i n g the circuit. T h e response of the w i n d turbine system to the fault is depicted i n F i g . 5.17 - F i g . 5.24. Large disturbances cause large initial fault currents b o t h i n the stator and rotor. Moreover the D C link w i l l experience overvoltages due to reduced voltage at the stator-side converter terminals. H i g h currents flowing through the converter may activate the crowbar protection i n order to preserve the converter which may then lead to the w i n d turbine being disconnected from the network. It has been recognized i n the literature that stability models of D F I G wind turbines may not adequately estimate the transients i n 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 w i t h the Transient Security Assessment T o o l ( T S A T ) [41], which is a good representative of commonly used stability simulation tools. The exact same case as tested w i t h O V N I could not be replicated due to the design of T S A T [41], which allows accurate computations i n the M V A range rather than i n 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 a n d consume 3 M V A r of reactive power. Default controls of a D F I G i n T S A T were modified to reflect similar conditions as i n the O V N I model. T h e D F I G model used i n T S A T is a user-defined model of type 3 that represents a rotor voltage-controlled D F I G w i t h flux transients. T h e shape and general trends of a l l the responses computed w i t h T S A T correspond to those calculated w i t h O V N I . The most significant difference i n the results is found i n the oscillatory behavior caused by the transient i n the converter's D C link. A s is commonly done i n transient stability modelling, the stator-side converter and the D C link are neglected and only the rotorside converter control is modelled. B y examining the results shown i n Figs. 5.25 and 5.26, one can conclude that this is probably an oversimplification imposed on the model. T h i s type of behavior, which is not well represented w i t h 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 i n this work, both implemented i n O V N I . T h i s 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 w i t h 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 t y p i c a l stability simulation tool. A similar comparison between any other E M T P - and stability-type tool w i 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 successfully replicated i n a few commercially available E M T P - t y p e simulators. However, . these simulators do not implement partitioning techniques. O u r O V N I simulator implements the concept of M A T E (or, i n the more general case, Multilevel M A T E described i n Chapter 3 ) , which allows parallel processing a n d efficient computation at the subsystem level. M u l t i l e v e l M A T E can be applied to simulations of very large power system networks. Stability tools use approximate modelling to cope w i t h the size of the network, while O V N I uses multilevel system partitioning a n d efficient modelling solutions to achieve this objective. Therefore, for the same objective of coping w i t h system size, Multilevel M A T E allows more detailed modelling than do the stability tools. T h e solutions proposed i n 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 t o o l we have replicated and simulated a n experimental setup from the literature. O u r comparison of responses to a step change i n w i n d speed has shown good agreement w i t h the experimental results. In addition, we have the same system simulated for a three-phase fault near the D F I G terminals. W e 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 w i n d 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 (d) -1000 -1500 _-2000 g-2500 < > -2500 a- 2000 O -3000 -3500 -4000 0.2 0.4 t(s) 0.6 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 0.8 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 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 111 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) 0 0.2 . 0.4 t (s) (b) 0.6 0.8 0 0.2 0.4 t (s) 0.6 Figure 5.20: Transient response of the D F I G w i n d 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 0.8 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI Figure 5.21: Transient response of the D F I G w i n d 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 Iw 600 o550 o I ll A A A / ^ ~ ~ 500 450 400 0.1 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 w i n d turbine D C - l i n k voltage to a threephase short circuit 114 Chapter 5. Modelling and Simulation of a DFIG Wind Turbine System in OVNI (a) 200 r (b) 0 1 > co >°" •S > 100 < ~ ds *CO . _°-4 V qs 0 0.2 0.4 > 200 0.8 0.6 0.8 AAA,- • V d1 I V . ' AA. \s CT100 > 0 -100 — q1 0.2 C 50 • 1 IP 1 CT -50 -100 0.4 (e) 1 0 > > > 0.6 r (C) 300 > I A A / W - — CO CT 0 •100 qs 2 < R IAAA/NIAA/V CT -P 0 qr C 5 • 'dr 0.2 0.4 t(s) 0.6 0.8 0.2 0.4 t(s) 0.6 Figure 5.23: Transient response of the D F I G w i n d 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 0.8 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 w i n d 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) (b) Generator terminal voltage (pu) Generator speed (Hz) 1.000 • 53.990 • 44.800 44.960 4^.120 45.280 5.440 45.600 44.800 44.9 45.120 Time (sec) (C) Generator active power (MW) 45'.280 45".440 45.600 Time (sec) Generator reactive power (MVAR) 2.800 (d) -2.000 V : : \: V ^ 44.800 44.960 : — 45.120 45.280 45.440 45.600 44.800 i 44.960 T 45.120 1 45.280 1 45.440 Time (sec) Time (sec) Figure 5.25: Transient response of the D F I G wind turbine to a three-phase short circuit obtained w i t h 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 45.600 Chapter 5. Modelling and Simulation, of a DFIG Wind Turbine System in OVNI (a) Output of WTGUDM block Line current magnitude (pu) 0.300 • (b) 0.042 / 44.800 44.960 45.120 45.440 TUt 45.600 Time (sec) Time (sec) Figure 5.26: Transient response of the D F I G w i n d turbine to a three-phase short circuit obtained w i t h the stability tool T S A T : (a) D F I G rotor dq voltages and currents (from top to b o t t o m - Vd , idr, qr and i ), (b) D F I G stator current v r qr 118 Chapter 6 Conclusion T h i s chapter summarizes the main contributions of the work presented i n the preceding chapters and provides suggestions for the continuation of this research. 6.1 Summary of Contributions T h i s 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 partitioning and functional modelling for the simulation of power system networks. T h e proposed concept by its nature allows a more efficient solution of the subsystems' equations, and provides a means for general modelling of power system components w i t h nodal and branch equations. T h e advantages of subsystem partitioning w i t h the M u l t i l e v e l M A T E concept introduced i n this work are presented i n terms of the computational speedup over the existing single-level M A T E . T h e greatest improvements are realized when simulating subsystems containing components of a changing nature. The advantages of modelling w i t h branch equations are demonstrated through the implementations of an ideal switch and a nonlinear control function. W i t h the Multilevel M A T E concept, the subsystem's admittance m a t r i x 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 w i t h iterations 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. • M o d e l l i n g of Power System Components w i t h M u l t i l e v e l 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. A l s o , the solution for calculating ungrounded subsystem reference voltages is given. We describe branch (sublink) implementation of phase-domain induction and synchronous machine models w i t h Multilevel M A T E . T h e advantage of branch implementation is evident from the fact that the changing nature of the machine's coupled branch equations does not affect the subsystem's admittance matrix. A l s o , the transformation of coupled branch equations into nodal equations is no longer required, which greatly simplifies the modelling approach. M u l t i l e v e l 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 M o d e l l i n g with M u l t i l e v e l M A T E t h r o u g h an E x a m p l e of a D o u b l y - F e d Induction Generator W i n d T u r b i n e System. We test our proposed approach to modelling and simulation on an example of a doubly fed induction generator w i n d turbine system. We describe modelling of the wind turbine system components and control i n detail. We have replicated an experimental setup from the literature [36] i n our simulation tool and tested it for two types of disturbances: decrease i n w i n d velocity and a three-phase fault in the connecting double-circuit transmission line. T h e results were successfully compared against the results i n the literature and against a traditional stability simulation tool. T h e comparison shows the advantages of using more detailed modelling, especially when control and protection devices play a major role i n the system's response. 6.2 Recommendation for Future Research T h i s thesis addresses a broad range of issues associated w i t h the modelling and simulation of power system networks. M a n y aspects of this work can be further investigated and developed, including the following: 1. Diakoptics versus Sparsity? T h i s 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 i n this thesis, but we noted that the efficiency of the multilevel partitioning approach depends on the sparse nature of the network. T h e actual comparison of the two approaches and the possibility of employing sparsity w i t h i n the M A T E concept is still to be investigated. 2. O p t i m a l network partitioning algorithm. O p t i m a l partitioning that takes into account the multilevel approach needs to be developed. T h e algorithm has to take into account the particular system's topology and characteristics. 120 Chapter 6. Conclusion 3. F i x e d point versus Newton-type iteration method? T h e iterative procedure for nonlinearities implemented i n 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. T h e possibility of implementing Newton's method for the iteration of nonlinear equations w i t h i n the M u l t i l e v e l M A T E concept needs to be further investigated. 4. E x p a n s i o n of the D F I G wind turbine system model. T h e w i n d turbine system implemented and tested i n O V N I does not include models of the protection 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 w i n d turbine farm into the power system network. T h e wind turbine model described i n this thesis could be used for stability analysis of large wind farms interconnected w i t h the traditional power system network. It could also be used for benchmarking the approximate wind-farm models used i n transient stability tools. 6. Implementation of latency and eigenvalue analysis into O V N I . M e t h ods have already been developed i n our Power System G r o u p to take advantage of the slow and fast nature of the subsystems as reflected i n the size of the time step used for their integration, eventually resulting i n even more efficient simulation speeds. Eigenvalue analysis provides a valuable insight for stability considerations of power system networks. These approaches need to be implemented into the current version of O V N I . 7. Development and implementation of power electronic models into O V N I . T h i s 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 i n O V N I that take advantage of the multilevel partitioning approach. T h e 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. C o u p l i n g 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 combination of software and hardware solutions. T h e software solution for O V N I ' s implementation w i t h Multilevel M A T E needs to be coupled w i t h 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 converter model for real time transients simulation. IEEE Transactions on Power Systems, 14(1):166-171, 1999. [4] H . W . D o m m e l . D i g i t a l computer solution of electromagnetic transients i n 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 . D o m m e l and W . S. Meyer. C o m p u t a t i o n 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 . D o m m e l , 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, C h i n a , 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, T h e University , of B r i t i s h C o l u m b i a , 2000. [9] J . A . H o l l m a n and J . R . M a r t i . Real time network simulation w i t h PC-cluster. IEEE Transactions on Power Systems, 18(2):563-569, 2003. [10] T . De R y b e l , J . A . H o l l m a n , 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 i n korperlichen Leitern m i t Anwendung auf die thierisch-elektrischen Versuche [Some laws concerning the distribution of electrical currents i n conductors w i t h applications to experiments on animal electricity]. Annalen der Physik und Chemie, 89(6):211-233, 1853. 13] L . Thevenin. Extension de l a l o i d ' O h m aux circuits electromoteurs complexes [Extension of O h m ' 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 theoremof dynamic electricity]. C. R. des Seances de lAcademie des Sciences, 97:159161, 1883. [15] I E E E recommended practice for industrial and commercial power systems analysis. IEEE Std 399-1990, 1990. [16] J . B . W a r d and H . W . Hale. Digital computer solution of power-flow problems. Transactions of AIEE, PAS-75:398-404, June 1956. [17] W . F . T i n n e y 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. C i t a t i o n classic - direct solutions of sparse network equations by o p t i m a l l y ordered triangular factorization. Current Contents - Engineering Technology and Applied Sciences, 18:12-12, 1979. [19] H . M . M a r k o w i t z . T h e 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) :500535, 1977. [21] I. Hajj. Sparsity considerations i n network solution by tearing. IEEE Transactions on Circuits and Systems, 27(5):357-366, 1980. [22] Chung-Wen H o , A . R u e h l i , and P. Brennan. T h e 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 . M u l t i - a r e a Thevenin Equivalent ( M A T E ) . A circuit concept. Research note to H . W . D o m m e l and L . Linares, August 1993. 123 Bibliography [24] P . Zhang, J . R . M a r t i , and H . W . D o m m e l . Network partitioning 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):222229, 1997. [26] F . A . M o r e i r a 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. K u n d u r . Power Sysem Stability and Control. T h e E P R I Power System E n g i neering Series. M c G r a w - H i l l , Inc., 1994. [28] H . W . D o m m e l . EMTP Theory Book. M i c r o t r a n Power System Analysis Corporation, 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 . D o m m e l , 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 modelling of control transfer functions. In Conference Proceedings of the 14th Power Systems Computation Conference, PSCC02, Seville, Spain, June 2002. [33] L . Linares, M . A r m s t r o n g , and J . R . M a r t i . Software implementation of controller representation i n the O V N I simulator. In Conference Proceedings of the 15th Power Systems Computation Conference, PSCC05, Liege, B e l g i u m , August 2005. [34] J . Mahseredjian, L . Dube, M i n g Zou, S. Dennetiere, and G . Joos. Simultaneous solution of control system equations i n E M T P . IEEE Transactions on Power Systems, 21(1):117-124, 2006. [35] L . Dube and H . W . D o m m e l . Simulation of control systems i n an electromagnetic transients program w i t h T A C S . In Proceedings of IEEE PICA Conference, pages 266-271, 1977. 124 Bibliography [36] R . Pena, J . C . Clare, a n d G . M . Asher. D o u b l y 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, M a y 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 w i n d turbine w i t h 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] S e u l - K i K i m , Eung-Sang K i m , Jae-Young Y o o n , and H o - Y o n g K i m . P S C A D / E M T D C based dynamic modeling and analysis of a variable speed w i n d 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 B o o k Company Inc., 1st edition, 1955. [41] Powertech Labs Inc. TSAT - Transient Security Assessment Tool, Version 5.1, 2005. [42] V . A k h m a t o v . 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. T h e E P R I Power System E n g i neering Series. M c G r a w - H i l l , Inc., 1994. [44] T . V . Cutsem a n d C . Vournas. Voltage Stability of Electric Power Systems. Power Electronics a n d 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 E l e c t r i c a l a n d Electronics E n g i neers, 1993. [46] S. B . L i p p m a n . C++ Primer. A d d i s o n Wesley Longman, Inc., 2nd edition, 1996. [47] J . A r r i l l a g a and N . R . Watson. Computer Modelling of Electrical Power Systems. John W i l e y a n d Sons, L t d . , 2nd edition, 2001. [48] E . Clarke. Circuit Analysis of A-C Power Systems. Wiley, 1950. [49] M . L u k i c . Feasibility of security assessment i n power systems using full timedomain solutions i n the E M T P . Master's thesis, T h e University of B r i t i s h C o l u m b i a , 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 Generators 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 using back-to-back P W M converters supplying an isolated load from a variable speed wind turbine. IEE Proceedings-Electric Power Applications, 143(5):380387, 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. Transactions on Circuits and Systems, 23(12):694-705, 1976. IEEE [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 ) : 8 4 - 9 2 , 2002. [64] Anonymous. Modified nodal analysis and the incorporation of multiphase, coupled networks. Engineering Science and Education Journal, l l ( 6 ) : 2 5 4 - 2 5 6 , 2002. [65] R . J . Koessler, S. P i l l u t l a , L . H . Trinh, and D . L . Dickmander. Integration of large w i n d farms into utility 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 . W o n g . Integration of large w i n d farms into utility grids (part 2 - performance issues). In IEEE Power Engineering Society General Meeting, 2003. [67] H . W . D o m m e l . 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 . D o m m e l . 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 networks was first introduced by Dommel in [4]. To demonstrate the use of the trapezoidal rule, the basic relationship of voltage and current for an inductance is considered: Mt) = (A.1) Rewriting the equation and integrating both sides, we obtain J t-At v {t)dt = L j di (t) L ' L (A.2) tk-At where t 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 v (t) between the points v (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 k L tk f / L v (t ) + v (t - At) . v (t)dt = area — —— • At L L k L k (A.3) tk-At A similar approach is applied for the capacitance to obtain its discrete-time equivalent. 129 Appendix A. Trapezoidal Integration Rule Appendix B Reference Frame Transformations The introduction of reference frame theory i n 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, i (t) + i (t) v (t) + a a + i (t) b = c v {t) b + v {t) c 0 = 0 (B.l) where the variables denote stator phase currents, voltages and flux linkages, respectively. A s a result of this redundancy, it is possible to transform the phase-variable representation into an equivalent two-phase representation. T h e transformation from three-phase to two-phase quantities proposed by Clarke [48], is written i n matrix form as _ i,p(t) _ 2 ' 1 cos(^ '2j ~ 3 0 3 COS (f ) sin ( f ) ia{t) h(t) (B.2) i (t) c T h e transformation is equally valid for the voltages and flux linkages. T h e stator current space vector is defined as the complex quantity h(t) = i (t) sa 131 +ji fi(t) a (B.3) Appendix B. Reference Frame Transformations Equation (B.2) can be written i n a more compact form by introducing a vector rotation factor, a — e T", as follows: J I (i) = H [ (t) + ai {t) + a\(t)] s ia (B.4) b The choice of the constant i n the transformations (B.2), (B.4) is somewhat arbitrary. T h e value of 2 / 3 has the advantage that magnitudes are preserved across the transformation. A sinusoidal phase current with a peak magnitude of I w i l l produce a current space vector w i t h a peak magnitude of I . 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 m m ' ia(t) ' ib(t) Q O ic(t) B.2 2 sin sin(f COS cos isp(t) (B.5) dq Transformation The stator current, voltage a n d flux-linkage space vectors are complex quantities defined i n a reference frame whose real axis is fixed to the magnetic axis of stator winding a. T h e corresponding quantities defined for the rotor circuit of a three-phase A C machine are similarly defined i n 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. W h i l e i n 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. T h i s results i n stationary current, voltage and flux-linkage space vectors. The current space vector i n the dq rotating reference frame can be written as h = isd + ji q = he" (B.6) 36 SQ which can be written i n m a t r i x form as Isq cos(0) - sin(0) sin(0) " cos(0) _ ^sct (B.7) The real component of the current space vector i n 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. T h e main advantage of the dq transformation is the elimination of position dependency from the machine electrical variables. The inverse dq transformation can be written i n m a t r i x form as ^sct " cos(0) sin(0) - sin(0) ' cos(0) isd 1sq (B.8) The dq transformation performed directly on the phase quantities has the following well-known form: id Isq s 2 ~ 3 cos (0) - sin (0) cos (0 - f ) -sin ( 0 - f ) 133 cos (0 - f ) -sin(0-f) ib ir (B.9) Appendix C Digital Pass-Band Filter for Filtering Stator-Flux To accurately calculate of the stator-flux vector position 9 f i n Section 5.6.3, it is necessary to filter out D C offsets i n the stator flux. Instantaneous values of stator flux linkages ip , ipi, and ip are filtered using the digital pass-band filter described here. A simple R L C parallel circuit i n 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. s a c Figure C . l : C i r c u i t configuration for filtering 60 H z component of ^ ( t ) The admittance of this circuit can be expressed as (C.l) In the resonant condition, the susceptance is equal to zero, which leads to the following well-known expression for the resonant frequency: 1 CO = UQ = (C.2) T h e resonant frequency is entirely specified by the choice of Lp and CF- T h e magnitude of the filtered flux linkage ^ _ w i l l vary as a function of the frequency of the input flux linkage \& , as depicted i n F i g . C . 2 . a F a 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 : V a r i a t i o n 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 Qo = LOBW — (C.3) LOQRFC'F W i t h the resonant frequency set to 50 H z (314 rad/s) for the particular example in Chapter 5 and w i t h the desired quality factor, the values of RF, L and CF can be determined from the above equations by freely choosing either one of them. F T h e circuit i n F i g . C . l can be discretized using the trapezoidal rule of integration. The equivalent discrete circuit is depicted i n F i g . C . 3 , where h {t) LF = - ^ ^ ( t - A t ) + h (t - At) LF h {t) = ^aAt ~ A t ) CF (CA) h (t CF - At) B y solving the circuit, the filtered stator-flux linkage of phase a is 1pa.F(t) + h (t) LF 135 + h {t) CF (C.5) Appendix C. Digital Pass-Band Filter for Filtering Stator^Flux Figure C . 3 : Equivalent circuit of the parallel R L C filter discretized w i t h the trapezoidal 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 Nm Rated current: 421.2 A Inertia: J = 63.87 kg m R = 0.029 ft R' = 0.022 ft X = 0.226 ft X' = 0.226 ft X - 13.04 ft 3 2 s r u lT m Table D . l : 2250 hp induction motor parameters 137 Appendix D. Test System Parameters D.2 Parameters of the 835 M V A Synchronous Generator from Chapter 4.4.2 R a t i n g : .835 MVA Line-to-line voltage: 26 kV Poles: 2 Speed: 3600 rpm Combined inertia of generator and turbine: J = 0.0658 • 10 J s R = 0.00243 fi Xi = 0.1538 fi X = 1.457 fi X = 1.457 fi R = 0.00144 fi R' = 0.00075 fi X' = 0.6578 fi X' = 0.1145 fi R' = 0.00681 fi R' = 0.01080 fi X' = 0.07602 fi X' = 0.06577 fi 6 s s q d kql fd lkql lfd kq2 kd lko2 lkd Table D . 2 : 835 M V A steam turbine generator parameters 138 2 Appendix D. D.3 Test. System Parameters Parameters of the 7.5 kW Doubly Fed Induction Generator from Chapter 5.10 R a t e d power = 7.5 k W Stator voltage = 415 V Rotor voltage = 440 V R a t e d stator current = 19 A R a t e d rotor current = 11 A Pole pairs =3 R a t e d speed = 970 rmp Base frequency = 50 H z N /N = 1.7 J = 7.5kgm stator connection = delta rotor connection = wye s r 2 Machine resistances and inductances per phase: R = 1.06 ft Rr = 0.80 ft L = 0.0664 H L = 0.0810 H L = 0.0320 H s s 0 r Control parameters: Stator-side converter Rotor-side converter K =0.12 tf =0.49 a„=0.9248 a =0.988 #,-=4.72 # =20 Oi=0.96 a =0.985 v w w i r ir Table D . 3 : 7.5 k W doubly fed induction generator w i n d 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
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 |
Aggregated Source Repository | 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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0100317/manifest