UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

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

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

Item Metadata


831-ubc_2006_fall_hollman_jorge.pdf [ 2.47MB ]
JSON: 831-1.0065611.json
JSON-LD: 831-1.0065611-ld.json
RDF/XML (Pretty): 831-1.0065611-rdf.xml
RDF/JSON: 831-1.0065611-rdf.json
Turtle: 831-1.0065611-turtle.txt
N-Triples: 831-1.0065611-rdf-ntriples.txt
Original Record: 831-1.0065611-source.json
Full Text

Full Text

Step by Step EigenvalueAnalysis with EMTP DiscreteTime SolutionsbyJORGE ARIEL HOLLMANIndustrial Engineer, Universidad Nacional del Comahue, 1996M.A.Sc., The University of British Columbia, 2000A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate Studies(Electrical and Computer Engineering)The University Of British ColumbiaOctober, 2006c© Jorge Ariel Hollman 2006AbstractThe present work introduces a methodology to obtain a discrete time statespace representation of an electrical network using the nodal [G] matrix ofthe Electromagnetic Transients Program (EMTP) solution. This is the firsttime the connection between the EMTP nodal analysis solution and a cor-responding state-space formulation is presented. Compared to conventionalstate space solutions, the nodal EMTP solution is computationally muchmore efficient. Compared to the phasor solutions used in transient stabilityanalysis, the proposed approach captures a much wider range of eigenvaluesand system operating states. A fundamental advantage of extracting the sys-tem eigenvalues directly from the EMTP solution is the ability of the EMTPto follow the characteristics of nonlinearities. The system’s trajectory can beaccurately traced and the calculated eigenvalues and eigenvectors correctlyrepresent the system’s instantaneous dynamics. In addition, the algorithmcan be used as a tool to identify network partitioning subsystems suitablefor real-time hybrid power system simulator environments, including the im-plementation of multi-time scale solutions.The proposed technique can be implemented as an extension to any EMTP-based simulator. Within our UBC research group, it is aimed at extendingthe capabilities of our real-time PC-cluster Object Virtual Network Integra-tor (OVNI) simulator.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xNotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Research Background . . . . . . . . . . . . . . . . . . . . . . . 51.2.1 Off-Line Power System Simulation . . . . . . . . . . . 51.2.2 Real-Time Power System Simulation . . . . . . . . . . 61.2.3 Large Scale Power System Stability Analysis . . . . . . 91.3 Research Claim and Contributions . . . . . . . . . . . . . . . 111.4 List of Publications . . . . . . . . . . . . . . . . . . . . . . . . 122 Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 The Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.1 Traditional Analysis Tools . . . . . . . . . . . . . . . . 172.2.2 Advantages of Step by Step Trajectory Analysis . . . . 27iiiTable of Contents3 Eigenvalue Analysis with EMTP Solution . . . . . . . . . . . 293.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2 EMTP Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2.1 Nodal Analysis . . . . . . . . . . . . . . . . . . . . . . 323.2.2 Numerical Integration Accuracy . . . . . . . . . . . . . 343.2.3 Discretization Rule Stability . . . . . . . . . . . . . . . 373.3 Discrete Time State Space Formulation of an Electrical Net-work from the EMTP Nodal Matrix Formulation . . . . . . . 403.4 EMTP Solution in State Space form . . . . . . . . . . . . . . 423.4.1 Discrete Time State Space Equation for an Inductance 443.4.2 Discrete Time State Space Equation of a Capacitance . 463.4.3 Discrete Time State Space Equation of a Resistance . . 483.4.4 Treatment of Series Connection RL and RC . . . . . . 493.4.5 Treatment of Series Branches . . . . . . . . . . . . . . 503.4.6 Treatment of Parallel Branches . . . . . . . . . . . . . 553.5 Nonlinear Networks . . . . . . . . . . . . . . . . . . . . . . . . 573.6 Continuous Time Eigenvalues from Discrete Time Domain So-lution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.6.1 Discrete to Continuous Time Mapping . . . . . . . . . 603.6.2 Time Discretization Considerations . . . . . . . . . . . 634 Application to the OVNI simulator . . . . . . . . . . . . . . . 664.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.2 OVNI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2.1 The MATE Concept . . . . . . . . . . . . . . . . . . . 684.2.2 Discrete State Space Formulation with MATE . . . . . 734.3 The Latency Concept . . . . . . . . . . . . . . . . . . . . . . . 764.3.1 Application of the New Discrete Time State Space For-mulation to Latency . . . . . . . . . . . . . . . . . . . 774.4 Automatic Time Step Selection . . . . . . . . . . . . . . . . . 795 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.1 Comparison of State Space Solutions between Continuous Timeand Discrete Time Domains . . . . . . . . . . . . . . . . . . . 815.1.1 RLC Series Circuit State Space Equation from Contin-uous Time Domain . . . . . . . . . . . . . . . . . . . . 825.1.2 RLC Series Circuit State Space Equation from the [G]EMTP Matrix . . . . . . . . . . . . . . . . . . . . . . 85ivTable of Contents5.2 Latency Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 935.3 Two-Area Resonant Circuit . . . . . . . . . . . . . . . . . . . 995.3.1 Case I - Poorly Damped . . . . . . . . . . . . . . . . . 1005.3.2 Case II - Improved Damping . . . . . . . . . . . . . . . 1065.4 Two-Area Case with Nonlinear Load . . . . . . . . . . . . . . 1105.5 Voltage Collapse in a Radial System . . . . . . . . . . . . . . 1125.6 11 Bus Two Area System . . . . . . . . . . . . . . . . . . . . . 1226 Conclusions and Recommendations . . . . . . . . . . . . . . . 1286.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1286.2 Recommendations for Future Work . . . . . . . . . . . . . . . 131Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134A The QR Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 144Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147vList of Tables2.1 Relevant modes for the two-area system case . . . . . . . . . . 25viList of Figures1.1 Relation among research projects at UBC power system group 82.1 Linear time invariant (LTI) dynamic system . . . . . . . . . . 182.2 Two-area 11 bus test system . . . . . . . . . . . . . . . . . . . 222.3 Phasor value of bus 8 voltage for the two-area test system,small-signal instability observed with detailed time domainsimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.4 Phasor value of bus 8 voltage for the two-area test system,small-signal simulated with FTD (no oscillation are capturedby this solution) . . . . . . . . . . . . . . . . . . . . . . . . . . 243.1 Inductor element . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 Inductor discrete time equivalent with trapezoidal rule . . . . 323.3 Magnitude frequency response of integration rules . . . . . . . 353.4 Phase frequency response of integration rules . . . . . . . . . . 363.5 State space representation of a multiple input/output system . 393.6 State space system diagram block . . . . . . . . . . . . . . . 433.7 Discrete time equivalent of an Inductor . . . . . . . . . . . . . 453.8 Discrete equivalent of a Capacitor . . . . . . . . . . . . . . . . 473.9 Discrete equivalent of a Resistor . . . . . . . . . . . . . . . . . 493.10 Discrete equivalent of a Series RL . . . . . . . . . . . . . . . . 503.11 Generic circuit . . . . . . . . . . . . . . . . . . . . . . . . . . 513.12 Aggregate treatment of parallel RLC branches . . . . . . . . 563.13 Piecewise linear inductance with two slopes . . . . . . . . . . . 583.14 Implementation of two slope piecewise linear inductance . . . 593.15 Continuous time and discrete time stability mapping for trape-zoidal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.16 Continuous time and discrete time stability mapping for back-ward Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.1 Generic description of a MATE link . . . . . . . . . . . . . . . 68viiList of Figures4.2 Generic circuit for discrete state space formulation with MATE 694.3 General update scheme for discrete state space formulationwith MATE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.4 Hybrid real-time/soft real-time simulator layout . . . . . . . . 754.5 Automatic EMTP time step selection scheme . . . . . . . . . . 805.1 RLC series circuit . . . . . . . . . . . . . . . . . . . . . . . . . 825.2 Discretized RLC series circuit . . . . . . . . . . . . . . . . . . 865.3 Latency test circuit . . . . . . . . . . . . . . . . . . . . . . . . 935.4 Circuit with slow and fast subareas . . . . . . . . . . . . . . . 995.5 Two-area resonant test case . . . . . . . . . . . . . . . . . . . 995.6 Node voltages for the two-area example of Fig. 5.5. Case I,poorly damped . . . . . . . . . . . . . . . . . . . . . . . . . . 1035.7 Discrete time eigenvalues for the two-area example of Fig. 5.5.Case I, poorly damped . . . . . . . . . . . . . . . . . . . . . . 1045.8 Reconstructed continuous time eigenvalues for the two-areaexample of Fig. 5.5. Case I, poorly damped . . . . . . . . . . 1055.9 Node voltages for the two-area example of Fig. 5.5. Case II,improved damping . . . . . . . . . . . . . . . . . . . . . . . . 1075.10 Discrete time eigenvalues for the two-area example of Fig. 5.5.Case II, improved damping . . . . . . . . . . . . . . . . . . . . 1085.11 Reconstructed continuous time eigenvalues for the two-areaexample of Fig. 5.5. Case II, improved damping . . . . . . . . 1095.12 Two-area circuit test case with nonlinear load . . . . . . . . . 1105.13 Nonlinear load - discrete time eigenvalues . . . . . . . . . . . . 1115.14 Nonlinear load - continuous time reconstructed eigenvaluesfrom discrete time solution . . . . . . . . . . . . . . . . . . . . 1125.15 Two-area system with nonlinear load . . . . . . . . . . . . . . 1135.16 A simple radial system . . . . . . . . . . . . . . . . . . . . . . 1145.17 Voltage collapse of a simple radial system . . . . . . . . . . . . 1165.18 Reconstructed continuous time eigenvalue trajectory over timefor a simple radial system with voltage collapse . . . . . . . . 1175.19 Complex plane projection of continuous time eigenvalues overtime line of simple radial system of Fig. 5.16 with voltagecollapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.20 Discrete time eigenvalue trajectory from EMTP solution forthe simple radial system of Fig. 5.16 with voltage collapse . . 1195.21 Zoom in of discrete time trajectory of Fig. 5.20 . . . . . . . . 120viiiList of Figures5.22 Discrete time eigenvalue 1 magnitude, |eig1|, over time . . . . 1215.23 Two-area 11 bus test system . . . . . . . . . . . . . . . . . . . 1235.24 Two-area 11 bus test system - Bus Voltages . . . . . . . . . . 1245.25 Two-area 11 bus test system - PV curve at Bus 9 . . . . . . . 1255.26 Two-area 11 bus test system - Discrete time Eigenvalues nor-mal operation point A . . . . . . . . . . . . . . . . . . . . . . 1265.27 Two-area 11 bus test system - Discrete time Eigenvalues atcollapse point B . . . . . . . . . . . . . . . . . . . . . . . . . . 127ixAcronymsATA Advanced Technology Attachment, a disk drive implementation that inte-grates the controller on the disk drive itselfDSP Digital Signal ProcessingEDF Electricite´ de FranceEMTDC Transients Simulation Software by Manitoba HVDC Research CentreEMTP Electromagnetic Transients ProgramEPRI Electric Power Research Institute FACTS Flexible AC Transmission Sys-temsFTD Fast Time DomainGW Gigawatt, one billion times the power unit of wattGbps Gigabits per second, a data transfer speed measurement for high-speednetworksHVDC High Voltage Direct CurrentIDE Intelligent Drive Electronics, an interface for mass storage deviceIREQ Hydro Quebec Research InstituteLTI Linear Time InvariantMAM Modified Arnoldi MethodMATE Multi Area Thevenin EquivalentOVNI Object Virtual Network IntegratorRTDNS Real-Time Distributed Network SimulatorPCI Peripheral Component Interconnect, a local bus standard developed by IntelCorporationxAcronymsRTNS Real-Time Network SolverSATA Often abbreviated Serial ATA or S-ATA, an evolution of the Parallel ATAphysical storage interfaceTNA Transient Network AnalyzerTSAT Transient Stability Assessment ToolUBC The University of British ColumbiaUILO UBC’s University - Industry Liaison OfficeIPO Initial Public OfferingATP Alternative Transient ProgramCDA Critical Damping AdjustmentJIIRP Joint Infrastructure Interdependencies Research ProjectI2Sim Infrastructures Interdependencies SimulatorxiNotation[G] admittance matrix[GB] admittance matrix of subnetwork [B]GB43 element (4, 3) of admittance matrix [B]gB43 admittance of branch nodes 4 and 3 in subnetwork [B][L] incidence matrix[A]d discrete time transition matrix[A]c continuous time transition matrixλ continuous time eigenvaluez discrete time eigenvaluex˙(t) future state in continuous time domainx(t+∆t) future state in discrete time domaina′ compact notation for a(t−∆t)a compact notation for a(t)vkm branch voltage of node k to node mh(t)km history of branch km at (t)h(t−∆t)km history of branch km at (t−∆t)_h(t)i compact notation for history of branch i at (t)_h(t−∆t)i history of branch i at (t−∆t)_h i history of branch i at (t)_h′i compact notation for history of branch i at (t−∆t)[k∧] diagonal matrix of branch nature coefficients[g∧] diagonal matrix of branch discretization rule coefficientsζ damping ratioζc damping ratio from continuous time domainf oscillation frequencyf cn undamped natural frequency from continuous time domain[P ] Participation matrix[φ] Right eigenvector[ϕ] Left eigenvectorfNy Nyquist frequencyxiiAcknowledgmentsWriting this thesis was a long process of hard work and dedication duringwhich many people have somehow left their imprint on my life. This ac-knowledgement tries to reflect my gratitude to some of them.This thesis could not have been written without the encouragement andsupport of my wife, Sandra Zappa-Hollman. Her infinite and unconditionallove, endless patience and support made easier even the most challengingtimes. Our daughters, Roc´ıo and Serena, are also part of this family effortof surviving in spite of having both parents writing doctoral dissertationsunder the same roof. Thank you Roc´ıo and Serena for reminding me withyour cheerfulness and love of what is really important in life.I cannot find words to fully capture how thankful I am to my parentsMar´ıa del Carmen and Esteban. Their support, teachings and care helpedme to reach this point in life. I couldn’t have done it without them. Thankyou dad for taking me during school holidays with you to the power substa-tions and Hydro plants, which sparked my passion for power systems. Thankyou mom for getting me my first 386-based PC in which I ran power systemsimulations 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 myM.A.Sc. and Ph.D. under their supervision.A special recognition goes to my sister Vero´nica, and my parents and brotherin law Lieschen, Hansi and Andre´s, who always believed in us and uncondi-tionally supported us all these years.I also want to extend my thanks to our friends Vero´nica, Mat´ıas, Lucas,Dante and Sof´ıa Salibia´n, who became our family in Canada and made ourdaily life a pleasant adventure; and to our good friends Paula Murcia, PabloRuiz, Jesu´s Calvin˜o-Fraga, Vero´nica, Javier and Mat´ıas Gazzarri, Lito andAlicia Caro, Jorma and Laura Neuvonen, Hazem and Samar Sabbagh, Gilland Jorge Palejko, Mar´ıa Isabel and Felipe Gajardo, Naoko and Je´re´mieSe´ror, Aya and Tomoaki Nakashima, and Daniel and Lise Lindenmeyer.Numerous people at UBC have provided me with support in different,xiiiAcknowledgmentsinvaluable ways: Dr. Jesu´s Calvin˜o-Fraga, Dr. Luis Linares, Dr. Carlos Ven-tura, Dr. Juri Jatskevich, Mazana Armstrong, Tom de Rybel and AlejandroCervantes-Larios, by being good friends, colleagues and sharing stimulatingconversations; Dr. Hermann Dommel, who on many occasions spent timeproviding helpful comments and suggestions; Dr. Michael Davies and Dr.Vijay Bhargava, who have given me their unconditional support, and lastbut not least, everybody in the UBC Power group.I especially want to thank Doris Metcalf, Cathleen Holtvogt, Anne Coatesand David Chu Chong for their continued professional assistance and friend-ship.This dissertation would not have been possible without the generous fi-nancial support provided by several sources, including The National Sciencesand Engineering Research Council of Canada through an ISP-2 ResearchGrant; Dr. Jose´ Mart´ı with appointments as Research Assistant; the Depart-ment of Electrical and Computer Engineering via part-time appointments asTeaching Assistant; Fundacio´n Repsol-YPF through the Estenssoro Award;The University of Bristish Columbia through the Kenneth George VansackerMemorial Prize; Consejo Nacional de Investigaciones Cient´ıficas y Te´cnicas;and Dr. Prabha Kundur and Dr. Ali Moshref from Powertech labs withpart-time research and consulting contracts. To all of them goes my appre-ciation.Finally, I would especially like to express my most sincere gratitude tomy mentor and supervisor, Dr. Jose´ Mart´ı. His constant availability, pa-tience, guidance, leadership and support, have encouraged me throughoutthese years. Furthermore, he has become a role model beyond the academicworld.Jorge Ariel HollmanxivTo Serena, Roc´ıo and SandraxvChapter 1Introduction1.1 MotivationSince Thomas Alva Edison set up the world’s first commercial electric powersystem in lower Manhattan in September 1882, power systems grids haveundergone a continuous and aggressive expansion which enhanced their op-eration complexity. Just taking into consideration the case of the UnitedStates-Canadian system, for instance, this amazing expansion of the gener-ation capacity grew from 171 GW in 1960 up to 943 GW by 20041. Duringthe late eighties and nineties, first in the United Kingdom, subsequently incountries such as Argentina, Chile, New Zealand, many European countries,Brazil, and finally in some states of the United States the power industry waschallenged to a new operational paradigm: the implementation of a dereg-ulated environment. This new operational paradigm does not accept thetraditional extra investment in network infrastructure to ensure conserva-tive secure operation margins. The traditional technical constraints togetherwith the newly introduced financial constraints have forced the utilities to1Information available from the Energy Information Administration web site,http://www.eia.doe.gov1Chapter 1. Introductionoperate the systems near their stability limits. The maximization of revenuecan only be achieved under the new deregulated scheme by not only knowingaccurately but also knowing moment by moment the actual safe operationallimits of the system. These two factors have created new interest in theanalysis of stability problems. The high degree of interconnection and therequirement of synchronicity among all the system components makes theoperation of power systems not an easy task. Furthermore, the unbundlingof generation, transmission and supply is less oriented towards the physicalnature of the synchronously interconnected power systems. In addition, thelatest trend for the new century power system layout is the fast developmentof distributed generation.It is easy to visualize that centralized control operation of the entire systemwill be impossible due to the size and complexity of the objective functionand the response time constraints. The United States-Canadian electric-ity system alone deals with more than 320,000 km of transmission lines at230 KV or higher level, 950 GW of generating capability distributed amongnearly 3,500 utility organizations. The August 2003 North American Eastcoast blackout which affected an area of 50 millon people with a loss of 62GW of load lends strong support to the validity of this statement.A centralized coordination paradigm will inevitably shift to a locally coor-dinated operation. The 21st century power system will make extensive useof distributed intelligent local processing nodes capable of making their owndecisions and able to interact with other related nodes and with centralized2Chapter 1. Introductionregional controls. Since all parts of a power network interact with each other,local control action may require extensive information regarding the externalsystem to which the node is connected. This information can be acquiredand transmitted to all the distributed nodes or, even more efficiently, it canbe simulated in real-time at each decision node. However, to implementsufficiently detailed external system representations at local nodes requiresconsiderable computational power. This computational power must be ca-pable of simulating external system behaviour in real-time at the local node.Consequently, the concept becomes practical only if this computational so-lution can be provided in a fast, compact, accurate and inexpensive manner.Clear statements calling for further research on the proposed topic can befound in the literature. According to S. Gissinger [1], “Power system opera-tion requires increasingly complex decision-making in order to find the rightcompromise between economy and security”. The importance of consideringtechnical issues as a variable of the economic equation is pointed out by M.Pereira [2]: “Today’s new environment brings new uncertainties and risks(...) Much talk about risk today is focused on financial risks but it is re-markable that many important risks are not measured in dollars”. As well,the relevance of an extensive and detailed representation of the power sys-tem for the stability analysis with modal analysis is addressed by P. Kundurin [3] “A mode of oscillation in one part of the system can excite units ina remote part due to mode coupling. Analysis requires detailed and samelevel of representation throughout the system”. The feasibility of distributed3Chapter 1. Introductioncomputation in power system stability studies is stressed by J.L. Jardim in[4]: “Parallel processing has been increasingly gaining acceptance as a wayof speeding up the solution of highly intensive computational problems. Itis particularly indicated for applications that exhibit a natural decouplingamong computation tasks”. In sum, the above citations reflect but just alimited portion of the innumerable calls for further exploration of very largepower systems stability analysis.In particular, a discrete time state space description of electrical networksfrom nodal equations is of fundamental interest for this research work. Someof these key research questions were inspired by Dr. Dommel in the TheoryBook [5], where he pondered “could such a closed-form solution be used inan EMTP?... where the method becomes almost unmanageable, or useless,is in networks with nonlinear elements”.Hence, it is clear from the bibliographical research and from the above in-troduced scenario that there are great opportunities and an imminent needfor the development of fast on-line power systems simulators applicationsseeking faster and more accurate analysis, especially those with integrationof stability analysis of power systems. This new brand of simulators withthe capability of identifying the natural areas of association of the powersystem components will allow for the implementation of a decentralized op-eration layout. Given the present development of The University of BritishColumbia’s (UBC) real-time power system simulator -Object Virtual Net-work Integrator- (OVNI) in the hardware and software, it is attractive to4Chapter 1. Introductionexplore, develop and implement an extension of its capabilities to the stabil-ity analysis of power systems.1.2 Research Background1.2.1 Off-Line Power System SimulationThe Electromagnetic Transients Program (EMTP) [6] has been the standardtool for digital off-line simulation of power system transients for more thanthirty years. The EMTP became an attractive alternative to the bulky andexpensive Transient Network Analyzers (TNA) due to its portability and re-duced cost. Much effort was invested during the last three decades to developaccurate and sophisticated models to represent the behaviour of the powersystem components. While off-line power system transients simulators suchas the EMTP[6], EMTDC [7] and Power System Blockset [8], successfullycope with the requirements of system planning and transient studies areas,the off-line simulation environment exhibits great limitations when sensitiveequipment is required to be tested in a closed loop with the modelled powersystem or when real-time decision is a must. This called for the developmentof a new breed of power systems simulators which could cope with these lim-itations. This new group is based on the concepts presented by H. Dommelin [6], [5], but it is capable of achieving real-time performance. Alternatively,some research groups developed hybrid solutions based on Digital SignalProcessing (DSP) [9] to cope with the requirements of real-time equipment5Chapter 1. Introductiontesting. Notwithstanding its limitations, the off-line simulation environmentkeeps improving and extending its capabilities of analysis of power systemsby incorporating new areas of knowledge, such as artificial intelligence, neu-ral networks, stochastic analysis, risk theory, stability analysis, sensitivityanalysis, and so forth.1.2.2 Real-Time Power System SimulationAfter almost a decade now, various research groups have been able to achievereal-time performance [9], [10], [11], [12]. The real-time research group of theUniversity of British Columbia has been continuously advancing the conceptof a full-size real-time power system simulator, which eventually resultedin the OVNI [13] power system simulator. The simulator’s architecture isbased on a cluster of PC workstations. Using partitioning techniques in thehardware implementation and network decoupling techniques in the softwaresolution has led to the development of a distributed solver architecture, wherea network of inexpensive workstations is used to achieve real-time simulationfor fast transients of large power systems. There are mainly five well-knownresearch groups working in the field of real-time power system simulatorsusing different approaches. Electricite´ de France (EDF) and Hydro Que-bec Research Institute (IREQ) together with TEQSIM and Mitsubishi havechosen the very expensive and not so portable supercomputer architecture,while Manitoba Research Group has developed a hybrid solution based ona convenient arrangement of DSPs. On the other hand, Mitsubishi started6Chapter 1. Introductionwith supercomputers and, after joint work with UBC, moved their researchto PC’s. UBC’s research group was the first to choose an inexpensive PCsolution and achieve real-time performance with this scheme.Recently [14] EDF presented results of a parallel approach based on sharedmemory architecture. This approach is still based on supercomputers andshows no extra gain in speed in proportion to the number of CPU’s in servicebecause the time needed for communication between the CPU’s increaseswith the degree of parallelization. Mitsubishi research group achieved aPC cluster system based on six interconnected machines using a MyricomMyrinet giga-bit network. This approach presents a poor performance forhard real-time simulation due to its round trip overhead [12], but it is promis-ing for soft real-time applications. In contrast, UBC’s real-time PC clusterarchitecture shows a constant communication time for each type of config-uration link between the subsystem nodes, independently of the number ofnode solvers included in the array.UBC’s Real-time Network Solver (RTNS) [15] created by J. Mart´ı andL. 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 implementationfor testing relays [16] developed by Jesu´s Calvin˜o; the Multi-layered tearingsolution software [10]; and the Real-Time Distributed Network Simulator,developed by J. Hollman (RTDNS) [17]. See Fig.1.1. These efforts are partof the OVNI project [18], [13], [19], [20], [21] to develop a very large size7Chapter 1. IntroductionEMTPRTNS(One-layer-tearing software)SRTNS for RelayTester(Single-PC simulator)MRTNS(Multi-layer-tearing software)RTDNS(PC-Cluster Simulator)OVNIHybrid SimulationEnvironmentHard or SoftReal-TimeOff-LineEnvironmentTransient and Stability AnalysisTransient AnalysisFigure 1.1: Relation among research projects at UBC power system groupreal-time power system simulator. By breaking the power network solu-tion into subparts while maintaining simultaneous solution of all the parts,the Real-Time Power Systems group at UBC has developed a PC-Clustersimulator made up of an extensible number of off-the-shelf inexpensive PCboards that is capable of simulating in real-time fast system phenomena forsystems of practically any size just by extending the number of PC boards8Chapter 1. Introductioninterconnected. RTDNS achieves one of the lowest latency2 network solu-tions available today. However, with the advent of the Multi Area TheveninEquivalent (MATE) [21] solution scheme, a new, more flexible hardware in-terconnect solution became necessary. Also, the current system is limitedin bandwidth and tied into the Intelligent Drive Electronics (IDE) hard-diskinterface. This technology is soon to disappear in favour of serial AdvancedTechnology Attachment (SATA). The new hardware system, called OVNI-NET [22], circumvents any of the normal PC standard ports and is builtdirectly on top of the industry-standard Peripheral Component Interconnect(PCI) bus. Also, the parallel connections between the nodes are replaced byhigh-speed 1.2 Gbps full-duplex serial links, allowing for a system of largerphysical size with a more manageable wiring solution. These enhancementsyield an even lower latency while significantly boosting the throughput. Theyalso 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 AnalysisAlthough sophisticated power systems stability analysis techniques are avail-able [23], [24], [25], [26], [27], [28], [29], a large number of simplificationsis a common characteristic among them in order to reduce the computa-tional burden, especially for the case of on-line stability analysis tools. As a2Latency in this context is the time interval between the instant at which an instructioncontrol unit issues a call for data and the instant at which the transfer of data is started.9Chapter 1. Introductionresult, this situation makes it very difficult to determine whether observedsystem instability is due to the system properties or to the assumptions inthe algorithms. Generally, recursive and extensive off-line simulations arerequired before the power system stability can be accurately assessed. Thelack of detailed representation of the system in fast on-line assessment alsobecomes a constraint at the moment of identifying corrective responses, sincenot all the variables of the system are considered at the same instant. Onthe other hand, while detailed system representation is achieved in off-lineanalysis programs, the inherent solution times do not permit to achieve real-time performance. Most of the currently available transient analysis softwareworks in the off-line environment. Nonetheless, some groups [30], [31], mo-tivated by the present industry transformation, have been able to reach softreal-time performance in the range of tens of seconds up to hundreds of min-utes, yet still based on a large number of previous off-line simulations in thetime domain. Some stability tools use an intermediate approach betweenthe phasorial and time domains. The so-called fast time domain simulationapproach, based only on time domain simulation, focuses on the evolutionof the system operating conditions driven by the subsequent slow dynam-ics while neglecting part of the fast dynamics present. In some cases theneglected dynamics, such as those associated one with high-voltage directcurrent (HVDC) and Flexible AC Transmission Systems (FACTS) devices,become fundamental factors of the possible corrective actions. Whereas thistype of tools is satisfactory for certain kinds of studies, such as mid-term10Chapter 1. Introductionvoltage stability, it is not capable of producing a correct stability assessmentof the analyzed power system in the case of small-signal instability.1.3 Research Claim and ContributionsIn light of the current need for fast and accurate analysis solution of powersystems, this thesis makes the following contribution: The description andimplementation of a new and original power system stability assessmentmethodology that identifies the system’s eigenvalue trajectories in a real-timeEMTP solution incorporating the effect of switching phenomena and nonlin-ear load behaviour.Some of the advantages of the new methodology can be summarized as fol-lows:• Trajectory tracking of nonlinear element eigenvalues moment by mo-ment.• Capability of identifying suitable network partitioning schemes for ap-plication of multi-step integration solution in a hybrid power systemsimulator environment.• Capability of identifying small-signal interactions between elements byanalyzing damping ratios of natural frequencies and participation in-dexes from discrete time eigenvalues.11Chapter 1. Introduction• Visualization of eigenvalues trajectories in discrete time for the purposeof assessing power system’s dynamic behaviour.• Automatic selection of discretization step from discrete time eigenvalueinformation.1.4 List of PublicationsThe following are publications in which the author of this Thesis is the prin-cipal author or co-author and which are in part related to the topic of thiswork.• T. De Rybel, J. Hollman, J. Mart´ı. OVNI-NET: A Flexible Cluster In-terconnect for the New OVNI Real-Time Simulator. Proceedings of the15th Power System Computation Conference. Lie´ge, Belgium, 2005.• J. Hollman, J. Mart´ı. Real-Time Network Simulation with PC-Clusters.IEEE Transactions on Power Systems,Vol. 18, No. 2, pp. 563-569,IEEE Transaction on Power Systems, 2003.• J. Mart´ı, L. Linares, J. Hollman, F. Moreira. OVNI: Integrated Soft-ware/Hardware Solution for Real-Time Simulation of Large Power Sys-tems. Proceedings of the 14th Power Systems Computation Conference.Seville, Spain, 7 journal pages, 2002.• F. Moreira, J. Hollman, L. Linares, J. Mart´ı. Network Decoupling byLatency Exploitation and Distributed Hardware Architecture. Proceed-12Chapter 1. Introductionings of the International Conference on Power System Transients. Riode Janeiro, Brazil, 2001.• J. Mart´ı, J. Hollman, J. Calvin˜o. Implementation of a Real-Time Dis-tributed Network Simulator with PC-Cluster. Proceedings of the Inter-national Conference on Parallel Computing in Electrical Engineering,Quebec, Canada, 2000.• J. Hollman. Real Time Distributed Network Simulation with PC-Clusters.M.A.Sc. Thesis. Vancouver, 2000.13Chapter 2Framework2.1 IntroductionThis chapter provides a summary of current representative state of the artanalysis tools in the power system stability environment, pointing out theirstrengths and limitations and the motivation for the development of newmethodologies based on the EMTP solution. We highlight the capabilitiesof the EMTP simulators to follow nonlinear systems dynamic trajectoriesincluding their instantaneous eigenvalue/eigenvector states.2.2 The ProblemMaintaining interconnected power systems within their operational limits hasbeen always challenging due to the highly dynamic and nonlinear nature ofits components, such as loads, generators and changing operation topologies.Historically, due to the complexity of the problem, the assessment of oper-ation limits has been conducted in an off-line environment. Due to thesecomputational limitations, the operational limits are computed well in ad-14Chapter 2. Frameworkvance of the time in which they are expected to occur. The off-line resultsare loaded into look-up tables that are accessed by the operators. Duringactual real-time operation, the system operators find the closest match be-tween the actual system condition and the study set loaded into the tables.The limits for that condition are then used as the actual operation limits ofthe system. The new operational paradigm introduced by de-regulation hasmoved system operation much closer to its limits and off-line conservativeand structured operation conditions are no longer considered adequate. Thetrend is now towards on-line assessment tools.The theoretical frame of traditional stability analysis presents an interest-ing paradox: it has a highly detailed representation of elements such asgenerator and load controllers, while at the same time it uses a very lowdetail representation of other system operating aspects and components, likeswitching operations and nonlinear loads. Differential equations are used forthe generators, while steady state phasor solutions are used for the electricalnetwork, including the large component of motor loads. In many situations,this traditional approach is insufficient. One example is the case of interareaoscillation, which is a particular case of small-signal stability. Small-signalstability is the ability of the power system to maintain synchronism undersmall disturbances. Such disturbances occur continuously in the system be-cause of small variations in loads and generation. Small-signal instability canbe of two forms: steady increase in generator rotor angle, or rotor oscillationsof increasing amplitude. With present power system configurations, small-15Chapter 2. Frameworksignal stability is largely a problem of insufficient damping of oscillations andit can be classified into the following modes: Local, Interarea, Control, andTorsional [3].Interarea modes are associated with the swinging of machines in one partof the system against machines in other parts. These electromechanicallycoupled elements can be geographically very far apart. The problem occurswhen the electromechanically closely coupled machines are interconnected byweak ties. Interarea oscillations can severely constrain power transfers acrosslarge distances and can lead to widespread system disturbances if cascadingoutages of transmission lines occur due to oscillatory power swings, like thoseduring the blackout in western North America on August 10, 1996. If therepresentation of the system is not detailed enough, those dominant modescan be overlooked. Here, dominance refers to largeness in modulus.For the case of interarea oscillations, the size of the computational problemand the aggregation degree of the reduced system are the main constraintsof the traditional stability analysis approach. Reduced-order techniques ofanalysis are based on eigenvalue analysis of only the dominant modes. Thisapproach provides excellent information of associated frequencies and damp-ing ratios of the dominant coupled modes; however, they are usually inten-sive and have to be used off-line. Modified Arnoldi [32], and spectral analysisbased on Prony’s algorithm [33], [34], are the preferred traditional techniquesfor interarea stability analysis of large power systems.16Chapter 2. Framework2.2.1 Traditional Analysis ToolsMethods such as Prony’s spectral analysis and Fast Time-Domain (FTD)simulations are well developed. These tools reconstruct the system dynamicsby using curve fitting techniques on quasi-steady state simulation results.For the calculation of eigenvalues the most used methods are the Arnoldimethod and the Modified Arnoldi Method (MAM). An overview of theserepresentative methods of power system analysis follows.Prony AnalysisProny analysis was developed by Gaspard Riche, Baron de Prony [35] in1795 to explain the expansion of various gases. Riche proposed fitting asum of exponentials to equally space data points, and extended the modelto interpolate at intermediate points. Prony analysis is both a signal analy-sis technique and a system identification method widely used in biomedicalmonitoring, geophysical analysis, speech processing and power system elec-tromechanical oscillation analysis. As compared to Fourier analysis, Pronyanalysis has the advantage of providing estimation of damping coefficients inaddition to frequency, phase and amplitude of the Fourier technique. Whilethe Fourier series fits a sum of undamped complex exponentials to the data,Prony analysis fits a sum of damped complex exponentials [36]. A brief de-scription of Prony analysis is provided next.Let us consider the case of a Linear Time Invariant (LTI) dynamic system17Chapter 2. Frameworksuch as the one illustrated in Fig. 2.1(LTI)x(t)u(t) y(t)Figure 2.1: Linear time invariant (LTI) dynamic systemHere,y(t) is the system’s outputx(t) is the state of the systemu(t) is the input to the systemThe evolution of the LTI system is expressed byx˙(t) = Ax(t) +Bu(t) (2.1)where A and B are constant matrices. If the input is removed (u(t) = 0) forall t, equation (2.1) can be rewritten asx˙(t) = Ax(t) (2.2)where A is a square matrix of size n x n whose eigenvalues are λi, lefteigenvectors are qi and right eigenvectors are pi. The solution to equation18Chapter 2. Framework(2.2) can be expressed as a sum of n components,x˙(t) =n∑i=1(qTi x0)pie(λit) (2.3)With u(t) = 0 for all t, our LTI system output equation is given byy(t) = Cx(t) (2.4)Prony analysis directly estimates the parameters of the eigenstructure de-scribed in equation (2.3) by fitting a sum of complex damped sinusoids toevenly spaced time sample values of the output, that is,y(t) =L∑i=1Aie(σit) cos(2pifit+ φi) (2.5)Where,Ai is the amplitude of the component iσi is the damping coefficient of component ifi is the frequency of component iL is the total number of damped exponential componentsy(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 time19Chapter 2. FrameworkUsing Euler’s theorem, the complex sinusoids can be represented as a sumof exponentials,cos(2pifit+ φi) =e(j2pifit+φi) + e−(j2pifit+φi)2=e(j2pifit)ejφi2+e−(j2pifit)e−jφi2(2.6)replacing (2.6) into (2.5) and letting t = kT the samples of y(t) are given byy[k] =L∑i=1Ci µki︸︷︷︸poles(2.7)Where,Ci =Ai2ejφi (2.8)µi = e(σi+j2pifi)T (2.9)Prony analysis computes Ci and µi in three steps [37], as follows:i. solve a linear prediction model, constructed from the observed data setii. find roots of the characteristic polynomial formed from the linear pre-diction coefficientsiii. solve the original set of linear equations to yield the estimates of theexponential amplitude and sinusoidal phaseWith Ci and µi known, the amplitude, frequency and phase damping coeffi-cients are computed from equations (2.8) and (2.9).20Chapter 2. FrameworkThe main limitation of identification techniques, such as Prony analysis, forreal-time applications is the prohibitively large computational time requiredfor the fitting process. The Prony algorithm involves the solution of anoverdetermined set of linear equations and roots of a high order polynomial,which are numerically intensive operations. In addition, identification tech-niques are sensitive to noise. It is known that Prony analysis behaves poorlywhen a signal is embedded in noise, leading to a significant misestimationof damping and frequency terms. Thus, Prony analysis requires care in thechoice of processing parameters [38]. Also Prony analysis requires that anoscillation should exist in the studied time interval. If not, the calculated re-sult contains an equivalent of a higher order eigenterm. This is critical for thecase of fast and highly nonlinear elements, in which the eigenvalue location ischanging over time. Notice that even in the case when the exact set of eigen-values is found, the Prony method is not capable of providing informationabout which of the network components are responsible for the particularpoorly damped modes. In contrast, the methodology presented in this workshows a direct association between network component and eigenmode.Fast Time Domain simulation and Transient Security AssessmentToolFast Time Domain (FTD) simulation [39] simplifies the model of a powersystem in an quasi-steady state simulation approach in which many of thefaster device dynamics are ignored and the simulation is focused on the ac-21Chapter 2. Frameworktion of slower devices such as ULTC’s, overexcitation limiters and relays.The FTD simulation technique was developed for fast dynamic voltage sta-bility assessment. Even though the method considers the slow devices of thesystem, it is not always suitable to identify small-signal instability dynamicsdue to the simplifications imposed in order to speed up the solution. As anexample, the illustrated two-area system in Fig. 2.2 is simulated using theindustry standard Transient Stability Assessment Tool (TSAT) [28] for thecase 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 usedin this simulation. The simulation is first performed using detailed time do-main analysis, see Fig. 2.3, and later using the FTD technique, see Fig. 2.4.The results show the 60Hz phasor value of the voltage at bus 8 for the samecontingency in both cases. The full time domain simulation shows that the12657 8 910 11 34G1G2 G4G3L7 L925 km 25 km10 km10 km110 km110 kmArea 1 Area 2400 MWFigure 2.2: Two-area 11 bus test systemsystem becomes unstable for the given conditions while the FTD techniquesimulation shows change in the phasor value of the voltage and therefore no22Chapter 2. FrameworkBus voltage magnitude (pu)Time in seconds0.00 3.87 7.74 11.61 15.48 19.350.000.280.560.841.121.40Figure 2.3: Phasor value of bus 8 voltage for the two-area test system, small-signal instability observed with detailed time domain simulationoscillation after the contingency is applied and cleared. The FTD techniqueis not capable of capturing this case of small-signal instability.In very large systems, the low frequency modes are the only ones capableof coupling through very large distances to create instabilities. The higherfrequency modes are naturally more damped and they can not propagatefar away from the source; consequently, they do not contribute to interareaoscillations. For instance, a complete eigenvalue analysis of the two areas ofFig. 2.2, using the Prony technique (Table 2.1), makes clear the presence oflocal and interarea modes. The term “local” is associated with the swingingof units at a particular generating station with respect to the rest of thesystem. The term “interarea” refers to the swinging of many units in one23Chapter 2. FrameworkBus voltage magnitude (pu)Time in seconds0.00 6.00 12.00 18.00 24.00 30.000.9700.9840.9981.0121.0261.040Figure 2.4: Phasor value of bus 8 voltage for the two-area test system, small-signal simulated with FTD (no oscillation are captured by this solution)part of the system against many units in another part of the system. Interareaoscillations are caused by two or more groups of closely coupled units beinginterconnected by weak (high impedance) ties. Using Prony’s method onthe voltage signals at the system generation nodes (1,2,3,4) all the existingmodes are fitted. As observed in Table 2.1, the local modes (except forthe 1.03 Hz mode at generators 3 and 4) have sufficient damping, but theinterarea mode has insufficient damping, thus leading to the small-signalinstability problem observed in the simulation. The location of the observedmodes suggests the corrective action, in this case increase the damping of the0.62Hz interarea mode at generators # 1, # 2, # 3, # 4. It is important toemphasize that system identification techniques such as Prony analysis are24Chapter 2. FrameworkTable 2.1: Relevant modes for the two-area system caseFreq[Hz]Damp % Real [1/s] Imag [rad] Type of observed modeNode 1 0.62 -2.144 0.0831 3.875 InterareaGen.# 1 0.97 6.702 -0.4103 6.109 Local between G1 & G21.2 -0.579 0.0459 7.942 Local between G1 & G2Node 2 0.62 -1.077 0.0421 3.910 InterareaGen.# 2 0.97 8.926 -0.5432 6.061 Local between G1 & G21.2 1.052 0.0421 6.061 Local between G1 & G2Node 3 0.62 -2.357 0.0925 3.924 InterareaGen.# 3 1.03 -1.884 0.1215 6.446 Local between G3 & G41.7 13.677 -0.1087 10.743 Local between G3 & G4Node 4 0.62 -1.344 0.0527 3.921 InterareaGen.# 4 1.03 3.442 -0.2335 6.820 Local between G3 & G41.7 -0.675 0.0734 10.870 Local between G3 & G4not able to provide a direct association between network components andmodes with poor damping.Arnoldi MethodThe Arnoldi method (1951) [32] was proposed as a means of reducing a ma-trix to the Hessenberg form. Arnoldi already suggested that the process couldalso provide good approximations to some eigenvalues if stopped before itscompletion. Later this concept was fully developed by Saad (1980) [40] mak-ing the Arnoldi method evolve into one of the most successful ones for non-Hermitian eigenproblems. This is in contrast to the so-called direct methodswhich have to be complete in order to provide any useful results. The Arnoldimethod has poor numerical stability due to extensive numerical cancellation,25Chapter 2. Frameworkas pointed out by Wilkinson [41], and requires re-orthogonalization in or-der to improve its numerical stability performance. However, the inclusionof re-orthogonalization makes the method computationally less efficient [41].The so-called Modified Arnoldi method (MAM) by EPRI [42] is probablythe most popular method for identifying the set of dominant eigenvalues ofpower systems. The Arnoldi algorithm can be summarized as follows: Givena square matrix A of order n, if n steps of Arnoldi’s method is carried outthen an orthogonal reduction to Hessenberg form is achieved,AV = V H (2.10)WhereH is the upper Hessenberg matrix of order n with positive subdiagonalelements and V is an orthogonal matrix. Matrices V and H are uniquelydetermined by the first column of V , a unit-norm vector that is called theinitial vector. When an initial vector fails afterm steps, the following relationholds:AVm − VmHm = fe∗m (2.11)Where the vector f is called the residual of the m-step Arnoldi factorization.The input required by the Arnoldi algorithm are the matrix A, the number ofsteps m and the initial vector v1 of norm 1. The method provides as outputthe orthogonal matrices Vm, Hm, the residual f and β, so that verifies theequation (2.11) for β = ‖f‖2.In practical situations, the number of steps m required to achieve a good26Chapter 2. Frameworkapproximation may be too large. Consequently, the method becomes expen-sive due to storage and computational load that grows in every step. Thisis critical for the case of real-time applications. One possible strategy toovercome this problem is the concept of restarting the algorithm after msteps trying not only to cope with the growth of computational load, butalso benefitting from the previously computed information to define the newinitial vector. The ARPACK library provides a Fortran implementation ofthe implicit restarted Arnoldi method, for both real and complex arithmetic.2.2.2 Advantages of Step by Step Trajectory AnalysisA simple but illustrative general non-power system scenario can be used toshow the advantage of discrete time trajectory analysis compared to steadystate analysis. Take the case of a marathon runner who has a race tar-get time of 3h 08min 00 s. Based on a steady state performance, in orderto achieve the desired time, the runner would have to be able to run the42.195km with a pace of 4min 27 s/km. For the sake of producing a moreaccurate estimation of the required race pace, a fade of 10%/km based onthe history of the runner will be considered and thus the runner would haveto cope with an initial pace of 4min02 s/km. This initial pace will graduallyincrease 10%/km until the finish line. However, on the race day, the runnerwill most probably find out that because of the highly nonlinear dynamics ofthe race, such as the elevation of the course, the temperature, the wind, andthe physical condition of the runner during the race day, she or he will have a27Chapter 2. Frameworkdifferent pace performance, despite of long training hours. If the performanceis 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 giveenough information to the runner to adjust her or his pace to achieve thetarget time.By analogy, the same is valid for more complex and large systems such aspower grids. Solutions based on steady state or quasi-steady state computa-tions will be less accurate to represent highly nonlinear characteristics. Thisadditional accuracy is the main advantage of having analysis tools that arecapable of tracing point by point the transitions of nonlinear elements in thenetwork. The EMTP is a proven tool to provide this type of information[6], [5], [43], [44], [45], [46] and many steady state based power system stud-ies are verified against EMTP transient simulations. Providing the EMTPwith the ability to trace eigenvalue trajectories extends the capabilities ofthis highly accurate tool and the ability to represent highly sensitive systemoperating conditions. In addition, the methodology developed in this Thesisdirectly maps the eigenvalues to those system components responsible fortheir existence, therefore, allowing for the most effective corrective actions.28Chapter 3Eigenvalue Analysis withEMTP Solution3.1 IntroductionThis chapter starts with an overview of the EMTP solution, which aims toprovide the reader with sufficient background to reach a full understandingof the proposed new methodology. It is followed by a description of theproposed 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 ofassessing dynamic behaviour.3.2 EMTP SolutionThe EMTP solution, based on H.W. Dommel’s algorithm [6], is probably themost widely accepted computer program for power systems transient studies29Chapter 3. Eigenvalue Analysis with EMTP Solution[47]. In addition it is used as an educational tool in power system universityprograms [48], [49], [50].The beauty and efficiency of the EMTP solution are rooted in the nature ofthe discretization process. As clearly described in [19], the EMTP solutionfirst defines the discretization at the element level instead of discretizing thefull set of differential equations that describe the network behaviour. Con-sequently, in the EMTP every component relationship between its currentand voltage is modelled with a discrete time equivalent. These differenceequations are represented by discrete time equivalents circuits consisting ofresistance and source combinations. The values depend on the integrationrule applied.For instance, for the case of an Inductor (Fig. 3.1) the relation between theLFigure 3.1: Inductor elementvoltage across it and its current is given byVL(t) = LdiL(t)dt(3.1)30Chapter 3. Eigenvalue Analysis with EMTP SolutionSolving for the voltage,VLdt = LdiL(t) (3.2)∫ tt−∆tVLdt = L∫ tt−∆tdiL(t) (3.3)∫ tt−∆tVLdt = L[iL(t)− iL(t−∆t)] (3.4)where ∆t is the discretization time step. Approximating the left hand sidewith the trapezoidal formula we obtainVL(t) + VL(t−∆t)2∆t = L[iL(t)− iL(t−∆t)] (3.5)Then the voltage across the inductor isVL(t) =2L∆t︸︷︷︸RequiviL(t) + [−VL(t−∆t)− 2L∆tiL(t−∆t)︸ ︷︷ ︸history] (3.6)Equation 3.6 can be represented by the equivalent discrete time circuit ofFig. 3.2.The trapezoidal integration rule is used in the EMTP with good accuracy,as illustrated in Appendix A. This rule, however, introduces numerical oscil-lations at discontinuities. This situation is averted in the EMTP by usingthe Critical Damping Adjustment algorithm (CDA) developed by Mart´ı andLin [51]. As described by Linares in [19], in systems with a large componentof power electronic devices the backward Euler integration rule can be an31Chapter 3. Eigenvalue Analysis with EMTP Solutioneh(t)vL+ -iL  2L/   tFigure 3.2: Inductor discrete time equivalent with trapezoidal ruleeffective alternative for Real-time simulation.If a constant discretization time step is applied, the nodal analysis matrixconsists of a constant matrix of conductances and a current source vector.If a variable discretization time step is used, the conductance matrix has tobe recalculated whenever the time-step is changed. The conductance matrixalso needs to be reevaluated if there is a change in the network topology.This is the case, for instance, of circuits with nonlinear characteristics orswitching electronic devices.3.2.1 Nodal AnalysisNodal analysis is based on Kirchoff’s current law, which states that the al-gebraic sum of currents entering a node is zero, leading to the formulation ofa system with the form[Y ][V ] = [J ] (3.7)32Chapter 3. Eigenvalue Analysis with EMTP Solutionwhere,Y is the admittance matrix of order N x NV is the vector of N node voltages referred to groundJ is the vector of N current sources injected into the nodesN is the number of nodes, excluding groundFor each independent voltage source, the voltage of the node at which thesource is connected is known. If there are N2 sources, the number of nodeswith unknown voltages isN1 = N −N2 (3.8)which can be partitioned as follows,Y11 Y12Y21 Y22V1V2 =J1J2 (3.9)thenY11V1 = J1 − Y12V2 (3.10)where,V1 is the vector of N1 nodes with unknown voltages referred to groundV2 is the vector of N2 nodes with voltages sources connected between thenode and ground (known voltages)Y11, Y12, Y21, Y22 are admittances submatrices of order N1 x N1, N1 x N2,N2 x N1, N2 x N2, respectivelyJ1 is the vector of current sources entering the N1 unknown voltages nodes33Chapter 3. Eigenvalue Analysis with EMTP SolutionJ2 is the vector of current sources entering the N2 known voltages nodesEquation 3.10 constitutes a system of N1 equations that can be solved forV1 by any numerical method. Currents in the elements are then obtainedusing the original elements’s branch equations. The admittance matrix Y isdefined as follows:• Diagonal elements of the Y matrix come from the sum of all admit-tances at the corresponding node.• Off-diagonal elements of the Y matrix come from the negative of theadmittances connecting the corresponding nodes.3.2.2 Numerical Integration AccuracyWhen analyzing the total distortion error introduced by solving a systemwith a discrete time solution, both the highest frequency represented and thedistortion error of the integration rule have to be taken into consideration.The Nyquist frequency is the bandwidth of a discrete time signal, and it isequal to half of the sampling frequency of that signal. In EMTP discrete timesimulation, the sampled signal will contain frequencies ranging from 0Hz tothe Nyquist frequency fNy,fNy =fsample2=12∆t(3.11)The smaller the integration time step, the higher Nyquist frequency. There34Chapter 3. Eigenvalue Analysis with EMTP Solution0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.500.ño 3Backward & Forward Eulerf(pu) = fsim(Hz)/fNy(Hz)He/H, Magnitude distortion (p.u.) Figure 3.3: Magnitude frequency response of integration rulesis no aliasing in the EMTP solution because the input to the system, thesources, have values generated at each ∆t of the system solution (no higherfrequency values are sampled in). However, as the frequencies in the net-work solution approach the Nyquist frequency, the accuracy of the solutiondecreases (“distortion error”). The distortion error introduced by the cho-sen integration rule increases as the frequency gets closer to the Nyquist fre-quency. Figures 3.3 and 3.4 show the magnitude and phase distortion [He/H],[ΦHe−Φexact] for four discretization rules. The magnitude distortion is shown35Chapter 3. Eigenvalue Analysis with EMTP Solutionin 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 curvesis shown in per unit of the Nyquist frequency, i.e. f(pu) = f(Hz) / fNy(Hz).For a given frequency in the simulation, the error will decrease if the time0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5−100−80−60−40−20020406080100TrapezoidalCalviño 3Backward EulerForward Eulerf(pu) = fsim(Hz)/fNy(Hz)ΦHe − Φexact    Phase distortion  [degrees]Figure 3.4: Phase frequency response of integration rulesstep is reduced (i.e., if the Nyquist frequency is “pushed away”). In powersystem transients simulation it is desirable to obtain an accurate represen-tation of high frequencies when fast transients, fast switching operations,HVDC converters, and TVR studies occur. In real-time simulation when the36Chapter 3. Eigenvalue Analysis with EMTP Solutioncomputer is synchronized to external hardware requiring a solution point ev-ery ∆t time units, there is a limit in the system size that can be solved suchthat the computer processing time plus the overhead to put each sample outinto the external hardware is less than that ∆t. A solution to this problem isto implement the simulation with a distributed multi-pc’s architecture, suchas the ones proposed and implemented in [17],[52],[22].3.2.3 Discretization Rule StabilityThe stability of the integration rule can be assessed through the TransferFunction in the Z-domain [53]. For the particular case of an Inductancediscretized with trapezoidal rule,i(t)− i(t−∆t) = ∆t2L+∆t2Lv(t−∆t) (3.12)Applying the Z-transform,I(z)− z−1I(z) = ∆t2LV (z) +∆t2Lz−1V (z) (3.13)Ye(z)Trap =I(z)V (z)=∆t2L(z + 1)(z − 1) (3.14)Trapezoidal with a p1 = 1 and z1 = −1 has the characteristic of a stableintegrator and a stable differentiator, with bounded oscillations at disconti-nuities for the latter case.37Chapter 3. Eigenvalue Analysis with EMTP SolutionSimilarly, for an inductor discretized with backward Euler the Z-transferfunction is given byYe(z)BE =∆tLz(z − 1) (3.15)Backward Euler with a p1 = 1 and a z1 = 0 has the behavior of a stableintegrator and a stable differentiator.Similarly, for the case of Calvin˜o 3 [54], which has a stable integrator anddifferentiator characteristic, we get the transfer function,Ye(z)C3 =∆tk0L1[1− k1z−1 + k2z−2 − k3z−3] (3.16)Figures 3.3 and 3.4 show the distortion introduced by the discretization ruleas 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 discretizedinductor with trapezoidal rule,Ye(w) =∆t2L(ejω∆t + 1)(ejω∆t − 1) (3.17)Similarly, for backward Euler in (3.15) we obtainYe(w)BE =∆tLejω∆t(ejω∆t − 1) (3.18)38Chapter 3. Eigenvalue Analysis with EMTP Solutionand for Calvin˜o 3 in (3.16) we obtainYe(w)C3 =∆tk0L1[1− k1e−jω∆t + k2e−jω2∆t − k3e−jω3∆t] (3.19)Notice from these figures that trapezoidal does not introduce phase distortionwhereas backward Euler introduces a strong phase distortion. Trapezoidalcan have stability problems at discontinuities, this is overcome by applyingthe concept of critical damping adjustment (CDA) [51]. Backward Euler hasno problems at discontinuities. The behaviour of trapezoidal, backward Eulerand other rules at discontinuities is analyzed in [51], where CDA combiningthe best characteristics of trapezoidal and backward Euler rules is introduced.3.5, as follows.internalstatesu1(t)u2(t)u2(t)um(t)y1(t)y2(t)y2(t)yp(t)x1(t),x2(t),x2(t),...xn(t)Figure 3.5: State space representation of a multiple input/output system39Chapter 3. Eigenvalue Analysis with EMTP Solution3.3 Discrete Time State Space Formulationof an Electrical Network from theEMTP Nodal Matrix FormulationThe concept of representing a discrete time system by a set of first orderdifference equations became a standard tool for the research engineer by theend of the 1950s. These techniques have since become generally known asstate space representations. Any Linear Time Invariant (LTI) system canbe described by a set of discrete time state space equations. A state spacerepresentation utilizes three types of variables to model the dynamics ofa system: the input, which are the external forces driving the system; theoutput, which are the quantities directly measurable by an external observer;and the states, which characterize the internal dynamics of the system. A setof state variables can be computed as a function of the present and past statesand the present and past system inputs. The definition implies that if weknow the state at time t we can then compute the energy stored in the systemat that instant. In a general framework, the state variables can be chosenas an arbitrary combination of inner system variables. This generalizationbuilds some distance between the state and its physical interpretation, but itmakes more evident that the choice of state variables is not unique. A statespace representation of a system with multiple inputs and outputs can bevisualized in a block diagram, see Fig. Here we assume m inputs, p outputs40Chapter 3. Eigenvalue Analysis with EMTP Solutionand n states. A system governed by n linear difference equations of orderone requires n internal states for a minimal state space representation. Usingforward Euler discretization, the values of the state variables at discrete timepoint (t) can be expressed as a linear combination of values of the state andinput variables at the previous discrete time point (t−∆t),x1(t)...xn(t) =a11 . . . a1n.... . ....an1 . . . annx1(t−∆t)...xn(t−∆t)+b11 . . . b1m.... . ....an1 . . . anmu1(t−∆t)...um(t−∆t)(3.20)while the output variables at discrete time (t) can also be expressed as linearcombinations of the values of the input and state variables at discrete time(t).y1(t)...yp(t) =c11 . . . c1n.... . ....cp1 . . . cpnx1(t)...xn(t)+d11 . . . d1m.... . ....dp1 . . . dpmu1(t)...um(t)(3.21)41Chapter 3. Eigenvalue Analysis with EMTP SolutionEven though, as discussed in the previous chapter, the forward Euler rule isnot stable for EMTP simulation, we will show subsequently that when adapt-ing the EMTP solution to state space form, for the chosen EMTP variables,forms (3.20) and (3.21) are obtained although the network components arediscretized using stable rules such as trapezoidal and backward Euler. For asingle input-output system (see figure 3.6), (3.20) and (3.21) becamex(t) = A x(t−∆t) +B u(t−∆t) (3.22)y(t) = C x(t) +D u(t) (3.23)System states change from moment to moment, and what defines the cur-rent value of any state of the system are the memory, current input and theconnectivity topology. Changes of the states can be originated by the inputor by the system’s dynamics. Changes due to input are forced, while changesdue to system dynamics reflect the nature of the system.3.4 EMTP Solution in State Space formH. Dommel in [6] proposed the following nodal solution form for the Electro-magnetic Transients Program (EMTP),[G] [v(t)]=[h(t)]+[u(t)](3.24)42Chapter 3. Eigenvalue Analysis with EMTP Solutionu(t)DBACZ-1y(t)x(t+ t) x(t)Figure 3.6: State space system diagram blockWhere [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 historysources. The system “output” is given by[v(t)]=[G]−1 [h(t)]+[G]−1 [u(t)](3.25)The updating formula for the history sources can be expressed as[h(t+∆t)]=[A] [h(t)]+[B] [u(t)](3.26)The formula (3.26) establishes the connection between the history sourceconcept of the EMTP solution and the previous state concept of the statespace 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 state43Chapter 3. Eigenvalue Analysis with EMTP Solutionspace 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 heredefines the branch history terms as states of the electrical network, but itretrieves all the necessary information from the EMTP nodal solution andthe case data input.If there is no input applied to the system, the changes in [h(t)] constitutethe homogeneous solution of the system. By excluding all external sourceswe can compute the transition matrix of the system due to the nature of thesystem. Finally, from [A] we can study the eigenvalues of the system andanalyze its stability behavior as will be shown in detail subsequently.3.4.1 Discrete Time State Space Equation for anInductanceConsider a single inductance discretized using the Trapezoidal integrationrule, (Fig. 3.7). This basic circuit element is used to model single-phase shuntreactors, discharge circuits, smoothing reactors in HVDC converter stations,the inductive part of Thevenin/Helmholtz equivalent circuits, the inductivepart of single-phase nominal pi circuits of balanced systems operation andeven the inductive equivalent part for loads of low frequencies.The voltage across the element isvkm(t) = vk(t)− vm(t) (3.27)44Chapter 3. Eigenvalue Analysis with EMTP Solutionig =   t / 2Lk mh(t)vkmk mLFigure 3.7: Discrete time equivalent of an Inductorand 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 astate variable.The current through the inductance isikm(t) = gvkm(t)− hkm(t) (3.29)while the current ikm(t−∆t) isikm(t−∆t) = gvkm(t−∆t)− hkm(t−∆t) (3.30)45Chapter 3. Eigenvalue Analysis with EMTP SolutionBy replacing (3.30) into (3.28), we obtain the state space equation of a selfinductance,hkm(t) = −gvkm(t−∆t)− gvkm(t−∆t) + hkm(t−∆t) (3.31)Rearranging,hkm(t)︸ ︷︷ ︸new state= hkm(t−∆t)︸ ︷︷ ︸old state+(−2g)vkm(t−∆t)︸ ︷︷ ︸forcing function(3.32)Compare this formula with (3.22) in the state space description. As men-tioned earlier, even though in state space (3.22) corresponds to forward Eulerdiscretization, in the EMTP (3.32) was obtained from trapezoidal discretiza-tion. The one-element state equation for the inductance can be extended tobranches connected in a network. The branch histories and branch voltagescan be easily described by nodal equations using the network connectivity in-formation contained in the incidence matrix. In this case the forcing functionv will be defined by the network solution.3.4.2 Discrete Time State Space Equation of aCapacitanceA single capacitance discretized using the trapezoidal integration rule isshown in Fig. 3.8. This model is typically used to represent series andshunt capacitors, shunt capacitances in nominal pi-circuit representations of46Chapter 3. Eigenvalue Analysis with EMTP Solutionig = 2C /   tk mh(t)vkmk mCFigure 3.8: Discrete equivalent of a Capacitortransmission lines, equipment in HVDC converter stations such as surge ca-pacitors and filters, capacitance potential transformers, capacitive voltagedividers, and stray capacitances of transformers, generators, especially forthe case of transient recovery voltage and lightning studies, among others.Again, we start with the update formula from the nodal equationhkm(t) = gvkm(t−∆t) + ikm(t−∆t) (3.33)The current through the capacitance isikm(t) = gvkm(t)− hkm(t) (3.34)while the current for (t−∆t), analogously to the case of the inductance is,ikm(t−∆t) = gvkm(t−∆t)− hkm(t−∆t) (3.35)47Chapter 3. Eigenvalue Analysis with EMTP SolutionBy replacing (3.35) into (3.33), the state space description of a single capac-itance is obtained,hkm(t) = gvkm(t−∆t) + gvkm(t−∆t)− hkm(t−∆t) (3.36)Rearranging,hkm(t)︸ ︷︷ ︸new state= −hkm(t−∆t)︸ ︷︷ ︸old state+2g vkm(t−∆t)︸ ︷︷ ︸forcing function(3.37)Which has the same form as (3.22) for the state space description. Analo-gously to the inductance case, the one-element state equation for the capac-itance can be extended to branches connected in a network, in which casethe forcing function v will be defined by the network solution.3.4.3 Discrete Time State Space Equation of aResistanceIn this case, Fig. 3.9, hkm(t) = 0 for all t. Even though resistive elements donot contribute with a state because they do not have branch history, theydo provide damping. To avoid a special case in the state space formulationthe R element can be discretized together with either an inductance or acapacitance branch as described in the following section.48Chapter 3. Eigenvalue Analysis with EMTP Solutiong = 1/Rk mh(t) = 0k mRFigure 3.9: Discrete equivalent of a Resistor3.4.4 Treatment of Series Connection RL and RCThe series connection of a resistor with an inductor (RL) or a resistor witha capacitor (RC) can be represented with a combined equivalent discrete el-ement. The combined element approach not only takes care of the dampingintroduced by the resistive element, but also improves the performance ofthe overall solution, since it reduces the number of nodes to solve and, con-sequently, the number of branches. For the case of an inductor in series witha resistor, the aggregation process is shown in Fig. 3.10.The general form of equation (3.32) is kept,hkm(t) = hkm(t−∆t)− 2gRLvkm(t−∆t) (3.38)where,gRL =1R + 2L∆t(3.39)49Chapter 3. Eigenvalue Analysis with EMTP Solutioneh(t)vL+ -iL  2L/   tgRLh(t)vRLiLLR RRL discrete equivalentk mk meh(t)vRL+ -k m RRL = R + 2L/   tFigure 3.10: Discrete equivalent of a Series RLAnalogously for the series RC,hkm(t) = −hkm(t−∆t) + 2gRCvkm(t−∆t) (3.40)gRC =1R + ∆t2C(3.41)3.4.5 Treatment of Series BranchesInitially we will consider pure L or pure C branches, see Fig. 3.11, to illustratethe concept.50Chapter 3. Eigenvalue Analysis with EMTP SolutionhL1hC2hL321gC2gL1gL3  L3  C2  L11 2Branch  3Branch  1Branch  2Figure 3.11: Generic circuitThe nodal equations for the generic circuit of Fig. 3.11 are[G][v] = [j] (3.42)G11 G12G21 G22v1v2 =j1j2In detail:gL1 + gL3 −gL3−gL3 gC2 + gL3v1(t)v2(t) =hL1(t) + hL3(t)hC2(t)− hL3(t)v1(t)v2(t) = [G−1]j1(t)j2(t)51Chapter 3. Eigenvalue Analysis with EMTP Solutionv1(t)v2(t) =R11 R12R21 R22j1(t)j2(t) (3.43)In the circuit of Fig. 3.11 we want to identify the transition matrix [A]dby 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 thehistory of the branch i is denoted as_hi. We know from equations (3.32) and(3.37) that for each branch element L or C, the branch histories are:_h1(t)_h2(t)_h3(t) =k∧1 0 00 k∧2 00 0 k∧3_h1(t−∆t)_h2(t−∆t)_h3(t−∆t)+g∧1 0 00 g∧2 00 0 g∧3_v 1(t−∆t)_v 2(t−∆t)_v 3(t−∆t)(3.44)In particular for our circuit,k∧1 = 1 g∧1 = −∆tL1k∧2 = −1 g∧2 = 4C2∆tk∧3 = 1 g∧3 = −∆tL3(3.45)Notice that k∧ = 0 and g∧ = 0 for the case of a purely resistive branch.52Chapter 3. Eigenvalue Analysis with EMTP SolutionThe relationship between branch and node voltages is given by the branch-node incidence matrix:_v 1(t−∆t)_v 2(t−∆t)_v 3(t−∆t) =1 00 11 −1v1(t−∆t)v2(t−∆t) (3.46)Substituting (3.43) in (3.46),_v 1(t−∆t)_v 2(t−∆t)_v 3(t−∆t) =1 00 11 −1R11 R12R21 R22j1(t−∆t)j2(t−∆t) (3.47)The nodal history sources in terms of the branch histories are given byj1(t−∆t)j2(t−∆t) =_h1(t−∆t) + _h3(t−∆t)_h2(t−∆t)−_h3(t−∆t)=1 0 10 1 −1_h1(t−∆t)_h2(t−∆t)_h3(t−∆t) (3.48)53Chapter 3. Eigenvalue Analysis with EMTP Solutionsubstituting (3.47), (3.48) in (3.44),_h1(t)_h2(t)_h3(t) =k∧1 0 00 k∧2 00 0 k∧3︸ ︷︷ ︸[k∧]_h1(t−∆t)_h2(t−∆t)_h3(t−∆t)+g∧1 0 00 g∧2 00 0 g∧3︸ ︷︷ ︸[g∧]1 00 11 −1︸ ︷︷ ︸[L]R11 R12R21 R22︸ ︷︷ ︸[G−1]1 0 10 1 −1︸ ︷︷ ︸[Lt]_h1(t−∆t)_h2(t−∆t)_h3(t−∆t)[_h(t)]=[k∧] [_h(t−∆t)]+[g∧] [L] [G]−1 [L]t [_h(t−∆t)](3.49)The discrete time transition matrix [A]d is then given by[A]d=[k∧]+[g∧] [L] [G−1] [Lt](3.50)which gives the connection between the [G] matrix based EMTP solution andthe [A] based state space description. This connection had not been proposedpreviuosly to this work. Some interesting properties can be observed:54Chapter 3. Eigenvalue Analysis with EMTP Solution• [A]d depends directly on [G−1]• [L] is the graph connectivity• [k∧] is a diagonal matrix that indicates the nature of the branch, i.e.,+1 for an L, −1 for a C or 0 for an R.• [g∧] is the diagonal matrix of the branch equivalent parameters as de-termined by the discretization rule and the time step ∆t.3.4.6 Treatment of Parallel BranchesFor the case of L and C elements in parallel, see Fig. 3.12, we can choosefrom two possible approaches. First we can keep the identity of each compo-nent and by augmenting the size of the [A]d matrix compute each element’seigenvalue. This alternative is convenient for the case in which it is neces-sary to identify, analyze and specify the corrective action for a single elementcontributing to a specific instability.Alternatively, we can produce an aggregation of elements of the same typeto obtain an equivalent Req, Leq and Ceq, see Fig. 3.12. By doing so we aredefining a new state and are able to compute its eigenvalue, which reflectsthe dynamic characteristic of the aggregate elements.Computationally, either option still has only one matrix inversion [G−1]. Theonly difference is the size of the resulting [A]d matrix.55Chapter 3. Eigenvalue Analysis with EMTP Solution2L1C21L3C4 L5 CeqLeq21Figure 3.12: Aggregate treatment of parallel RLC branches_h1(t)_h2(t)_h3(t)_h4(t)_h5(t)=k∧1 0 0 0 00 k∧2 0 0 00 0 k∧3 0 00 0 0 k∧4 00 0 0 0 k∧5_h1(t−∆t)_h2(t−∆t)_h3(t−∆t)_h4(t−∆t)_h5(t−∆t)+g∧1 0 0 0 00 g∧2 0 0 00 0 g∧3 0 00 0 0 g∧4 00 0 0 0 g∧5_v 1(t−∆t)_v 2(t−∆t)_v 3(t−∆t)_v 4(t−∆t)_v 5(t−∆t)(3.51)56Chapter 3. Eigenvalue Analysis with EMTP Solution_v 1(t−∆t)_v 2(t−∆t)_v 3(t−∆t)_v 4(t−∆t)_v 5(t−∆t)=1 −11 −11 −11 −11 −1v1(t−∆t)v2(t−∆t) (3.52)The general formulas are maintained and finally we can obtain [A]d withequation (3.50), as for the case of a single branch per node case,[A]d=[k∧]+[g∧] [L] [G−1] [Lt]3.5 Nonlinear NetworksOver the years the EMTP has used three different approaches to handlenonlinear elements:• current-source representation with time delay• compensation method• piecewise linear representationThe first one is no longer used, while the second and third ones are currentlyused. In this work the piecewise linear representation is chosen.In this representation, nonlinear elements are considered as made up ofpiecewise linear segments [5]. Depending on the particular component, a57Chapter 3. Eigenvalue Analysis with EMTP Solutionλ Saturationλ Saturationtan α2 = L2tan α1 = L1Figure 3.13: Piecewise linear inductance with two slopespiecewise segment can be relatively long (like the saturation characteristic ofpower transformers) or it can be as short as one segment per solution timestep. A change of piecewise segment in the EMTP solution corresponds toa change in the value of the component (around the particular operatingpoint), and therefore, to a new set of eigenvalue/eigenvectors. The dynamicsof these eigensets reflect the system behaviour in terms of trajectories (wherethe 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 L1and Lp in parallel (3.53), as shown in Fig. 3.14. The flux in Lp is always58Chapter 3. Eigenvalue Analysis with EMTP Solutioncomputed by integrating the voltage across the branch independent of theswitch position. In this representation a change in the operating region cor-responds to a change in the topology of the network. During the simulationLp will be connected whenever |λ| > λSaturation and disconnected as soon as|λ| < λSaturation. This is accurately enough in most cases to represent thesaturation characteristics of modern transformers.1L2=1L1+1Lp(3.53)LpL1λ1λpFigure 3.14: Implementation of two slope piecewise linear inductance3.6 Continuous Time Eigenvalues fromDiscrete Time Domain SolutionWhen eigenvalues are computed from the EMTP solution, the matrices in-volved contain only real numbers. Eigenvalue extraction routines, like the59Chapter 3. Eigenvalue Analysis with EMTP SolutionQR method (Appendix A) are much more efficient for real matrices than forcomplex matrices. The process of extracting eigenvalues from the discretetime EMTP solution and reconstructing the continuous time eigenvalues fromthe inverse discretization mapping (3.64) can actually be much more efficient,particularly for large systems, than the direct computation of complex eigen-values from steady state power system admittance matrices whose elementsare complex numbers.Notice also that in the important case of eigenvalue calculations with non-linearities in the circuit, the eigenvalues will change incrementally from onesolution step to the next. Using the previously calculated eigenvalue setas the starting point for the new set can also speed up considerably theeigenvalue calculation. This situation is particularly advantageous for theRestarted Arnoldi method.3.6.1 Discrete to Continuous Time MappingThe presented methodology allows for computing the discrete time eigen-values and eigenvectors of a given power system described with the EMTPsolution. Even though the eigenvalues obtained in the discrete time solutionare affected by the inaccuracies of the distortion introduced by the discretiza-tion rule (Fig.s 3.3 and 3.4), once the discrete time eigenvalues are obtained,it is possible to reconstruct the exact eigenvalues and eigenvectors of the orig-inal continuous time system. Similarly, the exact continuous time transitionmatrix can be built.60Chapter 3. Eigenvalue Analysis with EMTP SolutionConsider the homogeneous part of the continuous time state space equations,x˙(t) = Acx(t) (3.54)and its corresponding discrete time state space equivalent given by 3.22,x(t+∆t) = Adx(t) (3.55)Discretizing 3.54 with the trapezoidal integration rule,x˙(t+∆t) =2∆tx(t+∆t)− 2∆th(t+∆t) (3.56)whereh(t+∆t) = x(t) +∆t2x˙(t) (3.57)substituting 3.57 in 3.56,x˙(t+∆t) =2∆t{x(t+∆t)− x(t)} − x˙(t) (3.58)where x˙(t) = Acx(t) and x˙(t+∆t) = Acx(t+∆t),Acx(t+∆t) =2∆t{x(t+∆t)− x(t)} − Acx(t) (3.59)61Chapter 3. Eigenvalue Analysis with EMTP Solutionsolving for x(t+∆t),x(t+∆t) = x(t)2∆tI − Ac2∆tI + Ac(3.60)since,x(t+∆t) = Adx(t) (3.61)a relationship between Ad and Ac for the case of discretizing with Trapezoidal,can be established asAd =ΩI + AcΩI − Ac (3.62)where Ω = 2/∆t.The eigenvalues and eigenvectors of Ac are defined byAcvi = λivi (3.63)The continuous time eigenvalues λi of Ad can be computed immediately fromthose discrete time eigenvalues zi of Ad, and vice-versa, as follows,λi = Ωzi − 1zi + 1(3.64)zi = Ω1 + λi1− λi (3.65)Equation (3.65) corresponds to the classical Bilinear transformation in signalprocessing theory. In the case of the trapezoidal rule (Bilinear transforma-62Chapter 3. Eigenvalue Analysis with EMTP Solutiontion), the eigenvalues located anywhere in the entire left-hand side of thecomplex plane in the continuous time domain maps into the entire unit cir-cle in the discrete time domain as shown in Fig. 3.15.In the case when backward Euler is used as the integration rule, the relationbetween continuous and discrete time eigenvalues, as depicted in Fig. 3.16,can be expressed asλi = Ωzi − 1zi(3.66)In the case of backward Euler the left hand side of the complex plane ismapped into a subregion of the unit circle in the discrete time domain. Ingeneral, for a stable discretization rule the left half of the complex plane inthe continuous time maps into a subarea inside the unit circle in the discretetime domain.3.6.2 Time Discretization ConsiderationsThe linearization of differential equations into in difference equations is anapproximation that will produce accurate results when the highest frequencypresent in the system is 5 to 10 times below the Nyquist frequency, as dis-cussed in section 3.2.2. However, when using the discrete time formulationof a network to compute the eigenvalues some extra considerations have tobe in place.As discussed in section 3.6.1, the accuracy of the discrete to continuoustime mapping in LTI systems will be correct until the Nyquist frequency.63Chapter 3. Eigenvalue Analysis with EMTP SolutionReal RealImaginary ImaginaryλλzzComplex continuous-time plane                                    Complex discrete-time planeUnitcircleFigure 3.15: Continuous time and discrete time stability mapping for trape-zoidalReal RealImaginary ImaginaryλλzzComplex continuous-time plane                                    Complex discrete-time planeUnitcircleFigure 3.16: Continuous time and discrete time stability mapping for back-ward Euler64Chapter 3. Eigenvalue Analysis with EMTP SolutionHowever, in the case of systems with nonlinear elements the discrete to con-tinuous time mapping will be accurate only if the discretization time is smallenough to represent correctly the region change of the nonlinear elements.In addition to the limitations imposed by nonlinear elements, using smallerdiscretization times than required will not lead to further improvements inthe accuracy of the continuous time eigenvalues. However, a nonlinear scalingeffect will be present in the discrete time values due to the mapping charac-teristics of the discretization transformation. The smaller the discretizationtime, 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 aneigenvalue λ is given by eλt. Consequently,• A purely real eigenvalue corresponds to a non-oscillatory mode.– A positive real eigenvalue represents an exponential growing in-stability.– A negative real eigenvalue represents an exponential decaying mode.– A zero eigenvalue represents a constant mode.• Complex eigenvalues occur in conjugate pairs, and each pair corre-sponds to an oscillatory mode with an exponential amplitude with asbefore can be decaying in the stable case or exponentially growing inthe unstable case. In the critically damped case the amplitude becomesconstant.65Chapter 4Application to the OVNIsimulator4.1 IntroductionThis chapter starts with an overview of the OVNI simulator and its associatedMATE solution framework. Also a background of OVNI’s architecture forreal-time simulation is provided. Next follows a description of the procedureto obtain the state space formulation to compute the discrete time eigenvaluesfor large networks from MATE. The chapter concludes with two extendedapplications of the new methodology for the Latency concept [55], [56] andfor 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 ofassessing dynamic behaviour.3. Participation factors computation from discrete time eigenvalues.66Chapter 4. Application to the OVNI simulator4. Identification of segmentation schemes from discrete time eigenvalueinformation applied to Latency.5. Automatic selection of variable discretization steps from discrete timeeigenvalue information.4.2 OVNIThe Object Virtual Network Integrator (OVNI) is UBC’s power system real-time simulator. OVNI’s discrete time solution is based on the EMTP dis-cretization method. It was developed by Mart´ı and Linares to simulate largeelectric networks in real-time and its software and hardware architectures aredescribed in the following publications [10], [18], [13], [19], [21], [22] and [57],among others. The OVNI simulator is built around the Multi Area The´veninEquivalent (MATE) concept, which is described in the following section, andwas developed as an object oriented application written in C++. This char-acteristic facilitates the integration of new models and modules with OVNI.The OVNI solution achieved real-time performance in a single PC environ-ment allowing interaction with analog devices for protective relay testing ap-plications, as introduced by Calvin˜o in [16]. The OVNI PC-cluster solutiondeveloped by Hollman [17] showed the OVNI algorithm’s inherent advan-tages for parallelization allowing for real-time simulation of large networks.The introduction of the Multilevel MATE concept by Armstrong [57] pro-vides OVNI with the capability for efficient simultaneous solution of control67Chapter 4. Application to the OVNI simulator+ -iαz VαBAFigure 4.1: Generic description of a MATE linksystems and nonlinearities.4.2.1 The MATE ConceptThe MATE concept developed by Mart´ı [58], [13], [19] was introduced to copewith the dimension of large power system networks for real-time and on-linedynamic studies. It also provides a framework for the implementation of theLatency concept which will be summarized later in this chapter. MATE isused in OVNI to segment a complex network into independent subnetworksinterconnected by links. A link can be a transmission line, a voltage sourceor just simply an ideal link to allow the desired segmentation. The mostgeneral description of a link is shown in Fig. 4.1. To demonstrate the MATEconcept, 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 2GB23 is the element (2, 3) of the admittance matrix [B]68Chapter 4. Application to the OVNI simulatorhA1 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 α1zα1 is the impedance of branch α1hA2A2A1A3A4B3B2B1hB2hA4hA1gA23 gA34gA14gA21hB1gB13gB12gB32[ A ] [ B ]Link α1Link α2iα1iα2Zα2Zα1Figure 4.2: Generic circuit for discrete state space formulation with MATEFor the considered network, the general nodal equation of the system withthe links α1 and α2 removed is, A 00 B vAvB = hAhB (4.1)69Chapter 4. Application to the OVNI simulatorwhile the full system with the links α1 and α2 included is given byGA11 GA12 0 GA14 0 0GA21 GA22 GA23 0 0GA31 GA32 GA33 GA34 1 0GA41 0 GA34 GA44 0 1GB11 GB12 GB13 −1 0GB21 GB22 GB23 0 0GB31 GB32 GB33 0 −10 0 1 0 −1 0 0 −z1 00 0 0 1 0 0 −1 0 −z2vA1vA2vA3vA4vB1vB2vB3iα1iα2=hA1hA2hA3hA4hB1hB2hB300(4.2)70Chapter 4. Application to the OVNI simulatorwith its compact notation,A 0 p0 B qpt qt −zvAvBiα =hAhB0 (4.3)Where [p] and [q] are the injection matrices derived from the link branches(if the current comes out +1, if the current goes in −1).p =0 00 0+1 00 +1(4.4)q =−1 00 00 −1 (4.5)After manipulations, as explained in [21], equation (4.2) can be representedby (4.6), the MATE equivalent circuit which is implemented in the central71Chapter 4. Application to the OVNI simulatorsolver,1 a13 a141 a23 a241 a33 a341 a43 a441 −b11 −b131 −b21 −b231 −b31 −b33Zα11 Zα12Zα21 Zα22vA1vA2vA3vA4vB1vB2vB3iα1iα2=eA1eA2eA3eA4eB1eB2eB3eα1eα2(4.6)where,a = A−1p (4.7)b = B−1q (4.8)eA = A−1hA (4.9)eB = B−1hB (4.10)e = pteA + qteB (4.11)Z = pta+ qtb+ z (4.12)72Chapter 4. Application to the OVNI simulator4.2.2 Discrete State Space Formulation with MATEMATE allows for the implementation of the proposed new methodology toextract eigenvalues from the discrete time EMTP solution. The formulationof the state space description can be implemented at a subnetwork level,while the global analysis of the eigenvalues is performed at the core leveleither in a central or dedicated node of the hardware cluster. Each sub-network node contains all the information required to formulate the discretestate space equations for that subnetwork. Each subnetwork has to submit tothe central solver its discrete state space description, which is then relayedto the dedicated eigenvalue solver. Since the central node solver is oftencurrently used to preprocess the data input file, the process of constructingthe discrete state space formulation of the network can be achieved at once,together with the case preprocessing stage. Consequently, after the simula-tion is started, only those subnetworks in which there is a topology changeor there is a change of the nonlinear region will have to submit to the cen-tral solver the updated discrete state space formulation. This is particularlyefficient for the cases where the nonlinearities are isolated in subnetworks,since the computation size of the problem is reduced (see Fig. 4.3). Here thesubnetworks A2 and A3 are the only ones in need to submit its state spacedescription update to the dedicated eigenvalue solver. With all the systeminformation allocated in the dedicated eigenvalue solver node, the discretetransition matrix [A]d of the complete network can be computed. The con-tents of the segmented blocks’ transition matrices are located in the correct73Chapter 4. Application to the OVNI simulatorLinearsubnetworkA1LinearsubnetworkA2with topologychangenonlinearsubnetworkA3EigenvaluesolverUpdate is requiredonly when changeof regionUpdate is requiredonly when changeof topologyNo update requiredFigure 4.3: General update scheme for discrete state space formulation withMATEplace of the global transition matrix using the information of connectivitypresent 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 constructedfrom the data input case file. The subnetwork’s discrete transition matrices[AA]d and [AB]d are first constructed from the data input case, and later dur-ing the simulation are updated with the information from the subnetwork incase of change of topology or change of region of a nonlinear subnetwork.The new hardware network interconnect layout of OVNI [22] can smoothly74Chapter 4. Application to the OVNI simulatorrun simultaneously both hard and soft real-time environments. Thus, it ispossible to allocate all the required eigenvalue analysis of the electrical net-work within the soft real-time environment and interacting with the hardreal-time cluster when required, as depicted in Fig. 4.4. These interactionsexist whenever the topology of the network changes or whenever the dynam-ics of elements is updated in the stability analysis [59], [60].As an alternative to the proposed centralized computation of eigenvalueSubsystemIOVNIcentralsolverSubsystemIIISubsystemIIDedicatedEigenvalue solverReal-timeenvironmentSoft real-timeenvironmentInteraction requiredonly when there isa change in [A] iInteractionsrequired ateverytime steplinkslinkslinksFigure 4.4: Hybrid real-time/soft real-time simulator layouttrajectories, a segmented distributed computation can be implemented ex-75Chapter 4. Application to the OVNI simulatortending the MATE partitioning concept [13], [19] to the eigenvalue analysis.Computing the eigenvalues at the tearing blocks level can provide a bettercomputational performance in some cases for large scale systems by concen-trating the computation of eigenvalues within a much smaller structure. Inan extension of the Latency concept discussed in the next section, the sub-system eigenvalues can be determined by its MATE Thevenin equivalent, asseen from the links subsystem. This approach, however, will only producegood results when the Latency approach is effective, which is when there isa large time constant separation among subsystems. The Latency approachis discussed next.4.3 The Latency ConceptThe Latency concept3, introduced by Mart´ı and Linares in [55] and by Mor-eira in [56], aims at exploiting the use of different time steps in differentsubnetworks. By using larger time steps for slow subnetworks and smalltime steps for fast subnetworks considerable computational savings can beachieved without sacrifying accuracy. The Latency concept [20], [56] dependsfor its successful implementation on the accurate identification of areas withsignificantly different time constants. This is not by any means an easy task,and as quoted in [56] “it is one of the limitations of the Latency method...”.3In this work, like in previous power system references, the term “latency” is usedto mean “latent” or “dormant”. In Signal processing theory this concept corresponds tomultirate systems76Chapter 4. Application to the OVNI simulator4.3.1 Application of the New Discrete Time StateSpace Formulation to LatencyThe methodology presented in this work allows for the identification of dom-inant eigenvalues in the discrete time domain using the EMTP description.Each subnetwork can be described in the state space domain by equation(3.22) in which the eigenvalues of [A] determine directly the time-responseof that subnetwork. After the system eigenvalues and participation matrixare computed an efficient segmentation and correct time step selection canbe achieved.In order to produce a correct segmentation, the concept of the participationmatrix [3] is used. The Participation matrix [P ], is given byP = [p˜1p˜2 . . . p˜k] (4.13)withp˜i =φ1iϕi1φ2iϕi2...φniϕik(4.14)whereφki = kth entry of the right eigenvector φiϕik = kth entry of the left eigenvector ϕi.The dimensionless element Pki provides a measure of the relative partic-77Chapter 4. Application to the OVNI simulatoripation of the kth state in the ith mode. The participation matrix combinesright and left eigenvectors showing the association between state space vari-ables and modes. The relationship between states and modes can not easilybe extracted directly from the right and left eigenvectors due to their depen-dency on scaling. The Participation matrix overcomes this situation.The algorithm can be summarized as follows:i. obtain the discrete time state space formulation of the network from theEMTP solution, that is, compute [A]d with equation (3.50)ii. compute the discrete time eigenvalues of the system from [A]diii. reconstruct the continuous time eigenvalues with equation (3.64) in thecase of the trapezoidal integration ruleiv. compute the Participation matrix [P ] from equation (4.13)v. identify in [P ], clusters of elements sharing similar time constantsvi. apply the latency method to partition the system with the appropriatetime stepsThis procedure should be repeated whenever there is a change in the topologyof the system.In section 5.2 a detailed example is provided to illustrate this methodology.78Chapter 4. Application to the OVNI simulator4.4 Automatic Time Step SelectionWith the eigenvalue information obtained from the discrete state space for-mulation from the EMTP solution, it is possible to implement an automatictime step selection option in the EMTP. Assume the user specifies a ∆twhich is adequate for the maximum frequency of interest in the study. Theexact continuous time eigenvalues can be obtained from the discrete timeEMTP formulation as discussed in this Thesis. From the exact continuoustime eigenvalues the maximum frequency present in the system (below themaximum study frequency) can be determined. If this frequency is higherthan say a fifth of the Nyquist frequency then a smaller ∆t that satisfies thiscondition must be selected. The automatic selection of the simulation timestep must be performed every time a nonlinear element changes its operationregion or every time a topology change is applied. The general scheme isillustrated in Fig. 4.5.79Chapter 4. Application to the OVNI simulatorMaximum frequency of interest inthe study, specified by userTime step specified by user in datainput fileDiscrete time eigenvaluescomputation from EMTP solutionMaximum frequency is identifiedIf fmax system > fmax userIssue user notificationIs fmax system ofinterest?noyesRun simulation withfixed user time stepSelect best appropriate time stepRun simulationChange of topologyor nonlinear operationregion?noyesDiscrete time eigenvaluescomputation from EMTP solutionFigure 4.5: Automatic EMTP time step selection scheme80Chapter 5Validation5.1 Comparison of State Space Solutionsbetween Continuous Time and DiscreteTime DomainsFor comparison purposes and a detailed illustration of the proposed tech-nique, a state space formulation of a simple RLC series circuit, see Fig. 5.1,will be obtained using both the traditional approach from the continuoustime domain and the proposed new technique from the nodal EMTP matrix.This second order circuit is an interesting initial case, because as pointed outin [3], the performance of higher order systems is often viewed in terms of adominant second order system.81Chapter 5. ValidationRCLVin(t)i(t) vC(t)vL(t)Figure 5.1: RLC series circuit5.1.1 RLC Series Circuit State Space Equation fromContinuous Time DomainThe continuous time state space system equations are given byx˙(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 circuitvin = eL + eR + eC (5.3)whereeL = Ldidt; eR = R i ; eC =1C∫i(t) dt82Chapter 5. Validationi(t) is the current in the loop. Substitutingvin(t) = Ldi(t)dt+R i(t) +1C∫i(t) dt (5.4)Introducing the capacitor charge as a variable, i(t) = dq(t)/dt,vin(t) = Ld2q(t)dt2+Rdq(t)dt+1Cq(t) (5.5)If we select q(t) and i(t) as the state variables and define the input andoutput asx(t) =q(t)i(t)u(t) = vin(t) (5.6)y(t) = vL(t)then the state space equations of the RLC series circuit aredq(t)dtdi(t)dt︸ ︷︷ ︸rate of change= 0 1−1LC−RL︸ ︷︷ ︸[A]q(t)i(t)︸ ︷︷ ︸present states+01L︸︷︷︸[B]vin(t)︸ ︷︷ ︸input(5.7)vL(t)︸ ︷︷ ︸output=[−1C−R]︸ ︷︷ ︸[C]q(t)i(t)︸ ︷︷ ︸states+ 1︸︷︷︸[D]vin(t)︸ ︷︷ ︸input(5.8)83Chapter 5. ValidationIf 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 transitionmatrix [A]c and the eigenvalues [λ]c are[A]c= 0 1−800 −20 (5.9)[λ]c=−10.0000 + j26.4575−10.0000− j26.4575 (5.10)Another possible representation of the system can be obtained by se-lecting vR(t) and vC(t) as the state variables. In this case the state spaceequations for the RLC series circuit arex(t) =vR(t)vC(t)u(t) = vin(t) (5.11)y(t) = vL(t)dvR(t)dtdvC(t)dt︸ ︷︷ ︸rate of change=−RL −RL−1RC0︸ ︷︷ ︸[A]vR(t)vC(t)︸ ︷︷ ︸states+RL0︸︷︷︸[B]vin(t)︸ ︷︷ ︸input(5.12)84Chapter 5. ValidationvL(t)︸ ︷︷ ︸output=[−1 −1]︸ ︷︷ ︸[C]vR(t)vC(t)︸ ︷︷ ︸states+ 1︸︷︷︸[D]vin(t)︸ ︷︷ ︸input(5.13)The obtained continuous time transition matrix [A]c and the eigenvalues [λ]care[A]c=−20 −2040 0 (5.14)[λ]c=−10.0000 + j26.4575−10.0000− j26.4575 (5.15)We can see that the continuous time eigenvalues [λ]c are the same in bothcases (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 fromthe [G] EMTP MatrixThe EMTP equivalent discretized with the trapezoidal rule RLC series circuitis shown in Fig. 5.2. The nodal equations in matrix form are85Chapter 5. ValidationhC(t)hL(t)32gC equivgL equivVin(t)1RFigure 5.2: Discretized RLC series circuit1R−1R0−1R( 1R+ −∆t2L) −∆t2L0 ∆t2L(∆t2L+ 2C∆t)v1v2v3 =i1hL−hL + hC (5.16)Since v1 = vin and is a known voltage source, we reduce the order of thesystem by 1 by moving the node to the right hand side of the equation.Re-arranging the [G] matrix86Chapter 5. Validation( 1R+ ∆t2L) −∆t2L−1R−∆t2L(∆t2L+ 2C∆t) 0−1R0 1Rv2v3v1 =hL−hL + hCi1 (5.17) GAA GABGBA GBB vAvB = hAhB (5.18)( 1R + ∆t2L) −∆t2L−∆t2L(∆t2L+ 2C∆t)︸ ︷︷ ︸[G]v2v3 = hL−hL + hC︸ ︷︷ ︸[h]−−v1R0︸ ︷︷ ︸input(5.19)v2v3 =( 1R + ∆t2L) −∆t2L−∆t2L(∆t2L+ 2C∆t)−1  hL + v1R−hL + hC (5.20)The relation between branch and node voltages and history sources isgiven by the incidence matrix [L]:_v′1_v′2_v′3 =1 0−1 10 1v′2v′3 (5.21)h′2h′3 =1 −1 00 1 1h′R−h′LhC (5.22)87Chapter 5. ValidationSubstituting the nodal voltages (5.20) into (5.21),_v′1_v′2_v′3 =1 0−1 10 1( 1R + ∆t2L) −∆t2L−∆t2L(∆t2L+ 2C∆t)−1  h′L + v1R−h′L + h′C (5.23)For the considered RLC series circuit, [k∧] and [g∧] (3.45) are[k∧]=0 0 00 1 00 0 −1 (5.24)[g∧]=0 0 00 −∆tL00 0 4C∆t (5.25)while the branch histories, according to equation (3.49), can be obtained as_h1_h2_h3 =0 0 00 1 00 0 −1_h′1_h′2_h′3+0 0 00 −∆tL00 0 4C∆t1 00 11 −1( 1R + ∆t2L) −∆t2L−∆t2L(∆t2L+ 2C∆t)−1 1 0 10 1 −1_h′1_h′2_h′3(5.26)Finally, the discrete time transition matrix [A]d for the circuit is computed88Chapter 5. Validationfrom the nodal [G] matrix as indicated in (3.50),[A]d=0 0 00 1 00 0 −1+0 0 00 −∆tL00 0 4C∆t1 00 11 −1( 1R + ∆t2L) −∆t2L−∆t2L(∆t2L+ 2C∆t)−1 1 0 10 1 −1(5.27)Assuming the same circuit values and discretization time of 50µS:R L C V ∆t100Ω 5H 250µF 8.5V 50µSwe obtain the following discrete time transition matrix[A]d=0 0 01.99799−3 9.97998−1 −3.99599−61.99799−3 1.99799 9.99996−1 (5.28)The discrete time eigenvalues of this matrix are[z]d=9.98997−1 + j2.64310−39.98997−1 − j2.64310−30 (5.29)The eigenvalues of the transition matrix [A] are the poles of the corre-sponding linear time-invariant system. In order to validate the proposed89Chapter 5. Validationtechnique, it is desirable to evaluate the precision of the eigenvalues com-puted from the nodal [G] matrix against those obtained from the continuoustime domain. For the case of the trapezoidal integration rule, the continuoustime eigenvalues can be computed from the discrete time transition matrixobtained from the nodal [G] matrix with the proposed method using equation(3.64),λi =2∆t(zi − 1)(zi + 1)(5.30)where λi are the continuous time eigenvalues and zi are the discrete timeeigenvalues . For the RLC circuit, the reconstructed continuous time eigen-values from the discrete time eigenvalues, obtained with equation (3.64), witha ∆t of 100µs are[λ]creconstructed=2100µs−0.00050 + j0.00132−0.00050− j0.00132=−10.0000 + j26.4575−10.0000− j26.4575 (5.31)The continuous time eigenvalues of equation (5.10) and the reconstructedcontinuous time eigenvalues of equation (5.31) are exactly the same, thusvalidating the proposed technique. The example validates the discussion inchapter 3 that the exact eigenvalues of a continuous time system can be90Chapter 5. Validationexactly computed from the discrete time nodal method independently fromthe chosen integration time step. It is important to note that even when avery large discretization time step is selected, the continuous time eigenvaluescan still be recomputed “exactly”, with possibly some small error due totruncation. For instance, if we take ∆t = 1ms[λ]creconstructed=21ms−0.0050 + j0.01322−0.0050− j0.01322=−10.0000 + j26.4400−10.0000− j26.4400 (5.32)Even though the eigenvalues can still be reconstructed for large ∆t′s, inthe important case of a circuit with nonlinearities, large ∆t′s may move theoperating point too far in the nonlinear characteristic as we move from thesolution step to the next. In this case, we may miss intermediate criticalpoints in the operating characteristics.The real part of the eigenvalues gives the damping, while the imaginary partgives the frequency of oscillation. The circuit’s oscillation frequency f anddamping ratio ζ can be computed from the reconstructed continuous timeeigenvalues (5.31) as follows, letλ = σ ± jω (5.33)91Chapter 5. Validationthenf =ω2pi(5.34)andζ =−σ√σ2 + ω2(5.35)A positive damping ratio ζ provides damping to the natural frequency.If ζ is greater than 1, both eigenvalues are real and negative; if ζ is equal to1, both eigenvalues are equal to −wn; and if ζ is less than 1, the eigenvaluesare complex and conjugates.For the series RLC test case the reconstructed continuous oscillation fre-quency and damping ratio aref c = 4.21Hz (5.36)ζc = 0.35355 (5.37)The LC resonant frequency of the circuit f cn (“natural frequency”) isf cn =12pi√LC= 4.50Hz (5.38)and the damping ratio ζc isζc =R2√LC= 0.35355 (5.39)Notice that in a stable system due to the losses the oscillation frequency92Chapter 5. ValidationVgenL1C1 C2R1 L2State 1State 4State 3State 2Figure 5.3: Latency test circuit(5.34) is lower than the natural frequency (5.38).5.2 Latency CaseIn this section the simple test case shown in Fig. 5.3 from [56], is used todemonstrate the capabilities of the presented method for an efficient imple-mentation of Latency. The values of the circuit components areR1 L1 C1 L2 C20.1Ω 1µH 100µF 1µH 1µFFor validation purposes the eigenvalues will be also computed from the con-tinuous time description.Selecting the capacitor voltages and inductor currents as state variables, we93Chapter 5. Validationobtain the following state space equations describing the test circuit,iC1 = C1dvC1dt=⇒ x˙1 = 1C1 iC1vL1 = L1diL1dt=⇒ x˙2 = 1L1vL1iC2 = C2dvC2dt=⇒ x˙3 = 1C2 iC2vL2 = L2diL2dt=⇒ x˙4 = 1L2 iL2(5.40)Applying Kirchoff’s laws,x˙1 =1C1(vL1 − vL2)x˙2 =1L1(Vgen − vC1)x˙3 =1C2vL2x˙4 =1L2(vC1 − vC2 +R1vL2)(5.41)The continuous time [A] matrix for the test circuit is given byAc =0 1C10 −1C1−1L10 0 00 0 0 1C21L20 −1L2−R1L2(5.42)94Chapter 5. Validationfrom which the following eigenvalues are obtainedλc =−4.99950x104 + j1.00379x106 λ1−4.99950x104 − j1.00379x106 λ2−4.99801 + j9.94988x104 λ3−4.99801− j9.94988x104 λ4 (5.43)We next compute the eigenvalues and eigenvectors using the discrete timemethodology. The network’s [G] matrix from the discrete time description isG =∆t2L1−∆t2L10 0−∆t2L1( 1R1+ ∆t2L1+ 2C1∆t) −1R100 −1R1( 1R1+ ∆t2L2) −∆t2L20 0 −∆t2L2( ∆t2L2+ 2C2∆t)(5.44)the matrix [k∧] is given byk∧ =1 0 0 00 −1 0 00 0 1 00 0 0 −1(5.45)95Chapter 5. Validationwhile the matrix [g∧] is given byg∧ =−∆tL10 0 00 2C1∆t0 00 0 −∆tL200 0 0 4C2∆t(5.46)The incidence matrix [L]L =1 −1 0 00 1 0 00 0 1 −10 0 0 1(5.47)Applying (3.50), we obtain the discrete time transition matrix [A]d[A]d=[k∧]+[g∧] [L] [G−1] [Lt]Ad =9.99999x10−1 1.25000x10−15 1.24999x10−15 0−2.00000 9.99999x10−1 1.99999 1.25034x10−131.24999x10−15 −1.24999x10−15 9.99999x10−1 1.24999x10−13−1.24999x10−15 1.24999x10−15 −1.99999 9.99999x10−1(5.48)96Chapter 5. Validationfrom where the following discrete time eigenvalues are computed,λd =9.99999x10−1 + j4.97494x10−89.99999x10−1 + j4.97494x10−89.99999x10−1 + j5.01896x10−79.99999x10−1 − j5.01896x10−7 (5.49)Finally, the reconstructed eigenvalues from the discrete time domain givenby (3.64) are,λcreconstructed =−4.99752 + j9.94988x104−4.99752− j9.94988x104−4.99950x104 + j1.00379x106−4.99950x104 − j1.00379x106 (5.50)which matches up to the second decimal for the worse approximation thecontinuous time eigenvalues calculated directly from the continuous time so-lution (5.43). The participation matrix [P ] (4.13) can now be computed inorder to determine how Latency will be implemented. For the test circuit,97Chapter 5. Validationthe participation matrix [P ] computed from the discrete time domain isP crec. =λ3creconstructed λ4creconstructed4.9995x10−1 + j2.6108x10−5 4.9995x10−1 − j2.6108x10−5 x˙24.9495x10−1 + j7.5560x10−5 4.9495x10−1 − j7.5560x10−5 x˙14.9970x10−5 − j9.9168x10−7 4.9970x10−5 + j9.9168x10−7 x˙45.0474x10−3 − j1.0067x10−4 5.0474x10−3 + j1.0067x10−4 x˙3λ2creconstructed λ1creconstructed4.9470x10−5 − j7.5406x10−6 4.9470x10−5 + j7.5406x10−6 x˙25.0479x10−3 − j2.6137x10−4 5.0479x10−3 + j2.6137x10−4 x˙14.9995x10−1 + j2.4910x10−2 4.9995x10−1 − j2.4910x10−2 x˙44.9495x10−1 − j2.4641x10−2 4.9495x10−1 + j2.4641x10−2 x˙3(5.51)As can be seen, the eigenvalues λ1creconstructed and λ2creconstructed are stronglyassociated to the states of x˙3 and x˙4 (corresponding to C2 and L2), whileeigenvalues λ3creconstructed and λ4creconstructed are strongly associated to thestates of x˙1 and x˙2 (corresponding to C1 and L1).The oscillation frequencies given by (5.34) are 15.8kHz for the slow sub-network and 159kHz for the fast one. Consequently, the time constants ofλ1 and λ2 are very loosely coupled with λ3 and λ4 allowing for a suitable98Chapter 5. ValidationVgenState 1State 4State 3State 2L1C1 C2R1 L2Slow FastFigure 5.4: Circuit with slow and fast subareassegmentation into a fast subcircuit and a slow subcircuit to implement theLatency concept, as shown in Fig. 5.4.VgenState 1 State 4State 3State 2R1 L1C1 C2R2 L2SW1Rload1 2 3 4Figure 5.5: Two-area resonant test case5.3 Two-Area Resonant CircuitA two-area resonant circuit, Fig. 5.5, is presented here. The eigenvalueanalysis is performed with the discrete time state space algorithm, for the99Chapter 5. Validationcases of a poorly damped and an improved damped configuration.The selected values for the simulation are, for the poorly damped case I,R1 R2 C1 C2 L1 L2 Vgen ∆t3Ω 5Ω 10nF 50nF 350mH 10mH 230kV 50µsand for the improved damped case II,R1 R2 C1 C2 L1 L2 Vgen ∆t3Ω 5Ω 10nF 1µF 350mH 1mH 230kV 50µsThe source is applied at t = 0 with switch SW1 closed. At t = 25ms SW1 isopened and the two resonant circuits become isolated.5.3.1 Case I - Poorly DampedThe voltage outputs for the poorly damped caseI are shown in Fig. 5.6.Using the proposed methodology the trajectory of the eigenvalues of the testnetwork is computed. The discrete eigenvalues are presented in 5.7 and thereconstructed continuous time eigenvalues are presented in Fig. 5.8. Thereare two sets of eigenvalues which correspond to the circuit topology beforeand after the switching operation. The [G] discretized network matrix withswitch SW1 open is given by[Gopen] =1R1+2L1∆t−1R1+2L1∆t0 0− 1R1+2L1∆t1R1+2L1∆t+ 2C1∆t +1Ropen−1Ropen00 −1Ropen1Ropen+ 2C2∆t +1R2+2L2∆t−1R2+2L2∆t0 0 −1R2+2L2∆t1R2+2L2∆t+ 1Rload(5.52)100Chapter 5. Validationand with the switch SW1 closed,[Gclose] =1R1+2L1∆t−1R1+2L1∆t0 0− 1R1+2L1∆t1R1+2L1∆t+ 2C1∆t +1Rclose−1Rclose00 −1Rclose1Rclose+ 2C2∆t +1R2+2L2∆t−1R2+2L2∆t0 0 −1R2+2L2∆t1R2+2L2∆t+ 1Rload(5.53)the [k∧] matrix of branch nature coefficients is given byk∧ =1 0 0 0 00 −1 0 0 00 0 −1 0 00 0 0 1 00 0 0 0 0(5.54)while [g∧] matrix of branch discretization rule coefficients isg∧ =−2R1+2L1∆t0 0 0 00 4C1∆t0 0 00 0 4C2∆t0 00 0 0 −2R2+2L2∆t00 0 0 0 0(5.55)101Chapter 5. Validationand the incidence matrix [L] isL =1 −1 0 00 1 0 00 0 1 00 0 1 −10 0 0 1(5.56)The transition matrices [A] for the open and closed switch conditions areobtained with 3.50. Finally the discrete time eigenvalues areλdclosed =−1.419878x10−2 + j9.998990x10−1−1.419878x10−2 − j9.998990x10−1−9.999999x10−1−1.000000(5.57)λdopen =−1.049723x10−1 + j9.944750x10−1−1.049723x10−1 − j9.944750x10−19.999999x10−1−1.000000(5.58)while the reconstructed continuous time eigenvalues computed with (3.64)102Chapter 5. Validation00.0050.010.0150.020.0250.030.0350.040.0450.05−505x 105 Voltage Node 1time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−101x 106 Voltage Node 2time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−202x 104 Voltage Node 3time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−0.500.5 Voltage Node 4time [s][V]Figure 5.6: Node voltages for the two-area example of Fig. 5.5. Case I,poorly damped103Chapter 5. Validation−1−0.500.51−1−0.500.5100. [s] eig1 − L2eig2 − C2eig3 − C1eig4 − L1Figure 5.7: Discrete time eigenvalues for the two-area example of Fig. 5.5.Case I, poorly dampedareλcclosed =−5.072087x10−3 + j4.057204x104−5.072087x10−3 − j4.057204x104−1.200000x1012−4.000000x104(5.59)λcopen =−4.938281x10−3 + j4.444444x104−4.938281x10−3 − j4.444444x104−1.000044x10−7−4.000000x104(5.60)The continuous time damping ratios computed with (5.35) for both topolo-104Chapter 5. Validation−14−12−10−8−6−4−20x 1011−505x 10400. [s]eig1 − L2eig2 − C2eig3 − C1eig4 − L1Figure 5.8: Reconstructed continuous time eigenvalues for the two-area ex-ample of Fig. 5.5. Case I, poorly dampedgies areζclosed = 1.250x10−7ζopen = 1.111x10−7 (5.61)and the damped oscillation frequencies (5.34) arefclosed = 6.457kHzfopen = 7.073kHz (5.62)As can be observed, the analyzed circuit is poorly damped for the oscillation105Chapter 5. Validationfrequencies for both open switch and closed switch conditions.5.3.2 Case II - Improved DampingThe same circuit is once again analyzed with C2 and L2 changed in order toimprove the damping of the oscillation frequency. The new set of eigenvaluesand damping ratios for case II with improved damping are shown below.Discrete time eigenvalues,λdclosed =2.902763x10−1 + j9.569413x10−12.902763x10−1 − j9.569413x10−1−9.99999x10−1−1.00000(5.63)λdopen =2.857138x10−1 + j9.583134x10−12.857138x10−1 − j9.583134x10−1−9.99999x10−1−1.00000(5.64)while the reconstructed continuous eigenvalues areλcclosed =−4.444838x10−2 + j2.966627x104−4.444838x10−2 − j2.966627x104−1.009999x1012−4.000000x104(5.65)106Chapter 5. Validation00.0050.010.0150.020.0250.030.0350.040.0450.05−505x 105 Voltage Node 1time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−101x 106 Voltage Node 2time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−202x 104 Voltage Node 3time [s][V]00.0050.010.0150.020.0250.030.0350.040.0450.05−0.500.5 Voltage Node 4time [s][V]Figure 5.9: Node voltages for the two-area example of Fig. 5.5. Case II,improved damping107Chapter 5. Validation−1−0.500.51−1−0.500.5100. [s]eig1 − L2eig2 − C2eig3 − C1eig4 − L1Figure 5.10: Discrete time eigenvalues for the two-area example of Fig. 5.5.Case II, improved dampingλcopen =−4.444444x10−2 + j2.981423x104−4.444444x10−2 − j2.981423x1041.000044x10−7−4.000000x104(5.66)The new damping ratios areζclosed = 1.498x10−6ζopen = 1.490x10−6 (5.67)108Chapter 5. Validation−12−10−8−6−4−20x 1011−3−2−10123x 10400. [s]eig1 − L2eig2 − C2eig3 − C1eig4 − L1Figure 5.11: Reconstructed continuous time eigenvalues for the two-area ex-ample of Fig. 5.5. Case II, improved dampingand the new damped frequencies (5.34) arefclosed = 4.721kHzfopen = 4.745kHz (5.68)The improvement in damping due to the changes applied to L2 and C2 isclearly 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 associatedto L2 and C2 lay inside the unit circle in the discrete time domain, providingthe circuit with a better damping ratio compared to the first case analyzed.109Chapter 5. ValidationVgenState 1State 3State 2R1 L1C1 C2R2L2SW1Non linealLoadState 4Figure 5.12: Two-area circuit test case with nonlinear load5.4 Two-Area Case with Nonlinear LoadA particularly interesting area of study is the identification of transient eigen-value trajectories in systems with poorly damped eigenvalues that are beingheavily stressed. These systems can present very critical operational trajec-tories. This situation can be typical of overstressed power systems operatingunder 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. Assumethe load is defined as followsRLoad = iC2(VC2)η (5.69)with η = 0.18. This case presents discrete time eigenvalues that move outsidethe unit circle (unstable region) during the circuit’s energization (closingswitch SW1 at t = 0). The discrete time eigenvalues trajectories for this110Chapter 5. Validation−1.5 −1−0.5 00.5 11.5−2−101200.51RealImagtime [s]eig1 − L2eig2 − C2eig3 − C1eig4 − L1−1.5−1−0.500.511.5−1.5−1−0.500.511.5RealImagFigure 5.13: Nonlinear load - discrete time eigenvaluescase are plotted in Fig. 5.13, while the continuous time ones are plottedin Fig. 5.15. In the upper plot of Fig. 5.13, it can be observed that theeigenvalues are moving around their ultimate steady state location. Also,due to the damping of the circuit elements, the excursion around the steadystate location decreases as the simulation evolves. If some of the eigenvalueslocated close to the instability region (unit circle) are excited during thetransient, their trajectories can lead to an eventual system collapse. In thiscase static eigenvalue analysis which would show only the final steady statepoint would not be sufficient to predict the eventual system collapse.Computing the point by point eigenvalues trajectories from the EMTP111Chapter 5. Validation−15 −10−5 05x 1011−505x 10400.51RealImagtime [s]eig1 − L2eig2 − C2eig3 − C1eig4 − L1−14 −12 −10 −8 −6 −4 −2 0 2x 1011−505x 104RealImagFigure 5.14: Nonlinear load - continuous time reconstructed eigenvalues fromdiscrete time solutionformulation with the presented methodology brings out the actual nature ofthe system and identifies the system components the behaviour is associatedwith, thus facilitating mitigation and corrective actions.5.5 Voltage Collapse in a Radial SystemFigure 5.16 shows a simple radial system. The load is increased until voltagecollapse is reached. As pointed out in [3], the Voltage collapse is the processby which the sequence of events accompanying voltage instability leads to112Chapter 5. Validation00.−505x 105 Voltage Node 1time [s]voltage [V]−101x 106 Voltage Node 2time [s]voltage [V]−101x 106 Voltage Node 3time [s]voltage [V]−101x 106 Voltage Node 4time [s]voltage [V]Figure 5.15: Two-area system with nonlinear load113Chapter 5. ValidationVgenZlineV1 V2P + jQPcteExp. loadFigure 5.16: A simple radial systema low unacceptable voltage profile in a significant part of the power system.Several indices to predict voltage collapse are available [61]. One of them isthe state based “voltage drops” index, which relates the no-load voltage withthe on-load voltage. The voltage drops index is given byVV0=ejφgrid−φload22 cos(φgrid−φload2)(5.70)Where V is the voltage in on-load situation and V0 is the voltage at thesame point in no-load situation with φgrid as the representative power angleof 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 acceptable114Chapter 5. Validationoperation values. Five characteristic transition points are identified for theanalysis, A1, A2, B, C and D in which the system voltage drops. For thisexample, an exponential model was used to represent the load,P = P0(VV0)α(5.71)Q = Q0(VV0)β(5.72)With α = 0.1, β = 0.6. These coefficients are typical of an industrial loadaggregate with a large induction motors component. A bank of capacitors Cfor power factor correction is connected in parallel with the load. Equations(5.71) and (5.72) are RMS values over one cycle of the signal. The amountof 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.2The continuous time eigenvalue trajectories of the radial system of Fig.5.16 are computed from the discrete time domain with the presented method-ology and shown in Fig. 5.18. As can be observed, there is a distinctivechange of the system eigenvalues location in the vicinity of the progressive115Chapter 5. Validation1 2 3 4 5 6 7 8 9 10 11 1200. [s]Voltage Bus 2 (pu)A1 A2 B C D Figure 5.17: Voltage collapse of a simple radial system116Chapter 5. Validation−1 −0.8−0.6 −0.4−0.2 0−2−1012024681012time [s] Real Imag A1 B C D A2 Figure 5.18: Reconstructed continuous time eigenvalue trajectory over timefor a simple radial system with voltage collapsevoltage drops at node v2. Points A1, A2, B, C and D in Fig. 5.17 can beclearly 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 thecomplex plane in Fig. 5.19 shows the value of the eigenvalues at the criticalpoints of Fig. 5.17. This simple representation of specific points is not alwaysenough to understand the system dynamics due to overlapping of eigenvaluesat different time instants. It does provide, nonetheless, a visual clue to ex-117Chapter 5. Validation−45 −40 −35 −30 −25 −20 −15 −10 −5 0 5−2−1.5−1−0.500.511.52RealImagA1CBDCDBA2 A2A1Figure 5.19: Complex plane projection of continuous time eigenvalues overtime line of simple radial system of Fig. 5.16 with voltage collapsetreme locations of eigenvalues. Figure 5.20 shows the complete trajectory ofthe discrete time eigenvalues computed from the nodal EMTP equations. Thezoomed in area showed in Fig. 5.21 highlights the fact that by analyzing theeigenvalues trajectories it is possible to anticipate the voltage collapse (i.e.,by observing the speed of change and damping ratio of eigenvalues). Boththe continuous and discrete time eigenvalues trajectories are the consequenceof the system’s search for a new equilibrium operation point. Since the load118Chapter 5. Validation0.950.960.970.980.991−2−1.5−1−0.500.511.52x 10−3024681012RealImagtime [s] A1 A2 B C DFigure 5.20: Discrete time eigenvalue trajectory from EMTP solution for thesimple radial system of Fig. 5.16 with voltage collapsedemand is forcing the system to its limit and the system is not capable ofsatisfying it, the voltage is dropped and the system is able to find at least atemporal 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, providesus with additional information about the type of transition of the system, avery fast time constant associated to the change with higher damping thanthe original eigenvalue locations. Immediately after the system suffers the119Chapter 5. Validation0.99920.99940.99960.9998100.511.5x 10−33456789101112RealImag[s]A1 A2 BC D Figure 5.21: Zoom in of discrete time trajectory of Fig. 5.20voltage drop (A1, A2, B, C and D), the system eigenvalues return to theinitial locations unless a new voltage drop is suffered, in which case a newtransition is observed.The method proposed in this work can provide sensitive information aboutproximity of system voltage collapse with sufficient anticipation time to ini-tiate corrective actions. For example, Fig. 5.22 shows the trajectory of themagnitude of eigenvalue 1, |eig1| over time. This trajectory provides a goodindication of the proximity of a voltage drop. It can be noticed for instance120Chapter 5. Validation0 2 4 6 8 10 120.99880.9990.99920.99940.99960.99981time [s]|eig 1|A1 A2 BC DFigure 5.22: Discrete time eigenvalue 1 magnitude, |eig1|, over timethat eigenvalue 1 starts to show evident signs of voltage collapse in the prox-imity of 250ms, before the actual drop to point A1, 300ms before point B,200ms before point C and 450ms before the point of interest, D. This an-ticipation is enough to initiate automatic corrective actions as, for example,line automatic reclosing, or load shedding.121Chapter 5. Validation5.6 11 Bus Two Area SystemFigure 5.23 shows a Two Area System of 11 buses described in [28]. Thesystem represents two areas with two generators in each area. The connec-tion between areas is achieved by several parallel transmission lines betweenbuses 7, 8 and 9. In steady state operation the area composed by generators1 and 2 is exporting 400 MW to the other area composed by generators 3 and4. The system is simulated in the discrete time domain using the RTTP sim-ulator [16], [17], [54] in which the proposed methodology was implemented.The simulation represents the tap changing transformers and it does not takeinto account the generator controllers. The dynamic behaviour of the systemuntil 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 showsthe PV curve at bus 9 with two characteristic points A and B, normal initialoperation and collapse point due to load stress, respectively. Using the dis-crete time state space formulation from the EMTP solution presented in thiswork, the discrete time eigenvalues of the circuit in Fig. 5.23 are computed.Figure 5.26 shows the discrete time eigenvalues of the system for the char-acteristic point A while Figure 5.27 illustrates the discrete time eigenvaluesfor the characteristic point B. Similarly to the case of the radial systempresented in the previous section, it can be observed that the eigenvaluelocations change in the vicinity of a voltage collapse.122Chapter 5. ValidationSwV:1.075Ph:9.64V:0.872Ph:-26.36V:1.030Ph:20.20V:1.010Ph:6.92V:0.954Ph:-12.57V:1.030Ph:-83.94SettingsV:0.986Ph:-73.61V:1.029Ph:-78.40V:1.000Ph:-4.39V:1.088Ph:-81.74V:1.010Ph:-83.94Bus_3G1 Voltage=1.03Phase=20.2 G2V=1.01P=10.5G3V=1.03P=-6.8T1X=0.0167Tap=1.1291T2X=0.0167Tap=1.1291T3X=0.0167Tap=1.1082T5X=0.0167Tap=1.0810 kmR=0.0020X=0.020B=0.0170110 kmR=0.0075X=0.075B=0.1283R=0.0075X=0.075B=0.1283Q=-3.5Steady State Order=7Steady State Frequency=60.0Max Steps=10000Stop in Collapse=yCollapse Voltage=0.3Ignore Tolerance=noCompute RMS=1Tolerance=1.0e-5MVAR Base=10025 kmR=0.0025X=0.025B=0.04375R=0.0010X=0.010B=0.0175 Q=-2.0PQ[7]P=9.67Q=1110 kmR=0.018X=0.18B=0.19R=0.018X=0.18B=0.192510 kmR=0.001X=0.01B=0.0175PQ[9]P=1.767Q=125 kmR=0.0025X=0.025B=0.0437G4V=1.01P=-17.0R=0.0075X=0.075B=0.1283R=0.0075X=0.075B=0.1283R=0.018X=0.18B=0.1925Bus_1Bus_5Bus_2Bus_7Bus_11Bus_6Bus_8Bus_9Bus_10Bus_4Figure 5.23: Two-area 11 bus test system123Chapter 5. Validation0510152025300.850.90.9511. (s)Bus Voltage Magnitude (pu)Bus1 Bus2 Bus3 Bus4 Bus5 Bus6 Bus7 Bus8 Bus9 Bus10Bus11Figure 5.24: Two-area 11 bus test system - Bus Voltages124Chapter 5. Validation020040060080010001200140016001800200000.−9)pf=0.870P=176.700, V=0.986P=1583.783, V=0.566PL(MW)V(PU)A B Figure 5.25: Two-area 11 bus test system - PV curve at Bus 9125Chapter 5. Validation0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005−0.25−0.2−0.15−0.1−0.0500. Real Figure 5.26: Two-area 11 bus test system - Discrete time Eigenvalues normaloperation point A126Chapter 5. Validation0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 1.005−0.25−0.2−0.15−0.1−0.0500. Real Figure 5.27: Two-area 11 bus test system - Discrete time Eigenvalues atcollapse point B127Chapter 6Conclusions andRecommendations6.1 ConclusionsThe main contribution of this Thesis is the construction of the discrete timetransition matrix as a function of the branch parameters and connectivity. Tothe best of the author’s knowledge this has never been previously proposed.The present work expands the capabilities of leading transient analysis pro-grams for power systems analysis based on the EMTP solution to producenot only the time domain waveforms of voltages and currents in the sys-tem, but also to trace the system’s dynamic and stability behaviour. Thisis achieved by extracting the discrete time eigenvalue information from thesystem difference equations at each time step of the solution. This infor-mation can be used directly in the discrete time domain by observing thestability margins around the unit circle, or by re-mapping the discrete timeeigenvalues into the continuous time eigenvalues of the original physical sys-tem. This backwards mapping is exact because the distortion introduced by128Chapter 6. Conclusions and Recommendationsthe discretization rule in going from continuous to discrete time is appliedin reverse when going from discrete to continuous time. Notice that this isa specific property of the eigenvalues and it does not apply to voltages andcurrents in the network.By using the EMTP as the base solution, system nonlinearities (as for ex-ample the nonlinear nature of the load), which have a fundamental effect onthe system dynamics, can be easily taken into account in the calculation ofthe eigenvalues because they are part of the EMTP network representationat each time step of the system’s response.The ability of transients programs, like the EMTP, Microtran, or OVNI,together with the extension presented in this thesis to incorporate eigen-value/eigenvector solutions at each point of the system trajectory will allowfor a much more accurate assessment and corrective control measures undercritical system operating conditions like those leading to extensive blackouts.System collapse during these situations is driven by a complex interaction offast and slow system dynamics. Conventional quasi-steady state solutions,like those of traditional transient stability programs, do not have the capabil-ity of considering these interactions using computational tools that are fastenough for on-line system operation.In the context of the OVNI simulator, an important application of the tech-nique proposed in this thesis is the identification of network areas that canbe solved with large time steps and areas that require small time steps. Asshown, the necessary information can be obtained from the participation ma-129Chapter 6. Conclusions and Recommendationstrix [P ] computed directly from the nodal description. The exploitation ofthis system latency is a key aspect in the ability of OVNI to represent verylarge systems in real-time using network partitioning techniques.This dissertation has presented the theory of the proposed eigenvalue/eigenvectoranalysis technique, as well as illustrative examples of its ability to detectsystem conditions leading to instability due to the dynamic behaviour of thesystem eigenvalues under situations of nonlinear loads. This type of analy-sis was not available before within the EMTP. Traditionally, power systemstability information could only be obtained for the system’s steady stateoperating points, thus ignoring the fact that under situations of high stress,like those leading to voltage collapse and extensive blackouts, the system’strajectory may no lead to stable steady state points. With the new presentedapproach, it is possible, with the appropriate control strategies, to redirectthe system’s trajectory before the system becomes unstable and collapses.Current work to be reported in future publications includes the combinationof this technique with the multilevel MATE concept of [57] which achievesa very efficient representation of control systems corrective actions in theOVNI simulator for real-time on-line applications.The ability to trace dynamically the trajectory of the eigenvalues of nonlinearsystems is of particular importance for the following studies:i. design and operation of embedded adaptive controllersii. automatic selection of discretization steps130Chapter 6. Conclusions and Recommendationsiii. determination of suitable strategies to avert blackouts in power systemsworking close to their operational limitsiv. dynamic small-signal stability analysis with the EMTPAn additional advantage of the discrete time eigenvalue trajectory analy-sis when compared to traditional power systems stability tools is its abilityto keep the physical network mapped to the eigenvalue/eigenvector analysislayer, thus allowing for a fast and accurate identification of elements con-tributing to instability.6.2 Recommendations for Future WorkThe 21st century power system will make extensive use of distributed intel-ligent local processing nodes capable of making their own decisions and ableto interact with other related nodes and with centralized regional controls.Since all parts of a power network interact with each other, local controlaction may require extensive information regarding the external system towhich the node is connected. This information can be acquired and transmit-ted to all the distributed nodes or, even more efficiently, it can be simulatedin real-time at each decision node. The use of MATE concepts together witheigenvalue/eigenvector trajectory analysis will allow for efficient embeddedcontrol devices.131Chapter 6. Conclusions and RecommendationsOne area of particular interest is that of distributed generation with renew-able energy sources such as wind farms. These generation clusters are highlynonlinear and require advance on-site control devices to allow a healthy in-terconnection with the power grid. The proposed methodology together withembedded simulation capabilities can provide a solution to this problem.Current research work [62],[63] within the Joint Infrastructure Interdepen-dencies Research Program (JIIRP) (sponsored by Public Safety and Emer-gency Preparedness Canada (PSEPC) and the Natural Sciences and Engi-neering Research Council of Canada (NSERC)) is concerned with the iden-tification of critical interdependencies among multiple infrastructures (e.g.,power grid, water grid, etc.) and the dynamic evolution of these interde-pendencies over the time line. The dynamic behaviour of interdependenciesamong critical infrastructures can be mapped into a multiple-infrastructurestransportation matrix which can then be transformed into a state spaceequation from which eigenvalues and trajectories can be obtained with themethodology presented in this work.The following study areas are suggested for future work:i. Development of distributed intelligent control solutions based on embed-ded OVNI and discrete state space eigenvalue trajectory computation.ii. Integration of discrete state space eigenvalue analysis methodology withthe Latency concept.iii. Development of more detailed Induction motor models.132Chapter 6. Conclusions and Recommendationsiv. Implementation of the discrete state space eigenvalue analysis method-ology within UBC’s OVNI-NET simulator for stability analysis and de-termination of segmentation schemes.v. Implementation of discrete state space eigenvalue analysis methodologyto UBC-JIIRP’s I2Sim simulator for identification of trajectories of crit-ical interdependencies among Infrastructures.133Bibliography[1] S. Gissinger, P. Chaunmes, J.P. Antoine, A. Bihain, and M. Stubbe,“Advanced dispatcher training simulator,” Computer Applications inPower, IEEE, vol. 13, no. 2, pp. 25–30, 2000.[2] M.V.F. Pereira, M. McCoy, and H.M. Merril, “Managing risk in thenew 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 powersystem engineering series. McGraw-Hill, 1994.[4] C.A.Da Neto, D.M. Falcao, and J.L. Jardim, “A unified online securityassesment system,” in CIGRE session 2000. 2000, pp. 38–102, CIGRE.[5] H.W. Dommel, EMTP Theory Book, Microtran Power System AnalysisCorporation, Vancouver, 2nd edition, 1992.[6] H.W. Dommel, “Digital computer solution of electromagnetic transientin single and multiphase networks,” IEEE Transactions on Power Ap-paratus and Systems, vol. PAS-88, no. 4, pp. 388–399, 1969.134Bibliography[7] G.D. Irwin, R. Kuffel, and D.AWoodford, “EMTDC- high performanceelectromagnetic 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 IEEEPES winter meeting. 2000, IEEE.[9] P.G. McLaren, R. Kuffel, R. Wierckx, J. Giesbrecht, and L. Arendt, “Areal time digital simulator for testing relays,” IEEE Transactions onPower Delivery, vol. 7, pp. 207–213, 1992.[10] J.R. Mart´ı and L.R. Linares, “Real-time EMTP-based transients sim-ulation,” Power Systems, IEEE Transactions on, vol. 9, no. 3, pp.1309–1317, 1994.[11] M. Kesunovic, M. Aganic, V. Skendzic, J. Domaszewicz, J.K. Bladow,D.M. Hamai, and S.M. McKenna, “Transients computation for relaytesting 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 InternationalConference 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 Power135BibliographySystem Technology, 1998. Proceedings. POWERCON ’98. 1998 Inter-national Conference on, 1998, vol. 2, pp. 977–981 vol.2.[14] Fromonto, Levacher, and Strunz, “An efficient parallel computer archi-tecture for real time simulations of electromagnetic transients in powersystems,” 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] Jesu´s Calvin˜o Fraga, Implementation of a real time simulator for pro-tective relay testing, M.A.Sc. Thesis. University of British Columbia,Vancouver, 1999.[17] J. A. Hollman, Real time distributed network simulation with pc clusters,M.A.Sc. Thesis. University of British Columbia, Vancouver, 1999.[18] J.R. Mart´ı, L.R. Linares, Roberto Rosales, and H.W. Dommel, “OVNI:a full-size real time power system simulator,” in International Confer-ence on Digital Power System Simulators, Montreal, Quebec, Canada,1997, pp. 1309–1317, ICDS’97.[19] L.R. Linares, OVNI (Object Virtual Network Integrator) a new fastalgorithm for the simulation of very large electric networks in real time,Ph.D. Thesis. University of British Columbia, Vancouver, 2000.136Bibliography[20] F. Moreira, J.A. Hollman, L.R. Linares, and J. R. Mart´ı, “Networkdecoupling by latency exploitation and distributed hardware architec-ture,” in International Conference on Power System Transients, Rio deJaneiro, 2001, IPST’01.[21] J.R. Mart´ı, L.R. Linares, and J.A. Hollman, “OVNI: Integrated soft-ware/hardware solution for real-time simulation of large power sys-tems,” in 14th Power Systems Computer Conference, Seville, Spain,2002, PSCC’02.[22] T. De Rybel, J.A. Hollman, and J. Mart´ı, “OVNI-NET: A flexiblecluster interconnect for the new OVNI real-time simulator,” in 15thPower Systems Computer Conference, Liege, Belgium, 2005, PSCC’05.[23] B. Gao, G.K. Morison, and P. Kundur, “Towards the developmentof a systematic approach for voltage stability assessment of large-scalepower 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.137Bibliography[28] Powertech Labs Inc., “TSAT: Transient security assessment tool,” 2002.[29] K.S. Shim, H.K. Nam, S.G. Song, Y.G. Kim, and K.Y. Lee, “Applica-tion results of the eigen-sensitivity theory of augmented matrix to smallsignal stability analysis of large power systems,” in Power EngineeringSociety 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-timesecurity assessment of electrical power systems,” Power Systems, IEEETransactions on, vol. 11, no. 2, pp. 1112–1117, 1996, TY - JOUR.[31] J.L. Jardim, “Online dynamic security assessment: implementationproblems and potential use of artificial intelligence,” in Power Engi-neering Society Summer Meeting, 2000. IEEE, 2000, vol. 1, pp. 340–345vol. 1, TY - CONF.[32] W.E. Arnoldi, “The principle of minimized iterations in the solution ofthe 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 fornoisy data: Choosing the signal components and selecting the orderin exponential signal models,” Proceedings of the IEEE, vol. 72, pp.230–233, 1984.138Bibliography[34] M.A. Johnson, I.P. Zarafonitis, and M. Calligaris, “Prony analysisand power system stability-some recent theoretical and applications re-search,” in Power Engineering Society Summer meeting, 2000, vol. 3,pp. 1928–1923.[35] Baron (Gaspard Riche) de Prony, “Essai experimental et analytique,”J. E. Polytech., vol. 1, no. 2, pp. 24–76, 1795.[36] R. Kumaresan and D. W. Tufts, “Estimation the parameters of ex-ponentially damped sinusoids and pole-zero modeling in noise,” IEEETransactions 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, “Comparisonof prony and eigenanalysis for power system control design,” PowerSystems, IEEE Transactions on, vol. 8, no. 3, pp. 964–971, 1993.[39] P. Kundur, G.K. Morison, L. Wang, and H. Hamadanizadeh, “On-linedynamic security assessment of power systems,” in Fifth InternationalWorkshop on Electrical Power Control Centres, Hungary, 1999.[40] Y. Saad, “Variations of arnoldi’s method for computing eigenelementsof large unsymmetric matrices,” Linear Algebra and its applications,vol. 34, pp. 269–295, 1980.139Bibliography[41] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press,Oxford, 1965.[42] D.T. Wong, B. Rogers, G. andPorreta, and P. P. Kundur, “Eigenvalueanalysis 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 inthe EMTP: steady-state initialization,” IEEE Transactions on PowerSystems, vol. 10, no. 2, pp. 593–601, 1995.[44] M. Kezunovic, L. Kojovic, A. Abur, C.W. Fromen, D.R. Sevcik, andF. Phillips, “Experimental evaluation of EMTP-based current trans-former models for protective relay transient study,” IEEE Transactionson Power Delivery, vol. 9, no. 1, pp. 405–413, 1994.[45] B. Khodabakhchian and G.T. Vuong, “Modelling a mized residential-commercial load for simulations involving large disturbances,” IEEETransactions on Power Systems, vol. 12, no. 2, pp. 791–796, 1997.[46] N. Nagaoka, S. Tawada, and A. Ametani, “Magnetizing circuit model forEMTP including voltage dependence characteristic,” in Transmissionand Distribution Conference and Exhibition 2002, Asia Pacific, 2002,vol. 2, pp. 1265–1269.140Bibliography[47] W. Long, D. Cotcher, D. Ruiu, P. Adam, S. Lee, and R. Adapa,“EMTP- a powerful tool for analyzing power system transients,” Com-puter Applications in Power, IEEE, vol. 3, no. 3, pp. 36–41, 1990.[48] J.A. Martinez, “How to adapt the EMTP for classroom instruction,”Power Systems, IEEE Transactions on, vol. 7, no. 1, pp. 351–358, 1992.[49] J.A. Martinez, “Educational use of EMTP models for the study ofrotating machine transients,” Power Systems, IEEE Transactions on,vol. 8, no. 4, pp. 1392–1399, 1993.[50] Kim Chul-Hwan, Lee Myung-Hee, R.K. Aggarwal, and A.T. Johns, “Ed-ucational use of EMTP models for the study of rotating machine tran-sients,” Power Systems, IEEE Transactions on, vol. 15, no. 1, pp. 9–15,2000.[51] J.R. Mart´ı and J. Lin, “Supression of numerical oscillation in theEMTP,” Power Systems, IEEE Transactions on, vol. 4, no. 2, pp.739–747, 1989.[52] J.A. Hollman and J.R. Mart´ı, “Real time network simulation with pc-cluster,” Power Systems, IEEE Transactions on, vol. 18, no. 2, pp.563–569, 2003.[53] J.R. Mart´ı, “EECE 560 class notes,” Department of Electrical andComputer Engineering, The University of British Columbia, Vancouver.141Bibliography[54] Jesu´s Calvin˜o Fraga, Tunable Difference Equations for the Time-Domain Simulation of Power System Operating States, Ph.D. Thesis.University of British Columbia, Vancouver, 2003.[55] L.R. Linares and J.R. Mart´ı, “Sub-area latency in a real time power net-work simulator,” in International Conference on Power System Tran-sients, Lisbon, 1995, pp. 541–545, IPST’95.[56] Fernando Moreira, Latency Techniques in Power System TransientsSimulation, Ph.D. thesis, University of British Columbia, 2002.[57] M. Armstrong, J.R. Mart´ı, L.R. Linares, and P. Kundur, “MultilevelMATE for efficient simultaneous solution of control systems and non-linearities in the OVNI simulator,” Power Systems, IEEE Transactionson, vol. 21, no. 3, 2006.[58] J.R. Mart´ı, “MATE, a multi-area thev´enin equivalent concept,” In-ternal memorandun of the power systems group, UBC - Department ofElectrical and Computer Engineering, August 1993.[59] A. Brameller, Y. M. Hamam, and Ronald Norman Allan, Sparsity : itspractical application to systems analysis, Pitman, London ; New York,1976.[60] A. Brameller, Maldwyn Noel John, and M. R. Scott, Practical diakopticsfor electrical networks, Advanced engineering textbooks. Chapman andHall, London,, 1969.142Bibliography[61] Working group 38.02, “Indices predicting voltage collapse includingdynamic phenomena,” Task force 11, December 1994.[62] J.R. Mart´ı, J.A. Hollman, C. Ventura, and J. Jatskevich, “Design forsurvival: Real-time infrastructure coordination,” in International Work-shop on Complex Network and Infrastructure Protection, Rome, 2006,CNIP.[63] J.A. Hollman, J.R. Mart´ı, J. Jatskevich, and K.D. Srivastava, “Dynamicislanding of critical infrastructure a suitable strategy to survive and mit-igate critical events,” in International Workshop on Complex Networkand Infrastructure Protection, Rome, 2006, CNIP.[64] Gene Golub and Charles Van Loan, Matrix Computation, Johns Hop-kins Series in the Mathematical Sciences. The Johns Hopkins UniversityPress, Baltimore, Maryland 21211, third edition, 1989.143Appendix AThe QR MethodThe QR method allows the calculation of eigenvalues of real matrices. TheStability of a numerical eigenvalue problem depends on the matrix underconsideration. If the matrix is symmetric with symmetrically distributederror, then the calculated eigenvalues will approximate the actual eigenvalues,provided the eigenvalues are all simple. Otherwise, the numerical methodsmy fail to find all eigenvalues.The basis of the QR method for calculating the eigenvalues of A is the factthat an n x n real matrix can be written asA = QR (A.1)where Q is orthogonal and R is upper triangular. The method is efficientfor the calculation of all eigenvalues of a matrix. The construction of Qand R proceeds as follows. Matrices P1, P2, ...Pn−1 are constructed so thatPn−1Pn−2....P2P1A = R is upper triangular. These matrices can be chosenas orthogonal matrices and are called Householder matrices. Since the P ′sare orthogonal, the stability of the eigenvalue problem will not be worsened.144Appendix A. The QR MethodIf we letQT = Pn−1Pn−2...P2P1 (A.2)then we have QTA = R andQQTA = QRIA = QRA = QRThe method can be summarized with the following steps,1. Set A1 = A, Q1 = Q, and R1 = R2. First set A2 = R1Q1; then factor A2 as A2 = Q2R2.3. First set A3 = R2Q2; then factor A3 as A3 = Q3R3.4. Set Am = Rm−1Qm−1; then factor Am−1 as Am = QmRm.Matrix Am will tend toward a triangular or nearly triangular form. Thus theeigenvalues of Am will be easy to compute.If the eigenvalues can be ordered as |λ1| > |λ2| > |λ3| > · · · > |λn| > 0, thenas m increases the eigenvalues of Am approach the eigenvalues of A.The QR method can be used for an arbitrary real matrix, but works muchfaster on special matrices:• symmetric tri-diagonal,145Appendix A. The QR Method• Hessenberg,• symmetric band matrices.It is important to note that if the elements of A have widely varying magni-tudes, then A should be balanced before applying the QR algorithm. WhenA is balanced, the computed eigenvalues are often more accurate [64].146IndexArgentina, 1Arnoldigeneralized method, 17modified method, 26Bilinear Transformation, 63Brazil, 1Calvin˜o, Jesu´s, xiv, 7CDA, 31, 39Chile, 1continues time state space, 82continuous time eigenvalues, 90discrete eigenvalues, 89, 90discrete time state space, 4, 40capacitance, 46inductance, 44MATE, 73parallel branches, 55resistance, 48rlc series, 85series branches, 50trapezoidal, 44, 46Dommel, Hermann, xiv, 4, 5, 42Edison, Thomas Alva, 1Eigenvalue, 16, 55Electricite de France, 6EMTP, ii, 5fast time-domain, 17Gaspard Richie, Baron de Prony,17Gissinger, S., 3Helmholtz, 44Hollman, Jorge, 7HVDC, 10Hydro Quebec Research Institute,6147IndexI2Sim, 133interarea oscillation, 15, 24Jardim, J.L., 4JIIRP, 132, 133Kundur, Prabha, xiv, 3Linares, Luis, xiv, 7LTI, 40Mart´ı, Jose´, xiv, 7Mitsubishi, 6Modified Arnoldi, 16Moshref, Ali, xivNSERC, 132OVNI, ii, 4, 7, 67MATE, 9OVNI-NET, 9PC-cluster, iiPereira, M., 3Prony, 16PSEPC, 132QR method, 144RTDNS, 7RTNS, 7small-signal stability, 15TEQSIM, 6Thevenin, 44Transient Network Analyzer, 5transition matrix, 89trapezoidal, 90148


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items