Step by Step Eigenvalue Analysis with EMTP Discrete Time Solutions by JORGE ARIEL HOLLMAN Industrial Engineer, Universidad Nacional del Comahue, 1996 M.A.Sc., The University of British Columbia, 2000 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University Of British Columbia October, 2006 c© Jorge Ariel Hollman 2006 Abstract The present work introduces a methodology to obtain a discrete time state space representation of an electrical network using the nodal [G] matrix of the Electromagnetic Transients Program (EMTP) solution. This is the first time the connection between the EMTP nodal analysis solution and a cor- responding state-space formulation is presented. Compared to conventional state space solutions, the nodal EMTP solution is computationally much more efficient. Compared to the phasor solutions used in transient stability analysis, the proposed approach captures a much wider range of eigenvalues and system operating states. A fundamental advantage of extracting the sys- tem eigenvalues directly from the EMTP solution is the ability of the EMTP to follow the characteristics of nonlinearities. The system’s trajectory can be accurately traced and the calculated eigenvalues and eigenvectors correctly represent the system’s instantaneous dynamics. In addition, the algorithm can be used as a tool to identify network partitioning subsystems suitable for real-time hybrid power system simulator environments, including the im- plementation of multi-time scale solutions. The proposed technique can be implemented as an extension to any EMTP- based simulator. Within our UBC research group, it is aimed at extending the capabilities of our real-time PC-cluster Object Virtual Network Integra- tor (OVNI) simulator. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Research Background . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Off-Line Power System Simulation . . . . . . . . . . . 5 1.2.2 Real-Time Power System Simulation . . . . . . . . . . 6 1.2.3 Large Scale Power System Stability Analysis . . . . . . 9 1.3 Research Claim and Contributions . . . . . . . . . . . . . . . 11 1.4 List of Publications . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 The Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Traditional Analysis Tools . . . . . . . . . . . . . . . . 17 2.2.2 Advantages of Step by Step Trajectory Analysis . . . . 27 iii Table of Contents 3 Eigenvalue Analysis with EMTP Solution . . . . . . . . . . . 29 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 EMTP Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Nodal Analysis . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 Numerical Integration Accuracy . . . . . . . . . . . . . 34 3.2.3 Discretization Rule Stability . . . . . . . . . . . . . . . 37 3.3 Discrete Time State Space Formulation of an Electrical Net- work from the EMTP Nodal Matrix Formulation . . . . . . . 40 3.4 EMTP Solution in State Space form . . . . . . . . . . . . . . 42 3.4.1 Discrete Time State Space Equation for an Inductance 44 3.4.2 Discrete Time State Space Equation of a Capacitance . 46 3.4.3 Discrete Time State Space Equation of a Resistance . . 48 3.4.4 Treatment of Series Connection RL and RC . . . . . . 49 3.4.5 Treatment of Series Branches . . . . . . . . . . . . . . 50 3.4.6 Treatment of Parallel Branches . . . . . . . . . . . . . 55 3.5 Nonlinear Networks . . . . . . . . . . . . . . . . . . . . . . . . 57 3.6 Continuous Time Eigenvalues from Discrete Time Domain So- lution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.6.1 Discrete to Continuous Time Mapping . . . . . . . . . 60 3.6.2 Time Discretization Considerations . . . . . . . . . . . 63 4 Application to the OVNI simulator . . . . . . . . . . . . . . . 66 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 OVNI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.1 The MATE Concept . . . . . . . . . . . . . . . . . . . 68 4.2.2 Discrete State Space Formulation with MATE . . . . . 73 4.3 The Latency Concept . . . . . . . . . . . . . . . . . . . . . . . 76 4.3.1 Application of the New Discrete Time State Space For- mulation to Latency . . . . . . . . . . . . . . . . . . . 77 4.4 Automatic Time Step Selection . . . . . . . . . . . . . . . . . 79 5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.1 Comparison of State Space Solutions between Continuous Time and Discrete Time Domains . . . . . . . . . . . . . . . . . . . 81 5.1.1 RLC Series Circuit State Space Equation from Contin- uous Time Domain . . . . . . . . . . . . . . . . . . . . 82 5.1.2 RLC Series Circuit State Space Equation from the [G] EMTP Matrix . . . . . . . . . . . . . . . . . . . . . . 85 iv Table of Contents 5.2 Latency Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3 Two-Area Resonant Circuit . . . . . . . . . . . . . . . . . . . 99 5.3.1 Case I - Poorly Damped . . . . . . . . . . . . . . . . . 100 5.3.2 Case II - Improved Damping . . . . . . . . . . . . . . . 106 5.4 Two-Area Case with Nonlinear Load . . . . . . . . . . . . . . 110 5.5 Voltage Collapse in a Radial System . . . . . . . . . . . . . . 112 5.6 11 Bus Two Area System . . . . . . . . . . . . . . . . . . . . . 122 6 Conclusions and Recommendations . . . . . . . . . . . . . . . 128 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.2 Recommendations for Future Work . . . . . . . . . . . . . . . 131 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 A The QR Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 v List of Tables 2.1 Relevant modes for the two-area system case . . . . . . . . . . 25 vi List of Figures 1.1 Relation among research projects at UBC power system group 8 2.1 Linear time invariant (LTI) dynamic system . . . . . . . . . . 18 2.2 Two-area 11 bus test system . . . . . . . . . . . . . . . . . . . 22 2.3 Phasor value of bus 8 voltage for the two-area test system, small-signal instability observed with detailed time domain simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Phasor value of bus 8 voltage for the two-area test system, small-signal simulated with FTD (no oscillation are captured by this solution) . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1 Inductor element . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Inductor discrete time equivalent with trapezoidal rule . . . . 32 3.3 Magnitude frequency response of integration rules . . . . . . . 35 3.4 Phase frequency response of integration rules . . . . . . . . . . 36 3.5 State space representation of a multiple input/output system . 39 3.6 State space system diagram block . . . . . . . . . . . . . . . 43 3.7 Discrete time equivalent of an Inductor . . . . . . . . . . . . . 45 3.8 Discrete equivalent of a Capacitor . . . . . . . . . . . . . . . . 47 3.9 Discrete equivalent of a Resistor . . . . . . . . . . . . . . . . . 49 3.10 Discrete equivalent of a Series RL . . . . . . . . . . . . . . . . 50 3.11 Generic circuit . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.12 Aggregate treatment of parallel RLC branches . . . . . . . . 56 3.13 Piecewise linear inductance with two slopes . . . . . . . . . . . 58 3.14 Implementation of two slope piecewise linear inductance . . . 59 3.15 Continuous time and discrete time stability mapping for trape- zoidal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.16 Continuous time and discrete time stability mapping for back- ward Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.1 Generic description of a MATE link . . . . . . . . . . . . . . . 68 vii List of Figures 4.2 Generic circuit for discrete state space formulation with MATE 69 4.3 General update scheme for discrete state space formulation with MATE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4 Hybrid real-time/soft real-time simulator layout . . . . . . . . 75 4.5 Automatic EMTP time step selection scheme . . . . . . . . . . 80 5.1 RLC series circuit . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Discretized RLC series circuit . . . . . . . . . . . . . . . . . . 86 5.3 Latency test circuit . . . . . . . . . . . . . . . . . . . . . . . . 93 5.4 Circuit with slow and fast subareas . . . . . . . . . . . . . . . 99 5.5 Two-area resonant test case . . . . . . . . . . . . . . . . . . . 99 5.6 Node voltages for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.7 Discrete time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . . . . . . . . . . . . . . 104 5.8 Reconstructed continuous time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . . 105 5.9 Node voltages for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . . . . . . . . . . . . . . . . . . 107 5.10 Discrete time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . . . . . . . . . . . . . . 108 5.11 Reconstructed continuous time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . . 109 5.12 Two-area circuit test case with nonlinear load . . . . . . . . . 110 5.13 Nonlinear load - discrete time eigenvalues . . . . . . . . . . . . 111 5.14 Nonlinear load - continuous time reconstructed eigenvalues from discrete time solution . . . . . . . . . . . . . . . . . . . . 112 5.15 Two-area system with nonlinear load . . . . . . . . . . . . . . 113 5.16 A simple radial system . . . . . . . . . . . . . . . . . . . . . . 114 5.17 Voltage collapse of a simple radial system . . . . . . . . . . . . 116 5.18 Reconstructed continuous time eigenvalue trajectory over time for a simple radial system with voltage collapse . . . . . . . . 117 5.19 Complex plane projection of continuous time eigenvalues over time line of simple radial system of Fig. 5.16 with voltage collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.20 Discrete time eigenvalue trajectory from EMTP solution for the simple radial system of Fig. 5.16 with voltage collapse . . 119 5.21 Zoom in of discrete time trajectory of Fig. 5.20 . . . . . . . . 120 viii List of Figures 5.22 Discrete time eigenvalue 1 magnitude, |eig1|, over time . . . . 121 5.23 Two-area 11 bus test system . . . . . . . . . . . . . . . . . . . 123 5.24 Two-area 11 bus test system - Bus Voltages . . . . . . . . . . 124 5.25 Two-area 11 bus test system - PV curve at Bus 9 . . . . . . . 125 5.26 Two-area 11 bus test system - Discrete time Eigenvalues nor- mal operation point A . . . . . . . . . . . . . . . . . . . . . . 126 5.27 Two-area 11 bus test system - Discrete time Eigenvalues at collapse point B . . . . . . . . . . . . . . . . . . . . . . . . . . 127 ix Acronyms ATA Advanced Technology Attachment, a disk drive implementation that inte- grates the controller on the disk drive itself DSP Digital Signal Processing EDF Electricité de France EMTDC Transients Simulation Software by Manitoba HVDC Research Centre EMTP Electromagnetic Transients Program EPRI Electric Power Research Institute FACTS Flexible AC Transmission Sys- tems FTD Fast Time Domain GW Gigawatt, one billion times the power unit of watt Gbps Gigabits per second, a data transfer speed measurement for high-speed networks HVDC High Voltage Direct Current IDE Intelligent Drive Electronics, an interface for mass storage device IREQ Hydro Quebec Research Institute LTI Linear Time Invariant MAM Modified Arnoldi Method MATE Multi Area Thevenin Equivalent OVNI Object Virtual Network Integrator RTDNS Real-Time Distributed Network Simulator PCI Peripheral Component Interconnect, a local bus standard developed by Intel Corporation x Acronyms RTNS Real-Time Network Solver SATA Often abbreviated Serial ATA or S-ATA, an evolution of the Parallel ATA physical storage interface TNA Transient Network Analyzer TSAT Transient Stability Assessment Tool UBC The University of British Columbia UILO UBC’s University - Industry Liaison Office IPO Initial Public Offering ATP Alternative Transient Program CDA Critical Damping Adjustment JIIRP Joint Infrastructure Interdependencies Research Project I2Sim Infrastructures Interdependencies Simulator xi Notation [G] admittance matrix [GB] admittance matrix of subnetwork [B] GB43 element (4, 3) of admittance matrix [B] gB43 admittance of branch nodes 4 and 3 in subnetwork [B] [L] incidence matrix [A]d discrete time transition matrix [A]c continuous time transition matrix λ continuous time eigenvalue z discrete time eigenvalue ẋ(t) future state in continuous time domain x(t+∆t) future state in discrete time domain a′ compact notation for a(t−∆t) a compact notation for a(t) vkm branch voltage of node k to node m h(t)km history of branch km at (t) h(t−∆t)km history of branch km at (t−∆t) _ h(t)i compact notation for history of branch i at (t) _ h(t−∆t)i history of branch i at (t−∆t) _ h i history of branch i at (t) _ h′i compact notation for history of branch i at (t−∆t) [k∧] diagonal matrix of branch nature coefficients [g∧] diagonal matrix of branch discretization rule coefficients ζ damping ratio ζc damping ratio from continuous time domain f oscillation frequency f cn undamped natural frequency from continuous time domain [P ] Participation matrix [φ] Right eigenvector [ϕ] Left eigenvector fNy Nyquist frequency xii Acknowledgments Writing this thesis was a long process of hard work and dedication during which many people have somehow left their imprint on my life. This ac- knowledgement tries to reflect my gratitude to some of them. This thesis could not have been written without the encouragement and support of my wife, Sandra Zappa-Hollman. Her infinite and unconditional love, endless patience and support made easier even the most challenging times. Our daughters, Roćıo and Serena, are also part of this family effort of surviving in spite of having both parents writing doctoral dissertations under the same roof. Thank you Roćıo and Serena for reminding me with your cheerfulness and love of what is really important in life. I cannot find words to fully capture how thankful I am to my parents Maŕıa del Carmen and Esteban. Their support, teachings and care helped me to reach this point in life. I couldn’t have done it without them. Thank you dad for taking me during school holidays with you to the power substa- tions and Hydro plants, which sparked my passion for power systems. Thank you mom for getting me my first 386-based PC in which I ran power system simulations using the EMTP written by Dr. Dommel and using the Mart́ı line model. I never imagined then that I would be so blessed to do both my M.A.Sc. and Ph.D. under their supervision. A special recognition goes to my sister Verónica, and my parents and brother in law Lieschen, Hansi and Andrés, who always believed in us and uncondi- tionally supported us all these years. I also want to extend my thanks to our friends Verónica, Mat́ıas, Lucas, Dante and Sof́ıa Salibián, who became our family in Canada and made our daily life a pleasant adventure; and to our good friends Paula Murcia, Pablo Ruiz, Jesús Calviño-Fraga, Verónica, Javier and Mat́ıas Gazzarri, Lito and Alicia Caro, Jorma and Laura Neuvonen, Hazem and Samar Sabbagh, Gill and Jorge Palejko, Maŕıa Isabel and Felipe Gajardo, Naoko and Jérémie Séror, Aya and Tomoaki Nakashima, and Daniel and Lise Lindenmeyer. Numerous people at UBC have provided me with support in different, xiii Acknowledgments invaluable ways: Dr. Jesús Calviño-Fraga, Dr. Luis Linares, Dr. Carlos Ven- tura, Dr. Juri Jatskevich, Mazana Armstrong, Tom de Rybel and Alejandro Cervantes-Larios, by being good friends, colleagues and sharing stimulating conversations; Dr. Hermann Dommel, who on many occasions spent time providing helpful comments and suggestions; Dr. Michael Davies and Dr. Vijay Bhargava, who have given me their unconditional support, and last but not least, everybody in the UBC Power group. I especially want to thank Doris Metcalf, Cathleen Holtvogt, Anne Coates and David Chu Chong for their continued professional assistance and friend- ship. This dissertation would not have been possible without the generous fi- nancial support provided by several sources, including The National Sciences and Engineering Research Council of Canada through an ISP-2 Research Grant; Dr. José Mart́ı with appointments as Research Assistant; the Depart- ment of Electrical and Computer Engineering via part-time appointments as Teaching Assistant; Fundación Repsol-YPF through the Estenssoro Award; The University of Bristish Columbia through the Kenneth George Vansacker Memorial Prize; Consejo Nacional de Investigaciones Cient́ıficas y Técnicas; and Dr. Prabha Kundur and Dr. Ali Moshref from Powertech labs with part-time research and consulting contracts. To all of them goes my appre- ciation. Finally, I would especially like to express my most sincere gratitude to my mentor and supervisor, Dr. José Mart́ı. His constant availability, pa- tience, guidance, leadership and support, have encouraged me throughout these years. Furthermore, he has become a role model beyond the academic world. Jorge Ariel Hollman xiv To Serena, Roćıo and Sandra xv Chapter 1 Introduction 1.1 Motivation Since Thomas Alva Edison set up the world’s first commercial electric power system in lower Manhattan in September 1882, power systems grids have undergone a continuous and aggressive expansion which enhanced their op- eration complexity. Just taking into consideration the case of the United States-Canadian system, for instance, this amazing expansion of the gener- ation capacity grew from 171 GW in 1960 up to 943 GW by 20041. During the late eighties and nineties, first in the United Kingdom, subsequently in countries such as Argentina, Chile, New Zealand, many European countries, Brazil, and finally in some states of the United States the power industry was challenged to a new operational paradigm: the implementation of a dereg- ulated environment. This new operational paradigm does not accept the traditional extra investment in network infrastructure to ensure conserva- tive secure operation margins. The traditional technical constraints together with the newly introduced financial constraints have forced the utilities to 1Information available from the Energy Information Administration web site, http://www.eia.doe.gov 1 Chapter 1. Introduction operate the systems near their stability limits. The maximization of revenue can only be achieved under the new deregulated scheme by not only knowing accurately but also knowing moment by moment the actual safe operational limits of the system. These two factors have created new interest in the analysis of stability problems. The high degree of interconnection and the requirement of synchronicity among all the system components makes the operation of power systems not an easy task. Furthermore, the unbundling of generation, transmission and supply is less oriented towards the physical nature of the synchronously interconnected power systems. In addition, the latest trend for the new century power system layout is the fast development of distributed generation. It is easy to visualize that centralized control operation of the entire system will be impossible due to the size and complexity of the objective function and the response time constraints. The United States-Canadian electric- ity system alone deals with more than 320,000 km of transmission lines at 230 KV or higher level, 950 GW of generating capability distributed among nearly 3,500 utility organizations. The August 2003 North American East coast blackout which affected an area of 50 millon people with a loss of 62 GW of load lends strong support to the validity of this statement. A centralized coordination paradigm will inevitably shift to a locally coor- dinated operation. The 21st century power system will make extensive use of distributed intelligent local processing nodes capable of making their own decisions and able to interact with other related nodes and with centralized 2 Chapter 1. Introduction regional controls. Since all parts of a power network interact with each other, local control action may require extensive information regarding the external system to which the node is connected. This information can be acquired and transmitted to all the distributed nodes or, even more efficiently, it can be simulated in real-time at each decision node. However, to implement sufficiently detailed external system representations at local nodes requires considerable computational power. This computational power must be ca- pable of simulating external system behaviour in real-time at the local node. Consequently, the concept becomes practical only if this computational so- lution can be provided in a fast, compact, accurate and inexpensive manner. Clear statements calling for further research on the proposed topic can be found in the literature. According to S. Gissinger [1], “Power system opera- tion requires increasingly complex decision-making in order to find the right compromise between economy and security”. The importance of considering technical issues as a variable of the economic equation is pointed out by M. Pereira [2]: “Today’s new environment brings new uncertainties and risks (...) Much talk about risk today is focused on financial risks but it is re- markable that many important risks are not measured in dollars”. As well, the relevance of an extensive and detailed representation of the power sys- tem for the stability analysis with modal analysis is addressed by P. Kundur in [3] “A mode of oscillation in one part of the system can excite units in a remote part due to mode coupling. Analysis requires detailed and same level of representation throughout the system”. The feasibility of distributed 3 Chapter 1. Introduction computation in power system stability studies is stressed by J.L. Jardim in [4]: “Parallel processing has been increasingly gaining acceptance as a way of speeding up the solution of highly intensive computational problems. It is particularly indicated for applications that exhibit a natural decoupling among computation tasks”. In sum, the above citations reflect but just a limited portion of the innumerable calls for further exploration of very large power systems stability analysis. In particular, a discrete time state space description of electrical networks from nodal equations is of fundamental interest for this research work. Some of these key research questions were inspired by Dr. Dommel in the Theory Book [5], where he pondered “could such a closed-form solution be used in an EMTP?... where the method becomes almost unmanageable, or useless, is in networks with nonlinear elements”. Hence, it is clear from the bibliographical research and from the above in- troduced scenario that there are great opportunities and an imminent need for the development of fast on-line power systems simulators applications seeking faster and more accurate analysis, especially those with integration of stability analysis of power systems. This new brand of simulators with the capability of identifying the natural areas of association of the power system components will allow for the implementation of a decentralized op- eration layout. Given the present development of The University of British Columbia’s (UBC) real-time power system simulator -Object Virtual Net- work Integrator- (OVNI) in the hardware and software, it is attractive to 4 Chapter 1. Introduction explore, develop and implement an extension of its capabilities to the stabil- ity analysis of power systems. 1.2 Research Background 1.2.1 Off-Line Power System Simulation The Electromagnetic Transients Program (EMTP) [6] has been the standard tool for digital off-line simulation of power system transients for more than thirty years. The EMTP became an attractive alternative to the bulky and expensive Transient Network Analyzers (TNA) due to its portability and re- duced cost. Much effort was invested during the last three decades to develop accurate and sophisticated models to represent the behaviour of the power system components. While off-line power system transients simulators such as the EMTP[6], EMTDC [7] and Power System Blockset [8], successfully cope with the requirements of system planning and transient studies areas, the off-line simulation environment exhibits great limitations when sensitive equipment is required to be tested in a closed loop with the modelled power system or when real-time decision is a must. This called for the development of a new breed of power systems simulators which could cope with these lim- itations. This new group is based on the concepts presented by H. Dommel in [6], [5], but it is capable of achieving real-time performance. Alternatively, some research groups developed hybrid solutions based on Digital Signal Processing (DSP) [9] to cope with the requirements of real-time equipment 5 Chapter 1. Introduction testing. Notwithstanding its limitations, the off-line simulation environment keeps improving and extending its capabilities of analysis of power systems by incorporating new areas of knowledge, such as artificial intelligence, neu- ral networks, stochastic analysis, risk theory, stability analysis, sensitivity analysis, and so forth. 1.2.2 Real-Time Power System Simulation After almost a decade now, various research groups have been able to achieve real-time performance [9], [10], [11], [12]. The real-time research group of the University of British Columbia has been continuously advancing the concept of a full-size real-time power system simulator, which eventually resulted in the OVNI [13] power system simulator. The simulator’s architecture is based on a cluster of PC workstations. Using partitioning techniques in the hardware implementation and network decoupling techniques in the software solution has led to the development of a distributed solver architecture, where a network of inexpensive workstations is used to achieve real-time simulation for fast transients of large power systems. There are mainly five well-known research groups working in the field of real-time power system simulators using different approaches. Electricité de France (EDF) and Hydro Que- bec Research Institute (IREQ) together with TEQSIM and Mitsubishi have chosen the very expensive and not so portable supercomputer architecture, while Manitoba Research Group has developed a hybrid solution based on a convenient arrangement of DSPs. On the other hand, Mitsubishi started 6 Chapter 1. Introduction with supercomputers and, after joint work with UBC, moved their research to PC’s. UBC’s research group was the first to choose an inexpensive PC solution and achieve real-time performance with this scheme. Recently [14] EDF presented results of a parallel approach based on shared memory architecture. This approach is still based on supercomputers and shows no extra gain in speed in proportion to the number of CPU’s in service because the time needed for communication between the CPU’s increases with the degree of parallelization. Mitsubishi research group achieved a PC cluster system based on six interconnected machines using a Myricom Myrinet giga-bit network. This approach presents a poor performance for hard real-time simulation due to its round trip overhead [12], but it is promis- ing for soft real-time applications. In contrast, UBC’s real-time PC cluster architecture shows a constant communication time for each type of config- uration link between the subsystem nodes, independently of the number of node solvers included in the array. UBC’s Real-time Network Solver (RTNS) [15] created by J. Mart́ı and L. Linares in 1993 was the very first PC based real-time network solver. Stemming from this work several other research projects originated at UBC, examples of which are the RTNS extensions and hardware implementation for testing relays [16] developed by Jesús Calviño; the Multi-layered tearing solution software [10]; and the Real-Time Distributed Network Simulator, developed by J. Hollman (RTDNS) [17]. See Fig.1.1. These efforts are part of the OVNI project [18], [13], [19], [20], [21] to develop a very large size 7 Chapter 1. Introduction EMTP RTNS (One-layer-tearing software) SRTNS for Relay Tester (Single-PC simulator) MRTNS (Multi-layer-tearing software) RTDNS (PC-Cluster Simulator) OVNI Hybrid Simulation Environment Hard or Soft Real-Time Off-Line Environment Transient and Stability Analysis Transient Analysis Figure 1.1: Relation among research projects at UBC power system group real-time power system simulator. By breaking the power network solu- tion into subparts while maintaining simultaneous solution of all the parts, the Real-Time Power Systems group at UBC has developed a PC-Cluster simulator made up of an extensible number of off-the-shelf inexpensive PC boards that is capable of simulating in real-time fast system phenomena for systems of practically any size just by extending the number of PC boards 8 Chapter 1. Introduction interconnected. RTDNS achieves one of the lowest latency2 network solu- tions available today. However, with the advent of the Multi Area Thevenin Equivalent (MATE) [21] solution scheme, a new, more flexible hardware in- terconnect solution became necessary. Also, the current system is limited in bandwidth and tied into the Intelligent Drive Electronics (IDE) hard-disk interface. This technology is soon to disappear in favour of serial Advanced Technology Attachment (SATA). The new hardware system, called OVNI- NET [22], circumvents any of the normal PC standard ports and is built directly on top of the industry-standard Peripheral Component Interconnect (PCI) bus. Also, the parallel connections between the nodes are replaced by high-speed 1.2 Gbps full-duplex serial links, allowing for a system of larger physical size with a more manageable wiring solution. These enhancements yield an even lower latency while significantly boosting the throughput. They also allow the hardware to be used on non-PC platforms such as Sun, Apple, HP, and IBM workstations. 1.2.3 Large Scale Power System Stability Analysis Although sophisticated power systems stability analysis techniques are avail- able [23], [24], [25], [26], [27], [28], [29], a large number of simplifications is a common characteristic among them in order to reduce the computa- tional burden, especially for the case of on-line stability analysis tools. As a 2Latency in this context is the time interval between the instant at which an instruction control unit issues a call for data and the instant at which the transfer of data is started. 9 Chapter 1. Introduction result, this situation makes it very difficult to determine whether observed system instability is due to the system properties or to the assumptions in the algorithms. Generally, recursive and extensive off-line simulations are required before the power system stability can be accurately assessed. The lack of detailed representation of the system in fast on-line assessment also becomes a constraint at the moment of identifying corrective responses, since not all the variables of the system are considered at the same instant. On the other hand, while detailed system representation is achieved in off-line analysis programs, the inherent solution times do not permit to achieve real- time performance. Most of the currently available transient analysis software works in the off-line environment. Nonetheless, some groups [30], [31], mo- tivated by the present industry transformation, have been able to reach soft real-time performance in the range of tens of seconds up to hundreds of min- utes, yet still based on a large number of previous off-line simulations in the time domain. Some stability tools use an intermediate approach between the phasorial and time domains. The so-called fast time domain simulation approach, based only on time domain simulation, focuses on the evolution of the system operating conditions driven by the subsequent slow dynam- ics while neglecting part of the fast dynamics present. In some cases the neglected dynamics, such as those associated one with high-voltage direct current (HVDC) and Flexible AC Transmission Systems (FACTS) devices, become fundamental factors of the possible corrective actions. Whereas this type of tools is satisfactory for certain kinds of studies, such as mid-term 10 Chapter 1. Introduction voltage stability, it is not capable of producing a correct stability assessment of the analyzed power system in the case of small-signal instability. 1.3 Research Claim and Contributions In light of the current need for fast and accurate analysis solution of power systems, this thesis makes the following contribution: The description and implementation of a new and original power system stability assessment methodology that identifies the system’s eigenvalue trajectories in a real-time EMTP solution incorporating the effect of switching phenomena and nonlin- ear load behaviour. Some of the advantages of the new methodology can be summarized as fol- lows: • Trajectory tracking of nonlinear element eigenvalues moment by mo- ment. • Capability of identifying suitable network partitioning schemes for ap- plication of multi-step integration solution in a hybrid power system simulator environment. • Capability of identifying small-signal interactions between elements by analyzing damping ratios of natural frequencies and participation in- dexes from discrete time eigenvalues. 11 Chapter 1. Introduction • Visualization of eigenvalues trajectories in discrete time for the purpose of assessing power system’s dynamic behaviour. • Automatic selection of discretization step from discrete time eigenvalue information. 1.4 List of Publications The following are publications in which the author of this Thesis is the prin- cipal author or co-author and which are in part related to the topic of this work. • T. De Rybel, J. Hollman, J. Mart́ı. OVNI-NET: A Flexible Cluster In- terconnect for the New OVNI Real-Time Simulator. Proceedings of the 15th Power System Computation Conference. Liége, Belgium, 2005. • J. Hollman, J. Mart́ı. Real-Time Network Simulation with PC-Clusters. IEEE Transactions on Power Systems,Vol. 18, No. 2, pp. 563-569, IEEE Transaction on Power Systems, 2003. • J. Mart́ı, L. Linares, J. Hollman, F. Moreira. OVNI: Integrated Soft- ware/Hardware Solution for Real-Time Simulation of Large Power Sys- tems. Proceedings of the 14th Power Systems Computation Conference. Seville, Spain, 7 journal pages, 2002. • F. Moreira, J. Hollman, L. Linares, J. Mart́ı. Network Decoupling by Latency Exploitation and Distributed Hardware Architecture. Proceed- 12 Chapter 1. Introduction ings of the International Conference on Power System Transients. Rio de Janeiro, Brazil, 2001. • J. Mart́ı, J. Hollman, J. Calviño. Implementation of a Real-Time Dis- tributed Network Simulator with PC-Cluster. Proceedings of the Inter- national Conference on Parallel Computing in Electrical Engineering, Quebec, Canada, 2000. • J. Hollman. Real Time Distributed Network Simulation with PC-Clusters. M.A.Sc. Thesis. Vancouver, 2000. 13 Chapter 2 Framework 2.1 Introduction This chapter provides a summary of current representative state of the art analysis tools in the power system stability environment, pointing out their strengths and limitations and the motivation for the development of new methodologies based on the EMTP solution. We highlight the capabilities of the EMTP simulators to follow nonlinear systems dynamic trajectories including their instantaneous eigenvalue/eigenvector states. 2.2 The Problem Maintaining interconnected power systems within their operational limits has been always challenging due to the highly dynamic and nonlinear nature of its components, such as loads, generators and changing operation topologies. Historically, due to the complexity of the problem, the assessment of oper- ation limits has been conducted in an off-line environment. Due to these computational limitations, the operational limits are computed well in ad- 14 Chapter 2. Framework vance of the time in which they are expected to occur. The off-line results are loaded into look-up tables that are accessed by the operators. During actual real-time operation, the system operators find the closest match be- tween the actual system condition and the study set loaded into the tables. The limits for that condition are then used as the actual operation limits of the system. The new operational paradigm introduced by de-regulation has moved system operation much closer to its limits and off-line conservative and structured operation conditions are no longer considered adequate. The trend is now towards on-line assessment tools. The theoretical frame of traditional stability analysis presents an interest- ing paradox: it has a highly detailed representation of elements such as generator and load controllers, while at the same time it uses a very low detail representation of other system operating aspects and components, like switching operations and nonlinear loads. Differential equations are used for the generators, while steady state phasor solutions are used for the electrical network, including the large component of motor loads. In many situations, this traditional approach is insufficient. One example is the case of interarea oscillation, which is a particular case of small-signal stability. Small-signal stability is the ability of the power system to maintain synchronism under small disturbances. Such disturbances occur continuously in the system be- cause of small variations in loads and generation. Small-signal instability can be of two forms: steady increase in generator rotor angle, or rotor oscillations of increasing amplitude. With present power system configurations, small- 15 Chapter 2. Framework signal stability is largely a problem of insufficient damping of oscillations and it can be classified into the following modes: Local, Interarea, Control, and Torsional [3]. Interarea modes are associated with the swinging of machines in one part of the system against machines in other parts. These electromechanically coupled elements can be geographically very far apart. The problem occurs when the electromechanically closely coupled machines are interconnected by weak ties. Interarea oscillations can severely constrain power transfers across large distances and can lead to widespread system disturbances if cascading outages of transmission lines occur due to oscillatory power swings, like those during the blackout in western North America on August 10, 1996. If the representation of the system is not detailed enough, those dominant modes can be overlooked. Here, dominance refers to largeness in modulus. For the case of interarea oscillations, the size of the computational problem and the aggregation degree of the reduced system are the main constraints of the traditional stability analysis approach. Reduced-order techniques of analysis are based on eigenvalue analysis of only the dominant modes. This approach provides excellent information of associated frequencies and damp- ing ratios of the dominant coupled modes; however, they are usually inten- sive and have to be used off-line. Modified Arnoldi [32], and spectral analysis based on Prony’s algorithm [33], [34], are the preferred traditional techniques for interarea stability analysis of large power systems. 16 Chapter 2. Framework 2.2.1 Traditional Analysis Tools Methods such as Prony’s spectral analysis and Fast Time-Domain (FTD) simulations are well developed. These tools reconstruct the system dynamics by using curve fitting techniques on quasi-steady state simulation results. For the calculation of eigenvalues the most used methods are the Arnoldi method and the Modified Arnoldi Method (MAM). An overview of these representative methods of power system analysis follows. Prony Analysis Prony analysis was developed by Gaspard Riche, Baron de Prony [35] in 1795 to explain the expansion of various gases. Riche proposed fitting a sum of exponentials to equally space data points, and extended the model to interpolate at intermediate points. Prony analysis is both a signal analy- sis technique and a system identification method widely used in biomedical monitoring, geophysical analysis, speech processing and power system elec- tromechanical oscillation analysis. As compared to Fourier analysis, Prony analysis has the advantage of providing estimation of damping coefficients in addition to frequency, phase and amplitude of the Fourier technique. While the Fourier series fits a sum of undamped complex exponentials to the data, Prony analysis fits a sum of damped complex exponentials [36]. A brief de- scription of Prony analysis is provided next. Let us consider the case of a Linear Time Invariant (LTI) dynamic system 17 Chapter 2. Framework such as the one illustrated in Fig. 2.1 (LTI) x(t) u(t) y(t) Figure 2.1: Linear time invariant (LTI) dynamic system Here, y(t) is the system’s output x(t) is the state of the system u(t) is the input to the system The evolution of the LTI system is expressed by ẋ(t) = Ax(t) +Bu(t) (2.1) where A and B are constant matrices. If the input is removed (u(t) = 0) for all t, equation (2.1) can be rewritten as ẋ(t) = Ax(t) (2.2) where A is a square matrix of size n x n whose eigenvalues are λi, left eigenvectors are qi and right eigenvectors are pi. The solution to equation 18 Chapter 2. Framework (2.2) can be expressed as a sum of n components, ẋ(t) = n∑ i=1 (qTi x0)pie (λit) (2.3) With u(t) = 0 for all t, our LTI system output equation is given by y(t) = Cx(t) (2.4) Prony analysis directly estimates the parameters of the eigenstructure de- scribed in equation (2.3) by fitting a sum of complex damped sinusoids to evenly spaced time sample values of the output, that is, y(t) = L∑ i=1 Aie (σit) cos(2pifit+ φi) (2.5) Where, Ai is the amplitude of the component i σi is the damping coefficient of component i fi is the frequency of component i L is the total number of damped exponential components y(t) is the estimate of observed data for y(t) consisting of N samples y(t) = y[k] with k = 0, 1, 2, 3, ..., (N − 1) evenly spaced over the time 19 Chapter 2. Framework Using Euler’s theorem, the complex sinusoids can be represented as a sum of exponentials, cos(2pifit+ φi) = e(j2pifit+φi) + e−(j2pifit+φi) 2 = e(j2pifit)ejφi 2 + e−(j2pifit)e−jφi 2 (2.6) replacing (2.6) into (2.5) and letting t = kT the samples of y(t) are given by y[k] = L∑ i=1 Ci µ k i︸︷︷︸ poles (2.7) Where, Ci = Ai 2 ejφi (2.8) µi = e (σi+j2pifi)T (2.9) Prony analysis computes Ci and µi in three steps [37], as follows: i. solve a linear prediction model, constructed from the observed data set ii. find roots of the characteristic polynomial formed from the linear pre- diction coefficients iii. solve the original set of linear equations to yield the estimates of the exponential amplitude and sinusoidal phase With Ci and µi known, the amplitude, frequency and phase damping coeffi- cients are computed from equations (2.8) and (2.9). 20 Chapter 2. Framework The main limitation of identification techniques, such as Prony analysis, for real-time applications is the prohibitively large computational time required for the fitting process. The Prony algorithm involves the solution of an overdetermined set of linear equations and roots of a high order polynomial, which are numerically intensive operations. In addition, identification tech- niques are sensitive to noise. It is known that Prony analysis behaves poorly when a signal is embedded in noise, leading to a significant misestimation of damping and frequency terms. Thus, Prony analysis requires care in the choice of processing parameters [38]. Also Prony analysis requires that an oscillation should exist in the studied time interval. If not, the calculated re- sult contains an equivalent of a higher order eigenterm. This is critical for the case of fast and highly nonlinear elements, in which the eigenvalue location is changing over time. Notice that even in the case when the exact set of eigen- values is found, the Prony method is not capable of providing information about which of the network components are responsible for the particular poorly damped modes. In contrast, the methodology presented in this work shows a direct association between network component and eigenmode. Fast Time Domain simulation and Transient Security Assessment Tool Fast Time Domain (FTD) simulation [39] simplifies the model of a power system in an quasi-steady state simulation approach in which many of the faster device dynamics are ignored and the simulation is focused on the ac- 21 Chapter 2. Framework tion of slower devices such as ULTC’s, overexcitation limiters and relays. The FTD simulation technique was developed for fast dynamic voltage sta- bility assessment. Even though the method considers the slow devices of the system, it is not always suitable to identify small-signal instability dynamics due to the simplifications imposed in order to speed up the solution. As an example, the illustrated two-area system in Fig. 2.2 is simulated using the industry standard Transient Stability Assessment Tool (TSAT) [28] for the case of a three phase fault applied at bus 6, a fault clearing time of 100 ms, and a simulation time step of 10 ms. The trapezoidal integration rule is used in this simulation. The simulation is first performed using detailed time do- main analysis, see Fig. 2.3, and later using the FTD technique, see Fig. 2.4. The results show the 60Hz phasor value of the voltage at bus 8 for the same contingency in both cases. The full time domain simulation shows that the 1 2 65 7 8 9 10 11 3 4 G1 G2 G4 G3 L7 L9 25 km 25 km10 km 10 km 110 km 110 km Area 1 Area 2 400 MW Figure 2.2: Two-area 11 bus test system system becomes unstable for the given conditions while the FTD technique simulation shows change in the phasor value of the voltage and therefore no 22 Chapter 2. Framework B us v ol ta ge m ag ni tu de (p u) Time in seconds 0.00 3.87 7.74 11.61 15.48 19.350.00 0.28 0.56 0.84 1.12 1.40 Figure 2.3: Phasor value of bus 8 voltage for the two-area test system, small- signal instability observed with detailed time domain simulation oscillation after the contingency is applied and cleared. The FTD technique is not capable of capturing this case of small-signal instability. In very large systems, the low frequency modes are the only ones capable of coupling through very large distances to create instabilities. The higher frequency modes are naturally more damped and they can not propagate far away from the source; consequently, they do not contribute to interarea oscillations. For instance, a complete eigenvalue analysis of the two areas of Fig. 2.2, using the Prony technique (Table 2.1), makes clear the presence of local and interarea modes. The term “local” is associated with the swinging of units at a particular generating station with respect to the rest of the system. The term “interarea” refers to the swinging of many units in one 23 Chapter 2. Framework B us v ol ta ge m ag ni tu de (p u) Time in seconds 0.00 6.00 12.00 18.00 24.00 30.000.970 0.984 0.998 1.012 1.026 1.040 Figure 2.4: Phasor value of bus 8 voltage for the two-area test system, small- signal simulated with FTD (no oscillation are captured by this solution) part of the system against many units in another part of the system. Interarea oscillations are caused by two or more groups of closely coupled units being interconnected by weak (high impedance) ties. Using Prony’s method on the voltage signals at the system generation nodes (1,2,3,4) all the existing modes are fitted. As observed in Table 2.1, the local modes (except for the 1.03 Hz mode at generators 3 and 4) have sufficient damping, but the interarea mode has insufficient damping, thus leading to the small-signal instability problem observed in the simulation. The location of the observed modes suggests the corrective action, in this case increase the damping of the 0.62Hz interarea mode at generators # 1, # 2, # 3, # 4. It is important to emphasize that system identification techniques such as Prony analysis are 24 Chapter 2. Framework Table 2.1: Relevant modes for the two-area system case Freq [Hz] Damp % Real [1/s] Imag [rad] Type of observed mode Node 1 0.62 -2.144 0.0831 3.875 Interarea Gen.# 1 0.97 6.702 -0.4103 6.109 Local between G1 & G2 1.2 -0.579 0.0459 7.942 Local between G1 & G2 Node 2 0.62 -1.077 0.0421 3.910 Interarea Gen.# 2 0.97 8.926 -0.5432 6.061 Local between G1 & G2 1.2 1.052 0.0421 6.061 Local between G1 & G2 Node 3 0.62 -2.357 0.0925 3.924 Interarea Gen.# 3 1.03 -1.884 0.1215 6.446 Local between G3 & G4 1.7 13.677 -0.1087 10.743 Local between G3 & G4 Node 4 0.62 -1.344 0.0527 3.921 Interarea Gen.# 4 1.03 3.442 -0.2335 6.820 Local between G3 & G4 1.7 -0.675 0.0734 10.870 Local between G3 & G4 not able to provide a direct association between network components and modes with poor damping. Arnoldi Method The Arnoldi method (1951) [32] was proposed as a means of reducing a ma- trix to the Hessenberg form. Arnoldi already suggested that the process could also provide good approximations to some eigenvalues if stopped before its completion. Later this concept was fully developed by Saad (1980) [40] mak- ing the Arnoldi method evolve into one of the most successful ones for non- Hermitian eigenproblems. This is in contrast to the so-called direct methods which have to be complete in order to provide any useful results. The Arnoldi method has poor numerical stability due to extensive numerical cancellation, 25 Chapter 2. Framework as pointed out by Wilkinson [41], and requires re-orthogonalization in or- der to improve its numerical stability performance. However, the inclusion of re-orthogonalization makes the method computationally less efficient [41]. The so-called Modified Arnoldi method (MAM) by EPRI [42] is probably the most popular method for identifying the set of dominant eigenvalues of power systems. The Arnoldi algorithm can be summarized as follows: Given a square matrix A of order n, if n steps of Arnoldi’s method is carried out then an orthogonal reduction to Hessenberg form is achieved, AV = V H (2.10) WhereH is the upper Hessenberg matrix of order n with positive subdiagonal elements and V is an orthogonal matrix. Matrices V and H are uniquely determined by the first column of V , a unit-norm vector that is called the initial vector. When an initial vector fails afterm steps, the following relation holds: AVm − VmHm = fe∗m (2.11) Where the vector f is called the residual of the m-step Arnoldi factorization. The input required by the Arnoldi algorithm are the matrix A, the number of steps m and the initial vector v1 of norm 1. The method provides as output the orthogonal matrices Vm, Hm, the residual f and β, so that verifies the equation (2.11) for β = ‖f‖2. In practical situations, the number of steps m required to achieve a good 26 Chapter 2. Framework approximation may be too large. Consequently, the method becomes expen- sive due to storage and computational load that grows in every step. This is critical for the case of real-time applications. One possible strategy to overcome this problem is the concept of restarting the algorithm after m steps trying not only to cope with the growth of computational load, but also benefitting from the previously computed information to define the new initial vector. The ARPACK library provides a Fortran implementation of the implicit restarted Arnoldi method, for both real and complex arithmetic. 2.2.2 Advantages of Step by Step Trajectory Analysis A simple but illustrative general non-power system scenario can be used to show the advantage of discrete time trajectory analysis compared to steady state analysis. Take the case of a marathon runner who has a race tar- get time of 3h 08min 00 s. Based on a steady state performance, in order to achieve the desired time, the runner would have to be able to run the 42.195km with a pace of 4min 27 s/km. For the sake of producing a more accurate estimation of the required race pace, a fade of 10%/km based on the history of the runner will be considered and thus the runner would have to cope with an initial pace of 4min02 s/km. This initial pace will gradually increase 10%/km until the finish line. However, on the race day, the runner will most probably find out that because of the highly nonlinear dynamics of the race, such as the elevation of the course, the temperature, the wind, and the physical condition of the runner during the race day, she or he will have a 27 Chapter 2. Framework different pace performance, despite of long training hours. If the performance is analyzed km by km and adjusted based on the discrete considered points, it will be possible to compute a better estimation of the overall performance. Real-time embedded measurements of pace (such as a GPS device) will give enough information to the runner to adjust her or his pace to achieve the target time. By analogy, the same is valid for more complex and large systems such as power grids. Solutions based on steady state or quasi-steady state computa- tions will be less accurate to represent highly nonlinear characteristics. This additional accuracy is the main advantage of having analysis tools that are capable of tracing point by point the transitions of nonlinear elements in the network. The EMTP is a proven tool to provide this type of information [6], [5], [43], [44], [45], [46] and many steady state based power system stud- ies are verified against EMTP transient simulations. Providing the EMTP with the ability to trace eigenvalue trajectories extends the capabilities of this highly accurate tool and the ability to represent highly sensitive system operating conditions. In addition, the methodology developed in this Thesis directly maps the eigenvalues to those system components responsible for their existence, therefore, allowing for the most effective corrective actions. 28 Chapter 3 Eigenvalue Analysis with EMTP Solution 3.1 Introduction This chapter starts with an overview of the EMTP solution, which aims to provide the reader with sufficient background to reach a full understanding of the proposed new methodology. It is followed by a description of the proposed methodology to extract eigenvalues from the EMTP solution. The research contributions included in this chapter are: 1. Eigenvalue extraction from EMTP solution. 2. Eigenvalue trajectories visualization in discrete time for the purpose of assessing dynamic behaviour. 3.2 EMTP Solution The EMTP solution, based on H.W. Dommel’s algorithm [6], is probably the most widely accepted computer program for power systems transient studies 29 Chapter 3. Eigenvalue Analysis with EMTP Solution [47]. In addition it is used as an educational tool in power system university programs [48], [49], [50]. The beauty and efficiency of the EMTP solution are rooted in the nature of the discretization process. As clearly described in [19], the EMTP solution first defines the discretization at the element level instead of discretizing the full set of differential equations that describe the network behaviour. Con- sequently, in the EMTP every component relationship between its current and voltage is modelled with a discrete time equivalent. These difference equations are represented by discrete time equivalents circuits consisting of resistance and source combinations. The values depend on the integration rule applied. For instance, for the case of an Inductor (Fig. 3.1) the relation between the L Figure 3.1: Inductor element voltage across it and its current is given by VL(t) = L diL(t) dt (3.1) 30 Chapter 3. Eigenvalue Analysis with EMTP Solution Solving for the voltage, VLdt = LdiL(t) (3.2)∫ t t−∆t VLdt = L ∫ t t−∆t diL(t) (3.3)∫ t t−∆t VLdt = L[iL(t)− iL(t−∆t)] (3.4) where ∆t is the discretization time step. Approximating the left hand side with the trapezoidal formula we obtain VL(t) + VL(t−∆t) 2 ∆t = L[iL(t)− iL(t−∆t)] (3.5) Then the voltage across the inductor is VL(t) = 2L ∆t︸︷︷︸ Requiv iL(t) + [−VL(t−∆t)− 2L ∆t iL(t−∆t)︸ ︷︷ ︸ history ] (3.6) Equation 3.6 can be represented by the equivalent discrete time circuit of Fig. 3.2. The trapezoidal integration rule is used in the EMTP with good accuracy, as illustrated in Appendix A. This rule, however, introduces numerical oscil- lations at discontinuities. This situation is averted in the EMTP by using the Critical Damping Adjustment algorithm (CDA) developed by Mart́ı and Lin [51]. As described by Linares in [19], in systems with a large component of power electronic devices the backward Euler integration rule can be an 31 Chapter 3. Eigenvalue Analysis with EMTP Solution eh(t) v L + - i L 2L/ t Figure 3.2: Inductor discrete time equivalent with trapezoidal rule effective alternative for Real-time simulation. If a constant discretization time step is applied, the nodal analysis matrix consists of a constant matrix of conductances and a current source vector. If a variable discretization time step is used, the conductance matrix has to be recalculated whenever the time-step is changed. The conductance matrix also needs to be reevaluated if there is a change in the network topology. This is the case, for instance, of circuits with nonlinear characteristics or switching electronic devices. 3.2.1 Nodal Analysis Nodal analysis is based on Kirchoff’s current law, which states that the al- gebraic sum of currents entering a node is zero, leading to the formulation of a system with the form [Y ][V ] = [J ] (3.7) 32 Chapter 3. Eigenvalue Analysis with EMTP Solution where, Y is the admittance matrix of order N x N V is the vector of N node voltages referred to ground J is the vector of N current sources injected into the nodes N is the number of nodes, excluding ground For each independent voltage source, the voltage of the node at which the source is connected is known. If there are N2 sources, the number of nodes with unknown voltages is N1 = N −N2 (3.8) which can be partitioned as follows, Y11 Y12 Y21 Y22 V1 V2 = J1 J2 (3.9) then Y11V1 = J1 − Y12V2 (3.10) where, V1 is the vector of N1 nodes with unknown voltages referred to ground V2 is the vector of N2 nodes with voltages sources connected between the node and ground (known voltages) Y11, Y12, Y21, Y22 are admittances submatrices of order N1 x N1, N1 x N2, N2 x N1, N2 x N2, respectively J1 is the vector of current sources entering the N1 unknown voltages nodes 33 Chapter 3. Eigenvalue Analysis with EMTP Solution J2 is the vector of current sources entering the N2 known voltages nodes Equation 3.10 constitutes a system of N1 equations that can be solved for V1 by any numerical method. Currents in the elements are then obtained using the original elements’s branch equations. The admittance matrix Y is defined as follows: • Diagonal elements of the Y matrix come from the sum of all admit- tances at the corresponding node. • Off-diagonal elements of the Y matrix come from the negative of the admittances connecting the corresponding nodes. 3.2.2 Numerical Integration Accuracy When analyzing the total distortion error introduced by solving a system with a discrete time solution, both the highest frequency represented and the distortion error of the integration rule have to be taken into consideration. The Nyquist frequency is the bandwidth of a discrete time signal, and it is equal to half of the sampling frequency of that signal. In EMTP discrete time simulation, the sampled signal will contain frequencies ranging from 0Hz to the Nyquist frequency fNy, fNy = fsample 2 = 1 2∆t (3.11) The smaller the integration time step, the higher Nyquist frequency. There 34 Chapter 3. Eigenvalue Analysis with EMTP Solution 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Trapezoidal Calviño 3 Backward & Forward Euler f(pu) = f sim(Hz)/fNy(Hz) H e/ H, M ag ni tu de d ist or tio n (p. u.) Figure 3.3: Magnitude frequency response of integration rules is no aliasing in the EMTP solution because the input to the system, the sources, have values generated at each ∆t of the system solution (no higher frequency values are sampled in). However, as the frequencies in the net- work solution approach the Nyquist frequency, the accuracy of the solution decreases (“distortion error”). The distortion error introduced by the cho- sen integration rule increases as the frequency gets closer to the Nyquist fre- quency. Figures 3.3 and 3.4 show the magnitude and phase distortion [He/H], [ΦHe−Φexact] for four discretization rules. The magnitude distortion is shown 35 Chapter 3. Eigenvalue Analysis with EMTP Solution in per unit of the ideal response (He = rule’s response, H = exact response) while the phase distortion is shown in degrees. The frequency in those curves is shown in per unit of the Nyquist frequency, i.e. f(pu) = f(Hz) / fNy(Hz). For a given frequency in the simulation, the error will decrease if the time 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −100 −80 −60 −40 −20 0 20 40 60 80 100 Trapezoidal Calviño 3 Backward Euler Forward Euler f(pu) = f sim(Hz)/fNy(Hz) Φ H e − Φ e xa ct Ph as e di st or tio n [d eg ree s] Figure 3.4: Phase frequency response of integration rules step is reduced (i.e., if the Nyquist frequency is “pushed away”). In power system transients simulation it is desirable to obtain an accurate represen- tation of high frequencies when fast transients, fast switching operations, HVDC converters, and TVR studies occur. In real-time simulation when the 36 Chapter 3. Eigenvalue Analysis with EMTP Solution computer is synchronized to external hardware requiring a solution point ev- ery ∆t time units, there is a limit in the system size that can be solved such that the computer processing time plus the overhead to put each sample out into the external hardware is less than that ∆t. A solution to this problem is to implement the simulation with a distributed multi-pc’s architecture, such as the ones proposed and implemented in [17],[52],[22]. 3.2.3 Discretization Rule Stability The stability of the integration rule can be assessed through the Transfer Function in the Z-domain [53]. For the particular case of an Inductance discretized with trapezoidal rule, i(t)− i(t−∆t) = ∆t 2L + ∆t 2L v(t−∆t) (3.12) Applying the Z-transform, I(z)− z−1I(z) = ∆t 2L V (z) + ∆t 2L z−1V (z) (3.13) Ye(z)Trap = I(z) V (z) = ∆t 2L (z + 1) (z − 1) (3.14) Trapezoidal with a p1 = 1 and z1 = −1 has the characteristic of a stable integrator and a stable differentiator, with bounded oscillations at disconti- nuities for the latter case. 37 Chapter 3. Eigenvalue Analysis with EMTP Solution Similarly, for an inductor discretized with backward Euler the Z-transfer function is given by Ye(z)BE = ∆t L z (z − 1) (3.15) Backward Euler with a p1 = 1 and a z1 = 0 has the behavior of a stable integrator and a stable differentiator. Similarly, for the case of Calviño 3 [54], which has a stable integrator and differentiator characteristic, we get the transfer function, Ye(z)C3 = ∆t k0L 1 [1− k1z−1 + k2z−2 − k3z−3] (3.16) Figures 3.3 and 3.4 show the distortion introduced by the discretization rule as we approach the Nyquist frequency can be obtained as follows [53]. If we make z = ejω∆t in (3.14) we get the frequency response of the discretized inductor with trapezoidal rule, Ye(w) = ∆t 2L (ejω∆t + 1) (ejω∆t − 1) (3.17) Similarly, for backward Euler in (3.15) we obtain Ye(w)BE = ∆t L ejω∆t (ejω∆t − 1) (3.18) 38 Chapter 3. Eigenvalue Analysis with EMTP Solution and for Calviño 3 in (3.16) we obtain Ye(w)C3 = ∆t k0L 1 [1− k1e−jω∆t + k2e−jω2∆t − k3e−jω3∆t] (3.19) Notice from these figures that trapezoidal does not introduce phase distortion whereas backward Euler introduces a strong phase distortion. Trapezoidal can have stability problems at discontinuities, this is overcome by applying the concept of critical damping adjustment (CDA) [51]. Backward Euler has no problems at discontinuities. The behaviour of trapezoidal, backward Euler and other rules at discontinuities is analyzed in [51], where CDA combining the best characteristics of trapezoidal and backward Euler rules is introduced. 3.5, as follows. internal states u 1 (t) u 2 (t) u 2 (t) u m (t) y 1 (t) y 2 (t) y 2 (t) y p (t) x 1 (t),x 2 (t),x 2 (t),...x n (t) Figure 3.5: State space representation of a multiple input/output system 39 Chapter 3. Eigenvalue Analysis with EMTP Solution 3.3 Discrete Time State Space Formulation of an Electrical Network from the EMTP Nodal Matrix Formulation The concept of representing a discrete time system by a set of first order difference equations became a standard tool for the research engineer by the end of the 1950s. These techniques have since become generally known as state space representations. Any Linear Time Invariant (LTI) system can be described by a set of discrete time state space equations. A state space representation utilizes three types of variables to model the dynamics of a system: the input, which are the external forces driving the system; the output, which are the quantities directly measurable by an external observer; and the states, which characterize the internal dynamics of the system. A set of state variables can be computed as a function of the present and past states and the present and past system inputs. The definition implies that if we know the state at time t we can then compute the energy stored in the system at that instant. In a general framework, the state variables can be chosen as an arbitrary combination of inner system variables. This generalization builds some distance between the state and its physical interpretation, but it makes more evident that the choice of state variables is not unique. A state space representation of a system with multiple inputs and outputs can be visualized in a block diagram, see Fig. Here we assume m inputs, p outputs 40 Chapter 3. Eigenvalue Analysis with EMTP Solution and n states. A system governed by n linear difference equations of order one requires n internal states for a minimal state space representation. Using forward Euler discretization, the values of the state variables at discrete time point (t) can be expressed as a linear combination of values of the state and input variables at the previous discrete time point (t−∆t), x1(t) ... xn(t) = a11 . . . a1n ... . . . ... an1 . . . ann x1(t−∆t) ... xn(t−∆t) + b11 . . . b1m ... . . . ... an1 . . . anm u1(t−∆t) ... um(t−∆t) (3.20) while the output variables at discrete time (t) can also be expressed as linear combinations of the values of the input and state variables at discrete time (t). y1(t) ... yp(t) = c11 . . . c1n ... . . . ... cp1 . . . cpn x1(t) ... xn(t) + d11 . . . d1m ... . . . ... dp1 . . . dpm u1(t) ... um(t) (3.21) 41 Chapter 3. Eigenvalue Analysis with EMTP Solution Even though, as discussed in the previous chapter, the forward Euler rule is not stable for EMTP simulation, we will show subsequently that when adapt- ing the EMTP solution to state space form, for the chosen EMTP variables, forms (3.20) and (3.21) are obtained although the network components are discretized using stable rules such as trapezoidal and backward Euler. For a single input-output system (see figure 3.6), (3.20) and (3.21) became x(t) = A x(t−∆t) +B u(t−∆t) (3.22) y(t) = C x(t) +D u(t) (3.23) System states change from moment to moment, and what defines the cur- rent value of any state of the system are the memory, current input and the connectivity topology. Changes of the states can be originated by the input or by the system’s dynamics. Changes due to input are forced, while changes due to system dynamics reflect the nature of the system. 3.4 EMTP Solution in State Space form H. Dommel in [6] proposed the following nodal solution form for the Electro- magnetic Transients Program (EMTP), [ G ] [ v(t) ] = [ h(t) ] + [ u(t) ] (3.24) 42 Chapter 3. Eigenvalue Analysis with EMTP Solution u(t) D B A CZ-1 y(t)x(t+ t) x(t) Figure 3.6: State space system diagram block Where [v(t)] are the node voltages, [u(t)] are the external current sources, [G] is the matrix of nodal conductances and [h(t)] is the vector of history sources. The system “output” is given by [ v(t) ] = [ G ]−1 [ h(t) ] + [ G ]−1 [ u(t) ] (3.25) The updating formula for the history sources can be expressed as [ h(t+∆t) ] = [ A ] [ h(t) ] + [ B ] [ u(t) ] (3.26) The formula (3.26) establishes the connection between the history source concept of the EMTP solution and the previous state concept of the state space formulation. Vector [h(t)] in (3.24) constitutes the “system memory” (history sources). This memory is a mixture of system self dynamics h(t−∆t) and the input modulated dynamics u(t − ∆t) and it fits perfectly the state 43 Chapter 3. Eigenvalue Analysis with EMTP Solution space definition. The EMTP solution updates [h(t)] in terms of [h(t −∆t)] and external current sources (3.26), which follows exactly the form (3.20) of the state space formulation. To summarize, the method proposed here defines the branch history terms as states of the electrical network, but it retrieves all the necessary information from the EMTP nodal solution and the case data input. If there is no input applied to the system, the changes in [h(t)] constitute the homogeneous solution of the system. By excluding all external sources we can compute the transition matrix of the system due to the nature of the system. Finally, from [A] we can study the eigenvalues of the system and analyze its stability behavior as will be shown in detail subsequently. 3.4.1 Discrete Time State Space Equation for an Inductance Consider a single inductance discretized using the Trapezoidal integration rule, (Fig. 3.7). This basic circuit element is used to model single-phase shunt reactors, discharge circuits, smoothing reactors in HVDC converter stations, the inductive part of Thevenin/Helmholtz equivalent circuits, the inductive part of single-phase nominal pi circuits of balanced systems operation and even the inductive equivalent part for loads of low frequencies. The voltage across the element is vkm(t) = vk(t)− vm(t) (3.27) 44 Chapter 3. Eigenvalue Analysis with EMTP Solution i g = t / 2L k m h(t) v km k m L Figure 3.7: Discrete time equivalent of an Inductor and the update formula, hkm(t) = −gvkm(t−∆t)− ikm(t−∆t) (3.28) where h is the memory of the discretized inductance and can be chosen as a state variable. The current through the inductance is ikm(t) = gvkm(t)− hkm(t) (3.29) while the current ikm(t−∆t) is ikm(t−∆t) = gvkm(t−∆t)− hkm(t−∆t) (3.30) 45 Chapter 3. Eigenvalue Analysis with EMTP Solution By replacing (3.30) into (3.28), we obtain the state space equation of a self inductance, hkm(t) = −gvkm(t−∆t)− gvkm(t−∆t) + hkm(t−∆t) (3.31) Rearranging, hkm(t)︸ ︷︷ ︸ new state = hkm(t−∆t)︸ ︷︷ ︸ old state +(−2g)vkm(t−∆t)︸ ︷︷ ︸ forcing function (3.32) Compare this formula with (3.22) in the state space description. As men- tioned earlier, even though in state space (3.22) corresponds to forward Euler discretization, in the EMTP (3.32) was obtained from trapezoidal discretiza- tion. The one-element state equation for the inductance can be extended to branches connected in a network. The branch histories and branch voltages can be easily described by nodal equations using the network connectivity in- formation contained in the incidence matrix. In this case the forcing function v will be defined by the network solution. 3.4.2 Discrete Time State Space Equation of a Capacitance A single capacitance discretized using the trapezoidal integration rule is shown in Fig. 3.8. This model is typically used to represent series and shunt capacitors, shunt capacitances in nominal pi-circuit representations of 46 Chapter 3. Eigenvalue Analysis with EMTP Solution i g = 2C / t k m h(t) v km k m C Figure 3.8: Discrete equivalent of a Capacitor transmission lines, equipment in HVDC converter stations such as surge ca- pacitors and filters, capacitance potential transformers, capacitive voltage dividers, and stray capacitances of transformers, generators, especially for the case of transient recovery voltage and lightning studies, among others. Again, we start with the update formula from the nodal equation hkm(t) = gvkm(t−∆t) + ikm(t−∆t) (3.33) The current through the capacitance is ikm(t) = gvkm(t)− hkm(t) (3.34) while the current for (t−∆t), analogously to the case of the inductance is, ikm(t−∆t) = gvkm(t−∆t)− hkm(t−∆t) (3.35) 47 Chapter 3. Eigenvalue Analysis with EMTP Solution By replacing (3.35) into (3.33), the state space description of a single capac- itance is obtained, hkm(t) = gvkm(t−∆t) + gvkm(t−∆t)− hkm(t−∆t) (3.36) Rearranging, hkm(t)︸ ︷︷ ︸ new state = −hkm(t−∆t)︸ ︷︷ ︸ old state +2g vkm(t−∆t)︸ ︷︷ ︸ forcing function (3.37) Which has the same form as (3.22) for the state space description. Analo- gously to the inductance case, the one-element state equation for the capac- itance can be extended to branches connected in a network, in which case the forcing function v will be defined by the network solution. 3.4.3 Discrete Time State Space Equation of a Resistance In this case, Fig. 3.9, hkm(t) = 0 for all t. Even though resistive elements do not contribute with a state because they do not have branch history, they do provide damping. To avoid a special case in the state space formulation the R element can be discretized together with either an inductance or a capacitance branch as described in the following section. 48 Chapter 3. Eigenvalue Analysis with EMTP Solution g = 1/R k m h(t) = 0 k m R Figure 3.9: Discrete equivalent of a Resistor 3.4.4 Treatment of Series Connection RL and RC The series connection of a resistor with an inductor (RL) or a resistor with a capacitor (RC) can be represented with a combined equivalent discrete el- ement. The combined element approach not only takes care of the damping introduced by the resistive element, but also improves the performance of the overall solution, since it reduces the number of nodes to solve and, con- sequently, the number of branches. For the case of an inductor in series with a resistor, the aggregation process is shown in Fig. 3.10. The general form of equation (3.32) is kept, hkm(t) = hkm(t−∆t)− 2gRLvkm(t−∆t) (3.38) where, gRL = 1 R + 2L ∆t (3.39) 49 Chapter 3. Eigenvalue Analysis with EMTP Solution eh(t) v L + - i L 2L/ t g RL h(t) v RL i L LR R RL discrete equivalent k m k m eh(t) v RL + -k m R RL = R + 2L/ t Figure 3.10: Discrete equivalent of a Series RL Analogously for the series RC, hkm(t) = −hkm(t−∆t) + 2gRCvkm(t−∆t) (3.40) gRC = 1 R + ∆t 2C (3.41) 3.4.5 Treatment of Series Branches Initially we will consider pure L or pure C branches, see Fig. 3.11, to illustrate the concept. 50 Chapter 3. Eigenvalue Analysis with EMTP Solution h L1 h C2 h L3 21 g C2 g L1 g L3 L3 C2 L1 1 2 Branch 3 B ra n c h 1 B ra n c h 2 Figure 3.11: Generic circuit The nodal equations for the generic circuit of Fig. 3.11 are [G][v] = [j] (3.42) G11 G12 G21 G22 v1 v2 = j1 j2 In detail:gL1 + gL3 −gL3 −gL3 gC2 + gL3 v1(t) v2(t) = hL1(t) + hL3(t) hC2(t)− hL3(t) v1(t) v2(t) = [G−1] j1(t) j2(t) 51 Chapter 3. Eigenvalue Analysis with EMTP Solution v1(t) v2(t) = R11 R12 R21 R22 j1(t) j2(t) (3.43) In the circuit of Fig. 3.11 we want to identify the transition matrix [A]d by establishing a relationship in terms of the nodal solution. In what follows, to simplify the notation, the voltage of branch i is denoted as _ v i while the history of the branch i is denoted as _ hi. We know from equations (3.32) and (3.37) that for each branch element L or C, the branch histories are: _ h1(t) _ h2(t) _ h3(t) = k∧1 0 0 0 k∧2 0 0 0 k∧3 _ h1(t−∆t) _ h2(t−∆t) _ h3(t−∆t) + g∧1 0 0 0 g∧2 0 0 0 g∧3 _ v 1(t−∆t) _ v 2(t−∆t) _ v 3(t−∆t) (3.44) In particular for our circuit, k∧1 = 1 g∧1 = −∆tL1 k∧2 = −1 g∧2 = 4C2∆t k∧3 = 1 g∧3 = −∆tL3 (3.45) Notice that k∧ = 0 and g∧ = 0 for the case of a purely resistive branch. 52 Chapter 3. Eigenvalue Analysis with EMTP Solution The relationship between branch and node voltages is given by the branch- node incidence matrix: _ v 1(t−∆t) _ v 2(t−∆t) _ v 3(t−∆t) = 1 0 0 1 1 −1 v1(t−∆t) v2(t−∆t) (3.46) Substituting (3.43) in (3.46), _ v 1(t−∆t) _ v 2(t−∆t) _ v 3(t−∆t) = 1 0 0 1 1 −1 R11 R12 R21 R22 j1(t−∆t) j2(t−∆t) (3.47) The nodal history sources in terms of the branch histories are given by j1(t−∆t) j2(t−∆t) = _h1(t−∆t) + _h3(t−∆t)_ h2(t−∆t)− _ h3(t−∆t) = 1 0 1 0 1 −1 _ h1(t−∆t) _ h2(t−∆t) _ h3(t−∆t) (3.48) 53 Chapter 3. Eigenvalue Analysis with EMTP Solution substituting (3.47), (3.48) in (3.44), _ h1(t) _ h2(t) _ h3(t) = k∧1 0 0 0 k∧2 0 0 0 k∧3 ︸ ︷︷ ︸ [k∧] _ h1(t−∆t) _ h2(t−∆t) _ h3(t−∆t) + g∧1 0 0 0 g∧2 0 0 0 g∧3 ︸ ︷︷ ︸ [g∧] 1 0 0 1 1 −1 ︸ ︷︷ ︸ [L] R11 R12 R21 R22 ︸ ︷︷ ︸ [G−1] 1 0 1 0 1 −1 ︸ ︷︷ ︸ [Lt] _ h1(t−∆t) _ h2(t−∆t) _ h3(t−∆t) [ _ h(t) ] = [ k∧ ] [ _ h(t−∆t) ] + [ g∧ ] [ L ] [ G ]−1 [ L ]t [ _ h(t−∆t) ] (3.49) The discrete time transition matrix [A]d is then given by [ A ]d = [ k∧ ] + [ g∧ ] [ L ] [ G−1 ] [ Lt ] (3.50) which gives the connection between the [G] matrix based EMTP solution and the [A] based state space description. This connection had not been proposed previuosly to this work. Some interesting properties can be observed: 54 Chapter 3. Eigenvalue Analysis with EMTP Solution • [A]d depends directly on [G−1] • [L] is the graph connectivity • [k∧] is a diagonal matrix that indicates the nature of the branch, i.e., +1 for an L, −1 for a C or 0 for an R. • [g∧] is the diagonal matrix of the branch equivalent parameters as de- termined by the discretization rule and the time step ∆t. 3.4.6 Treatment of Parallel Branches For the case of L and C elements in parallel, see Fig. 3.12, we can choose from two possible approaches. First we can keep the identity of each compo- nent and by augmenting the size of the [A]d matrix compute each element’s eigenvalue. This alternative is convenient for the case in which it is neces- sary to identify, analyze and specify the corrective action for a single element contributing to a specific instability. Alternatively, we can produce an aggregation of elements of the same type to obtain an equivalent Req, Leq and Ceq, see Fig. 3.12. By doing so we are defining a new state and are able to compute its eigenvalue, which reflects the dynamic characteristic of the aggregate elements. Computationally, either option still has only one matrix inversion [G−1]. The only difference is the size of the resulting [A]d matrix. 55 Chapter 3. Eigenvalue Analysis with EMTP Solution 2 L 1 C 2 1 L 3 C 4 L5 CeqLeq 2 1 Figure 3.12: Aggregate treatment of parallel RLC branches _ h1(t) _ h2(t) _ h3(t) _ h4(t) _ h5(t) = k∧1 0 0 0 0 0 k∧2 0 0 0 0 0 k∧3 0 0 0 0 0 k∧4 0 0 0 0 0 k∧5 _ h1(t−∆t) _ h2(t−∆t) _ h3(t−∆t) _ h4(t−∆t) _ h5(t−∆t) + g∧1 0 0 0 0 0 g∧2 0 0 0 0 0 g∧3 0 0 0 0 0 g∧4 0 0 0 0 0 g∧5 _ v 1(t−∆t) _ v 2(t−∆t) _ v 3(t−∆t) _ v 4(t−∆t) _ v 5(t−∆t) (3.51) 56 Chapter 3. Eigenvalue Analysis with EMTP Solution _ v 1(t−∆t) _ v 2(t−∆t) _ v 3(t−∆t) _ v 4(t−∆t) _ v 5(t−∆t) = 1 −1 1 −1 1 −1 1 −1 1 −1 v1(t−∆t) v2(t−∆t) (3.52) The general formulas are maintained and finally we can obtain [A]d with equation (3.50), as for the case of a single branch per node case, [ A ]d = [ k∧ ] + [ g∧ ] [ L ] [ G−1 ] [ Lt ] 3.5 Nonlinear Networks Over the years the EMTP has used three different approaches to handle nonlinear elements: • current-source representation with time delay • compensation method • piecewise linear representation The first one is no longer used, while the second and third ones are currently used. In this work the piecewise linear representation is chosen. In this representation, nonlinear elements are considered as made up of piecewise linear segments [5]. Depending on the particular component, a 57 Chapter 3. Eigenvalue Analysis with EMTP Solution λ Saturation λ Saturation tan α 2 = L 2 tan α 1 = L 1 Figure 3.13: Piecewise linear inductance with two slopes piecewise segment can be relatively long (like the saturation characteristic of power transformers) or it can be as short as one segment per solution time step. A change of piecewise segment in the EMTP solution corresponds to a change in the value of the component (around the particular operating point), and therefore, to a new set of eigenvalue/eigenvectors. The dynamics of these eigensets reflect the system behaviour in terms of trajectories (where the system is going towards and how fast) and stability. Figure 3.53 illustrates how a nonlinear inductor is represented with two slopes [5]. This model can be implemented in discrete time simulation with L1 and Lp in parallel (3.53), as shown in Fig. 3.14. The flux in Lp is always 58 Chapter 3. Eigenvalue Analysis with EMTP Solution computed by integrating the voltage across the branch independent of the switch position. In this representation a change in the operating region cor- responds to a change in the topology of the network. During the simulation Lp will be connected whenever |λ| > λSaturation and disconnected as soon as |λ| < λSaturation. This is accurately enough in most cases to represent the saturation characteristics of modern transformers. 1 L2 = 1 L1 + 1 Lp (3.53) L p L 1 λ 1 λ p Figure 3.14: Implementation of two slope piecewise linear inductance 3.6 Continuous Time Eigenvalues from Discrete Time Domain Solution When eigenvalues are computed from the EMTP solution, the matrices in- volved contain only real numbers. Eigenvalue extraction routines, like the 59 Chapter 3. Eigenvalue Analysis with EMTP Solution QR method (Appendix A) are much more efficient for real matrices than for complex matrices. The process of extracting eigenvalues from the discrete time EMTP solution and reconstructing the continuous time eigenvalues from the inverse discretization mapping (3.64) can actually be much more efficient, particularly for large systems, than the direct computation of complex eigen- values from steady state power system admittance matrices whose elements are complex numbers. Notice also that in the important case of eigenvalue calculations with non- linearities in the circuit, the eigenvalues will change incrementally from one solution step to the next. Using the previously calculated eigenvalue set as the starting point for the new set can also speed up considerably the eigenvalue calculation. This situation is particularly advantageous for the Restarted Arnoldi method. 3.6.1 Discrete to Continuous Time Mapping The presented methodology allows for computing the discrete time eigen- values and eigenvectors of a given power system described with the EMTP solution. Even though the eigenvalues obtained in the discrete time solution are affected by the inaccuracies of the distortion introduced by the discretiza- tion rule (Fig.s 3.3 and 3.4), once the discrete time eigenvalues are obtained, it is possible to reconstruct the exact eigenvalues and eigenvectors of the orig- inal continuous time system. Similarly, the exact continuous time transition matrix can be built. 60 Chapter 3. Eigenvalue Analysis with EMTP Solution Consider the homogeneous part of the continuous time state space equations, ẋ(t) = Acx(t) (3.54) and its corresponding discrete time state space equivalent given by 3.22, x(t+∆t) = Adx(t) (3.55) Discretizing 3.54 with the trapezoidal integration rule, ẋ(t+∆t) = 2 ∆t x(t+∆t)− 2 ∆t h(t+∆t) (3.56) where h(t+∆t) = x(t) + ∆t 2 ẋ(t) (3.57) substituting 3.57 in 3.56, ẋ(t+∆t) = 2 ∆t {x(t+∆t)− x(t)} − ẋ(t) (3.58) where ẋ(t) = Acx(t) and ẋ(t+∆t) = Acx(t+∆t), Acx(t+∆t) = 2 ∆t {x(t+∆t)− x(t)} − Acx(t) (3.59) 61 Chapter 3. Eigenvalue Analysis with EMTP Solution solving for x(t+∆t), x(t+∆t) = x(t) 2 ∆t I − Ac 2 ∆t I + Ac (3.60) since, x(t+∆t) = Adx(t) (3.61) a relationship between Ad and Ac for the case of discretizing with Trapezoidal, can be established as Ad = ΩI + Ac ΩI − Ac (3.62) where Ω = 2/∆t. The eigenvalues and eigenvectors of Ac are defined by Acvi = λivi (3.63) The continuous time eigenvalues λi of A d can be computed immediately from those discrete time eigenvalues zi of A d, and vice-versa, as follows, λi = Ω zi − 1 zi + 1 (3.64) zi = Ω 1 + λi 1− λi (3.65) Equation (3.65) corresponds to the classical Bilinear transformation in signal processing theory. In the case of the trapezoidal rule (Bilinear transforma- 62 Chapter 3. Eigenvalue Analysis with EMTP Solution tion), the eigenvalues located anywhere in the entire left-hand side of the complex plane in the continuous time domain maps into the entire unit cir- cle in the discrete time domain as shown in Fig. 3.15. In the case when backward Euler is used as the integration rule, the relation between continuous and discrete time eigenvalues, as depicted in Fig. 3.16, can be expressed as λi = Ω zi − 1 zi (3.66) In the case of backward Euler the left hand side of the complex plane is mapped into a subregion of the unit circle in the discrete time domain. In general, for a stable discretization rule the left half of the complex plane in the continuous time maps into a subarea inside the unit circle in the discrete time domain. 3.6.2 Time Discretization Considerations The linearization of differential equations into in difference equations is an approximation that will produce accurate results when the highest frequency present in the system is 5 to 10 times below the Nyquist frequency, as dis- cussed in section 3.2.2. However, when using the discrete time formulation of a network to compute the eigenvalues some extra considerations have to be in place. As discussed in section 3.6.1, the accuracy of the discrete to continuous time mapping in LTI systems will be correct until the Nyquist frequency. 63 Chapter 3. Eigenvalue Analysis with EMTP Solution Real Real Imaginary Imaginary λ λ z z Complex continuous-time plane Complex discrete-time plane Unit circle Figure 3.15: Continuous time and discrete time stability mapping for trape- zoidal Real Real Imaginary Imaginary λ λ z z Complex continuous-time plane Complex discrete-time plane Unit circle Figure 3.16: Continuous time and discrete time stability mapping for back- ward Euler 64 Chapter 3. Eigenvalue Analysis with EMTP Solution However, in the case of systems with nonlinear elements the discrete to con- tinuous time mapping will be accurate only if the discretization time is small enough to represent correctly the region change of the nonlinear elements. In addition to the limitations imposed by nonlinear elements, using smaller discretization times than required will not lead to further improvements in the accuracy of the continuous time eigenvalues. However, a nonlinear scaling effect will be present in the discrete time values due to the mapping charac- teristics of the discretization transformation. The smaller the discretization time, the closer to the unit circle the discrete time eigenvalues will be located. In the continuous time domain, the behaviour of the term associated to an eigenvalue λ is given by eλt. Consequently, • A purely real eigenvalue corresponds to a non-oscillatory mode. – A positive real eigenvalue represents an exponential growing in- stability. – A negative real eigenvalue represents an exponential decaying mode. – A zero eigenvalue represents a constant mode. • Complex eigenvalues occur in conjugate pairs, and each pair corre- sponds to an oscillatory mode with an exponential amplitude with as before can be decaying in the stable case or exponentially growing in the unstable case. In the critically damped case the amplitude becomes constant. 65 Chapter 4 Application to the OVNI simulator 4.1 Introduction This chapter starts with an overview of the OVNI simulator and its associated MATE solution framework. Also a background of OVNI’s architecture for real-time simulation is provided. Next follows a description of the procedure to obtain the state space formulation to compute the discrete time eigenvalues for large networks from MATE. The chapter concludes with two extended applications of the new methodology for the Latency concept [55], [56] and for the automatic discretization step selection in the EMTP. The research contributions included in this chapter are: 1. Eigenvalue extraction from OVNI solution. 2. Eigenvalue trajectories visualization in discrete time for the purpose of assessing dynamic behaviour. 3. Participation factors computation from discrete time eigenvalues. 66 Chapter 4. Application to the OVNI simulator 4. Identification of segmentation schemes from discrete time eigenvalue information applied to Latency. 5. Automatic selection of variable discretization steps from discrete time eigenvalue information. 4.2 OVNI The Object Virtual Network Integrator (OVNI) is UBC’s power system real- time simulator. OVNI’s discrete time solution is based on the EMTP dis- cretization method. It was developed by Mart́ı and Linares to simulate large electric networks in real-time and its software and hardware architectures are described in the following publications [10], [18], [13], [19], [21], [22] and [57], among others. The OVNI simulator is built around the Multi Area Thévenin Equivalent (MATE) concept, which is described in the following section, and was developed as an object oriented application written in C++. This char- acteristic facilitates the integration of new models and modules with OVNI. The OVNI solution achieved real-time performance in a single PC environ- ment allowing interaction with analog devices for protective relay testing ap- plications, as introduced by Calviño in [16]. The OVNI PC-cluster solution developed by Hollman [17] showed the OVNI algorithm’s inherent advan- tages for parallelization allowing for real-time simulation of large networks. The introduction of the Multilevel MATE concept by Armstrong [57] pro- vides OVNI with the capability for efficient simultaneous solution of control 67 Chapter 4. Application to the OVNI simulator + - iα z Vα BA Figure 4.1: Generic description of a MATE link systems and nonlinearities. 4.2.1 The MATE Concept The MATE concept developed by Mart́ı [58], [13], [19] was introduced to cope with the dimension of large power system networks for real-time and on-line dynamic studies. It also provides a framework for the implementation of the Latency concept which will be summarized later in this chapter. MATE is used in OVNI to segment a complex network into independent subnetworks interconnected by links. A link can be a transmission line, a voltage source or just simply an ideal link to allow the desired segmentation. The most general description of a link is shown in Fig. 4.1. To demonstrate the MATE concept, let us consider the system illustrated in Fig. 4.2. Where [A] is the admittance matrix of subnetwork [A] gA12 is the branch admittance between nodes 1 and 2 GB23 is the element (2, 3) of the admittance matrix [B] 68 Chapter 4. Application to the OVNI simulator hA1 is the current injected from ground to node 1 of subnetwork [A] VA3 is the voltage of node 3 to ground in subnetwork [A] iα1 is the current in link branch α1 zα1 is the impedance of branch α1 h A2 A 2 A 1 A 3 A 4 B 3 B 2 B 1 h B2 h A4 h A1 g A23 g A34 g A14 g A21 h B1 g B13 g B12 g B32 [ A ] [ B ] Link α 1 Link α 2 iα 1 iα 2 Zα 2 Zα 1 Figure 4.2: Generic circuit for discrete state space formulation with MATE For the considered network, the general nodal equation of the system with the links α1 and α2 removed is, A 0 0 B vA vB = hA hB (4.1) 69 Chapter 4. Application to the OVNI simulator while the full system with the links α1 and α2 included is given by GA11 GA12 0 GA14 0 0 GA21 GA22 GA23 0 0 GA31 GA32 GA33 GA34 1 0 GA41 0 GA34 GA44 0 1 GB11 GB12 GB13 −1 0 GB21 GB22 GB23 0 0 GB31 GB32 GB33 0 −1 0 0 1 0 −1 0 0 −z1 0 0 0 0 1 0 0 −1 0 −z2 vA1 vA2 vA3 vA4 vB1 vB2 vB3 iα1 iα2 = hA1 hA2 hA3 hA4 hB1 hB2 hB3 0 0 (4.2) 70 Chapter 4. Application to the OVNI simulator with its compact notation, A 0 p 0 B q pt qt −z vA vB iα = hA hB 0 (4.3) Where [p] and [q] are the injection matrices derived from the link branches (if the current comes out +1, if the current goes in −1). p = 0 0 0 0 +1 0 0 +1 (4.4) q = −1 0 0 0 0 −1 (4.5) After manipulations, as explained in [21], equation (4.2) can be represented by (4.6), the MATE equivalent circuit which is implemented in the central 71 Chapter 4. Application to the OVNI simulator solver, 1 a13 a14 1 a23 a24 1 a33 a34 1 a43 a44 1 −b11 −b13 1 −b21 −b23 1 −b31 −b33 Zα11 Zα12 Zα21 Zα22 vA1 vA2 vA3 vA4 vB1 vB2 vB3 iα1 iα2 = eA1 eA2 eA3 eA4 eB1 eB2 eB3 eα1 eα2 (4.6) where, a = A−1p (4.7) b = B−1q (4.8) eA = A −1hA (4.9) eB = B −1hB (4.10) e = pteA + q teB (4.11) Z = pta+ qtb+ z (4.12) 72 Chapter 4. Application to the OVNI simulator 4.2.2 Discrete State Space Formulation with MATE MATE allows for the implementation of the proposed new methodology to extract eigenvalues from the discrete time EMTP solution. The formulation of the state space description can be implemented at a subnetwork level, while the global analysis of the eigenvalues is performed at the core level either in a central or dedicated node of the hardware cluster. Each sub- network node contains all the information required to formulate the discrete state space equations for that subnetwork. Each subnetwork has to submit to the central solver its discrete state space description, which is then relayed to the dedicated eigenvalue solver. Since the central node solver is often currently used to preprocess the data input file, the process of constructing the discrete state space formulation of the network can be achieved at once, together with the case preprocessing stage. Consequently, after the simula- tion is started, only those subnetworks in which there is a topology change or there is a change of the nonlinear region will have to submit to the cen- tral solver the updated discrete state space formulation. This is particularly efficient for the cases where the nonlinearities are isolated in subnetworks, since the computation size of the problem is reduced (see Fig. 4.3). Here the subnetworks A2 and A3 are the only ones in need to submit its state space description update to the dedicated eigenvalue solver. With all the system information allocated in the dedicated eigenvalue solver node, the discrete transition matrix [A]d of the complete network can be computed. The con- tents of the segmented blocks’ transition matrices are located in the correct 73 Chapter 4. Application to the OVNI simulator Linear subnetwork A1 Linear subnetwork A2 with topology change nonlinear subnetwork A3 Eigenvalue solver Update is required only when change of region Update is required only when change of topology No update required Figure 4.3: General update scheme for discrete state space formulation with MATE place of the global transition matrix using the information of connectivity present within the class elm t [19] For the case of the sample network of Fig. 4.2, the incidence matrices [LA], [LB] and the coefficients matrices [k∧A], [k∧B], [g∧A] and [g∧B] are constructed from the data input case file. The subnetwork’s discrete transition matrices [AA] d and [AB] d are first constructed from the data input case, and later dur- ing the simulation are updated with the information from the subnetwork in case of change of topology or change of region of a nonlinear subnetwork. The new hardware network interconnect layout of OVNI [22] can smoothly 74 Chapter 4. Application to the OVNI simulator run simultaneously both hard and soft real-time environments. Thus, it is possible to allocate all the required eigenvalue analysis of the electrical net- work within the soft real-time environment and interacting with the hard real-time cluster when required, as depicted in Fig. 4.4. These interactions exist whenever the topology of the network changes or whenever the dynam- ics of elements is updated in the stability analysis [59], [60]. As an alternative to the proposed centralized computation of eigenvalue Subsystem I OVNI central solver Subsystem III Subsystem II Dedicated Eigenvalue solver Real-time environment Soft real-time environment Interaction required only when there is a change in [A] i Interactions required at every time step links links links Figure 4.4: Hybrid real-time/soft real-time simulator layout trajectories, a segmented distributed computation can be implemented ex- 75 Chapter 4. Application to the OVNI simulator tending the MATE partitioning concept [13], [19] to the eigenvalue analysis. Computing the eigenvalues at the tearing blocks level can provide a better computational performance in some cases for large scale systems by concen- trating the computation of eigenvalues within a much smaller structure. In an extension of the Latency concept discussed in the next section, the sub- system eigenvalues can be determined by its MATE Thevenin equivalent, as seen from the links subsystem. This approach, however, will only produce good results when the Latency approach is effective, which is when there is a large time constant separation among subsystems. The Latency approach is discussed next. 4.3 The Latency Concept The Latency concept3, introduced by Mart́ı and Linares in [55] and by Mor- eira in [56], aims at exploiting the use of different time steps in different subnetworks. By using larger time steps for slow subnetworks and small time steps for fast subnetworks considerable computational savings can be achieved without sacrifying accuracy. The Latency concept [20], [56] depends for its successful implementation on the accurate identification of areas with significantly different time constants. This is not by any means an easy task, and as quoted in [56] “it is one of the limitations of the Latency method...”. 3In this work, like in previous power system references, the term “latency” is used to mean “latent” or “dormant”. In Signal processing theory this concept corresponds to multirate systems 76 Chapter 4. Application to the OVNI simulator 4.3.1 Application of the New Discrete Time State Space Formulation to Latency The methodology presented in this work allows for the identification of dom- inant eigenvalues in the discrete time domain using the EMTP description. Each subnetwork can be described in the state space domain by equation (3.22) in which the eigenvalues of [A] determine directly the time-response of that subnetwork. After the system eigenvalues and participation matrix are computed an efficient segmentation and correct time step selection can be achieved. In order to produce a correct segmentation, the concept of the participation matrix [3] is used. The Participation matrix [P ], is given by P = [p̃1p̃2 . . . p̃k] (4.13) with p̃i = φ1iϕi1 φ2iϕi2 ... φniϕik (4.14) where φki = k th entry of the right eigenvector φi ϕik = k th entry of the left eigenvector ϕi. The dimensionless element Pki provides a measure of the relative partic- 77 Chapter 4. Application to the OVNI simulator ipation of the kth state in the ith mode. The participation matrix combines right and left eigenvectors showing the association between state space vari- ables and modes. The relationship between states and modes can not easily be extracted directly from the right and left eigenvectors due to their depen- dency on scaling. The Participation matrix overcomes this situation. The algorithm can be summarized as follows: i. obtain the discrete time state space formulation of the network from the EMTP solution, that is, compute [A]d with equation (3.50) ii. compute the discrete time eigenvalues of the system from [A]d iii. reconstruct the continuous time eigenvalues with equation (3.64) in the case of the trapezoidal integration rule iv. compute the Participation matrix [P ] from equation (4.13) v. identify in [P ], clusters of elements sharing similar time constants vi. apply the latency method to partition the system with the appropriate time steps This procedure should be repeated whenever there is a change in the topology of the system. In section 5.2 a detailed example is provided to illustrate this methodology. 78 Chapter 4. Application to the OVNI simulator 4.4 Automatic Time Step Selection With the eigenvalue information obtained from the discrete state space for- mulation from the EMTP solution, it is possible to implement an automatic time step selection option in the EMTP. Assume the user specifies a ∆t which is adequate for the maximum frequency of interest in the study. The exact continuous time eigenvalues can be obtained from the discrete time EMTP formulation as discussed in this Thesis. From the exact continuous time eigenvalues the maximum frequency present in the system (below the maximum study frequency) can be determined. If this frequency is higher than say a fifth of the Nyquist frequency then a smaller ∆t that satisfies this condition must be selected. The automatic selection of the simulation time step must be performed every time a nonlinear element changes its operation region or every time a topology change is applied. The general scheme is illustrated in Fig. 4.5. 79 Chapter 4. Application to the OVNI simulator Maximum frequency of interest in the study, specified by user Time step specified by user in data input file Discrete time eigenvalues computation from EMTP solution Maximum frequency is identified If f max system > f max user Issue user notification Is f max system of interest? no yes Run simulation with fixed user time step Select best appropriate time step Run simulation Change of topology or nonlinear operation region? no yes Discrete time eigenvalues computation from EMTP solution Figure 4.5: Automatic EMTP time step selection scheme 80 Chapter 5 Validation 5.1 Comparison of State Space Solutions between Continuous Time and Discrete Time Domains For comparison purposes and a detailed illustration of the proposed tech- nique, a state space formulation of a simple RLC series circuit, see Fig. 5.1, will be obtained using both the traditional approach from the continuous time domain and the proposed new technique from the nodal EMTP matrix. This second order circuit is an interesting initial case, because as pointed out in [3], the performance of higher order systems is often viewed in terms of a dominant second order system. 81 Chapter 5. Validation R C L V in (t) i(t) v C (t) v L (t) Figure 5.1: RLC series circuit 5.1.1 RLC Series Circuit State Space Equation from Continuous Time Domain The continuous time state space system equations are given by ẋ(t) = A x(t) +B u(t) (5.1) y(t) = C x(t) +D u(t) (5.2) Applying Kirchoff Voltage Law to the circuit vin = eL + eR + eC (5.3) where eL = L di dt ; eR = R i ; eC = 1 C ∫ i(t) dt 82 Chapter 5. Validation i(t) is the current in the loop. Substituting vin(t) = L di(t) dt +R i(t) + 1 C ∫ i(t) dt (5.4) Introducing the capacitor charge as a variable, i(t) = dq(t)/dt, vin(t) = L d2q(t) dt2 +R dq(t) dt + 1 C q(t) (5.5) If we select q(t) and i(t) as the state variables and define the input and output as x(t) = q(t) i(t) u(t) = vin(t) (5.6) y(t) = vL(t) then the state space equations of the RLC series circuit are dq(t)dt di(t) dt ︸ ︷︷ ︸ rate of change = 0 1 −1 LC −R L ︸ ︷︷ ︸ [A] q(t) i(t) ︸ ︷︷ ︸ present states + 0 1 L ︸︷︷︸ [B] vin(t)︸ ︷︷ ︸ input (5.7) vL(t)︸ ︷︷ ︸ output = [ −1 C −R ] ︸ ︷︷ ︸ [C] q(t) i(t) ︸ ︷︷ ︸ states + 1︸︷︷︸ [D] vin(t)︸ ︷︷ ︸ input (5.8) 83 Chapter 5. Validation If we assign the following values to the circuit shown in Fig. 5.1, R = 100Ω, L = 5H, C = 250µF , V = 8.5V , the obtained continuous time transition matrix [A]c and the eigenvalues [λ]c are [ A ]c = 0 1 −800 −20 (5.9) [ λ ]c = −10.0000 + j26.4575 −10.0000− j26.4575 (5.10) Another possible representation of the system can be obtained by se- lecting vR(t) and vC(t) as the state variables. In this case the state space equations for the RLC series circuit are x(t) = vR(t) vC(t) u(t) = vin(t) (5.11) y(t) = vL(t) dvR(t)dt dvC(t) dt ︸ ︷︷ ︸ rate of change = −RL −RL −1 RC 0 ︸ ︷︷ ︸ [A] vR(t) vC(t) ︸ ︷︷ ︸ states + RL 0 ︸︷︷︸ [B] vin(t)︸ ︷︷ ︸ input (5.12) 84 Chapter 5. Validation vL(t)︸ ︷︷ ︸ output = [ −1 −1 ] ︸ ︷︷ ︸ [C] vR(t) vC(t) ︸ ︷︷ ︸ states + 1︸︷︷︸ [D] vin(t)︸ ︷︷ ︸ input (5.13) The obtained continuous time transition matrix [A]c and the eigenvalues [λ]c are [ A ]c = −20 −20 40 0 (5.14) [ λ ]c = −10.0000 + j26.4575 −10.0000− j26.4575 (5.15) We can see that the continuous time eigenvalues [λ]c are the same in both cases (5.10), (5.15). The transition matrix [A]c is not unique, (5.9) (5.14), but the eigenvalues [λ]c are. 5.1.2 RLC Series Circuit State Space Equation from the [G] EMTP Matrix The EMTP equivalent discretized with the trapezoidal rule RLC series circuit is shown in Fig. 5.2. The nodal equations in matrix form are 85 Chapter 5. Validation h C (t) h L (t) 32 g C equiv g L equiv V in (t) 1 R Figure 5.2: Discretized RLC series circuit 1 R −1 R 0 −1 R ( 1 R + −∆t 2L ) −∆t 2L 0 ∆t 2L (∆t 2L + 2C ∆t ) v1 v2 v3 = i1 hL −hL + hC (5.16) Since v1 = vin and is a known voltage source, we reduce the order of the system by 1 by moving the node to the right hand side of the equation. Re-arranging the [G] matrix 86 Chapter 5. Validation ( 1 R + ∆t 2L ) −∆t 2L −1 R −∆t 2L (∆t 2L + 2C ∆t ) 0 −1 R 0 1 R v2 v3 v1 = hL −hL + hC i1 (5.17) GAA GAB GBA GBB vA vB = hA hB (5.18) ( 1R + ∆t2L) −∆t2L −∆t 2L (∆t 2L + 2C ∆t ) ︸ ︷︷ ︸ [G] v2 v3 = hL −hL + hC ︸ ︷︷ ︸ [h] − −v1R 0 ︸ ︷︷ ︸ input (5.19) v2 v3 = ( 1R + ∆t2L) −∆t2L −∆t 2L (∆t 2L + 2C ∆t ) −1 hL + v1R −hL + hC (5.20) The relation between branch and node voltages and history sources is given by the incidence matrix [L]: _ v′1 _ v′2 _ v′3 = 1 0 −1 1 0 1 v′2 v′3 (5.21) h′2 h′3 = 1 −1 0 0 1 1 h′R −h′L hC (5.22) 87 Chapter 5. Validation Substituting the nodal voltages (5.20) into (5.21), _ v′1 _ v′2 _ v′3 = 1 0 −1 1 0 1 ( 1R + ∆t2L) −∆t2L −∆t 2L (∆t 2L + 2C ∆t ) −1 h′L + v1R −h′L + h′C (5.23) For the considered RLC series circuit, [k∧] and [g∧] (3.45) are [ k∧ ] = 0 0 0 0 1 0 0 0 −1 (5.24) [ g∧ ] = 0 0 0 0 −∆t L 0 0 0 4C ∆t (5.25) while the branch histories, according to equation (3.49), can be obtained as _ h1 _ h2 _ h3 = 0 0 0 0 1 0 0 0 −1 _ h′1 _ h′2 _ h′3 + 0 0 0 0 −∆t L 0 0 0 4C ∆t 1 0 0 1 1 −1 ( 1R + ∆t2L) −∆t2L −∆t 2L (∆t 2L + 2C ∆t ) −1 1 0 1 0 1 −1 _ h′1 _ h′2 _ h′3 (5.26) Finally, the discrete time transition matrix [A]d for the circuit is computed 88 Chapter 5. Validation from the nodal [G] matrix as indicated in (3.50), [ A ]d = 0 0 0 0 1 0 0 0 −1 + 0 0 0 0 −∆t L 0 0 0 4C ∆t 1 0 0 1 1 −1 ( 1R + ∆t2L) −∆t2L −∆t 2L (∆t 2L + 2C ∆t ) −1 1 0 1 0 1 −1 (5.27) Assuming the same circuit values and discretization time of 50µS: R L C V ∆t 100Ω 5H 250µF 8.5V 50µS we obtain the following discrete time transition matrix [ A ]d = 0 0 0 1.99799−3 9.97998−1 −3.99599−6 1.99799−3 1.99799 9.99996−1 (5.28) The discrete time eigenvalues of this matrix are [ z ]d = 9.98997−1 + j2.64310−3 9.98997−1 − j2.64310−3 0 (5.29) The eigenvalues of the transition matrix [A] are the poles of the corre- sponding linear time-invariant system. In order to validate the proposed 89 Chapter 5. Validation technique, it is desirable to evaluate the precision of the eigenvalues com- puted from the nodal [G] matrix against those obtained from the continuous time domain. For the case of the trapezoidal integration rule, the continuous time eigenvalues can be computed from the discrete time transition matrix obtained from the nodal [G] matrix with the proposed method using equation (3.64), λi = 2 ∆t (zi − 1) (zi + 1) (5.30) where λi are the continuous time eigenvalues and zi are the discrete time eigenvalues . For the RLC circuit, the reconstructed continuous time eigen- values from the discrete time eigenvalues, obtained with equation (3.64), with a ∆t of 100µs are [ λ ]c reconstructed = 2 100µs −0.00050 + j0.00132 −0.00050− j0.00132 = −10.0000 + j26.4575 −10.0000− j26.4575 (5.31) The continuous time eigenvalues of equation (5.10) and the reconstructed continuous time eigenvalues of equation (5.31) are exactly the same, thus validating the proposed technique. The example validates the discussion in chapter 3 that the exact eigenvalues of a continuous time system can be 90 Chapter 5. Validation exactly computed from the discrete time nodal method independently from the chosen integration time step. It is important to note that even when a very large discretization time step is selected, the continuous time eigenvalues can still be recomputed “exactly”, with possibly some small error due to truncation. For instance, if we take ∆t = 1ms [ λ ]c reconstructed = 2 1ms −0.0050 + j0.01322 −0.0050− j0.01322 = −10.0000 + j26.4400 −10.0000− j26.4400 (5.32) Even though the eigenvalues can still be reconstructed for large ∆t′s, in the important case of a circuit with nonlinearities, large ∆t′s may move the operating point too far in the nonlinear characteristic as we move from the solution step to the next. In this case, we may miss intermediate critical points in the operating characteristics. The real part of the eigenvalues gives the damping, while the imaginary part gives the frequency of oscillation. The circuit’s oscillation frequency f and damping ratio ζ can be computed from the reconstructed continuous time eigenvalues (5.31) as follows, let λ = σ ± jω (5.33) 91 Chapter 5. Validation then f = ω 2pi (5.34) and ζ = −σ√ σ2 + ω2 (5.35) A positive damping ratio ζ provides damping to the natural frequency. If ζ is greater than 1, both eigenvalues are real and negative; if ζ is equal to 1, both eigenvalues are equal to −wn; and if ζ is less than 1, the eigenvalues are complex and conjugates. For the series RLC test case the reconstructed continuous oscillation fre- quency and damping ratio are f c = 4.21Hz (5.36) ζc = 0.35355 (5.37) The LC resonant frequency of the circuit f cn (“natural frequency”) is f cn = 1 2pi √ LC = 4.50Hz (5.38) and the damping ratio ζc is ζc = R 2 √ L C = 0.35355 (5.39) Notice that in a stable system due to the losses the oscillation frequency 92 Chapter 5. Validation Vgen L1 C1 C2 R1 L2 State 1 State 4 State 3 State 2 Figure 5.3: Latency test circuit (5.34) is lower than the natural frequency (5.38). 5.2 Latency Case In this section the simple test case shown in Fig. 5.3 from [56], is used to demonstrate the capabilities of the presented method for an efficient imple- mentation of Latency. The values of the circuit components are R1 L1 C1 L2 C2 0.1Ω 1µH 100µF 1µH 1µF For validation purposes the eigenvalues will be also computed from the con- tinuous time description. Selecting the capacitor voltages and inductor currents as state variables, we 93 Chapter 5. Validation obtain the following state space equations describing the test circuit, iC1 = C1 dvC1 dt =⇒ ẋ1 = 1C1 iC1 vL1 = L1 diL1 dt =⇒ ẋ2 = 1L1vL1 iC2 = C2 dvC2 dt =⇒ ẋ3 = 1C2 iC2 vL2 = L2 diL2 dt =⇒ ẋ4 = 1L2 iL2 (5.40) Applying Kirchoff’s laws, ẋ1 = 1 C1 (vL1 − vL2) ẋ2 = 1 L1 (Vgen − vC1) ẋ3 = 1 C2 vL2 ẋ4 = 1 L2 (vC1 − vC2 +R1vL2) (5.41) The continuous time [A] matrix for the test circuit is given by Ac = 0 1 C1 0 −1 C1 −1 L1 0 0 0 0 0 0 1 C2 1 L2 0 −1 L2 −R1 L2 (5.42) 94 Chapter 5. Validation from which the following eigenvalues are obtained λc = −4.99950x104 + j1.00379x106 λ1 −4.99950x104 − j1.00379x106 λ2 −4.99801 + j9.94988x104 λ3 −4.99801− j9.94988x104 λ4 (5.43) We next compute the eigenvalues and eigenvectors using the discrete time methodology. The network’s [G] matrix from the discrete time description is G = ∆t 2L1 −∆t 2L1 0 0 −∆t 2L1 ( 1 R1 + ∆t 2L1 + 2C1 ∆t ) −1 R1 0 0 −1 R1 ( 1 R1 + ∆t 2L2 ) −∆t 2L2 0 0 −∆t 2L2 ( ∆t 2L2 + 2C2 ∆t ) (5.44) the matrix [k∧] is given by k∧ = 1 0 0 0 0 −1 0 0 0 0 1 0 0 0 0 −1 (5.45) 95 Chapter 5. Validation while the matrix [g∧] is given by g∧ = −∆t L1 0 0 0 0 2C1 ∆t 0 0 0 0 −∆t L2 0 0 0 0 4C2 ∆t (5.46) The incidence matrix [L] L = 1 −1 0 0 0 1 0 0 0 0 1 −1 0 0 0 1 (5.47) Applying (3.50), we obtain the discrete time transition matrix [A]d [ A ]d = [ k∧ ] + [ g∧ ] [ L ] [ G−1 ] [ Lt ] Ad = 9.99999x10−1 1.25000x10−15 1.24999x10−15 0 −2.00000 9.99999x10−1 1.99999 1.25034x10−13 1.24999x10−15 −1.24999x10−15 9.99999x10−1 1.24999x10−13 −1.24999x10−15 1.24999x10−15 −1.99999 9.99999x10−1 (5.48) 96 Chapter 5. Validation from where the following discrete time eigenvalues are computed, λd = 9.99999x10−1 + j4.97494x10−8 9.99999x10−1 + j4.97494x10−8 9.99999x10−1 + j5.01896x10−7 9.99999x10−1 − j5.01896x10−7 (5.49) Finally, the reconstructed eigenvalues from the discrete time domain given by (3.64) are, λcreconstructed = −4.99752 + j9.94988x104 −4.99752− j9.94988x104 −4.99950x104 + j1.00379x106 −4.99950x104 − j1.00379x106 (5.50) which matches up to the second decimal for the worse approximation the continuous time eigenvalues calculated directly from the continuous time so- lution (5.43). The participation matrix [P ] (4.13) can now be computed in order to determine how Latency will be implemented. For the test circuit, 97 Chapter 5. Validation the participation matrix [P ] computed from the discrete time domain is P crec. = λ3 c reconstructed λ4 c reconstructed 4.9995x10−1 + j2.6108x10−5 4.9995x10−1 − j2.6108x10−5 ẋ2 4.9495x10−1 + j7.5560x10−5 4.9495x10−1 − j7.5560x10−5 ẋ1 4.9970x10−5 − j9.9168x10−7 4.9970x10−5 + j9.9168x10−7 ẋ4 5.0474x10−3 − j1.0067x10−4 5.0474x10−3 + j1.0067x10−4 ẋ3 λ2 c reconstructed λ1 c reconstructed 4.9470x10−5 − j7.5406x10−6 4.9470x10−5 + j7.5406x10−6 ẋ2 5.0479x10−3 − j2.6137x10−4 5.0479x10−3 + j2.6137x10−4 ẋ1 4.9995x10−1 + j2.4910x10−2 4.9995x10−1 − j2.4910x10−2 ẋ4 4.9495x10−1 − j2.4641x10−2 4.9495x10−1 + j2.4641x10−2 ẋ3 (5.51) As can be seen, the eigenvalues λ1 c reconstructed and λ2 c reconstructed are strongly associated to the states of ẋ3 and ẋ4 (corresponding to C2 and L2), while eigenvalues λ3 c reconstructed and λ4 c reconstructed are strongly associated to the states of ẋ1 and ẋ2 (corresponding to C1 and L1). The oscillation frequencies given by (5.34) are 15.8kHz for the slow sub- network and 159kHz for the fast one. Consequently, the time constants of λ1 and λ2 are very loosely coupled with λ3 and λ4 allowing for a suitable 98 Chapter 5. Validation Vgen State 1 State 4 State 3 State 2 L1 C1 C2 R1 L2 Slow Fast Figure 5.4: Circuit with slow and fast subareas segmentation into a fast subcircuit and a slow subcircuit to implement the Latency concept, as shown in Fig. 5.4. Vgen State 1 State 4 State 3State 2 R1 L1 C1 C2 R2 L2 SW1 Rload 1 2 3 4 Figure 5.5: Two-area resonant test case 5.3 Two-Area Resonant Circuit A two-area resonant circuit, Fig. 5.5, is presented here. The eigenvalue analysis is performed with the discrete time state space algorithm, for the 99 Chapter 5. Validation cases of a poorly damped and an improved damped configuration. The selected values for the simulation are, for the poorly damped case I, R1 R2 C1 C2 L1 L2 Vgen ∆t 3Ω 5Ω 10nF 50nF 350mH 10mH 230kV 50µs and for the improved damped case II, R1 R2 C1 C2 L1 L2 Vgen ∆t 3Ω 5Ω 10nF 1µF 350mH 1mH 230kV 50µs The source is applied at t = 0 with switch SW1 closed. At t = 25ms SW1 is opened and the two resonant circuits become isolated. 5.3.1 Case I - Poorly Damped The voltage outputs for the poorly damped caseI are shown in Fig. 5.6. Using the proposed methodology the trajectory of the eigenvalues of the test network is computed. The discrete eigenvalues are presented in 5.7 and the reconstructed continuous time eigenvalues are presented in Fig. 5.8. There are two sets of eigenvalues which correspond to the circuit topology before and after the switching operation. The [G] discretized network matrix with switch SW1 open is given by [Gopen] = 1 R1+2 L1 ∆t −1 R1+2 L1 ∆t 0 0 − 1 R1+2 L1 ∆t 1 R1+2 L1 ∆t + 2C1∆t + 1 Ropen −1 Ropen 0 0 −1Ropen 1 Ropen + 2C2∆t + 1 R2+2 L2 ∆t −1 R2+2 L2 ∆t 0 0 −1 R2+2 L2 ∆t 1 R2+2 L2 ∆t + 1Rload (5.52) 100 Chapter 5. Validation and with the switch SW1 closed, [Gclose] = 1 R1+2 L1 ∆t −1 R1+2 L1 ∆t 0 0 − 1 R1+2 L1 ∆t 1 R1+2 L1 ∆t + 2C1∆t + 1 Rclose −1 Rclose 0 0 −1Rclose 1 Rclose + 2C2∆t + 1 R2+2 L2 ∆t −1 R2+2 L2 ∆t 0 0 −1 R2+2 L2 ∆t 1 R2+2 L2 ∆t + 1Rload (5.53) the [k∧] matrix of branch nature coefficients is given by k∧ = 1 0 0 0 0 0 −1 0 0 0 0 0 −1 0 0 0 0 0 1 0 0 0 0 0 0 (5.54) while [g∧] matrix of branch discretization rule coefficients is g∧ = −2 R1+ 2L1 ∆t 0 0 0 0 0 4C1 ∆t 0 0 0 0 0 4C2 ∆t 0 0 0 0 0 −2 R2+ 2L2 ∆t 0 0 0 0 0 0 (5.55) 101 Chapter 5. Validation and the incidence matrix [L] is L = 1 −1 0 0 0 1 0 0 0 0 1 0 0 0 1 −1 0 0 0 1 (5.56) The transition matrices [A] for the open and closed switch conditions are obtained with 3.50. Finally the discrete time eigenvalues are λdclosed = −1.419878x10−2 + j9.998990x10−1 −1.419878x10−2 − j9.998990x10−1 −9.999999x10−1 −1.000000 (5.57) λdopen = −1.049723x10−1 + j9.944750x10−1 −1.049723x10−1 − j9.944750x10−1 9.999999x10−1 −1.000000 (5.58) while the reconstructed continuous time eigenvalues computed with (3.64) 102 Chapter 5. Validation 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 505 x 10 5 Vo lta ge N od e 1 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 101 x 10 6 Vo lta ge N od e 2 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 202 x 10 4 Vo lta ge N od e 3 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 0. 50 0. 5 Vo lta ge N od e 4 tim e [s] [V] Figure 5.6: Node voltages for the two-area example of Fig. 5.5. Case I, poorly damped 103 Chapter 5. Validation −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.01 0.02 0.03 0.04 0.05 Real Imag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 Figure 5.7: Discrete time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped are λcclosed = −5.072087x10−3 + j4.057204x104 −5.072087x10−3 − j4.057204x104 −1.200000x1012 −4.000000x104 (5.59) λcopen = −4.938281x10−3 + j4.444444x104 −4.938281x10−3 − j4.444444x104 −1.000044x10−7 −4.000000x104 (5.60) The continuous time damping ratios computed with (5.35) for both topolo- 104 Chapter 5. Validation −14 −12 −10 −8 −6 −4 −2 0 x 1011 −5 0 5 x 104 0 0.01 0.02 0.03 0.04 0.05 Real Imag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 Figure 5.8: Reconstructed continuous time eigenvalues for the two-area ex- ample of Fig. 5.5. Case I, poorly damped gies are ζclosed = 1.250x10 −7 ζopen = 1.111x10 −7 (5.61) and the damped oscillation frequencies (5.34) are fclosed = 6.457kHz fopen = 7.073kHz (5.62) As can be observed, the analyzed circuit is poorly damped for the oscillation 105 Chapter 5. Validation frequencies for both open switch and closed switch conditions. 5.3.2 Case II - Improved Damping The same circuit is once again analyzed with C2 and L2 changed in order to improve the damping of the oscillation frequency. The new set of eigenvalues and damping ratios for case II with improved damping are shown below. Discrete time eigenvalues, λdclosed = 2.902763x10−1 + j9.569413x10−1 2.902763x10−1 − j9.569413x10−1 −9.99999x10−1 −1.00000 (5.63) λdopen = 2.857138x10−1 + j9.583134x10−1 2.857138x10−1 − j9.583134x10−1 −9.99999x10−1 −1.00000 (5.64) while the reconstructed continuous eigenvalues are λcclosed = −4.444838x10−2 + j2.966627x104 −4.444838x10−2 − j2.966627x104 −1.009999x1012 −4.000000x104 (5.65) 106 Chapter 5. Validation 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 505 x 10 5 Vo lta ge N od e 1 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 101 x 10 6 Vo lta ge N od e 2 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 202 x 10 4 Vo lta ge N od e 3 tim e [s] [V] 0 0. 00 5 0. 01 0. 01 5 0. 02 0. 02 5 0. 03 0. 03 5 0. 04 0. 04 5 0. 05 − 0. 50 0. 5 Vo lta ge N od e 4 tim e [s] [V] Figure 5.9: Node voltages for the two-area example of Fig. 5.5. Case II, improved damping 107 Chapter 5. Validation −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 0.01 0.02 0.03 0.04 0.05 Real Imag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 Figure 5.10: Discrete time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping λcopen = −4.444444x10−2 + j2.981423x104 −4.444444x10−2 − j2.981423x104 1.000044x10−7 −4.000000x104 (5.66) The new damping ratios are ζclosed = 1.498x10 −6 ζopen = 1.490x10 −6 (5.67) 108 Chapter 5. Validation −12 −10 −8 −6 −4 −2 0 x 1011 −3 −2 −1 0 1 2 3 x 104 0 0.01 0.02 0.03 0.04 0.05 Real Imag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 Figure 5.11: Reconstructed continuous time eigenvalues for the two-area ex- ample of Fig. 5.5. Case II, improved damping and the new damped frequencies (5.34) are fclosed = 4.721kHz fopen = 4.745kHz (5.68) The improvement in damping due to the changes applied to L2 and C2 is clearly seen in both the computed discrete time domain eigenvalues (Fig. 5.10) and in the re-constructed continuous time ones (Fig. 5.11). In this case, after the switching operation, the new location of the eigenvalues associated to L2 and C2 lay inside the unit circle in the discrete time domain, providing the circuit with a better damping ratio compared to the first case analyzed. 109 Chapter 5. Validation Vgen State 1 State 3State 2 R1 L1 C1 C2 R2 L2 SW1 Non lineal Load State 4 Figure 5.12: Two-area circuit test case with nonlinear load 5.4 Two-Area Case with Nonlinear Load A particularly interesting area of study is the identification of transient eigen- value trajectories in systems with poorly damped eigenvalues that are being heavily stressed. These systems can present very critical operational trajec- tories. This situation can be typical of overstressed power systems operating under deregulated electricity markets. To illustrate the behaviour of eigenvalues trajectories of a nonlinear circuit, consider the presence of a nonlinear load in the circuit of Fig. 5.12. Assume the load is defined as follows RLoad = iC2(VC2) η (5.69) with η = 0.18. This case presents discrete time eigenvalues that move outside the unit circle (unstable region) during the circuit’s energization (closing switch SW1 at t = 0). The discrete time eigenvalues trajectories for this 110 Chapter 5. Validation −1.5 −1 −0.5 0 0.5 1 1.5 −2 −1 0 1 2 0 0.5 1 RealImag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 −1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 RealImag Figure 5.13: Nonlinear load - discrete time eigenvalues case are plotted in Fig. 5.13, while the continuous time ones are plotted in Fig. 5.15. In the upper plot of Fig. 5.13, it can be observed that the eigenvalues are moving around their ultimate steady state location. Also, due to the damping of the circuit elements, the excursion around the steady state location decreases as the simulation evolves. If some of the eigenvalues located close to the instability region (unit circle) are excited during the transient, their trajectories can lead to an eventual system collapse. In this case static eigenvalue analysis which would show only the final steady state point would not be sufficient to predict the eventual system collapse. Computing the point by point eigenvalues trajectories from the EMTP 111 Chapter 5. Validation −15 −10 −5 0 5 x 1011 −5 0 5 x 104 0 0.5 1 RealImag tim e [s] eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1 −14 −12 −10 −8 −6 −4 −2 0 2 x 1011 −5 0 5 x 104 Real Imag Figure 5.14: Nonlinear load - continuous time reconstructed eigenvalues from discrete time solution formulation with the presented methodology brings out the actual nature of the system and identifies the system components the behaviour is associated with, thus facilitating mitigation and corrective actions. 5.5 Voltage Collapse in a Radial System Figure 5.16 shows a simple radial system. The load is increased until voltage collapse is reached. As pointed out in [3], the Voltage collapse is the process by which the sequence of events accompanying voltage instability leads to 112 Chapter 5. Validation 0 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 1 − 505 x 10 5 Vo lta ge N od e 1 tim e [s] voltage [V] 0 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 1 − 101 x 10 6 Vo lta ge N od e 2 tim e [s] voltage [V] 0 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 1 − 101 x 10 6 Vo lta ge N od e 3 tim e [s] voltage [V] 0 0. 1 0. 2 0. 3 0. 4 0. 5 0. 6 0. 7 0. 8 0. 9 1 − 101 x 10 6 Vo lta ge N od e 4 tim e [s] voltage [V] Figure 5.15: Two-area system with nonlinear load 113 Chapter 5. Validation Vgen Zline V1 V2 P + jQ P cteExp. load Figure 5.16: A simple radial system a low unacceptable voltage profile in a significant part of the power system. Several indices to predict voltage collapse are available [61]. One of them is the state based “voltage drops” index, which relates the no-load voltage with the on-load voltage. The voltage drops index is given by V V0 = ej φgrid−φload 2 2 cos( φgrid−φload 2 ) (5.70) Where V is the voltage in on-load situation and V0 is the voltage at the same point in no-load situation with φgrid as the representative power angle of the grid impedances and φload the representative power angle of the loads. Typically threshold values for the | V V0 | index is 0.6 (φgrid = 85◦, cos(φload) = 0.95). Figure 5.17 shows the voltage at the load as it collapses below acceptable 114 Chapter 5. Validation operation values. Five characteristic transition points are identified for the analysis, A1, A2, B, C and D in which the system voltage drops. For this example, an exponential model was used to represent the load, P = P0 ( V V0 )α (5.71) Q = Q0 ( V V0 )β (5.72) With α = 0.1, β = 0.6. These coefficients are typical of an industrial load aggregate with a large induction motors component. A bank of capacitors C for power factor correction is connected in parallel with the load. Equations (5.71) and (5.72) are RMS values over one cycle of the signal. The amount of load is increased as indicated in the table below, Time Interval Pt Qt [s] p.u. p.u. [0-2] 0.5 0.2 [2-3] 0.7 0.3 [3-4] 1.8 0.1 [4-5] 1.8 0.12 [5-6] 1.9 0.1 [6-12] 1.55 0.2 The continuous time eigenvalue trajectories of the radial system of Fig. 5.16 are computed from the discrete time domain with the presented method- ology and shown in Fig. 5.18. As can be observed, there is a distinctive change of the system eigenvalues location in the vicinity of the progressive 115 Chapter 5. Validation 1 2 3 4 5 6 7 8 9 10 11 12 0 0.2 0.4 0.6 0.8 1 1.1 time [s] Vo lta ge B us 2 (p u) A1 A2 B C D Figure 5.17: Voltage collapse of a simple radial system 116 Chapter 5. Validation −1 −0.8 −0.6 −0.4 −0.2 0 −2 −1 0 1 2 0 2 4 6 8 10 12 tim e [s] Real Imag A1 B C D A2 Figure 5.18: Reconstructed continuous time eigenvalue trajectory over time for a simple radial system with voltage collapse voltage drops at node v2. Points A1, A2, B, C and D in Fig. 5.17 can be clearly mapped to the eigenvalue trajectories shown in Figs. 5.18 to 5.22. The projection of the continuous time eigenvalues over the time line in the complex plane in Fig. 5.19 shows the value of the eigenvalues at the critical points of Fig. 5.17. This simple representation of specific points is not always enough to understand the system dynamics due to overlapping of eigenvalues at different time instants. It does provide, nonetheless, a visual clue to ex- 117 Chapter 5. Validation −45 −40 −35 −30 −25 −20 −15 −10 −5 0 5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Real Im ag A1 C B D C D B A2 A2 A1 Figure 5.19: Complex plane projection of continuous time eigenvalues over time line of simple radial system of Fig. 5.16 with voltage collapse treme locations of eigenvalues. Figure 5.20 shows the complete trajectory of the discrete time eigenvalues computed from the nodal EMTP equations. The zoomed in area showed in Fig. 5.21 highlights the fact that by analyzing the eigenvalues trajectories it is possible to anticipate the voltage collapse (i.e., by observing the speed of change and damping ratio of eigenvalues). Both the continuous and discrete time eigenvalues trajectories are the consequence of the system’s search for a new equilibrium operation point. Since the load 118 Chapter 5. Validation 0.95 0.96 0.97 0.98 0.99 1 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x 10−3 0 2 4 6 8 10 12 RealImag tim e [s] A1 A2 B C D Figure 5.20: Discrete time eigenvalue trajectory from EMTP solution for the simple radial system of Fig. 5.16 with voltage collapse demand is forcing the system to its limit and the system is not capable of satisfying it, the voltage is dropped and the system is able to find at least a temporal more stable operation point. This new point, (A1 or A2, or B, etc) closer to the origin of the complex plane in the discrete domain, provides us with additional information about the type of transition of the system, a very fast time constant associated to the change with higher damping than the original eigenvalue locations. Immediately after the system suffers the 119 Chapter 5. Validation 0.9992 0.9994 0.9996 0.9998 1 0 0.5 1 1.5 x 10−3 3 4 5 6 7 8 9 10 11 12 RealImag [s] A1 A2 B C D Figure 5.21: Zoom in of discrete time trajectory of Fig. 5.20 voltage drop (A1, A2, B, C and D), the system eigenvalues return to the initial locations unless a new voltage drop is suffered, in which case a new transition is observed. The method proposed in this work can provide sensitive information about proximity of system voltage collapse with sufficient anticipation time to ini- tiate corrective actions. For example, Fig. 5.22 shows the trajectory of the magnitude of eigenvalue 1, |eig1| over time. This trajectory provides a good indication of the proximity of a voltage drop. It can be noticed for instance 120 Chapter 5. Validation 0 2 4 6 8 10 12 0.9988 0.999 0.9992 0.9994 0.9996 0.9998 1 time [s] |ei g 1 | A1 A2 B C D Figure 5.22: Discrete time eigenvalue 1 magnitude, |eig1|, over time that eigenvalue 1 starts to show evident signs of voltage collapse in the prox- imity of 250ms, before the actual drop to point A1, 300ms before point B, 200ms before point C and 450ms before the point of interest, D. This an- ticipation is enough to initiate automatic corrective actions as, for example, line automatic reclosing, or load shedding. 121 Chapter 5. Validation 5.6 11 Bus Two Area System Figure 5.23 shows a Two Area System of 11 buses described in [28]. The system represents two areas with two generators in each area. The connec- tion between areas is achieved by several parallel transmission lines between buses 7, 8 and 9. In steady state operation the area composed by generators 1 and 2 is exporting 400 MW to the other area composed by generators 3 and 4. The system is simulated in the discrete time domain using the RTTP sim- ulator [16], [17], [54] in which the proposed methodology was implemented. The simulation represents the tap changing transformers and it does not take into account the generator controllers. The dynamic behaviour of the system until it reaches the steady state operation can be seen in Fig. 5.24. The system is stressed in steps to obtain the PV curves. Figure 5.25 shows the PV curve at bus 9 with two characteristic points A and B, normal initial operation and collapse point due to load stress, respectively. Using the dis- crete time state space formulation from the EMTP solution presented in this work, the discrete time eigenvalues of the circuit in Fig. 5.23 are computed. Figure 5.26 shows the discrete time eigenvalues of the system for the char- acteristic point A while Figure 5.27 illustrates the discrete time eigenvalues for the characteristic point B. Similarly to the case of the radial system presented in the previous section, it can be observed that the eigenvalue locations change in the vicinity of a voltage collapse. 122 Chapter 5. Validation Sw V: 1. 07 5 Ph : 9. 64 V: 0. 87 2 Ph :-2 6. 36 V: 1. 03 0 Ph :2 0. 20 V: 1. 01 0 Ph : 6. 92 V: 0. 95 4 Ph :-1 2. 57 V: 1. 03 0 Ph :-8 3. 94 Se tti ng s V: 0. 98 6 Ph :-7 3. 61 V: 1. 02 9 Ph :-7 8. 40 V: 1. 00 0 Ph :- 4. 39 V: 1. 08 8 Ph :-8 1. 74 V: 1. 01 0 Ph :-8 3. 94 Bu s_ 3 G 1 Vo lta ge =1 .0 3 Ph as e= 20 .2 G 2 V= 1. 01 P= 10 .5 G 3 V= 1. 03 P= -6 .8 T1 X= 0. 01 67 Ta p= 1. 12 91 T2 X= 0. 01 67 Ta p= 1. 12 91 T3 X= 0. 01 67 Ta p= 1. 10 82 T5 X= 0. 01 67 Ta p= 1. 08 10 k m R =0 .0 02 0 X= 0. 02 0 B= 0. 01 70 11 0 km R =0 .0 07 5 X= 0. 07 5 B= 0. 12 83 R =0 .0 07 5 X= 0. 07 5 B= 0. 12 83 Q= -3 .5 St ea dy S ta te O rd er =7 St ea dy S ta te F re qu en cy =6 0. 0 M ax S te ps =1 00 00 St op in C ol la ps e= y Co lla ps e Vo lta ge =0 .3 Ig no re T ol er an ce =n o Co m pu te R M S= 1 To le ra nc e= 1. 0e -5 M VA R B as e= 10 0 25 k m R =0 .0 02 5 X= 0. 02 5 B= 0. 04 37 5 R =0 .0 01 0 X= 0. 01 0 B= 0. 01 75 Q =- 2. 0 PQ [7] P= 9. 67 Q= 1 11 0 km R =0 .0 18 X= 0. 18 B= 0. 19 R =0 .0 18 X= 0. 18 B= 0. 19 25 10 k m R =0 .0 01 X= 0. 01 B= 0. 01 75 PQ [9] P= 1. 76 7 Q= 1 25 k m R =0 .0 02 5 X= 0. 02 5 B= 0. 04 37 G 4 V= 1. 01 P= -1 7. 0 R =0 .0 07 5 X= 0. 07 5 B= 0. 12 83 R =0 .0 07 5 X= 0. 07 5 B= 0. 12 83 R =0 .0 18 X= 0. 18 B= 0. 19 25 Bu s_ 1 Bu s_ 5 Bu s_ 2 Bu s_ 7 Bu s_ 11 Bu s_ 6 Bu s_ 8 Bu s_ 9 Bu s_ 10 Bu s_ 4 Figure 5.23: Two-area 11 bus test system 123 Chapter 5. Validation 0 5 10 15 20 25 30 0. 850. 9 0. 951 1. 051. 1 1. 151. 2 1. 25 Ti m e (s) Bus Voltage Magnitude (pu) Bu s1 Bu s2 Bu s3 Bu s4 Bu s5 Bu s6 Bu s7 Bu s8 Bu s9 Bu s1 0 Bu s1 1 Figure 5.24: Two-area 11 bus test system - Bus Voltages 124 Chapter 5. Validation 0 20 0 40 0 60 0 80 0 10 00 12 00 14 00 16 00 18 00 20 00 0 0. 2 0. 4 0. 6 0. 81 (B US −9 ) pf =0 .8 70 P= 17 6. 70 0, V =0 .9 86 P= 15 83 .7 83 , V =0 .5 66 PL (M W ) V(PU) A B Figure 5.25: Two-area 11 bus test system - PV curve at Bus 9 125 Chapter 5. Validation 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 Im ag Real Figure 5.26: Two-area 11 bus test system - Discrete time Eigenvalues normal operation point A 126 Chapter 5. Validation 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 Im ag Real Figure 5.27: Two-area 11 bus test system - Discrete time Eigenvalues at collapse point B 127 Chapter 6 Conclusions and Recommendations 6.1 Conclusions The main contribution of this Thesis is the construction of the discrete time transition matrix as a function of the branch parameters and connectivity. To the best of the author’s knowledge this has never been previously proposed. The present work expands the capabilities of leading transient analysis pro- grams for power systems analysis based on the EMTP solution to produce not only the time domain waveforms of voltages and currents in the sys- tem, but also to trace the system’s dynamic and stability behaviour. This is achieved by extracting the discrete time eigenvalue information from the system difference equations at each time step of the solution. This infor- mation can be used directly in the discrete time domain by observing the stability margins around the unit circle, or by re-mapping the discrete time eigenvalues into the continuous time eigenvalues of the original physical sys- tem. This backwards mapping is exact because the distortion introduced by 128 Chapter 6. Conclusions and Recommendations the discretization rule in going from continuous to discrete time is applied in reverse when going from discrete to continuous time. Notice that this is a specific property of the eigenvalues and it does not apply to voltages and currents in the network. By using the EMTP as the base solution, system nonlinearities (as for ex- ample the nonlinear nature of the load), which have a fundamental effect on the system dynamics, can be easily taken into account in the calculation of the eigenvalues because they are part of the EMTP network representation at each time step of the system’s response. The ability of transients programs, like the EMTP, Microtran, or OVNI, together with the extension presented in this thesis to incorporate eigen- value/eigenvector solutions at each point of the system trajectory will allow for a much more accurate assessment and corrective control measures under critical system operating conditions like those leading to extensive blackouts. System collapse during these situations is driven by a complex interaction of fast and slow system dynamics. Conventional quasi-steady state solutions, like those of traditional transient stability programs, do not have the capabil- ity of considering these interactions using computational tools that are fast enough for on-line system operation. In the context of the OVNI simulator, an important application of the tech- nique proposed in this thesis is the identification of network areas that can be solved with large time steps and areas that require small time steps. As shown, the necessary information can be obtained from the participation ma- 129 Chapter 6. Conclusions and Recommendations trix [P ] computed directly from the nodal description. The exploitation of this system latency is a key aspect in the ability of OVNI to represent very large systems in real-time using network partitioning techniques. This dissertation has presented the theory of the proposed eigenvalue/eigenvector analysis technique, as well as illustrative examples of its ability to detect system conditions leading to instability due to the dynamic behaviour of the system eigenvalues under situations of nonlinear loads. This type of analy- sis was not available before within the EMTP. Traditionally, power system stability information could only be obtained for the system’s steady state operating points, thus ignoring the fact that under situations of high stress, like those leading to voltage collapse and extensive blackouts, the system’s trajectory may no lead to stable steady state points. With the new presented approach, it is possible, with the appropriate control strategies, to redirect the system’s trajectory before the system becomes unstable and collapses. Current work to be reported in future publications includes the combination of this technique with the multilevel MATE concept of [57] which achieves a very efficient representation of control systems corrective actions in the OVNI simulator for real-time on-line applications. The ability to trace dynamically the trajectory of the eigenvalues of nonlinear systems is of particular importance for the following studies: i. design and operation of embedded adaptive controllers ii. automatic selection of discretization steps 130 Chapter 6. Conclusions and Recommendations iii. determination of suitable strategies to avert blackouts in power systems working close to their operational limits iv. dynamic small-signal stability analysis with the EMTP An additional advantage of the discrete time eigenvalue trajectory analy- sis when compared to traditional power systems stability tools is its ability to keep the physical network mapped to the eigenvalue/eigenvector analysis layer, thus allowing for a fast and accurate identification of elements con- tributing to instability. 6.2 Recommendations for Future Work The 21st century power system will make extensive use of distributed intel- ligent local processing nodes capable of making their own decisions and able to interact with other related nodes and with centralized regional controls. Since all parts of a power network interact with each other, local control action may require extensive information regarding the external system to which the node is connected. This information can be acquired and transmit- ted to all the distributed nodes or, even more efficiently, it can be simulated in real-time at each decision node. The use of MATE concepts together with eigenvalue/eigenvector trajectory analysis will allow for efficient embedded control devices. 131 Chapter 6. Conclusions and Recommendations One area of particular interest is that of distributed generation with renew- able energy sources such as wind farms. These generation clusters are highly nonlinear and require advance on-site control devices to allow a healthy in- terconnection with the power grid. The proposed methodology together with embedded simulation capabilities can provide a solution to this problem. Current research work [62],[63] within the Joint Infrastructure Interdepen- dencies Research Program (JIIRP) (sponsored by Public Safety and Emer- gency Preparedness Canada (PSEPC) and the Natural Sciences and Engi- neering Research Council of Canada (NSERC)) is concerned with the iden- tification of critical interdependencies among multiple infrastructures (e.g., power grid, water grid, etc.) and the dynamic evolution of these interde- pendencies over the time line. The dynamic behaviour of interdependencies among critical infrastructures can be mapped into a multiple-infrastructures transportation matrix which can then be transformed into a state space equation from which eigenvalues and trajectories can be obtained with the methodology presented in this work. The following study areas are suggested for future work: i. Development of distributed intelligent control solutions based on embed- ded OVNI and discrete state space eigenvalue trajectory computation. ii. Integration of discrete state space eigenvalue analysis methodology with the Latency concept. iii. Development of more detailed Induction motor models. 132 Chapter 6. Conclusions and Recommendations iv. Implementation of the discrete state space eigenvalue analysis method- ology within UBC’s OVNI-NET simulator for stability analysis and de- termination of segmentation schemes. v. Implementation of discrete state space eigenvalue analysis methodology to UBC-JIIRP’s I2Sim simulator for identification of trajectories of crit- ical interdependencies among Infrastructures. 133 Bibliography [1] S. Gissinger, P. Chaunmes, J.P. Antoine, A. Bihain, and M. Stubbe, “Advanced dispatcher training simulator,” Computer Applications in Power, IEEE, vol. 13, no. 2, pp. 25–30, 2000. [2] M.V.F. Pereira, M. McCoy, and H.M. Merril, “Managing risk in the new power business,” Computer Applications in Power, IEEE, vol. 13, no. 2, pp. 18–24, 2000. [3] P. Kundur, Power System Stability and Control, The EPRI power system engineering series. McGraw-Hill, 1994. [4] C.A.Da Neto, D.M. Falcao, and J.L. Jardim, “A unified online security assesment system,” in CIGRE session 2000. 2000, pp. 38–102, CIGRE. [5] H.W. Dommel, EMTP Theory Book, Microtran Power System Analysis Corporation, Vancouver, 2nd edition, 1992. [6] H.W. Dommel, “Digital computer solution of electromagnetic transient in single and multiphase networks,” IEEE Transactions on Power Ap- paratus and Systems, vol. PAS-88, no. 4, pp. 388–399, 1969. 134 Bibliography [7] G.D. Irwin, R. Kuffel, and D.AWoodford, “EMTDC- high performance electromagnetic transient simulation,” Electro Soft Journal, 1981. [8] G. Sybille and et al., “Theory and application of power system blockset, a matlab/simulink-based simulation tools for power systems,” in IEEE PES winter meeting. 2000, IEEE. [9] P.G. McLaren, R. Kuffel, R. Wierckx, J. Giesbrecht, and L. Arendt, “A real time digital simulator for testing relays,” IEEE Transactions on Power Delivery, vol. 7, pp. 207–213, 1992. [10] J.R. Mart́ı and L.R. Linares, “Real-time EMTP-based transients sim- ulation,” Power Systems, IEEE Transactions on, vol. 9, no. 3, pp. 1309–1317, 1994. [11] M. Kesunovic, M. Aganic, V. Skendzic, J. Domaszewicz, J.K. Bladow, D.M. Hamai, and S.M. McKenna, “Transients computation for relay testing in real-time,” IEEE Transactions on Power Delivery, vol. 9, pp. 1298–1307, 1994. [12] Y. Fuyimoto, Y. Bin, H. Taoka, H. Tezuka, S. Sumimoto, and Y. Ishika, “Real-time power system simulator on a pc cluster,” in International Conference on Power System Technology, Budapest, 1999, IPST’99. [13] J.R. Mart́ı, L.R. Linares, J. Calvino, H.W. Dommel, and J. Lin, “OVNI: an object approach to real-time power system simulators,” in Power 135 Bibliography System Technology, 1998. Proceedings. POWERCON ’98. 1998 Inter- national Conference on, 1998, vol. 2, pp. 977–981 vol.2. [14] Fromonto, Levacher, and Strunz, “An efficient parallel computer archi- tecture for real time simulations of electromagnetic transients in power systems,” in 31st North American Power Sysmposium, 1999, pp. 490– 496. [15] L.R. Linares, A Real time simulator for power electric networks, M.A.Sc. Thesis. University of British Columbia, Vancouver, 1993. [16] Jesús Calviño Fraga, Implementation of a real time simulator for pro- tective relay testing, M.A.Sc. Thesis. University of British Columbia, Vancouver, 1999. [17] J. A. Hollman, Real time distributed network simulation with pc clusters, M.A.Sc. Thesis. University of British Columbia, Vancouver, 1999. [18] J.R. Mart́ı, L.R. Linares, Roberto Rosales, and H.W. Dommel, “OVNI: a full-size real time power system simulator,” in International Confer- ence on Digital Power System Simulators, Montreal, Quebec, Canada, 1997, pp. 1309–1317, ICDS’97. [19] L.R. Linares, OVNI (Object Virtual Network Integrator) a new fast algorithm for the simulation of very large electric networks in real time, Ph.D. Thesis. University of British Columbia, Vancouver, 2000. 136 Bibliography [20] F. Moreira, J.A. Hollman, L.R. Linares, and J. R. Mart́ı, “Network decoupling by latency exploitation and distributed hardware architec- ture,” in International Conference on Power System Transients, Rio de Janeiro, 2001, IPST’01. [21] J.R. Mart́ı, L.R. Linares, and J.A. Hollman, “OVNI: Integrated soft- ware/hardware solution for real-time simulation of large power sys- tems,” in 14th Power Systems Computer Conference, Seville, Spain, 2002, PSCC’02. [22] T. De Rybel, J.A. Hollman, and J. Mart́ı, “OVNI-NET: A flexible cluster interconnect for the new OVNI real-time simulator,” in 15th Power Systems Computer Conference, Liege, Belgium, 2005, PSCC’05. [23] B. Gao, G.K. Morison, and P. Kundur, “Towards the development of a systematic approach for voltage stability assessment of large-scale power systems,” Power Systems, IEEE Transactions on, vol. 11, no. 3, pp. 1314–1324, 1996, TY - JOUR. [24] P. Kundur and G.K. Morison, “Developments in power system analysis,” IEEE Canadian Review, pp. 5–9, 1995. [25] Powertech Labs Inc., “RAP: Reliability assesment package,” 1998. [26] Powertech Labs Inc., “PSAPAC: Power system analysis package,” 1998. [27] Powertech Labs Inc., “VSAT: Voltage security assesment tool,” 2002. 137 Bibliography [28] Powertech Labs Inc., “TSAT: Transient security assessment tool,” 2002. [29] K.S. Shim, H.K. Nam, S.G. Song, Y.G. Kim, and K.Y. Lee, “Applica- tion results of the eigen-sensitivity theory of augmented matrix to small signal stability analysis of large power systems,” in Power Engineering Society Summer Meeting, 2000. IEEE, 2000, vol. 3, pp. 1835–1838 vol. 3. [30] C.G. Groom, K.W. Chan, R.W. Dunn, and A.R. Daniels, “Real-time security assessment of electrical power systems,” Power Systems, IEEE Transactions on, vol. 11, no. 2, pp. 1112–1117, 1996, TY - JOUR. [31] J.L. Jardim, “Online dynamic security assessment: implementation problems and potential use of artificial intelligence,” in Power Engi- neering Society Summer Meeting, 2000. IEEE, 2000, vol. 1, pp. 340–345 vol. 1, TY - CONF. [32] W.E. Arnoldi, “The principle of minimized iterations in the solution of the matrix eigenvalue problem,” Quart. Appl. Math., vol. 9, pp. 17–29, 1951. [33] R. Kumaresan, D.W.; Tufts, and L.L Scharf, “A prony method for noisy data: Choosing the signal components and selecting the order in exponential signal models,” Proceedings of the IEEE, vol. 72, pp. 230–233, 1984. 138 Bibliography [34] M.A. Johnson, I.P. Zarafonitis, and M. Calligaris, “Prony analysis and power system stability-some recent theoretical and applications re- search,” in Power Engineering Society Summer meeting, 2000, vol. 3, pp. 1928–1923. [35] Baron (Gaspard Riche) de Prony, “Essai experimental et analytique,” J. E. Polytech., vol. 1, no. 2, pp. 24–76, 1795. [36] R. Kumaresan and D. W. Tufts, “Estimation the parameters of ex- ponentially damped sinusoids and pole-zero modeling in noise,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. ASSP-30, pp. 833–840, 1982. [37] B. Porat, A Course in Digital Signal Processing, John Wiley and Sons, 1st edition edition, 1997. [38] C.E. Grund, J.J. Paserba, J.F. Hauer, and S.L Nilsson, “Comparison of prony and eigenanalysis for power system control design,” Power Systems, IEEE Transactions on, vol. 8, no. 3, pp. 964–971, 1993. [39] P. Kundur, G.K. Morison, L. Wang, and H. Hamadanizadeh, “On-line dynamic security assessment of power systems,” in Fifth International Workshop on Electrical Power Control Centres, Hungary, 1999. [40] Y. Saad, “Variations of arnoldi’s method for computing eigenelements of large unsymmetric matrices,” Linear Algebra and its applications, vol. 34, pp. 269–295, 1980. 139 Bibliography [41] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press, Oxford, 1965. [42] D.T. Wong, B. Rogers, G. andPorreta, and P. P. Kundur, “Eigenvalue analysis of very large power systems,” IEEE Trans. on Power Systems, vol. PWRS-3, 1998. [43] B.K. Perkins, J.R. Mart́ı, and H.W. Dommel, “Nonlinear elements in the EMTP: steady-state initialization,” IEEE Transactions on Power Systems, vol. 10, no. 2, pp. 593–601, 1995. [44] M. Kezunovic, L. Kojovic, A. Abur, C.W. Fromen, D.R. Sevcik, and F. Phillips, “Experimental evaluation of EMTP-based current trans- former models for protective relay transient study,” IEEE Transactions on Power Delivery, vol. 9, no. 1, pp. 405–413, 1994. [45] B. Khodabakhchian and G.T. Vuong, “Modelling a mized residential- commercial load for simulations involving large disturbances,” IEEE Transactions on Power Systems, vol. 12, no. 2, pp. 791–796, 1997. [46] N. Nagaoka, S. Tawada, and A. Ametani, “Magnetizing circuit model for EMTP including voltage dependence characteristic,” in Transmission and Distribution Conference and Exhibition 2002, Asia Pacific, 2002, vol. 2, pp. 1265–1269. 140 Bibliography [47] W. Long, D. Cotcher, D. Ruiu, P. Adam, S. Lee, and R. Adapa, “EMTP- a powerful tool for analyzing power system transients,” Com- puter Applications in Power, IEEE, vol. 3, no. 3, pp. 36–41, 1990. [48] J.A. Martinez, “How to adapt the EMTP for classroom instruction,” Power Systems, IEEE Transactions on, vol. 7, no. 1, pp. 351–358, 1992. [49] J.A. Martinez, “Educational use of EMTP models for the study of rotating machine transients,” Power Systems, IEEE Transactions on, vol. 8, no. 4, pp. 1392–1399, 1993. [50] Kim Chul-Hwan, Lee Myung-Hee, R.K. Aggarwal, and A.T. Johns, “Ed- ucational use of EMTP models for the study of rotating machine tran- sients,” Power Systems, IEEE Transactions on, vol. 15, no. 1, pp. 9–15, 2000. [51] J.R. Mart́ı and J. Lin, “Supression of numerical oscillation in the EMTP,” Power Systems, IEEE Transactions on, vol. 4, no. 2, pp. 739–747, 1989. [52] J.A. Hollman and J.R. Mart́ı, “Real time network simulation with pc- cluster,” Power Systems, IEEE Transactions on, vol. 18, no. 2, pp. 563–569, 2003. [53] J.R. Mart́ı, “EECE 560 class notes,” Department of Electrical and Computer Engineering, The University of British Columbia, Vancouver. 141 Bibliography [54] Jesús Calviño Fraga, Tunable Difference Equations for the Time- Domain Simulation of Power System Operating States, Ph.D. Thesis. University of British Columbia, Vancouver, 2003. [55] L.R. Linares and J.R. Mart́ı, “Sub-area latency in a real time power net- work simulator,” in International Conference on Power System Tran- sients, Lisbon, 1995, pp. 541–545, IPST’95. [56] Fernando Moreira, Latency Techniques in Power System Transients Simulation, Ph.D. thesis, University of British Columbia, 2002. [57] M. Armstrong, J.R. Mart́ı, L.R. Linares, and P. Kundur, “Multilevel MATE for efficient simultaneous solution of control systems and non- linearities in the OVNI simulator,” Power Systems, IEEE Transactions on, vol. 21, no. 3, 2006. [58] J.R. Mart́ı, “MATE, a multi-area thev́enin equivalent concept,” In- ternal memorandun of the power systems group, UBC - Department of Electrical and Computer Engineering, August 1993. [59] A. Brameller, Y. M. Hamam, and Ronald Norman Allan, Sparsity : its practical application to systems analysis, Pitman, London ; New York, 1976. [60] A. Brameller, Maldwyn Noel John, and M. R. Scott, Practical diakoptics for electrical networks, Advanced engineering textbooks. Chapman and Hall, London,, 1969. 142 Bibliography [61] Working group 38.02, “Indices predicting voltage collapse including dynamic phenomena,” Task force 11, December 1994. [62] J.R. Mart́ı, J.A. Hollman, C. Ventura, and J. Jatskevich, “Design for survival: Real-time infrastructure coordination,” in International Work- shop on Complex Network and Infrastructure Protection, Rome, 2006, CNIP. [63] J.A. Hollman, J.R. Mart́ı, J. Jatskevich, and K.D. Srivastava, “Dynamic islanding of critical infrastructure a suitable strategy to survive and mit- igate critical events,” in International Workshop on Complex Network and Infrastructure Protection, Rome, 2006, CNIP. [64] Gene Golub and Charles Van Loan, Matrix Computation, Johns Hop- kins Series in the Mathematical Sciences. The Johns Hopkins University Press, Baltimore, Maryland 21211, third edition, 1989. 143 Appendix A The QR Method The QR method allows the calculation of eigenvalues of real matrices. The Stability of a numerical eigenvalue problem depends on the matrix under consideration. If the matrix is symmetric with symmetrically distributed error, then the calculated eigenvalues will approximate the actual eigenvalues, provided the eigenvalues are all simple. Otherwise, the numerical methods my fail to find all eigenvalues. The basis of the QR method for calculating the eigenvalues of A is the fact that an n x n real matrix can be written as A = QR (A.1) where Q is orthogonal and R is upper triangular. The method is efficient for the calculation of all eigenvalues of a matrix. The construction of Q and R proceeds as follows. Matrices P1, P2, ...Pn−1 are constructed so that Pn−1Pn−2....P2P1A = R is upper triangular. These matrices can be chosen as orthogonal matrices and are called Householder matrices. Since the P ′s are orthogonal, the stability of the eigenvalue problem will not be worsened. 144 Appendix A. The QR Method If we let QT = Pn−1Pn−2...P2P1 (A.2) then we have QTA = R and QQTA = QR IA = QR A = QR The method can be summarized with the following steps, 1. Set A1 = A, Q1 = Q, and R1 = R 2. First set A2 = R1Q1; then factor A2 as A2 = Q2R2. 3. First set A3 = R2Q2; then factor A3 as A3 = Q3R3. 4. Set Am = Rm−1Qm−1; then factor Am−1 as Am = QmRm. Matrix Am will tend toward a triangular or nearly triangular form. Thus the eigenvalues of Am will be easy to compute. If the eigenvalues can be ordered as |λ1| > |λ2| > |λ3| > · · · > |λn| > 0, then as m increases the eigenvalues of Am approach the eigenvalues of A. The QR method can be used for an arbitrary real matrix, but works much faster on special matrices: • symmetric tri-diagonal, 145 Appendix A. The QR Method • Hessenberg, • symmetric band matrices. It is important to note that if the elements of A have widely varying magni- tudes, then A should be balanced before applying the QR algorithm. When A is balanced, the computed eigenvalues are often more accurate [64]. 146 Index Argentina, 1 Arnoldi generalized method, 17 modified method, 26 Bilinear Transformation, 63 Brazil, 1 Calviño, Jesús, xiv, 7 CDA, 31, 39 Chile, 1 continues time state space, 82 continuous time eigenvalues, 90 discrete eigenvalues, 89, 90 discrete time state space, 4, 40 capacitance, 46 inductance, 44 MATE, 73 parallel branches, 55 resistance, 48 rlc series, 85 series branches, 50 trapezoidal, 44, 46 Dommel, Hermann, xiv, 4, 5, 42 Edison, Thomas Alva, 1 Eigenvalue, 16, 55 Electricite de France, 6 EMTP, ii, 5 fast time-domain, 17 Gaspard Richie, Baron de Prony, 17 Gissinger, S., 3 Helmholtz, 44 Hollman, Jorge, 7 HVDC, 10 Hydro Quebec Research Institute, 6 147 Index I2Sim, 133 interarea oscillation, 15, 24 Jardim, J.L., 4 JIIRP, 132, 133 Kundur, Prabha, xiv, 3 Linares, Luis, xiv, 7 LTI, 40 Mart́ı, José, xiv, 7 Mitsubishi, 6 Modified Arnoldi, 16 Moshref, Ali, xiv NSERC, 132 OVNI, ii, 4, 7, 67 MATE, 9 OVNI-NET, 9 PC-cluster, ii Pereira, M., 3 Prony, 16 PSEPC, 132 QR method, 144 RTDNS, 7 RTNS, 7 small-signal stability, 15 TEQSIM, 6 Thevenin, 44 Transient Network Analyzer, 5 transition matrix, 89 trapezoidal, 90 148
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Step by step eigenvalue analysis with EMTP discrete...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Step by step eigenvalue analysis with EMTP discrete time solutions Hollman, Jorge 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Step by step eigenvalue analysis with EMTP discrete time solutions |
Creator |
Hollman, Jorge |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | The present work introduces a methodology to obtain a discrete time state space representation of an electrical network using the nodal [G] matrix of the Electromagnetic Transients Program (EMTP) solution. This is the first time the connection between the EMTP nodal analysis solution and a corresponding state-space formulation is presented. Compared to conventional state space solutions, the nodal EMTP solution is computationally much more efficient. Compared to the phasor solutions used in transient stability analysis, the proposed approach captures a much wider range of eigenvalues and system operating states. A fundamental advantage of extracting the system eigenvalues directly from the EMTP solution is the ability of the EMTP to follow the characteristics of nonlinearities. The system's trajectory can be accurately traced and the calculated eigenvalues and eigenvectors correctly represent the system's instantaneous dynamics. In addition, the algorithm can be used as a tool to identify network partitioning subsystems suitable for real-time hybrid power system simulator environments, including the implementation of multi-time scale solutions. The proposed technique can be implemented as an extension to any EMTP-based simulator. Within our UBC research group, it is aimed at extending the capabilities of our real-time PC-cluster Object Virtual Network Integrator (OVNI) simulator. |
Extent | 2766613 bytes |
Subject |
EMTP Power systems Real time simulation Eigenvalue analysis Discrete time state space OVNI PC-cluster Latency Multi-time scale solutions |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2006-10-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0065505 |
URI | http://hdl.handle.net/2429/67 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2006_fall_hollman_jorge.pdf [ 2.64MB ]
- Metadata
- JSON: 24-1.0065505.json
- JSON-LD: 24-1.0065505-ld.json
- RDF/XML (Pretty): 24-1.0065505-rdf.xml
- RDF/JSON: 24-1.0065505-rdf.json
- Turtle: 24-1.0065505-turtle.txt
- N-Triples: 24-1.0065505-rdf-ntriples.txt
- Original Record: 24-1.0065505-source.json
- Full Text
- 24-1.0065505-fulltext.txt
- Citation
- 24-1.0065505.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0065505/manifest