UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Step by step eigenvalue analysis with EMTP discrete time solutions Hollman, Jorge 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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 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 EMTPbased 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.  ii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables  vi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Notation  x  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii  Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Dedication  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv  1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . 1.2 Research Background . . . . . . . . . . . . . . . . . 1.2.1 Off-Line Power System Simulation . . . . . 1.2.2 Real-Time Power System Simulation . . . . 1.2.3 Large Scale Power System Stability Analysis 1.3 Research Claim and Contributions . . . . . . . . . 1.4 List of Publications . . . . . . . . . . . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . 1 . 1 . 5 . 5 . 6 . 9 . 11 . 12  2 Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The Problem . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Traditional Analysis Tools . . . . . . . . . . . . 2.2.2 Advantages of Step by Step Trajectory Analysis  . . . . .  . . . . .  . . . . .  . . . . .  iii  . . . . . . . .  14 14 14 17 27  Table of Contents 3 Eigenvalue Analysis with EMTP Solution . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 EMTP Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Nodal Analysis . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Numerical Integration Accuracy . . . . . . . . . . . . . 3.2.3 Discretization Rule Stability . . . . . . . . . . . . . . . 3.3 Discrete Time State Space Formulation of an Electrical Network from the EMTP Nodal Matrix Formulation . . . . . . . 3.4 EMTP Solution in State Space form . . . . . . . . . . . . . . 3.4.1 Discrete Time State Space Equation for an Inductance 3.4.2 Discrete Time State Space Equation of a Capacitance . 3.4.3 Discrete Time State Space Equation of a Resistance . . 3.4.4 Treatment of Series Connection RL and RC . . . . . . 3.4.5 Treatment of Series Branches . . . . . . . . . . . . . . 3.4.6 Treatment of Parallel Branches . . . . . . . . . . . . . 3.5 Nonlinear Networks . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Continuous Time Eigenvalues from Discrete Time Domain Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Discrete to Continuous Time Mapping . . . . . . . . . 3.6.2 Time Discretization Considerations . . . . . . . . . . . 4 Application to the OVNI simulator . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 OVNI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 The MATE Concept . . . . . . . . . . . . . . . . . . 4.2.2 Discrete State Space Formulation with MATE . . . . 4.3 The Latency Concept . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Application of the New Discrete Time State Space Formulation to Latency . . . . . . . . . . . . . . . . . . 4.4 Automatic Time Step Selection . . . . . . . . . . . . . . . .  . . . . . .  40 42 44 46 48 49 50 55 57 59 60 63 66 66 67 68 73 76  . 77 . 79  5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Comparison of State Space Solutions between Continuous Time and Discrete Time Domains . . . . . . . . . . . . . . . . . . . 5.1.1 RLC Series Circuit State Space Equation from Continuous Time Domain . . . . . . . . . . . . . . . . . . . . 5.1.2 RLC Series Circuit State Space Equation from the [G] EMTP Matrix . . . . . . . . . . . . . . . . . . . . . . iv  29 29 29 32 34 37  81 81 82 85  Table of Contents 5.2 Latency Case . . . . . . . . . . . . . 5.3 Two-Area Resonant Circuit . . . . . 5.3.1 Case I - Poorly Damped . . . 5.3.2 Case II - Improved Damping . 5.4 Two-Area Case with Nonlinear Load 5.5 Voltage Collapse in a Radial System 5.6 11 Bus Two Area System . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  93 99 100 106 110 112 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 2.1 Linear time invariant (LTI) dynamic system . . . . . . . . . 2.2 Two-area 11 bus test system . . . . . . . . . . . . . . . . . . 2.3 Phasor value of bus 8 voltage for the two-area test system, small-signal instability observed with detailed time domain simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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) . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15  Inductor element . . . . . . . . . . . . . . . . . . . . . . . . Inductor discrete time equivalent with trapezoidal rule . . . Magnitude frequency response of integration rules . . . . . . Phase frequency response of integration rules . . . . . . . . . State space representation of a multiple input/output system State space system diagram block . . . . . . . . . . . . . . Discrete time equivalent of an Inductor . . . . . . . . . . . . Discrete equivalent of a Capacitor . . . . . . . . . . . . . . . Discrete equivalent of a Resistor . . . . . . . . . . . . . . . . Discrete equivalent of a Series RL . . . . . . . . . . . . . . . Generic circuit . . . . . . . . . . . . . . . . . . . . . . . . . Aggregate treatment of parallel RLC branches . . . . . . . Piecewise linear inductance with two slopes . . . . . . . . . . Implementation of two slope piecewise linear inductance . . Continuous time and discrete time stability mapping for trapezoidal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Continuous time and discrete time stability mapping for backward Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . .  8  . 18 . 22  . 23  . 24 . . . . . . . . . . . . . .  30 32 35 36 39 43 45 47 49 50 51 56 58 59  . 64 . 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 4.3 General update scheme for discrete state space formulation with MATE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Hybrid real-time/soft real-time simulator layout . . . . . . . . 4.5 Automatic EMTP time step selection scheme . . . . . . . . . . 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19  5.20 5.21  RLC series circuit . . . . . . . . . . . . . . . . . . . . . . . . Discretized RLC series circuit . . . . . . . . . . . . . . . . . Latency test circuit . . . . . . . . . . . . . . . . . . . . . . . Circuit with slow and fast subareas . . . . . . . . . . . . . . Two-area resonant test case . . . . . . . . . . . . . . . . . . Node voltages for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . . . . . . . . . . . . . . . . . Discrete time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . . . . . . . . . . . . . Reconstructed continuous time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped . . . . . . . . . Node voltages for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . . . . . . . . . . . . . . . . . Discrete time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . . . . . . . . . . . . . Reconstructed continuous time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping . . . . . . . Two-area circuit test case with nonlinear load . . . . . . . . Nonlinear load - discrete time eigenvalues . . . . . . . . . . . Nonlinear load - continuous time reconstructed eigenvalues from discrete time solution . . . . . . . . . . . . . . . . . . . Two-area system with nonlinear load . . . . . . . . . . . . . A simple radial system . . . . . . . . . . . . . . . . . . . . . Voltage collapse of a simple radial system . . . . . . . . . . . Reconstructed continuous time eigenvalue trajectory over time for a simple radial system with voltage collapse . . . . . . . Complex plane projection of continuous time eigenvalues over time line of simple radial system of Fig. 5.16 with voltage collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discrete time eigenvalue trajectory from EMTP solution for the simple radial system of Fig. 5.16 with voltage collapse . Zoom in of discrete time trajectory of Fig. 5.20 . . . . . . . viii  . . . . .  69 74 75 80 82 86 93 99 99  . 103 . 104 . 105 . 107 . 108 . 109 . 110 . 111 . . . .  112 113 114 116  . 117  . 118 . 119 . 120  List of Figures 5.22 5.23 5.24 5.25 5.26  Discrete time eigenvalue 1 magnitude, |eig1 |, over time . . . Two-area 11 bus test system . . . . . . . . . . . . . . . . . . Two-area 11 bus test system - Bus Voltages . . . . . . . . . Two-area 11 bus test system - PV curve at Bus 9 . . . . . . Two-area 11 bus test system - Discrete time Eigenvalues normal operation point A . . . . . . . . . . . . . . . . . . . . . 5.27 Two-area 11 bus test system - Discrete time Eigenvalues at collapse point B . . . . . . . . . . . . . . . . . . . . . . . . .  ix  . . . .  121 123 124 125  . 126 . 127  Acronyms ATA Advanced Technology Attachment, a disk drive implementation that integrates the controller on the disk drive itself DSP Digital Signal Processing EDF Electricit´e de France EMTDC Transients Simulation Software by Manitoba HVDC Research Centre EMTP Electromagnetic Transients Program EPRI Electric Power Research Institute FACTS Flexible AC Transmission Systems 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] [GB ] GB43 gB43 [L] [A]d [A]c λ z x(t) ˙ x(t + ∆t) a a vkm h(t)km h(t − ∆t)km  admittance matrix admittance matrix of subnetwork [B] element (4, 3) of admittance matrix [B] admittance of branch nodes 4 and 3 in subnetwork [B] incidence matrix discrete time transition matrix continuous time transition matrix continuous time eigenvalue discrete time eigenvalue future state in continuous time domain future state in discrete time domain compact notation for a(t − ∆t) compact notation for a(t) branch voltage of node k to node m history of branch km at (t) 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)  hi  history of branch i at (t)  hi [k∧ ] [g∧ ] ζ ζc f fnc [P ] [φ] [ϕ] fN y  compact notation for history of branch i at (t − ∆t) diagonal matrix of branch nature coefficients diagonal matrix of branch discretization rule coefficients damping ratio damping ratio from continuous time domain oscillation frequency undamped natural frequency from continuous time domain Participation matrix Right eigenvector Left eigenvector 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 acknowledgement 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, Roc´ı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 Roc´ı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 Mar´ı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 substations 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´onica, and my parents and brother in law Lieschen, Hansi and Andr´es, who always believed in us and unconditionally supported us all these years. I also want to extend my thanks to our friends Ver´onica, Mat´ıas, Lucas, Dante and Sof´ıa Salibi´an, who became our family in Canada and made our daily life a pleasant adventure; and to our good friends Paula Murcia, Pablo Ruiz, Jes´ us Calvi˜ no-Fraga, Ver´onica, Javier and Mat´ıas Gazzarri, Lito and Alicia Caro, Jorma and Laura Neuvonen, Hazem and Samar Sabbagh, Gill and Jorge Palejko, Mar´ıa Isabel and Felipe Gajardo, Naoko and J´er´emie S´eror, 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´ us Calvi˜ no-Fraga, Dr. Luis Linares, Dr. Carlos Ventura, 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 friendship. This dissertation would not have been possible without the generous financial support provided by several sources, including The National Sciences and Engineering Research Council of Canada through an ISP-2 Research Grant; Dr. Jos´e Mart´ı with appointments as Research Assistant; the Department of Electrical and Computer Engineering via part-time appointments as Teaching Assistant; Fundaci´on 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´ecnicas; and Dr. Prabha Kundur and Dr. Ali Moshref from Powertech labs with part-time research and consulting contracts. To all of them goes my appreciation. Finally, I would especially like to express my most sincere gratitude to my mentor and supervisor, Dr. Jos´e Mart´ı. His constant availability, patience, 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, Roc´ı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 operation complexity. Just taking into consideration the case of the United States-Canadian system, for instance, this amazing expansion of the generation 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 deregulated environment. This new operational paradigm does not accept the traditional extra investment in network infrastructure to ensure conservative secure operation margins. The traditional technical constraints together with the newly introduced financial constraints have forced the utilities to 1 Information 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 electricity 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 coordinated 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 capable of simulating external system behaviour in real-time at the local node. Consequently, the concept becomes practical only if this computational solution 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 operation 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 remarkable that many important risks are not measured in dollars”. As well, the relevance of an extensive and detailed representation of the power system 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 introduced 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 operation layout. Given the present development of The University of British Columbia’s (UBC) real-time power system simulator -Object Virtual Network 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 stability analysis of power systems.  1.2 1.2.1  Research Background 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 reduced 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 limitations. 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, neural 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´e de France (EDF) and Hydro Quebec 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 promising for soft real-time applications. In contrast, UBC’s real-time PC cluster architecture shows a constant communication time for each type of configuration 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´ us Calvi˜ no; 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)  MRTNS (Multi-layer-tearing software)  Transient Analysis  SRTNS for Relay  RTDNS  Tester (Single-PC simulator)  (PC-Cluster Simulator)  OVNI  Transient and Stability Analysis Hybrid Simulation Environment Hard or Soft Real-Time  Off-Line Environment  Figure 1.1: Relation among research projects at UBC power system group  real-time power system simulator. By breaking the power network solution 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 solutions available today. However, with the advent of the Multi Area Thevenin Equivalent (MATE) [21] solution scheme, a new, more flexible hardware interconnect 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 OVNINET [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 available [23], [24], [25], [26], [27], [28], [29], a large number of simplifications is a common characteristic among them in order to reduce the computational burden, especially for the case of on-line stability analysis tools. As a 2  Latency 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 realtime performance. Most of the currently available transient analysis software works in the off-line environment. Nonetheless, some groups [30], [31], motivated 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 minutes, 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 dynamics 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 nonlinear load behaviour. Some of the advantages of the new methodology can be summarized as follows: • Trajectory tracking of nonlinear element eigenvalues moment by moment. • Capability of identifying suitable network partitioning schemes for application 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 indexes 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 principal 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 Interconnect for the New OVNI Real-Time Simulator. Proceedings of the 15th Power System Computation Conference. Li´ege, 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 Software/Hardware Solution for Real-Time Simulation of Large Power Systems. 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. Proceed12  Chapter 1. Introduction ings of the International Conference on Power System Transients. Rio de Janeiro, Brazil, 2001. • J. Mart´ı, J. Hollman, J. Calvi˜ no. Implementation of a Real-Time Distributed Network Simulator with PC-Cluster. Proceedings of the International 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 operation 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 between 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 interesting 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 because 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, small15  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 damping ratios of the dominant coupled modes; however, they are usually intensive 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 analysis technique and a system identification method widely used in biomedical monitoring, geophysical analysis, speech processing and power system electromechanical 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 description 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)  u(t)  y(t)  x(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  x(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  x(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, n  (qiT x0 )pi e(λi t)  x(t) ˙ =  (2.3)  i=1  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 described in equation (2.3) by fitting a sum of complex damped sinusoids to evenly spaced time sample values of the output, that is, L  Ai e(σi t) cos(2πfi t + φi )  y(t) =  (2.5)  i=1  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(2πfi t + φi ) =  e(j2πfi t+φi ) + e−(j2πfi t+φi ) e(j2πfi t) ejφi e−(j2πfi t) e−jφi = + 2 2 2 (2.6)  replacing (2.6) into (2.5) and letting t = kT the samples of y(t) are given by L  Ci µki  y[k] = i=1  (2.7)  poles  Where, Ai jφi e 2  (2.8)  µi = e(σi +j2πfi )T  (2.9)  Ci =  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 prediction 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 coefficients 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 techniques 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 result 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 eigenvalues 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 ac21  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 stability 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 domain 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 7 1  5 25 km  8  9  110 km  6  10  110 km 10 km  10 km  11  3  25 km  G3  G1 2  4 L7  L9  G2  G4 400 MW  Area 1  Area 2  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 1.40  Bus voltage magnitude (pu)  1.12  0.84  0.56  0.28  0.00  0.00  3.87  7.74 11.61 Time in seconds  15.48  19.35  Figure 2.3: Phasor value of bus 8 voltage for the two-area test system, smallsignal 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 1.040  Bus voltage magnitude (pu)  1.026  1.012  0.998  0.984  0.970  0.00  6.00  12.00 18.00 Time in seconds  24.00  30.00  Figure 2.4: Phasor value of bus 8 voltage for the two-area test system, smallsignal 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  Node 1 Gen.# 1 Node 2 Gen.# 2 Node 3 Gen.# 3 Node 4 Gen.# 4  Freq Damp % [Hz] 0.62 -2.144 0.97 6.702 1.2 -0.579 0.62 -1.077 0.97 8.926 1.2 1.052 0.62 -2.357 1.03 -1.884 1.7 13.677 0.62 -1.344 1.03 3.442 1.7 -0.675  Real [1/s]  Imag [rad]  0.0831 -0.4103 0.0459 0.0421 -0.5432 0.0421 0.0925 0.1215 -0.1087 0.0527 -0.2335 0.0734  3.875 6.109 7.942 3.910 6.061 6.061 3.924 6.446 10.743 3.921 6.820 10.870  Type of observed mode Interarea Local between G1 Local between G1 Interarea Local between G1 Local between G1 Interarea Local between G3 Local between G3 Interarea Local between G3 Local between G3  & G2 & G2 & G2 & G2 & G4 & G4 & G4 & 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 matrix 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] making the Arnoldi method evolve into one of the most successful ones for nonHermitian 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 order 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)  Where H 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 after m steps, the following relation holds: AVm − Vm Hm = f e∗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 expensive 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 target time of 3 h 08 min 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 4 min 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 4 min02 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 computations 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 studies 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. Consequently, 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  30  diL (t) dt  (3.1)  Chapter 3. Eigenvalue Analysis with EMTP Solution Solving for the voltage,  VL dt = LdiL (t) t  (3.2)  t  VL dt = L  diL (t)  (3.3)  VL dt = L[iL (t) − iL (t − ∆t)]  (3.4)  t−∆t  t−∆t  t t−∆t  where ∆t is the discretization time step. Approximating the left hand side with the trapezoidal formula we obtain VL (t) + VL (t − ∆t) ∆t = L[iL (t) − iL (t − ∆t)] 2  (3.5)  Then the voltage across the inductor is  VL (t) =  2L 2L iL (t) + [−VL (t − ∆t) − iL (t − ∆t)] ∆t ∆t Requiv  (3.6)  history  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 oscillations 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  2L/ t iL  +eh(t) vL  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 algebraic sum of currents entering a node is zero, leading to the formulation of a system with the form [Y ][V ] = [J]  32  (3.7)  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  Y21      Y12  V1  J1    =   Y22 V2 J2  (3.9)  then Y11 V1 = J1 − Y12 V2  (3.10)  where, V1 is the vector of N 1 nodes with unknown voltages referred to ground V2 is the vector of N 2 nodes with voltages sources connected between the node and ground (known voltages) Y11 , Y12 , Y21 , Y22 are admittances submatrices of order N 1 x N 1, N 1 x N 2, N 2 x N 1, N 2 x N 2, respectively J1 is the vector of current sources entering the N 1 unknown voltages nodes 33  Chapter 3. Eigenvalue Analysis with EMTP Solution J2 is the vector of current sources entering the N 2 known voltages nodes Equation 3.10 constitutes a system of N 1 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 admittances 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 0 Hz to the Nyquist frequency fN y ,  fN y =  fsample 1 = 2 2∆t  (3.11)  The smaller the integration time step, the higher Nyquist frequency. There  34  Chapter 3. Eigenvalue Analysis with EMTP Solution 1.6  1.4  He/H, Magnitude distortion (p.u.)  Backward & Forward Euler 1.2  1  0.8  Trapezoidal 0.6  Calviño 3 0.4  0.2  0  0  0.05  0.1  0.15  0.2  0.25  f(pu) = f  0.3  0.35  0.4  0.45  0.5  (Hz)/f (Hz)  sim  Ny  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 network solution approach the Nyquist frequency, the accuracy of the solution decreases (“distortion error”). The distortion error introduced by the chosen integration rule increases as the frequency gets closer to the Nyquist frequency. 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) / fN y (Hz). For a given frequency in the simulation, the error will decrease if the time 100  Backward Euler  Phase distortion [degrees]  80  60  40  Calviño 3  20  Trapezoidal  exact  0  −40  Φ  He  −Φ  −20  −60  Forward Euler  −80  −100  0  0.05  0.1  0.15  0.2  0.25  f(pu) = f  0.3  0.35  0.4  0.45  0.5  (Hz)/f (Hz)  sim  Ny  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 representation 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 every ∆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, ∆t ∆t + v(t − ∆t) 2L 2L  (3.12)  ∆t ∆t −1 V (z) + z V (z) 2L 2L  (3.13)  i(t) − i(t − ∆t) =  Applying the Z-transform, I(z) − z −1 I(z) =  Ye (z)Trap =  I(z) ∆t (z + 1) = V (z) 2L (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 discontinuities 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 z L (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˜ no 3 [54], which has a stable integrator and differentiator characteristic, we get the transfer function,  Ye (z)C3 =  ∆t 1 −1 k0 L [1 − k1 z + k2 z −2 − k3 z −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 (ejω∆t + 1) 2L (ejω∆t − 1)  (3.17)  Similarly, for backward Euler in (3.15) we obtain  Ye (w)BE  ∆t ejω∆t = L (ejω∆t − 1)  38  (3.18)  Chapter 3. Eigenvalue Analysis with EMTP Solution and for Calvi˜ no 3 in (3.16) we obtain  Ye (w)C3 =  ∆t 1 k0 L [1 − k1 e−jω∆t + k2 e−jω2∆t − k3 e−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.  y1(t)  u1(t) u2(t) u2(t)  internal states  y2(t) y2(t)  x1(t),x2(t),x2(t),...xn(t) yp(t)  um(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)   a11  .   .  ..  =  ..       an1 xn (t)   b11  .  ..   an1  . . . a1n   x1 (t − ∆t)   .  .. ..  + . ..  .     xn (t − ∆t) . . . ann   . . . b1m   u1 (t − ∆t)    ..  .. ...   .  .     . . . anm 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) c11  .   .  ..  =  ..       yp (t) cp1  d11  .  ..   dp1      . . . c1n   x1 (t)   .  . . . ..   .  .   . +   . . . cpn xn (t)   . . . d1m   u1 (t)   .  ..  ..  .  . .   .    um (t) . . . dpm  41  (3.21)  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 adapting 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 current 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 Electromagnetic Transients Program (EMTP),  G  v(t) = h(t) + u(t)  42  (3.24)  Chapter 3. Eigenvalue Analysis with EMTP Solution  D  x(t+ t)  u(t)  x(t)  y(t)  Z-1  B  C  A  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 −1  v(t) = G  −1  h(t) + G  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 π 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) 44  (3.27)  Chapter 3. Eigenvalue Analysis with EMTP Solution g = t / 2L  L m  k  k  i  m h(t) vkm  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)  45  (3.30)  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) = hkm (t − ∆t) + (−2g)vkm (t − ∆t) new state  old state  (3.32)  f orcing f unction  Compare this formula with (3.22) in the state space description. As mentioned earlier, even though in state space (3.22) corresponds to forward Euler discretization, in the EMTP (3.32) was obtained from trapezoidal discretization. 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 information 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 π-circuit representations of 46  Chapter 3. Eigenvalue Analysis with EMTP Solution g = 2C / t  C k  m  k  i  m h(t) vkm  Figure 3.8: Discrete equivalent of a Capacitor  transmission lines, equipment in HVDC converter stations such as surge capacitors 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)  47  (3.35)  Chapter 3. Eigenvalue Analysis with EMTP Solution By replacing (3.35) into (3.33), the state space description of a single capacitance is obtained,  hkm (t) = gvkm (t − ∆t) + gvkm (t − ∆t) − hkm (t − ∆t)  (3.36)  Rearranging,  hkm (t) = −hkm (t − ∆t) + 2g vkm (t − ∆t) new state  old state  (3.37)  f orcing f unction  Which has the same form as (3.22) for the state space description. Analogously to the inductance case, the one-element state equation for the capacitance 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  R k  m  m  k  h(t) = 0  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 element. 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, consequently, 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) − 2gRL vkm (t − ∆t)  (3.38)  where, gRL =  1 R+  49  2L ∆t  (3.39)  Chapter 3. Eigenvalue Analysis with EMTP Solution L  R  R  2L/ t iL  k  m  +m  k  eh(t) vL  gRL RRL = R + 2L/ t  iL  +-  k  m  eh(t) h(t) vRL vRL RL discrete equivalent  Figure 3.10: Discrete equivalent of a Series RL  Analogously for the series RC,  hkm (t) = −hkm (t − ∆t) + 2gRC vkm (t − ∆t)  gRC =  3.4.5  1 ∆t R + 2C  (3.40)  (3.41)  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 hL3 1  Branch 3 1  2  2  gL3  L1  C2  Branch 2  Branch 1  L3  gL1  hL1  gC2  hC2  Figure 3.11: Generic circuit  The nodal equations for the generic circuit of Fig. 3.11 are  [G][v] = [j]   G11  G21  (3.42)      G12  v1  j1    =   G22 v2 j2  In detail:           −gL3  v1 (t)  hL1 (t) + hL3 (t)  gL1 + gL3  =   hC2 (t) − hL3 (t) v2 (t) −gL3 gC2 + gL3     j (t) v1 (t) −1  1   = [G ]   j2 (t) v2 (t)  51  Chapter 3. Eigenvalue Analysis with EMTP Solution       v1 (t) R11 R12  j1 (t)  =   v2 (t) R21 R22 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 h i . We know from equations (3.32) and (3.37) that for each branch element L or C, the branch histories are:        0   h 1 (t − ∆t)  h 1 (t) k∧1 0       h (t) =  0 k   0  ∧2  2     h 2 (t − ∆t) +      h 3 (t) 0 0 k∧3 h 3 (t − ∆t)    0   v 1 (t − ∆t) g∧1 0     0 g   0  ∧2    v 2 (t − ∆t)    0 0 g∧3 v 3 (t − ∆t)  (3.44)  In particular for our circuit, g∧1 =  −∆t L1  k∧2 = −1  g∧2 =  4C2 ∆t  k∧3 =  g∧3 =  −∆t L3  k∧1 =  1  1  (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 branchnode incidence matrix:       v 1 (t − ∆t) 1 0   v1 (t − ∆t)     v (t − ∆t) = 0 1      2    v2 (t − ∆t)    1 −1 v 3 (t − ∆t)   (3.46)  Substituting (3.43) in (3.46),       v 1 (t − ∆t) 1 0      R11 R12  j1 (t − ∆t)  v (t − ∆t) = 0 1      2      R21 R22   j2 (t − ∆t) 1 −1 v 3 (t − ∆t)     (3.47)  The nodal history sources in terms of the branch histories are given by     j1 (t − ∆t)  h 1 (t − ∆t) + h 3 (t − ∆t)   =   j2 (t − ∆t) h 2 (t − ∆t) − h 3 (t − ∆t)    h (t − ∆t)   1   1 0 1    =    h 2 (t − ∆t)   0 1 −1  h 3 (t − ∆t)  53  (3.48)  Chapter 3. Eigenvalue Analysis with EMTP Solution substituting (3.47), (3.48) in (3.44),       0   h 1 (t − ∆t)  h 1 (t) k∧1 0        h (t) =  0 k  0  ∧2  2     h 2 (t − ∆t) +      h 3 (t) 0 0 k∧3 h 3 (t − ∆t) [k∧ ]         0  1 0   g∧1 0    R11 R12   0 g  0 1    0 ∧2       R21 R22  0 0 g∧3 1 −1 [G−1 ] [g∧ ]  [L]    h (t − ∆t)  1  1 0 1       h 2 (t − ∆t)    0 1 −1  h 3 (t − ∆t)     [Lt ]  −1  h(t) = k∧  h(t − ∆t) + g∧  L  G  t  L  h(t − ∆t)  (3.49)  The discrete time transition matrix [A]d is then given by d  A  = 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 determined 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 component 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 necessary 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  1  L1  C2  1  L3  C4  L5  2  Leq  Ceq  2  Figure 3.12: Aggregate treatment of parallel RLC branches        h 1 (t) k∧1 0 0 0 0 h 1 (t − ∆t)            h 2 (t)  0 k∧2 0   0 0       h 2 (t − ∆t)       h (t) =  0   h (t − ∆t) + 0 k 0 0 ∧3  3    3             h 4 (t)  0 0 0 k∧4 0   h 4 (t − ∆t)      h 5 (t) 0 0 0 0 k∧5 h 5 (t − ∆t)    g∧1 0 0 0 0 v 1 (t − ∆t)        0 g∧2 0   0 0   v 2 (t − ∆t)       0   v (t − ∆t) 0 g 0 0 ∧3 3           0 0 0 g∧4 0   v 4 (t − ∆t)    0 0 0 0 g∧5 v 5 (t − ∆t)  56  (3.51)  Chapter 3. Eigenvalue Analysis with EMTP Solution   1 v (t − ∆t)    1     v 2 (t − ∆t) 1        v (t − ∆t) = 1    3        v 4 (t − ∆t) 1    1 v 5 (t − ∆t)    −1    −1   v1 (t − ∆t)  −1   v2 (t − ∆t)  −1  −1  (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, d  A  3.5  = k∧ + g∧  L  G−1  Lt  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  tan α2 = L2  λ Saturation tan α1 = L1  λ Saturation  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 corresponds 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 1 1 = + L2 L1 Lp  (3.53)  Lp λp L1 λ1 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 involved 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 eigenvalues from steady state power system admittance matrices whose elements are complex numbers. Notice also that in the important case of eigenvalue calculations with nonlinearities 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 eigenvalues 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 discretization 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 original 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,  x(t) ˙ = Ac x(t)  (3.54)  and its corresponding discrete time state space equivalent given by 3.22, x(t + ∆t) = Ad x(t)  (3.55)  Discretizing 3.54 with the trapezoidal integration rule,  x(t ˙ + ∆t) =  2 2 x(t + ∆t) − h(t + ∆t) ∆t ∆t  (3.56)  where h(t + ∆t) = x(t) +  ∆t x(t) ˙ 2  (3.57)  substituting 3.57 in 3.56,  x(t ˙ + ∆t) =  2 {x(t + ∆t) − x(t)} − x(t) ˙ ∆t  (3.58)  where x(t) ˙ = Ac x(t) and x(t ˙ + ∆t) = Ac x(t + ∆t), Ac x(t + ∆t) =  2 {x(t + ∆t) − x(t)} − Ac x(t) ∆t  61  (3.59)  Chapter 3. Eigenvalue Analysis with EMTP Solution solving for x(t + ∆t), 2  I − Ac I + Ac ∆t  x(t + ∆t) = x(t) ∆t 2  (3.60)  since, x(t + ∆t) = Ad x(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 Ac vi = λi vi  (3.63)  The continuous time eigenvalues λi of Ad can be computed immediately from those discrete time eigenvalues zi of Ad , and vice-versa, as follows, zi − 1 zi + 1  (3.64)  1 + λi 1 − λi  (3.65)  λi = Ω  zi = Ω  Equation (3.65) corresponds to the classical Bilinear transformation in signal processing theory. In the case of the trapezoidal rule (Bilinear transforma62  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 circle 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 discussed 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  Imaginary  Imaginary  Unit circle  λ  z Real  Real z  λ  Complex continuous-time plane  Complex discrete-time plane  Figure 3.15: Continuous time and discrete time stability mapping for trapezoidal  Imaginary  Imaginary Unit circle  λ  z Real  Real z  λ  Complex continuous-time plane  Complex discrete-time plane  Figure 3.16: Continuous time and discrete time stability mapping for backward Euler  64  Chapter 3. Eigenvalue Analysis with EMTP Solution However, in the case of systems with nonlinear elements the discrete to continuous 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 characteristics 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 instability. – 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 corresponds 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 realtime simulator. OVNI’s discrete time solution is based on the EMTP discretization 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´evenin Equivalent (MATE) concept, which is described in the following section, and was developed as an object oriented application written in C++. This characteristic facilitates the integration of new models and modules with OVNI. The OVNI solution achieved real-time performance in a single PC environment allowing interaction with analog devices for protective relay testing applications, as introduced by Calvi˜ no in [16]. The OVNI PC-cluster solution developed by Hollman [17] showed the OVNI algorithm’s inherent advantages for parallelization allowing for real-time simulation of large networks. The introduction of the Multilevel MATE concept by Armstrong [57] provides OVNI with the capability for efficient simultaneous solution of control 67  Chapter 4. Application to the OVNI simulator  iα +A  B z  Vα  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  Link α1  A3  Zα1  gA23  iα1  gA34 Link α2  A4  [A]  A2  Zα2 hA2  hA4 gA21  B1 gB13  B3  [B]  hB1 gB12  iα2 gB32  gA14  B2  A1 hB2  hA1  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   vA   h A  =    0 B vB hB  69  (4.1)  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  GA21 GA22 GA23  0  GA31 GA32 GA33 GA34  1  GA41  0  0  GA34 GA44 GB11 GB12 GB13  −1  GB21 GB22 GB23  0  GB31 GB32 GB33  0  0  0  1  0  −1  0  0  0  0  1  0  0  0                             70  −1  0  −z1   0    vA1       vA2        vA3       vA4        = vB1       vB2        vB3        iα1      iα2     0    0    1     0    0     −1    0    −z2  hA1    hA2    hA3    hA4     hB1    hB2     hB3    0    0  (4.2)  Chapter 4. Application to the OVNI simulator with its compact notation,            A 0 p   0 B q   pt q t −z    vA   v  B  iα    hA   = h   B   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).     0  0   0 0  p=   +1 0  0 +1            (4.4)     −1 0  q= 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,    a13  1   a23 1    1 a33     1 a43    1 −b11    1 −b21     1 −b31    Zα11   Zα21  a14     a24     a34      a44     −b13     −b23      −b33     Zα12    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−1 p  (4.7)  b = B −1 q  (4.8)  eA = A−1 hA  (4.9)  eB = B −1 hB  (4.10)  e = pt eA + q t eB  (4.11)  Z = pt a + q t b + 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 subnetwork 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 simulation 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 central 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 contents of the segmented blocks’ transition matrices are located in the correct 73  Chapter 4. Application to the OVNI simulator No update required  Update is required only when change of topology  Linear subnetwork A1  Eigenvalue solver  Linear subnetwork A2 with topology change  Update is required only when change of region  nonlinear subnetwork A3  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 during 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 network 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 dynamics of elements is updated in the stability analysis [59], [60]. As an alternative to the proposed centralized computation of eigenvalue Soft real-time environment Interaction required only when there is a change in [A] i  Real-time environment Subsystem I links links  Subsystem II  links  Dedicated Eigenvalue solver  OVNI central solver  Subsystem III  Interactions required at every time step  Figure 4.4: Hybrid real-time/soft real-time simulator layout  trajectories, a segmented distributed computation can be implemented ex75  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 concentrating the computation of eigenvalues within a much smaller structure. In an extension of the Latency concept discussed in the next section, the subsystem 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 Moreira 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...”. 3  In 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 dominant 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 = [p1 p2 . . . pk ]  with    (4.13)     φ1i ϕi1   φ2i ϕi2  pi =  ..  .   φ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 partic77  Chapter 4. Application to the OVNI simulator ipation of the k th state in the ith mode. The participation matrix combines right and left eigenvectors showing the association between state space variables and modes. The relationship between states and modes can not easily be extracted directly from the right and left eigenvectors due to their dependency 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 formulation 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 fmax system > fmax user Issue user notification  no  Is fmax system of interest? yes Select best appropriate time step  Run simulation with fixed user time step  Run simulation  no  Change of topology or nonlinear operation region?  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 technique, 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  vL(t)  R  L  i(t )  C  Vin(t)  vC(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  x(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 ;  82  eC =  1 C  i(t) dt  Chapter 5. Validation i(t) is the current in the loop. Substituting di(t) 1 + R i(t) + dt C  vin (t) = L  i(t) dt  (5.4)  Introducing the capacitor charge as a variable, i(t) = dq(t)/dt,  vin (t) = L  d2 q(t) dq(t) 1 +R + q(t) 2 dt dt C  (5.5)  If we select q(t) and i(t) as the state variables and define the input and output as  q(t) x(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          1  q(t) 0  +   vin (t)   1 −R i(t) input L L  0 =  −1 LC  rate of change  [A]  present states  (5.7)  [B]   q(t)  + 1 vin (t) −R  i(t) [D] input   vL (t) = output  −1 C  [C] states  83  (5.8)  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  c  A c  λ    1   0 =   −800 −20   −10.0000 + j26.4575 =   −10.0000 − j26.4575  (5.9)  (5.10)  Another possible representation of the system can be obtained by selecting vR (t) and vC (t) as the state variables. In this case the state space equations for the RLC series circuit are     vR (t) x(t) =   vC (t) u(t) = vin (t)  (5.11)  y(t) = vL (t)      dvR (t)  dt     dvC (t) dt     =  rate of change      −R  L  −R v (t) L  R    R  L  +   vin (t)  vC (t) 0 0 input  −1 RC  states  [A]  84  [B]  (5.12)  Chapter 5. Validation     vR (t) vL (t) = −1 −1   + 1 vin (t) v (t) [D] input C output  (5.13)  [C] states  The obtained continuous time transition matrix [A]c and the eigenvalues [λ]c are   c  A c  λ   −20 −20 =   40 0   −10.0000 + j26.4575 =   −10.0000 − j26.4575  (5.14)  (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  hL(t)  R  1  2  3 gL equiv  Vin(t)  gC equiv  hC(t)  Figure 5.2: Discretized RLC series circuit  1 R     −1 R  0      0 i1  v1             −∆t = ( R1 + −∆t ) v h 2 L      2L 2L     ∆t ∆t 2C ( 2L + ∆t ) v3 −hL + hC 2L −1 R  (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          ( R1  +  ∆t ) 2L  −∆t 2L −1 R    +  ∆t ) 2L  −∆t 2L        −1 R  hL    v2        ( ∆t + 2C ) 0    v3  =  −hL + hC 2L ∆t    1 0 v1 i1 R       hA   GAA GAB   vA    =    GBA GBB vB hB   1 ( R  −∆t 2L        −v1 hL  v2     R    =  −  + 2C ) v −h + h 0 3 L C ∆t  −∆t 2L  ( ∆t 2L [G]  [h]        (5.17)  (5.18)  (5.19)  input      −1  v1 1 ∆t −∆t v2  ( R + 2L )   hL + R  2L  =    −∆t ∆t 2C −h + h v3 ( + ) L C 2L 2L ∆t  (5.20)  The relation between branch and node voltages and history sources is given by the incidence matrix [L]:           v 1   1 0      v2  v  = −1 1   2        v3 v3 0 1  (5.21)     h     R   1 −1 0 h 2        = −hL    0 1 1 h3 hC  (5.22)  87  Chapter 5. Validation Substituting the nodal voltages (5.20) into (5.21),       −1  v 1   1 0 1 ∆t v1 −∆t     ( R + 2L )   hL + R  2L v  = −1 1      2    −∆t ∆t 2C     ( 2L + ∆t ) −hL + hC 2L v3 0 1  (5.23)  For the considered RLC series circuit, [k∧ ] and [g∧ ] (3.45) are  k∧  g∧    0 0 0     =  0 1 0    0 0 −1   0 0 0   −∆t =  0  0 L   4C 0 0 ∆t  (5.24)  (5.25)  while the branch histories, according to equation (3.49), can be obtained as     0  1 0   h 1  0 0 0  h 1  0 0          h  = 0 1 0  h  + 0 −∆t 0  0 1   2      2  L         0 0 4C h3 0 0 −1 1 −1 h3 ∆t   −1   h   1 ∆t −∆t 1  1 0 1 + ) (     R 2L 2L      h 2   ∆t 2C −∆t 0 1 −1   ( 2L + ∆t ) 2L h3           (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),  d  A          0  1 0  0 0 0  0 0        −∆t 0  0 1  = 0 1 0  + 0 L        0 0 4C 0 0 −1 1 −1 ∆t  −1   ∆t −∆t 1 ( R + 2L )  1 0 1  2L     −∆t ∆t 2C ( + ) 0 1 −1 2L 2L ∆t  (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  d  A   0  0  0      −3 9.97998−1 −3.99599−6  = 1.99799    1.99799−3 1.99799 9.99996−1  (5.28)  The discrete time eigenvalues of this matrix are  d  z    −1 + j2.64310−3 9.98997      −1 −3 = 9.98997 − j2.64310    0  (5.29)  The eigenvalues of the transition matrix [A] are the poles of the corresponding 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 computed 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 (zi − 1) ∆t (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 eigenvalues from the discrete time eigenvalues, obtained with equation (3.64), with a ∆t of 100µs are     c  =  λ reconstructed  2 −0.00050 + j0.00132   100µs −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 −0.0050 + j0.01322   1ms −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ω 91  (5.33)  Chapter 5. Validation then f=  ω 2π  (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 frequency and damping ratio are f c = 4.21Hz  (5.36)  ζ c = 0.35355  (5.37)  The LC resonant frequency of the circuit fnc (“natural frequency”) is fnc =  1 √ = 4.50Hz 2π LC  (5.38)  and the damping ratio ζ c is R  ζc = 2  L C  = 0.35355  (5.39)  Notice that in a stable system due to the losses the oscillation frequency 92  Chapter 5. Validation  State 2  State 4  L1  Vgen  R1  L2  C1  State 1  C2  State 3  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 implementation 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 continuous 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 dvdtC1 =⇒ x˙ 1 =  1 i C1 C1  vL1 = L1 didtL1 =⇒ x˙ 2 =  1 v L1 L1  iC2 = C2 dvdtC2 =⇒ x˙ 3 =  1 i C2 C2  vL2 = L2 didtL2 =⇒ x˙ 4 =  1 i L2 L2  (5.40)  Applying Kirchoff’s laws, x˙ 1 =  1 (v C1 L1  x˙ 2 =  1 (Vgen L1  x˙ 3 =  1 v C2 L2  x˙ 4 =  1 (v L2 C1  − vL2 ) − vC1 )  (5.41)  − vC2 + R1 vL2 )  The continuous time [A] matrix for the test circuit is given by      0   −1  L c A = 1  0    1 C1  0  0  0  0  0  1 L2  0  −1 L2  94  −1 C1    0     1 C2    −R1 L2  (5.42)  Chapter 5. Validation from which the following eigenvalues are obtained    −4.99950x104 + j1.00379x106 λ1      4 6  −4.99950x10 − j1.00379x10 λ2  c   λ =  4  −4.99801 + j9.94988x10 λ3    −4.99801 − j9.94988x104  (5.43)  λ4  We next compute the eigenvalues and eigenvectors using the discrete time methodology. The network’s [G] matrix from the discrete time description is    ∆t 2L1  −∆t 2L1     −∆t ( 1 + ∆t +  2L R1 2L1 G= 1  0 −1  R1  0 0  2C1 ) ∆t  0  0  −1 R1  0  ( R11 +  ∆t ) 2L2  −∆t 2L2  −∆t 2L2 ∆t ( 2L + 2           (5.44)  2C2 ) ∆t  the matrix [k∧ ] is given by      1 0 0 0   0 −1 0 0  k∧ =   0 0 1 0   0 0 0 −1  95           (5.45)  Chapter 5. Validation while the matrix [g∧ ] is given by    −∆t L1     0  g∧ =    0  0  0  0  2C1 ∆t  0  0  −∆t L2  0  0  0    0     0    (5.46)  4C2 ∆t  The incidence matrix [L]     1 −1 0 0      0 1 0 0    L=     0 0 1 −1    0 0 0 1  (5.47)  Applying (3.50), we obtain the discrete time transition matrix [A]d d  A  = k∧ + g∧  L  G−1  Lt     9.99999x10−1     −2.00000 d A =   1.24999x10−15  −1.24999x10−15  1.25000x10−15  1.24999x10−15  0  9.99999x10−1  1.99999  1.25034x10−13  −1.24999x10−15  9.99999x10−1  1.24999x10−13  1.24999x10−15  −1.99999  9.99999x10−1          (5.48)  96  Chapter 5. Validation from where the following discrete time eigenvalues are computed,    9.99999x10−1 + j4.97494x10−8     9.99999x10−1 + j4.97494x10−8 d λ =   9.99999x10−1 + j5.01896x10−7           (5.49)  9.99999x10−1 − j5.01896x10−7  Finally, the reconstructed eigenvalues from the discrete time domain given by (3.64) are,    −4.99752 +  λcreconstructed  j9.94988x104     −4.99752 − j9.94988x104 =   −4.99950x104 + j1.00379x106           (5.50)  −4.99950x104 − j1.00379x106  which matches up to the second decimal for the worse approximation the continuous time eigenvalues calculated directly from the continuous time solution (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  λ3 creconstructed  λ4 creconstructed  4.9995x10−1 + j2.6108x10−5 4.9995x10−1 − j2.6108x10−5 x˙ 2 c Prec. = 4.9495x10−1 + j7.5560x10−5 4.9495x10−1 − j7.5560x10−5 x˙ 1  4.9970x10−5 − j9.9168x10−7 4.9970x10−5 + j9.9168x10−7 x˙ 4 5.0474x10−3 − j1.0067x10−4 5.0474x10−3 + j1.0067x10−4 x˙ 3  λ2 creconstructed  λ1 creconstructed  4.9470x10−5 − j7.5406x10−6 4.9470x10−5 + j7.5406x10−6 x˙ 2 5.0479x10−3 − j2.6137x10−4 5.0479x10−3 + j2.6137x10−4 x˙ 1 4.9995x10−1 + j2.4910x10−2 4.9995x10−1 − j2.4910x10−2 x˙ 4 4.9495x10−1 − j2.4641x10−2 4.9495x10−1 + j2.4641x10−2 x˙ 3 (5.51)  As can be seen, the eigenvalues λ1 creconstructed and λ2 creconstructed are strongly associated to the states of x˙ 3 and x˙ 4 (corresponding to C2 and L2 ), while eigenvalues λ3 creconstructed and λ4 creconstructed are strongly associated to the states of x˙ 1 and x˙ 2 (corresponding to C1 and L1 ). The oscillation frequencies given by (5.34) are 15.8kHz for the slow subnetwork 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 State 2  State 4  L1  Vgen  L2  R1 C1  State 1  C2  Slow  State 3  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. State 1 1  SW1  2 R1  Vgen  State 4 3  L1  4 R2  L2  State 3  State 2 C1  Rload  C2  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  3Ω  5Ω 10nF  C2  L1  L2  Vgen  ∆t  50nF  350mH  10mH  230kV  50µs  and for the improved damped case II, R1  R2  3Ω 5Ω  C1  C2  L1  L2  Vgen  ∆t  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    1 L R1 +2 ∆t1      − 1 L1  R1 +2 ∆t [Gopen ] =   0    0  −1 L R1 +2 ∆t1 1 L R1 +2 ∆t1  +  2C1 ∆t  +  1 Ropen  −1 Ropen  1 Ropen  +  0  0  −1 Ropen  0  2C2 ∆t  +  −1 L R2 +2 ∆t2  0  100  1 L R2 +2 ∆t2  −1 L R2 +2 ∆t2 1 L R2 +2 ∆t2  +  (5.52)  1 Rload             Chapter 5. Validation  and with the switch SW1 closed,    1 L R1 +2 ∆t1      − 1 L1  R1 +2 ∆t [Gclose ] =   0    0  −1 L R1 +2 ∆t1 1 L R1 +2 ∆t1  +  2C1 ∆t  +  1 Rclose 1 Rclose  −1 Rclose  +  0  0  −1 Rclose  0  2C2 ∆t  +  1 L R2 +2 ∆t2  −1 L R2 +2 ∆t2  0  −1 L R2 +2 ∆t2 1 L R2 +2 ∆t2  +  (5.53)  the [k∧ ] matrix of branch nature coefficients is given by     1 0 0 0 0    0 −1 0 0 0   k∧ =   0 0 −1 0 0   0 1 0  0 0  0 0 0 0 0              (5.54)  while [g∧ ] matrix of branch discretization rule coefficients is        g∧ =         −2 2L R1 + ∆t1  0  0  0  0  4C1 ∆t  0  0  0  0  4C2 ∆t  0  0  0  0  −2 2L R2 + ∆t2  0  0  0  0  101  0    0    0    0    0  (5.55)  1 Rload             Chapter 5. Validation and the incidence matrix [L] is     1 −1 0 0      0 1 0 0       L= 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    −2  λdclosed  −1   −1.419878x10 + j9.998990x10   −1.419878x10−2 − j9.998990x10−1  =  −9.999999x10−1   −1.000000   (5.57)   −1  λdopen           −1   −1.049723x10 + j9.944750x10   −1.049723x10−1 − j9.944750x10−1  =  9.999999x10−1   −1.000000           (5.58)  while the reconstructed continuous time eigenvalues computed with (3.64)  102  103  [V] [V] 0  0.5  −2  0  2  −1  0  1  −5  0  −0.5  [V] [V]  5  0  0.01  0.015  0.01  0.015  0.005  0.01  0.015  Voltage Node 4  0.005  Voltage Node 3  4 0.005 0 x 10  0  0.015  Voltage Node 2  0.01  Voltage Node 1  6 0.005 0 x 10  5  x 10  0.02 0.025 time [s]  0.02 0.025 time [s]  0.02 0.025 time [s]  0.02 0.025 time [s]  0.03  0.03  0.03  0.03  0.035  0.035  0.035  0.035  0.04  0.04  0.04  0.04  0.045  0.045  0.045  0.045  0.05  0.05  0.05  0.05  Chapter 5. Validation  Figure 5.6: Node voltages for the two-area example of Fig. 5.5. Case I, poorly damped  Chapter 5. Validation  0.05  time [s]  0.04  0.03 eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  0.02  0.01  0 1 1  0.5 0.5  0 0  Imag  −0.5  −0.5 −1  Real  −1  Figure 5.7: Discrete time eigenvalues for the two-area example of Fig. 5.5. Case I, poorly damped are     −3  λcclosed  4  −5.072087x10 + j4.057204x10    −5.072087x10−3 − j4.057204x104  =  −1.200000x1012   −4.000000x104   (5.59)   −3  λcopen           4   −4.938281x10 + j4.444444x10   −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  0.05  0.04  eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  time [s]  0.03  0.02  0.01  0 5 0 −2 −4  0  4  x 10  −6 −8  Imag  −10 −5  −12 −14  11  x 10  Real  Figure 5.8: Reconstructed continuous time eigenvalues for the two-area example 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,    −1  λdclosed  −1   2.902763x10 + j9.569413x10   2.902763x10−1 − j9.569413x10−1  =  −9.99999x10−1   −1.00000           (5.63)   −1 −1 2.857138x10 + j9.583134x10      2.857138x10−1 − j9.583134x10−1    =    −1 −9.99999x10     −1.00000  (5.64)    λdopen  while the reconstructed continuous eigenvalues are    −2  λcclosed  4  −4.444838x10 + j2.966627x10    −4.444838x10−2 − j2.966627x104  =  −1.009999x1012   −4.000000x104 106           (5.65)  107  [V] [V] 0  0.5  −2  0  2  −1  0  1  −5  0  −0.5  [V] [V]  5  0  0  0.01  0.01  0.005  0.01  Voltage Node 4  0.005  Voltage Node 3  4 0.005 0 x 10  Voltage Node 2  0.01  Voltage Node 1  6 0.005 0 x 10  5  x 10  0.015  0.015  0.015  0.015  0.02 0.025 time [s]  0.02 0.025 time [s]  0.02 0.025 time [s]  0.02 0.025 time [s]  0.03  0.03  0.03  0.03  0.035  0.035  0.035  0.035  0.04  0.04  0.04  0.04  0.045  0.045  0.045  0.045  0.05  0.05  0.05  0.05  Chapter 5. Validation  Figure 5.9: Node voltages for the two-area example of Fig. 5.5. Case II, improved damping  Chapter 5. Validation  0.05  0.04  time [s]  0.03 eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  0.02  0.01  0 1 0.5 0  1  Imag  0.5  −0.5  0 −1  −0.5 −1  Real  Figure 5.10: Discrete time eigenvalues for the two-area example of Fig. 5.5. Case II, improved damping    −2  λcopen  4  −4.444444x10 + j2.981423x10    −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  108  (5.67)  Chapter 5. Validation  0.05 eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  0.04  time [s]  0.03  0.02  0.01  0 3 2  0  1  −2 0  4  x 10  −4 −6  −1 Imag  11  −8  −2  −10 −3  x 10 Real  −12  Figure 5.11: Reconstructed continuous time eigenvalues for the two-area example 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  State 1 SW1 R1  Vgen  R2  L1 L2  State 3  State 2 C1  C2  State 4 Non lineal Load  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 eigenvalue trajectories in systems with poorly damped eigenvalues that are being heavily stressed. These systems can present very critical operational trajectories. 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  time [s]  1  0.5  0 2  1  0  −1  −2  −1.5  −1  Imag  −0.5  0  0.5  Real  1.5  1  eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  1.5 1  1.5  0.5  1 0  0.5 0  −0.5  −0.5  −1  Imag  −1 −1.5  −1.5  Real  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  time [s]  1 eig1 − L2 eig2 − C2 eig3 − C1 eig4 − L1  0.5  0 5 0 −5  4  x 10  −15  5  0  −5  −10  11  x 10  Imag  Real  5  4  x 10  0  Imag  −5 −14  −12  −10  −8  −6  −4  −2  Real  0  2  11  x 10  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  0  0  1  0  1  Figure 5.15: Two-area system with nonlinear load −1  0  1  −1  voltage [V]  113 voltage [V]  −1  voltage [V]  −5  voltage [V]  5  5  0  6 0 x 10  6 0 x 10  6 0 x 10  x 10  0.1  0.1  0.1  0.1  0.2  0.2  0.2  0.2  0.3  0.3  0.3  0.3  0.4  0.4  0.4  0.4  0.5 time [s]  0.5 time [s]  0.5 time [s]  0.5 time [s]  0.6  0.6  0.6  0.6  0.7  0.7  0.7  0.7  0.8  0.8  0.8  0.8  0.9  Voltage Node 4  0.9  Voltage Node 3  0.9  Voltage Node 2  0.9  Voltage Node 1  1  1  1  1  Chapter 5. Validation  Chapter 5. Validation V1  Vgen  V2  P + jQ  Zline  Exp. load  Pcte  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 φgrid −φload  2 V ej = φ −φ V0 2 cos( grid 2 load )  (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 | VV0 | 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, V V0  α  P = P0  V V0  β  Q = Q0  (5.71) (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 RM S 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 methodology 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.1 1  A1 A2 Voltage Bus 2 (pu)  0.8  0.6  B C 0.4  D  0.2  0  1  2  3  4  5  6  7  8  9  10  time [s]  Figure 5.17: Voltage collapse of a simple radial system  116  11  12  Chapter 5. Validation  D  12  10  C  time [s]  8  6  B  A2  A1 4  2  0 2 1 0  Imag  −1  −1  −2  −0.8  −0.6  −0.4  −0.2  Real  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  0  Chapter 5. Validation 2  D B A1 C  1.5  1  0.5  A2  Imag  A2 0  −0.5  C −1  A1 −1.5  B D −2 −45  −40  −35  −30  −25  −20  −15  −10  −5  0  5  Real  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  D C 12  B  time [s]  10  A2  8  A1 6 4 2 1 0.99  0 2  1.5 −3  x 10  0.98 1  0.5  0  0.97 −0.5  −1  0.96 −1.5  Imag  −2  0.95  Real  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  D  12 11 10 9  [s]  8  C  B  7 6 5  A2  A1  4 3 1.5 −3  1 0.9998  1 0.9996  x 10  0.5  Imag  0.9994 0.9992  0  Real  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 initiate 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 1  0.9998  A1  0.9996  B  |eig 1|  D  A2  0.9994  0.9992  0.999 C  0.9988  0  2  4  6  8  10  12  time [s]  Figure 5.22: Discrete time eigenvalue 1 magnitude, |eig1 |, over time  that eigenvalue 1 starts to show evident signs of voltage collapse in the proximity 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 anticipation 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 connection 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 simulator [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 discrete 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 characteristic 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  Sw  123  V: 1.010 Ph: 6.92  V=1.01 P=10.5  Bus_2  X=0.0167 Tap=1.1291  T1  V: 1.075 Ph: 9.64  Bus_1  G2  Phase=20.2 Voltage=1.03  G1  V: 1.030 Ph: 20.20  PQ[7] P=9.67 Q=1  X=0.0167 Tap=1.1291  T2  25 km R=0.0025 X=0.025 B=0.04375  Bus_5  V: 1.000 Ph: -4.39  Q=-2.0  R=0.0010 X=0.010 B=0.0175  R=0.0020 X=0.020 B=0.0170  10 km  Bus_6  V: 0.954 Ph:-12.57  Bus_7  R=0.0075 X=0.075 B=0.1283  R=0.0075 X=0.075 B=0.1283  R=0.0075 X=0.075 B=0.1283  R=0.0075 X=0.075 B=0.1283  110 km  V: 0.872 Ph:-26.36  Bus_8  R=0.018 X=0.18 B=0.1925  R=0.018 X=0.18 B=0.1925  R=0.018 X=0.18 B=0.19  110 km  V: 0.986 Ph:-73.61  Q=-3.5  10 km  V: 1.029 Ph:-78.40  R=0.001 X=0.01 B=0.0175  PQ[9] P=1.767 Q=1  Bus_9  R=0.0025 X=0.025 B=0.0437  25 km  V: 1.010 Ph:-83.94  X=0.0167 Tap=1.08  T5  Bus_10  V: 1.030 Ph:-83.94  V=1.03 P=-6.8  Bus_3  Tolerance=1.0e-5 Ignore Tolerance=no Steady State Order=7 Steady State Frequency=60.0 Max Steps=10000 Stop in Collapse=y Compute RMS=1 MVAR Base=100 Collapse Voltage=0.3  Settings  G4  X=0.0167 Tap=1.1082  Bus_11 T3  V=1.01 P=-17.0  Bus_4  V: 1.088 Ph:-81.74  G3  Chapter 5. Validation  Figure 5.23: Two-area 11 bus test system  Bus Voltage Magnitude (pu)  124  0.85 0  0.9  0.95  1  1.05  1.1  1.15  1.2  1.25  5  10  15 Time (s)  20  25  30  Bus1 Bus2 Bus3 Bus4 Bus5 Bus6 Bus7 Bus8 Bus9 Bus10 Bus11  Chapter 5. Validation  Figure 5.24: Two-area 11 bus test system - Bus Voltages  125 0 0  0.2  0.4  0.6  0.8  1  200  pf=0.870  A  400  600  P=176.700, V=0.986  800  1000 1200 PL(MW)  1400  P=1583.783, V=0.566  1600  1800  B  (BUS−9)  2000  Chapter 5. Validation  Figure 5.25: Two-area 11 bus test system - PV curve at Bus 9 V(PU)  Chapter 5. Validation  0.25  0.2  0.15  0.1  Imag  0.05  0  −0.05  −0.1  −0.15  −0.2  −0.25 0.96  0.965  0.97  0.975  0.98  0.985  0.99  0.995  1  1.005  Real  Figure 5.26: Two-area 11 bus test system - Discrete time Eigenvalues normal operation point A  126  Chapter 5. Validation  0.25  0.2  0.15  0.1  Imag  0.05  0  −0.05  −0.1  −0.15  −0.2  −0.25 0.96  0.965  0.97  0.975  0.98  0.985  0.99  0.995  1  1.005  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 programs for power systems analysis based on the EMTP solution to produce not only the time domain waveforms of voltages and currents in the system, 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 information 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 system. 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 example 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 eigenvalue/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 capability 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 technique 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 ma129  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 analysis 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 analysis 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 contributing to instability.  6.2  Recommendations for Future Work  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 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. 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 renewable energy sources such as wind farms. These generation clusters are highly nonlinear and require advance on-site control devices to allow a healthy interconnection 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 Interdependencies Research Program (JIIRP) (sponsored by Public Safety and Emergency Preparedness Canada (PSEPC) and the Natural Sciences and Engineering Research Council of Canada (NSERC)) is concerned with the identification of critical interdependencies among multiple infrastructures (e.g., power grid, water grid, etc.) and the dynamic evolution of these interdependencies 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 embedded 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 methodology within UBC’s OVNI-NET simulator for stability analysis and determination of segmentation schemes. v. Implementation of discrete state space eigenvalue analysis methodology to UBC-JIIRP’s I2Sim simulator for identification of trajectories of critical 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 Apparatus and Systems, vol. PAS-88, no. 4, pp. 388–399, 1969.  134  Bibliography [7] G.D. Irwin, R. Kuffel, and D.A Woodford, “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 simulation,” 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 International Conference on, 1998, vol. 2, pp. 977–981 vol.2. [14] Fromonto, Levacher, and Strunz, “An efficient parallel computer architecture 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´ us Calvi˜ no Fraga, Implementation of a real time simulator for protective 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 Conference 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 architecture,” 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 software/hardware solution for real-time simulation of large power systems,” 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, “Application 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 Engineering 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 research,” 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 exponentially 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 transformer 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 residentialcommercial 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,” Computer 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, “Educational use of EMTP models for the study of rotating machine transients,” 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 pccluster,” 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´ us Calvi˜ no 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 network simulator,” in International Conference on Power System Transients, 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 nonlinearities in the OVNI simulator,” Power Systems, IEEE Transactions on, vol. 21, no. 3, 2006. [58] J.R. Mart´ı, “MATE, a multi-area the´ venin equivalent concept,” Internal 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 Workshop 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 mitigate critical events,” in International Workshop on Complex Network and Infrastructure Protection, Rome, 2006, CNIP. [64] Gene Golub and Charles Van Loan, Matrix Computation, Johns Hopkins 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−1 Pn−2 ....P2 P1 A = 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−1 Pn−2 ...P2 P1  (A.2)  then we have QT A = R and QQT A = 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 = R1 Q1 ; then factor A2 as A2 = Q2 R2 . 3. First set A3 = R2 Q2 ; then factor A3 as A3 = Q3 R3 . 4. Set Am = Rm−1 Qm−1 ; then factor Am−1 as Am = Qm Rm . 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 magnitudes, 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  rlc series, 85  Arnoldi  series branches, 50  generalized method, 17  trapezoidal, 44, 46  modified method, 26  Dommel, Hermann, xiv, 4, 5, 42  Bilinear Transformation, 63  Edison, Thomas Alva, 1  Brazil, 1  Eigenvalue, 16, 55 Electricite de France, 6  Calvi˜ no, Jes´ us, xiv, 7  EMTP, ii, 5  CDA, 31, 39 Chile, 1  fast time-domain, 17  continues time state space, 82 Gaspard Richie, Baron de Prony,  continuous time eigenvalues, 90  17 discrete eigenvalues, 89, 90  Gissinger, S., 3  discrete time state space, 4, 40 Helmholtz, 44  capacitance, 46  Hollman, Jorge, 7  inductance, 44  HVDC, 10  MATE, 73  Hydro Quebec Research Institute,  parallel branches, 55  6  resistance, 48 147  Index I2Sim, 133  RTDNS, 7  interarea oscillation, 15, 24  RTNS, 7  Jardim, J.L., 4  small-signal stability, 15  JIIRP, 132, 133  TEQSIM, 6  Kundur, Prabha, xiv, 3  Thevenin, 44 Transient Network Analyzer, 5  Linares, Luis, xiv, 7  transition matrix, 89  LTI, 40  trapezoidal, 90 Mart´ı, Jos´e, 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 148  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0065505/manifest

Comment

Related Items