@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Electrical and Computer Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Moreira, Fernando Augusto"@en ; dcterms:issued "2009-10-06T21:41:35Z"@en, "2002"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """This thesis work presents a methodology for exploiting time domain latency in EMTP-type programs for the simulation of electromagnetic transients phenomena. Latency exploitation is related to the capability of numerically solving the differential equations governing the behaviour of electric networks with different integration steps. In many situations, the efficiency of numerical simulation is compromised by the requirements of the fastest network components that force the choice of a very small integration step to accurately track the highest frequency transients developing on the network. Through an efficient and accurate network partitioning technique, the electric network is divided into distinct subsystems that present different dominant eigenvalues. The subnetworks can then be simulated according to their own time step requirements. The effect of the subnetwork couplings are taken into account by resynchronizing the global solution at certain time instants. For this reason, the ratio between the large time step AT and the small time step At is always an integer number. The thesis includes a detailed description of how the interface between the fast and slow subnetworks is adjusted, in order to guarantee maximum accuracy and efficiency. The latency approach followed in this thesis work results in an expanded form of the trapezoidal and backward Euler integration rules, which are then tested for accuracy and stability. The analyses performed show that the expanded rules are as accurate and stable as the conventional rules for single time step simulation. Lumped networks and networks containing transmission lines, modelled with the travelling wave equations, are simulated with the proposed latency technique. Different line-models are implemented and the results compared to those obtained from conventional EMTP simulations, where a single integration step is used for the complete network solution."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/13663?expand=metadata"@en ; dcterms:extent "5096497 bytes"@en ; dc:format "application/pdf"@en ; skos:note "LATENCY TECHNIQUES IN POWER SYSTEM TRANSIENTS SIMULATION by Fernando Augusto Moreira Electrical Engineer, Universidade de Sao Paulo, Brazil, 1994 M . A . S c , Electrical Engineering, Universidade de Sao Paulo, Brazil, 1997 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF D O C T O R OF P H I L O S O P H Y in THE F A C U L T Y OF G R A D U A T E STUDIES (Department of Electrical & Computer Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 2002 © Fernando Augusto Moreira, 2002 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of E L R T K \\ ^ \\ f\\\\ft. W E The University of British Columbia Vancouver, Canada Date ^ T ^ v A DE-6 (2/88) A B S T R A C T This thesis work presents a methodology for exploiting time domain latency in EMTP-type programs for the simulation of electromagnetic transients phenomena. Latency exploitation is related to the capability of numerically solving the differential equations governing the behaviour of electric networks with different integration steps. In many situations, the efficiency of numerical simulation is compromised by the requirements of the fastest network components that force the choice of a very small integration step to accurately track the highest frequency transients developing on the network. Through an efficient and accurate network partitioning technique, the electric network is divided into distinct subsystems that present different dominant eigenvalues. The subnetworks can then be simulated according to their own time step requirements. The effect of the subnetwork couplings are taken into account by resynchronizing the global solution at certain time instants. For this reason, the ratio between the large time step AT and the small time step At is always an integer number. The thesis includes a detailed description of how the interface between the fast and slow subnetworks is adjusted, in order to guarantee maximum accuracy and efficiency. The latency approach followed in this thesis work results in an expanded form of the trapezoidal and backward Euler integration rules, which are then tested for accuracy and stability. The analyses performed show that the expanded rules are as accurate and stable as the conventional rules for single time step simulation. Lumped networks and networks containing transmission lines, modelled with the travelling wave equations, are simulated with the proposed latency technique. Different line-models are implemented and the results compared to those obtained from conventional EMTP simulations, where a single integration step is used for the complete network solution. ii Portuguese version: Este trabalho de tese apresenta uma metodologia para explorar latencia no dominio do tempo, em programas do tipo EMTP, para a simulacao de transitorios eletromagneticos. A exploracao da latencia esta relacionada com a capacidade de resolver numericamente as equacoes diferenciais que governam o comportamento de redes eletricas com diferentes passos de integracao. Em muitos casos, a eficiencia da simulacao numerica e comprometida pelos requisitos dos componentes mais rapidos da rede que forcam a escolha de um passo de integracao muito pequeno para a reproducao precisa dos transitorios de frequencia mais altas que se desenvolvem na rede. Atraves de uma tecnica de particao de rede eficiente e precisa, a rede eletrica e dividida em subsistemas distintos que apresentam autovalores dominantes diferentes. As subredes podem entao ser simuladas de acordo com seus proprios requisitos de passos de integracao. Os efeitos dos acoplamentos entre as subredes sao levados em consideracao atraves da re-sincronizacao da solucao global em certos instantes de tempo. Por esse motivo, a razao entre o passo de integracao longo AT e o passo de integracao curto At e sempre um numero inteiro. Esta tese inclui a descricao em detalhes de como a interface entre as subredes rapida e lenta e ajustada, de forma a garantir precisao e eficiencia maximas. A tecnica de latencia seguida neste trabalho de tese resulta em uma forma expandida das regras de integracao trapezoidal e \"backward Euler\", que sao entao testadas para precisao e estabilidade. As analises realizadas mostram que as regras expandidas sao tao estaveis e precisas quanto as regras convencionais para simulacao com um unico passo de integracao. Redes com parametros concentrados e redes contendo linhas de transmissao, modeladas com as equacoes de propagacao de ondas, sao simuladas com a tecnica de latencia proposta. Diferentes modelos de linha sao implementados e os resultados sao comparados com aqueles obtidos de simulacoes convencionais com o EMTP, nas quais um unico passo de integracao e usado para a solucao da rede inteira. iii C O N T E N T S Abstract ii Abstract - Portuguese version i i i Contents iv List of Tables vi List of Figures vii Acknowledgements x Dedication xii 1. Introduction and Overview 1 1.1. Background 1 1.2. Thesis Motivation 3 1.3. Thesis Organization 4 2. Literature Review 6 2.1. Network Partitioning 6 2.2. The Origins of Latency Exploitation 7 2.3. The Waveform Relaxation Method 8 2.4. The Variable Time Step Method 10 2.5. Exploiting Latency Using Implicit Integration Rules 12 2.6. Coupling Electromagnetic and Transient Stability Solution Methods 15 2.7. Latency Exploitation for Power Electronic Systems 16 3. Time-Domain Simulation of Electromagnetic Transients: A Brief Overview 18 3.1. Network Discretization in Time-Domain 18 3.1.1. The Trapezoidal Integration Rule 19 3.1.2. The Backward Euler Integration Rule 21 3.2. Accuracy of Integration Rules 24 3.2.1. Trapezoidal Integration Rule 24 3.2.2. Backward Euler Integration Rule 27 3.3. Stability of Integration Rules 30 3.3.1. Trapezoidal Integration Rule 30 3.3.2. Backward Euler Integration Rule 31 iv 4. Proposed Methodology for Latency Exploitation and its Main Characteristics 32 4.1. The Latency Simulation Method 32 4.1.1. History Source Accumulation for Resynchronized Solution 35 4.1.2. Thevenin Equivalent Interpolation for Fast-Varying Subnetwork Solution 39 4.2. Stability and Accuracy of the Proposed Latency Technique 44 4.2.1. Accuracy of the Expanded Integration Rules 46 4.2.2. Stability of the Expanded Integration Rules 61 4.3. Eigenanalysis for Latency Suitability 63 5. Case Studies 74 5.1. Latency Exploitation in Networks Consisting Solely of Lumped Elements 74 5.1.1. Simulation A 74 5.1.2. Simulation B 81 5.1.3. Simulation C 88 5.2. Latency Exploitation in Networks Containing Transmission Lines Modelled as Ideal Lossless Lines 90 5.3. Latency Exploitation in Networks Containing Transmission Lines Modelled with the cp-Line Model 94 6. Conclusions and Future Work 100 6.1. Conclusions 100 6.2. Recommendations for Future Work 102 Bibliography .....105 Appendix A 109 Appendix B 114 Appendix C 116 LIST O F T A B L E S 5.1 Number of FLOP S for the test case A considering a total simulation time of 1 Oms 81 5.2 Number of FLOPS for the test case B considering a total simulation time of 1.0ms 84 5.3 Number of FLOPS for the lossless line model considering a total simulation time of 1.5ms 93 5.4 Number of FLOPS for the cp-line model considering a total simulation time of 5ms 99 v i LIST OF FIGURES 2.1 Trapezoidal integration rule 12 2.2 Backward Euler integration rule 13 3.1 Inductance represented in continuous-time and its equivalent representation in discrete-time with the trapezoidal rule 20 3.2 Capacitance represented in continuous-time and its equivalent representation in discrete-time with the trapezoidal rule 21 3.3 Inductance represented in continuous-time and its equivalent representation in discrete-time with the backward Euler rule 22 3.4 Capacitance represented in continuous-time and its equivalent representation in discrete-time with the backward Euler rule 23 3.5 Frequency response in amplitude and phase of the trapezoidal integration rule 26 3.6 Frequency response in amplitude and phase of the backward Euler integration rule 29 4.1 Schematic diagram of the network decoupling 33 4.2 Timeline of the simulation with the large integration step AT containing several small integration steps At 34 4.3 Expanded trapezoidal rule for history source accumulation 36 4.4 5 Expanded backward Euler rule for history source accumulation 38 4.5 Network decoupling with the detailed slow subcircuit replaced by its Thevenin equivalent 40 4.6 Thevenin equivalent representation of slow and fast subnetworks at linking nodes 41 4.7 Thevenin equivalent history source linear interpolation for AT = 4At 44 4.8 Frequency domain representation of inductance and capacitance using the expanded trapezoidal integration rule for latency exploitation 49 4.9 Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 2At (partial and complete responses) 50 4.10 Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 2At (complete response) 51 vii 4.11 Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 4At (partial and complete responses) 52 4.12 Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 4At (complete response) 54 4.13 Frequency domain representation of inductance and capacitance using the expanded backward Euler integration rule for latency exploitation 58 4.14 Frequency response in magnitude and phase of the expanded backward Euler rule for AT = 4At (partial and complete responses) 59 4.15 Frequency response in magnitude and phase of the expanded backward Euler rule for AT = 4At (complete response) 60 4.16 Network studied for latency suitability 64 4.17 Modified network studied for latency suitability 68 5.1 Lumped circuit for latency exploitation (case A) 75 5.2 Voltage across the inductor for the three methods proposed 76 5.3 Current through the inductor for the three methods proposed 76 5.4 Voltage across the inductor up to 2ms and including the initialized value 78 5.5 Voltage across the capacitor for the three methods proposed 79 5.6 Current through the capacitor for the three methods proposed 79 5.7 Current through the capacitor up to 4ms and including the initialized value 80 5.8 Lumped circuit for latency exploitation (case B) 81 5.9 Voltage across capacitor C. (\"slow\") 82 5.10 Voltage across capacitor C2 (\"fast\") 83 5.11 Voltage across capacitor C. (\"slow\") for R = O.OOin 85 5.12 Voltage across capacitor C 2 (\"fast\") for R = 0.001Q 85 5.13 Voltage across capacitor C i (\"slow\") with the backward Euler rule 86 5.14 Voltage across capacitor C2 (\"fast\") with the backward Euler rule 87 5.15 Lumped circuit for latency exploitation (case C) (modified from case B) 88 5.16 Voltage across inductor L3 89 5.17 Network for latency exploitation in single-phase transmission lines 90 5.18 Voltage across the inductor located on the source side of the network 91 5.19 Voltage across the capacitor 92 viii 5.20 Circuit proposed for latency testing with the cp-line model 95 5.21 Transmission line divided into five segments for cp-line modelling 96 5.22 Voltage to the right of the first short transmission line segment 97 5.23 Voltage across the 350mH inductors in the external subnetwork 97 5.24 Voltage to the right of the first short transmission line segment 98 A. 1 M A T E example network 109 C l Separation of basic effects in the z-line model 118 ix A C K N O W L E D G E M E N T S Many people have contributed to the accomplishment of this Ph.D. thesis work. It is without doubt the biggest achievement in my professional career so far. My sincere apologies to anyone I may miss. A special gratitude to my supervisor, Professor lose R. Marti, for his dedicated guidance during the development of this thesis. His excellent teaching skills, combined with the unconditional support of his graduate students and their respective research have been fundamental for the successful conclusion of this work. Prof. Marti follows the work of his students very closely and with great enthusiasm. I am sure that having him as my supervisor was the right choice. Special thanks to Professor Hermann W. Dommel for his sharing of knowledge and experience. He has always been willing to contribute with his wisdom whenever students require some orientation. My gratitude and praise also to Prof. Dommel for opening the doors to the exciting field of digital simulation of electromagnetic transients in power systems when he created the EMTP. I am also thankful to Dr. Luis Linares, with whom I shared many hours of research and joy. The quality of his work and his broad knowledge were an inspiration for me. Special thanks to my friends Benedito D. Bonatto, Richard Rivas, and Antonio Carlos S. de Lima for their interest and discussion of ideas for my thesis. Many thanks also to them and to Luciana A . B. Bonatto and Martha R. P. Sanabria for their companionship during this time in Vancouver. Thanks to the Power Group members of U B C who were and still are my friends. They have always treated me with the respect and politeness expected from such a world-class institution. I wish them all the best on their research and future endeavours. Big thanks also to my friends Luiz Alexandre P. de Freitas and Fulvia de Freitas for their kindness with me and the happy times we spent together. I gratefully acknowledge the financial support of \"Fundacao Coordenacao de Aperfeicoamento de Pessoal de Nivel Superior\" (CAPES) from Brazil. In particular I am very thankful to Professor Sandoval Carneiro Jr. from the Federal University of Rio de Janeiro, who always demonstrated great interest in my research and was very supportive of it. I am also very grateful to Professor Luiz Cera Zanetta Jr. from the University of Sao Paulo. It was with him that I started to study the \"scary\", at first sight, field of electromagnetic transients in power systems. I met him for the first time in June, 1993 because I was interested in a position he had announced for an undergraduate research assistant. I went on to present my graduation paper having Prof. Zanetta as my advisor in 1994 and my M.A.Sc. thesis in 1997. He very much motivated me to pursue my Ph.D. degree abroad. Thanks to the marvellous country of Canada, the beautiful province of British Columbia, and the wonderful city of Vancouver and its people. I will certainly miss the quality of life I enjoyed during the years I lived here. The cultural diversity of this country and city were very important for my informal education. I consider Canada now to be my second nation, and I will treasure in my heart the memories of these years for the rest of my life. Thanks to my beloved country Brazil. And thanks to its people, very simple, very modest, but with a big heart. I believe in the potential of my country. We have a surging economy, an abundance of natural resources, and hard-working people. I believe I will still see this country as competitive as the most developed in the world. And I will be doing my part to help this goal be achieved. Finally, I wish to thank my brother Jose Eduardo Moreira, and my parents Jose Roberto Moreira and Wany Augusta Moreira for all their caring and help. M y brother and father also contributed with technical discussions, and my mother was present at moments of complicated personal problems. Thanks to a special friend Maria L. Meneses, and to my godparents, aunt and cousins. Without their unconditional love and support my dream would have never come true. xi To my parents Jose Roberto Moreira and Wany Augusta Moreira \"It's not your blue blood, your pedigree or your college degree. It's what you do with your life that counts.\" Millard Fuller (1935-), founder and president of Habitat for Humanity International (HFHI). xiii CHAPTER 1. INTRODUCTION AND OVERVIEW 1.1. BACKGROUND The Electromagnetic Transients Program (EMTP) [1,2] has been the standard tool for the digital simulation of power system transients for just over thirty years now. Prior to that, the Transient Network Analyzer (TNA) was largely used for transients simulation. The T N A extended to transient conditions the idea of the steady-state analyzer, which, for many years, was used to investigate load flow, short-circuit, and perform stability studies on a.c. systems [3]. Since the investigation of electromagnetic transients requires the consideration of very short time scales, the T N A had to be capable of assuring that all components possessed the high frequency transient characteristics of the power system equipments they represented. With the introduction of the EMTP in the late 1960's by H.W. Dommel [1], the T N A started to lose ground as the preferred tool for power system transient analyses, and by now its utilization is restricted to particular studies. The EMTP is a fully digital network simulator and has several advantages over the TNA. Firstly, EMTP-type of programs (DCG-EPRI, ATP, MICROTRAN, EMTDC, etc.) can be run on costless off-the-shelf computer workstations, whereas the TNA, as an analog simulator, is much more expensive and requires significantly more space. Another point of limitation for the TNA is that it is adequate for the analysis of transient phenomena only up to about 10kHz, since in the presence of higher frequencies, difficult aspects of modelling like parasitic capacitances have a significant importance. A major advantage of the EMTP is that very accurate and sophisticated models for the system components can be developed and implemented. For transmission lines, for example, it is not possible with scaled down analog models to simulate the distributed nature of the line parameters [4]. In fact, the EMTP attracted so much attention and was so widely used by power engineers in 1 the 1980's that individuals and groups started to contribute with newer and more elaborate models for transmission lines and cables [5,6,7,8], transformers [9,10,11], synchronous machines [12,13], and non-linear elements [14,15]. The development and implementation of new models in the several EMTP versions continue through the present days. Although well established as the major tool for electromagnetic transient studies in power systems, the EMTP has an important disadvantage when compared to the TNA: it is incapable of performing real-time simulation (at least in its traditional versions). When a very sensitive equipment needs to be tested under similar conditions to the ones it will face in the actual power system, real-time simulation is particularly indicated. Therefore, the equipment can be tested in a closed loop scenario with the digital network simulator representing the real system and the equipment interfaced with the simulator through digital/analog converters. Examples of such equipments are protective relays, and High Voltage Direct Current (HVDC) and Flexible Alternating Current Transmission System (FACTS) control systems. The real-time simulation requirement provided the motivation for the development of real-time digital simulators that combine the real-time simulation capability of the T N A with the flexibility and sophistication of the EMTP. Over almost a decade now, thanks to the continuous advance in hardware capability, there are a number of real-time digital simulators available [4,16,17,18]. The real-time research group at the Department of Electrical and Computer Engineering of the University of British Columbia has been continuously advancing the concept of a full-size real-time power system simulator, which has eventually resulted in the Object Virtual Network Integrator (OVNI) power system simulator [19]. OVNI is aimed at representing full-size power system networks for real-time and on-line simulations using off-the-shelf Pentium-based personal workstations. This requires a very efficient network solution, which takes advantage of network partitioning techniques. 2 The motivation for the development of this thesis work was provided by the possible exploitation of a particular partitioning technique, based on the widely separated time constants that electric networks may present. 1.2. THESIS M O T I V A T I O N One of the natural ways to separate solution regions is to recognize that some subsystems or components have long time constants, and relatively large integration steps At are sufficient for an accurate solution, while some other subsystems or components have short time constants and require much smaller At's. In these cases it is wasteful to solve the slow system at each integration step of the fast system solution. However, EMTP-type of programs does not offer the possibility of using different time steps according to the necessities of the subsystems. The same discretization step size At is used for all parts of the network and therefore, due to the overall required accuracy, the common step size is determined by the solution needs of the \"fastest\" subsystem in the network. There have been some successful reports on the exploitation of different discretization steps taking advantage of the separate dynamics in a network, and, even further, on the combination of different solution methods, such as those used for electromagnetic transients and transient stability. However, no attempt has ever been made to generalize the concept of network simulation with dual or multiple time steps for implementation exclusively on EMTP-type of programs. The capability of simulating a network with different time steps in different subparts of it is termed \"latency exploitation\". Latency is viewed as a network property where some signals remain constant or change very slowly in comparison to the rest of the network. The part of the network with such behaviour is said to be latent, as in \"dormant\". 3 Latency exploitation may be particularly important in the field of power system transients for the simulation of power electronic converters and controllers. With relatively new switching devices such as the gate turn-off thyristors (GTO's) and insulated gate bipolar transistors (IGBT's) very high switching frequencies are obtained, and small integration steps are needed for an accurate simulation. In contrast, the power network may be subjected to much lower frequencies and a larger integration step may provide the accuracy desired. Also, an important new line model, the z-line [8], in principle similar to the cp-line model, can ideally be suited for latency exploitation since it requires the partitioning of transmission lines into several segments for an accurate simulation. A very small time step, equivalent to the travelling time of the short line segments should be used for the solution update of the internal nodes of the line, while the subnetworks connected to both the sending and receiving ends of the line could be simulated with a much larger time step. In these situations, gains in the order of ten to a hundred times can be achieved applying the techniques developed in this thesis work. In the next chapter, a review of the available literature on latency exploitation will be performed. The main purpose of this thesis is to present in detail the methodology that allows a successful network simulation using distinct time steps in different parts of it, and to apply the latency technique in case studies of networks containing lumped and distributed-parameter elements. 1.3. THESIS ORGANIZATION Chapter 2 reviews the technical literature on the topic of numerical circuit simulation with multiple time steps dating back to the simulation of Very Large Scale Integration (VLSI) components in the 1980's. In power networks, examples presented include simulations that 4 couple solution methods for electromagnetic transients and electromechanical transients, and approaches in which the integration step is adjusted during the simulation to reflect the change in the dynamics of the network with time. Also, the partitioning method [19] that provides the solution framework in which this thesis is built upon will be introduced. Chapter 3 presents a thorough review on the time-domain integration rules that are used in EMTP-type of programs. Topics especially discussed are the stability and accuracy of these rules. In chapter 4, the latency technique and its implementation in a time-domain EMTP-type of solution are presented in detail. The expanded integration rules obtained in the latency approach are analyzed in terms of their stability and accuracy. Case studies are presented in chapter 5. The latency technique is applied to networks containing exclusively lumped elements and networks containing transmission lines, modelled with the travelling wave equations, as well. Cases with transmission lines modelled as ideal lossless lines and with the constant-parameter line (cp-line) model are presented. The cp-line model is tested under the latency environment for a future implementation of the z-line model exploiting latency. The z-line model [8] is a full frequency-dependent phase-domain transmission line model that eliminates the necessity of synthesizing the frequency dependent transformation matrix that decouples modal propagation. The z-line model will be discussed in more detail when appropriate. Finally, chapter 6 presents the conclusions of the thesis, emphasizing the original contributions of this work. This is a new and promising field of investigation and many suggestions for future research are discussed as well. 5 CHAPTER 2. LITERATURE REVIEW 2.1. N E T W O R K P A R T I T I O N I N G As previously mentioned, this thesis work focuses on \"latency exploitation\", i.e., the capability of simulating a network with distinct time steps for its different unknown variables. Although it is theoretically possible to numerically solve an electric network in time-domain applying different integration steps for individual component elements, it would be desirable to apply latency in the context of a block partitioned matrix, where the complete network is divided into subnetworks coupled through links. Each subnetwork could then be simulated with its own time step according to its required simulation accuracy. Information would then be exchanged from one subnetwork to the others through the links. This approach guarantees computational efficiency and a straight forward implementation procedure. The network partitioning scheme followed in this thesis work [19] is based on the concept of Diakoptics, introduced by G. Kron in the 1950's [20] for the solution of large power system networks. Despite its advantages, Diakoptics did not fulfil its potential of becoming a universal tool for the study of large and complex power systems, due especially to the success of Sparsity Techniques. While Sparsity Techniques take advantage of the very sparse nature of the admittance matrix, Diakoptics splits the network into dense subsystems connected by a few links. In the work reported in [21,22,19], Diakoptics is extended to the Multi-Area Thevenin Equivalent (MATE) concept. M A T E formulates the main ideas of Diakoptics in a conceptually clear and computationally efficient solution framework. Its characteristics will be explained in more detail in the upcoming chapters and in appendix A. Once established the solution framework in which latency exploitation will be based, the attention now can be turned specifically to the historical developments in latency. 6 2.2. THE ORIGIN OF LATENCY EXPLOITATION The limitations of a single and fixed-size step simulator for electric networks were first reported by researchers in the area of VLSI simulation, as well as the initial attempts to overcome these limitations [23,24,25,26]. The SPICE circuit simulator has played a key role in electronics circuit simulation since the 1970's. However, when VLSI circuits were simulated, the efficiency of SPICE reached a bottleneck. Although very accurate, SPICE and its optimized versions could not simulate VLSI circuits in detail in a time and cost effective manner [23]. This prompted various industrial and academic research groups to pursue advanced techniques for faster circuit simulation without sacrificing its accuracy. It was on this context that researchers turned their attention to latency exploitation. Before continuing on to the contributions previously achieved in latency it is important to clarify and differentiate the nomenclature used in these works. First of all, it is important to differentiate between time-domain latency and iteration latency [23]. During circuit simulation, latent sections of the circuit have smaller changes in their state variables than other more active sections. They may not even change at all during a certain time range. Latent sections can therefore be simulated with larger time steps than active sections without incurring in any significant loss of accuracy. Approaches which use larger time steps for latent sections are exploiting time-domain latency. On the other hand, i f a uniform time step is used, latency can still be exploited to solve the nonlinear equations of the circuit. Latent sections require fewer iterations for the solution of the nonlinear equations. Approaches using the increased convergence rate of latent sections are exploiting iteration latency. The latter form of latency, however, provides very limited gains since usually a latent section can save only one or two iterations [23]. In this thesis work, latency is applied only in its time-domain form. Furthermore, since this is a latency implementation for an EMTP type of simulation, nonlinear 7 elements may usually be approximated by its well defined piecewise linear representations [2]. Iteration latency exploitation for the solution of nonlinear equations is out of the scope of this thesis. Some authors prefer to differentiate the term latency from the term multirate behaviour. In [25], the authors define latency as long periods over which the value of a given signal remains constant. Latency here may apply to a single variable and can be exploited simply by varying the time step through some mechanism that detects the fact that the given signal value is not changing appreciably, and also detects when the latency approach is no longer valid. Latency is viewed as a subset of a more general property of waveforms called multirate behaviour. It is only at the multirate behaviour level that the authors refer to signal values changing at different rates relative to one another, over the same interval of time. In this thesis work, it is this property that is exploited, although the name latency exploitation will be maintained. Here, the terms latency and multirate behaviour are indistinguishable. As a matter of fact, the variable step size approach for transient simulation will not be explored in this work, since it may be accurate only for some specific conditions, as it will be later explained. Some contributions in this area, however, will be cited when opportune. 2.3. THE WAVEFORM RELAXATION METHOD The initial attempts to exploit latency in VLSI circuit simulation were based on the waveform relaxation method [27]. Some successful contributions using this and similar methods, such as the nonlinear relaxation, are reported in [25,26]. The waveform relaxation method is an iterative method that can exploit time-domain latency to attain significant speed gains with SPICE like accuracy. It is based on the Gauss-Seidel and Gauss-Jacobi relaxation methods used for solving large systems of algebraic equations, expanding them for the numerical integration of 8 a system of differential-algebraic equations (DAE). The D A E system is broken into subsystems (like in MATE) , which are solved independently. However, unlike M A T E , the waveform relaxation method does not provide a tie between the subsystems (the links in MATE) and therefore, delays are introduced in the solution of the system. Each subsystem uses the previous iterate waveforms as \"guesses\" about the behaviour of the state and algebraic variables in other subsystems [28]. The waveforms are then exchanged between subsystems, which are resolved with improved information about the other subsystems. In real-time simulation, however, due to the requirements of sustained continuous simulation, all parts of the network have to be solved together, i.e., iterations through At delays among subsystems are not allowed. According to [23], the waveform relaxation method has proven to be successful only for special classes of circuits and perform poorly for many other types of circuits. Iterative matrices methods present slow rates of convergence and lack robustness when compared to direct methods that use some form of implicit integration method. The waveform relaxation algorithm was explored for the first time in the context of power systems in [29], and it was further investigated for transient stability simulation on a parallel processing scenario in [28]. The authors claimed that the characteristics that prevented the waveform relaxation method from achieving its full potential for VLSI circuit simulation, such as the presence of floating capacitors, are not encountered in power system problems. One of the attractions of the waveform relaxation method is that it is well suited for parallel simulation [26,28], since each processor can be assigned the task of simulating a particular subsystem. Vast improvements in simulation time can be achieved i f the system partitioning is optimized, grouping together the latent subsystems under the responsibility of a small number of processors. Direct methods exploiting latency, however, can also be well suited for parallel processing systems under the M A T E context [21]. The fundamental ideas regarding the implementation of latency simulation using direct methods in a parallel processing environment were outlined in [30] and will be revisited on the recommendation for future work 9 section. The latency algorithm developed in this thesis falls under the domain of the direct method solution techniques, as implemented in the EMTP. The choice of this approach is based on the accuracy and robustness characteristics offered for almost all circuits. 2.4. THE VARIABLE TIME STEP METHOD As mentioned before, time-domain latency may be applied in the context of variable time steps. The idea is to dynamically adjust the step size during the simulation. If a component is identified as being latent, then its integration step may be increased. Problems arise i f the latency condition is not detected properly. For example, i f a component is identified prematurely as being latent, then any small error in its value will be propagated to the other components producing errors in their solutions [25]. If a component is thought to be latent but is actually changing very slowly, the resulting solution may be deceiving. Applying a variable, but single time step for power systems simulation may be efficient in certain cases, as described in [31,32]. However, it lacks the generality required from a real-time simulator which is used in a closed loop circuit scenario, usually targeting equipment testing. The variable step method may be well suited for simulating dynamic systems which are primarily slow response systems, but exhibit fast decaying transients. In this case, the simulation may be initialized with a small time step to track the fast transients. If these fast transients decay after a short period of time, the step size may be increased in order to track only the remaining slow transients. Computational time is thus saved in this case by adjusting the time step size. There is, however, a problem i f the fast dynamics are again excited during the simulation, although this limitation may be overcome by developing a \"wake-up\" condition [25] (opposite to the \"dormant\" condition), and dynamically adjusting once again the step size. 10 A more serious limitation of the variable step method occurs i f the system presents fast transients sustained for a large portion of the simulation interval. As reported in [33], this is the case, for example, when the power system modelled includes fast switching devices, such as FACTS controllers or converters necessary to intertie DC lines into an A C system. Under these conditions, although only a few of the state variables are affected by the fast dynamics, the time step must remain small, thus defeating the goal of increased computational efficiency. There is also an inherent condition that prevents the variable step size method from being of significant advantage for real-time simulators. This has to do with the synchronized solution output from the simulator. By definition, a real-time simulator has to be capable of performing the calculations for the simulation of a time step At faster than the elapsed real time At itself (it cannot be equal due to the time overhead incurred with the use of D/A and A/D converters and other hardware overheads). However, the solution time per time step of the simulator depends on the efficiency of the code and the hardware capability. Therefore, the solution time cannot be adjusted to meet the demands of a variable integration step size. In general, it is very important to keep the simulation time step constant in a real-time simulator, so that real-time calculation is always guaranteed. Although the proposed algorithm for latency exploitation developed in this thesis work has not been implemented under a real-time simulation scenario, the purpose is to lay the foundation for this next step. An effort has been made to exploit the latency characteristics of electric power networks in a way that facilitates its adaptation to a real-time simulator without compromising accuracy. Therefore, it has been decided that the best approach to exploit latency is on its time-domain form and using dual or multiple time steps with implicit integration methods. 11 2.5. EXPLOITING LATENCY USING IMPLICIT INTEGRATION RULES The implicit integration method normally used in EMTP-type of programs is the trapezoidal rule, which is numerically A-stable, i.e., the region of absolute stability includes the entire left-half plane of the complex plane C. Therefore, i f the trapezoidal rule is applied to a stable linear system of equations, the step size At can only influence the accuracy of the simulation, and not its numerical stability. Figure 2.1 shows the trapezoidal rule applied for the numerical integration of a function of the type / (x, t) - • f(x,t) Figure 2.1: Trapezoidal integration rule The value of x at t = tn+i is given by = Xn + \"^r [/\"fa , 0+ f(Xn+\\ > )] (2.1) The other implicit integration method used in certain conditions in EMTP-type of programs is the backward Euler rule. This rule is also A-stable and, although it may not be generally as 12 accurate as the trapezoidal rule, it presents some characteristics related to the damping of numerical oscillations that make it attractive for use in time-domain simulations [34]. These characteristics will be further explored in this thesis. Figure 2.2 shows the backward Euler rule applied for the numerical integration of the same function as in figure 2.1 f(x,t) f(Vo) At Figure 2.2: Backward Euler integration rule The value of x at t = tn+i is given by =*«+A'[/(X+i>'»+i)] (2.2) Latency exploitation will be tested with both the trapezoidal and the backward Euler rules of integration. Implicit integration methods are more suitable than explicit techniques (e.g., Euler, Runge-Kutta) for handling numerical simulations of systems with increased stiffness, which is associated with the range of time constants in the system models. Stiffness is measured by the ratio of the largest to smallest eigenvalues of the system. Explicit integration methods have weak numerical stability and usually require small step sizes in order to prevent instability. EMTP-type 13 of programs are well known for their capacity of modelling power systems in great detail. As the modelling detail of a system increases, so does its stiffness, and since stiff systems are potential candidates for latency exploitation, the application of implicit integration methods is justified. Latency exploitation with implicit integration methods was applied to VLSI circuit simulation as early as the mid 1980's [35]. However, as explained in [23], the contribution from the latent subcircuits when solving only for the active ones were poorly computed. The contribution was either neglected completely (as i f the circuits were totally decoupled) or computed with polynomial extrapolation of the last obtained values. Both approaches are severely limited since most latent subcircuits are not completely dormant and polynomial extrapolations are inherently inaccurate for predicting the exponential curves typical in electric networks. In [23], an approach called Overdetermined Polynomial Integration (ODPM) is used for evaluating the contribution of the latent subcircuit to the active ones. However, this method is used for the integration of nonlinear equations, which require the calculation of the Jacobian matrix. In this thesis work, only linear systems are considered and the contribution of the latent subcircuits is evaluated in a much simpler form. In the context of power systems, latency exploitation employing implicit integration methods have been reported in [33,36] for transient stability studies and in [37,38] for electromagnetic transient studies. In [33], the trapezoidal integration rule was used for numerical integration of the fast time scales. However, in order to account for the contribution of the latent subcircuits, explicit integration methods were still used for the slow time scales. Linear interpolation was then performed between calculated values. In a subsequent paper [36], the authors claim to have improved their approach by using implicit integration methods for both fast and slow time scales. Although this reference does present important contributions to the subject of multirate simulation, no explanation is provided on how the contribution of the latent subcircuits is taken into account when using only implicit integration methods. 14 As far as electromagnetic transients simulation is concerned, the possibility of latency exploitation was studied in [37]. This reference provided a good background for the start of this thesis work. The authors, however, never pursued further developments in this area or the implementation of the latency technique in an EMTP-type simulator. This thesis aims at improving some of the methods adopted in [37] in order to overcome the limitations presented in this reference. Also, in [37], only the possibility of exploiting latency was validated, while important issues such as stability and accuracy of the simulation were not deeply investigated. In [38] some aspects of numerical stability and accuracy of the implemented latency technique were presented. However, the contribution of the slow time scales was neglected when solving for the fast time scales. 2.6. COUPLING ELECTROMAGNETIC AND TRANSIENT STABILITY SOLUTION METHODS As mentioned in chapter 1, there have been some successful attempts to couple solution methods for electromagnetic and eletromechanical transients. A comprehensive formulation as well as simulation results of this approach were presented in the late 1980's [39,40]. The authors proposed this approach to overcome the limitations in ac/dc system transient studies, where the modelling of one system, a detailed dc system simulation or an ac transient stability program (TSP), is compromised by the restricted representation of the other system. In [39,40] the solution of the network is done in a non-simultaneous way, with the EMTP solution occurring prior to the TSP solution for a given large TSP time step. In [41], the method is further improved to allow the suspension of the EMTP simulation after the quasi steady-state dc system conditions have been restored. The slower system dynamics can then be accommodated by the TSP quasi steady-state simulation alone. Although the interface of the fast and slow subsystems, 15 represented by the parts of the network solved respectively by the EMTP and TSP, is relatively similar to the approach followed in this thesis, the method proposed is not suitable for real-time simulation. 2.7. LATENCY EXPLOITATION FOR POWER ELECTRONIC SYSTEMS Finally, it is worth mentioning the attempts to exploit multirate behaviour in power electronic systems. In [42], a simple switched mode power supply converter, consisting of three stages, namely, a diode rectifier, a switched-mode dc-dc converter, and a dc load, was considered for latency exploitation. The diode rectifier operates at 50 or 60 Hz and the dc load has no switching. However, since the dc-dc converter may switch from 20 kHz up to 1 MHz, the simulation time step in a traditional time-domain simulation (the smallest At in the system) is significantly restricted by this part of the circuit. A decoupling technique based on the transmission line modelling (TLM) method was used. The T L M technique consists of basic transmission line models that can be used to represent an inductor, a capacitor, or a combination of both. However, the T L M models for lumped inductances and capacitances are approximations, and the EMTP approach of more exact representation of these components is preferred in this thesis work. In [43], a decoupling technique is introduced for the detailed simulation of a power electronics circuit enclosing several widely different time constants, which vary from the switching period in the order of microseconds up to the self-heating of the components in the order of seconds and minutes. This technique, however, is based on the extrapolation of sampled nodes in the fast subcircuit and therefore is inappropriate for the stability robustness sought in this work. 16 In the following chapter, a brief summary of the numerical integration rules commonly used in time-domain simulation of electromagnetic transients will be presented. 17 CHAPTER 3. TIME-DOMAIN SIMULATION OF ELECTROMAGNETIC TRANSIENTS: A BRIEF OVERVIEW As previously mentioned, the trapezoidal integration rule is the most common implicit integration method used in EMTP-type of programs. However, in certain conditions where the damping of numerical oscillations introduced by the trapezoidal integration rule is required, the backward Euler integration rule is also employed for time-domain simulations, at least for a certain number of time intervals. Before presenting in this thesis an expanded integration method and its characteristics for latency exploitation, it is important to briefly review the trapezoidal and the backward Euler rules, as well as, their accuracy and stability. This chapter is dedicated to this task, while chapter 4 will present the methodology proposed for efficient latency exploitation. 3.1. N E T W O R K D I S C R E T I Z A T I O N I N T I M E D O M A I N Since closed-form analytical solutions of large systems in time-domain, as in the case of power systems, are virtually impossible to derive, the best strategy to obtain the system solution is to discretize time and obtain a successive step-by-step solution of the system. Whenever this approach is followed, issues regarding the time step size and the choice of the discretization rule need to be evaluated in order to guarantee an accurate and stable solution. From the analyses of simple inductances and capacitances, it is easy to represent these elements in a discrete-time form for both the trapezoidal and backward Euler rules. 18 3.1.1. THE TRAPEZOIDAL INTEGRATION RULE Equation (2.1) is repeated below for reference: *„+i = x „ +^[/(*fl>0 + / ( W n + i ) ] (2-1) A. INDUCTANCE In continuous-time, the relationship between voltage and current in an inductance is given by the expression: . . di, (t) Integrating both sides with respect to time, from t-At to t, and after some manipulations, eq. (3.1) can be written as h (0=iL ( ' - A 0 + f ^ h ( 0 + v L ('-A0] (3-2) Equation (3.2) is now in the same format as eq. (2.1). For network analysis it is convenient to re-write eq. (3.2) as where hL (t) = —vL (t-At)+iL (t-At) is the history term, which depends on the voltage and current values at the previous time-step only. 19 Figure 3.1 shows the relationship between current and voltage for the continuous-time inductance and the discrete-time inductance. i(t) m , i ( t ) k • v(t) continuous-time L G = At / 2L -A/W m hL(t) v(t) discrete-time L Figure 3.1: Inductance represented in continuous-time and its equivalent representation in discrete-time with the trapezoidal rule In discrete-time, the inductance L can be represented by a constant conductance G = in parallel with a current source liL(t). B . C A P A C I T A N C E In continuous-time, the relationship between voltage and current in a capacitor is given by the expression: ic(t) = C dvc (t) dt (3.4) Integrating both sides with respect to time, from t-At to t, and after some manipulations, eq. (3.4) can be written as At vc (0 = v c (t - At)+—[ic (t)+ic (t - At)] (3.5) 20 Equation (3.5) is also in the same format as eq. (2.1). Equation (3.5) can be rewritten as 2C J'c(0 = -7J-Vc(0+Ac(0 (3.6) 2C where hc (f) = —-^-v c (r —Ar)-? c (t-At) is the history term, which depends on the voltage and current values at the previous time step only. Figure 3.2 shows the relationship between current and voltage for the continuous-time capacitance and the discrete-time capacitance. i(t) k • o m , i ( t ) k • o v(t) continuous-time C G = 2C/ At - A A A / -m hc(t) v(t) discrete-time C Figure 3.2: Capacitance represented in continuous-time and its equivalent representation in discrete-time with the trapezoidal rule In discrete-time, the capacitance C can also be represented by a constant conductance 2C G in parallel with a current source hc(t). At 3.1.2. THE BACKWARD EULER INTEGRATION RULE Equation (2.2) is repeated below for reference: (2.2) 21 A . I N D U C T A N C E Integrating both sides of eq. (3.1) with respect to time, from t-At to t, using the backward Euler rule of integration; and after some manipulations, eq. (3.7) is obtained. iL(t) = iL(t-At)+^-vL(t) (3.7) Equation (3.7) is in the same format as eq. (2.2). It is convenient to re-write eq. (3.7) as (3.8) where hL(t) — iL(t-At) is the history term, which depends on the current value at the previous time step only. Figure 3.3 shows the relationship between current and voltage for the continuous-time inductance and the discrete-time inductance. i(t) m , i ( t ) k • o v(t) continuous-time L G = A t / L -AMr m hL(t) V(t) discrete-time L Figure 3.3: Inductance represented in continuous-time and its equivalent representation in discrete-time with the backward Euler rule In discrete-time, the inductance L can be represented by a constant conductance G = — in parallel with a current source hL(t). 22 B. C A P A C I T A N C E Integrating both sides of eq. (3.4) with respect to time, from t-At to t, using the backward Euler rule of integration; and after some manipulations, eq. (3.9) is obtained. vc=vc(t-At)+^ic(t) (3.9) Equation (3.9) is also in the same format as eq. (2.2). It is convenient to re-write eq. (3.9) as 'c(0 = £ v c ( 0 + A c ( 0 (3-10) c where A c (/)==——vc (f — Af) is the history term, which depends on the voltage value at the previous time step only. Figure 3.4 shows the relationship between current and voltage for the continuous-time capacitance and the discrete-time capacitance. i(t) k • o m , i(t) k • o v(t) continuous-time C G = C/ At A W m hc(t) v(t) discrete-time C Figure 3.4: Capacitance represented in continuous-time and its equivalent representation in discrete-time with the backward Euler rule In discrete-time the capacitance C can also be represented by a constant conductance C G = -— in parallel with a current source hc(t). AJ 23 / 3.2. ACCURACY OF INTEGRATION RULES Solving first order differential equations numerically using integration rules such as the trapezoidal and backward Euler rules is never completely exact. Since, at each time-step, the area under the curve is approximated by a trapezium (trapezoidal rule) or by a rectangle (backward Euler rule), a truncation error is always incurred. Ideally, i f the time step could be reduced to a value in the order of the numerical precision of the computer, the truncation error would be insignificant. However, execution time would then increase beyond practical levels. The solution is to find out a compromise between the level of accuracy required and a practical simulation time. Fortunately, the truncation error decreases very fast with the reduction in time step (the truncation error is in the order of At 3), and enough accuracy can most often be guaranteed. Once an analog signal is sampled in time, as is the case in digital computers, integration rules act as digital filters, eliminating frequency components of the signal that are beyond the Nyquist frequency, which is given by 4 = i (3.1D where At is the time-step adopted for network discretization. Therefore, it is important and more accurate to assess the performance of integration rules in the frequency domain as it will be performed in the next subsections. 3.2.1. TRAPEZOIDAL INTEGRATION RULE In continuous-time, eq. (3.1) gives the voltage-current relationship in an inductance, and its respective admittance in the frequency-domain is given by 24 In discrete-time, using the trapezoidal rule, the voltage-current relationship is given by eq. (3.2), which can be re-written as h{t)-h(t~&) = fjrVL{t)+^vL{t-At) (3.13) Applying an input of the form vL(t)=eJ0\" and assuming an output of the form iL {t) = Ye{co)ejax, where Ye(co) is the admittance in the frequency-domain, lead to Ye(cv) = At eJaAt+l 2L ejoAt-\\ (3.14) Accordingly, the impedance Ze ( -0.2 -0.4 -0.6 -0.8 -1 -i 1 1 1— _i_ 0 0.2 0.4 0.6 0.8 1.0 frequency (per unit of Myquist frequency) Figure 3.5: Frequency response in amplitude and phase of the trapezoidal integration rule 26 It is clearly verified that for frequencies well below the Nyquist frequency (1/10 to 1/5 of it), the distortion factor is very close to unity, and therefore results are very accurate in this range. On the other hand, the distortion factor grows very quickly as frequency approaches the Nyquist frequency, and the simulation results become totally inaccurate. One important aspect to mention is that the trapezoidal rule does not present any phase distortion in the frequency spectrum up to the Nyquist frequency. 3.2.2. B A C K W A R D E U L E R I N T E G R A T I O N R U L E In discrete-time, using the backward Euler rule, the voltage-current relationship in an inductance is given by eq. (3.5), which can be rewritten as iL(t)-iL(t-At) = ^ -vL(t) (3.20) Performing the same steps as in the trapezoidal case: Ye(eo) = At e jaAt (3.21) Equation (3.21) can also be written as Ye(co) = '&} }Lj + K2LJ j tan f coAt^ v 2 , (3.22) From the form of Ye(a>): (3.23) 27 where Re and Le(aj) are given by Re = — (constant resistance) tan Le(co) = L-coAt (same as trapezoidal) (3.24) (3.25) with the distortion factor tan Ke{aS) = —r f coAt^ coAt (same as trapezoidal) (3.26) The distortion factor of the equivalent inductance for the backward Euler rule is the same as in trapezoidal. However, the backward Euler rule adds a resistance to the equivalent circuit in frequency-domain of a discrete-time inductance. In the case of a capacitance this resistance would be Re = At 2C (3.27) The presence of this resistance causes additional losses to the discretized elements, in comparison to their continuous-time representation. The effect on network simulation using the backward Euler integration rule is an increase in the damping characteristics of the network. This is related to the use of the backward Euler rule for the damping of numerical oscillations introduced by the trapezoidal rule [34]. Once again, it is desirable to obtain the frequency response of the equivalent admittance Ye (co) in comparison to the exact continuous-time admittance Y (ca). 28 1 1 f requency (per unit of Nyquist f requency) f requency (per unit of Nyquist f requency) Figure 3.6: Frequency response in amplitude and phase of the backward Euler integration rule The backward Euler rule presents a linearly varying phase distortion as observed in figure 3.6, thus making it usually inappropriate for sustained long simulations. However, combining the better accuracy of the trapezoidal rule with the stronger damping capabilities of the backward Euler rule results in a very efficient approach. This method is called Critical Damping Adjustment (CDA), and is thoroughly described in [34]. 29 3.3. STABILITY OF INTEGRATION RULES In this section, the stability of the trapezoidal and backward Euler integration rules will be analyzed. The z-transform is applied to the difference eqs. (3.13) and (3.20) and the poles and zeros of the resulting transfer functions are calculated. 3.3.1. TRAPEZOIDAL INTEGRATION RULE Applying the z-transform to eq. (3.13), the following transfer function is obtained: Ye{z) = V{z) ( At V + l z -1 (3.29) This is the same as eq. (3.14) substituting eJ ' by z. Equation (3.29), which gives the transfer function of an integrator, has a pole at z = 1 and a zero at z = - 1 , and the system is stable as both an integrator and differentiator. However, the pole at z = -1 for the equivalent differentiator introduces numerical bounded oscillations at discontinuities. Although the numerical oscillation is bounded for each discontinuity, simulations for fast and sustained switching circuits, such as H V D C converters and pulse width modulation (PWM) inverters, would quickly become unstable, since the amplitude of the numerical oscillation would double at each subsequent discontinuity. 30 3.3.2. B A C K W A R D E U L E R I N T E G R A T I O N R U L E Applying the z-transform to eq. (3.20), the following transfer function is obtained: V[z) [L )z-l This is the same as eq. (3.21) substituting ejaAt by z. Equation (3.30) has a pole at z = 1 and a zero at z — 0, and the system is again stable as both an integrator and differentiator. The pole at z = 0 for the equivalent differentiator causes the system to be critically damped, eliminating the numerical oscillations generated by the trapezoidal rule. This is the property exploited in C D A [34]. In the next chapter, the latency simulation technique proposed in this thesis work will be fully explained, and the resulting expanded trapezoidal and backward Euler rules will be checked for their accuracy and stability properties. 31 CHAPTER 4. PROPOSED METHODOLOGY FOR LATENCY EXPLOITATION AND ITS MAIN CHARACTERISTICS In this chapter, the numerical characteristics of the method proposed for accurate and efficient latency exploitation will be described. The modifications performed on the trapezoidal and backward Euler integration rules for a natural adaptation to an EMTP-type simulator are fully explained. Also the stability and accuracy of the adapted integration rules are verified and an eigenvalue analysis is proposed for the verification of electric networks that are suitable for latency exploitation 4.1 T H E L A T E N C Y S I M U L A T I O N M E T H O D As previously mentioned, latency is exploited in this thesis work in the context of a block partitioned matrix where the complete network is divided into subnetworks coupled through links. The partitioning technique is based on the recognition that some subsystems or components have big time constants, and relatively large time steps are sufficient for an accurate solution, while some other subsystems have small time constants and require much smaller time steps. This partitioning technique fits nicely into a latency exploitation environment. However, the interface between the fast and slow subsystems has to be cautiously adjusted in order to guarantee an accurate and efficient simulation. Incorrect approaches in the coupling of the slow and fast subsystem solutions can lead to strong numerical oscillations and divergence of long-term simulation. Assuming that a generic electric network has been partitioned into a fast and a slow subsystem, figure 4.1 shows a schematic diagram for visualization of the network decoupling. 32 fast subsystem slow subsystem Figure 4.1: Schematic diagram of the network decoupling Each subsystem contains any combination of voltage and current sources, active and passive elements, and transmission lines. They are linked through a resistance in this figure, although any other element could have been used as a link. The M A T E concept [19] guarantees an exact solution of the complete network i f both subsystems are first solved separately as i f they were completely decoupled, and the node voltages are subsequently updated taking into account the current flowing through the link. However, the results match exactly the ones obtained from traditional EMTP-type of programs only i f the same integration step is used in the transient simulation of both subsystems. Another important property of the M A T E technique extensively used in the latency approach is the network reduction into a Thevenin equivalent, easily obtained through matrix manipulation. In appendix A , the M A T E formulation is clearly stated. Even i f different integration steps are used, it is still possible to guarantee very accurate results using latency i f each of the subsystems can be characterized as being either slow or fast, i.e., i f the dominant eigenvalues in each of the subnetworks are fairly independent. A complete generalization of the latency methodology for application in a generic network is out of the scope of this thesis, although some important contributions are outlined in section 4.3. 33 In the technique proposed in this, thesis work, the large time step should always be a multiple of the small time step. Defining AT as the large time step (for the slow-varying subsystem) and At as the small time step (for the fast-varying subsystem), this requirement can be expressed as AT = nAt (4.1) where n is a positive integer larger than unity Figure 4.2 shows the timeline of the simulation at a general time interval when latency is exploited. k ^ T (k+1)AT (k+2)AT kAT + At kAT-•-2 At kAT+ {n-1)At (k+1)A T4At (k+1) AT+2 At (k+1)z •\\T+(n-1)At n is the sarre as in equation (4.1) k is a non-negative integer Figure 4.2: Timeline of the simulation with the large integration step AT containing several small integration steps At In order to guarantee a stable solution, both the slow-varying subsystem and the fast-varying subsystem are solved together at instants multiple of AT. The novelty of the method is that for the small time steps At's located between two consecutive AT's interval, only the fast-varying subsystem is solved: Two important issues arise here that are fundamental to the success of the latency approach proposed: 1) How to consider the contribution of the latent subcircuits when solving only for the active ones at small time step At intervals; 2) How to consider the contribution of the fast-varying subsystem through a complete AT interval when it is time to solve for the complete network. 34 The approach followed in this thesis work tries to find a compromise between maximum accuracy and efficiency of the simulation. 4.1.1. HISTORY SOURCE ACCUMULATION FOR RESYNCHRONIZED SOLUTION Whenever it is time to solve for the complete network (at every instant multiple of the large time step AT), the history sources for both the slow-varying and the fast-varying subnetworks should be updated. For the slow-varying elements, the previous solution was calculated at t-AT, and therefore the history sources are naturally updated from this solution. However, during the large time interval AT, the fast-varying elements have been solved a number of times, the exact quantity depending on the ratio AT/At. It is, therefore, more accurate to consider all the solutions obtained for the fast subnetwork within the interval AT when evaluating the new history sources of the fast components. If only the last solution is used in the update of the history sources, they may turn out to be of an ill-defined value and thus inadequate for reliable calculations. This problem is solved by applying an integration rule (trapezoidal or backward Euler in this thesis work) to each of the fast-varying components within the large interval AT. A. TRAPEZOIDAL INTEGRATION RULE Figure 4.3 illustrates the procedure when using the trapezoidal integration rule: 35 vL(t) kAT kAT+At kAT+2At kAT4{n-1)At (k+1)AT t Figure 4.3: Expanded trapezoidal rule for history source accumulation The expanded trapezoidal rule is developed here for an inductor and a capacitor respectively. Starting from eqs. (3.1) and (3.4): ic(t) = C diL (t) dt dvc(t) dt (3.1) (3-4) Integrating these equations over the large time step AT: t J* vL(t)dt = L[iL(t)-iL(t-AT)] (4.2) 1-&T j ic(t)dt = C[vc(t)-vc(t-AT)] (4.3) t-AT 36 Applying the trapezoidal rule in order to evaluate the integrals in eqs. (4.2) and (4.3) results in At r - \" _ l J vL(t)dt = —[vL(t)+vL(t-AT)] + At^vL(t-AT + kAt) (4.4) t-AT 1 *=' J ic(t)dt = — [ic(t)+ic(t-AT)] + At^ic(t-AT + kAt) (4.5) t-AT 2 k=\\ n-1 where n = At The difference equations giving the voltage-current relationships for an inductor and a capacitor in discrete-time can then be written as ^(0 = ^ (0+^(0 (4-6) 2 C Equations (4.6) and (4.7) are identical to eqs. (3.3) and (3.6), respectively. The history terms though are modified and are given by At At hL(t) = iL(t-AT)+—vL(t-AT)+—^vL(t-AT + kAt) (4.8) 2C , • hc(t) = vc(t-AT)-ic(t-AT)-2^ic(t-AT + kAt) (4.9) Equations (4.8) and (4.9) give the history sources for an inductor and a capacitor respectively, with all the information gathered from these elements within a large integration step AT. This avoids the quasi-randomness of only taking into account the single values of the last calculated solution of the fast-varying subnetwork. Although, in principle, it may seem that a 37 significant increase in computational burden is incurred, the reality is that in the code implementation, the summation part of the history terms in eqs. (4.8) and (4.9) is updated at the end of each small time step At. As a consequence, the number of numerical operations does not increase when the history terms of the fast-varying elements have to be finally updated at the time instants when the complete network has to be solved together. B. B A C K W A R D E U L E R I N T E G R A T I O N R U L E Figure 4.4 illustrates the history source accumulation procedure for the backward Euler integration rule. t kAT kAT+At kAT+2A» kAT+(n-1)At (k+1)AT t Figure 4.4: Expanded backward Euler rule for history source accumulation As for the case of the trapezoidal rule, the expanded backward Euler rule is developed here for an inductor and a capacitor respectively. Applying the backward Euler integration rule in order to evaluate the left-hand terms of eqs. (4.2) and (4.3) results in 38 n-1 J vL(t)dt = AtvL(t)+Al^vL(t-AT + kAt) (4.10) t-AT k=\\ n-\\ j ic(t)dt = Atic(t)+At^ic (t-AT + kAt) (4.11) t-AT k=\\ The difference equations giving the voltage-current relationships for an inductor and a capacitor in discrete-time can then be written as ^(0=^(0+^(0 (4-12) 'c(0 = £ - v c ( 0 + A c ( 0 (4-13) And the history terms are given by At \"_1 hL(t) = iL(t-AT) + — ^ vL (t-AT + kAt) (4.14) L k=\\ hc(t) = -^-vc(t-AT)-^ic(t-AT + kAt) (4.15) 4.1.2. THEVENIN EQUIVALENT INTERPOLATION FOR FAST-VARYING SUBNETWORK SOLUTION As described in chapter 2, one of the most important requirements of a simulator exploiting latency is to accurately and efficiently evaluate the contribution from the latent subnetwork when solving only for the fast-varying elements. Since polynomial extrapolations are inherently inaccurate for predicting curves in electric networks, attention in this work has been turned to the possibility of using interpolation instead, in the manner described as follows. 39 Referring again to the schematic diagram in figure 4.1, it is readily observed that the link current can be exactly calculated i f both subsystems are first transformed into their respective Thevenin equivalents. Indeed, this is one of the most important properties of M A T E [19]. The Thevenin equivalent voltage source represents the combined effects of the voltage, current, and history sources existent in the original discretized network. After a particular solution instant when the complete network has been solved, there will be a number of solution steps (depending on the ratio AT/At) where only the fast subsystem needs to be solved. Thus, the representation of the slow subsystem as a Thevenin equivalent for the interval AT seems to be ideal. Figure 4.5 shows the Thevenin equivalent substituting the detailed representation of the latent subnetwork. slow subsystem f a s t subsystem Figure 4.5: Network decoupling with the detailed slow subcircuit replaced by its Thevenin equivalent If the Thevenin equivalent is maintained constant throughout the simulation interval AT, the latent subnetwork would be represented as completely dormant, resulting in a possibly inaccurate solution. If a varying Thevenin equivalent is applied at each small time step At, accuracy may be significantly increased. Extending the idea of a Thevenin equivalent representation to the active subnetwork as well, allows the complete network to be represented as shown in figure 4.6. 40 R t h slow AAAr th slow slow subsystem ARJinJ< AAAr 'link R t h fast AAAr '08 ruv, thfast fast subsystem Figure 4.6: Thevenin equivalent representation of slow and fast subnetworks at linking nodes As discussed in detail in appendix A , voltages eaA and eaB can be evaluated and the current through the link calculated as i f the slow and fast subnetworks were completely decoupled. The internal node voltages are then updated taking into account the link current. Referring to appendix A , eq. (A.4) shows the nodal and link equations of the partitioned matrix depicted in figure A. 1. A A A A B] B2 £3 ax A GAU GA\\2 0 GAU 0 0 0 0 VAX K A GA2l G All G A l l 0 0 0 0 0 VA2 kAl A 0 GA31 GA33 GA34 0 0 0 1 VA3 hA3 A GAM 0 GAA1 GA44 0 0 0 0 VAA = kA4 5, 0 0 0 0 GBU GB\\1 GB\\3 -1 VB\\ K B2 0 0 0 0 GB1\\ GB11 GB23 0 VB1 kB2 B, 0 0 0 0 GB1\\ G B i l G B i i 0 VB3 hB3 «, 0 0 1 0 -1 0 0 ~z«i 0 (A.4) Rewriting eq. (A.4) using partitioned matrices results in 41 A B a A A 0 P VA K B 0 B q VB K a P' q' —z ia 0 (A.5) where p and q are the connectivity matrices By following up the derivation presented from eqs. (A.7) to (A.9), eq. (A.5) can be rewritten A B a A 1 0 a VA eA B 0 1 b VB eB a 0 0 za la (A.10) where a = A~]p b = B~lq e, = A~xhA A A (4.16) eB=B-% Za = p'a + q'b + z ea=p'eA+q'eB From eq. (A. 10), it is possible to verify that the internal node voltages [v^ ] and [vB ] can be calculated once the link current ia is obtained. It is also clear that for the evaluation of ia, the Thevenin equivalent representation shown in figure 4.6 is accurate. Comparing eq. (A. 10) with figure 4.6, the following equivalences may be derived: 42 RthsloW =Pa Rt¥ast = q'b (4.17) e, aA = PeA e. = qet B As soon as the global solution is obtained at, for example, t = AT, one can go one step further and calculate eaA at t = 2AT. With the values of eaA at t = AT and t = 2AT, the contribution of the slow subnetwork when solving only for the fast subnetwork may be obtained through an interpolation process of the values eaA (AT) and eaA (2AT). The Thevenin equivalent impedance Rti,siow is maintained constant throughout the interpolation process, since its value depends only on the time step AT used for the slow subnetwork discretization, as verified from the first expression in eq. (4.17). Therefore, the interpolation is performed on the value of the Thevenin equivalent voltage source. With the assumption that the latent subnetwork is varying slowly through the interval AT (otherwise the integration step AT itself would have to be decreased), a linear interpolation for the Thevenin equivalent voltage source is executed between the two consecutive values eaA (AT) and eaA (2AT). A question that arises here is how the equivalent voltage source can be calculated at a yet unknown solution time t + AT. The answer is by simply updating the history sources of the slow subnetwork before incrementing the solution step. As soon as the solution of the slow-varying subnetwork is obtained at t, the history sources for the solution at t + AT can be evaluated. And since the values of independent voltage and current sources are a function of time, the Thevenin equivalent voltage source for the slow-varying subnetwork at t + AT can be evaluated even before the next time instant for the slow subnetwork solution is reached. Figure 4.7 shows how the interpolation procedure works, assuming that AT = 4At. 43 exact value | t t + At t + 2At t + 3At t + AT t Figure 4.7: Thevenin equivalent history source linear interpolation for AT = 4At. 4.2. S T A B I L I T Y A N D A C C U R A C Y OF T H E PROPOSED L A T E N C Y T E C H N I Q U E In order to check for the stability and accuracy of the proposed latency technique for network transients simulation, the same approach adopted in sections 3.2 and 3.3 should be followed here. For a proper analysis of the stability and accuracy of the latency technique, it is important to realize that the network is discretized slightly differently according to the solution time instant. If the time instant corresponds to an integer multiple of the large time step AT, the network must be solved in detail, i.e., both the slow and fast subnetworks must have their solutions updated. Accordingly, at all other time instants, only the fast subnetwork must be solved. For the network elements located in the slow subnetwork and considering a time instant when the complete network must be solved in detail, the discretization is performed exactly in the same way as i f latency were not exploited, albeit with a larger time step. Therefore, no further analysis of accuracy and stability is needed. On the other hand, for the network elements 44 history source exact value located in the fast subnetwork and the same time instant, the discretization is performed following the history source accumulation procedure for resynchronized solution. This approach results in the expanded trapezoidal and backward Euler rules previously described. These expanded rules should be checked for accuracy and stability. Switching to a time instant when only the fast subnetwork must be solved, the network elements contained within this part of the network are discretized also in the same way as i f latency were not exploited. The history terms depend only on the previous fast time step and no further analysis of accuracy and stability is needed. For the solution of the fast subnetwork, the history sources of the slow subnetwork elements are not explicitly determined, since the concept of the Thevenin equivalent voltage source interpolation is applied. Stability then is not an issue in this case, as the linear interpolation for the Thevenin equivalent voltage source is performed between two known values at consecutive large time intervals t-AT and t. Another way of visualizing this is to associate the fact that the linear interpolation corresponds to the application of the trapezoidal integration rule for network discretization. Since the trapezoidal rule is stable, so is the Thevenin equivalent interpolation performed on the slow subnetwork. The accuracy of the linear interpolation procedure is dependent on the correct assumption that the latent subnetwork is varying slowly through the large interval AT. The expanded trapezoidal and backward Euler integration rules are therefore the only techniques in this latency approach that need to be further evaluated for accuracy and stability. 45 4.2.1. ACCURACY OF THE EXPANDED INTEGRATION RULES Both the trapezoidal and backward Euler expanded integration rules will be tested for accuracy. A. EXPANDED TRAPEZOIDAL INTEGRATION RULE The case of an inductance will be investigated. The derivation and the resulting distortion factor are the same for a capacitance, however, as was the case for the conventional trapezoidal integration rule. Equations (4.6) and (4.8) can also be written as vL(t) = ^ iL(t)+ehL(t) (4.18) ehL (t) = -—iL (t-AT)-vL (t-At)-2^vL (t-AT + kAt) (4.19) Af k=\\ Substitution of eq. (4.19) into eq. (4.18) yields 2 T ? T \"_1 vL (0 = T 7 ^ (O—TTk {t-^)-^ (t-AT)-2^vL (t-AT + kAt) (4.20) Applying an input of the form vL (t) = eJC0t and assuming the output to be of the form iL (t) — Ye(co)eJ0}t, where Ye(cd) is the admittance in the frequency-domain, leads to ejM+eM,-ATh2^eJ^'-AT+kAt)= — Ye(a))ejM-^Ye(co)eM'-AT) (4.21) k=\\ A/ &t Cancellation of the common factor ejM and multiplication by ejwAT results in 46 e J o A T + l J«*T _ i + -2L «-i JcokAl JaAT . (4.22) 7e(£y) may also be written as Ye(co)^Ye^(co)+Ye(2)(co) (2), (4.23) where 7 ^ ) = — v ' 2L JmnAt + 1 jconAt -1 (4.24) and 7e ( 2 )(/y) = At_ ~2~L n-\\ JatkAt JojnAt (4.25) Both (<2>) and 7e ( 2 ) ( cos CO k — { In j k=\\ I 2 J n and (4.30) The frequency-dependent distortion factor Ke^ (co) can then be obtained as sin Ke{2)(co) = n coAT ^coAT^\"-1 cos k=\\ CO f n\\AT k (4.31) The summation term in the denominator of eq. (4.31) takes into account all the contributions of the small time step At solutions within the large time step AT. Figure 4.8 shows the 48 correspondence between the frequency-domain representation of an inductance and a capacitance in continuous-time and in discrete-time. circuit element in continuous-time circuit element in discrete-time with the trapezoidal rule L e ( 1 ,(o} I_e(2)(c4 c<1>(o} Figure 4.8: Frequency domain representation of inductance and capacitance using the expanded trapezoidal integration rule for latency exploitation. Le® (co) and Le(2) (co) are given respectively by eqs. (4.27) and (4.30). Meanwhile, ,(2) Ce®(co) and Ce(2) (co) can be calculated in similar manner: , 2) tan ( coAT\\ Ce®(co)=nC-' coAT^ [ 2 , (4.32) sin ( coAT^ Ce{2)(co) = nC f coAT V -* J cos k=\\ CO V J AT (4.33) 49 The effect of latency exploitation in network simulation accuracy can be made clearer by some practical examples. The frequency response of the equivalent admittance Ye(co) of the inductance Le(co) using latency can be compared to the exact continuous-time admittance Y(co). 1 1 - + -1 Ye(co)_ jmLe (co) _ jCoLe® (co) jcoLe™ (co) _ Ke® (CD) + Ke{2) (co) _ 1 Y(0))~ _ J _ J~ ~ ~ jCOL jCOL Ke{x)(co)Ke{2)(co) Ke(co) (4.34) Figure 4.9 shows the frequency response of eq. (4.34) in amplitude and phase, assuming AT = 2At. 10 3 C D E 1/Ke<1>(0) 1/Ke<2>(oj I I 1/K e (o) 0.2 0.4 0.6 0.8 0.95 4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 f N y (AT)=0 .5 f N y (A t ) frequency (per unit of Nyquist f requency for a time step At) / / „ 1/Ke<1>(r4 \\ \\ :1/Ke(2)(f4 = 1/Ke(f4 0 0.5 1.0 frequency (per unit of Nyquist f requency for a time step At) Figure 4.9: Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 2At (partial and complete responses) It is interesting to notice that even though the Nyquist frequency of the frequency response component associated to Ke® (co) is located at half the global Nyquist frequency, both the 50 frequency response of the magnitude and phase of eq. (4.34) follow the same pattern of the conventional integration rule. This guarantees the accuracy needed for the fast-varying elements at the time instants when the complete network must be solved in detail. Figure 4.10 shows once again the frequency response of eq. (4.34), focusing only on the complete response however, for ease of comparison with the normal trapezoidal integration rule. 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 f requency (per unit of Nyquist f requency) f requency (per unit of Nyquist f requency) Figure 4.10: Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 2At (complete response) Comparing figure 4.10 with figure 3.5, which gives the frequency response of the conventional trapezoidal integration rule, it is readily observed that the history source accumulation procedure for latency exploitation has not introduced any difference to the frequency response. As an extra example, figure 4.11 presents the frequency response of eq. (4.34) when adopting AT = 4At. 51 Figure 4.11: Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 4At (partial and complete responses) Once again the combined frequency response of the expanded integration rule (l/Ke{co)) in magnitude and phase follows the frequency response of the conventional trapezoidal integration rule. However, in this case, the result obtained for the expanded trapezoidal rule is not analytically easily verified for f = 0.5pu, since both subcomponents (l/Ke^'Xco) and l/Ke^2\\co)) tend to infinite for this value of frequency. The correct result though, can be verified by a detailed analysis of eq. (4.34), which is rewritten here for completeness. 1 Ke{l)(co)+Ke{2)(co) - w K J (4.34) Ke(co) Ke{l)(co)Ke{2)(co) where Ke^ (co) and Ke^ (co) are respectively given by eqs. (4.28) and (4.31). For a value of n = 4, eqs. (4.28) and (4.31) can be rewritten, respectively as 52 tan f coAT \\ Ke{x)(co) = A-f coAT} v 2 j (4.35) f coAT ^ sin Ke(2)(co) = 4 coAT sin j[cos (-coAt)+1 + cos (coAt)] j[l + 2 cos (ffi-Ar)] (4.36) Substituting AT by nAt and after some manipulations, eq. (4.34) results in 1 coAt Ke(co) 2 1 2cos(ffi(Af). s\\n(2coAt) sin(2coAt) tan(2ffltAf) (4.37) As coAt tends to 7t/2, the functions sin(2coAt), tan(2coAt), and cos(coAf) tend to zero, and it is necessary to evaluate the limit of the right-hand side of eq. (4.37) in this case. lim 1 coAt lim VIA, Ke(co) 2 °^72A, 1 ^2cos(coAt)^ 1 sin(2coAt) sin(2ffiA/) tan(2coAt) (4.38) The individual limits of the right-hand side terms of eq. (4.38) are lim 2 At sin (2coAt) • — -f-oo (4.39) lim '2At tan(2coAt) 2cos(coAt) lim = lim '!& sin (2coAt) '2At sin (coAt) (4.40) (4.41) The individual limits of eqs. (4.39) and (4.40) cancel out since sin(oc) = tan(a) for a small. Therefore, eq. (4.38) can now be rewritten as 53 1 coAt K l i m / — - = = — °^7iMKe{(o\\ 2 4 This is exactly the value depicted in figure 4.11 for Ke(co) when co = 7t/2At = 0.5cuNyquist-Figure 4.12 shows again the frequency response of eq. (4.34), focusing only on the complete response for ease of comparison to the normal trapezoidal rule. frequency (per unit of Nyquist frequency) frequency (per unit of Nyquist frequency) Figure 4.12: Frequency response in magnitude and phase of the expanded trapezoidal rule for AT = 4At (complete response) Once again, the frequency response of the expanded trapezoidal rule in both magnitude and phase for the case when AT = 4At matches exactly the frequency response of the conventional trapezoidal integration rule with a single integration step. It is now convenient to analyze the case when latency is exploited using the backward Euler integration rule for history source accumulation instead of the trapezoidal rule. 54 B . E X P A N D E D B A C K W A R D E U L E R I N T E G R A T I O N R U L E One more time the case of an inductance will be investigated. The derivation and the resulting distortion factor are the same for a capacitance, however, as was the case for the conventional backward Euler integration rule. The only change in case of a capacitance occurs in the value of the constant resistance, as was also the case for the conventional rule. Equations (4.12) and (4.14) can also be written as ^(0=^(0+^(0 ehL(t)^-^-iL(t-AT)-^vL(t-AT + kAt) At k=\\ (4.43) (4.44) Substitution of eq. (4.44) into eq. (4.43) results in n-1 vL(0 = ^ iL(0-^h(t-*T)-Y,vL(t-AT + kAt) (4.45) Applying an input of the form vL(t) = eJ'ax and assuming the output to be of the form iL (t) = Ye{co)ejox, where Ye(cd) is the admittance in the frequency-domain, leads to +YeJ^-AT+kAt) =—Ye(oj)e^-±Ye(co)ej<'-^ ^ At K J At (4.46) k=\\ Cancellation of the common factor eJM and multiplication by e j a A T results in n-1 Jwk&t Ye(co)- A/_ L jca\\T JwAT -1 A H L k=\\ eJa*T-l (4.47) 55 Ye(co) may also be written as Ye(co) = Ye{1)(co)+Yel2)(co) (4.48) where Ye^(co)-. At e J a A T - l (4.49) and 7e ( 2 ) (co) At n-\\ x< k=\\ JcokAl (4.50) Both Ye^ (co) and Ye(2) (co) depend on the ratio AT/At. However, the number of terms in ,(2) Ye^ (co) increases as n grows, exactly as shown for the trapezoidal rule. Ye^ (co) can always be written as YeV(a>) =—+ At At 1 2L 2L . 4 (COnAt^ J tan v 2 j (4.51) with i?e^ = — (constant resistance) ( - ^ in case of a capacitance) Af . (4.52) tan ( conAt ^ and Le(l)(co) = L-tan f coAt^ = L-tan ' coAT^ v 2 y f coAT^ v 2 \" y = nL-f coAT^ v ^ y (4.53) The frequency-dependent distortion factor Ke® (co) can then be obtained 56 tan ( coAT\\ Ke{l)(co) = n-' coAT ' { 2 (4.54) The distortion factor Ke^ (co) obtained for the expanded backward Euler integration rule is the same as the one obtained for the expanded trapezoidal integration rule. For the equivalent integration rules with single integration steps this relationship is also verified. Careful examination of eq. (4.50) results in the observation that Ye^ (co) for the expanded backward Euler rule is exactly the same as Ye^ (co) for the expanded trapezoidal rule. Ye^(co), Le{2)(co), and Ke{2) (co) are then given by eqs. (4.29), (4.30), and (4.31) respectively, repeated below for completeness. Ye{2)(co)--n-1 2J y cos CO At 2L jsm COnAt (4.29) V - j . f coAT ^ sin Le{2)(co) = nL-( coAT^ n-1 f f \"1 AT' > COS CO k — { 2 ; k=\\ I 2 ; n (4.30) Ke{2)(co) = n sin coAT f coAT^ 2, cos k=\\ CO L V f , n\\AT k (4.31) Figure 4.13 shows the correspondence between the frequency-domain representation of an inductance and a capacitance in continuous-time and in discrete-time. 57 circuit element in continuous-time circuit element in discrete-time with the backward Euler rule Le(1>(o} U2>(cc) Re = 2L_ t Re = At/2C A A A C<1>(o} Figure 4.13: Frequency domain representation of inductance and capacitance using the expanded backward Euler integration rule for latency exploitation Le-' (co) and Le- ' (co) are given respectively by eqs. (4.53) and (4.30). Ce-' (co) and Ce^ (co) are also unchanged when compared to the expanded trapezoidal rule and their values are given by eqs. (4.32) and (4.33) respectively. Once again, the effect of latency exploitation in network simulation accuracy employing the backward Euler rule will be made clearer by a practical example. The frequency response of the equivalent admittance Ye(co) using latency can be compared to the exact continuous-time admittance Y(co). 1 ( 1 1 , 1 , 1 Ye(co)_Re jcoLe(co) _ Re jcoLe{x)(co) jcoLe{2)(co) Y(co) jcoL jcoL 58 Ye (co) _ jcoL Ke(l) (co)+Ke{2) (co) Y(co) ~~Re~ Ke(])(co)Ke(2)(co) (4.55) Figure 4.14 shows the frequency response of eq. (4.55) in amplitude and phase, assuming AT = 4At. S. 2 1/K e ( 1 ) (o} i 1/K<1)(c4| 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 frequency (per unit of Nyquist frequency for a time step At) frequency (per unit of Nyquist frequency for a time step At) Figure 4.14: Frequency response in magnitude and phase of the expanded backward Euler rule for AT — 4At (partial and complete responses) The frequency response pattern of the magnitude and phase of eq. (4.55) follows the one from the normal backward Euler integration rule with single integration step, as was already the case for the trapezoidal rule. Therefore, accuracy is also guaranteed for the fast-varying elements at the time instants when the complete network must be solved in detail. The same limit analysis performed for the expanded trapezoidal rule for AT = 4At can be repeated here to verify that the magnitude and phase values for f = 0.5pu are indeed correct. Figure 4.15 shows once again the frequency response of eq. (4.55), focusing only on the complete response. This facilitates the comparison with the regular backward Euler integration rule. 59 0 0.2 0.4 0.6 0.8 1.0 f requency (per unit of Nyquist f requency) - - • 71/2 0 0.2 0.4 0.6 0.8 1 0 f requency (per unit of Nyquist f requency) Figure 4.15: Frequency response in magnitude and phase of the expanded backward Euler rule for AT = 4At (complete response) Comparing figure 4.15 with figure 3.6, which gives the frequency response of the conventional backward Euler integration rule, it is easily observed that the history source accumulation procedure for latency exploitation with the backward Euler rule is exactly the same as the frequency response of the regular rule. The phase distortion is maintained and therefore the damping capabilities of the backward Euler rule are preserved when applying the history source accumulation technique. In the upcoming section the stability of both expanded integration rules will be verified. 60 4.2.2. STABILITY OF THE EXPANDED INTEGRATION RULES Stability of the expanded trapezoidal and backward Euler integration rules will be verified in this section. Not only will the two integration rules be proven as being stable, but also the poles and zeros of the equivalents integrator and differentiator will be shown to be the same as in the regular single time step rules, independent of the ratio AT/At. A. EXPANDED TRAPEZOIDAL INTEGRATION RULE The frequency dependent admittance Ye(co) obtained when an inductance is discretized with the expanded trapezoidal integration rule is repeated below, in terms of the small time step At. jamAt + 1 jomAt -1 + -2L n-1 jcokAt jconAt -1 (4.56) Replacing eJ0At byz in eq. (4.56) results in Ye{co)-2L n-1 z\"+i+25y k=\\ zn-\\ (4.57) Equation (4.57) can be expanded, yielding Ye{co) = At (z + l ) ( z + z 2 ) ( z + z 3 ) _ At \"z+r 2L _ ( z - l ) ( z + z 2 ) ( z + z 3 ) ~2L (4.58) 61 This is the same transfer function obtained for the single time step case of the trapezoidal rule, given by eq. (3.29). Equation (4.58) has a pole at z = 1 and a zero at z = -1 and the transfer function is stable as both an integrator and differentiator. However, like in the single time step case, the pole at z = -1 for the equivalent differentiator introduces numerical bounded oscillations at discontinuities. B. EXPANDED BACKWARD EULER INTEGRATION RULE The frequency dependent admittance Ye(co) obtained when an inductance is discretized with the expanded backward Euler integration rule is repeated below, in terms of the small time step At. Ye(o))- At L JwnAt jwnAt At L n-1 k=\\ jwkAt jwnAt (4.59) Replacing eJoAl by z in eq. (4.59) results in n-1 Z +2JZ k=l zn-\\ (4.60) Equation (4.60) can be expanded resulting in At z (z + z2)(z + z3) _ At z T ( z - l ) ( z + z 2 ) ( z + z 3 ) - Z - l -(4.61) 62 This is the same transfer function obtained for the single time step case of the backward Euler integration rule, given by eq. (3.30). Equation (4.61) has a pole at z = l and a zero at z = 0, and the transfer function is again stable as both an integrator and a differentiator. The pole at z = 0 causes the system to be critically damped as a differentiator, avoiding numerical oscillations. This property suggests the possibility of creating a C D A version for latency exploitation. Its implementation, however, is beyond the scope of this thesis work. 4.3. EIGENANALYSIS F O R L A T E N C Y SUITABILITY Before proceeding on to presenting the results of network simulation with latency exploitation, it is important to describe one of the limitations that this method has. In this section, it will be shown that networks suitable for latency exploitation have to be chosen with care, particularly in cases where the partitioning is performed within a lumped area of a network. It is often not true that individual system components can be characterized in an absolute way as being either slow or fast [37]. The eigenvectors associated with such modes give an indication of by how much these modes participate in each state variable. When modelling a system by its state equations, it is straight forward to compute the system's eigenvalues. Latency may then be exploited i f the state variables and eigenvalues of the slow subnetwork are virtually independent from the state variables and eigenvalues of the fast subnetwork. A brief review of state-space representation is performed in appendix B together with a systematic procedure for assigning state variables and writing the dynamical equations for generic lumped linear R L C networks which may contain independent voltage and current sources. In this way, it will be eventually possible to measure the association between the state variables and the modes. This measurement will be of significant importance in determining whether a network may be decoupled and solved with distinct integration steps. 63 The lumped network shown in figure 4.16 will be tested for latency suitability. This network is a good candidate for latency exploitation, since inductor L i and capacitor Ci have a resonant frequency ten times smaller than inductor L2 and capacitor C2. , C 2 = 1|oF Figure 4.16: Network studied for latency suitability Following the state variables identification procedure described in appendix B, the normal tree is identified by the heavy lines in figure 4.16. The voltages across the capacitors in the normal tree and the current through the inductors in the links are assigned as state variables, therefore resulting in four state variables, also depicted in figure 4.16. xi = v c i x2 — ia Xi ~ VC2 X4 ~ lL2 The state space equations may be written as 64 • -r d V a • - — • lc\\ ~<~yi ^ ^ X> ~ C l a V n k dt 1 x, = — V • _ r d v C 2 • 1C2 ~ ^2 J j L ^ X 3 — lC2 c 2 (4.62) dt di 1 VL2 = L 2 ^ T = > X 4 =T~ VL2 at L2 Applying Kirchhoff s current and voltage laws: xx — (x2 x 4 ) x2=Y(~X^+V^ -f-i . _ J _ x 3 — —-x 4 c 2 X 4 ~~~ f^ X| X 3 ^ X 4 ) Z 2 (4.63) Rearranging eqs. (4.63) in matrix form and considering only the response to zero input: 0 1 0 1 X , x 2 1 0 0 0 X , x 3 x4_ 0 0 0 1 _x4 1 0 1 R _L2 L2 (4.64) Replacing the parameters by their actual values shown in figure 4.16, the following genvalues are obtained: 65 \\ =(-0.05 + yi.0038)xl0 6 \\ =(-0.05- jl.0038 ) x l 0 6 ^ =(-0.000004998 + y0.0995)xl0 6 A4 = (-0.000004998-y0.0995)xl0 6 (4.65) The simple determination of the eigenvalues of the circuit is not enough to establish a clear definition on whether latency may be applied to the system. The problem is that the rate of change of each state variable is a linear combination of all the eigenvalues. The best alternative to identify the relationship between the state variables and the modes is to evaluate what is called the participation matrix P [44], which combines the right and left eigenvectors as a measure of the association between the state variables and the modes. with Pn (4.66) where Pu ~u¥n~ Pli = JC, = 1 . ~~p7lc\\ = k dt => x 2 1 lC2 = c2 dvC2 dt -=> x 3 I . _ ~771C2 VL2 = L2 diL2 dt =>x4 1 lC3 = C 3 dvCi dt => x 5 1 . = 4 diLi dt =>x6 1 (4.71) Applying Kirchhoffs current and voltage laws: •*i ( x 2 x * ) C, *2 = — ( - ^ 1 + V , ) X 3 — ( X 4 X 6 ) x 4 = — (x, — x3 — i?,x 4) L 2 . _ J _ •\"•5 — Q X6 x 6 = — (x 3 — x 5 — R2x6) (4.72) 69 Rearranging eq. (4.72) in matrix form and considering only the response to zero input: 0 1 ~c~x 0 1 ~Cx 0 0 X, 1 0 0 0 0 0 X, x 2 0 0 0 1 c2 0 1 x 2 x 4 x5 1 T2 0 1 T2 0 0 x 4 x 5 0 0 0 o • 0 1 0 0 1 0 1 R2 A . (4.73) Replacing the network parameters by their real values shown in figure 4.17, the following eigenvalues are obtained: \\ =(-0.049999 + yl.416872)xl0 6 ^ =(-0.049999 -jl .416872)xl0 6 /L = (-0.019417+ /0.116982)xl0 6 ^ V ' (4.74) /t4 = (-0.019417-y0.116982)xl0 6 A5 =(-0.030584 + y0.051016)xl0 6 \\ = (-0.030584-y'0.051016)xl06 As performed in the previous case, the participation matrix P should be evaluated. Starting from the matrix of right eigenvectors: 70 4 A, 0.0035 + ./0.0004 0.0035-70.0004 -0.0367 + 70.0846 X , -0.0002 + ;0.0025 -0.0002-70.0025 -0.7547-70.1883 x 2 -0.7017-7O.O8OI -0.7017 + 7O.O8OI -0.0025 + 70.0328 xi 0.0743-70.4951 0.0743 + 70.4951 0.1641 + 70.4051 x 4 0.0035 + 7O.OOO4 0.0035 -7O.OOO4 0.0315-70.0196 x5 -0.0743 + 70.4951 -0.0743 -70.4951 0.1679 + 70.4060 x6 K -0.0367-7'0.0846 0.0276-70.0243 0.0276 + yO.0243 x \\ -0.7547 + 7O.I883 0.5900 + 7O.I88I 0.5900-70.1881 x 2 -0.0025-70.0328 -0.0120-70.0505 -0.0120 + 70.0505 x 3 0.1641-70.4051 0.5503 -7O.O274 0.5503 + yO.0274 x 4 0.0315 + 70.0196 -0.0514-70.0765 -0.0514 + 70.0765 *5 0.1679-7'0.4060 0.5474-7'0.0283 0.5474 + y0.0283 x6 The matrix of left eigenvectors is given by 0.3503-70.0527 0.3503 + 7'0.0527 -3.5338-7'4.9570 -3.5338 + 74.9570 -2.6416 + 70.9305 -2.6416-70.9305 -0.0005-7O.OO25 -0.6971 + 70.1046 \\ -0.0005+ 70.0025 -0.6971-70.1046 ^ -0.3636 + 70.3624 -0.0177-70.0126 A$ -0.3636-70.3624 -0.0177 + 70.0126 A4 0.3625 + 70.3005 -0.0059 + 70.0390 A5 0.3625-70.3005 -0.0059-70.0390 \\ 0.0567 + 70.4965 0.0567-70.4965 0.2849 + 70.0453 0.2849-70.0453 0.3958 + 70.1373 0.3958 -70.1373 0.3485-70.0523 0.3485 + 70.0523 0.0343 + 72.4453 0.0343-72.4453 1.4486 + 7'6.8553 1.4486-76.8553 -0.0567-70.4965 -0.0567 + 70.4965 0.2867 + yO.0435 0.2867-70.0435 0.3940 + yO.1358 0.3940-70.1358 4 A, (4.76) And finally, the participation matrix P can be computed using eq. (4.67): 71 0.001256 + y0.000045 0.001256 -./0.000045 0.000006 + 7O.OOOOOI 0.000006-7O.OOOOOI 0.497500 + yO.017555 0.497500-/0.017555 0.250003 - 70.008822 0.250003 + y0.008822 0.001244 + JU000044 0.001244-70.000044 0.249991-70.008823 0.249991 + 7'0.008823 K 0.549100 + 70.117230 0.342619 + 70.205052 0.000457 + 70.000548 0.028401-70.122851 0.048940-70.076261 0.030483-70.123719 0.549100-70.117230 0.342619-70.205052 0.000457-70.000548 0.028401 + 70.122851 0.048940 + 70.076261 0.030483 + 7'0.123719 -0.050356-70.090034 0.157375-70.245457 0.002043 + 70.000167 0.221596-70.064686 0.449817 + 7'0.463155 0.219526-70.063144 -0.050356 + 70.090034 0.157375 + y0.245457 0.002043-7'0.000167 0.221596 + 70.064686 0.449817-70.463155 0.219526 + 70.063144 (4.77) From a careful analysis of the participation matrix, it can be observed that no clear independence between fast and slow behaviour can be observed for some state variables. This is the case, for example, of state variable which is equally dependent on the pair of eigenvalues X\\ and A2 (very fast) and on the pair X5 and Xs (slow). Latency exploitation in the way proposed in figure 4.17 would then not guarantee very accurate results as will be shown in chapter 5. In power systems, however, the separation of time constants is usually more explicit than in general lumped R L C networks. This is the case, for example, of an H V D C converter station, where, for an accurate simulation, the switching frequency of the valves requires an integration step much smaller than the power network connected to the converter. Synchronous generators eletromechanical transients are much slower than electromagnetic transients and the decoupling of time constants is again well defined in this case. It should also be mentioned that lossy frequency dependent transmission lines and power transformers restrict the propagation of high frequency components usually to adjacent buses to the fault location. Therefore, while small time steps may be necessary for an accurate transient simulation near the fault location, larger time steps may be used for the feeding networks. 72 The eigenanalyses presented in this section have been performed based on a continuous-time description of the electrical networks. However, time-domain simulation is based on a discrete-time description. Indeed, the eigenvalues of the discrete-time system are not exactly the same as the eigenvalues of the continuous-time system. They are rather related by the well known bilinear transformation, corresponding to the application of the trapezoidal rule for the discretization of the differential equations [45]. 2 z -1 4= T~M (4-78) where the set of At are the continuous-time and the set of z, are the discrete-time eigenvalues. In case the backward Euler is used for discretization of the differential equations, the transformation relating continuous-time and discrete-time eigenvalues is given by 1 z -1 4 = 7 - — ( 4 - 7 9 ) At z.. However, as also shown in [45], the eigenvectors of the continuous-time and discrete-time systems are the same. This results from the fact that the state space matrix is also mapped from continuous-time to discrete-time according to either the bilinear transformation, given by eq. (4.78) or the transformation given by eq. (4.79), i f the differential equations are discretized with the trapezoidal or the backward Euler rule, respectively. Since the eigenvectors of the continuous-time system are equal to the eigenvectors of the discrete-time system, the participation matrix P remains unchanged and the approach described in this section is valid. 73 CHAPTER 5. CASE STUDIES The methodology proposed for latency exploitation, as well as the motivation to pursue a better simulation performance in terms of processing time have been explained in the previous chapters. It is now essential to present simulation results for some typical circuit configurations employing latency. The results are compared with a typical EMTP simulation where a single time step, necessary to track the fastest frequency components present in the network, is adopted. This comparison, eventually, validates the latency technique applied. In section 5.1, latency is tested on networks consisting solely of lumped elements. Section 5.2 shows the results for networks containing transmission lines, modelled as ideal lossless lines. Finally, in section 5.3 the results obtained for networks modelled with the cp-line model are presented. 5.1. L A T E N C Y E X P L O I T A T I O N IN N E T W O R K S CONSISTING S O L E L Y O F L U M P E D E L E M E N T S In this section, latency is employed for the time-domain simulation of electric networks containing lumped elements only. Three different networks have been chosen, and each one will be described in detail. 5.1.1. S I M U L A T I O N A The first circuit simulated with two distinct integration steps is a simple lumped circuit proposed in [38]. Figure 5.1 shows the circuit under test. 74 L = 100mH R, = 40QQ fast / i slow A / W r\\j) Vs = cos(ol) C_= lOCpF Figure 5.1: Lumped circuit for latency exploitation (case A) With the element parameters indicated in figure 5.1, the two poles of the system are negative real numbers. This condition is necessary for this particular circuit i f latency is to be exploited, since the only other possibility would be to have a pair of complex conjugate poles. Clearly, latency would be of no advantage in this case since both state-variables would be varying at the same speed. However, the case of a non-oscillatory response with two distinct real poles can illustrate the benefits of latency exploitation. The voltage source is applied at t = 0 and three different simulation options are tested, all using the trapezoidal integration rule. The time constant associated with the fast part of the circuit is x fast = L/ = 250 us and the time constant associated with the slow part is (1) Standard procedure using a small time step for the complete network solution in order to guarantee maximum accuracy: At = 25us (normal EMTP solution); (2) Dual step sizes: At = 25us for the fast part of the circuit and AT = 250us for the slow part (latency approach); ^ = ^ = 2500 us. 75 (3) Large time step for the complete network solution: AT = 250LIS (normal EMTP solution). Figure 5.2 shows the voltage across the inductor for the three methods proposed, while figure 5.3 shows the current through the inductor for the same three methods. Method 1: Smal l t ime-step for the complete circuit > 1 © 0.5 cn S o o > - 0 . 5 a> 0.5 cn CO = 0 o -0.5 -0.5 i i p i i r 0 1 2 3 4 5 6 7 8 9 10 Method 2: Apply ing latency i i i i i 0 1 2 3 4 5 6 7 8 9 10 Method 3: La rge t ime-step for the comple te circuit • i i i i i i • ! i ! ! 1 i i i 4 5 6 time (ms) 10 Figure 5.2: Voltage across the inductor for the three methods proposed Method 1: Smal l time step for the complete circuit 0 1 2 3 4 5 Method 2: Applying latency Method 3: Large time step for the complete circuit 4 5 6 time (ms) Figure 5.3: Current through the inductor for the three methods proposed 76 This is an opportune moment to introduce a conceptually simple modification performed in the way the transient simulation is initialized when exploiting latency. Although simple, this modification is very important for a correct implementation of the latency methodology. Usually, in EMTP-type simulators, steady-state solution, in case there is one, or simple zero initialization is performed at t = 0 and the first transient solution is obtained at / = At. When latency is exploited, this approach has to be modified due to the requirements of subnetwork solutions resynchronization at t = AT. If the approach followed by EMTP were adopted, both slow and fast subnetworks would be initialized at t = 0. However, the first transient solution for the fast subnetwork would be calculated at t = At, while for the slow subnetwork the first solution would be calculated at t = AT = nAt. Therefore, the first solutions would already be desynchronized and, potentially, inaccurate results could develop. To avoid this desynchronization, the initialization is performed at / = -At for the fast subnetwork and at t — —AT for the slow subnetwork. Each subnetwork has then its first transient solution calculated at / = 0 and the solutions are synchronized from the very beginning. Following their own integration steps, the complete network solution would be resynchronized at t = AT,2AT, and so on, as it should be. Figure 5.4 shows once again the voltage across the inductor for the three methods proposed including the initialized value at a negative time. For sake of comparison, methods 1 and 3 are also initialized at t — -At and t = -AT, respectively. 77 Method 1: Small time step for the complete circuit CD CJ) o > 1 0.5 0 -0.5 -0.5 0 0.5 1.0 Method 2: Applying latency 0.5 1.0 time (ms) 1.5 2.0 0 0.5 1.0 1.5 2.0 Method 3: Large time step for the complete circuit 2.0 Figure 5.4: Voltage across the inductor up to 2ms and including the initialized value From figure 5.4 it is clearly visible that method 3 cannot track the fast surge on the voltage across the inductance at the start of the simulation. While methods 1 and 2 give very similar results (almost 1 .OV) at t = 0, method 3 predicts this value erroneously (approximately 0.6V). This is just a simple example of what an unsuitable time step may cause to the results of a simulation. After the initial transient, method 3 gives similar results to methods 1 and 2. However, i f the mode most related to the inductance is re-excited, simulation results from method 3 will be compromised once again. The current through the inductor shown in figure 5.3 also presents some difference at the beginning of the simulation, when comparing methods 1 and 2 with method 3. The latency method, however, provides results that are virtually indistinguishable from the base method (method 1) where a very small time step is employed for the complete circuit. Figures 5.5 and 5.6 show the voltage across and the current through the capacitor, respectively. 78 Method 1: Small time step for the complete circuit 9 10 Method 2: Applying latency ~i 1 1 1 1 1 1 1 r Method 3: Large time step for the complete circuit time (ms) Figure 5.5: Voltage across the capacitor for the three methods proposed The voltage across the capacitor does not present fast transients at any instant of the simulation. A l l three methods, therefore, give virtually the same results. Method 1: Small time step for the complete network y—^1 i i i i i i i i 0 1 2 3 4 5 6 7 8 9 10 Method 2: Applying latency S 1 1 I 1 1 1 i i ' : : : : : : : i i i i i i i i 0 1 2 3 4 5 6 7 8 9 10 Method 3: Large time step for the complete network y- ~ J I I I I I i i : : : : ! i ^ ^ - - L ^ j ; i \\ J ^ . — \" i i i i i i i i i 0 1 2 '3 4 5 6 7 8 9 10 time (ms) Figure 5.6: Current through the capacitor for the three methods proposed 79 The current through the capacitor will also be shown with the initialized value at a negative time. Figure 5.7 presents the simulation results in this case. Method 1: Small time step for the complete circuit „ 2 < r 1 c £ 0 *_ o -1 i — T — 1 • i y \\ \\ \\ i t i i I i -0.5 0 1 2 3 4 p Method 2: Applying latency < r 1 c £ 0 i_ 3 O -1 i i i _ J / \\ \\ 1 i i i i i -0.5 0 1 2 3 Method 3: Large time step for the complete circuit < E 3 o Figure 5.7: Current through the capacitor up to 4ms and including the initialized value Although the current through the capacitance presents a fast rate of change initially, method 2, when compared to the results given by method 1, guarantees more accurate results than method 3. After the current reaches its positive peak, all three methods provide accurate results. Table 5.1 shows the number of floating point operations (FLOPS) performed for a total simulation time of 1 Oms for each of the three methods proposed. 80 Table 5.1: Number of FLOPS for test case A considering a total simulation time of 10ms Type of simulation performed FLOPS Method 1 19,653 Method 2 17,127 Method 3 2,013 The decrease in the number of numerical operations by approximately 15% with the method exploiting latency may not seem very significant. However, this is a very small circuit and the slow subnetwork consists of one single node, while the fast subnetwork contains two nodes and therefore demands the most computational efforts. 5.1.2. S I M U L A T I O N B Latency will now be exploited in the lumped circuit initially shown in figure 4.16, when the discussion revolved around the eigenanalysis of generic lumped networks with the purpose of determining latency suitability. Figure 5.8 shows the circuit under consideration once again. L, = 1 MH f 1.0 V 60 Hz slow V =0.1 CI •AAA-t = 0 = 100 uF Slow L, = 1 uH r 'fast Fast C 2 = 1 uF Figure 5.8: Lumped circuit for latency exploitation (case B) 81 This circuit was initially proposed in [37] and, as the eigenanalysis previously indicated, seems to be a good choice for latency exploitation. Vf a s t is the voltage across the \"fast\" capacitor (l | iF), while Vsiow is the voltage across the \"slow\" capacitor (IOOLLF). A voltage source with a frequency of 60Hz and an amplitude of 1 .OV has been connected for a long time, while the ideal switch closes at t - 0. As in case A, three different simulation methods were proposed and tested using the trapezoidal integration rule, as follows: (1) Standard procedure using a small time step for the complete network solution in order to guarantee maximum accuracy: At = 0.2u.s (normal EMTP solution). (2) Dual step sizes: At = 0.2u,s for the fast part of the circuit and AT = 2.0\\xs for the slow part (latency approach); (3) Large time step for the complete network solution: AT = 2.0(is (normal EMTP solution). Figure 5.9 shows the voltage across the \"slow\" capacitor for the three methods proposed, while figure 5.10 shows the voltage across the \"fast\" capacitor for the same three methods. Method 1: Small time step for the complete network > 1 CD CO S 0.96 o > 0.92 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Method 2: Applying latency > 1 CD c? 0.96 o > 0.92 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Method 3: Large time step for the complete network > 1 CD cn 3 0.96 o > 0.92 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time (ms) Figure 5.9: Voltage across capacitor Ci (\"slow\") 82 Method 1: Small time-step for the complete circuit 2 1.5 1 0.5 0 0 0.5 1 1.5 Method 3: Largel time-step for the complete circuit 0.5 , 1 time (x10\"1ms) 1.5 Figure 5.10: Voltage across capacitor C2 (\"fast\") The voltage across the slow capacitor, as shown in figure 5.9, is accurately predicted by the three different methods. Since the participation factors of the \"fast\" eigenvalues on the \"slow\" state variables are very small when compared to the participation factors of the \"slow\" eigenvalues, it is possible to increase the integration step by a factor of ten and still obtain accurate results. However, the voltage across the fast capacitor, as shown in figure 5.10, is predicted very poorly by method 3. Since the \"fast\" eigenvalues present a very high oscillating frequency, a time step of 2.0us imposes a Nyquist rate too close to the mentioned frequency and, therefore, the results are not accurate enough anymore. Method 2, on the other hand, guarantees very accurate results for the voltage across capacitor C2, to the point where the simulation results obtained with methods 1 and 2 become virtually indistinguishable. Table 5.2 shows the number of FLOPS performed for a total simulation time of 1.0ms for each of the three methods proposed. 83 Table 5.2: Number of FLOPS for test case B considering a total simulation time of 1.0ms Type of simulation performed FLOPS Method 1 405,085 Method 2 257,123 Method 3 40,585 In this case a significant decrease in the number of numerical operations was obtained. From this particular example it is straight forward to evaluate the maximum theoretical decrease in the number of FLOPS when latency is adopted. Since each subnetwork is composed of two nodes, it is virtually impossible to obtain a decrease greater than 50% of the number of FLOPS from method 1. As the reduction obtained is around 40% for a ratio of step sizes equivalent to 10, the latency approach is very efficient. Note that method 3 provides the fewest number of FLOPS, as expected. However, the simulation results may be very inaccurate. From the simulation results presented in figs. 5.9 and 5.10, it is possible to verify that the eigenvalues with the most significant contribution to the voltage across the \"fast\" capacitor result in a very high damping coefficient. The fast transients die out very quickly and the only frequencies remaining are the resonant frequency 6% = l /^/L,C, and the source frequency. The latency behaviour should be verified for the case when the fast transients are very slowly damped. This can be achieved by reducing the value of resistor R by two orders of magnitude (R = 0.001 Q.) in the circuit shown in figure 5.8. Once again, the same three methods will be employed. Figure 5.11 shows the voltage across capacitor C-, while figure 5.12 shows the voltage across capacitor Ci in this case. 84 Method 1: Small time step for the complete circuit 0.2 0.3 0.4 0.5 0.6 0.7 Method 2: Appliying latency 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Method 3: Large time step for the complete circuit 0.2 0.3 0.4 0.5 0.6 time (ms) 0.7 0.8 0.9 Figure 5.11: Voltage across capacitor Ci (\"slow\") for R = 0.001Q Method 1: Small time-step for the complete circuit 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Method 2: Applying latency 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 time (x10\" 1ms) 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Method 3: Large time-step for the complete circuit Figure 5.12: Voltage across capacitor C 2 (\"fast\") for R = 0.001Q Even though the fast transient is very lightly damped, method 2 is capable of accurately predicting the voltage across both the slow and fast capacitors, while method 3 clearly introduces a considerable distortion factor to the voltage across capacitor C2. These examples show that latency exploitation with the trapezoidal integration rule is a very accurate and efficient alternative to the traditional simulation when only a single integration step is used for the whole network. The backward Euler integration rule will now be tested under the same conditions as the ones used for the trapezoidal rule. Figures 5.13 and 5.14 show the voltage across the \"slow\" and \"fast\" capacitors respectively, for the original value of resistor R (R = 0.1Q). Method 1: Small time step for the complete circuit CD CO CD D) CO o > 0.2 0.3 0.4 0.5 0.6 0.7 Method 2: Applying latency 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Method 3: Large time step for the complete circuit 0.4 0.5 time (ms) Figure 5.13:Voltage across capacitor C- (\"slow\") with the backward Euler rule 86 Method 1: Samll time step for the complete circuit 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Method 2: Applying latency 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Method 3: Large time step for the corripte circuit CS 1.5 CD r? 1 CO o > 0.5 _1_ 0.1 0.2 0.3 0.4 0.5 0.6 time (x1Cr1ms) 0.7 0.8 0.9 Figure 5.14: Voltage across capacitor C 2 (\"fast\") with the backward Euler rule As expected, the backward Euler rule introduces extra damping to the circuit under simulation. This is observed in figs. 5.13 and 5.14 when compared to figs. 5.9 and 5.10, respectively. Even in the case that a single and small integration step is used, the damping is very noticeable. For the voltage across the \"slow\" capacitor C i , method 2 presents a much faster damping than method 1. This is not a consequence of the latency technique itself though. The reason is simply that AT is ten times larger than At used in method 1. Therefore, the equivalent resistances in the frequency domain model of the backward Euler discretized elements are also increased by a factor of 10 (see eqs. (3.24) and (3.27)). Indeed, method 2 shows the same damping characteristics as method 3 for the voltage across the \"slow\" capacitor. Meanwhile, methods 1 and 2 predict accurately the voltage across the fast capacitor, albeit with increased damping when compared to the results obtained with the trapezoidal integration rule. Method 3, once again, cannot follow the fast transients at C2. 87 5.1.3. S I M U L A T I O N C The next example will illustrate a case where latency application does not provide very accurate results, due to the impossibility of decoupling the eigenvalues and the state variables. It will be shown, however, that the results obtained with the latency approach are still more accurate than those obtained i f the network were simulated with a single large integration step. The circuit is the same as the one presented in figure 4.17 and is repeated below in figure 5.15. In this circuit, an additional L C cell is added so that the resonant frequency of this cell is the same as in the slow subnetwork of the previous circuit analyzed. It is assumed that the new cell has a \"slow\" behaviour and therefore can be simulated with the same large integration step of 2.0us. i.ov / p . \\ 60 Hz V W / V s l o w l L 1 = 1 uH R, = 0.1 Ci ^ =0.1 ci w v — C , = 1 0 0 uF Slow L, = 1 uM ' f a s t Fast t = 0 C 2 = 1 uF — Fast •4 1-3 = 1 MH slow 2 C , = 100 uF Slow Figure 5.15: Lumped circuit for latency exploitation (case C) (modified from case B) The analysis of the participation matrix for this network presented in chapter 4 has indicated that the voltage across the inductor L 3 , which is supposedly located in a \"slow\" part of the network, is actually very dependent upon the eigenvalues with the highest oscillating frequency, associated with the elements L2 and C2. Figure 5.16 shows the voltage across the inductor L3 for the same three methods used for the simulation of test case B; i.e., At = 0.2us and AT = 2.0us. 88 Method 1: Small time step for the complete circuit -0.5 ' 1 ' ' 0 100 200 300 Method 3: Large time step for the complete circuit 0 100 200 300 time (ps) Figure 5.16: Voltage across inductor L3 Figure 5.16 clearly shows the strong dependence of the voltage across inductor L3 on the fast modes. The assumption that L3 is located in a \"slow\" part of the network is not very accurate. However, latency can still be used in this case to track the average value of the signal. Even i f the high frequency components were damped slower than the current example shows, the latency technique would still be capable of following the slow modes associated with this variable accurately since it acts as a low-pass filter. The decision on whether to exploit latency in this case may be dependent on the simulation needs of the user, and a proper judgement has to be made. The worst that can happen i f latency is exploited when there is no clear definition on the dominant modes for some state variables is not making use of all the latency potentiality. However, results will be more accurate with the latency approach than with a single and large time step because the latency algorithm is absolutely stable and accuracy always increases in those parts of the network where smaller steps are used. 89 Up to now, only electric networks consisting exclusively of lumped elements have been considered. In power systems however, transmission lines usually have to be modelled for transient simulation. The representation of the lines with a characteristic impedance and propagation time is generally more accurate than lumped based models. The application of latency in such environments will be analyzed in the following section. 5.2. LATENCY EXPLOITATION IN NETWORKS CONTAINING TRANSMISSION LINES MODELLED AS IDEAL LOSSLESS LINES The latency technique proposed in this thesis work will now be applied to networks containing ideal lossless single-phase transmission lines. A lossless transmission line means that the line does not introduce any attenuation or distortion between the signals at the sending and receiving ends of the line. The line has a travelling time T solely dependent on its length and the voltage and the current at the receiving end of the line at any given instant t are equal to the voltage and current at the sending end at the instant t-T. Figure 5.17 shows the network studied. This network is derived from a simulation example proposed in [37]. 150 Km R, = 1 £ 2 L,=25mH R 3 = 1 Q AA/V 1 V 60 Hz Slow 1 L 3 = 25 Figure 5.17: Network for latency exploitation in single-phase transmission lines 90 The eigenvalues associated with the fast part of the network have a very high oscillating frequency, around 30kHz. The eigenvalues associated with the two slow subnetworks are negative real numbers with a much smaller magnitude than the complex conjugate ones. The fast transients of this network are excited by applying an initial charge of 1.0V to the capacitor of the fast load. One more time three different simulation procedures are tested. (1) Standard procedure using a small time-step for the complete network solution in order to guarantee maximum accuracy: At = 1 .Ous (normal EMTP solution). (2) Dual step sizes: At = l.Ous for the fast part of the circuit and AT = 20.Ous for the slow part (latency approach); (3) Large time step for the complete network solution: AT = 20.Ous (normal EMTP solution). Figure 5.18 shows the voltage across the inductor on the source side of the network for the three simulation methods proposed, while figure 5.19 shows the voltage across the capacitor for the same three methods. Method 1: Small time step for the complete circuit 0.5 1 Method 2: Applying latency 1.5 0.5 1 1.5 Method 3: Large time step for the complete circuit Figure 5.18: Voltage across the inductor located on the source side of the network 91 Method 1: Small time step for the complete circuit _i_ _ i_ _i_ 0.2 0.3 0.4 Method 2: Applying latency i r~ 0.5 0.1 0.2 0.3 0.4 0.5 Method 3: Large time step for the complete circuit 0.2 0.3 time (ms) 0.5 Figure 5.19: Voltage across the capacitor The latency approach (method 2) is capable of reproducing very accurately the voltage across the capacitor (according to figure 5.19), which is located in the fast subnetwork. Method 3, however, fails to provide accurate results since a time step of 20us is unsuitable to track the fast transients. Regarding the voltage across the slow inductor (figure 5.18), methods 2 and 3 give accurate results for the slow transients, initiated by the connection of the voltage source to the network. Nevertheless, since the transmission lines are modelled as ideal lossless lines, the fast transients excited by the initial charge to the capacitance travel down the 150km-long line and reach the slow part of the network at t = 0.5ms. These fast transients are, therefore, injected into the slow subnetwork. As a result, both methods 2 and 3 are incapable of reproducing the high frequency components that appear on the voltage across the inductor at t = 0.5ms. The situation is similar to the lumped circuits previously analysed, where presumably slow components may also present significant high frequency transients due to the non-decoupling characteristics between the 92 eigenvalues and the state variables. According to figure 5.18, however, method 2 at least guarantees a very similar time constant for the transient envelope to the one observed in method 1 for the oscillations initiated at t = 0.5ms. With method 3, the transient envelope fails to attenuate as fast as with methods 1 and 2. This observation corroborates the previous statement regarding the fact that the latency approach should always provide more accurate results than using a single and large integration step. Also, while in the case of. lumped R L C networks the fast modes may truly influence significantly the response of presumably slow components, the same is not true when- the network is decoupled by long transmission lines. In the more realistic case of frequency dependent lossy transmission lines, the fast transients are damped out during their propagation since the frequency dependent transmission line resistance increases as frequency increases. This is the reason why, as previously mentioned, the high frequency oscillations usually do not propagate far into the network from the point where a transient event occurs. Table 5.3 shows the number of FLOPS performed for a total simulation time of 1.5ms for each of the three methods proposed. Table 5.3: Number of FLOPS for the lossless line model considering a total simulation time of 1.5ms Type of simulation performed FLOPS Method 1 144,223 Method 2 105,903 Method 3 7,328 A decrease of about 30% in the number of FLOPS is achieved for this network when exploiting latency according to method 2. In principle, it appears that this reduction could be more significant. However, even though the fast subnetwork contains only two nodes compared 93 to four nodes in the slow subnetworks, there are two energy storing elements in the fast subnetwork (there are two such elements in the slow subnetworks as well) and the decrease in the number of FLOPS is limited to 30% due to this constrain. 5.3. LATENCY EXPLOITATION IN NETWORKS CONTAINING TRANSMISSION LINES MODELLED WITH THE cp-LINE MODEL The constant-parameter line (cp-line) model [2] is an improvement to the ideal lossless line model. In the cp-line model, the transmission line is no longer considered as lossless. Instead, the losses are lumped at specified locations of the line. In the EMTP, the total line resistance is lumped in three different locations: Vi of it is lumped at the middle of the line and % at each end of the line. For the simulations presented in this section, the total line resistance is divided by two and lumped only at both ends of the line. The purpose of testing the latency technique with transmission lines modelled with the cp-line model is to build a framework for the future latency exploitation with the z-line model for transmission lines. The z-line model [8], as described in appendix C, is essentially a model that eliminates the necessity of synthesizing the frequency dependent transformation matrix in multiphase lines when taking into account the frequency dependence of the line parameters. For this purpose, the frequency dependent lossy section is lumped at different points of the line. The frequency independent section is then modelled as ideal lossless lines with the travelling wave equations. This model combines the simplicity of the cp-line model with the accuracy of the complex frequency dependent line models available in the EMTP where transformation matrices relating phase and modal quantities are, in general, frequency dependent [2]. One of the severe limitations of the z-line model is that the original transmission line has to be partitioned into a considerable number of segments for the lumping of the lossy section. The 94 number of segments must increase as the accuracy desired for higher frequencies increases. In turn, very small integration steps are required, resulting in possibly very time consuming simulations. Latency exploitation may overcome this limitation since the networks connected to both the sending and receiving ends of the transmission line do not have to be simulated with the same small time step required for the simulation of the short internal segments within the line. Figure 5.20 shows the network studied. 3 Q 350mH r J W V i Y Y Y V 230kV (rXj 3 Q 350mH ^ A / V ^ Y Y Y Y . 230kV 1 1 1 100nF i 1 100nF 15km 15km A' t = 0 350mH 230kV yXj) Figure 5.20: Circuit proposed for latency testing with the cp-line model The single phase transmission lines are 15km long, and each line is divided into five 3km-long segments. The 3km segments are modelled with the cp-line model, in the way previously described. The main characteristic of the z-line model is maintained, since each transmission line is split into several segments, therefore limiting the global time step for EMTP-like simulators to at most the travelling time of each short segment. Figure 5.21 presents the transmission line partitioning implemented for this case study. A short-circuit is applied at t = 0 between the two 15km-long transmission lines, as indicated in figure 5.20. This network contains a total of sixteen nodes and nine of them are internal nodes of the transmission lines. The remaining seven nodes are contained within the slow subnetworks. The total internal resistance of each 15km-long transmission line is assumed to be 10Q. 95 15km -AAA Q } W \\ r -R | i n e /10 3km R | i n e /5 3km R | i n e /5 3km R | i n e /5 3km R | i n e /5 3km R | i n e /10 - A A A r - r l )—V\\A/- f~ l )-^ VvV-n / - A A A r - Q M A A r O > A A A r -Figure 5.21: Transmission line divided into five segments for cp-line modelling Two simulation methods are proposed in this case as follows: (1) Standard procedure using a small time step for the complete circuit for maximum accuracy. The time step adopted is equivalent to the travelling time of the 3km-long segments: At = lOps (normal EMTP solution); (2) Dual step size: At = IOLLS for the solution updating of the internal nodes of the transmission lines and AT = 50LLS (equal to the travelling time of the original 15km-long lines) for the external subnetworks connected to the transmission lines (latency approach). A simulation method employing a AT = 50u.s for the complete network solution is impractical in this case, since the time step used for the transmission line simulation cannot be larger than the travelling time of IOJIS required by the short line segments. Figure 5.22 shows the voltage to the right of the first short transmission line segment while figure 5.23 shows the voltage across the 350mH inductors in the external subnetwork for the two methods proposed. 96 Method 1: Small time step for the complete network 2 3 Method 2: Applying latency Figure 5.22: Voltage to the right of the first short transmission line segment Method 1: Small time step for the complete network -100 O) CO 200 100 -100 Method 2: Applying latency time (ms) Figure 5.23: Voltage across the 350mH inductors in the external subnetwork Both methods give the same results for the voltage across the inductors, as a time step of 5 Ous is more than sufficient to simulate the slow dynamics of this part of the network. Method 2 is also very accurate for predicting the voltage of the internal node of the transmission line. The total internal impedance of each 15km-long transmission line is now increased to 50Q to represent an increase in losses due to high frequency transients. Figure 5.24 shows the voltage to the right of the first short transmission line segment. Method 1: Small time step for the complete netw ork cn CO Method 2: Applying latency CD CO s o > time (ms) Figure 5.24: Voltage to the right of the first short transmission line segment (internal line resistance increased to 50Q) With the line resistance increased to 50Q. the high frequency oscillations are attenuated much faster than for a line resistance of \\0£l. The latency method (method 2) continues to provide accurate results. Table 5.4 shows the number of FLOPS performed for a total simulation time of 5ms for each of the two simulation methods proposed. 98 Table 5.4: Number of FLOPS for the cp-line model considering a total simulation time of 5ms Type of simulation performed FLOPS Method 1 166,336 Method 2 121,462 Once again, method 2 considerably reduces the number of floating point operations when compared to method 1. The decrease in the number of FLOPS would be even more significant i f the external networks connected to the ends of the lines were represented in much more detail. 99 CHAPTER 6. CONCLUSIONS AND FUTURE WORK 6.1. C O N C L U S I O N S A new methodology for latency exploitation in time-domain EMTP-type simulators has been presented. The latency concept in this thesis work is related to the possibility of using different integration steps for the numerical solution of electric networks. The network is partitioned according to the M A T E technique, and the resulting subnetworks are discretized with time steps convenient to the accurate solution of each particular subnetwork. The purpose of the latency solution is to allow a decrease in total simulation time, eventually facilitating the achievement of real-time simulation for large power system networks. The latency simulation method has been described in detail, with emphasis on the two proposed techniques for increased efficiency and accuracy of the simulation: the history source accumulation of the fast subnetwork elements for resynchronized solution, and the Thevenin equivalent voltage source interpolation for the slow subnetwork representation. The stability and accuracy of the modified integration rules obtained have been compared to the conventional integration rules. Both stability and accuracy are not affected when latency is exploited. An eigenvalue analysis is proposed for the verification of latency suitability of generic networks. The methodology is based on the evaluation of the participation matrix, which measures the degree of coupling between each eigenvalue and state-variable of the system. Lumped networks as well as networks containing transmission lines, modelled with the travelling waves approach, have been simulated employing the latency method proposed. The transmission lines have been modelled as ideal lossless lines and with the cp-line model. The results obtained when latency is exploited are in good agreement with those obtained when the 100 entire network is simulated with a small integration step. The simulation results have shown that exploiting latency guarantees more accurate results than employing a single and large integration step, even for networks where there is no clear decoupling between the fast and slow eigenvalues and state-variables. The number of numerical operations is significantly reduced, even for small networks, raising high expectations for the possible time savings incurred when latency is applied to large electric systems. Further gains may be forecasted with the code optimization. Some concern may be raised regarding the possible occurrence of aliasing when different parts of a network are simulated with distinct integration steps. Aliasing is not an issue in this latency implementation since the fast transients developing in the fast subnetwork are not injected into the slow subnetwork. For the time instants when only the fast subnetwork is solved, the link current connecting the subnetworks (see figure 4.6) is never used to update the internal slow variables. The fast transients are therefore filtered out from the point of view of the slow subnetwork. The main contributions of this thesis work may be summarized as follows: 1. Development of a latency exploitation technique for time domain transient simulation. Of particular importance are the history source accumulation technique and the Thevenin equivalent interpolation procedure. These two methods allow a compromise between accuracy and efficiency of the simulation, and guarantee a smooth interface between the slow and fast variables. The latency technique developed has also been checked for accuracy and stability. 2. Implementation of the latency technique within the M A T E framework, which propitiates a natural transition to a real-time simulator. In fact, latency is one of the applications for which M A T E was initially conceived and the successful latency implementation also justifies the development of the M A T E technique. 101 3. Achievement of latency exploitation using non-iterative solutions. The robust characteristics of EMTP-type simulators are thus maintained. 4. Evaluation of the participation matrix for latency exploitation suitability in generic R L C networks. Since the state-variables of a network are a linear combination of all its eigenvalues, it is generally not true that the network may be partitioned and solved with different integration steps. If there is not a strong physical evidence of the possibility of exploiting latency for time domain simulation, the evaluation of the participation matrix may be a good alternative in order to verify the degree of coupling between the system eigenvalues and the state variables. If this analysis shows that the eigenvalues of the distinct subnetworks are fairly independent, then latency exploitation should provide very accurate results. 5. Simulation results for some case studies show that latency simulation is a viable alternative for EMTP-type simulators. Combining accurate results with a decrease in computational burden may prove to be a valuable resource to the future of electromagnetic transients simulation. 6.2. R E C O M M E N D A T I O N S F O R F U T U R E W O R K This thesis work provides a first glimpse on the possible application of latency techniques for the time domain simulation of electromagnetic transients. Although the thesis describes in detail a methodology for latency exploitation and the results obtained corroborate its efficiency and accuracy, the author's perspective is that work on latency is far from being exhausted. Some of the main recommendations for future work related to this thesis project are listed below. 102 1. Development of a generalized latency methodology that allows an automatic and efficient network partitioning. The reverse application of the participation matrix as proposed in this thesis may prove to be a good starting point in order to reach this objective. Different network partitions may be performed and the most suitable for latency exploitation would be the one where the participation matrix indicates that state-variables predominantly fast or slow have been decoupled. 2. Implementation of the latency technique into a real-time simulator. This thesis work proposes a latency technique in a framework suitable for real-time simulation, following the M A T E methodology, which is implemented into OVNI, UBC' s real time simulator. The real implementation of the latency technique into OVNI, however, was out of the scope of the thesis. Latency may ideally be suited for implementation into the distributed OVNI version, where a network of inexpensive personal computers is assembled. The simulation task may be assigned to the different workstations in the most efficient way by allowing the fastest subnetwork to be solved by a dedicated workstation. The slower subnetworks are then solved by the remaining units, the exact number depending on the demands of the network and of the simulation. Further information on this topic may be obtained in [19] and [30]. 3. Extension of the dual time step method proposed in this thesis work to a multiple time step case. As long as each consecutive integration step (starting from the largest) is an integer multiple of the next step, the extension is straight forward. There would be solution instants when only the fastest components are solved, instants when the fastest and the second fastest components are solved, and so on, until a general resynchronization instant, when the complete network is solved. This would guarantee stable results, according to the theory described in this thesis and would make the present code more robust. Another more interesting and challenging 103 approach would be to verify the possibility of latency exploitation when the integration steps do not follow the pattern above described of consecutive multiples. 4. Adaptation of the latency methodology to large networks where, according to the network partitioning chosen, different solution methods may be combined. Usually, for large power system networks, the transient effect of fast surges is restricted to the area relatively adjacent to where the transient event occurred. The solution for the rest of the network may be accurately obtained through a steady-state solution using phasors. Therefore, latency may seem to be a good alternative for performing phasorial and transient solutions in the same simulation. The subnetwork containing the buses most affected by the fast transients would be modelled for an EMTP-type of simulation, therefore resulting in a fast subnetwork, while the other buses would be solved in steady-state and grouped in a slow subnetwork. 5. Application of latency to the simulation of power electronic converters and controllers. Power electronics based systems are also good candidates for latency exploitation. Especially with the development of voltage and current source inverters using gate turn-off thyristors (GTO's) and insulated gate bipolar transistors (IGBT's), very high switching frequencies are obtained, and small integration steps are needed for an accurate simulation. In contrast, the power network may be subjected to much lower frequencies and a larger integration step may provide the accuracy desired. 104 B I B L I O G R A P H Y [I] H. W. Dommel, \"Digital Computer Solution of Electromagnetic Transients in Single and Multiphase Networks\", IEEE Transactions on Power Apparatus and Systems, vol.88, no.4, pp.388-399, April 1969. [2] H. W. Dommel, EMTP Theory Book, Microtran Power System Analysis Corporation, Vancouver, B C , 2 n d Edition, May 1992. Latest update: April 1996. [3] Allan Greenwood, \"Electrical Transients in Power Systems\", 2 n d Edition, Ed. lohn Wiley & Sons, New York, 1991. [4] J. R. Marti, L. R. Linares, \"Real-time EMTP-based Transients Simulation\", IEEE Transactions on Power Systems, vol.9, no.3, pp.1309-1317, August 1994. [5] J. R. Marti, \"Accurate Modelling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations\", IEEE Transactions on Power Apparatus and Systems, vol.101, no.l, pp.147-155, January 1982. [6] L. Marti, \"Simulation of Transients in Underground Cables with Frequency-Dependent Modal Transformation Matrices\", IEEE Transactions on Power Delivery, vol.3, no.3, pp.1099-1110, July 1988. [7] T. Noda, N . Nagaoka, A . Ametani, \"Phase Domain Modeling of Frequency-Dependent Transmission Lines by Means of an A R M A Model\", IEEE Transactions on Power Delivery, vol. 11, no. 1, pp.401 -411, January 1996. [8] F. Castellanos, J. R. Marti, \"Full Frequency-Dependent Phase-Domain Transmission Line Model\", IEEE Transactions on Power Systems, vol.12, no.3, pp.1331-1339, August 1997. [9] V . Brandwajn, H. W. Dommel, I. I. Dommel, \"Matrix Representations of Three-Phase n-Winding Transformers for Steady-State and Transient Studies\", IEEE Transactions on Power Apparatus and Systems, vol.101, no.6, pp.1369-1378, June 1982. [10] A . Morched, L. Marti, J. Ottevangers, \" A High Frequency Transformer Model for the EMTP\", IEEE Transactions on Power Delivery, vol.8, no.3, pp.1615-1626, July 1993. [II] S. Chimklai, J. R. Marti, \"Simplified Three-Phase Transformer Model for Electromagnetic Transient Studies\", IEEE Transactions on Power Delivery, vol.10, no.3, pp.1316-1325, July 1995. [12] A . M . Gole, R. W. Menzies, H. M . Turanli, D. A . Woodford, \"Improved Interfacing of Electrical Machines Models to Electromagnetic Transients Programs\", IEEE Transactions on Power Apparatus and Systems, vol.104, no.9, pp.2367-2373, September 1985. [13] N . J. Bacalao, P. de Arizon, R. O. Sanchez, \" A Model for the Synchronous Machine Using Frequency Response Measurements\", IEEE Transactions on Power Systems, vol.10, no.l , pp.457-464, February 1995. 105 [14] IEEE Working Group on Surge Arrester Modeling, \"Modeling of Metal Oxide Surge Arresters\", IEEE Transactions on Power Delivery, vol.7, no.l, pp.302-309, January 1992. [15] E. J. Tarasiewicz, A . S. Morched, A . Narang, E P. Dick, \"Frequency Dependent Eddy Current Models for Nonlinear Iron Cores\", IEEE Transactions on Power Systems, vol.8, no.2, pp.588-597, May 1993. [16] P. G. McLaren, R. Kuffel, R. Wierckx, J. Giesbrecht, L. Arendt, \" A Real Time Digital Simulator for Testing Relays\", IEEE Transactions on Power Delivery, vol.7, no.l, pp. 207-213,January 1992. [17] M . Kezunovic, M . Aganagic, V . Skendzic, J. Domaszewicz, J. K . Bladow, D. M . Hamai, S. M . McKenna, \"Transients Computation for Relay Testing in Real-Time\", IEEE Transactions on Power Delivery, vol.9, no.3, pp. 1298-1307, July 1994. [18] O. Devaux, L. Levacher, O. Huet, \"An Advanced and Powerful Real-Time Digital Transient Network Analyser\", IEEE Transactions on Power Delivery, vol.13, no. 2, pp.421-426, April 1998. [19] J. R. Marti, L. R. Linares, J. A . Hollman, F. A . Moreira, \"OVNI: Integrated Software/Hardware Solution for Real-Time Simulation of Large Power Systems\", in Proceedings 2002 PSCC (Power Systems Computation Conference), Session 33, Paper 4, Seville, Spain. [20] G. Kron, \"Tensorial Analysis of Integrated Transmission Systems, Part III - The Primitive Division\", AIEE Transactions, vol.71, Part 3, pp.814-821, 1952. [21] Q. Wang, J. R. Marti, \" A Waveform Relaxation Technique for Steady State Initialization of Circuits with Nonlinear Elements and Ideal Diodes\", IEEE Transactions on Power Delivery, vol.11, no.3, pp.1437-1443, July 1996. [22] J. R. Marti, L. R. Linares, J. Calvino, H . W. Dommel, J. Lin, \"OVNI: An Object Approach to Real-time Power System Simulators\", in Proceedings 1998 POWERCON (International Conference on Power System Technology), Beijing, China, pp.977-981. [23] P. F. Cox, R. G. Burch, P. Yang, D. E. Hocevar, \"New Implicit Integration Method for Efficient Latency Exploitation in Circuit Simulation\", IEEE Transactions on Computer-Aided Design, vol.8, no. 10, pp.1051-1064, October 1989. [24] P. F. Cox, R. G. Burch, D. E. Hocevar, P. Yang, B. D. Epler, \"Direct Circuit Simulation Algorithms for Parallel Processing\", IEEE Transactions on Computer-Aided Design, vol.10, no.6, pp.714-725, June 1991. [25] R. A . Saleh, A . R. Newton, \"The Exploitation of Latency and Multirate Behavior Using Nonlinear Relaxation for Circuit Simulation\", IEEE Transactions on Computer-Aided Design, vol.8, no.12, pp.1286-1298, December 1989. 106 [26] A . I. Zecevic, N . Gacic, \" A Partitioning Algorithm for the Parallel Solution of Differential-Algebraic Equations by Waveform Relaxation\", IEEE Transactions on Circuits and Systems-I: Fundamental Theory and Applications, vol.46, no.4, pp.421-434, April 1999. [27] E. Lelarasmee, A . E. Ruehli, A . L. Sangiovanni-Vincentelli, \"The Waveform Relaxation Method for Time-Domain Analysis of Large Integrated Circuits\", IEEE Transactions on Computer-Aided Design, vol.1, pp.131-145, July 1982. [28] M . L. Crow, M . Ilic, \"The Parallel Implementation of the Waveform Relaxation Method for Transient Stability Simulations\", IEEE Transactions on Power Systems, vol.5, no.3, pp.922-932, August 1990. [29] M . Ilic-Spong, M . L. Crow, M . A . Pai, \"Transient Stability Simulation by Waveform Relaxation Methods\", IEEE Transactions on Power Systems, vol.2, no.4, November 1987. [30] F. A . Moreira, J. A . Hollman, L. R. Linares, J. R. Marti, \"Network Decoupling by Latency Exploitation and Distributed Hardware Architecture\", in Proceedings 2001 IPST (International Conference on Power System Transients), Rio de Janeiro, Brazil, pp.317-321. [31] T. Kato, K. Ikeuchi, \"Variable Order and Variable Step-Size Integration Method for Transient Analysis Programs\", IEEE Transactions on Power Systems, vol.6, no.l , pp.206-213, February 1991. [32] A . Kurita, H . Okubo, K. Oki, S. Agematsu, D. B. Klapper, N . W. Miller, W. W. Price, J. J. Sanchez-Gasca, K. A . Wirgau, T. D. Younkins, \"Multiple Time-Scale Power System Dynamic Simulation\", IEEE Transactions on Power Systems, vol.8, no.l , pp.216-223, February 1993. [33] M . L. Crow, J. G. Chen, \"The Multirate Method for Simulation of Power System Dynamics\", IEEE Transactions on Power Systems, vol.9, no.3, pp. 1684-1690, August 1994. [34] J. R. Marti, 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. [35] K . A . Sakallah, S. W. Director, \"SAMSON2: An Event Driven VLSI Circuit Simulator\", IEEE Transactions on Computer-Aided Design\", vol.4, pp.668-685, October 1985. [36] M . L. Crow, J. G. Chen, \"The Multirate Simulation of FACTS Devices in Power System Dynamics\", IEEE Transactions on Power Systems, vol.11, no.l , pp.376-382, February 1996. [37] A . Semlyen, F. de Leon, \"Computation of Electromagnetic Transients Using Dual or Multiple Time Steps\", IEEE Transactions on Power Systems, vol.8, no.3, pp. 1274-1281, August 1993. [38] L. R. Linares, J. R. Marti, \"Sub-area Latency in a Real Time Power Network Simulator\", in Proceedings 1995 IPST (International Conference on Power System Transients), Lisbon, Portugal, pp.541-545. 107 [39] J. Reeve, R. Adapa, \" A New Approach to Dynamic Analysis of A C Networks Incorporating Detailed Modeling of DC Systems. Part I: Principles and Implementation\", IEEE Transactions on Power Delivery, vol.3, no.4, pp.2005-2011, October 1988. [40] R. Adapa, J. Reeve, \" A New Approach to Dynamic Analysis of A C Networks Incorporating Detailed Modeling of DC Systems. Part II: Application to Interaction of DC and Weak A C Systems\", IEEE Transactions on Power Delivery, vol.3, no.4, pp.2012-2019, October 1988. [41] M . Sultan, J. Reeve, R. Adapa, \"Combined Transient and Dynamic Analysis of H V D C and FACTS Systems\", IEEE Transactions on Power Delivery, vol.13, no.4, pp. 1271-1277, October 1998. [42] K . K. Fung, S. Y . R. Hui, \"Fast Simulation of Multistage Power Electronic Systems with Widely Separated Operating Frequencies\", IEEE Transactions on Power Electronics, vol.11, no.3, pp.405-412, May 1996. [43] C. Demi, P. Turkes, \"Fast Simulation Technique for Power Electronic Circuits with Widely Different Time Constants\", IEEE Transactions on Industry Applications, vol.35, no.3, pp.657-662, May/June 1999. [44] Prabha Kundur, \"Power System Stability and Control\", Ed. McGraw-Hill, New York, 1994. [45] M . Poller, \"Discrete Time State Space Analysis of Electrical Networks\", in Proceedings 2001 IPST (International Conference on Power System Transients), Rio de Janeiro, Brazil, pp. 305-310. [46] C. W. Ho, A . E. Ruehli, P. A . Brennan, \"The Modified Nodal Approach to Network Analysis\", IEEE Transactions on Circuits and Systems, vol.22, pp.504-509, June 1975. [47] Chi-Tsong Chen, \"Linear System Theory and Design\", Ed. Holt, Rinehart and Winston, New York, 1984. [48] T. C. Yu , J.R. Marti, \"z-Cable Model for Frequency Dependent Modelling of Cable Transmission Systems\", in Proceedings 2001 IPST (International Conference on Power System Transients), Rio de Janeiro, Brazil, pp.55-60. 108 APPENDIX A. THE MATE CONCEPT The M A T E (\"Multi-Area Thevenin Equivalent\") concept extends the main ideas of Diakoptics [20] in a conceptually clear and computationally efficient solution framework. For this purpose it embodies also the concepts of Multinode Thevenin Equivalents [2] and Modified Nodal Analysis [46]. The main idea in M A T E , as in Diakoptics, is that a large network can be split into smaller parts (subsystems) by tearing apart some of its branches. M A T E can be clearly formulated using partitioned matrices and, even though it applies to a general circuit solution, the case of an EMTP discrete-time solution is used here to explain its formulation [19]. Figure A . l shows a simple example network, consisting of two subsystems [A] and [B] connected by the link a i . Figure A . 1: M A T E example network If the link a in figure A . l did not exist, there would be two independent subnetworks, [A] and [B], with corresponding nodal equations: 109 (A.1) (A.2) Putting these equations in a \"common container\": 'A 0~ V 0 B VB_ (A.3) Including the link equation and expanding submatrices [A] and [B], eq. (A.3) can be written as A A A A 5, B2 B, A GA\\\\ GA\\2 0 GA14 0 0 0 0 VAX K\\ A2 GA2\\ GA22 GA2i 0 0 0 0 0 VA2 hA2 A 0 GA32 GA13 GA34 0 0 0 1 VA1 HA3 A GA4l 0 GA43 GA44 0 0 0 0 VA4 - kA4 Bi 0 0 0 0 GBU GB\\2 GB\\3 -1 VBl K B2 0 0 0 0 GB2\\ GB22 GB2i 0 VB2 hB2 0 0 0 0 G B l l GBi2 G B i i 0 VB3 h 0 0 1 0 -1 0 0 0 The diagonal elements of matrices [A] and [B] are built by adding all the conductances connected to the respective node, while the off-diagonal elements are the negative of the conductances connecting the two involved nodes. Rewriting eq. (A.4) using partitioned matrices: A B a A A 0 P VA K B 0 q VB hB a P' —z 0 (A.5) where p and q are given by 110 p = q = - l o o (A.6) Inverting matrix A and applying it to row A in eq. (A.5): A' VA A 0 P VB = A-1 ia VA 1 0 a VB = eA a- Ax p eA = A~ hA (A.7a) (A.7b) (A.7c) Inverting matrix B and applying it to row B in eq. (A.5): VA 0 B q VB B~l VA 0 1 b VB = eB b = £-> q eB = B~ hB (A.8a) (A.8b) (A.8c) And the last row in eq. (A.5) represents the link branch: VA p' q' —z VB = 0 (A.9a) 111 But [vA] and [vB] are determined in eqs. (A.7b) and (A.8b) respectively, in terms of [ia], [eA] and [eB]. Substituting in eq. (A.9a): e A - < p' q' —z 0 (A.9b) Recovering the original set of variables [v^ ] and [vB ]: VA 0 0 za VB = Za — p'a + q'b + z ea=P'eA+q'eB (A.9c) (A.9d) Using partitioned matrices once again to put eqs. (A.7b), (A.8b) and (A.9c) together: A B a A 1 0 a VA eA B 0 1 b VB eB a 0 0 za la (A. 10) From the last row in eq. (A. 10) the link current [ia] can be readily calculated: i =5-a Z„ (A.11) And the internal node voltages [v^ ] and [vB ] are updated by substituting the value of [ia ] in eq. (A. 10): 112 [ v J=[ e J-4 'J M = [eB]~b[ia] (A.12) (A. 13) The expansion of the M A T E concept for subnetworks connected by more than one link or for networks decoupled in more than two subnetworks is straight forward. If the subnetworks are connected through more than one link, there would be as many link currents as link connections. For example, i f the subnetworks A and B were connected through two links, [/„ ] and [ea ] would become: & a 2 The node voltages would still be updated in the same way as before, except that now [a] and [b] become matrices with the number of columns exactly the same as the number of links. In case M A T E decoupling generates, say, three subnetworks, there would be three decoupled block diagonal nodal equations and the number of branch equations would increase either to two or three depending i f the new subnetwork is connected to just one or both of the pre-existing subnetworks, respectively. 113 APPENDIX B. THE STATE VARIABLE DESCRIPTION The state variable description of a system is an alternative to the input-output description, which is applicable only when the system is initially relaxed. The input-output description for a linear, causal, time-invariant, and relaxed at time to system is given by where u is the input of the system; y is the output of the system; G is the impulse response matrix of the system. Equation (B.l) however, does not provide any information about the internal behaviour of the system. If the system is not initially relaxed, the output y will not only depend on the input u, but also on the initial conditions at t = to. The set of initial conditions required to determine the output y uniquely is called the state of the system at time t = to. In fact, the complete behaviour of the system is uniquely defined, including its state for all t>t0. The set of equations that describes the unique relations between the input, output, and state of a system is called a dynamical equation. If a system contains p inputs and q outputs, an n-dimensional linear and time-invariant dynamical equation is of the form: where x is the state of the system Matrices A , B, C, and D are respectively of \"nxn\", \"nxp\", \"qxn\", and \"qxp\" dimensions. Equation (B.2a) is called the state equation and eq. (B.2b) is called the output equation. This set of equations is called the state variable description of a system. From its analysis, general properties such as controllability, observability, and stability of the system; as well as the exact responses of the system due to some excitation may be determined. (B.l) x(t)-Ax(t)+Bu (t) y(t) = Cx(t)+Du (t) (B.2b) (B.2a) 114 In the case where the systems are electric networks, a systematic procedure for assigning state variables and writing the dynamical equations may be introduced. Most of the times it is not necessary to choose all the inductor currents and capacitor voltages as state variables. This is particularly true in the case of large networks where not all capacitor voltages or inductor currents may be independent. This procedure is extracted from [47]. It is well known that i f all the inductor currents and the capacitor voltages of an R L C network are known, then the behaviour of the network is uniquely determinable for any input. The procedure described here requires the definition of the concepts of tree, link, and cutset of a network. A tree of a network is defined as any connected graph (connection of branches) containing all the nodes of the network but not any loop. Every branch in a given tree is called a tree branch. Every branch not in the tree is called a link. A cutset of a connected network is any minimum set of branches such that the removal of all the branches in this set causes the remaining network to be unconnected. With respect to any fixed tree, every link and some tree branches form a unique cutset called a fundamental cutset. With these definitions in mind, the procedure for assigning state variables is as follows: 1) Choose a tree called a normal tree. The branches of the normal tree are chosen in the following priority order: voltage sources, capacitors, resistors, inductors, and current sources. Hence, a normal tree consists of all the voltage sources, the maximum number of permissible capacitors, the resistors, and finally the minimum number of inductors. Usually, it does not contain any current source. 2) Assign the charges or voltages of the capacitors in the normal tree and the flux or current of the inductors in the links as state variables. The voltages of the capacitors in the links and the currents of the inductors in the normal tree do not need to be chosen as state variables. 3) Express the branch variables (branch voltage and current) of all the resistors, the capacitors in the links, and the inductors in the normal tree in terms of the state variables and the inputs by applying the 1s t or 2 n d Kirchhoff laws to the fundamental loops or cutsets of these branches. 4) Apply the 1s t or 2 n d Kirchhoff laws to the fundamental loop or cutset of every branch that is assigned as a state variable. 115 APPENDIX C. THE Z-LINE MODEL Traditional EMTP multiphase transmission line and cable models use modal decomposition to decouple the physical system of conductors into a mathematically-equivalent one. This modal decomposition has the advantage that each natural decoupled mode has its own propagation velocity. In the case of constant parameter line models, the transformation matrix that relates phase and mode quantities are real and constant, and the modal-based models work very well. However, when taking into account the frequency dependence of the line parameters, the transformation matrix also becomes frequency dependent. Usually, for a symmetric system of overhead conductors, the transformation matrix can be assumed real and constant, and the results are still very reliable. In the case of underground cables and asymmetric overhead lines though, this assumption may result in very inaccurate results. To overcome this limitation, frequency dependent transformation matrices are used. These matrices are synthesized with rational functions in the frequency domain. This approach may become computationally expensive and is prone to numerical instability of the frequency dependent modal domain functions. In order to avoid the need for a frequency dependent transformation matrix, F. Castellanos and J. R. Marti proposed the z-line model [8]. This model can be formulated directly in phase coordinates and avoids the use of modal decomposition matrices. The theory of the z-line model is now reviewed as formulated in [8]. The modeling of transmission lines is based on the travelling wave equations, which can be expressed for a given frequency co as ~ = [ZY]V and ^ = [YZ]l (C.l) ax ax [ZY] and [YZ] are full matrices that couple the propagation of voltage and current waves in all phases. The elements of the frequency dependent series impedance matrix Z can be described as Z.. (CQ) = RV(CQ)+JCQ($ (a>)+E?) (C.2) Ry (co) includes the resistance of the conductor and ground return; 116 L'™ (at) is the internal inductance associated with the flux inside the conductor and ground return: Lfy1 is the external inductance due to the flux outside the conductor. This last term is frequency independent and matrix Z can be written as [Z] = [zhss] + jco[£f«] (C.3) ^ z l ™ ( ^ = Rv(co)+jcoLf(c^ The shunt admittance matrix [Y] can be written (assuming shunt conductance G = 0) as [7] = 7 « [ C ] (C.4) Matrices [if*1 ] and [C] only depend on the geometry of the system and matrix [YZ] can be written as [YZ] = ja>[C][z^ ] - co2 [C][IT ] ( C 5 ) Equation (C.5) shows that [YZ] can be expressed as the combination of two main components: the first one related to the losses and ground effect corrections and the second one related to the ideal propagation. These components may be separated as shown in figure C . l . 117 correction for losses and internal flux ideal line, only external flux ground Figure C l : Separation of basic effects in the z-line model The ideal line can be modelled directly in phase domain and the losses can be lumped; therefore avoiding the necessity of a transformation matrix. However, the losses can only be lumped i f the line is split into several z-line cells such as the one represented in figure C l . The number of \"cells\" required for an accurate simulation must increase as frequency increases. This is a limitation of the z-line model since in traditional EMTP simulation, the integration step cannot be larger than the minimum travelling time of the lines within the network. The reader interested in more details about the z-line model, including the fitting procedure for the synthesis of the frequency dependent lumped lossy section, is referred to [8] and to [48]. In [48] the concept of the z-line model has been extended for the modelling of underground cables as well. 118 "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2002-11"@en ; edm:isShownAt "10.14288/1.0065858"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Electrical and Computer Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Latency techniques in power system transients simulation"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/13663"@en .