Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Tarazona Gomez, Javier Orlando
Doctor of Philosophy - PhD
University of British Columbia
Traditional power systems have high inertia, and disturbances can be represented with relatively large discretization steps. However, power electronic interfaced generation requires time steps, at least two orders of magnitude smaller to represent the switching operations. In addition, the increased penetration of these devices reduces the total inertia of the system, and a single-time step size transient stability (TS) solution that considers the details of the power electronic devices would take two orders of magnitude longer (that is, electromagnetic transients (EMT) simulators time step sizes). Hybrid TS-EMT simulators aim at using large time steps for the traditional (slow) part of the system and small time steps for the power electronic devices. However, in the current state-of-the-art, the proposed hybrid solutions are limited in a number of ways. Most of them introduce a one TS time step interfacing delay between solutions, which reduces accuracy and can create numerical oscillations. Other approaches use the travelling time of transmission lines to tear the system and interface the solutions, but this limits the flexibility of connecting power electronic components anywhere in the system. The hybrid solution presented in this work combines a Shifted Frequency Analysis (SFA) solution that is more accurate than traditional phasor solutions to follow transient stability trajectories in lower inertia systems with the traditional EMT solution for the fast power electronic converters. In order to interface the SFA and EMT solutions, a new hybrid protocol has been developed that overcomes the disadvantages of the current hybrid solutions. The two subsystems are joined using the Multi-Area Thévenin Equivalent (MATE) solution framework, and introducing the novel concept of employing a second parallel EMT simulation. Illustrative examples demonstrate the validity and accuracy of the SFA-EMT multirate protocol. The new protocol achieves significant computational savings when compared with an all-EMT solution. A simulator operating with the proposed protocol will allow improved accuracy when conducting power system transient studies with high penetration of renewable generation without considerably increasing the computational time.
HYBRID SHIFTED FREQUENCY ANALYSIS-ELECTROMAGNETIC TRANSIENTS MULTIRATE SIMULATION OF POWER SYSTEMS by Javier Orlando Tarazona Gomez B.Eng., Universidad Central de Venezuela, 1986 M.Sc., University of Manchester, 1990 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2023 © Javier Orlando Tarazona Gomez, 2023 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: Hybrid Shifted Frequency Analysis-Electromagnetic Transients Multirate Simulation of Power Systems submitted by Javier Orlando Tarazona Gomez in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Electrical and Computer Engineering Examining Committee: José R. Martí, Professor, Electrical and Computer Engineering, UBC Supervisor Hermann Dommel, Professor, Electrical and Computer Engineering, UBC Supervisory Committee Member Juri Jatskevich, Professor, Electrical and Computer Engineering, UBC Supervisory Committee Member Zhen Jane Wang, Professor, Electrical and Computer Engineering, UBC University Examiner Ryozo Nagamune, Professor, Mechanical Engineering, UBC University Examiner iii Abstract Traditional power systems have high inertia, and disturbances can be represented with relatively large discretization steps. However, power electronic interfaced generation requires time steps, at least two orders of magnitude smaller to represent the switching operations. In addition, the increased penetration of these devices reduces the total inertia of the system, and a single-time step size transient stability (TS) solution that considers the details of the power electronic devices would take two orders of magnitude longer (that is, electromagnetic transients (EMT) simulators time step sizes). Hybrid TS-EMT simulators aim at using large time steps for the traditional (slow) part of the system and small time steps for the power electronic devices. However, in the current state-of-the-art, the proposed hybrid solutions are limited in a number of ways. Most of them introduce a one TS time step interfacing delay between solutions, which reduces accuracy and can create numerical oscillations. Other approaches use the travelling time of transmission lines to tear the system and interface the solutions, but this limits the flexibility of connecting power electronic components anywhere in the system. The hybrid solution presented in this work combines a Shifted Frequency Analysis (SFA) solution that is more accurate than traditional phasor solutions to follow transient stability trajectories in lower inertia systems with the traditional EMT solution for the fast power electronic converters. In order to interface the SFA and EMT solutions, a new hybrid protocol has been developed that overcomes the disadvantages of the current hybrid solutions. The two subsystems are joined iv using the Multi-Area Thévenin Equivalent (MATE) solution framework, and introducing the novel concept of employing a second parallel EMT simulation. Illustrative examples demonstrate the validity and accuracy of the SFA-EMT multirate protocol. The new protocol achieves significant computational savings when compared with an all-EMT solution. A simulator operating with the proposed protocol will allow improved accuracy when conducting power system transient studies with high penetration of renewable generation without considerably increasing the computational time. v Lay Summary Current power system networks are undergoing a reduction of inertia due to the coexistence of significant renewable energy resources with traditional generators. Shifted Frequency Analysis (SFA) is a simulator that performs significantly better than traditional tools to study the stability of low inertia systems. Electromagnetic Transients (EMT) simulators are used to simulate fast devices, like renewable energy generators. This research has developed a new hybrid simulator that combines the advantages of the SFA and EMT solutions. The work introduces a novel interfacing technique that avoids the limitations of previously proposed hybrid solutions. The thesis demonstrates the proposed protocol's mathematical validity, accuracy, and numerical stability. vi Preface Part of the research results presented in this dissertation have been published in a scientific journal and have appeared in a conference proceedings. Additional results will be submitted for peer review. In all publications, I have been responsible for deriving all the mathematical formulations, developing and implementing the models proposed, carrying out the computer studies, analyzing the results, and writing the first draft of the publication. My supervisor, Dr. José Martí, has provided guidance and instructive comments and observations during the whole research process and in the preparation of the manuscripts and their revisions. Other co-authors have also contributed to the research by providing their feedback and discussions, editorial comments, and by reviewing the manuscripts. The simulator developed in this research was made from an existing software created by A. T. J. Martí for the individual EMT and SFA solvers. The enhancements to the modeling capabilities of those initial solvers as well as the development of the integrated hybrid SFA-EMT simulator was carried out by myself. The first article that has resulted from this doctoral work is indicated below and was based on part of the results presented in Chapters 4 and 6 of this dissertation. J. O. Tarazona, A. T. J. Martí, J. R. Martí, and F. A. Moreira, “Shifted frequency analysis-EMTP multirate simulation of power systems,” Electric Power Systems Research, vol. 197, p. 107292, 2021, doi: https://doi.org/10.1016/j.epsr.2021.107292. vii This article was also presented at the 2021 International Conference on Power Systems Transients (IPST 2021) in Belo Horizonte, Brazil, June 07-10. viii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ....................................................................................................................... viii List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................. xiii Nomenclature ............................................................................................................................ xvii List of Abbreviations ............................................................................................................... xviii Acknowledgements .................................................................................................................... xix Dedication ................................................................................................................................... xxi Chapter 1: Introduction ................................................................................................................1 1.1 Background ..................................................................................................................... 1 1.1.1 Power Systems Transients .......................................................................................... 1 1.1.2 Computation and Analysis of Electromagnetic Transients ......................................... 2 1.1.3 Power System Stability ............................................................................................... 3 1.1.4 Power System Stability Analysis ................................................................................ 3 1.1.5 Transient Stability-Electromagnetic Transients Hybrid Simulators ........................... 4 1.1.6 Dynamic Phasors ........................................................................................................ 6 1.1.7 Shifted Frequency Analysis ........................................................................................ 7 1.1.8 Hybrid Simulators Using Dynamic Phasors ............................................................... 8 1.2 Motivation and Research Objectives .............................................................................. 9 ix Chapter 2: Shifted Frequency Analysis of Power Systems ......................................................13 2.1 Theoretical Foundations................................................................................................ 14 2.1.1 The Dynamic Phasor ................................................................................................. 14 2.1.2 The Analytic Signal .................................................................................................. 15 2.1.3 The SFA Transformation .......................................................................................... 17 2.2 SFA Models of Basic Network Elements ..................................................................... 19 2.2.1 Dynamic Phasors in an Electric Circuit .................................................................... 19 2.2.2 Continuous-time SFA Equivalent Circuits of Basic R-L-C Elements ...................... 21 2.2.3 Discrete-time SFA Equivalents of M-phase Coupled Resistances ........................... 23 2.2.4 Discrete-time SFA Equivalents of M-phase Coupled Inductances ........................... 24 2.2.5 Discrete-time SFA Equivalents of M-phase Coupled R-L Branches ....................... 24 2.2.6 Discrete-time SFA Equivalents of M-phase Coupled Capacitances ......................... 25 2.2.7 Structure of the R-L-C Discrete-time SFA Models .................................................. 26 2.2.8 SFA models of Other Circuit Elements .................................................................... 27 2.3 The SFA Solution ......................................................................................................... 27 2.3.1 General Procedure ..................................................................................................... 27 2.3.2 SFA for Electromagnetic Transients ......................................................................... 29 2.3.3 Nonstandard SFA Solutions ...................................................................................... 30 2.4 Electric Power in Complex and Shifted Frequency Domains ...................................... 31 Chapter 3: Simulation of Power Systems in a Partitioned Environment ...............................38 3.1 Multi-Area Thevenin Equivalent Conceptual Framework ........................................... 39 3.2 MATE Solution Algorithm ........................................................................................... 42 3.3 Use of MATE in EMT-EMT Multirate Simulation ...................................................... 44 x 3.4 EMT-EMT MATE Multirate Interfacing Protocol ....................................................... 48 Chapter 4: SFA-EMT Multirate Simulation .............................................................................50 4.1 Interfacing SFA and EMT Simulators in MATE .......................................................... 50 4.2 Single Rate SFA-EMT Simulation ............................................................................... 52 4.2.1 Links Solution ........................................................................................................... 53 4.2.2 Subsystems Solutions ............................................................................................... 55 4.2.3 Single Rate SFA-EMT Solution Algorithm .............................................................. 56 4.3 Multirate SFA-EMT Simulation ................................................................................... 58 4.3.1 Updating the EMT Subsystem .................................................................................. 58 4.3.2 Updating the SFA Subsystem ................................................................................... 60 4.3.3 Antialiasing Processing of EMT Solution Variables ................................................ 62 4.3.3.1 Standard SFA-EMT MATE Multirate Protocol ............................................... 63 4.3.3.2 Alternative SFA-EMT MATE Multirate Protocol ............................................ 67 4.4 Selection of Boundaries Between Subsystems ............................................................. 72 Chapter 5: Description of simulators employed .......................................................................75 5.1 General Characteristics of the Hybrid Simulator .......................................................... 75 5.2 The Link Solver ............................................................................................................ 76 5.3 The SFA and EMT Solvers ........................................................................................... 78 Chapter 6: Case Studies ..............................................................................................................80 6.1 Case 1: IEEE 39 Bus Test System Under Fault Conditions ......................................... 80 6.1.1 Case 1a: Hybrid SFA-EMT Single Rate Solution .................................................... 82 6.1.2 Case 1b: Hybrid SFA-EMT Multirate Solution ........................................................ 86 6.2 Case 2: Modified CIGRE HVDC Benchmark System ............................................... 100 xi 6.2.1 Case 2a: Hybrid SFA-EMT Single Rate Solution .................................................. 103 6.2.2 Case 2b: Hybrid SFA-EMT Multirate Solution ...................................................... 107 6.2.3 Case 2c: Hybrid SFA-EMT Multirate Solution in Normal Operation .................... 111 Chapter 7: Conclusion ...............................................................................................................115 7.1 Contributions............................................................................................................... 115 7.2 Future Work ................................................................................................................ 121 References ...................................................................................................................................125 Appendix .....................................................................................................................................137 Data of the IEEE 39-Bus Test System .................................................................................... 137 xii List of Tables Table 6.1 Error in first peak of fault current from bus 28 (phase "b") for various antialiasing methods. ........................................................................................................................................ 93 Table 6.2 Execution times for simulations of Case 1b ................................................................. 96 Table 6.3 Estimated computational gain of the SFA-EMT solution with respect of an all EMT solution, for various relative sizes of the SFA and EMT subsystems and various rate ratios of the multirate solution. ......................................................................................................................... 98 xiii List of Figures Figure 2.1 Amplitude spectra of a signal in different representations for the SFA process: a) Real bandpass signal, b) Analytic signal and c) SFA dynamic phasor. ................................................ 18 Figure 2.2 One-port circuit to study the concept of power in complex domain, normal time and SFA domains ................................................................................................................................. 32 Figure 3.1 Sample system for MATE explanation ....................................................................... 40 Figure 3.2 Interpolation of Thévenins of the slow subsystem to calculate the link solutions required to update the fast subsystem at various solution times. .................................................. 45 Figure 3.3 Amplitude spectrum of the Fourier transform of the Thévenin of the fast subsystem, with indication of the range of frequencies that would produce aliasing in the slow subsystem solution. ......................................................................................................................................... 47 Figure 3.4 MATE Interfacing protocol for an EMT-EMT multirate solution .............................. 49 Figure 4.1 System split into two subsystems with “n” links for SFA-EMT simulation with MATE ........................................................................................................................................... 53 Figure 4.2 Standard SFA-EMT MATE multirate protocol ........................................................... 62 Figure 5.1 General flow chart of the integrated hybrid SFA-EMT simulator. ............................. 77 Figure 6.1 IEEE 39-bus system and SFA and EMT subsystems for case study 1. ....................... 81 Figure 6.2 Frequency of generator 8 in the SFA subsystem (single rate solution). ...................... 83 Figure 6.3 Electrical power of generator 8 in the SFA subsystem (single rate solution). ............ 83 Figure 6.4 Voltage at bus 27 (phase “a”) in the SFA subsystem (single rate solution). ............... 84 Figure 6.5 Fault current from bus 28 (phase “a”) in the EMT subsystem (single rate solution). . 84 xiv Figure 6.6 Fault current from bus 28 (phase “b”) in the EMT subsystem (single rate solution). . 85 Figure 6.7 Fault current from bus 28 (phase “c”) in the EMT subsystem (single rate solution). . 85 Figure 6.8 Voltage at bus 28 (phase “c”) in the EMT subsystem (single rate solution). .............. 86 Figure 6.9 Frequency of generator 8 in the SFA subsystem (multirate solution). ........................ 87 Figure 6.10 Angle difference between generators 8 and 1 in the SFA subsystem (multirate solution). ....................................................................................................................................... 88 Figure 6.11 Electric power of generator 8 in the SFA subsystem (multirate solution). ............... 88 Figure 6.12 Voltage at bus 27 (phase “a”) in the SFA subsystem (multirate solution). ............... 89 Figure 6.13 Fault current from bus 28 (phase “a”) in the EMT subsystem (multirate solution). . 90 Figure 6.14 Fault current from bus 28 (phase “b”) in the EMT subsystem (multirate solution). . 90 Figure 6.15 Fault current from bus 28 (phase “c”) in the EMT subsystem (multirate solution). . 91 Figure 6.16 Voltage at bus 28 (phase “a”) in the EMT subsystem (multirate solution). .............. 91 Figure 6.17 Voltage at bus 28 (phase “b”) in the EMT subsystem (multirate solution). .............. 92 Figure 6.18 Voltage at bus 28 (phase “c”) in the EMT subsystem (multirate solution). .............. 92 Figure 6.19 Estimated computational gain of the SFA-EMT multirate solution with respect of an all EMT solution, as a function of the size of the SFA subsystem and the rate ratio of the multirate solution (This estimation does not consider the execution time of the link solver for the multirate solution). ........................................................................................................................ 99 Figure 6.20 Modified CIGRE HVDC benchmark system and SFA and EMT subsystems for case study 2 ......................................................................................................................................... 101 Figure 6.21 Voltage at bus 1 (phase “b”) in the SFA subsystem (single rate solution). ............. 104 Figure 6.22 Current from AC system into bus 1 (phase “a”) in the EMT subsystem (single rate solution). ..................................................................................................................................... 104 xv Figure 6.23 DC voltage at rectifier terminal in the EMT subsystem (single rate solution). ....... 105 Figure 6.24 DC voltage at inverter terminal in the EMT subsystem (single rate solution). ....... 105 Figure 6.25 DC current at rectifier terminal in the EMT subsystem (single rate solution). ....... 106 Figure 6.26 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (single rate solution). ..................................................................................................................................... 106 Figure 6.27 Voltage at bus 1 (phase “b”) in the SFA subsystem (multirate solution). ............... 108 Figure 6.28 Current from AC system into bus 1 (phase “b”) in the EMT subsystem (multirate solution). ..................................................................................................................................... 108 Figure 6.29 DC voltage at rectifier terminal in the EMT subsystem (multirate solution). ......... 109 Figure 6.30 DC voltage at inverter terminal in the EMT subsystem (multirate solution). ......... 109 Figure 6.31 DC current at rectifier terminal in the EMT subsystem (multirate solution). ......... 110 Figure 6.32 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (multirate solution)...................................................................................................................................................... 110 Figure 6.33 Voltage at bus 1 (phase “b”) in the SFA subsystem (steady state multirate solution)...................................................................................................................................................... 112 Figure 6.34 Current from AC system into bus 1 (phase “b”) in the EMT subsystem (steady state multirate solution). ...................................................................................................................... 112 Figure 6.35 DC voltage at rectifier terminal in the EMT subsystem (steady state multirate solution). ..................................................................................................................................... 113 Figure 6.36 DC voltage at inverter terminal in the EMT subsystem (steady state multirate solution). ..................................................................................................................................... 113 Figure 6.37 DC current at rectifier terminal in the EMT subsystem (steady state multirate solution). ..................................................................................................................................... 114 xvi Figure 6.38 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (steady state multirate solution). ...................................................................................................................... 114 xvii Nomenclature In this dissertation complex quantities are represented with barred characters and real quantities with normal characters. The exception to this rule is the SFA transformation that is denoted by a plain capital T. Vectors and matrices are denoted by bold characters and scalar quantities are represented with normal type characters. Most of the time dynamic phasors or complex signals in the SFA domain are represented with the subscript “m”, but in a few occasions they are denoted by uppercase letters without a subscript. Throughout the dissertation the meaning of all variables is properly explained, however, in this section only a few general and basic variables are included, since their meaning is assumed in the text. I : Identity matrix 0 : Null matrix xviii List of Abbreviations AC Alternating Current CDA Critical Damping Adjustment CIG Converter Interfaced Generation DC Direct Current DP Dynamic Phasor EMT Electromagnetic Transients EMTP Electromagnetic Transients Program FACTS Flexible Alternating Current Transmission Systems FAST Frequency Adaptive Simulation of Transients FDNE Frequency Dependent Network Equivalent GAM Generalized Averaging Method HVDC High Voltage Direct Current MATE Multi-Area Thévenin Equivalent MMC Modular Multilevel Converter rms Root mean square SFA Shifted Frequency Analysis S-G Savitzky-Golay TICP Total Instantaneous Complex Power TIAP Total Instantaneous Active Power TIRP Total Instantaneous Reactive Power TS Transient Stability xix Acknowledgements I want to express my sincere and deep gratitude to my supervisor Dr. José R. Martí, for sharing his knowledge, insights and intuitions with me, for his constant and unconditional guidance, support, inspiration, motivation, patience, kindness, understanding and above all for his trust. I am also grateful to my research committee members, Dr. Hermann W. Dommel and Dr. Juri Jatskevich for their help and cooperation during this research work. I thank also Dr. Luis Linares for sharing his knowledge and ideas and for his support and guidance. In addition, I want to thank Dr. José R. Martí for providing funding for this research through the Discovery Grant of the Natural Sciences and Engineering Research Council (NSERC) of Canada. Likewise, I am grateful to the University of British Columbia for its financial support during my studies, especially during the pandemic. Special thanks go to Andrea Martí, for adjusting, thoroughly documenting and allowing me to use the program she developed for her MASc project, and for sharing with me her solid and extensive Python knowledge. Special thanks also for Dr.Verónica Galván for her support and help while she was a visiting PhD student at UBC and later from Mexico. Likewise, I want to specially thank Dr. Fernando Moreira, for his continuous and unconditional help, motivation and support while he was a visiting professor at UBC and later from Brazil. I also want to thank Dr. Yankan Song for trying my first version of my hybrid protocol and providing feedback on its performance. xx I also want to thank my friends, fellow colleagues and graduate students in the Electric Power and Energy Systems research group and in the i2SIM group at UBC: Andrea, Ivelina, Verónica, Ehssan, Yang, Arash, César, Jiafeng, Alireza, Felix, Yingwey, Navid, Milad, Ahmed, Merhdad, Matin, Wonbae, Zaeed, and many more not included in this not at all exhaustive short list, for many stimulating discussions, for their help and support and for making my program and my stay at UBC an enjoyable and memorable experience that I will never forget. And last, but not least, I want to thank my wife Carmen and my sons Javier and Arturo for their love and support, for their patience and commitment with this project and for being my inspiration and motivation to accomplish this life goal. I also want to thank my brothers and sisters, William, Miguel, Beatriz and Luisa and their families, for their love, support and motivation from the distance. May God bless you all! Javier O. Tarazona G. xxi Dedication To Carmen, Javier and Arturo 1 Chapter 1: Introduction 1.1 Background 1.1.1 Power Systems Transients In general, a power system transient is the process that a power system undergoes when is forced, by an external or internal event or disturbance, to change from a certain operating state, usually steady state, to another state. During a transient, the voltages, currents and other parameters in the system suffer temporary changes. Power system transients can be broadly classified into Electromagnetic Transients and Electromechanical Transients. Electromagnetic transients are initiated by the abrupt change in the operating conditions of the system due to external electromagnetic fields, equipment failure or changes in the configuration of the network such as switching and loading. Electromagnetic transients in power systems result from a combination of transients deriving from the response of lumped-parameters elements and traveling waves in distributed-parameter elements, either type can predominate over the other depending on the characteristics and conditions of the system and time of observation [1]. In general, they are fast phenomena and decay rapidly in the system. Electromechanical transients involve mechanical oscillations among the turbine generators of the system. These oscillations are relative to the synchronous speed and can last for seconds, minutes and possibly longer, depending on the case [2]. Thus, these transients are slower than electromagnetic transients. Electromechanical transients are caused by disturbances or events 2 that involve a mismatch, temporary or permanent, between the generation and consumption of power, or a rearrangement in the distribution of power within the network. The majority of these transients lie in the realm of what is known as power system stability. As opposed to mismatches between generation and load, EMT transients are due to fast changes in the topology of the circuit, like in power electronics commutation. 1.1.2 Computation and Analysis of Electromagnetic Transients Electromagnetic transients (EMT) analysis refers to the evaluation of the electromagnetic stresses to which a piece of equipment or a group of equipment of the system is subjected to. This requires the calculation of the temporary voltages and currents caused by switching operations, faults, lightning surges and other intended and unintended disturbances [3]. This kind of studies are used in the design stage of the system and its different pieces of equipment or for trouble-shooting purposes in the operation of the system. The preferred tools worldwide to perform this type of studies for power systems are the EMTP-type of simulators (EMTP stands for ElectroMagnetic Transients Program), such as PSCAD/EMTDC [4], EMTP® [5], ATP [6] and MicroTran [7]. Due to the objectives of these studies and the nature of the phenomena involved, this type of tools allows for the use of detailed models for the different types of equipment so as to fully capture all the frequencies involved in the phenomenon to be analyzed, and thereby simulate with good accuracy the response of the system or equipment under study. This requires the solution of a large number of differential equations and the use of small time-steps, which in turn means that the computational burden associated with these studies is high and, therefore, to 3 maintain reasonable solution times the size of the system is usually reduced to the smallest size that would render sufficiently accurate results for the type of study under consideration. 1.1.3 Power System Stability “Power system stability is the ability of an electric power system, for a given initial operating condition, to regain a state of operating equilibrium after being subjected to a physical disturbance, with most system variables bounded so that practically the entire system remains intact” [8]. It can be classified into 3 main categories: Rotor Angle Stability, Voltage Stability and Frequency Stability. Rotor Angle Stability is related to the ability of the synchronous machines of a power system to remain in synchronism after being subjected to a disturbance. It can be further classified into: Large-disturbance Rotor Angle Stability (Transient Stability) and Small-disturbance (Small-signal Rotor Angle Stability) [8], [9]. Voltage Stability is related to the ability of a power system to maintain steady voltages at all buses after being subjected to a disturbance [8], [9]. Frequency stability is concerned with the ability of a power system to maintain steady frequency after a severe system disturbance that results in a significant imbalance between load and generation [8]. 1.1.4 Power System Stability Analysis In stability studies the fast and mostly local transients associated to the disturbances are not usually of interest from a system point of view. The interest is not the integrity of one piece of equipment or a group of equipment, but rather the integrity of the system as a whole. Thus, in this type of analysis fully detailed models of the electrical components are not commonly 4 employed and the evaluations are conducted by assuming that the voltages and currents in the system are more or less at power frequency (60 or 50 Hz) and therefore, the electrical network is represented with steady-state voltage and current phasors [2]. Since electromechanical oscillations are slower phenomena than electromagnetic transients, the time span of interest in power system stability analysis is usually much longer and much larger integration time steps can be used. For transient stability the time interval of interest is typically 3 to 5 s, but it can go up to 20 s for large interconnected systems with dominant inter-area swings [8], [9]; the integration time step is in the order of the milliseconds. This type of analysis requires a complete description of the power network, so the size of the system to be modeled is usually quite large, which in turn means that the computational load is very heavy. However, the use of simplified models, together with the use large integration time-steps, helps in keeping the computational burden manageable while rendering sufficiently accurate results for this kind of studies. Commonly used tools worldwide to perform this type of studies are PSS/E [10], PowerFactory [11], EUROSTAG [12], TSAT [13] and PSASP [14]. 1.1.5 Transient Stability-Electromagnetic Transients Hybrid Simulators Computational tools for transient stability and electromagnetic transients studies had traditionally been used separately. However, with the increase in the complexity of power systems, due to the introduction of new technologies and its controls, and an increasing number of interconnections, the situation began to change. It was with the advent of the High Voltage Direct Current (HVDC) transmission that the need first arose for the incorporation of the kind of models available in EMTP-type of programs in the traditional tools used for stability studies. However, 5 the model simplifications and the large integration steps used by TS-tools precluded that possibility; thus, the need for a hybrid simulation framework between a transient stability type of simulation and an EMTP type of simulation became evident. The detailed modeling of an HVDC link within a power system stability simulator was first proposed and implemented by Hefferman et al. in 1981 [15]. The DC part was simulated using a state space formulation and the AC part was simulated with a conventional stability program. Later, Reeve and Adapa [16] developed a similar simulator for an HVDC link application, but this time, the detailed system (HVDC link plus part of the AC network) was simulated with the EMTP [17] while the external system (the rest of the AC network) was modeled again with a conventional transient stability program. After that, several publications have been presented in the literature dealing with different variations and approaches for the implementation of hybrid TS-EMT simulators and their applications. In [18] a comprehensive review of the literature on hybrid TS-EMT simulators is presented. In [19], [20], [21], [22] and [23] a Frequency Dependent Network Equivalent (FDNE) is proposed to represent the external AC system (TS part) in the detailed system (EMT part), to take into account the harmonic distortion at the interface. Hybrid TS-EMT simulators for real-time applications are presented in [24], [25], [22] and [26]. The hybrid simulators presented in [27], [23], [28] and [29] use interfacing protocols that iterate between the TS and EMT models to find the solution of the whole system, unlike all previous hybrid simulators based on direct interaction protocols, either serial or parallel, that do not iterate between the two models and therefore can introduce up to a one time step delay in the multirate solution. Some authors [30] 6 and [31], have applied a hybrid TS-EMT framework for the simulation of offshore wind power connected through a multi-terminal voltage sourced converters-high voltage direct current (VSC-HVDC) and have evaluated variations and improvements in the interaction protocols and the exchange of data between the simulators. In [32] a TS-EMT simulator that permits the switching from a parallel to a serial interaction protocol and combines a commercial EMT simulator with a three-sequence based TS simulator, was applied in a fault-induced delayed voltage recovery study. All hybrid TS-EMT simulators mentioned previously are based on conventional TS programs, either positive-sequence or three-sequence based, that use conventional or standard AC phasors. 1.1.6 Dynamic Phasors Dynamic phasors (DP) is a modeling technique originally based on the Generalized Averaging Method (GAM) [33]. The GAM consists in representing the variables in the circuit by time-dependent Fourier series evaluated over a sliding window and retaining only the most important harmonics of those series. The corresponding complex Fourier coefficients are called dynamic phasors. This averaging method can be applied to any periodically (or almost periodically) driven circuit [33] and allows the simulation of electromagnetic dynamics with adjustable accuracy depending on how many harmonics are retained for the analysis. DP (GAM-type) based simulations allow larger time steps than standard time domain simulations [34] since the DP models need to capture only up to the highest harmonic of interest rather than the highest frequency of the circuit. The DP modeling technique has been used to model induction and synchronous machines [35], [36], asymmetrical faults [34], power electronics components [36], 7 [37] and many other types of power systems components. Power system simulators have also been developed based on this type of modeling [36], [38] and [39]. The GAM approach is an option for modeling power systems dynamics that is sort of a midway between conventional transient stability modeling and standard detailed time domain modeling. Building equivalent circuits of the individual power system components is not straightforward with GAM-type DP modeling, this is more readily done, however, with Shifted Frequency Analysis, as it is discussed below. 1.1.7 Shifted Frequency Analysis Shifted Frequency Analysis (SFA) is based on the knowledge that following an event in a power system, after the electromagnetic transients have died out, the remaining transients (the electromechanical ones) manifest themselves as bandpass signals around the fundamental frequency of the system (50/60 Hz) [40]. That is, the frequency spectrum of the variables in the system is composed of a band around the system frequency. By shifting by 50/60 Hz the frequency spectrum of the variables through the use of a time domain rotational transformation, the resulting spectrum of the transformed variables is centered around zero frequency, thereby allowing the use of large time steps to trace the envelopes of the waveforms without compromising the accuracy of the solution [40], [41]. The transformed variables in the SFA domain can be considered as fundamental frequency dynamic phasors [40], but the methods to solve for these variables are based on frequency translations and time discretization. By discretizing the differential equations of the system in the same fashion as with the EMTP [17], equivalent circuits can be obtained for the individual system components and the network 8 solution can be found in a similar way as with the EMTP [40], [42], [43]. SFA models have been developed for synchronous machines [44],[45], induction machines[46], transformers and transmission lines [47] and for power electronic devices such as modular multi-level converters (MMC’s) [48]. This framework for modeling and solving the network is known as shifted frequency analysis [40]. The SFA dynamic phasors evolve or converge to the standard AC phasors when the system reaches steady state, so classical AC phasor analysis can be seen as a particular case of SFA after the transients have died out. 1.1.8 Hybrid Simulators Using Dynamic Phasors Due to the accuracy and computational benefits offered by the dynamic phasors, in the last few years several hybrid simulators that make use of dynamic phasors have been published in the literature. For example, in [49] and [50] simulators interfacing a dynamic phasor solver with an EMT solver are presented. Both simulators require the use of transmission lines as the interface between the solvers in order to allow the operation in multirate mode. [49] uses a constant parameter transmission line for that purpose and [50] can use both, constant parameters and frequency dependent transmission lines for the interface. The simulator in [50] is applied to study sub- and super-synchronous interactions in large scale wind farms. [51] presents the application in a practical HVDC project, of a TS-EMT simulator where both solvers are interfaced with a lump parameter pi-circuit model of a transmission line based on dynamic phasors (GAM-type). In [52] and [53] three solvers, one TS program, a dynamic phasor (SFA-type) solver and an EMT solver are combined into one single simulator TS-DP-EMT, the part of the system modeled with the DP solver serves as a buffer zone between the other two parts. The simulator in [52] is for real-time applications but requires balanced operating conditions in the 9 system. The simulator in [53] is more general and uses a 3-sequence TS program and is not limited to balanced operating conditions. It presents an application to the study of a practical HVDC project. Other research in hybrid simulators has centred on developing an integrative approach to model TS and EMT tasks. Instead of coupling two different simulation frameworks, the solution is based on the concept called Frequency Adaptive Simulation of Transients (FAST), to cover the application spectrum of typical EMT and TS programs, but the TS part is done using a dynamic phasor of the SFA-type [54]. In this approach, the entire system is first solved for fast transients in EMT mode with a very small time-step size, which is not very efficient, and then solved in TS mode for the electromechanical transients in the entire system, with a much larger time step size and using dynamic phasors. In [55] a multirate co-simulator that interfaces a FAST solver with an EMT solver is presented. In this implementation the FAST solver is called BFAST, because it employs what they call base-frequency dynamic phasors, which use a combination GAM- and SFA-type dynamic phasors. As with other simulators, the multirate operation of this BFAST-EMT co-simulator is achieved through the use of transmission lines at the partitioning locations to decouple the two solvers. 1.2 Motivation and Research Objectives Since the early 80’s the power engineering community has needed a simulator that, as efficiently as a TS simulator and as accurately as an EMT simulator, could model large-scale power systems to perform system wide functional analyses such as stability analysis. Nowadays, with the ever-increasing penetration of power electronics interfaced devices and distributed energy resources 10 in the modern power grid, the need for such a simulation framework has become more urgent. Conventional synchronous generators are being replaced by converter interfaced generation, which reduces the system’s inertia and affects the dynamic phenomena, making the system to experience higher frequency deviations and higher rates of change of frequency under disturbances, among other effects [56]. Because of this, TS programs are becoming less and less adequate to study the dynamic phenomena in modern power systems, due to the approximations involved in the modeling with such programs. An SFA simulator has the capability of using large time steps of the same order of magnitude as the TS simulators and since SFA solves the differential equations associated with the electrical network, it can trace more accurately the trajectory of the dynamic system states, even in operating conditions close to collapse, and therefore it is expected to perform better than a TS simulator in low inertia systems, and in fact, recent comparative studies [57], [58] [59] showed that SFA delivers better accuracy than a TS simulator for both, classic and low inertia systems. On the other hand, the option of using pure EMT simulation to study the stability of modern large-scale power systems, though the most accurate, it is still not computationally efficient due to the size of the systems and the proliferation of HVDC and power electronics interfaced devices, which entails the use of very small time steps. Based on the above, interfacing an SFA solver and an EMT seems a step on the right direction to achieve a hybrid simulator more accurate than hybrid TS-EMT simulators and more efficient than an EMT simulator. Such an SFA-EMT simulator would have these additional benefits. 11 Since the SFA solution algorithm is the same as of an EMTP-type of simulator, prior to actually solving the system, it is possible to get a Thévenin equivalent for the next SFA solution instant and this feature can be exploited through the use of the conceptual solution framework known as Multi-Area Thévenin Equivalent (MATE), as in [60] and [61], to allow for the direct and simultaneous interfacing of the SFA and EMT solutions, avoiding the problem of one TS time step delay that most TS-EMT simulators have, except the ones that can iterate between the TS and EMT solutions, like the simulators presented in [29], [27] and [23]. This direct way of interfacing the SFA and EMT solutions does not require transmission lines to decouple the two solutions, as is the case in some hybrid simulators that use dynamic phasors, for example [49], [50] and [55], which imposes an important restriction because the maximum time step that can be used in the dynamic phasor solution is limited by the travel time of the line. Since SFA solves for the dynamic phasors in the entire external system, variations in magnitude, angle and frequency are naturally taken into account and therefore no assumptions are made regarding the frequency of the equivalent impedance of the SFA subsystem, unlike many traditional TS-EMT simulators that evaluate the TS subsystem equivalent impedances at the fundamental frequency. The SFA solution is naturally in abc coordinates, so it can handle any unbalance with the same ease as an EMT simulator, however, many TS-EMT simulators use TS solvers that are single-phase, positive sequence-based, and this deteriorates the accuracy of the solution when there are heavy unbalances at the interface buses, as reported in [32], for example. 12 An SFA-EMT simulator based on MATE does not require the EMT time step to be an integer multiple of the SFA time step, unlike most of hybrid TS-EMT simulators and also most of the hybrid simulators that use dynamic phasors. Based on the previous considerations, the objective of this research was set to develop a simulation framework with the capability to perform transient stability type of analyses through the use of a hybrid simulation environment that employs SFA modeling for most of the system and EMT modeling for particular loads, components, or subsystems that require such a detailed, device-level type of modeling. It is expected that this hybrid framework will achieve higher accuracy than traditional hybrid TS-EMT simulators and a better computational efficiency than an EMT simulator. This macro objective was broken down into the following two steps: Objective 1: Develop a general interfacing scheme between the SFA and EMT simulation frameworks that allows their integration into an accurate and efficient hybrid multirate simulation framework. Objective 2: Develop a hybrid SFA-EMT simulator that operates with the derived interfacing scheme and allows the simultaneous simulation of the electromechanical transients in a power system and the electromagnetic transients in some parts of it. 13 Chapter 2: Shifted Frequency Analysis of Power Systems When studying the dynamic behavior of a power system due to its electromechanical transients, it is observed that after the high frequencies associated to the electromagnetic transients have decayed, the signals in the system, voltages and currents, have frequency spectra that are centered around the fundamental frequency of the system, with half bandwidths less than that frequency. This type of signals is called bandpass signals [62]. In time domain, these signals are 50 or 60 Hz sinusoids with a slow modulation in amplitude and phase or frequency. When the system is stable, as it returns to steady state the bandwidth around the system frequency shrinks more and more, eventually vanishing when the steady state is reached. In the simulation of the electromechanical transients to assess the stability of the system, the power engineer is interested in evaluating the modulations in amplitude and phase and the instantaneous frequency of the signals in the system, rather than the instantaneous values in time of those signals. Since the variations in amplitude, phase and frequency of those bandpass voltages and currents are much slower than the system frequency, from a computational efficiency standpoint, it is desirable to capture in simulation only those modulations of the signals, because that would enable the use of larger time integration steps. In fact, this is what conventional transient stability (TS) simulators do with their quasi-steady state approach, whereby differential equations are used for the electromechanical solution and steady state solutions are used for the network, through a standard AC phasor representation of it. 14 EMTP-type of simulators are very accurate to study the stability of a power system, however, from a computational efficiency point of view, the major difficulty of those simulators is that, because they are based on a three phase instantaneous waveform representation of the system variables, the EMT solution also tracks the AC carrier (i.e. the 50 or 60 Hz component), which limits the maximum time step that can be used to simulate the system. The slower a signal varies, the larger the integration step required to capture it in simulation; therefore, to enable the efficient use of an EMTP-type of simulator, the idea is to transform the original signal into an equivalent one that varies more slowly and solve for that signal. This is in essence what is done in shifted frequency analysis and in the following it is shown how. 2.1 Theoretical Foundations 2.1.1 The Dynamic Phasor A power system bandpass signal 𝑢ሺ𝑡ሻ with its frequency spectrum centered on the system frequency 𝜔଴ can be represented [63], [64] in terms of two low-pass signals, 𝑢ூሺ𝑡ሻ and 𝑢ொሺ𝑡ሻ modulating two in quadrature sinusoidal carrier signals (same frequency and 90° apart in phase) 𝑐𝑜𝑠𝜔଴𝑡 and 𝑠𝑖𝑛𝜔଴𝑡, as in (2.1): 𝑢ሺ𝑡ሻ ൌ 𝑢ூሺ𝑡ሻ𝑐𝑜𝑠𝜔଴𝑡 െ 𝑢ொሺ𝑡ሻ𝑠𝑖𝑛𝜔଴𝑡 (2.1) The low-pass signals 𝑢ூሺ𝑡ሻ and 𝑢ொሺ𝑡ሻ are referred to as the in-phase and in quadrature components of 𝑢ሺ𝑡ሻ, respectively. 15 The dynamic phasor 𝑈ഥሺ𝑡ሻ for the signal 𝑢ሺ𝑡ሻ is defined as: 𝑈ഥሺ𝑡ሻ ൌ 𝑢ூሺ𝑡ሻ ൅ 𝑗𝑢ொሺ𝑡ሻ (2.2) 𝑈ഥሺ𝑡ሻ, is also known as the complex envelope of 𝑢ሺ𝑡ሻ in signal processing and is as an analogous complex-valued low-pass representation of the original bandpass signal 𝑢ሺ𝑡ሻ. 2.1.2 The Analytic Signal Another representation of 𝑢ሺ𝑡ሻ is through its analytic signal 𝑧̅ሺ𝑡ሻ, defined as: 𝑧̅ሺ𝑡ሻ ൌ 𝑢ሺ𝑡ሻ ൅ 𝑗𝐻ሾ𝑢ሺ𝑡ሻሿ (2.3) where 𝐻ሾ⋅ሿ denotes the Hilbert transform [65] and is given by 𝐻ሾ𝑢ሺ𝑡ሻሿ ൌ 1𝜋𝑡 ∗ 𝑢ሺ𝑡ሻ ൌ 1𝜋න 𝑢ሺ𝜏ሻ𝑡 െ 𝜏 𝑑𝜏ஶିஶ (2.4) The Hilbert transform is referred to as the quadrature function of 𝑢ሺ𝑡ሻ and when added to it as the imaginary part of the analytic function, turns it into a rotating vector in the complex plane [66]. The analytic signal 𝑧̅ሺ𝑡ሻ is an analogous complex-valued bandpass representation of the original real-valued bandpass signal 𝑢ሺ𝑡ሻ. It has no negative frequency components and its Fourier transform is twice the positive frequency part of the Fourier transform of 𝑢ሺ𝑡ሻ [67], [66], [68]. Hence, this offers an alternative method to obtain the analytic signal by operating in the frequency domain on the Fourier transform of the signal 𝑢ሺ𝑡ሻ. 16 If 𝑢ሺ𝜔ሻ is the Fourier transform of 𝑢ሺ𝑡ሻ and 𝑠ሺ𝜔ሻ is a unit step function in the frequency domain, then 𝑧̅ሺ𝑡ሻ ൌ 𝐹ିଵሾ2𝑠ሺ𝜔ሻ𝑢ሺ𝜔ሻሿ (2.5) where 𝐹ିଵሾ⋅ሿ denotes inverse Fourier transform. Since 12 𝛿ሺ𝑡ሻ ൅ 𝑗2𝜋𝑡 ൌ 𝐹ିଵሾ𝑠ሺ𝜔ሻሿ (2.6) it follows that 𝑧̅ሺ𝑡ሻ ൌ 2 ൤12 𝛿ሺ𝑡ሻ ൅ 𝑗2𝜋𝑡൨ ∗ 𝑢ሺ𝑡ሻ ൌ 𝑢ሺ𝑡ሻ ൅ 𝑗 ൬ 1𝜋𝑡൰ ∗ 𝑢ሺ𝑡ሻ (2.7) Substituting the definition of the Hilbert transform, equation (2.4), in (2.7) renders 𝑧̅ሺ𝑡ሻ ൌ 𝑢ሺ𝑡ሻ ൅ 𝑗𝐻ሾ𝑢ሺ𝑡ሻሿ (2.8) Since equations (2.3) and (2.8) coincide, the above proves that the Fourier transform of the analytic signal contains only positive frequency component and therefore, it can either be obtained, a) operating in the frequency domain on the Fourier transform of the real signal by removing its negative frequency components, duplicating the result and then finding the inverse Fourier transform of it, or it can be obtained, b) operating in the time domain on the real signal by adding an imaginary part to it equal to its Hilbert transform. The importance of the analytic signal is that it retains the real signal 𝑢ሺ𝑡ሻ in its real part, its magnitude is the envelope of 𝑢ሺ𝑡ሻ and the derivative of its angle is the instantaneous frequency of 𝑢ሺ𝑡ሻ [68]. The relationship between the dynamic phasor 𝑈ഥሺ𝑡ሻ and the analytic signal 𝑧̅ሺ𝑡ሻ can be found as follows. Substituting (2.1) in (2.3) renders 𝑧̅ሺ𝑡ሻ ൌ 𝑢ூሺ𝑡ሻ𝑐𝑜𝑠𝜔଴𝑡 െ 𝑢ொሺ𝑡ሻ𝑠𝑖𝑛𝜔଴𝑡 ൅ 𝑗𝐻ൣ𝑢ூሺ𝑡ሻ𝑐𝑜𝑠𝜔଴𝑡 െ 𝑢ொሺ𝑡ሻ𝑠𝑖𝑛𝜔଴𝑡൧ ⇒ 17 𝑧̅ሺ𝑡ሻ ൌ 𝑢ூሺ𝑡ሻ𝑐𝑜𝑠𝜔଴𝑡 െ 𝑢ொሺ𝑡ሻ𝑠𝑖𝑛𝜔଴𝑡 ൅ 𝑗ሺ𝑢ூሺ𝑡ሻ𝑠𝑖𝑛𝜔଴𝑡 ൅ 𝑢ொሺ𝑡ሻ𝑐𝑜𝑠𝜔଴𝑡ሻ ⇒ 𝑧̅ሺ𝑡ሻ ൌ ቀ𝑢ூሺ𝑡ሻ ൅ 𝑗𝑢ொሺ𝑡ሻቁ 𝑒௝ఠబ௧ ൌ 𝑈ഥሺ𝑡ሻ𝑒௝ఠబ௧ (2.9) 2.1.3 The SFA Transformation Equation (2.9) shows that the dynamic phasor can be obtained from the analytic signal multiplying it by 𝑒ି௝ఠబ௧. This can be seen as the analytic signal being mapped into the dynamic phasor through a transformation defined as follows [41], [42]: 𝑇 ൌ 𝑒ି௝ఠబ௧ (2.10) 𝑇ିଵ ൌ 𝑒௝ఠబ௧ (2.11) The signal in the normal time domain would be related to the signal in the transformed or “modal domain” (denoted with subscript “m”) as follows: 𝑧௠̅ሺ𝑡ሻ ൌ 𝑈ഥሺ𝑡ሻ ൌ 𝑇𝑧̅ሺ𝑡ሻ (2.12) From Fourier analysis, it is known that multiplying a signal by 𝑒ି௝ఠబ௧, causes its Fourier transform to be shifted 𝜔଴ frequency units to the left. Thus, for an analytic bandpass signal with its frequency spectrum centered on 𝜔଴, the transformed signal is a low-pass signal with its frequency spectrum centered at zero frequency. This low-pass signal, the dynamic phasor, changes much less rapidly than the analytic signal and, therefore, larger time steps can be used to capture this signal in simulation, because the maximum frequency associated to the transformed signal is much lower. The above transformation is the one used for Shifted Frequency Analysis (SFA) and the transformed signals in the shifted frequency domain or SFA domain, are the corresponding dynamic phasors of the original real-valued time signals. 18 The whole SFA process to derive the dynamic phasor from the real bandpass signal is exemplified graphically in Figure 2.1. Figure 2.1 Amplitude spectra of a signal in different representations for the SFA process: a) Real bandpass signal, b) Analytic signal and c) SFA dynamic phasor. Since the dynamic phasor is calculated by multiplying the analytic signal by the SFA transformation and its magnitude is 1, it follows that the magnitude of the dynamic phasor is equal to the magnitude of the analytic signal and hence equal to the envelope of the real signal 𝑢ሺ𝑡ሻ. The situation is different regarding the phase of the dynamic phasor, because the transformation demodulates the analytic signal by removing the sinusoidal carrier. It can be easily shown that the derivative of the phase of the dynamic phasor is the instantaneous deviation of the instantaneous frequency of the real signal 𝑢ሺ𝑡ሻ from the fundamental electrical frequency of the system 𝜔଴. The SFA transformation is a complex transformation; therefore, the application of it and its inverse involve math operations in the complex plane. However, by separating the complex 19 variables into their real and imaginary parts, these operations can be done in the vector space 𝑅ଶ, using real arithmetic. In this case, the SFA transformation becomes the 2x2 matrix 𝑹 indicated below. 𝑹 ൌ ൤ 𝑐𝑜𝑠𝜔଴𝑡 𝑠𝑖𝑛𝜔଴𝑡െ𝑠𝑖𝑛𝜔଴𝑡 𝑐𝑜𝑠𝜔଴𝑡൨ (2.13) 𝑹ି𝟏 ൌ ൤𝑐𝑜𝑠𝜔଴𝑡 െ𝑠𝑖𝑛𝜔଴𝑡𝑠𝑖𝑛𝜔଴𝑡 𝑐𝑜𝑠𝜔଴𝑡 ൨ (2.14) In this research project and in this thesis all the operations with the SFA transformation have been done in the complex plane using complex arithmetic. Also, throughout this dissertation, unless otherwise stated, the term dynamic phasor (DP) will refer exclusively to SFA-type DP’s or DP’s in the SFA domain. 2.2 SFA Models of Basic Network Elements 2.2.1 Dynamic Phasors in an Electric Circuit In SFA the idea is to solve for the dynamic phasors in simulation, but before moving to that and proceeding to present the SFA models of basic circuit elements it is useful to point out some properties of the dynamic phasors, analytic signals and linear systems. Property 1: Let the response of a linear time-invariant system (LTI) to an input 𝑣ሺ𝑡ሻ, be 𝑢ሺ𝑡ሻ. If this system is now excited with a signal that is the Hilbert transform of 𝑣ሺ𝑡ሻ, then the response of the system to that signal will be the Hilbert transform of 𝑢ሺ𝑡ሻ. Also, because of linearity, if the system is 20 excited with an analytic signal the response will also be an analytic signal. This can be easily proved by noting that the response of an LTI system is the convolution of the input signal with the impulse response of the system [62] and by considering the definition of the Hilbert transform and the properties of the convolution integral. Property 2: Kirchhoff’s voltage and current laws, KVL and KCL, respectively, hold in a linear circuit excited by a complex signal and, in particular, by an analytic signal. Since the dynamic phasors are obtained by multiplying the analytic signals by the SFA transformation, it is also true that KVL and KCL hold for the dynamic phasors in the SFA domain. Based on the two properties above, it follows that if the forcing functions of a linear circuit are bandpass analytic signals; then, by transforming those forcing functions and the circuit to the SFA domain, the solution of the circuit in that domain will consist of low-pass dynamic phasors everywhere in the circuit. This solution can then be found by formulating and solving the differential equations that describe the dynamic behavior of the circuit in the SFA domain, which will have a different structure, or form, than the differential equations in normal time domain because now they have to reflect relationships between dynamic phasors. The solutions of the differential equations of power systems are usually found numerically by discretizing them and solving them in discrete-time domain. In the case of EMTP-like programs, this is done using the trapezoidal rule and also the backward Euler rule. By discretizing the differential equations of the circuit elements in the same fashion as in an EMT program, EMTP-21 like equivalent circuits can be derived for the different network elements and the solution in discrete-time can be found in SFA in a similar way as in an EMTP-type of program [40], [42] [43]. In the next two sections the derivation of these equivalent circuits is presented. In the first section the continuous-time differential equations in the SFA domain of the various circuit elements are determined and then, in the following section, the difference equations in discrete-time SFA domain for the various circuit elements are obtained by discretizing, with either the trapezoidal or the backward Euler rules, the continuous-time differential equations in the SFA domain. Detailed derivations of these equivalents for the single-phase components are shown in [43]. Here they are extended for the multiphase case. 2.2.2 Continuous-time SFA Equivalent Circuits of Basic R-L-C Elements Notation. - As a reminder on notation, in this dissertation complex quantities are represented with barred characters and real quantities with normal characters. Vectors and matrices are indicated with bold type characters and scalar quantities are represented with normal type characters. Dynamic phasors, or complex signals in the SFA domain, are represented with the subscript “m”. To derive the continuous-time equivalents for the basic circuit elements, it is useful to employ the property of the SFA transformation when applied to the derivative of a signal, being described as follows. 22 Let 𝒛തሺ𝑡ሻ be a vector of M complex signals. The application of the SFA transformation 𝑇 ൌ𝑒ି௝ఠబ௧, to the derivative of 𝒛തሺ𝑡ሻ renders: 𝑇ത ൤ 𝑑𝑑𝑡 𝒛ത ሺ𝑡ሻ൨ ൌ 𝑇ത ൤ 𝑑𝑑𝑡 𝒛തሺ𝑡ሻ𝑒ି௝ఠ೚௧𝑒௝ఠబ௧൨ ൌ 𝑇ത ൤ 𝑑𝑑𝑡 ൫𝒛ത௠ሺ𝑡ሻ𝑒௝ఠబ௧൯൨ ⇒ 𝑇ത ൤ 𝑑𝑑𝑡 𝒛ത ሺ𝑡ሻ൨ ൌ 𝑒ି௝ఠబ௧ ൤𝑒௝ఠబ௧ 𝑑𝑑𝑡 𝒛ത௠ሺ𝑡ሻ ൅ 𝑗𝜔଴𝑒௝ఠబ௧𝒛ത௠ሺ𝑡ሻ൨ ⇒ 𝑇ത ൤ 𝑑𝑑𝑡 𝒛ത ሺ𝑡ሻ൨ ൌ 𝑑𝑑𝑡 𝒛ത௠ሺ𝑡ሻ ൅ 𝑗𝜔଴𝒛ത௠ሺ𝑡ሻ (2.15) M-phase Coupled Resistances The normal time domain relationship between the complex vectors of currents and voltages in M coupled resistances is 𝒗ഥሺ𝑡ሻ ൌ 𝑹ଙሺ̅𝑡ሻ (2.16) Applying the SFA transformation (2.10) to (2.16) renders 𝒗ഥ௠ሺ𝑡ሻ ൌ 𝑹ଙ௠̅ሺ𝑡ሻ (2.17) Therefore, the equivalent circuit of M coupled resistances in the SFA domain is just the same as in normal time domain. M-phase Coupled Inductances The normal time domain relationship between the complex vectors of currents and voltages in M coupled inductances is 𝒗ഥሺ𝑡ሻ ൌ 𝑳 𝑑𝑑𝑡 ଙሺ̅𝑡ሻ (2.18) Applying the SFA transformation (2.10) to (2.18) and using property (2.15) renders 23 𝒗ഥ௠ሺ𝑡ሻ ൌ 𝑳 𝑑𝑑𝑡 ଙ௠̅ሺ𝑡ሻ ൅ 𝑗𝜔଴𝑳ଙ௠̅ሺ𝑡ሻ (2.19) Hence, the equivalent circuit of M coupled inductances in the SFA domain is the M coupled inductances in series with M coupled inductive reactances evaluated at 𝜔଴. Series Connection of M-phase Coupled Resistances and Inductances By combining (2.17) and (2.19) it follows that in the SFA domain the differential equation for M coupled resistances in series with M coupled inductances is given by 𝒗ഥ௠ሺ𝑡ሻ ൌ 𝑳 𝑑𝑑𝑡 ଙ௠̅ሺ𝑡ሻ ൅ ሺ𝑹 ൅ 𝑗𝜔଴𝑳ሻଙ௠̅ሺ𝑡ሻ (2.20) M-phase Coupled Capacitances The normal time domain relationship between the complex vectors of currents and voltages in M coupled capacitances is ଙሺ̅𝑡ሻ ൌ 𝑪 𝑑𝑑𝑡 𝒗ഥሺ𝑡ሻ (2.21) Applying the SFA transformation (2.10) to (2.20) and using property (2.15) renders ଙ௠̅ሺ𝑡ሻ ൌ 𝑪 𝑑𝑑𝑡 𝒗ഥ௠ሺ𝑡ሻ ൅ 𝑗𝜔଴𝑪𝒗ഥ௠ሺ𝑡ሻ (2.22) Therefore, the equivalent circuit of M coupled capacitances in the SFA domain is the M coupled capacitances in parallel with M coupled capacitive susceptances evaluated at 𝜔଴. 2.2.3 Discrete-time SFA Equivalents of M-phase Coupled Resistances For M coupled resistances the difference equation, with either the trapezoidal or backward Euler rules, has the same form as the continuous-time differential equation (2.17), that is 24 𝒗ഥ௠ሺ𝑡ሻ ൌ 𝑹ଙ௠̅ሺ𝑡ሻ (2.23) 2.2.4 Discrete-time SFA Equivalents of M-phase Coupled Inductances The discrete-time equivalent of M coupled inductances in the SFA domain is given by the following equations, obtained discretizing equation (2.19) with the trapezoidal and backward Euler rules. a) Trapezoidal ଙ௠̅ሺ𝑡ሻ ൌ 𝒀௅்௥𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅௅௠்௥ ሺ𝑡ሻ (2.24) 𝒀௅்௥ ൌ ൤ 2Δ𝑡 𝑳 ൅ 𝑗𝜔଴𝑳൨ିଵ ൌ 12Δ𝑡 ൅ 𝑗𝜔଴ 𝑳ିଵ (2.25) ଙ௛̅௅௠்௥ ሺ𝑡ሻ ൌ 𝒀௅்௥𝒗ഥ௠ሺ𝑡 െ Δ𝑡ሻ ൅2Δ𝑡 െ 𝑗𝜔଴2Δ𝑡 ൅ 𝑗𝜔଴ଙ௠̅ሺ𝑡 െ Δ𝑡ሻ (2.26) b) Backward Euler ଙ௠̅ሺ𝑡ሻ ൌ 𝒀௅஻ா𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅௅௠஻ா ሺ𝑡ሻ (2.27) 𝒀௅஻ா ൌ ൤ 1Δ𝑡 𝑳 ൅ 𝑗𝜔଴𝑳൨ିଵ ൌ 11Δ𝑡 ൅ 𝑗𝜔଴ 𝑳ିଵ (2.28) ଙ௛̅௅௠஻ா ሺ𝑡ሻ ൌ 11 ൅ 𝑗𝜔଴Δ𝑡 ଙ௠̅ሺ𝑡 െ Δ𝑡ሻ (2.29) 2.2.5 Discrete-time SFA Equivalents of M-phase Coupled R-L Branches The discrete-time equivalent in the SFA domain of M coupled resistances in series with M coupled inductances is given by the following equations, obtained discretizing equation (2.20) with the trapezoidal and backward Euler rules. 25 a) Trapezoidal ଙ௠̅ሺ𝑡ሻ ൌ 𝒀ோ௅்௥𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅ோ௅௠்௥ ሺ𝑡ሻ (2.30) 𝒀ோ௅்௥ ൌ ൤𝑹 ൅ 2Δ𝑡 𝑳 ൅ 𝑗𝜔଴𝑳൨ିଵ (2.31) ଙ௛̅ோ௅௠்௥ ሺ𝑡ሻ ൌ 𝒀ோ௅்௥𝒗ഥ௠ሺ𝑡 െ Δ𝑡ሻ ൅ 𝒀ோ௅்௥ ൬െ𝑹 ൅ 2Δ𝑡 𝑳 െ 𝑗𝜔଴𝑳൰ ଙ௠̅ሺ𝑡 െ Δ𝑡ሻ (2.32) b) Backward Euler ଙ௠̅ሺ𝑡ሻ ൌ 𝒀ோ௅஻ா𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅ோ௅௠஻ா ሺ𝑡ሻ (2.33) 𝒀ோ௅஻ா ൌ ൤𝑹 ൅ 1Δ𝑡 𝑳 ൅ 𝑗𝜔଴𝑳൨ିଵ (2.34) ଙ௛̅௅௠஻ா ሺ𝑡ሻ ൌ 𝒀ோ௅஻ா 1Δ𝑡 𝑳ଙ௠̅ሺ𝑡 െ Δ𝑡ሻ (2.35) 2.2.6 Discrete-time SFA Equivalents of M-phase Coupled Capacitances The discrete-time equivalent of M coupled capacitances in the SFA domain is given by the following equations, obtained discretizing equation (2.22) with the trapezoidal and backward Euler rules. a) Trapezoidal ଙ௠̅ሺ𝑡ሻ ൌ 𝒀஼்௥𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅஼௠்௥ ሺ𝑡ሻ (2.36) 𝒀஼்௥ ൌ 2Δ𝑡 𝑪 ൅ 𝑗𝜔଴𝑪 ൌ ൬ 2Δ𝑡 ൅ 𝑗𝜔଴൰ 𝑪 (2.37) ଙ௛̅஼௠்௥ ሺ𝑡ሻ ൌ െଙ௠̅ሺ𝑡 െ Δ𝑡ሻ െ 𝒀஼்௥∗𝒗ഥ௠ሺ𝑡 െ Δ𝑡ሻ (2.38) b) Backward Euler ଙ௠̅ሺ𝑡ሻ ൌ 𝒀஼஻ா𝒗ഥ௠ሺ𝑡ሻ ൅ ଙ௛̅஼௠஻ா ሺ𝑡ሻ (2.39) 26 𝒀஼஻ா ൌ 1Δ𝑡 𝑪 ൅ 𝑗𝜔଴𝑪 ൌ ൬ 1Δ𝑡 ൅ 𝑗𝜔଴൰𝑪 (2.40) ଙ௛̅஼௠஻ா ሺ𝑡ሻ ൌ െ 1Δ𝑡 𝑪𝒗ഥ௠ሺ𝑡 െ Δ𝑡ሻ (2.41) 2.2.7 Structure of the R-L-C Discrete-time SFA Models By inspecting the difference equations for the various R-L-C elements in SFA presented in the previous sub-sections, it can be observed that they can be used to represent the respective elements by an equivalent circuit in Norton form composed of a complex admittance matrix in parallel with a vector of complex history current sources. The structure of these equivalent circuits coincides with the structure of the equivalent circuits for the corresponding elements in an EMT-type of solution, except that for the latter the equivalents contain only real quantities. Thus, the SFA solution algorithm is in essence the same as for an EMT solution, except that the calculations have to be done using complex arithmetic. It is interesting to note that all the difference equations presented above for the R-L-C elements reduce to the classical AC phasor analysis relationships when the system reaches steady state, because the dynamic phasors become static AC phasors with constant amplitude and phase. This can be quickly confirmed for the case of the coupled capacitances, by substituting: ଙ௠̅ሺ𝑡ሻ ൌଙ௠̅ሺ𝑡 െ Δ𝑡ሻ and 𝒗ഥ௠ሺ𝑡ሻ ൌ 𝒗௠ሺ𝑡 െ Δ𝑡ሻ in the history equations (2.38) and (2.41), for trapezoidal and backward Euler, respectively, and checking that equations (2.36) and (2.39) both convert to: ଙ௠̅ሺ𝑡ሻ ൌ 𝑗𝜔଴𝑪𝒗௠ሺ𝑡ሻ (for the coupled L’s and the R-L’s elements the verification has to be done using the Thévenin form of the difference equations, which are not shown here, as only the Norton form of these equations is presented). Thus, in steady state, SFA is a discrete-time 27 solution of classical AC phasor analysis; the latter therefore, can be seen as a particular case of SFA after all transients have died out. 2.2.8 SFA models of Other Circuit Elements The models in the SFA domain of other network elements such as transformers, transmission lines, sources, generators and loads are discussed in Chapter 5, where the simulators employed in this research are presented. 2.3 The SFA Solution 2.3.1 General Procedure As explained earlier in the chapter, in SFA the idea is to solve for the dynamic phasors in the shifted frequency domain and this is done using the same discrete-time solution method as the one employed in EMTP-type of programs using a suitable numerical integration rule and employing a proper time integration interval. To analyze a power system with SFA, the network is transformed to the shifted frequency domain by replacing each circuit element with its discrete-time SFA model and by deriving the analytic signals of the external real forcing functions of the system and transforming them to the shifted domain by applying the transformation 𝑇 ൌ 𝑒ି௝ఠబ௧ to them. The next step is to formulate the system equations in nodal form in the shifted domain by forming the nodal admittance matrix of the system and the vector of total injected currents, which would include both the contributions of the external sources and the history sources. By solving for the node voltages, 28 the solution of the system is found. This solution includes the time-varying magnitudes and angles of the dynamic phasors in the system. The magnitudes trace the envelopes of the real signals in the system and the angles provide information regarding the instantaneous frequency deviations of those signals from the fundamental frequency, which are parameters the power system analyst needs to study the dynamics of the system under the influence of the electromechanical transients. In addition, the SFA solution can be transformed back into normal time domain to get the actual real-valued solution of the system. This is done by interpolating the signals in the shifted domain, applying the inverse transformation 𝑇ିଵ ൌ 𝑒௝ఠబ௧ to them and taking the real part of the complex-valued analytic signals obtained. Since the solution of the system in the SFA domain is composed of low-pass dynamic phasors, the discrete-time solution of the circuit can be found using much larger time steps than would be required if solving for the bandpass analytic signals or the real bandpass signals. For example, assuming the frequency spectrum of the real signals in the network lie within a band of ±5 Hz around the frequency of the system, 50 or 60 Hz, this implies that the maximum frequency of those signals would be 65 Hz. On the other hand, the spectrum of the dynamic phasors will also lie within a band of ±5 Hz, but around 0 Hz, thus the maximum frequency would be 5 Hz. Considering the limitations imposed by the Shannon - Nyquist theorem and allowing for an error of 3% with the trapezoidal rule [69], an EMT-type of simulation would require a time step of: 1.5 ms to capture the real signals, whereas to capture the dynamic phasors an SFA simulation would require a time step of 20 ms, which is more than 13 times larger! This shows how SFA can boost the computational efficiency of EMT simulation when used to study the dynamics 29 associated to the electromechanical transients in a power system, without compromising the accuracy of the solution. 2.3.2 SFA for Electromagnetic Transients So far in the presentation, the application suggested for SFA has been the simulation of the bandpass signals associated to the electromechanical transients in a power system. However, with SFA the differential equations of the network are solved and, in fact, in the derivations of the SFA models for the circuit elements in Section 2.3 no assumption was made regarding the type of transients to be analyzed or even the type of complex signals involved. The SFA transformation is just a mathematical transformation that maps the signals in normal time domain to signals in another domain, the frequency shifted domain, therefore the electromagnetic transients can also be studied in the shifted frequency domain. This has been suggested in [70]. Of course, this would not be computationally advantageous because, as a discrete-time solution, the SFA solution is governed by the Shannon - Nyquist Sampling theorem [71], [72], which establishes a theoretical bound for the maximum frequency that can be captured in simulation. That theoretical maximum is known as the Nyquist frequency (fNy) is given by 𝑓ே௬ ൌ ଵଶ୼௧. The SFA time step would have to be set to capture the highest frequency associated to the electromagnetic transients and hence would be of the same order of magnitude as the time step for an EMT solution. However, the point to be made here is that SFA can indeed be used to capture fast changing transients in a system by using a properly small time step and this is something that, under certain conditions, may be convenient to do, as will become evident in Chapters 4 and 6. 30 Having said that, it is important to clarify what happens when an SFA simulation is set up to study the dynamics of the system under the influence of the electromechanical transients using a large time step, appropriate for low frequencies, and a disturbance occurs in the system that generates electromagnetic transients which have high frequencies associated to it. In this case, recalling again the Sampling theorem [71], [72], any frequencies of the electromagnetic transients above the Nyquist frequency would not appear in the discrete-time solution and, consequently, most of the electromagnetic transients would be “filtered out” by the SFA simulation, or said in another way, would not be present in the SFA solution. 2.3.3 Nonstandard SFA Solutions Elaborating further from the previous sub-section, nothing prevents the use of the SFA transformation on signals that are not analytic or even plain real signals. The SFA models for the circuit elements do not make any assumptions regarding the type of signals involved and the differential and difference equations reproduce accurately the behavior of the circuit elements in the transformed domain, provided that an appropriate time step is used to capture the frequencies of interest; thus, from a pure mathematical point of view there are no restrictions to the type of applications in which to use the SFA solution, only computational efficiency and pragmatic considerations would be involved. The application of the SFA transformation directly on a real signal does not reduce the maximum frequency contained in the transformed signal, indeed it increases it, because the Fourier transform of a real-valued signal 𝑢ሺ𝑡ሻ is Hermitian, which means it has the property that 𝑢ሺെ𝜔ሻ ൌ 𝑢∗ሺ𝜔ሻ and hence the frequency spectrum is double sided. Even if the signal is a 31 bandpass signal, in the shifted frequency domain it will still require a time step of the same order of magnitude as in normal time domain. The reason is that, though the SFA transformation would shift the band around ω0 to the origin, the band around -ω0 would be shifted to -2ω0! Nevertheless, SFA can be used that way, though with no computational advantage, unless it is done in combination with other actions as indicated in [41]. 2.4 Electric Power in Complex and Shifted Frequency Domains When analyzing power systems, electrical engineers are most of the time interested in knowing how the electrical power and the energy moves in the network, therefore, since this research is concerned with conducting power system analysis using analytic signals and solving the circuit in the shifted frequency domain, it is important to study the concept of power in both, normal time domain and in the shifted frequency domain. General One-Port Circuit Let us consider a simple, one-port circuit in normal time domain that absorbs a complex current 𝚤ሺ̅𝑡ሻ and has a complex voltage ?̅?ሺ𝑡ሻ across its terminals, as indicated in Figure 2.2. At this time, for the purpose of the discussion, these complex current and voltage need not be analytic. 32 Figure 2.2 One-port circuit to study the concept of power in complex domain, normal time and SFA domains Let us extend the concept of instantaneous power from real circuits and the concept of complex power from standard AC circuit analysis by defining the Total Instantaneous Complex Power (TICP) absorbed by a circuit excited by complex variables, ?̅?௧ሺ𝑡ሻ, as. ?̅?௧ሺ𝑡ሻ ൌ ?̅?ሺ𝑡ሻ𝚤∗̅ሺ𝑡ሻ (2.42) Now let us decompose ?̅?௧ሺ𝑡ሻ, ?̅?ሺ𝑡ሻ and 𝚤ሺ̅𝑡ሻ into their real and imaginary components. Let us then assume that ?̅?ሺ𝑡ሻ ൌ 𝑣௥ሺ𝑡ሻ ൅ 𝑗𝑣௜ሺ𝑡ሻ (2.43) 𝚤ሺ̅𝑡ሻ ൌ 𝑖௥ሺ𝑡ሻ ൅ 𝑗𝑖௜ሺ𝑡ሻ (2.44) ?̅?௧ሺ𝑡ሻ ൌ 𝑝௧ሺ𝑡ሻ ൅ 𝑗𝑞௧ሺ𝑡ሻ (2.45) Substituting (2.43) - (2.45) into (2.42) and dropping the dependence of time for convenience in the notation, renders ?̅?௧ ൌ 𝑝௧ ൅ 𝑗𝑞௧ ൌ 𝑣௥𝑖௥ ൅ 𝑣௜𝑖௜ ൅ 𝑗ሾ𝑣௜𝑖௥ െ 𝑣௥𝑖௜ሿ (2.46) 33 As it can be observed, 𝑝௧, the real part of ?̅?௧, has two components, one that is due to the interaction of the real parts of the voltage and the current and another that is due to the interaction of the corresponding imaginary parts, so this variable can be defined as the Total Instantaneous Active Power (TIAP) absorbed by the circuit and includes both, the instantaneous active power due to the real parts of ?̅? and 𝚤,̅ and the corresponding instantaneous active power due to their imaginary parts. Similarly, 𝑞௧, the imaginary part of ?̅?௧, has also two components, one that results from the interplay of the imaginary part of the voltage with the real part of the current and another that corresponds to the negative of the interplay between the real part of the voltage and the imaginary part of the current; thus, for the purpose of this research, this variable can be interpreted and defined as the Total Instantaneous Reactive Power (TIRP) absorbed by the circuit and includes the instantaneous reactive powers that result from the interplays of imaginary part with real part and vice versa of voltage and current, respectively. Now, let us evaluate the TICP in the shifted frequency domain. For this purpose, the circuit variables are transformed to the SFA domain using the equations below ?̅?௠ሺ𝑡ሻ ൌ ?̅?ሺ𝑡ሻ𝑇 ൌ ?̅?ሺ𝑡ሻ𝑒ି௝ఠబ௧ (2.47) 𝚤௠̅ሺ𝑡ሻ ൌ 𝚤ሺ̅𝑡ሻ𝑇 ൌ 𝚤ሺ̅𝑡ሻ𝑒ି௝ఠబ௧ (2.48) The TICP in the SFA domain ?̅?௧,௠ሺ𝑡ሻ is then ?̅?௧,௠ሺ𝑡ሻ ൌ ?̅?௠ሺ𝑡ሻ𝚤௠̅∗ ሺ𝑡ሻ ൌ ?̅?ሺ𝑡ሻ𝑒ି௝ఠబ௧൫𝚤ሺ̅𝑡ሻ𝑒ି௝ఠబ௧൯∗ ൌ ?̅?ሺ𝑡ሻ𝚤∗̅ሺ𝑡ሻ ൌ ?̅?௧ሺ𝑡ሻ (2.49) 34 Thus, the TICP is the same in normal time domain and the SFA domain and, in consequence, we can say that the SFA transformation is total instantaneous complex power invariant. Circuit in AC Steady State Let us now evaluate the TICP when the circuit is in sinusoidal steady state. To do that, let us assume that the complex voltage and current in the circuit of Figure 2.2, in normal time domain and SFA domain, are given by ?̅?ሺ𝑡ሻ ൌ 𝑉௣ cosሺ𝜔଴𝑡 ൅ 𝜙ሻ ൅ 𝑗𝑉௣ sinሺ𝜔଴𝑡 ൅ 𝜙ሻ ൌ 𝑉௣𝑒௝ሺఠబ௧ାథሻ (2.50) ?̅?௠ሺ𝑡ሻ ൌ 𝑉௣𝑒௝థ (2.51) 𝚤ሺ̅𝑡ሻ ൌ 𝐼௣ cosሺ𝜔଴𝑡 ൅ 𝜃ሻ ൅ 𝑗𝐼௣ sinሺ𝜔଴𝑡 ൅ 𝜃ሻ ൌ 𝐼௣𝑒௝ሺఠబ௧ାఏሻ (2.52) 𝚤௠̅ሺ𝑡ሻ ൌ 𝐼௣𝑒௝ఏ (2.53) It is easier to evaluate the TICP in the SFA domain. Using (2.49), (2.51) and (2.53) we have ?̅?௧ሺ𝑡ሻ ൌ ?̅?௧,௠ሺ𝑡ሻ ൌ ?̅?௠ሺ𝑡ሻ𝚤௠̅∗ ሺ𝑡ሻ ൌ 𝑉௣𝑒௝థ𝐼௣𝑒ି௝ఏ ൌ 𝑉௣𝐼௣𝑒௝ሺథିఏሻ (2.54) Hence 𝑝௧ሺ𝑡ሻ ൌ 𝑉௣𝐼௣cos ሺ𝜙 െ 𝜃ሻ (2.55) 𝑞௧ሺ𝑡ሻ ൌ 𝑉௣𝐼௣sin ሺ𝜙 െ 𝜃ሻ (2.56) If the expression for ?̅?௧ is evaluated in normal time domain, it is interesting to note the expressions for 𝑝௧ and 𝑞௧, prior to obtaining the compact expressions in (2.55) and (2.56), are. 𝑝௧ሺ𝑡ሻ ൌ 12𝑉௣𝐼௣ሾcosሺ2𝜔଴𝑡 ൅ 𝜙 ൅ 𝜃ሻ ൅ cosሺ𝜙 െ 𝜃ሻ െ cosሺ2𝜔଴𝑡 ൅ 𝜙 ൅ 𝜃ሻ ൅ cosሺ𝜙 െ 𝜃ሻሿ 35 𝑞௧ሺ𝑡ሻ ൌ 12𝑉௣𝐼௣ሾsinሺ2𝜔଴𝑡 ൅ 𝜙 ൅ 𝜃ሻ ൅ sinሺ𝜙 െ 𝜃ሻ െ sinሺ2𝜔଴𝑡 ൅ 𝜙 ൅ 𝜃ሻ ൅ sinሺ𝜙 െ 𝜃ሻሿ From the above equations it can be readily observed that the two components of 𝑝௧, one due to the real parts of voltage and current and the other due to the imaginary parts of voltage and current, they both have a double frequency oscillating term, as expected; however, those terms cancel each other out. A similar thing happens to the two components of 𝑞௧. As it can be seen from equations (2.54) – (2.56), in a circuit in sinusoidal steady state, ?̅?௧, 𝑝௧ and 𝑞௧ are all constant! Recalling the expressions for the complex power S̅, the active power P and the reactive power Q from standard AC circuit analysis (equations (2.57) – (2.59) below), it is interesting to note that these three parameters are equal, respectively, to half of the corresponding total instantaneous powers in the complex domain, either in normal time domain or in the SFA domain. 𝑃 ൌ 𝑉௥௠௦𝐼௥௠௦ cosሺ𝜙 െ 𝜃ሻ ൌ 𝑉௣𝐼௣2 cosሺ𝜙 െ 𝜃ሻ (2.57) 𝑄 ൌ 𝑉௥௠௦𝐼௥௠௦ sinሺ𝜙 െ 𝜃ሻ ൌ 𝑉௣𝐼௣2 sinሺ𝜙 െ 𝜃ሻ (2.58) 𝑆̅ ൌ 𝑃 ൅ 𝑗𝑄 ൌ 𝑉௣𝐼௣2 𝑒௝ሺథିఏሻ (2.59) The expression for 𝑞௧ in equation (2.56) offers an alternative perspective to see the reactive power Q in AC circuits as half of the total instantaneous reactive power in an SFA dynamic phasor domain. This topic offers an interesting opening for further research that seems worth exploring but, at this moment, it is considered beyond the scope of this research project. 36 Circuit with Dynamic Phasors Let us assume that the circuit in Figure 2.2 has been subjected to a disturbance, but that the associated electromagnetic transients have already died out and there remain only the slow electromechanical transients. Under these operating conditions, the variables in the circuit have become bandpass signals with their frequency spectra centered around the fundamental frequency ω0 and can, therefore, be represented by analytic signals in normal time domain and by dynamic phasors in the SFA domain. Let us further assume that the voltage across the terminals of the circuit and current absorbed by it can be described by the following equations: ?̅?ሺ𝑡ሻ ൌ 𝑉௣ሺ𝑡ሻ cos൫𝜔଴𝑡 ൅ 𝜙ሺ𝑡ሻ൯ ൅ 𝑗𝑉௣ሺ𝑡ሻ sin൫𝜔଴𝑡 ൅ 𝜙ሺ𝑡ሻ൯ ൌ 𝑉௣ሺ𝑡ሻ𝑒௝൫ఠబ௧ାథሺ௧ሻ൯ (2.60) ?̅?௠ሺ𝑡ሻ ൌ 𝑉௣ሺ𝑡ሻ𝑒௝థሺ௧ሻ (2.61) 𝚤ሺ̅𝑡ሻ ൌ 𝐼௣ሺ𝑡ሻ cos൫𝜔଴𝑡 ൅ 𝜃ሺ𝑡ሻ൯ ൅ 𝑗𝐼௣ሺ𝑡ሻ sin൫𝜔଴𝑡 ൅ 𝜃ሺ𝑡ሻ൯ ൌ 𝐼௣ሺ𝑡ሻ𝑒௝൫ఠబ௧ାఏሺ௧ሻ൯ (2.62) 𝚤௠̅ሺ𝑡ሻ ൌ 𝐼௣ሺ𝑡ሻ𝑒௝ఏሺ௧ሻ (2.63) As long as the voltage and the current can be described by the above equations, the derivations previously made for the total instantaneous powers of the circuit in AC sinusoidal steady state also apply for this case; thus, under these conditions ?̅?௧, 𝑝௧ and 𝑞௧ are given by the following expressions: ?̅?௧ሺ𝑡ሻ ൌ ?̅?௧,௠ሺ𝑡ሻ ൌ 𝑉௣ሺ𝑡ሻ𝐼௣ሺ𝑡ሻ𝑒௝൫థሺ௧ሻିఏሺ௧ሻ൯ (2.64) 𝑝௧ሺ𝑡ሻ ൌ 𝑉௣ሺ𝑡ሻ𝐼௣ሺ𝑡ሻcos ሺ𝜙ሺ𝑡ሻ െ 𝜃ሺ𝑡ሻሻ (2.65) 37 𝑞௧ሺ𝑡ሻ ൌ 𝑉௣ሺ𝑡ሻ𝐼௣ሺ𝑡ሻsin ሺ𝜙ሺ𝑡ሻ െ 𝜃ሺ𝑡ሻሻ (2.66) Now ?̅?௧, 𝑝௧ and 𝑞௧ are no longer constant, but vary slowly according to the change in time of the amplitudes and phases of the voltage and current. Half of these parameters would also describe the evolution in time of the complex, active and reactive powers, S̅, P and Q, respectively, normally used in standard AC circuit analysis 38 Chapter 3: Simulation of Power Systems in a Partitioned Environment When simulating a system as a whole, to study its dynamic behavior, the numerical solution of its differential equations has to be found using a single integration interval or step size Δt, which can be variable. However, there are several reasons why, under certain circumstances, it would be convenient or required to partition a system into subsystems and solve it as a set of interconnected subsystems. For instance, it may occur that for a certain system some parts of it exhibit much faster responses while other parts change more slowly. This type of systems is said to have latency; so, simulating such a system as a single unit may imply a waste of computational resources. This may be irrelevant for small systems but can be significant for large scale systems. For large systems of this type it may then be convenient to improve the computational efficiency, by partitioning the system into subsystems and using different step sizes to solve each subsystem with a dual rate or with a multirate type of solution. Not all multirate solutions aim at exploiting latency, there are cases in which the type of study and the modeling capabilities and nature of the simulation programs involved, require the partitioning of the system to solve some parts of it with one type of simulator and other parts with another type of simulator. This is the case of hybrid TS-EMT simulators for the study of the transient stability of a system and it is also the case for this research, in which we want to interface an SFA simulator with an EMT simulator. These simulators have similar structure and solution methods but different domains of solution, the SFA works in a complex plane and the EMT works in a real plane. In other cases, the partitioning of the system may be convenient to employ parallel computing when solving it, again for reasons of computational efficiency. Still 39 in other instances, known as hardware-in-loop applications, the simulation may involve the simultaneous operation of software and hardware, like in the case of testing actual protection devices in a real-time simulator. In these cases, the system object of the simulation is already physically partitioned; therefore, the simulation has to run on an actually partitioned system. Depending on the type of system, the application study, the simulators involved and the way in which the splitting of the system is accomplished, the multirate simulation of a partitioned system requires the implementation of a specific solution protocol to define the way in which the solutions of the different subsystems interface and it may involve approximations in the handling of the subsystems solutions or the inclusion of real or artificial interfacing devices. For example, some schemes may require iterations, interpolation, prediction or extrapolation of variables and some others may also require the use of transmission lines to decouple the solutions through the travel time in the lines. In this research the Multi-Area Thévenin Equivalent (MATE) solution framework [73] has been used as it provides a means to find an exact, direct and simultaneous solution of a partitioned system in a close-form. In this chapter this framework is presented and various concepts associated to the multirate simulation of power systems are explained. 3.1 Multi-Area Thevenin Equivalent Conceptual Framework A detailed presentation of MATE and its formulation in EMT solutions is available in [73], [74] and [75]. Here, a summary of the procedure is presented to establish a frame of reference for the discussions that follow. The MATE solution framework provides an effective means for 40 partitioning a large system into subsystems connected by link branches. The subsystems are solved independently and the whole system solution is found by solving for the links’ interactions. In the general formulation of MATE, the link branches can be any type of element but, for the application in this research, they are considered as resistances with no mutual coupling. The number of subsystems in MATE can be N in general, however, to explain the concept, the system in Figure 3.1 composed of two subsystems will be considered. Furthermore, although the concepts in MATE apply to a general circuit solution, since this research is concerned with EMT and SFA solutions which share the same solution algorithm, in the following explanation a discrete-time EMT solution is used. Figure 3.1 Sample system for MATE explanation The system in Figure 3.1 can be described by a set of modified nodal equations [74], as follows ቎𝒀ଵ 𝟎 𝑪ଵ𝟎 𝒀ଶ 𝑪ଶ𝑪ଵ் 𝑪ଶ் െ𝑹𝑳቏ ൥𝑽ଵ𝑽ଶ𝑰𝑳൩ ൌ ൥𝒉ଵ𝒉ଶ𝟎൩ (3.1) 41 In MATE, Yj, Vj and hj are, respectively for subsystem j, the admittance matrix, the vector of nodal voltages and the vector of injected currents. The latter includes the history sources in the EMT model, and the actual sources, but not the currents from the links. Matrix Cj is the incidence or connection matrix that registers the connectivity of the nodes of subsystem j to the links. IL is the vector of link currents and RL is a diagonal matrix with the link resistances. Operating on (3.1) to find IL renders ቎𝑰 𝟎 𝒀ଵି ଵ𝑪ଵ𝟎 𝑰 𝒀ଶିଵ𝑪ଶ𝟎 𝟎 𝒁𝑻቏ ൥𝑽ଵ𝑽ଶ𝑰𝑳൩ ൌ ൥𝒆ଵ𝒆ଶ𝑬𝑳൩ (3.2) where I is the identity matrix and, 𝒆௝ ൌ 𝒀௝ି ଵ𝒉௝ (3.3) 𝒆்௛,௝ ൌ 𝑪௝் 𝒆௝ (3.4) 𝑬௅ ൌ 𝒆்௛,ଵ ൅ 𝒆்௛,ଶ (3.5) 𝒁்௛,௝ ൌ 𝑪௝் 𝒀௝ି ଵ𝑪௝ (3.6) 𝒁் ൌ 𝒁்௛,ଵ ൅ 𝒁்௛,ଶ ൅ 𝑹௅ (3.7) From the inspection of (3.3) it can be seen that ej is the vector of node voltages of subsystem j when solved independently of the rest of the subsystems, that is, with all the links open-circuited. eTh,j and ZTh,j are, respectively, the vector of Thévenin voltages and the matrix of Thévenin impedances of subsystem j as seen from the links. Equation (3.4) shows that the connection matrix picks from the vector of open-circuit voltages, the voltages of the nodes connected to the links to render the vector of Thévenin voltages of the respective subsystem. A similar task is performed by the connection matrix in (3.6) to produce the matrix of Thévenin impedances. EL is 42 the sum of the Thévenin voltages of the subsystems and ZT is the sum of the Thévenin impedances of the subsystems plus the resistances of the links. In the general case with more than two subsystems, EL and ZT will include, respectively, the sum of the Thévenin voltages and the Thévenin impedances of all subsystems. The set of equations for IL in (3.2) states that the currents through the links can be calculated by substituting all subsystems with their respective Thévenin equivalents and applying the following equation 𝑰௅ ൌ 𝒁்ିଵ𝑬௅ (3.8) After finding IL the voltages in all subsystems can be found from the corresponding equations obtained from (3.2), which for subsystem j is 𝑽௝ ൌ 𝒆௝ െ 𝒀௝ି ଵ𝑪௝𝑰௅ (3.9) This equation shows how the solution voltages for subsystem j are found by correcting or updating the open circuit solution with the injection of the link currents into that subsystem. Again, the connection matrix selects which nodes of the subsystem will receive link currents. From the equations above a general solution algorithm of MATE can be derived [73]. 3.2 MATE Solution Algorithm For a solution in single rate mode, that is, using the same step size Δt for all subsystems, the standard MATE solution procedure to solve a system at a given time t, is the following: i. Solve for the link currents applying equation (3.8) and using the Thévenin equivalents of all subsystems calculated at the end of the previous solution time. 43 ii. Solve for the node voltages in all subsystems by injecting the link currents, that is, by applying equation (3.9) for all subsystems. This will update the open circuit solutions of all subsystems and generate the solution of the whole system with all subsystems connected with the links. iii. For each subsystem evaluate the current injections hj that would be required for the next solution time at t = t + t. iv. Evaluate the Thévenin equivalents of all subsystems for next solution time at t = t + t, by solving the subsystems in open circuit (with the links disconnected) using the current injections hj evaluated in step (iii) and applying equations (3.3) and (3.4). If there are topology changes, the new Thévenin impedances are found using equation (3.6). v. Advance the solution time by Δt, to t = t + Δt and go back to step (i) until t = tmax. Notice that as long as the same integration step is used for all subsystems, the MATE solution of the partitioned system is exact, that is, it is the same as solving the entire system with the same step size. This is because there are no approximations in the derivation of MATE. The Thévenin equivalents used in the determination of the link currents with all subsystems connected, are an exact representation of each subsystem. Therefore, even though the system is solved in a partitioned way, it is as if the entire system was being solved as a single unit. An important feature of the above MATE solution method is that the Thévenin equivalents of each subsystem, that will be used to find the solution of the system at the present solution time, were evaluated right after the solution of the subsystems were found at the previous time step 44 (check steps (iii), (iv) and (i) above). This means that at every solution time, for each subsystem there is a Thévenin equivalent available which was calculated at the previous time instant. 3.3 Use of MATE in EMT-EMT Multirate Simulation In a multirate simulation of a power system, the subsystems that comprise it are solved at different rates and their solutions are jointly interfaced for a whole system solution. The MATE procedure as presented in the previous section applies for an EMT-EMT simulator operating in single rate mode, but it can be easily modified to use it in an EMT-EMT simulator operating in multirate mode. In [60], [61], [76] and [77] successful MATE multirate EMT-EMT simulations are presented. The applications in [60], [76] and [77] required the time steps to be multiples, but in [61] that restriction was removed; that is, the time steps can be arbitrary. This latter approach is the one presented here and the one used in this research. In MATE, to solve for the link currents needed to update the solution of a subsystem, a Thévenin equivalent for each of the subsystems must be available. If the system is being solved in multirate mode with no restrictions on the integration steps used to solve each subsystem; then, in general, each subsystem will require to be solved or updated at a different rate and hence at specific times and, at those specific times, there must be available Thévenin equivalents of the other subsystems which are running at different rates. However, this does not represent a problem to MATE since in an EMT solution the Thévenin required for the solution at the next time step of a particular subsystem is calculated at the present time step solution of that subsystem; therefore, by linear interpolation, a Thévenin equivalent of that subsystem can be made available for the updating of the other subsystems that need updating before that subsystem 45 and so on for the other subsystems. This process of interpolating the Thévenins voltages in presented graphically in Figure 3.2 for a case of two subsystems (SSs), slow (using a large Δt) and fast (using a small Δt). Figure 3.2 Interpolation of Thévenins of the slow subsystem to calculate the link solutions required to update the fast subsystem at various solution times. By updating the subsystems in a serial manner, beginning with the one with the fastest rate (lowest time step), followed by the second fast and ending with the one with the lowest rate (largest time step) the whole simulation can be accomplished. Every time a subsystem begins updating, it keeps updating until it is time for another subsystem to begin updating. Each subsystem will update when its turn comes and subsystems running at the same rate will update at the same time, hence their solutions can be parallelized for more efficiency in the computation. This is the general idea of the MATE protocol in multirate mode and, as it can be seen, the time steps need not be multiples, they can indeed be arbitrary. However, some specifics have to be considered to ensure the accuracy of the solution, as indicated below. 46 In this multirate solution method when a subsystem is being updated, the solutions of the subsystems advancing at a slower pace (larger time step) have to be, using the terminology of signal processing, oversampled (interpolated) and the solutions of subsystems running at a higher rate (smaller time step) have to be decimated [61]. The interpolation of slower solutions is straightforward and the interpolated values of the Thévenin voltages do not need additional conditioning or processing, they are ready to use. Note that these interpolations will not occur at uniformly distributed intervals, instead they will happen at any time as required by the updating subsystem, so in this sense, these interpolations differ from the standard oversampling process presented in signal processing textbooks. In the decimation of solutions advancing at faster rates, the plain interpolation of the Thévenin voltages will not be enough. In fact, it may actually degrade the accuracy of the multirate solution, because the fast solution may be ill-conditioned by the high frequencies that it may contain, which could not be accurately captured by the slower subsystem solution. In the decimation process the Shannon - Nyquist Sampling theorem [64], [71] gets into play. In plain terms, this theorem states that the maximum observable frequency in a sampled signal (theoretical), known as the Nyquist frequency (fNy) is given by the inverse of twice the sampling interval (Δt). Any frequency higher than the Nyquist frequency present in the signal being sampled will produce a phenomenon known as aliasing in the sampled signal, which will distort the sampled signal in comparison with the original signal object of the sampling. The implications of this for our case are as follows. 47 Let us assume that the system is partitioned into two subsystems, fast (using a small Δt) and slow (using a large Δt). Every time the slow SS is updated using the Thévenin equivalents from the fast SS, it is actually sampling (down sampling) those Thévenin voltages of the fast SS with a sampling interval equal to its own large time step (ΔtSlow); therefore, any frequencies contained in the Thévenin voltages coming from the fast SS that are higher than the Nyquist frequency of the slow SS, defined by its own time step (ΔtSlow), will produce aliasing in the solution of the slow subsystem, which will deteriorate it and in turn the whole multirate solution. This situation is explained graphically in Figure 3.3. Figure 3.3 Amplitude spectrum of the Fourier transform of the Thévenin of the fast subsystem, with indication of the range of frequencies that would produce aliasing in the slow subsystem solution. To prevent aliasing, prior to down sampling the Thévenin voltages coming from the fast subsystem, those Thévenin voltages have to be low-pass filtered (conditioned) to remove any frequencies above the Nyquist frequency of the slow SS. Thus, the decimation process includes a first stage of antialiasing processing (low-pass filtering) and after that, the smoothed out or low-48 pass filtered Thévenins are down sampled. These down sampled Thévenins will be the ones used in the updating of the slow subsystem. Additional discussions on antialiasing processing will be provided in the next chapter. 3.4 EMT-EMT MATE Multirate Interfacing Protocol Based on the previous considerations, for an EMT-EMT multirate solution of a system split into two subsystems (SSs), slow or SS1 and fast or SS2, using a large and a small time step, respectively, the MATE multirate interfacing protocol is as follows [61]. At the beginning of the simulation, the SSs are solved with the links open and their Thévenin equivalents (TEs) for the first time step solution of each SS are obtained. Then, the solution of the fast SS for the next small time steps begins. Using an interpolated value of the TE of SS1 and the TE of the SS2, the link currents are evaluated using equation (3.8) of the standard MATE solution. These currents are then injected into the fast SS to update it using equation (3.9). A new TE of the fast SS is evaluated by solving it in open circuit and using equations (3.3) and (3.4). Time is then advanced one small time step and the solution of the fast SS continues using the previous process, until it is time to update or solve the slow SS. At this point, SS2 has to provide its TE for the updating of the SS1. Since the slow SS will be receiving a down-sampled version of the TE of the fast SS, as discussed previously, in order to avoid aliasing errors, the TEs of the fast SS are first low-pass filtered or smoothed out before down sampling them [61]. Using the TE of SS1 and the low-pass filtered TE of SS2, the link currents are evaluated using equation (3.8). These currents are then injected into the slow SS to update it applying equation (3.9). A new TE of the slow SS for the solution at the next large time step is evaluated by solving it with 49 the links open and using equations (3.3) and (3.4). After this, the solution of SS2 resumes and the whole process is repeated until the maximum simulation time is reached. This procedure is shown graphically in Figure 3.4. Figure 3.4 MATE Interfacing protocol for an EMT-EMT multirate solution 50 Chapter 4: SFA-EMT Multirate Simulation In this chapter the main contribution of this research is presented, which is the implementation of an SFA-EMT hybrid simulation framework. In Section 4.1 the approach to accomplishing the hybrid SFA-EMT simulation is presented. SFA-EMT simulation can be done in both, single rate or multirate mode, so both simulation modes are presented, the former in Section 4.2 and the latter in Section 4.3. In search for optimizing the accuracy of the multirate SFA-EMT solution, various antialiasing processing schemes were studied in this project. One of the implemented antialiasing schemes is based on energy balance considerations and makes use of the original contribution of this research with regards to the concept of total instantaneous complex power in circuits fed by complex signals, as presented in Chapter 2, Section 2.4. The antialiasing methods are presented in a sub-section of Section 4.3. Section 4.4 presents a discussion about the selection of the boundaries between subsystems. 4.1 Interfacing SFA and EMT Simulators in MATE The SFA transformation allows the use of much larger time steps than would otherwise be required in a normal EMT simulation, but it turns the circuit variables into complex-valued ones. In a system partitioned into an SFA SS (slow SS) and an EMT SS (fast SS), intuitively, it would seem that the real-valued EMT solution should be interfaced with the real part of the SFA SS solution. However, this cannot be done directly, because the Thévenin equivalent of the SFA SS is complex-valued, both the voltages and the impedances. This causes the real and imaginary parts in normal time domain to be coupled at the Thévenin equivalent and this in turn, requires that the combined MATE solution at the links level to update the subsystems be carried out in 51 terms of complex numbers. This implies that for every SFA and EMT solution in the SFA-EMT multirate simulation, a complex Thévenin of the EMT SS is required and this means that for every single time step of the EMT solution, from the very first one to the very last, including any transient periods in between, in addition to the real EMT Thévenin, an accurate estimation of an imaginary EMT Thévenin is also necessary. Hence, to keep the relationships for the link solutions in terms of complex variables it is necessary to keep track of the response of the EMT SS in terms of both real and imaginary parts to combine with the SFA SS solution. This can be implemented efficiently and accurately by performing two simultaneous and parallel simulations of the EMT SS; one for the real part and one for the imaginary part. The solution for the real part is the actual solution of the system, while the solution for the imaginary part maintains the analyticity of the corresponding complex-plane solution of SFA and allows the direct interfacing of the real EMT solution with the complex-plane SFA solution, taking full advantage of the benefits of the SFA transformation of using large time steps and fully exploiting the strengths of MATE of interfacing the subsystems without a time step delay. This is a major contribution of this research as this concept has not been previously presented in the literature. It is important to mention that the option of numerically calculating the Hilbert transform of the real EMT Thévenins to obtain the imaginary EMT Thévenin at every solution time is not an accurate option. It can work in steady state, because the signals are sinusoidal, but not during the transient periods right after a disturbance. 52 The proposed method of interfacing the SFA and EMT solutions with MATE can be carried out in both, single or multirate modes. Although this research in mainly concerned with the SFA-EMT multirate simulation, before presenting the multirate mode of the protocol, the operation of the single rate mode will be explained because studying it facilitates the understanding of the multirate mode of the proposed protocol. 4.2 Single Rate SFA-EMT Simulation Let us consider a system partitioned into two subsystems connected with “n” links, one identified as Subsystem 1 (SS1) or SFA Subsystem (SFA SS) and the other identified as SS2 or EMT SS, as indicated in Figure 4.1. Applying the concept proposed in the previous section, the EMT SS actually comprises two separate subsystems, one for the real part of the EMT solution and a duplicate of it that runs in parallel for the imaginary part of the EMT solution. For ease of presentation, in the rest of this thesis the EMT SS will be treated as one single SS, with the understanding that whatever is done by that SS is done by both, the real and the imaginary parts of it. The SFA SS is to be solved with an SFA simulator and the EMT SS with an EMT simulator, both solutions use the same integration step Δt and the whole system solution will be found using MATE. 53 Figure 4.1 System split into two subsystems with “n” links for SFA-EMT simulation with MATE 4.2.1 Links Solution Let us assume that the Thévenin equivalents of both SSs for the solution at time t, have already been determined at the end of the previous time solution and now it is time to solve for the links at time t. For SS1 at time t, and with the links connected, we have that 𝒆ത்௛ଵ௠ሺ𝑡ሻ ൌ 𝒁ഥ்௛ଵଙ௅̅௠ሺ𝑡ሻ ൅ 𝒗ഥ்ଵ௠ሺ𝑡ሻ (4.1) where 𝒁ഥ்௛ଵ is the matrix of complex SFA Thévenin impedances of SS1 as seen from the links, and 𝒆ത்௛ଵ௠, 𝒗ഥ்ଵ௠ and ଙ௅̅௠, are the SFA vectors of the Thévenin sources and terminal voltages of SS1 and the currents through the links in the SFA domain, respectively (directions of link 54 currents is as in Figure 4.1). As there are n links, the dimension of these vectors is n and the dimension of the matrix is 𝑛𝑥𝑛. For updating the SSs, it is required to solve for the link currents in normal time domain. Thus, unshifting the SFA quantities in (4.1) by applying the inverse SFA transformation yields 𝒆ത்௛ଵሺ𝑡ሻ ൌ 𝒁ഥ்௛ଵଙ௅̅ሺ𝑡ሻ ൅ 𝒗ഥ்ଵሺ𝑡ሻ (4.2) where 𝒆ത்௛ଵ, 𝒗ഥ்ଵ and ଙ௅̅, are the vectors of Thévenin sources and terminal voltages of SS1 and the currents through the links, respectively, in normal time domain. Let us now consider both, the real and imaginary solutions of SS2 at time t with the links connected. Let 𝒆்௛ଶ,ோ, 𝒆்௛ଶ,ூ, 𝒗்ଶ,ோ, 𝒗்ଶ,ூ, 𝒊௅,ோ and 𝒊௅,ூ be, respectively, the vectors of real and imaginary Thévenin voltages, the voltages at the terminals of the EMT SS and the currents injected from the links into it. Let the matrices of Thévenin resistances of the EMT SS, as seen from the links, be 𝑹்௛ଶ,ோ and 𝑹்௛ଶ,ூ, for the real and imaginary solutions, respectively. Then for SS2, we have 𝒗்ଶ,ோሺ𝑡ሻ ൌ 𝒆்௛ଶ,ோሺ𝑡ሻ ൅ 𝑹்௛ଶ,ோ𝒊௅,ோሺ𝑡ሻ (4.3) 𝒗்ଶ,ூሺ𝑡ሻ ൌ 𝒆்௛ଶ,ூሺ𝑡ሻ ൅ 𝑹்௛ଶ,ூ𝒊௅,ூሺ𝑡ሻ (4.4) Separating (4.2) into real and imaginary parts, combining with (4.3) and (4.4) and considering the voltage drops of the link resistances renders ሺ𝒁்௛ଵ,ோ൅𝑹௟ ൅ 𝑹்௛ଶ,ோሻ𝒊௅,ோሺ𝑡ሻ െ 𝒁்௛ଵ,ூ𝒊௅,ூሺ𝑡ሻ ൌ 𝒆ோሺ𝑡ሻ (4.5) 𝒁்௛ଵ,ூ𝒊௅,ோሺ𝑡ሻ ൅ ሺ𝒁்௛ଵ,ோ൅𝑹௟ ൅ 𝑹்௛ଶ,ூሻ𝒊௅,ூሺ𝑡ሻ ൌ 𝒆ூሺ𝑡ሻ (4.6) 55 where 𝑹௟ is the 𝑛𝑥𝑛 diagonal matrix of link resistances, 𝒁்௛ଵ,ோ and 𝒁்௛ଵ,ூ are the real and imaginary parts of the 𝑛𝑥𝑛 complex-valued matrix 𝒁ഥ்௛ଵ, respectively, and 𝒆ோ and 𝒆ூ are given by 𝒆ோሺ𝑡ሻ ൌ 𝒆்௛ଵ,ோሺ𝑡ሻ െ 𝒆்௛ଶ,ோሺ𝑡ሻ (4.7) 𝒆ூሺ𝑡ሻ ൌ 𝒆்௛ଵ,ூሺ𝑡ሻ െ 𝒆்௛ଶ,ூሺ𝑡ሻ (4.8) Making 𝒊ሺ𝑡ሻ ൌ ቈ𝒊௅,ோሺ𝑡ሻ𝒊௅,ூሺ𝑡ሻ቉ (4.9) 𝒆ሺ𝑡ሻ ൌ ൤𝒆ோሺ𝑡ሻ𝒆ூሺ𝑡ሻ൨ (4.10) 𝒁 ൌ ൤𝒁்௛ଵ,ோ൅𝑹௟ ൅ 𝑹்௛ଶ,ோ െ𝒁்௛ଵ,ூ𝒁்௛ଵ,ூ 𝒁்௛ଵ,ோ൅𝑹௟ ൅ 𝑹்௛ଶ,ூ൨ (4.11) then, from (4.5) - (4.11) we have that 𝒁𝒊ሺ𝑡ሻ ൌ 𝒆ሺ𝑡ሻ (4.12) 𝒊ሺ𝑡ሻ ൌ 𝒁ିଵ𝒆ሺ𝑡ሻ (4.13) Equation (4.13) provides the real and imaginary parts of the link currents at time t. 4.2.2 Subsystems Solutions From the solution of (4.13) 𝒊௅,ோሺ𝑡ሻ and 𝒊௅,ூሺ𝑡ሻ, the real and imaginary link currents necessary to update the subsystems solutions can be extracted. EMT SS The real and imaginary EMT solutions can be calculated injecting the link currents and using the following equations, equivalent to equation (3.9) of Chapter 3. 56 𝒗ଶ,ோሺ𝑡ሻ ൌ 𝒆ଶ,ோሺ𝑡ሻ െ 𝑮ଶ,ோିଵ𝑪ଶ𝒊௅,ோሺ𝑡ሻ (4.14) 𝒗ଶ,ூሺ𝑡ሻ ൌ 𝒆ଶ,ூሺ𝑡ሻ െ 𝑮ଶ,ூିଵ𝑪ଶ𝒊௅,ூሺ𝑡ሻ (4.15) where, 𝒆ଶ,ோሺ𝑡ሻ, 𝒆ଶ,ூሺ𝑡ሻ, 𝒗ଶ,ோሺ𝑡ሻ, 𝒗ଶ,ூሺ𝑡ሻ, 𝑮ଶ,ோ and 𝑮ଶ,ூ are, for the EMT SS, the vectors of open circuit voltages for the solution at time t, the vectors of solution voltages at time t and the conductance matrices, for the real and imaginary solutions, respectively. 𝑪ଶ is the connectivity matrix of SS2. SFA SS To update the SFA SS, the real and imaginary link currents have to be combined in complex-number form using equation (4.16), transformed to the SFA domain using equation (4.17) and then injected into the SFA SS employing equation (4.18). ଙ௅̅ሺ𝑡ሻ ൌ 𝒊௅,ோሺ𝑡ሻ ൅ 𝑗𝒊௅,ூሺ𝑡ሻ (4.16) ଙ௅̅௠ሺ𝑡ሻ ൌ 𝑇ଙ௅̅ሺ𝑡ሻ ൌ 𝑒ି௝ఠబ௧ଙ௅̅ሺ𝑡ሻ (4.17) 𝒗ഥଵ௠ሺ𝑡ሻ ൌ 𝒆തଵ௠ሺ𝑡ሻ െ 𝒀ഥଵି ଵ𝑪ଵଙ௅̅௠ሺ𝑡ሻ (4.18) where ଙ௅̅௠ሺ𝑡ሻ is the vector of link currents at time t in the SFA domain. 𝒆തଵ௠ሺ𝑡ሻ and 𝒗ഥଵ௠ሺ𝑡ሻ are for SS1 the SFA vectors of open circuit voltages for the solution at time t and the solution voltages at time t, respectively. 𝒀ഥଵ and 𝑪ଵ are admittance matrix and the connectivity matrix of SS1. 4.2.3 Single Rate SFA-EMT Solution Algorithm Using the same step size Δt for the SFA and EMT SSs, the MATE solution procedure to solve the entire system at a given time t, is the following: 57 i. Solve for the link currents applying equation (4.13) and using the Thévenin equivalents of the SFA and EMT SSs calculated at the end of the previous solution time. ii. Solve for the node voltages in the subsystems by injecting the link currents into them, that is, by applying equations (4.14 – 4.15) for the EMT SS and equation (4.18) for the SFA SS. This will update the open circuit solutions of the subsystems and generate the solution of the whole system with the links connected. iii. Evaluate the current injections 𝒉ഥଵ for SS1 and 𝒉ଶ,ோ and 𝒉ଶ,ூ for SS2, that would be required for the next solution time at t = t + t. iv. Evaluate the Thévenin equivalents of both subsystems for next solution time at t = t + t, by solving the subsystems in open circuit (with the links disconnected) using the current injections evaluated in step (iii) and applying the equations below. EMT SS: 𝒆ଶ,ோ ൌ 𝑮ଶ,ோିଵ𝒉ଶ,ோ (4.19) 𝒆ଶ,ூ ൌ 𝑮ଶ,ூିଵ𝒉ଶ,ூ (4.20) 𝒆்௛ଶ,ோ ൌ 𝑪ଶ்𝒆ଶ,ோ (4.21) 𝒆்௛ଶ,ூ ൌ 𝑪ଶ்𝒆ଶ,ூ (4.22) 𝑹்௛ଶ,ோ ൌ 𝑪ଶ்𝑮ଶ,ோିଵ𝑪ଶ (4.23) 𝑹்௛ଶ,ூ ൌ 𝑪ଶ்𝑮ଶ,ூିଵ𝑪ଶ (4.24) SFA SS: 𝒆തଵ௠ ൌ 𝒀ഥଵି ଵ𝒉ഥଵ (4.25) 𝒆ത்௛ଵ௠ ൌ 𝑪ଵ்𝒆തଵ௠ (4.26) 𝒁ഥ்௛ଵ ൌ 𝑪ଵ்𝒀ഥଵି ଵ𝑪ଵ (4.27) 58 v. Advance the solution time by Δt, to t = t + Δt and go back to step (i) until t = tmax. 4.3 Multirate SFA-EMT Simulation The SFA-EMT multirate interfacing protocol will use similar equations as the ones derived for the single rate protocol in the previous section and will use a similar approach as the one employed for the EMT-EMT MATE multirate protocol described in Section 3.4. Let us consider the same system as the one indicated in Figure 4.1, but this time the SFA SS will be solved with a large time step equal to ΔTSFA and the EMT SS with a small time step equal to ΔtEMT. For this system, the updating of the SSs is as indicated below. 4.3.1 Updating the EMT Subsystem Let us assume that the SFA SS (SS1) has already been updated at t1 and is now time to update the EMT SS (SS2) at t, where t1 < t < t1 + ΔTSFA. As indicated in Section 3.3, every time a SS gets updated, a Thévenin equivalent (TE) of that SS is generated for the solution at the next time step. That is, after SS1 has been updated at time t1, a vector of SFA equivalent Thévenin sources is calculated for the next solution at time t1 + ΔTSFA. Therefore, the TE of the SFA SS required to update SS2 at time t can be found by linearly interpolating in the SFA domain between the two equivalent Thévenin sources at time t1 and at time t1 + ΔTSFA [78]. To update SS2 at time t, first the link solution is found applying equation (4.13) but using the interpolated TE of SS1 at t, that is, equations (4.1), (4.2) and (4.5) to (4.8) are evaluated employing the vector of interpolated Thévenins of the SFA SS. After that, the solution of the SS2 can be obtained by injecting the real and imaginary link currents extracted from the solution of 59 (4.13) into the respective real and imaginary parts of the EMT SS and applying equations (4.14) and (4.15). Now the real and imaginary current injections for the EMT SS, that would be required for the next solution time at t = t + tEMT are evaluated and then, applying equations (4.19) to (4.24) the Thévenin equivalents of SS2 at t = t + tEMT for the real and imaginary SS2 solutions are calculated. The solution time is then advanced to t + ΔtEMT and the above process is repeated, every time using a new interpolated TE of SS1, until is time to update the SFA SS, that is, until t would become greater than t1 + ΔTSFA [78]. Regarding the Thévenin equivalents used in MATE for the determination of the link currents either in the single rate or the multirate modes, it is worth mentioning that they are derived from the discretization of the differential equations of each of the subsystems and therefore they capture the frequency response of each subsystem for a range of frequencies determined by the integration step used. In particular, the TE of the SFA SS employed in equation (4.13) for the determination of the link current to update the EMT SS can capture the frequency response of the SFA subsystem for a range of frequencies around the fundamental which would be larger or smaller depending on the time step being used. For example, in a 60 Hz system, for time step Δt = 4 ms in SFA, the Nyquist frequency fNy is equal to 𝑓ே௬ ൌ ଵଶ୼௧ ൌ 125 𝐻𝑧; hence, assuming the use of the trapezoidal rule, frequencies up to fNy/5 (that is, േ25 Hz in SFA) can be captured with 3% accuracy [69]. This would correspond to a range in normal time domain from 35 to 85 Hz. This range can be increased by reducing the time step. If a time step of 1.67 ms is used instead, the range of frequencies that could be accurately captured in SFA would become േ60 Hz, which would correspond in normal time domain to a range from 0 Hz to 120 Hz! It follows from this, 60 that transients within the SFA SS or flowing into it through the boundaries and with frequencies within the range determined by the SFA time step, can be accurately captured with this SFA-EMT multirate solution. This is relevant in studies involving harmonic and sub-harmonic resonances, as well as in stability studies with large frequency excursions [78]. 4.3.2 Updating the SFA Subsystem Let us assume now that the EMT SS (SS2) has already been updated at t2 and is now time to update the SFA SS (SS1) at t, where t2 < t < t2 + ΔtEMT. As mentioned in the previous chapter (see Section 3.4 or 3.5), when updating the slow SS (the SFA SS in this case), a low-pass filtered or smoothed out version of the TE of the fast SS (the EMT SS2 in this case) must be used to evaluate the link currents to avoid aliasing errors in the solution. Therefore, to update the SFA SS, all the EMT Thévenin sources obtained in the EMT solutions (real and imaginary) after the last updating of SS1 are combined into complex form, shifted with the SFA transformation and then low-pass filtered or smoothed out by a certain antialiasing processing and down sampled. In this research, with the purpose of maximizing the accuracy of the multirate solution, several antialiasing processing methods were studied which are presented and discussed in the next section. For the present explanation it suffices to say that after the antialiasing and the down sampling have been performed, a smoothed out complex-valued vector of Thévenin voltages of the EMT SS at t is generated. Since this Thévenin vector is in the SFA domain, it is brought back to normal time domain by applying the inverse SFA transformation. This low-pass filtered Thévenin of the EMT SS in normal time domain will be the one used in the updating the SFA SS as indicated below. 61 Even though the solution of SS1 is in the SFA domain, it is easier to solve for the link currents in normal time domain. Hence, the link solution for the SFA SS is also found applying equation (4.13) but using the low-pass filtered TE of SS2, that is, equations (4.3) to (4.8) are calculated employing the real and the imaginary parts of the smoothed out Thévenins of the EMT SS at t. After that, the real and imaginary link currents are extracted from the solution of (4.13), combined into complex form, transformed to the SFA domain and injected into SS1 to update its solution. This is carried out applying equations (4.16) to (4.18). Now the complex current injections for the SFA SS that would be required for the next solution time at t = t + tSFA are evaluated and then, applying equations (4.25) and (4.27), the Thévenin equivalents for the solution of SS1 at t = t + tSFA are calculated. The solution time is then advanced to t2 + ΔtEMT and the updating of the EMT SS is resumed and will continue until is time to update again SS1. The updating of both SSs will continue in that fashion until the maximum simulation time is reached [78]. Note that, as opposed to previous MATE multirate implementations (e.g., [60], [76], [77]), this procedure does not require the large time step to be an integer multiple of the small time step. The above SFA-EMT multirate protocol will be identified in the rest of the thesis as the Standard SFA-EMT MATE multirate protocol to distinguish it from a variant of it that will be discussed later. This standard SFA-EMT MATE multirate protocol is shown graphically in Figure 4.2. 62 Figure 4.2 Standard SFA-EMT MATE multirate protocol 4.3.3 Antialiasing Processing of EMT Solution Variables In order to maximize the accuracy of the multirate solution, several ways of conducting the antialiasing process on the EMT variables were evaluated in this research. These were accomplished by combining three basic low-pass-filtering or smoothing methods with two ways of carrying out the SFA-EMT MATE multirate protocol when updating the SFA SS. The low-pass filtering methods used were: a) Averaging, b) Linear least-square fitting and c) A Savitzky-Golay low-pass filter. The approaches to implement the SFA-EMT MATE multirate protocol were: a) The standard one, that uses low-pass filtered Thévenin voltages in the Thévenin equivalents of the EMT SS employed in the link solution to update the SFA SS, as explained in Section 4.3, and b) An approach that uses alternative methods to find the link solution necessary to update the SFA SS. One of the methods considered in the latter, incorporates energy balance 63 constraints in the derivation of the link current to update the SFA SS. All these antialiasing procedures are described in the following. Let us bear in mind that the antialiasing processing is to be utilized for the updating of the SFA SS at time t; thus, the discussions that follow will be restricted to that part of the protocol only. 4.3.3.1 Standard SFA-EMT MATE Multirate Protocol The standard format of the SFA-EMT MATE multirate protocol, corresponds to the one described in Section 4.3 and shown in Figure 4.2. In this form of the protocol, the SSs are modeled with Thévenin equivalents for the link solution of each SS. Hence, to update the SFA SS the low-pass filtering is performed on the vectors of Thévenin voltages that come from the EMT solution (real and imaginary). As mentioned in Section 4.3, the low-pass filtering is conducted in the SFA domain, therefore the real and imaginary vectors of EMT Thévenin voltages are combined into complex form and transformed to the SFA domain. When the low-pass filtering is done, the vectors of low-pass filtered values of EMT Thévenin voltages are transformed back to normal time domain where they are combined with the corresponding matrices of Thévenin resistances (real and imaginary) to form the Thévenin equivalents of the EMT SS (real and imaginary) to use in the link solution to update the SFA SS. This applies to all three smoothing out methods presented below. i. Averaging In this approach, all the EMT Thévenin sources obtained in the EMT solutions (real and imaginary) after the last updating of SS1 are considered for the average, including an interpolated one evaluated at t. The average is performed as in [60], [79]; however, since 64 in this case the time steps can be arbitrary, they need not be multiples, the averaging is carried out by evaluating the integral below, as in [79]. 𝒆ത்௛ଶ௠,௔௩௘ሺ𝑡ሻ ൌ 1Δ𝑇ௌி஺ න 𝒆ത்௛ଶሺ𝑢ሻ𝑒ି௝ఠబ௨𝑑𝑢௧௧ି୼்ೄಷಲ (4.28) where, 𝒆ത்௛ଶሺ𝑢ሻ ൌ 𝒆்௛ଶ,ோሺ𝑢ሻ ൅ 𝑗𝒆்௛ଶ,ூሺ𝑢ሻ (4.29) 𝒆ത்௛ଶ௠,௔௩௘ሺ𝑡ሻ is the SFA domain average value of the vectors of complex EMT Thévenin sources, 𝒆ത்௛ଶ, over a period of time equal to ΔTSFA. Since this vector results from a plain average it is constant over the entire ΔTSFA interval. This average vector corresponds to an equivalent fundamental frequency phasor of the EMT Thévenin sources. In filter terminology, this averaging procedure is known as moving average filter. Moving average filters are known to introduce a delay in the filtered signal, equal to half the averaging interval [61]. However, the averaging procedure described above does not introduce delay in steady state though, because it is performed on analytic signals transformed to the SFA domain. When the analytic signal is a complex exponential at ω0; that is, the real and imaginary parts are sinusoids at ω0, the transformed signal is a complex number with constant magnitude and phase. Therefore, no matter what section of the period of the complex exponential is being averaged, the average is always the same; but when that average value is brought back to normal time domain at time t, that unshifted value will correspond exactly to the value of the complex exponential at time t with no delay whatsoever! Given that this method extracts the fundamental frequency 65 phasor, not the dynamic phasor, when the system is not in steady state it is expected to introduce some inaccuracy in the solution; especially if the EMT SS is undergoing significant frequency deviations under the influence of the electromechanical transients. Despite of that, the performance of this method was good in the case study where it was evaluated, as can be seen in Chapter 6, Sub-Section 6.1.2 For comparison and because in the method described above the averaging is done only on the EMT Thévenin voltages starting from the last update of the SFA SS, it was decided to test doing the averaging of the EMT Thévenin voltages over a full period of the fundamental and also averaging over a full fundamental period, only the real EMT Thévenin voltages. Notice that, for the last case, equation (4.28) just becomes the known integral for the complex Fourier coefficient C1 of a real signal. ii. Linear Least-square Fitting Because a constant value in the SFA domain corresponds to a fundamental frequency phasor in normal time domain, the averaging of the EMT Thévenin voltages will render only an equivalent fundamental frequency phasor of the EMT solution for the interval under consideration; hence, in order to extract the dynamic phasor from the EMT solution a standard Linear Least-Squares (LLS) regression method (straight line approximation) was used to accomplish the smoothing out process. In this approach, all the EMT Thévenin sources obtained in the EMT solutions (real and imaginary) after the last updating of SS1, including an interpolated one evaluated at t, are 66 considered for the LLS fitting procedure. This procedure was applied in the SFA domain to the complex vectors of EMT Thévenins obtained after combining the EMT solutions in complex form and shifting them. After the fitting is finished, the down sampling is carried out by evaluating the fitted line at t and the value thus obtained is transformed to normal time domain with the inverse SFA transformation to solve for the links using equation (4.13). The LLS fitting procedure does not introduce any delay between the smoothed signal and the original signal, because the smoothed value is obtained evaluating the fitted line at the end of the interval. This process is repeated every time the SFA SS has to be updated. iii. Savitzky-Golay Low-pass Filter Savitzky-Golay (S-G) filters are based on a local least-squares approximation of a set of input samples, to a polynomial of specified order. S-G filters are appropriate to match the waveform of an oversampled signal corrupted by noise. They have extremely flat passbands with moderate attenuation in their stopbands and have the important property that they do not introduce a phase delay; thus, the components of the signal are not distorted while some high-frequency noise is attenuated, though not completely eliminated [80]. Following the recommendations in [61] and with the purpose of testing a more elaborate filter than the two previously discussed, we also tested this filter to check the accuracy of the multirate solution. 67 The S-G filter requires the samples to be equally spaced and the number of samples to be an odd number. Hence, in this case the highest odd number of EMT Thévenins obtained in the EMT solutions since the last update of the SFA SS are used. These Thévenins are combined in complex form and shifted to the SFA domain where the S-G low-pass filtering is applied. The filtered data is then evaluated at t to carry out the down sampling. The resultant vector is brought back to normal time domain with inverse SFA transformation to be used in equation (4.13) for the link solution to update the SFA SS. This process is repeated every time SS1 has to be updated. All proposed filters were tested and compared for Case 1b in Chapter 6. 4.3.3.2 Alternative SFA-EMT MATE Multirate Protocol As mentioned earlier, an alternative approach of implementing the SFA-EMT MATE multirate protocol when updating the SFA SS was also considered. This approach consisted in making use the EMT link currents or other variables of the EMT solution, to find the link solution necessary to update the SFA SS, instead of utilizing the method of the standard protocol applying equation (4.13). Two methods were considered for this approach, namely: i. Equivalent EMT Link Current In this method, the SFA SS is updated using a link current equal to a smoothed value of the prior EMT link currents. To accomplish this, the SS2 link currents both, real and imaginary, for the time interval starting from the last update of SS1 until t, the time to update the SFA SS, are combined into complex form and shifted to the SFA domain 68 applying the SFA transformation. Once in the SFA domain, those currents are low-pass filtered using any of the three methods described before; namely, Averaging, LLS fitting or S-G low-pass filter. The smoothed value evaluated at t is injected into SS1 to update it making use of equation (4.18). Notice that for this calculation equation (4.13) is not utilized. This process is repeated every time SS1 has to be updated. ii. Power Balancing Link Current In the standard single rate MATE, the balance of the power/energy exchanged between subsystems is achieved automatically through the solution itself. This is because the Thévenin equivalents are exact and all subsystems are updated at the same time using the same link currents. This is not necessarily the case for the general MATE multirate, when the time steps are not multiples and antialiasing processing of the fast solution is required because, when the link currents to update a subsystem are found, only that subsystem’s own Thévenin equivalent is exact, the other Thévenin equivalent is either interpolated or smoothed and, therefore, is approximate. Hence, every time a subsystem is solved, whatever power is received/supplied by that subsystem at that time is not supplied/received by the other subsystem, but by an equivalent circuit of it which is approximate and vice versa. As there are no explicit power/energy balance constraints in the MATE multirate solution, energy exchange balance is not guaranteed and may not be achieved. In view of this, the idea arose of incorporating explicit energy balance constraints in the updating of the SFA SS to check if this would improve the accuracy of the multirate 69 solution in comparison with the standard formulation of the SFA-EMT MATE multirate based on pure Thévenin equivalents. To include such constraints, the following considerations have to be made. First, the granularities of the SFA and the EMT solutions are different. Due to its larger time step, the SFA solution points are more spaced in time than the EMT solution points. Everything that happens in between SFA solution points is unavailable to the SFA SS, it gets news from the EMT solution only at its solution points and this implies that any power/energy balance constraints need to be considered only at the SFA solution points. Second, let us consider two cases here. a) The particular case, where the Nyquist frequency of the SFA SS defined by its time step is higher than the highest frequency of the Thévenins of the EMT solution. In this case there is no need for low-pass filter the EMT solution, since there is no risk of aliasing; so, the EMT Thévenins are ready to be used in the SFA solution, except that they may need to be interpolated at the SFA updating time. Energy balance is not an issue under this condition. b) The general case, where the Nyquist frequency of the SFA SS is lower than the highest frequency of the Thévenins of the EMT solution. In this case the total energy absorbed/supplied by the EMT SS between SFA solution points may include energy associated to frequencies higher than the Nyquist frequency of the SFA SS, therefore this total energy cannot be used as an energy balance constraint; or said in another way, the instantaneous power determined by the EMT terminal voltages and the currents flowing through the links into the EMT SS at the SFA updating time would not work as a power balance constraint to 70 update the SFA SS, because this would cause inaccuracies in the solution (“aliasing”) as the energy content of the frequencies present in the SFA solution would be altered to accommodate alien energies associated to higher frequencies that do not exist in its solution but that are, nevertheless, present in the EMT solution. An effective way to avoid this problem and incorporate power/energy exchanged balance constraints in the solution is by smoothing out both, the instantaneous power at the terminals of the EMT SS and the terminal voltages, and determining a current that matches that smoothed power with the smoothed voltage and using that current to update the SFA SS. By balancing the power at the solution times of the SFA SS, the associated energy is also balanced. This can be done using the following procedure. Determination of Smoothed Link Current for Power Balance Using the equations developed in Chapter 2, Section 2.4 and extending them to multi-port circuits, the vector of total instantaneous complex power 𝒔ത௧௢௧ଶሺ𝑡ሻ absorbed by the EMT SS at any solution time t, is given by 𝒔ത௧௢௧ଶሺ𝑡ሻ ൌ 𝒑௧௢௧ଶሺ𝑡ሻ ൅ 𝑗𝒒௧௢௧ଶሺ𝑡ሻ (4.30) 𝒑௧௢௧ଶሺ𝑡ሻ ൌ 𝒗்ଶ,ோሺ𝑡ሻ ⊙ 𝒊௅ଶ,ோሺ𝑡ሻ ൅ 𝒗்ଶ,ூሺ𝑡ሻ ⊙ 𝒊௅ଶ,ூሺ𝑡ሻ (4.31) 𝒒௧௢௧ଶሺ𝑡ሻ ൌ 𝒗்ଶ,ூሺ𝑡ሻ ⊙ 𝒊௅ଶ,ோሺ𝑡ሻ െ 𝒗்ଶ,ோሺ𝑡ሻ ⊙ 𝒊௅ଶ,ூሺ𝑡ሻ (4.32) Where 𝒑௧௢௧ଶ and 𝒒௧௢௧ଶ are the vectors of total instantaneous active and reactive powers absorbed by the EMT SS. These vectors contained the combined effect of the active or reactive power of both, the real and imaginary parts of the EMT solution for the terminal voltages and the link currents. The number of elements of these vectors is equal to the number or links. 𝒗்ଶ,ோ, 𝒗்ଶ,ூ, 𝒊௅ଶ,ோ and 𝒊௅ଶ,ூ are, respectively, the real and imaginary 71 solutions for the vector of terminal voltages and the vector of link currents of SS2 at a certain time t. ⊙ denotes elementwise multiplication of vectors. Equations (4.30) - (4.32) are evaluated in normal time domain; however, they would render the same results if they were evaluated in the SFA domain because the SFA transformation is total instantaneous complex power invariant (see Chapter 2, Section 2.4). Note that the subscript “tot” for the powers does not refer to the fact that those powers account for the contributions from all links, which they do; instead, the subscript in this context specifically refers to the fact that those powers include the contributions from both, the real and the imaginary parts of the voltages and currents. Now the vector of terminal voltages in the SFA domain of SS2 at any solution time t, 𝒗ഥ்ଶ௠ሺ𝑡ሻ, can be determined applying equations (4.33) and (4.34) below 𝒗ഥ்ଶሺ𝑡ሻ ൌ 𝒗்ଶ,ோሺ𝑡ሻ ൅ 𝑗𝒗்ଶ,ூሺ𝑡ሻ (4.33) 𝒗ഥ்ଶ௠ሺ𝑡ሻ ൌ 𝑇𝒗ഥ்ଶሺ𝑡ሻ ൌ 𝑒ି௝ఠబ௧𝒗ഥ்ଶሺ𝑡ሻ (4.34) 𝒔ത௧௢௧ଶ and 𝒗ഥ்ଶ௠ have to be evaluated for all EMT solutions from the one that occurred right after the last update of SS1 until the latest prior to the new updating of SS1. These sets of solutions of 𝒔ത௧௢௧ଶ and 𝒗ഥ்ଶ௠ are then smoothed out using any of the low-pass filtering methods described in Sub-Sub-Section 4.3.3.1 and evaluated at t to obtain the corresponding smoothed values 𝒔ത௧௢௧ଶ,௟௣௙ሺ𝑡ሻ and 𝒗ഥ்ଶ௠,௟௣௙ሺ𝑡ሻ. The vector of link currents ଙ௅̅ଵ௠,௣௕, with power balance constraint, to update SS1 can be evaluated with the following equation 72 ଙ௅̅ଵ௠,௣௕ሺ𝑡ሻ ൌ ቀ𝒔ത௧௢௧ଶ,௟௣௙ሺ𝑡ሻ ⊢ 𝒗ഥ்ଶ௠,௟௣௙ሺ𝑡ሻቁ∗ (4.35) where ⊢ denotes element wise division of vectors. By injecting this current to solve SS1 we are performing an antialiasing processing with a power balance constraint, since ଙ௅̅ଵ௠,௣௕ has been derived from a smoothed instantaneous power and a smoothed voltage, which do not contain contributions from high frequency components above the Nyquist frequency of the SFA SS. It is worth mentioning however, that this power balance constraint is approximate, since the actual power absorbed or supplied by the SFA SS will only be known when that subsystem is solved and it may actually differ from the one used to derive the current, so there could be a small power mismatch. 4.4 Selection of Boundaries Between Subsystems The use of hybrid phasors/EMT modeling does not result out of the desire to take advantage of the existence of some natural latency in the system, but out of necessity; that is, a need for accuracy and a need for computational efficiency. Accuracy is needed because detailed device level modeling (EMT-type) is required for certain components of the system, to evaluate the impact of their fast/non-linear dynamics on the stability of the system. For very large power systems, system wide functional modeling (phasors-type) is needed to achieve computational efficiency because an all-EMT solution would be too computationally onerous. Thus, selecting the boundary between subsystems in a hybrid simulation is a compromise between computational efficiency and the accuracy of the solution. Ideally for maximum computational 73 efficiency, the size of the EMT subsystem should be as small as needed for an accurate solution. The effect of this minimization on the accuracy will depend on the particular hybrid method used and the ratio between the large- and small-time steps. On the other hand, for maximum accuracy, the EMT subsystem should be made much larger than just the region or the equipment of interest, so that the boundaries are sufficiently far from the source of harmonics and fast transients, but this would be detrimental to the computational efficiency [78]. Some hybrid TS-EMT techniques, (e.g., [23]), use an FDNE in the Thevenin equivalent used to represent the TS subsystem in the EMT subsystem, which makes it valid over a wide frequency range and allows the TS-EMT boundary to be moved closer to the source of transients without compromising the accuracy [18]. This approach, however, only considers the effect of the AC system on the EMT system, but not the effect of the EMT system on the AC system. In terms of boundaries, these techniques would require, a minimum buffer zone around the core part of interest in the EMT subsystem. Other TS-EMT methods (e.g., [32], [81]), that use multi-port equivalents to represent the TS in the EMT subsystem and employ three sequence-based TS simulators (symmetrical components), would require, a larger buffer region than the previous ones, since the impedances of the equivalents are evaluated at a single frequency. Hybrid TS-EMT methods, as in [26], that use simplified single-port Norton equivalents to represent the TS subsystem in the EMT subsystem and positive sequence-based TS simulators (e.g., [26], [29]) would require a buffer region of a still bigger size to achieve satisfactory accuracy. On the other side, the extraction of the fundamental frequency variables at the boundary, which all hybrid methods need to do, gets easier and the accuracy improves, as the buffer region gets larger 74 because the voltages and currents at the boundary get closer to sinusoidal and balanced. Thus, having at least a minimal buffer region is in general desirable [78]. The proposed method in this research uses multi-port Thévenin equivalents to represent all three phases of the dynamic phasors of SFA subsystem as seen by the EMT subsystem and, as explained previously (see Sub-Section 4.3.1), these equivalents are valid for a frequency range around the fundamental frequency. These aspects of the proposed method ameliorate the problem of the buffer zone and good accuracy is obtained with boundaries relatively close to the source of transients and with large multirate ratios [78]. 75 Chapter 5: Description of simulators employed MATE is a powerful and versatile solution framework but its implementation with commercially available simulators is not easy, as it requires the solution of the subsystems to be found in two steps, first the subsystems are solved independently disconnected from each other and second, their actual solutions are found by injecting the corresponding link current to each of the subsystems. Since the internal algorithms of commercial simulators are not accessible to the user, implementing MATE using existing industry-grade simulators is very difficult. An alternative method is when the software already has the capability to provide the Thévenin equivalents of the system to handle certain network components, like non-linear elements. Examples of these simulators is MICROTRAN [7]; but, since the feature is specific to handling non-linear elements, its application to other purposes like a hybrid simulator is limited for research purposes. For these reasons, in this research project the SFA and the EMT simulators to implement the hybrid multirate protocol, together with the link solver were all developed in house. In this chapter the main features of the hybrid simulator developed for the research are presented. 5.1 General Characteristics of the Hybrid Simulator The hybrid simulator has been developed as an integrated piece of software and is coded in Python 3. It consists of three main solvers: the SFA solver, the EMT solver and the link solver. The latter executes the interfacing protocol and the code to handle the input/output of data. The SFA and EMT solvers share the same code, which can be set as a different object in the program, 76 depending on whether that code is to be used as an SFA or an EMT simulator. In Figure 5.1 a general flowchart of the simulator is presented. The code for the SFA and EMT solvers used in this research was developed by extending an existing single-phase version of the SFA and EMT solvers developed in a previous research project by Andrea Martí [43]. 5.2 The Link Solver The link solver developed in this project contains the implementation of the hybrid SFA-EMT protocol in either multirate or single rate mode, as described in sections 4.3 and 4.2, respectively. The activities of the protocol are executed sequentially, which corresponds to a serial implementation of the protocol. In multirate mode the two EMT solutions can be parallelized and in single rate mode the SFA and the EMT’s solutions can be parallelized as well, however the implementation for this work does not include parallelization. The link solver supports all the antialiasing processing methods described in Sub-Section 4.3.3. It supports the use of Critical Damping Adjustment (CDA) [82] in the EMT solutions to prevent numerical oscillations when using trapezoidal, however, in the SFA solution that feature is not available. To avoid numerical oscillations in the SFA solution the SFA solver uses the backward Euler integration rule. 77 Figure 5.1 General flow chart of the integrated hybrid SFA-EMT simulator. StartInput data of subsystems (SS’s)Steady state solution to initialize SS’s.tSFA = 0 + TSFA, temt = 0 + temtEvaluate Thevenins of SS’s.Interpolate Thevenin of SS1Solve for link currents to update EMT SS real and imaginary (r&i)Inject link currents into SS2 (r&i) to update EMT solutiontemt = temt + temtEvaluate Thevenin eq. of SS2 r&iDo any required topology changes in SS2 r&iTime to update SS1?Interpolate Thevenin of SS2 r&iModified MATE?Do antialiasing processing to determine low‐pass filtered Thevenin eq. of EMT SSDo antialiasing processing to determine low‐pass filtered equiv. current source of EMT SSSolve for link current to update SFA SSInject current into SS1 to update SFA solutiontSFA = tSFA + TSFAEvaluate Thevenin eq. of SS1End simulation?Begin transient solutionEndYes NoNoNoYesYes78 5.3 The SFA and EMT Solvers The SFA and EMT solvers use the same solution algorithm, except that the SFA solution is in the complex domain and the EMT solution in the real domain. They can both use, either the trapezoidal or backward Euler integration rules for the solution. The EMT solver can use CDA to avoid numerical oscillations at discontinuities, like switchings, when using the trapezoidal integration rule. The SFA solver does not support this feature. The simulation can be carried out, either in units or in per unit. The modeling capabilities of these solvers in the present implementation regarding circuit elements are indicated below:  Single-phase and coupled three-phase lumped resistors, inductors and capacitors. Transmission lines can be modeled as single or cascaded π-circuits.  Ideal current or voltage sources.  Ideal single-phase transformers. Three-phase transformers are modeled as three-phase banks of single-phase units. Modeling follows the guidelines in section 6.8 of [69].  Ideal single-, double- or three-phase switches.  Ideal single-, double- or three-phase to ground faults.  Three-phase loads, modeled as constant impedance.  Six-pulse diode bridges. This capability is available in the EMT solver only.  Generators are modeled in both solvers with a classical generator model which corresponds to a voltage source behind a constant transient reactance. This model is commonly used in transient stability studies, it is a simplified one and is described by the second order differential equation known as the swing equation [57], [83]. 79 The swing equation only accounts for the electromechanical interactions in the generator, therefore in this model the dynamics of the field circuit and the electromagnetic interactions between the stator and the rotor are neglected. As a result, the voltage source of the model displays a constant magnitude and a variable angle during the transients. This angle is defined by the relative position of the rotor with respect to a synchronously rotating reference frame [57], [83]. The swing equation is solved in the same way by the SFA and EMT solvers, the user can choose to use either the trapezoidal or the backward Euler rule. In the SFA solver, the three-phase instantaneous electrical active power of the generator is calculated in the normal time domain by unshifting from the SFA domain the voltage and current of its three phases. The solvers only model the generator itself, governors, exciters or power system stabilizers are not considered in the present implementation. 80 Chapter 6: Case Studies The purpose of this chapter is to present the test cases considered to validate and show the capabilities the SFA-EMT interfacing protocol derived in this research, for different types of systems. Two separate cases were analyzed, one, the IEEE 39 bus test system, to represent a conventional power system and the other, the CIGRE benchmark system, to represent a power system with HVDC and power electronics components. Within each of these cases, various sub-cases were also studied in order to show the accuracy of the proposed protocol and to assess the performance of the proposed antialiasing methods. 6.1 Case 1: IEEE 39 Bus Test System Under Fault Conditions To assess the performance of the proposed interfacing protocol, the IEEE 39-bus test system [84] was simulated. This system was partitioned into two subsystems, a large one identified as SFA Subsystem (SS) or Subsystem 1 (SS1), and a small one, identified as the EMT SS or SS2, as shown in the one-line diagram of Figure 6.1. This system does not exhibit any natural latency between the dynamic responses of the two subsystems selected; thus, the selection of the boundary between subsystems was done considering the closeness to the disturbance while maintaining the EMT SS small. In other words, it was done thinking of a compromise between accuracy and computational efficiency [78]. SS1 was simulated with the SFA solver using a time step of 3.96 ms and SS2 was simulated with the EMT solver with a time step of 0.1 ms. The simulation consisted of applying a solid 3-phase fault on the line between buses 26 and 28, very close to bus 28, at t = 50 ms, which was cleared 81 after 8 cycles (60 Hz) by tripping that line. The dynamic behaviour of the system was recorded up to 2 s. Generators were modelled with the classical model, lines with π-sections and loads were represented as constant impedances. In the Appendix the system data used, modified from [84], is presented. Figure 6.1 IEEE 39-bus system and SFA and EMT subsystems for case study 1. 82 The purpose of the study is to evaluate the electromechanical transients in SS1 and to capture the electromagnetic transients in SS2 using the proposed interfacing protocol. To evaluate the performance of the protocol, the results are compared with a reference solution conducted by simulating the entire system as a single unit with an EMT solver using a single small time-step of 0.1 ms and using the trapezoidal rule with Critical Damping Adjustment (CDA) [82]. To prevent numerical oscillations in the SFA-EMT solutions, in all simulations with the hybrid protocol, the SFA solver employs the backward Euler integration rule and the EMT solver uses the trapezoidal integration rule with CDA. For Case 1 two sub-cases were considered. Case 1a is for the single rate solution and Case 1b for the multirate solution. 6.1.1 Case 1a: Hybrid SFA-EMT Single Rate Solution For this case the hybrid multirate protocol was used in single rate mode employing a time step of 0.1 ms for both the SFA and the EMT solutions. Results for the SFA SS are presented in Figure 6.2 to Figure 6.4 and for the EMT SS in Figure 6.5 to Figure 6.8. Since the interest in this case is to verify that the hybrid single rate solution can match the EMT only solution to the fastest transients, for greater resolution of the curves in all these graphs only a fraction of the total simulation is shown. 83 Figure 6.2 Frequency of generator 8 in the SFA subsystem (single rate solution). Figure 6.3 Electrical power of generator 8 in the SFA subsystem (single rate solution). 84 Figure 6.4 Voltage at bus 27 (phase “a”) in the SFA subsystem (single rate solution). Figure 6.5 Fault current from bus 28 (phase “a”) in the EMT subsystem (single rate solution). 85 Figure 6.6 Fault current from bus 28 (phase “b”) in the EMT subsystem (single rate solution). Figure 6.7 Fault current from bus 28 (phase “c”) in the EMT subsystem (single rate solution). 86 Figure 6.8 Voltage at bus 28 (phase “c”) in the EMT subsystem (single rate solution). As it can be observed from the inspection of the previous figures, the hybrid SFA-EMT solution in single rate mode coincide with the full EMT solution, the maximum percentage error is 0.06% and occurred in the electric power of generator 8, all other variables have much less error. This validates the accuracy of the proposed interface model between an SFA solution and an EMT solution in single rate mode for this type of systems. 6.1.2 Case 1b: Hybrid SFA-EMT Multirate Solution For all cases in this sub-section, as mentioned previously, SS1 was simulated with the SFA solver, using a time step of 3.96 ms and SS2 with the EMT solver and a time step of 0.1 ms, which corresponds to a rate ratio of 39.6. The reason for choosing this peculiar value of the time 87 step for the SFA solution was to show that the proposed method can work for time steps that are not integer multiples and can indeed be arbitrary. The results for the SFA subsystem are presented in Figure 6.9 to Figure 6.12. For the SFA solution the interest is on the electromechanical transients of the system and it can be easily appreciated that those transients are accurately captured by the hybrid multirate solution proposed. The error in the frequency for generator 8 is less than 0.01%. In Figure 6.12 it can be observed how the magnitude of the dynamic phasor of the SFA solution accurately captures the envelope of the voltage time signal. Figure 6.9 Frequency of generator 8 in the SFA subsystem (multirate solution). 88 Figure 6.10 Angle difference between generators 8 and 1 in the SFA subsystem (multirate solution). Figure 6.11 Electric power of generator 8 in the SFA subsystem (multirate solution). 89 Figure 6.12 Voltage at bus 27 (phase “a”) in the SFA subsystem (multirate solution). The results for the EMT subsystem are shown in Figure 6.13 to Figure 6.18 below. Since for this subsystem we are interested in the electromagnetic transients, only a time window of about 200 ms in presented. From these results it can be observed that the hybrid multirate solution captures with good accuracy the electromagnetic transients in the EMT subsystem. The errors in the first peak of the fault currents are for phases “a”, “b” and “c”, respectively, 3.3 %, 5.5 % and 2.0 %, these are associated to the inception of the fault. The maximum errors in the voltages, which are all associated to the clearance of the fault are, respectively, 3.2 %, 7.9 % and 1.7 %, for phases “a”, “b” and “c”. If desired, these errors can be reduced by either, moving the boundary between subsystems farther from the area of the disturbance, that is, by providing a larger buffer zone, or by reducing the time step for the solution of the SFA subsystem. 90 Figure 6.13 Fault current from bus 28 (phase “a”) in the EMT subsystem (multirate solution). Figure 6.14 Fault current from bus 28 (phase “b”) in the EMT subsystem (multirate solution). 91 Figure 6.15 Fault current from bus 28 (phase “c”) in the EMT subsystem (multirate solution). Figure 6.16 Voltage at bus 28 (phase “a”) in the EMT subsystem (multirate solution). 92 Figure 6.17 Voltage at bus 28 (phase “b”) in the EMT subsystem (multirate solution). Figure 6.18 Voltage at bus 28 (phase “c”) in the EMT subsystem (multirate solution). 93 Antialiasing Processing The antialiasing processing method employed in this simulation for the updating of the SFA subsystem was a Savitzky-Golay (S-G) low-pass filter using a polynomial of order 2 with the standard MATE procedure. This method was selected after testing all the methods presented in Chapter 4, Sub-Section 4.3.3 and verifying that it offered the best performance in terms of the accuracy of the multirate solution with respect to the reference EMT solution. The comparison was made using the error in the first peak of phase “b” of the fault current from bus 28. In Table 6.1 the results are presented. Table 6.1 Error in first peak of fault current from bus 28 (phase "b") for various antialiasing methods. AAPMError (%) ‐6.2 ‐7.4 ‐7.7 ‐5.8 ‐5.5AAPMError (%) ‐6.2 ‐5.8 ‐5.4 ‐6.3 ‐5.7 ‐5.5Legend:AAPM: Antialiasing processing methodAVG: AveragingAVG_1P: Averaging over 1 period of the fundamentalFFP_RP: Fundamental frequency phasor extraction from the real part solutionLLSF: Linear least‐square fittingSGF: Savitzky ‐Golay low‐pass filter All methods as explained in Chapter 4, Sub‐Section 4.3.3 LLSF SGFThevenin EquivalentLink Solution to Update SS1Standard SFA‐EMT MATE Multirate ProtocolSGFAlternative SFA‐EMT MATE Multirate ProtocolLink Solution to Update SS1Equivalent Link Current Power Balancing Link CurrentAVG LLSF SGF AVG LLSFAVG AVG_1P FFP_RP94 Various observations can be made from the error measurements in Table 6.1. First of all, the S-G filter has the best performance of all antialiasing methods in both, the standard and the alternative implementations of the SFA-EMT MATE Multirate protocol. This means that when the SFA-EMT protocol uses the S-G filter as the antialiasing method, the multirate solution has the best accuracy with respect to the reference EMT solution, when compared with the multirate solutions obtained using any of the other antialiasing methods. On the other hand, in general all antialiasing methods have comparable performances, there are no radical differences between them. This indicates that for this simulation they all seem good methods for smoothing out the data, which could be related to the fact that neither of them introduces delays in the processing; nevertheless, in this case the high frequency transients decay very rapidly and what stays longer is the exponentially decaying DC component, which has mainly low frequencies associated to it. Regarding the alternative approach that utilizes a smoothed value of the EMT link currents to update the SFA SS, the accuracy of the SFA-EMT multirate solution is essentially the same as with the standard implementation of the hybrid multirate protocol, which was to be expected since the link currents are determined using the Thévenin equivalents of the EMT subsystem. The S-G filter had a slightly better performance for this approach, but it was preferred to keep the hybrid interfacing protocol with the standard implementation. 95 Regarding the alternative approach that employs a power balancing link current to update the SFA SS, the performance of the antialiasing processing with respect to the homologous methods with the standard implementation is very similar and three conclusions can be extracted from that. First, this result validates the method derived for carrying out the antialiasing processing with a power balance constraint. Second, it confirms the applicability of the proposed concept of total instantaneous complex power in circuits with complex signals as derived in Chapter 2. Third, it proves that, the standard implementation of the SFA-EMT MATE multirate protocol, based on Thévenin equivalents, implicitly enforces power exchange balance in the multirate solution. This indicates that with the proposed SFA-EMT MATE multirate interfacing protocol there is no need to incorporate explicit power exchange balance constraints in the antialiasing processing, which adds another strength to the MATE multirate solution approach. Regarding the averaging or moving averaging approach derived in this research, even though it just uses the samples for a time span equal to the time step of the SFA subsystem, it outperforms the approach of averaging over a full period of the fundamental frequency or extracting the fundamental frequency phasor from the real part solution. Thus, this method is computationally more efficient and it allows to extract an equivalent fundamental frequency phasor that considers both, the real and the imaginary EMT solutions, without requiring more elaborate procedures such as Fourier analysis or least-square fitting techniques. However, since in our case the interest is to extract the dynamic phasor from the EMT solution, as modern power systems are expected to exhibit larger frequency excursions in their dynamic behavior due to their lowering inertia, other methods like linear least-square fitting or S-G filters would be preferred as they are not restricted to an equivalent fundamental frequency phasor output. 96 Computational Gain Evaluation In order to evaluate the computational gain provided by the hybrid multirate simulation with respect to the reference simulation using an all EMT solution, the execution times of both simulations were evaluated for Case 1b. The results are presented in Table 6.2 below. Table 6.2 Execution times for simulations of Case 1b It can be readily observed from these results that only a moderate gain in computational efficiency was achieved, around 32%, for a significant rate ratio of almost 40. Three factors contribute to limit the computational gain for the hybrid multirate SFA-EMT simulation. First the SFA solution is executed in the complex domain; thus, for the same system it involves roughly double the number of operations than an EMT solution with the same time step. To compensate for this, a more efficient SFA implementation would be required, for example the ones presented in [85] and [86] which are based on symmetrical components. Second, for the EMT solution the interaction protocol proposed requires two identical simulations, one for the real part and one for the imaginary part; hence, for the same system, the number of operations is duplicated on the EMT side in the hybrid solution. Because the proposed hybrid protocol uses MATE, this second factor can be readily compensated by executing the two EMT simulations of the protocol in parallel in both stages of the process, for the open circuit solution to determine the Thévenin equivalents of the EMT subsystems and for updating their solutions by injecting the tEMT TSFAReference: all EMT 0.1 161Hybrid SFA‐EMT 0.1 3.96 39.6 122SimulationTime step [ms]Rate Ratio Execution time [s]97 link currents. However, the parallel processing of the two EMT solutions is not done in the present implementation of the hybrid simulator. And third, in the multirate solution the execution time of the link solver has to be considered, and since all link solutions required in the protocol have to be conducted in complex numbers this further increases the computational overheads associated to operation of the link solver. The impact of this factor could be minimized by optimizing the code of the link solver for computational efficiency, however, this has not been done for the present implementation of the simulator. Considering the above factors, significant computational gains can be expected for large power systems where the SFA subsystem can be made much bigger than the EMT subsystem, with the EMT solutions parallelized and for cases where the EMT solution requires much smaller time steps than the ones used for SFA, that is, for cases with higher rate ratios. In order to have an approximate theoretical estimation of the computational gain in such conditions, Table 6.3 presents the ratio of the execution time of an all EMT solution with respect to the execution time of the hybrid SFA-EMT multirate solution, which is a measure of how many times faster is the SFA-EMT solution in comparison with an all EMT solution and therefore serves as an indicator of the computational gain. 98 Table 6.3 Estimated computational gain of the SFA-EMT solution with respect of an all EMT solution, for various relative sizes of the SFA and EMT subsystems and various rate ratios of the multirate solution. The evaluation was made for various sizes of the SFA subsystem and for different values of the rate ratio of the SFA-EMT solution ሺ𝑟𝑎𝑡𝑒 𝑟𝑎𝑡𝑖𝑜 ൌ ୼்ೄಷಲ୼௧ಶಾ೅ሻ. The data in the table is also presented graphically in Figure 6.19. The values in the table and in the graph were calculated assuming that the execution time of the SFA solver is exactly twice the execution time of an EMT solution for the same system and with the same time step, and neglecting the execution time of the link solver for the multirate solution. SFA SS EMT SS50 50 1.0 1.9 1.0 2.0 1.0 2.075 25 1.9 3.6 1.9 3.8 2.0 3.980 20 2.3 4.3 2.4 4.6 2.5 4.990 10 4.2 7.4 4.6 8.5 4.9 9.695 5 7.2 11.4 8.4 14.5 9.5 18.399 1 16.8 20.2 25.1 33.6 40.1 66.999.5 0.5 20.1 22.3 33.4 40.2 66.8 100.3Note: This estimation does not consider the execution times of the l ink solver for the SFA‐EMT protocol.Size of Subsystems in % of System's SizeSFA‐EMT Protocol's Rate Ratio = SFA / tEMT50 100 400EMT Sols. Not Paral.EMT Sols. Parallel.EMT Sols. Not Paral.EMT Sols. Parallel.EMT Sols. Not Paral.EMT Sols. Parallel.99 Figure 6.19 Estimated computational gain of the SFA-EMT multirate solution with respect of an all EMT solution, as a function of the size of the SFA subsystem and the rate ratio of the multirate solution (This estimation does not consider the execution time of the link solver for the multirate solution). Although the values in Table 6.3 and in Figure 6.19 are approximate and the actual execution times of the SFA-EMT multirate protocol would be higher, they do show that the proposed SFA-EMT multirate solution can be significantly faster than an all EMT solution, when the SFA subsystem is much larger than the EMT subsystem, the rate ratio of the hybrid multirate protocol is high and the EMT solutions are parallelized. It is important to note that these conditions are precisely the ones associated with the simulation of actual power systems, which are very large and may span over various countries or even continents. In those systems, a large percentage of 100 the system would be modeled in SFA and for the rest it, an EMT model would be used. On the other hand, typical integration steps that would be used for the SFA solution to study the transient stability of the system would be in the order of 3-5 ms, and the time steps necessary in the EMT solution to simulate the transient behaviour of power electronics-interfaced devices would be in the order of 10-50 s; so, rate ratios would be in the range of 60-500. Therefore, in the simulation of actual power systems, a simulator operating with the proposed SFA-EMT multirate protocol is expected to be considerably faster than an all EMT simulator. 6.2 Case 2: Modified CIGRE HVDC Benchmark System With the purpose to evaluate the performance of the proposed interfacing protocol in a system that contains power electronics-based devices, a modified version of the CIGRE HVDC benchmark system [87] as indicated in Figure 6.20 was simulated. More specifically, the objective of this case study was to test the protocol in a system containing power electronic devices, whose operation involves internal recurrent or periodic switchings; because, under these conditions the topology of the EMT subsystem for the real and the imaginary parts would not be the same at every instant of time; so, the impact of that on the protocol needed to be evaluated. The modifications to the benchmark system were as follows. Only the AC system on the rectifier side, the rectifier station and the DC line were modeled. The inverter station was modeled with an equivalent resistive load of 247.5 Ω. The 3-winding 345/211.42/211.42 kV, 1196 MVA converter transformer was modeled as two 345/211.42 kV, 598 MVA, three-phase banks of three single-phase ideal transformers, one with a grounded star – delta connection and the other with a grounded star – ungrounded star connection. The thyristors of the 6-pulse bridges of the rectifier 101 station were modeled as diodes; that is, with no firing control. The rest of the data of that system is as reported in [87]. Figure 6.20 Modified CIGRE HVDC benchmark system and SFA and EMT subsystems for case study 2 The system was partitioned into two subsystems, the equivalent AC system was identified as the SFA Subsystem (SS) or Subsystem 1 (SS1) and the HVDC part, that includes the rectifier converter station, the filters, capacitive compensators and all the DC part was identified as the EMT SS or SS2. This is shown in the diagram of Figure 6.20. As it can be readily observed, the interface between subsystems is right on the rectifier converter station, the device whose dynamic behavior is to be studied; so, there is no buffer zone between the SFA SS and the source 102 of transients and, therefore, there is no latency between the dynamic responses of the two subsystems. SS1 was simulated with the SFA solver, using a time step of 0.49 ms and SS2 was simulated with the EMT solver with a time step of 0.05 ms (rate ratio: 9.8). The simulation consisted of the startup of the rectifier station from a blocking condition of its valves; that is, in steady state at initial conditions (t = 0 s) everything was connected but all valves were kept in blocking condition. For the transient simulation (t = 0+) the blocking condition of the valves was removed and they were allowed to operate normally. The dynamic behaviour of the system was recorded for 300 ms. The purpose of the study is mainly to capture the electromagnetic transients in SS2 using the proposed interfacing protocol. In this case there are no electromechanical transients in SS1 as the generator is an ideal source. To evaluate the performance of the protocol, the results are compared with a reference solution conducted by simulating the entire system as a single unit with an EMT solver using a single small time-step of 0.05 ms and using the trapezoidal rule with Critical Damping Adjustment (CDA) [82]. As mentioned in the previous section for Case 1, in all SFA-EMT simulations with the hybrid protocol, the SFA solver uses the backward Euler integration rule and the EMT solver uses the trapezoidal integration rule with CDA. This action effectively prevents any numerical oscillations from arising in the hybrid SFA-EMT solutions. 103 The sub-cases cases presented below are as follows. Case 2a is for the single rate solution, Case 2b is the multirate solution with the rate ratio as indicated above and Case 2c is a multirate solution with a rate ratio of 80.2, to study just the dynamics of the system when the startup process has ended and the rectifier station is operating normally. 6.2.1 Case 2a: Hybrid SFA-EMT Single Rate Solution For this case the hybrid multirate protocol was used in single rate mode employing a time step of 0.05 ms for both the SFA and the EMT solutions. The purpose of presenting this case is to show how the hybrid SFA-EMT single rate solution matches the EMT only solution to the fastest transients. Results for the SFA SS are presented in Figure 6.21 and for the EMT SS in Figure 6.22 to Figure 6.26. Though 300 ms were simulated, for better resolution of the graphs, only the first 100 ms or less are shown. As it can be seen from these results, there is a very good match between the reference solution and the hybrid single rate one. Only at one single point the error reached 1% in the rectifier DC voltage, but in general, for the whole simulation and for all variables, the error is less than 0.5%. These results confirm the accuracy of the proposed interface model between an SFA solution and an EMT solution for systems with power electronics-based devices and in single rate mode. 104 Figure 6.21 Voltage at bus 1 (phase “b”) in the SFA subsystem (single rate solution). Figure 6.22 Current from AC system into bus 1 (phase “a”) in the EMT subsystem (single rate solution). 105 Figure 6.23 DC voltage at rectifier terminal in the EMT subsystem (single rate solution). Figure 6.24 DC voltage at inverter terminal in the EMT subsystem (single rate solution). 106 Figure 6.25 DC current at rectifier terminal in the EMT subsystem (single rate solution). Figure 6.26 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (single rate solution). 107 6.2.2 Case 2b: Hybrid SFA-EMT Multirate Solution As mentioned earlier, for this case the SFA subsystem was simulated using a time step of 0.49 ms (490 μs) and the EMT subsystem using a time step of 0.05 ms (50 μs), which corresponds to a rate ratio of 9.8 (not an integer!). The reason for choosing this low rate ratio was to have an adequate level of accuracy for the EMT solution during the transient period, since for this partitioned system, there is no buffer zone between SS1 and SS2 and the only variables of interest in the SFA SS are just the variables at the boundary with the EMT SS. The antialiasing processing method used for this case was a Savitzky-Golay low-pass filter with a polynomial of order 2. The standard MATE procedure was used. Results for the SFA SS are presented in Figure 6.27 where it can be observed how the magnitude of the dynamic phasor of the SFA solution accurately captures the envelope of the voltage time signal after the transient period, the error is 0.2 %. Note that for the SFA solution the interest is not on the electromagnetic transients of the system. The results for the EMT SS are shown in Figure 6.28 to Figure 6.32. From the graphs of the EMT solution it can be observed that it reproduces with good accuracy the transient period of the startup of the rectifier. The error in the peak current (phase “a”) that flows from the AC system is 2.5 % (see Figure 6.28). For the DC voltage at the rectifier the maximum error in the transient period is 5.4 % but at the peak is just 0.8 % (see Figure 6.29). The errors at the peak values of the inverter voltage, the current at the rectifier terminal and the current in valve V 1-2 are, respectively, 1.4 %, 3.5 % and 3.5 % (Figure 6.30 to Figure 6.32). 108 Figure 6.27 Voltage at bus 1 (phase “b”) in the SFA subsystem (multirate solution). Figure 6.28 Current from AC system into bus 1 (phase “b”) in the EMT subsystem (multirate solution). 109 Figure 6.29 DC voltage at rectifier terminal in the EMT subsystem (multirate solution). Figure 6.30 DC voltage at inverter terminal in the EMT subsystem (multirate solution). 110 Figure 6.31 DC current at rectifier terminal in the EMT subsystem (multirate solution). Figure 6.32 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (multirate solution). 111 These results validate the application of the proposed hybrid interface model between an SFA solution and an EMT solution for the multirate simulation of systems with power electronics-based components. They also show, however, that in order to obtain higher accuracy during the transient period and to enable the use of high rate ratios, the partition of the system must provide a significant buffer zone between the subsystems, which means that the boundary of the EMT subsystem should be placed at some distance from the source of transients. These results also show that the accuracy is excellent when the rectifier has reached normal operating conditions. This situation is examined in more detail in Case 2c below. 6.2.3 Case 2c: Hybrid SFA-EMT Multirate Solution in Normal Operation The purpose of this case was to assess the accuracy of the hybrid multirate protocol for the period after the startup process has ended and the rectifier station is operating normally. The conditions for this case were the same as for Case 2b, except that this time a much larger time step of 4.01 ms was used for the SFA solution. The time step for the EMT subsystem was kept in 0.05 ms; therefore, the multirate solution for this case had a rate ratio of 80.2. Results for the SFA SS are presented in Figure 6.33 and for the EMT SS in Figure 6.34 to Figure 6.38. As it can be observed in those graphs, the multirate solution matches very well the reference EMT solution during the steady state operation of the rectifier station. For these conditions, the maximum errors in the AC quantities were less than 1.25 % and less than 0.6 % in the DC quantities. This confirms the accuracy of the proposed hybrid multirate interfacing protocol with high rate ratios for normal operating conditions in systems with power electronics-based devices. 112 Figure 6.33 Voltage at bus 1 (phase “b”) in the SFA subsystem (steady state multirate solution). Figure 6.34 Current from AC system into bus 1 (phase “b”) in the EMT subsystem (steady state multirate solution). 113 Figure 6.35 DC voltage at rectifier terminal in the EMT subsystem (steady state multirate solution). Figure 6.36 DC voltage at inverter terminal in the EMT subsystem (steady state multirate solution). 114 Figure 6.37 DC current at rectifier terminal in the EMT subsystem (steady state multirate solution). Figure 6.38 Current in valves V 1-2, V 1-4 and V 1-6 in the EMT subsystem (steady state multirate solution). 115 Chapter 7: Conclusion The research project object of this dissertation was focused on developing a protocol to directly interface, in multirate mode, an SFA dynamic phasor solution with an EMT solution in an accurate and efficient way and to develop a simulator that operated with such protocol. As a result of this research, several contributions have been made and several areas of future work have been identified. In the following all these are presented. 7.1 Contributions Contribution 1: As part of this research a new hybrid protocol to interface a dynamic phasor solution based on shifted frequency analysis (SFA) and an electromagnetic transient (EMT) solution has been developed and validated. This protocol integrates both solutions in a multirate simulation framework and has the following characteristics:  It performs two simultaneous parallel EMT simulations of the EMT subsystem, one to interface with the real part and one to interface with the imaginary part of the complex-plane SFA solution. The solution for the real part is the actual solution of the system, while the solution for the imaginary part maintains the mathematical property of analyticity of the corresponding complex-plane solution in normal time domain and is, therefore, consistent with the use of large time steps in the SFA solution. This is a major contribution of the work presented in this dissertation, which has not been suggested prior to this research. 116  It allows for direct simultaneous multirate interfacing of hybrid SFA/EMT simulators within the MATE solution framework, without iterations and without a time delay to decouple the solutions. No decoupling elements, like transmission lines, are required either. With MATE, at each EMT time step there is available an interpolated value of the SFA solution, which increases the accuracy of the hybrid multirate solution and decreases the possibility of numerical instability.  In this interfacing method, the SFA subsystem is represented, using MATE, as a multi-port Thévenin equivalent that results from the discretization of the differential equations and therefore, can capture the frequency response of the SFA subsystem for a range of frequencies around the fundamental, which is determined by the time step used in the SFA solution. This allows the simulation of scenarios involving harmonic and sub-harmonic resonances and is also important in stability studies where large frequency excursions are present.  The proposed framework is very general in the multirate sense as it allows for arbitrary time steps in each subsystem, and it is not restricted to these time steps being integer multiples. This provides generality in choosing the most adequate time step size according to the needs of each subsystem, thus further increasing the flexibility and efficiency of the simulation.  It is based in a-b-c coordinates and therefore needs no assumptions regarding balanced operating conditions in the system. This protocol represents the first introduction of a new concept to interface dynamic phasors and EMT multirate systems in a direct and simultaneous combined solution. 117 The validity of the proposed multirate interfacing method has been demonstrated using the IEEE 39-bus system for a case of transient stability. The results show that the proposed method exhibits very good accuracy for capturing the slow electromechanical transients in the SFA solution and, vice-versa, for capturing the fast electromagnetic transients in the EMT solution. The multirate nature of the solution provided, for this case, modest computational savings when compared with a full-system EMT solution. This was due to the fact that, to significantly outperform an EMT solution, an optimized implementation of the SFA-EMT protocol is required with the two EMT solutions parallelized. Nevertheless, actual power systems are very large and, therefore, the part of the system that would be modeled with SFA would be a large percentage of the system and, with integration steps for the SFA solution in the range of milliseconds and in the range of microseconds for the EMT solution, the rate ratio of an SFA-EMT multirate solution would easily fall in the range of 100-500. Therefore, for actual power systems considerable computational gains can be expected with the proposed hybrid protocol. The validity of the proposed multirate interfacing method has also been demonstrated in a system with power electronics-based components. On this occasion a modified version of the CIGRE HVDC benchmark system was used. The purpose this time was solely to test the protocol in a system with components that exhibit fast switching/commutation events in normal operation, since under these conditions the two EMT solutions that the protocol uses may not have the same topology at all times. The results show that the proposed method exhibits good accuracy for capturing the transient behaviour of the EMT subsystem for the startup or initialization of the rectifier converter station, but only a modest rate ratio (ΔTSFA/ΔtEMT) of 9.8 was reached due to the fact that there was no buffer zone between the two subsystems. For the operation in normal 118 conditions of the converter rectifier station, however, the accuracy of the solution was very good and permitted the increase of the rate ratio up to 80.2 while maintaining good accuracy. For both, the IEEE and the CIGRE test systems, the single rate version of the protocol has been also validated, showing an accuracy as good as the full EMT solution. The results are almost identical. Contribution 2: In search for optimizing the accuracy of the multirate solution, in this research several methods of conditioning the fast EMT solutions to avoid aliasing in the multirate solution, also known as antialiasing methods, were explored, implemented and tested. Findings and tests performed are summarized below.  Tests of antialiasing methods were carried out along two different axes. On one axis, three main types of methods for smoothing the EMT solution were tested; namely, plain average or moving average, linear least-square (LLS) fitting and a Savitzky-Golay (S-G) low-pass filter using a polynomial of second order. On the other axis, two ways of carrying out the SFA-EMT multirate protocol were considered, the standard one, in which a Thévenin equivalent of the EMT subsystem is used in the updating the SFA subsystem, and an alternative implementation of the protocol that uses other methods to find the link solution to update the SFA subsystem. In the standard formulation of the protocol two additional averaging methods were considered; namely, averaging over a full period of the fundamental frequency and averaging over a full period of the fundamental frequency only the real part of the EMT solution. In the alternative 119 formulation, two methods to find the link solution for the SFA subsystem were evaluated, one method was founded on smoothing out the EMT link currents and the other was based on determining a link current to balance the power exchange between the two subsystems.  All smoothing methods developed here are unique in the sense that they are performed over a time period equal to the time step of the SFA solution and this makes them more general, because they are not restricted to the frequency of the EMT subsystem being equal to the fundamental.  All antialiasing methods developed for this research do not introduce delays in the smoothed signal, which means that they will not distort the phase of that signal and, therefore, they will not affect the accuracy of the multirate solution. This is a very important feature for the direct interface protocol proposed here.  For the system considered, the antialiasing processing methods that employ an S-G filter showed the best performance in terms of the accuracy of the SFA-EMT multirate solution. Also, it is important to point out that, for this hybrid protocol the antialiasing method based on plain averaging, though it does not introduce delays and its performance is not radically worse than the other two smoothing methods, S-G filters and, in second place LLS fitting are to be preferred over plain averaging because they offer the possibility of extracting the dynamic phasor from the EMT solution, while plain average will extract only an equivalent fundamental frequency phasor from the EMT solution which would be a limiting feature, especially in stability simulations of low inertia systems. 120  In the alternative formulation of the SFA-EMT multirate protocol, it could be confirmed that doing the antialiasing processing directly on the EMT link currents gives approximately the same performance as the standard formulation of the protocol which does the antialiasing on the EMT Thévenin voltages.  The incorporation of a power exchange balance constraint in the antialiasing processing proved to be a good and effective approach to carry out the antialiasing task, but the performance is approximately the same as working with the Thévenins in the standard formulation of the SFA-EMT protocol and, since it involves additional computations, it is not computationally attractive. However, this test proved that the proposed SFA-EMT MATE multirate protocol implicitly accounts for power exchange balance in the multirate solution, which constitutes a positive feature of this protocol and also adds to the strengths of the MATE solution framework.  Having a second EMT solution in the protocol allowed for an easy implementation of the antialiasing method with a power exchange balance constraint and also offered an opportunity to put into practice the novel concepts of Total Instantaneous Complex Power (TICP), Total Instantaneous Active Power (TIAP) and Total Instantaneous Reactive Power (TIRP) proposed in Section 2.4 of this dissertation for circuits with complex signals. These are useful and interesting concepts as the SFA transformation is TICP, TIAP and TIRP invariant. This means that the values of these instantaneous powers are equal in both, the normal time domain and the SFA domain, and this is independent of the state in which the circuit is on and it is also independent of the normal time domain signals being analytic or not. 121 Contribution 3: As a result of this research project, an interfacing software that integrates an SFA solver and an EMT solver and that operates with the proposed hybrid SFA-EMT multirate protocol was developed. The software is based on the MATE solution framework and uses the Python programming language. It can employ several antialiasing processing methods that result from the combination of the standard and alternative formulations of the protocol and three types of smoothing methods, namely, averaging, linear least-square fitting and Savitzky-Golay low-pass filters. The development also included expanding the modeling capabilities of the existing SFA and EMT solvers to incorporate: six-pulse diode bridges (EMT solver only), ideal transformers and three-phase coupled branches. 7.2 Future Work In the course of this research and during the processes of development of the simulator and the validation of the protocol, various areas of improvement and of future research were identified. Based on that, in the following some ideas for future work are proposed. a) Since SFA solves for the differential equations in the system, the SFA solution in a hybrid SFA-EMT simulator can capture not only the fundamental frequency in the SFA subsystem, but also a full range of frequencies around it, which would be larger or smaller depending on the SFA time step being used (see last paragraph of Section 4.3.1). Therefore, the proposed SFA-EMT multirate protocol lends itself naturally to the simulation of sub-harmonic and harmonic interactions between the EMT and the SFA subsystems. Thus, future research could be directed to explore this area of application of 122 the protocol to study the effect on the AC grid of the injection of harmonics and subharmonics from distributed energy resources. That would allow the confirmation of this type of application for the proposed protocol, together with the identification of the conditions and limitations associated to it and also any adaptations or enhancements that the protocol may require. b) The proposed hybrid interfacing scheme has three factors that limit its performance against an all EMT solution with regards to computational efficiency: a) it requires a link solver which has to carry out all link solutions in complex numbers, b) the SFA solution uses complex arithmetic and involves, for the same time step, roughly twice the number of operations as an EMT solution and c) the proposed interface requires two simultaneous EMT solutions. The present implementation of the protocol is serial; however, this could be improved by parallelizing the two EMT solutions. On the other hand, the code of the link solver of the simulator can be optimized to improve the computational efficiency of the multirate simulation. Additionally, the code of the SFA solver can be improved for faster processing; in particular, a different implementation of the solution algorithm for the SFA solver could be employed to render it more efficient, using for example a three-sequence based approach as suggested in [85] and [86]. All these actions would make the simulator more efficient, which added to its application to larger systems would permit more precise evaluations of the time savings that can be achieved with this hybrid multirate protocol. c) The SFA and EMT solvers that integrate the hybrid simulator are in their first stage of development and, in consequence, their modeling capabilities are limited at present which restricts the type of systems to which the interfacing protocol can be applied and 123 the operating conditions that can be simulated. To overcome this, the models library of these solvers should be enhanced by incorporating, for example: detailed models of generators and their controls, a variety of models of power electronics devices and their controls, detailed models of transformers with saturation, distributed parameters transmission lines, non-linear elements and different types of loads. d) Developing the modeling capabilities of the EMT solver as suggested above could take some time; so, a faster way to expose the hybrid multirate protocol to a wider range of systems and operating conditions would be to implement the protocol using an industry-grade EMT simulator instead of the current EMT solver. Using commercially available simulators with MATE is very difficult because the internal algorithms are not accessible to the user and MATE requires the solution of the subsystems to be found in two steps, first in open circuit and then injecting the link current. Nevertheless, in this research this possibility was explored and, from the preliminary analyses, it was concluded that this option could prove to be feasible by using a modified MATE formulation and embedding the protocol as an user-defined-model into a commercial EMT simulator. This would require an EMT time step delay in the hybrid multirate solution. Note that is not an SFA time step. At first sight this does not seem so inconvenient, but the impact of this delay in the accuracy of the solution would have to be properly assessed. Hence, further research on this direction would be necessary to thoroughly evaluate the feasibility of this option. e) With the SFA solution the interest is to exploit its capability to use large time steps; however, when there are generators in the SFA subsystem, large time steps in the order of several milliseconds may require iterations in the SFA solution to ensure convergence of the generator models. At present the simulator does not have this capability which limits 124 the maximum size of the SFA time step that can be used. On the other hand, incorporating this capability its challenging for the hybrid multirate interfacing protocol, because if this is not done properly it could deteriorate its computational performance; thus, further research is necessary to identify an effective and efficient way of incorporating this feature in the hybrid protocol. 125 References [1] W. D. Humpage and K.-P. Wong, “Electromagnetic transient analysis in EHV power networks,” Proceedings of the IEEE, vol. 70, no. 4, pp. 379–402, Mar. 1982, doi: 10.1109/PROC.1982.12315. [2] H. W. Dommel, “Techniques for analyzing electromagnetic transients,” IEEE Computer Applications in Power, vol. 10, no. 3, pp. 18–21, Jul. 1997, doi: 10.1109/67.595285. [3] H. W. Dommel and W. S. Meyer, “Computation of electromagnetic transients,” Proceedings of the IEEE, vol. 62, no. 7, pp. 983–993, Jul. 1974, doi: 10.1109/PROC.1974.9550. [4] Manitoba HVDC Research Centre and RTDS Technologies Inc., “PSCAD.” https://www.pscad.com/ (accessed Sep. 23, 2022). [5] EMTP Alliance, “EMTP®.” https://www.emtp.com/ (accessed Sep. 23, 2022). [6] ATP User’s Group, “ATP-EMTP.” http://www.emtp.org/ (accessed Sep. 23, 2022). [7] MicroTran Power System Analysis Corp., “MicroTran.” http://www.microtran.com/ (accessed Sep. 23, 2022). [8] P. Kundur et al., “Definition and classification of power system stability,” IEEE Transactions on Power Systems, vol. 19, no. 3, pp. 1387–1401, Aug. 2004, doi: 10.1109/TPWRS.2004.825981. [9] L. Grigsby Ed., Power System Stability and Control. Boca Raton: CRC Press, 2007. [10] Siemens AG, “PSS®E.” https://new.siemens.com/global/en/products/energy/energy-automation-and-smart-grid/pss-software/pss-e.html (accessed Sep. 24, 2022). 126 [11] DIgSILENT GmbH, “PowerFactory.” https://www.digsilent.de/en/ (accessed Sep. 24, 2022). [12] Tractebel Engineering, “EUROSTAG.” https://www.eurostag.be/ (accessed Sep. 24, 2022). [13] Powertech Labs Inc., “TSAT.” https://www.dsatools.com/ (accessed Sep. 24, 2022). [14] W. Zhongxi and Z. Xiaoxin, “Power system analysis software package (PSASP)-an integrated power system analysis tool,” in Proc. POWERCON ’98, International Conference on Power System Technology, Institute of Electrical and Electronics Engineers Inc., 1998, pp. 7–11. doi: 10.1109/ICPST.1998.728602. [15] M. D. Heffernan, K. S. Turner, J. Arrillaga, and C. P. Arnold, “Computation of a.c.-d.c. system disturbances — Parts I, II and III,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 11, pp. 4341–4363, Nov. 1981, doi: 10.1109/TPAS.1981.316825. [16] J. Reeve and R. Adapa, “A new approach to dynamic analysis of ac networks incorporating detailed modeling of dc systems. Parts I and II,” IEEE Transactions on Power Delivery, vol. 3, no. 4, pp. 2005–2019, Oct. 1988, doi: 10.1109/61.194011. [17] H. W. Dommel, “Digital computer solution of electromagnetic transients in single- and multiphase networks,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-88, no. 4, pp. 388–399, Apr. 1969, doi: 10.1109/TPAS.1969.292459. [18] V. Jalili-Marandi, V. Dinavahi, K. Strunz, J. A. Martinez, and A. Ramirez, “Interfacing techniques for transient stability and electromagnetic transient programs - IEEE task force on Interfacing Techniques for Simulation Tools,” IEEE Transactions on Power Delivery, vol. 24, no. 4, pp. 2385–2395, Oct. 2009, doi: 10.1109/TPWRD.2008.2002889. 127 [19] G. W. Anderson, C. P. Arnold, N. R. Watson, and J. Arrillaga, “A new hybrid AC-DC transient stability program,” in Proc. International Conference on Power System Transients (IPST 95), Lisbon, Portugal, Sep. 1995, pp. 535–540. [20] M. Sultan, J. Reeve, and R. Adapa, “Combined transient and dynamic analysis of HVDC and FACTS systems,” IEEE Transactions on Power Delivery, vol. 13, no. 4, pp. 1271–1277, Oct. 1998, doi: 10.1109/61.714495. [21] Y. Zhang, A. M. Gole, W. Wu, B. Zhang, and H. Sun, “Development and analysis of applicability of a hybrid transient simulation platform combining TSA and EMT elements,” IEEE Transactions on Power Systems, vol. 28, no. 1, pp. 357–366, Feb. 2013, doi: 10.1109/TPWRS.2012.2196450. [22] Xi Lin, A. M. Gole, and Ming Yu, “A wide-band multi-port system equivalent for real-time digital power system simulators,” IEEE Transactions on Power Systems, vol. 24, no. 1, pp. 237–249, Feb. 2009, doi: 10.1109/TPWRS.2008.2007000. [23] X. Zhang, A. J. Flueck, and S. Abhyankar, “Implicitly coupled electromechanical and electromagnetic transient analysis using a frequency-dependent network equivalent,” IEEE Transactions on Power Delivery, vol. 32, no. 3, pp. 1262–1269, Jun. 2017, doi: 10.1109/TPWRD.2015.2501387. [24] K. H. Su, W. Chan, L. A. Snider, and J. -C. Soumagne, “Advancements on the integration of electromagnetic transients simulator and transient simulator,” in Proc. International Conference on Power System Transients (IPST 05), Jun. 2005, pp. 1–6. [25] H. T. Su, K. W. Chan, and L. A. Snider, “Parallel interaction protocol for electromagnetic and electromechanical hybrid simulation,” IEE Proceedings - Generation, Transmission and Distribution, vol. 152, no. 3, p. 406, 2005, doi: 10.1049/ip-gtd:20049019. 128 [26] P. Zadkhast, “Practical challenges of hybrid simulations interfacing with real-time digital simulators,” in Proc. International Conference on Power Systems Transients, Jun. 2019, pp. 1–6. [27] S. Abhyankar, “Development of an implicitly coupled electromechanical and electromagnetic transients simulator for power systems,” PhD Dissertation, Illinois Institute of Technology, Chicago, Illinois, U.S.A., 2011. [28] F. Plumier, “Co-simulation of electromagnetic transients and phasor models of electric power systems,” PhD Dissertation, Institut Montefiore, Liege, Belgium, 2015. [29] F. Plumier, P. Aristidou, C. Geuzaine, and T. van Cutsem, “Co-simulation of electromagnetic transients and phasor models: A relaxation approach,” IEEE Transactions on Power Delivery, vol. 31, no. 5, pp. 2360–2369, Oct. 2016, doi: 10.1109/TPWRD.2016.2537927. [30] A. A. van der Meer, R. L. Hendriks, M. Gibescu, and W. L. Kling, “Interfacing methods for combined stability and electromagnetic transient simulations applied to VSC-HVDC,” in Proc. International Conference on Power System Transients (IPST 11), Delft, The Netherlands, 2011, pp. 1–9. [31] A. A. van der Meer, M. Gibescu, M. A. van der Meijden, W. L. Kling, and J. A. Ferreira, “Advanced Hybrid Transient Stability and EMT Simulation for VSC-HVDC Systems,” IEEE Transactions on Power Delivery, vol. 30, no. 3, pp. 1057–1066, Jun. 2015, doi: 10.1109/TPWRD.2014.2384499. [32] Q. Huang and V. Vittal, “Application of electromagnetic transient-transient stability hybrid simulation to FIDVR study,” IEEE Transactions on Power Systems, vol. 31, no. 4, pp. 2634–2646, Jul. 2016, doi: 10.1109/TPWRS.2015.2479588. 129 [33] S. R. Sanders, J. M. Noworolski, X. Z. Liu, and G. C. Verghese, “Generalized averaging method for power conversion circuits,” IEEE Trans Power Electron, vol. 6, no. 2, pp. 251–259, Apr. 1991, doi: 10.1109/63.76811. [34] A. M. Stankovic and T. Aydin, “Analysis of asymmetrical faults in power systems using dynamic phasors,” IEEE Transactions on Power Systems, vol. 15, no. 3, pp. 1062–1068, 2000, doi: 10.1109/59.871734. [35] A. M. Stankovic, B. C. Lesieutre, and T. Aydin, “Modeling and analysis of single-phase induction machines with dynamic phasors,” IEEE Transactions on Power Systems, vol. 14, no. 1, pp. 9–14, Feb. 1999, doi: 10.1109/59.744460. [36] T. Demiray, G. Andersson, and L. Busarello, “Evaluation study for the simulation of power system transients using dynamic phasor models,” in 2008 IEEE/PES Transmission and Distribution Conference and Exposition: Latin America, IEEE, Aug. 2008, pp. 1–6. doi: 10.1109/TDC-LA.2008.4641716. [37] E. Zhijun, K. W. Chan, and D. Z. Fang, “Hybrid simulation of power systems with dynamic phasor SVC transient model,” in 2007 7th International Conference on Power Electronics and Drive Systems, IEEE, Nov. 2007, pp. 1670–1675. doi: 10.1109/PEDS.2007.4487933. [38] T. H. Demiray, “Simulation of power system dynamics using dynamic phasor models,” PhD Dissertation, Swiss Federal Institute of Technology, Zurich, Switzerland, 2008. [39] M. Kulasza, “Generalized dynamic phasor-based simulation for power systems,” MSc Dissertation, University of Manitoba, Winnipeg, MB, Canada, 2014. [40] P. Zhang, J. R. Marti, and H. W. Dommel, “Shifted-frequency analysis for EMTP simulation of power-system dynamics,” IEEE Transactions on Circuits and Systems I: 130 Regular Papers, vol. 57, no. 9, pp. 2564–2574, Sep. 2010, doi: 10.1109/TCSI.2010.2043992. [41] J. Marti, “Shifted frequency analysis (SFA) for EMTP simulation of fundamental frequency dynamics,” Power Systems Lab., Internal Report, University of British Columbia, Vancouver, BC, Canada, 2005. [42] J. R. Marti, H. W. Dommel, B. D. Bonatto, and A. F. R. Barreto, “Shifted frequency analysis (SFA) concepts for EMTP modelling and simulation of power system dynamics,” in 2014 Power Systems Computation Conference, IEEE, Aug. 2014, pp. 1–8. doi: 10.1109/PSCC.2014.7038487. [43] A. T. J. Marti, “The application of shifted frequency analysis in power system transient stability studies,” M.A.Sc. Dissertation, Univ. British Columbia, Vancouver, BC, Canada, 2018. [44] P. Zhang, J. R. Marti, and H. W. Dommel, “Synchronous machine modeling based on shifted frequency analysis,” IEEE Transactions on Power Systems, vol. 22, no. 3, pp. 1139–1147, Aug. 2007, doi: 10.1109/TPWRS.2007.901288. [45] Y. Huang, M. Chapariha, F. Therrien, J. Jatskevich, and J. R. Marti, “A constant-parameter voltage-behind-reactance synchronous machine model based on shifted-frequency analysis,” IEEE Transactions on Energy Conversion, vol. 30, no. 2, pp. 761–771, Jun. 2015, doi: 10.1109/TEC.2014.2367098. [46] Peng Zhang, J. R. Marti, and H. W. Dommel, “Induction machine modeling based on shifted frequency analysis,” IEEE Transactions on Power Systems, vol. 24, no. 1, pp. 157–164, Feb. 2009, doi: 10.1109/TPWRS.2008.2005715. 131 [47] P. Zhang, “Shifted frequency analysis for EMTP simulation of power system dynamics,” PhD Dissertation, University of British Columbia, Vancouver, BC, Canada, 2009. [48] D. Shu, V. Dinavahi, X. Xie, and Q. Jiang, “Shifted frequency modeling of hybrid modular multilevel converters for simulation of MTDC grid,” IEEE Transactions on Power Delivery, vol. 33, no. 3, pp. 1288–1298, Jun. 2018, doi: 10.1109/TPWRD.2017.2749427. [49] K. Mudunkotuwa, S. Filizadeh, and U. Annakkage, “Development of a hybrid simulator by interfacing dynamic phasors with electromagnetic transient simulation,” IET Generation, Transmission & Distribution, vol. 11, no. 12, pp. 2991–3001, Aug. 2017, doi: 10.1049/iet-gtd.2016.1616. [50] Y. Li, D. Shu, F. Shi, Z. Yan, Y. Zhu, and N. Tai, “A multi-rate co-simulation of combined phasor-domain and time-domain models for large-scale wind farms,” IEEE Transactions on Energy Conversion, vol. 35, no. 1, pp. 324–335, Mar. 2020, doi: 10.1109/TEC.2019.2936574. [51] D. Shu, X. Xie, V. Dinavahi, C. Zhang, X. Ye, and Q. Jiang, “Dynamic phasor based interface model for EMT and transient stability hybrid simulations,” IEEE Transactions on Power Systems, vol. 33, no. 4, pp. 3930–3939, Jul. 2018, doi: 10.1109/TPWRS.2017.2766269. [52] H. Konara, “Interfacing a transient stability model to a real-time electromagnetic transient simulation using dynamic phasors,” PhD Dissertation, The University of Manitoba, Winnipeg, MB, Canada, 2019. [53] D. Shu, Z. Ouyang, and Z. Yan, “Multirate and mixed solver based cosimulation of combined transient stability, shifted-frequency phasor, and electromagnetic models: A 132 Practical LCC HVDC Simulation Study,” IEEE Transactions on Industrial Electronics, vol. 68, no. 6, pp. 4954–4965, Jun. 2021, doi: 10.1109/TIE.2020.2988228. [54] F. Gao and K. Strunz, “Frequency-adaptive power system modeling for multiscale simulation of transients,” IEEE Transactions on Power Systems, vol. 24, no. 2, pp. 561–571, May 2009, doi: 10.1109/TPWRS.2009.2016587. [55] J. Rupasinghe, S. Filizadeh, A. M. Gole, and K. Strunz, “Multi‐rate co‐simulation of power system transients using dynamic phasor and EMT solvers,” The Journal of Engineering, vol. 2020, no. 10, pp. 854–862, Oct. 2020, doi: 10.1049/joe.2020.0099. [56] V. N. Sewdien et al., “Effects of increasing power electronics on system stability: results from MIGRATE questionnaire,” in 2018 International Conference and Utility Exhibition on Green Energy for Sustainable Development (ICUE), IEEE, Oct. 2018, pp. 1–9. doi: 10.23919/ICUE-GESD.2018.8635602. [57] G. Nakti, J. Dinkelbach, M. Mirz, and A. Monti, “Comparative assessment of shifted frequency modeling in transient stability analysis using the open source simulator DPsim,” in 2022 Open Source Modelling and Simulation of Energy Systems (OSMSES), IEEE, Apr. 2022, pp. 1–6. doi: 10.1109/OSMSES54027.2022.9769135. [58] J. Dinkelbach, G. Nakti, M. Mirz, and A. Monti, “Simulation of low inertia power systems based on shifted frequency analysis,” Energies (Basel), vol. 14, no. 7, p. 1860, Mar. 2021, doi: 10.3390/en14071860. [59] A. T. J. Marti and J. Jatskevich, “Transient Stability Analysis Using Shifted Frequency Analysis (SFA),” in 2018 Power Systems Computation Conference (PSCC), IEEE, Jun. 2018, pp. 1–7. doi: 10.23919/PSCC.2018.8442806. 133 [60] F. A. Moreira, J. R. Marti, L. C. Zanetta, and L. R. Linares, “Multirate simulations with simultaneous-solution using direct integration methods in a partitioned network environment,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 53, no. 12, pp. 2765–2778, Dec. 2006, doi: 10.1109/TCSI.2006.882821. [61] V. A. Galván-Sánchez, J. R. Martí, E. S. Bañuelos-Cabral, J. Sotelo-Castañón, J. L. García-Sánchez, and J. A. Gutiérrez-Robles, “An asynchronous MATE-multirate method for the modeling of electric power systems,” Electrical Engineering, vol. 103, no. 2, pp. 993–1007, Apr. 2021, doi: 10.1007/s00202-020-01128-x. [62] H. Hsu, Signals and Systems. New York: McGraw-Hill, 1995. [63] S. Henschel, “Analysis of electromagnetic and electromechanical power system transients with dynamic phasors,” PhD Dissertation, University of British Columbia, Vancouver, BC, Canada, 1999. [64] P. Denbigh, System Analysis & Signal Processing. Harlow, England: Addison-Wesley, 1998. [65] S. L. Hahn, Hilbert Transforms in Signal Processing. Boston: Artech House, 1996. [66] D. Gabor, “Theory of communication. Part 1: The analysis of information,” Journal of the Institution of Electrical Engineers - Part III: Radio and Communication Engineering, vol. 93, no. 26, pp. 429–441, Nov. 1946, doi: 10.1049/ji-3-2.1946.0074. [67] R. Bracewell, The Fourier Transform and its Applications, Second Edition. New York: McGraw-Hill, Inc., 1986. [68] L. Cohen, Time-Frequency Analysis. New York: Prentice-Hall, 1995. [69] H. Dommel, EMTP Theory Book, Second editon. Vancouver, BC, Canada: Microtran Power Systems Analysis Corporation, 1986. 134 [70] V. Venkatasubramanian, “Tools for dynamic analysis of the general large power system using time-varying phasors,” International Journal of Electrical Power & Energy Systems, vol. 16, no. 6, pp. 365–376, 1994, doi: https://doi.org/10.1016/0142-0615(94)90023-X. [71] P. Lynn and W. Fuerst, Introductory Digital Signal Processing with Computer Applications, Second edition. New York: John Wiley & Sons, 1998. [72] R. W. Hamming, Digitial Filters, Third Edition. New Jersey: Prentice Hall, 1989. [73] L. R. Linares-Rojas, “OVNI (Object Virtual Network Integrator) : a new fast algorithm for the simulation of very large electric networks in real time,” University of British Columbia, Vancouver, 2000. [74] J. R. Marti, L. R. Linares, J. A. Hollman, and Moreira F.A., “OVNI: Integrated software/hardware solution for real-time simulation of large power systems,” in 14th Power Systems Computation Conference (PSSC), Sevilla, Spain, Jun. 2002. [75] J. R. Marti, L. R. Linares, J. Calvino, H. W. Dommel, and J. Lin, “OVNI: an object approach to real-time power system simulators,” in POWERCON ’98. 1998 International Conference on Power System Technology. Proceedings (Cat. No.98EX151), Beijing, China: IEEE, Aug. 1998, pp. 977–981. doi: 10.1109/ICPST.1998.729230. [76] V. A. Galvan, J. R. Marti, H. W. Dommel, J. A. Gutierrez-Robles, and J. L. Naredo, “MATE multirate modelling of power electronic converters with mixed integration rules,” in 2016 Power Systems Computation Conference (PSCC), Genoa, Italy: IEEE, Jun. 2016, pp. 1–7. doi: 10.1109/PSCC.2016.7540945. 135 [77] F. A. Moreira and J. R. Marti, “Latency techniques for time-domain power system transients simulation,” IEEE Transactions on Power Systems, vol. 20, no. 1, pp. 246–253, Feb. 2005, doi: 10.1109/TPWRS.2004.841223. [78] J. O. Tarazona, A. T. J. Martí, J. R. Martí, and F. A. Moreira, “Shifted frequency analysis-EMTP multirate simulation of power systems,” Electric Power Systems Research, vol. 197, p. 107292, 2021, doi: https://doi.org/10.1016/j.epsr.2021.107292. [79] S. D. Pekarek et al., “An efficient multirate simulation technique for power-electronic-based systems,” IEEE Transactions on Power Systems, vol. 19, no. 1, pp. 399–409, Feb. 2004, doi: 10.1109/TPWRS.2003.821452. [80] R. W. Schafer, “On the frequency-domain properties of Savitzky-Golay filters,” in 2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE), Palo Alto, CA: IEEE, Jan. 2011, pp. 54–59. doi: 10.1109/DSP-SPE.2011.5739186. [81] D. Shu, X. Xie, Q. Jiang, Q. Huang, and C. Zhang, “A novel interfacing technique for distributed hybrid simulations combining EMT and transient stability models,” IEEE Transactions on Power Delivery, vol. 33, no. 1, pp. 130–140, Feb. 2018, doi: 10.1109/TPWRD.2017.2690145. [82] J. R. Marti and J. Lin, “Suppression of numerical oscillations in the EMTP power systems,” IEEE Transactions on Power Systems, vol. 4, no. 2, pp. 739–747, May 1989, doi: 10.1109/59.193849. [83] P. Kundur, Power System Stability and Control. San Francisco, CA, USA: McGraw-Hill, Inc. [84] IEEE PES Task force on benchmark systems for stability controls, “Benchmark systems for small-signal stability analysis and control, Tech. Rep. PES-TR18,” Aug. 2015. 136 [85] Y. Song, Y. Chen, Z. Yu, S. Huang, and C. Shen, “Efficient shifted frequency analysis for power system in symmetrical component domain,” in International Conference on Power Systems Transients (IPST2019), Perpignan, France, Jun. 2019, pp. 1–6. [86] S. Gao, Y. Song, Y. Chen, Z. Yu, and Z. Tan, “Shifted frequency‐based electromagnetic transient simulation for AC power systems in symmetrical component domain,” IET Renewable Power Generation, Feb. 2022, doi: 10.1049/rpg2.12418. [87] M. Szechtman, T. Wess, and C. V. Thio, “First benchmark model for HVDC control studies,” Electra, vol. 135, pp. 56–73, Apr. 1991. 137 Appendix Data of the IEEE 39-Bus Test System Tables A1 and A2 below present the data used for the IEEE 39-bus system, that has been modified from [84]. The rest of the data, transmission lines and transformers, is as indicated in [84]. Table A1 Generator data for the IEEE 39-bus system used in Case Study 1. 1 1000.00 100.0 0.3 6 1.52 677.87 16.5 0.3 5 1.53 650.00 16.5 0.3 5 1.54 632.00 16.5 0.3 5 1.55 508.00 16.5 0.3 5 1.56 650.00 16.5 0.3 5 1.57 560.00 16.5 0.3 5 1.58 540.00 16.5 0.3 5 1.59 830.00 16.5 0.3 5 1.510 250.00 16.5 0.3 5 1.5(a): per unit of the generator base (b): per unit of torque baseGenerator #Rated Power [MVA]Rated Voltage [kV]X'd [p.u.a]Inertia Constant (H) [s]Damping Coefficient (K) [p.u.b]138 Table A2 Power flow data for the IEEE 39-bus system used for the initialization of Case Study 1 Bus V Angle P Q P Q# [pu] [Deg] [MW] [MVAR] [MW] [MVAR]3 1.0273 257.60 1.904 1.0118 400.00 147.207 1.0074 187.00 67.208 1.0062 417.60 141.3012 1.0125 6.00 70.4015 1.0201 256.00 122.4016 1.0315 263.20 25.8018 1.0300 126.40 24.0020 0.9953 502.40 82.4021 1.0260 219.20 92.0023 1.0278 198.00 67.6824 1.0344 246.88 ‐73.7625 1.0444 179.20 37.7626 1.0487 111.20 13.6027 1.0373 224.80 60.4028 1.0470 164.80 22.0829 1.0449 226.80 21.5230 1.0100 ‐4.0704 200.00 8.3331 0.9820 0 488.34 160.73 7.36 3.6832 0.9841 0.7238 520.00 157.9133 0.9972 1.1941 505.60 90.2634 1.0123 0.3701 406.40 139.5435 1.0200 1.6891 480.00 102.7836 1.0250 3.7142 408.00 3.5437 1.0150 0.4308 432.00 ‐7.1538 1.0150 4.6056 664.00 ‐29.3939 1.0150 ‐9.3307 800.00 ‐3.20 883.20 200.00Generation Load
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Martí, José R.
Hybrid shifted frequency analysis-electromagnetic transients multirate simulation of power systems
Text