The Application of Shifted FrequencyAnalysis in Power System TransientStability StudiesbyAndrea T.J. Mart´ıB.A.Sc., The University of British Columbia, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)February 2018c© Andrea T.J. Mart´ı 2018AbstractPower system engineers use transient stability computer simulation programs tomodel the power grid’s behaviour when large disruptions cause the grid to deviatefrom its 60-Hz operating frequency. These programs must be able to capture thesefrequency dynamics around 60 Hz while being computationally efficient, as an exten-sive number of simulations are typically run for a given scenario. In the traditionalpower grid, the large mechanical inertia of the synchronous generators stabilizes thenetwork during disturbances and maintains the system frequency close to the 60 Hzoperating frequency. In the modern grid, however, the increase of renewable energysources lowers the grid’s inertia and larger frequency deviations can occur. Thephasor solution method employed in the traditional programs solves the networkassuming a constant 60-Hz frequency. When deviations from 60 Hz are prominent,the Electromagnetic Transients Program (EMTP) is used as an alternative to thephasor solution to capture these fluctuations. The EMTP models the electricalnetwork based on the differential equations of the network components, which al-lows the tracing of the network waveforms. However, this discretization requiressmall time-steps, which makes the solution method computationally expensive. TheShifted Frequency Analysis (SFA) method discussed in this work is an alternativeto the traditional phasor solution and to the EMTP solution. In this work, a gener-alized SFA-based program is written and used for transient stability analysis. SFAis a discrete-time solution method, like the EMTP, but uses a frequency-shiftingtransformation to bring the solution domain down to 0 Hz. Because of this trans-formation, SFA can capture network dynamics around 60 Hz using large time-steps,making it suitable for transient stability analysis studies.iiLay SummaryPower system engineers use transient stability computer simulation programs tomodel the power grid’s behaviour when large disruptions cause the grid to deviatefrom its 60 Hz operating frequency. These programs must be able to capture the fre-quency dynamics around 60 Hz while being computationally efficient, as an extensivenumber of simulations are typically run for a given scenario. With the increase of re-newable energy sources connected to the grid, the power system is more susceptibleto network changes than the traditional power grid, whose high-inertial synchronousgenerators kept the grid close to 60 Hz. The traditional simulation programs usethe constant 60 Hz frequency in their models and the alternative programs who cancapture frequency dynamics can be computationally slow for the 60-Hz fluctuations.This thesis proposes the solution method, Shifted Frequency Analysis (SFA), as amodel to capture frequency dynamics around 60 Hz computationally efficiently forpower system transient stability studies.iiiPrefaceThe research work in this thesis to the best of my knowledge is original and unpub-lished, except where references are made to previous works. The work is done bythe author under the supervision of Prof. Juri Jatskevich.The research involved in Chapters 4, 6 and 7 has been accepted in the Power SystemsComputation Conference (PSCC): A. Mart´ı, and J. Jatskevich, “Transient StabilityAnalysis using Shifted Frequency Analysis (SFA),” in Power Systems ComputationConference., Dublin, Ireland, June 2018.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 The Role of AC Synchronization and Power Equilibrium For PowerGrid Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Transient Solution Methods Used in Transient Stability Studies . . . 61.3.1 The Traditional Phasor Solution Method . . . . . . . . . . . 61.3.2 The EMTP Solution Method . . . . . . . . . . . . . . . . . . 91.3.3 Time-varying Phasors and the Introduction of the Shifted Fre-quency Analysis Solution Method . . . . . . . . . . . . . . . . 111.4 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.5 Research Contributions . . . . . . . . . . . . . . . . . . . . . . . . . 181.6 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18vTABLE OF CONTENTS2 Shifted Frequency Analysis: A Discrete-Time Solution Methodand Solution Framework . . . . . . . . . . . . . . . . . . . . . . . . . . 202.1 Shifted Frequency Analysis Solution Framework . . . . . . . . . . . . 202.1.1 A Single-Frequency Signal . . . . . . . . . . . . . . . . . . . . 212.1.2 The Time-Frequency Relationship . . . . . . . . . . . . . . . 212.1.3 The SFA Solution Framework . . . . . . . . . . . . . . . . . . 212.2 The SFA Solution Method . . . . . . . . . . . . . . . . . . . . . . . . 232.2.1 The Choice of Time-Step in the SFA Solution Method . . . . 262.2.2 A Summary of the SFA Solution Method . . . . . . . . . . . 262.3 The SFA Equivalent Network Elements . . . . . . . . . . . . . . . . 272.4 The Discretization of the SFA Network Elements . . . . . . . . . . . 292.4.1 Equivalent SFA Network Elements Using the Trapezoidal Dis-cretization Rule . . . . . . . . . . . . . . . . . . . . . . . . . . 302.4.2 Equivalent SFA Network Elements Using the Backward EulerDiscretization Rule . . . . . . . . . . . . . . . . . . . . . . . . 322.5 A Comparison of the Traditional Phasor, the EMTP, and the SFANetwork Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.6 An RL Circuit Energized with an AC Source Using the TraditionalPhasor, the EMTP and SFA . . . . . . . . . . . . . . . . . . . . . . . 363 The Numerical Discretization Methods Used in SFA . . . . . . . . 393.1 The Evaluation of a Numerical Discretization Method Using the Fre-quency Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2 The System’s Frequency Response in the Physical Domain and in theSFA Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.2.1 The Frequency Response in the Physical Domain . . . . . . . 403.2.2 The Frequency Response in the SFA Domain . . . . . . . . . 423.2.3 The Continuous-Time Frequency Response in the SFA Domain 423.3 The Frequency Response of the Trapezoidal Discretization Rule inthe SFA Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3.1 Approach One: Applying the Single-Frequency Exponential . 433.3.2 Approach Two: Applying the Z-Domain Transfer Function . 443.3.3 The Trapezoidal Rule’s Accuracy in the SFA Domain . . . . 463.3.4 The Numerical Accuracy of the Trapezoidal Rule: Magnitudeand Phase Responses . . . . . . . . . . . . . . . . . . . . . . . 473.3.5 The Distortion of the Inductor Using the Trapezoidal Rule . 50viTABLE OF CONTENTS3.4 The Frequency Response of the Backward Euler Discretization Rulein the SFA Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4.1 Applying the Z-Domain Transfer Function . . . . . . . . . . . 513.4.2 The Backward Euler Rule’s Accuracy in the SFA Domain . . 523.4.3 The Numerical Accuracy of the Backward Euler Rule: Mag-nitude and Phase Responses . . . . . . . . . . . . . . . . . . . 533.4.4 The Distortion of the Inductor Using the Backward Euler Rule 543.4.5 A Comparison of The Trapezoidal and Backward Euler RuleAround 0 Hz (DC) . . . . . . . . . . . . . . . . . . . . . . . . 593.4.6 A Summary of the Accuracy of the Two Discretization Rules 603.5 The Discretization Rule’s Stability in the SFA Domain . . . . . . . . 613.6 The Response to a Step Input in the SFA Domain . . . . . . . . . . 623.7 An RL Circuit Energized with a DC Source Using the TraditionalPhasor, the EMTP and SFA . . . . . . . . . . . . . . . . . . . . . . . 694 Transient Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . 714.1 A Physical Representation . . . . . . . . . . . . . . . . . . . . . . . . 724.1.1 The Equal Area Criteria and The Critical Clearing Time . . 734.1.2 The Approximated Synchronous Machine Model . . . . . . . 744.2 A Mathematical Representation . . . . . . . . . . . . . . . . . . . . . 744.2.1 The Coupling of the Electromechanical and Electrical Net-work Equations in the Solution Methods . . . . . . . . . . . . 764.2.2 The Discretization of the Electromechanical Equation . . . . 774.3 Implementation of SFA in Three Transient Stability Case Studies . . 784.3.1 Fault Current Zero-Crossing Detection Algorithm in TSFA . 825 Single-Generator Infinite-Bus Case Study . . . . . . . . . . . . . . . 855.1 Critical Clearing Times . . . . . . . . . . . . . . . . . . . . . . . . . 865.2 A Three-Phase Fault at the Generator’s Terminals . . . . . . . . . . 875.2.1 A Comparison of Using the Actual Frequency with the Ap-proximate Frequency . . . . . . . . . . . . . . . . . . . . . . . 935.2.2 A Comparison of the Discretization Methods in the TSFASolution During the Three-Phase Fault. . . . . . . . . . . . . 946 Three-Bus Network Case Study . . . . . . . . . . . . . . . . . . . . . 966.1 Critical Clearing Times . . . . . . . . . . . . . . . . . . . . . . . . . 98viiTABLE OF CONTENTS6.2 A Three-Phase Fault at the Terminals of Transformer Two . . . . . 1006.2.1 A Comparison of Using the Actual Velocity with the Approx-imate Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . 1046.2.2 A Comparison of the Discretization Rules in the TSFA Solu-tion During the Three-Phase Fault. . . . . . . . . . . . . . . . 1066.2.3 The DC Offset Captured in the Fault Current in SFA . . . . 1086.3 Generator One Disconnected From the Network . . . . . . . . . . . . 1097 Thirty-Nine Bus Network Case Study . . . . . . . . . . . . . . . . . 1117.1 Critical Clearing Times . . . . . . . . . . . . . . . . . . . . . . . . . 1127.2 A Three-Phase Fault Applied to Bus 16 . . . . . . . . . . . . . . . . 1157.3 Generator Nine Disconnected from the Network . . . . . . . . . . . . 1198 Conclusion and Future Research . . . . . . . . . . . . . . . . . . . . . 1208.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1208.2 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131viiiList of TablesTable 2.1 The branch equivalents for the traditional phasor, the EMTP, andthe SFA methods using the trapezoidal discretization rule for theSFA and the EMTP. . . . . . . . . . . . . . . . . . . . . . . . . . 35Table 2.2 The branch equivalents for the traditional phasor, the EMTP, andthe SFA methods using the backward Euler discretization rule forthe EMTP and SFA. . . . . . . . . . . . . . . . . . . . . . . . . . 35Table 3.1 A comparison of the EMTP and SFA frequency response accuracyratio’s when the trapezoidal rule is behaving as an integrator . . 47Table 3.2 The discontinuities at the Nyquist frequency of different time-stepswhen using the trapezoidal rule. . . . . . . . . . . . . . . . . . . . 47Table 3.3 A comparison of the EMTP’s and SFA’s distortion of the inductorwhen using the trapezoidal rule. . . . . . . . . . . . . . . . . . . . 51Table 3.4 A comparison of the EMTP and SFA frequency response accuracyratio for the backward Euler rule behaving as an integrator. . . . 53Table 3.5 A comparison of the EMTP’s and SFA’s distortion of the inductorwhen using the backward Euler rule. . . . . . . . . . . . . . . . . 59Table 3.6 The response of a step input using the trapezoidal rule, comparingthe SFA solution with the EMTP solution. . . . . . . . . . . . . . 63Table 3.7 The response of a step input using the backward Euler rule, com-paring the SFA solution with the EMTP. . . . . . . . . . . . . . . 63Table 5.1 Network parameters for the SGIB network (Figure 5.1). . . . . . 85Table 5.2 The CCTs for the three-phase fault on the generator’s terminalsof the SGIB system in Figure 5.1 in 60-Hz cycles and milliseconds,and the critical clearing angle in degrees. . . . . . . . . . . . . . . 86Table 5.3 A comparison of the error in frequency among TSAT (at ∆t =8 ms)and TSFA (at ∆t =5 ms and ∆t =8 ms) with respect to the EMTP. 89ixLIST OF TABLESTable 5.4 A comparison of the error in frequency between TSAT, TSAT*,TPhasor and TPhasor* with respect to the EMTP reference . . . 93Table 6.1 Network parameters for the three-bus network (Figure 6.1) . . . . 97Table 6.2 The CCTs for three-phase faults at different locations of the three-bus system in Figure 6.1 in 60-Hz cycles and the critical clearingangle in degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99Table 6.3 A comparison of the error in frequency between TSAT and TSFAwith respect to the EMTP reference solution. . . . . . . . . . . . 100Table 6.4 A comparison of the error in frequency between the different TSFAsolutions with respect to the EMTP reference solution. . . . . . . 105Table 6.5 The four different conditions for the phasor solution in Figure 6.8. 105Table 6.6 A comparison of the error in frequency between the different TPha-sor solutions with respect to the EMTP reference solution. . . . . 106Table 7.1 The CCTs for three-phase faults at different locations of the three-bus system in Figure 7.1 in 60-Hz cycles and the critical clearingangle in degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113Table 7.2 A comparison of the error in frequency between TSAT and TSFAwith respect to the EMTP reference solution. . . . . . . . . . . . 117Table 7.3 A comparison of the error in frequency between the different TSFAsolutions with respect to the EMTP reference solution. . . . . . . 117Table A1 Generator data for the IEEE 39-bus case study (Modified from [1]).131Table A2 Transformer data for the IEEE 39-bus case study. . . . . . . . . . 131Table A3 Power flow results for the IEEE 39-bus case study. . . . . . . . . 132Table A4 Transmission line data for the IEEE 39-bus case study. . . . . . . 133xList of FiguresFigure 1.1 A classical equivalent representation of the power system showsthat the electrical power delivered depends on the voltages, anglesand network reactance. . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.2 AC power transmission limits for (a) steady-state conditions and(b) transient conditions. . . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.1 The process to create a complex signal with frequency side-bandsaround 0 Hz from a sinusoidal signal with frequency side-bandsaround ωo, such that a sinusoidal signal rotating around 60 Hzcan now rotate around 0 Hz. . . . . . . . . . . . . . . . . . . . . . 24Figure 2.2 The SFA solution method: a real continuous-time signal is shiftedto the SFA world, discretized and solved for in the discrete-timeSFA domain, interpolated into a continuous-time signal in the SFAdomain, and reverse transformed back to the physical domain.The imaginary portion is dropped and the real part of the signalis kept to create the physical continuous-time solution. . . . . . . 25Figure 2.3 Numerical discretization with the trapezoidal rule, where the areaunder the curve is Area = v(t)+v(t−∆t)2 (∆t). . . . . . . . . . . . . . 30Figure 2.4 Numerical discretization with the backward Euler rule, where thearea under the curve is Area = v(t)∆t. . . . . . . . . . . . . . . . 32Figure 2.5 Four representations of an RL circuit energized by an AC source:(a) the time-domain, (b) the phasor domain, (c) the EMTP, and(d) the SFA domain. . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 2.6 The node voltages and current through the inductor of the RLcircuit in Figure 2.5 comparing the traditional phasor, the EMTP,and the SFA solution methods. . . . . . . . . . . . . . . . . . . . 38xiLIST OF FIGURESFigure 3.1 An LTI system’s response to an input signal that is the single-frequency exponential [2]. . . . . . . . . . . . . . . . . . . . . . . 41Figure 3.2 The frequency responses of an inductor and a capacitor whosevalues are L=1 and C=1 are the frequency responses of a differ-entiator and an integrator when current is the input and voltageis the output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.3 The accuracy of the trapezoidal rule acting as a differentiator(magnitude and phase response). . . . . . . . . . . . . . . . . . . 48Figure 3.4 The accuracy of the trapezoidal rule acting as an integrator (mag-nitude and phase response). . . . . . . . . . . . . . . . . . . . . . 49Figure 3.5 The accuracy of the trapezoidal rule acting as a differentiator(around 0% error, which occurs when∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1). . . . . . . 50Figure 3.6 The accuracy of the backward Euler rule acting as a differentiator(magnitude and phase response). . . . . . . . . . . . . . . . . . . 55Figure 3.7 The accuracy of the backward Euler rule acting as an integrator(magnitude and phase response). . . . . . . . . . . . . . . . . . . 56Figure 3.8 The accuracy of the backward Euler rule acting as a differentiator(around 0% error, which occurs when∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1). . . . . . . 57Figure 3.9 A comparison of the accuracy of the trapezoidal and backwardEuler rules acting as an integrator around 0 Hz. . . . . . . . . . 59Figure 3.10 A comparison of the magnitude of the accuracy of the trapezoidaland backward Euler frequency responses acting as a differentiatoraround 60 Hz using ∆t = 8 ms and ∆t = 10 ms. . . . . . . . . . . 61Figure 3.11 A DC current source of I=1 A energizes an inductor of L=1 mH. . 64Figure 3.12 The voltage across an inductor in response to a 1A DC current inthe EMTP: (a) the response given when using the trapezoidal ruleand (b) the response given when using the backward Euler rule. . 65Figure 3.13 The voltage across an inductor in response to a 1A DC currentusing the trapezoidal rule in SFA: (a) the response given in theshifted domain and (b) the response given in the physical domain. 67Figure 3.14 The voltage across an inductor in response to a 1A DC currentusing the backward Euler rule in SFA: (a) the response given in theshifted domain and (b) the response given in the physical domain. 68Figure 3.15 An RL circuit energized with a DC voltage source. . . . . . . . . 69xiiLIST OF FIGURESFigure 3.16 The RL circuit discretized with the EMTP (left) and the SFA(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 3.17 A comparison of the trapezoidal and backward Euler rules in re-sponse to a DC source, using the EMTP as the reference solution. 70Figure 4.1 The rotor angle of the generator is the relative angle between thestator winding and the rotor’s field winding, taking the stator asthe reference frame and using the simplified model of the syn-chronous generator (Section 4.1.2). . . . . . . . . . . . . . . . . . 72Figure 4.2 The TSFA solution algorithm. . . . . . . . . . . . . . . . . . . . . 84Figure 5.1 The single-line diagram of the Single-Generator Infinite-Bus (SGIB)network, where the electrical network is represented as an infinitebus. A three-phase balanced fault is applied to the terminals ofthe generator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85Figure 5.2 A comparison of the TSFA and TSAT solutions compared with theEMTP solution for the electromechanical variables in the systemof Figure 5.1 during a three-phase fault at the terminals of thegenerator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Figure 5.3 A comparison of the EMTP, TSFA using ∆t = 5 ms and ∆t =8 ms and the phasor solution using ∆t = 5 ms and ∆t = 8 ms forthe electromechanical variables during the three-phase fault in thesystem of Figure 5.1. . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 5.4 The voltage at the faulted bus (terminals of Generator 1) duringthe 3-phase fault. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 5.5 The fault current during a 3-phase fault at the terminals of Gener-ator 1 (Top: using ∆t = 8 ms for TSFA, Bottom: using ∆t = 1 msfor TSFA). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 5.6 A comparison of using the grid frequency, ω, with the approximatefrequency, ωo, in the SFA solution and in the phasor solution com-pared to the EMTP for the three-phase fault on the SGIB. . . . . 93Figure 5.7 A comparison of the trapezoidal and backward Euler discretiza-tion rules on the rotor angle, frequency, electrical power, and faultvoltage for the three-phase fault condition. Top left: Generator’srotor angle, Top right: Generator’s frequency, Bottom left: Gen-erator’s electrical power, Bottom right: Fault voltage. . . . . . . . 95xiiiLIST OF FIGURESFigure 6.1 The one-line diagram of the three-bus network, where two gener-ators feed a constant impedance load. . . . . . . . . . . . . . . . . 96Figure 6.2 A comparison of the EMTP, TSFA, and TSAT solutions for theelectromechanical variables of Generator 1 in the system of Figure6.1 during a three-phase fault at the terminals of Transformer 2. 101Figure 6.3 A comparison of the EMTP, TSFA, and TSAT solutions for theelectromechanical variables of Generator 2 in the system of Figure6.1 during a three-phase fault at the terminals of Transformer 2. 102Figure 6.4 The voltage at the faulted bus during a 3-phase fault at the ter-minals of Transformer 2. . . . . . . . . . . . . . . . . . . . . . . . 103Figure 6.5 The fault current during a 3-phase fault at the terminals of Trans-former 2(using a ∆t = 8 ms for the SFA solution). . . . . . . . . . . . . . 103Figure 6.6 The fault current during a 3-phase fault at the terminals of Trans-former 2(using a ∆t = 1 ms for the SFA solution). . . . . . . . . . . . . . 103Figure 6.7 A comparison of using the grid frequency, ω, with the approxi-mate frequency, ωo, in the SFA solution in the electromechanicalequation for the three-phase fault. . . . . . . . . . . . . . . . . . . 104Figure 6.8 A comparison of using the grid frequency, ω, with the approximatefrequency, ωo, in the phasor solution for the three-phase fault. . . 106Figure 6.9 A comparison of the trapezoidal and backward Euler discretiza-tion rules inSFA on the rotor angle, frequency, electrical power, and fault volt-age during the three-phase fault at the terminals of Transformer2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure 6.10 A three-phase fault is cleared when the voltage crosses zero tocapture the maximum DC offset. The SFA solution, using boththe trapezoidal and backward Euler rules, and TPhasor* solutionare compared to the EMTP. . . . . . . . . . . . . . . . . . . . . . 108Figure 6.11 The frequency of Generator 2 after disconnecting Generator 1 fromthe three-bus network in Figure 6.1 using the backward Euler rule(top) and both rules (bottom). . . . . . . . . . . . . . . . . . . . 110Figure 6.12 The insert of Figure 6.11. Left: without the trapezoidal rule,Right: with the trapezoidal rule. . . . . . . . . . . . . . . . . . . 110xivLIST OF FIGURESFigure 6.13 A comparison of the traditional phasor solution using the correctedelectromechanical equation. . . . . . . . . . . . . . . . . . . . . . 110Figure 7.1 The IEEE 39-bus test system. . . . . . . . . . . . . . . . . . . . . 111Figure 7.2 A comparison of the EMTP, TSFA, and TSAT solutions for theelectromechanical variables of Generator 5 in the system of Figure7.1 during a three-phase fault on Bus 16. . . . . . . . . . . . . . . 116Figure 7.3 A comparison of the SFA solution using the actual frequency ωinstead of the approximated frequency ωo. . . . . . . . . . . . . . 117Figure 7.4 The voltage at the faulted bus during a 3-phase fault at Bus 16 inFigure 7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118Figure 7.5 The fault current during a 3-phase fault at Bus 16 in Figure 7.1. 118Figure 7.6 The frequency of Generator 5 after losing Generator 9 in the 39-bus Network of Figure 7.1 using the backward Euler rule (top) andboth the backward Euler and trapezoidal rules (bottom). Usingthe trapezoidal rule gives numerical oscillations. . . . . . . . . . . 119xvList of AbbreviationsCCT Critical Clearing TimeEMTP Electromagnetic Transients Programp.u. Per UnitSFA Shifted Frequency AnalysisSGIB Single-Generator Infinite-BusTPhasor Transient PhasorTSAT Transient Security Assessment ToolTSFA Transient Shifted Frequency AnalysisxviNomenclatureIn this thesis, three types of representations for the electrical variables are used:instantaneous time-varying scalars, steady-state phasors, and time-varying phasors.To differentiate among the representations, the following nomenclature is used (un-less used in another format in the literature):• Instantaneous time-varying scalars are written in italics. E.g., v(t) and i(t).• Phasors are written in italics with a bar, but independent of time. E.g., V¯and I¯• Time-varying phasors are written in italics with a bar, as functions of time.E.g., V¯ (t) and I¯(t)xviiAcknowledgementsI would like to express my deepest gratitude to my supervisor Prof. Juri Jatskevich.His gentle kindness and patience during my journey in the graduate program allowedme to not only grow in technical knowledge, but to mature as a person. His tech-nical advice and guidance in this research allowed the ideas to strengthen and bedynamic. His ingrained teachings to take things one step at a time is a life-longlesson. I am grateful for his kindness, not only in this research, but in life.I would also like to express my gratitude to Prof. KD Srivastava, whose wisdomguided me at the beginning of my journey, teaching me to reach out to all areasin the world and society beyond the scope of engineering. Our long conversationswould leave me in dumbfounded awe and with the notion that even as a small entity,I have the commitment and desire to do as much good as possible in this world.I am thankful to my friends in the power group for every conversation we have hadover the years has opened my eyes to new beliefs and understandings about our rolein the community, and our purpose in this world. I am grateful for your kindnessand deep desires to respect knowledge and freedom. Thank you for teaching methat all ideas and feelings are integral.To my Papa, Mama and Catherine. The values they have instilled in me throughtheir own lives have taught me more than anything else. Their strengths to overcomeany obstacles, optimisms in humanity and quests for the truth resonant throughouteach of their lives, providing me with a little spark of courage to live my own lifeas such. Their unwavering patience, wisdom and love for me and their belief in menever falter and each time I fall, they pick me up with trust that I will go fromcrawling to walking a little farther each time. Thank you for imagining far beyondwhat is physically seen.xviiiDedicationTo my Papa, Mama and Catherine. You are my reason for being.xixChapter 1Introduction1.1 Research MotivationFrequency fluctuations around the grid frequency are prominent on the electricpower grid when large disturbances perturb the network. Traditional transient sta-bility simulation programs assume that these fluctuations are close to grid frequency,fo, which is 60 Hz in North America and the solution methods used in these pro-grams use the constant 60-Hz frequency when modelling the electric network [3].This assumption was acceptable when modelling the traditional power grid, whereheavy, high-inertial synchronous generators convert resources, such as coal, oil, hy-dro and natural gas, into electricity. Physically, these generators store kinetic energyin their rotating components, which is proportional to their mass and their frequencyof rotationEk =12Iω2, (1.1)where Ek is the amount of stored kinetic energy (J), I is the generator’s inertiaconstant (kg·m2) and ω is the frequency of rotation ( rads ), which in this case is thegrid frequency: ω = ωo = 2pifo [4]. The stored kinetic energy contributes to thegenerator’s ability to react and to adapt to external forces. The heavy hydroelectricand thermoelectric generators store large amounts of kinetic energy and have a largeconstant of inertia, which gives them the ability to remain relatively stable duringdisturbances. When there are no disturbances on the grid, the generators rotatein synchronicity with each other at ωo. When a disturbance occurs on the grid,the generators will deviate from their steady-state frequency, or velocity, but areinertially reluctant to do so, as mathematically captured in Newton’s second law forrotational bodiesTnet = Iα, (1.2)11.1. Research Motivationwhich, for power systems, is used in the following form [4]Tnet = Id2δdt2, (1.3a)Tnet = Idωdt, (1.3b)Pnet = ωoIdωdt. (1.3c)where in (1.3a), I is the machine’s inertia constant (kg·m2), δ is the machine’s rotorangle (radians), Tnet is the net difference between the mechanical input torque andthe electrical output torque (Nm) and in (1.3b), ω is the machine’s velocity ( rads ).In (1.3c), net torque is converted to net power using the grid frequency ωo, as itis net power that is measured in the power network. These equations show thatwhen there is a power imbalance, the rate of frequency change is determined by themachine’s inertia constant. Due to the large inertia in these synchronous genera-tors, when such a disturbance changes the power balance, the frequency can deviate≈ ±0.1-0.5 Hz from the grid frequency[5]. Even so, larger frequency excursions fromthe grid frequency can occur, such as in the 2003 Italian blackout, when the gridfrequency decreased 2.5 Hz, from 50 Hz to below 47.5 Hz, before the frequency con-trollers interfered [6].On the modern grid, there is an increase of renewable energy plants, such as windfarms and solar photovoltaic systems contributing to power generation. The powerelectronic controllers, which couple these renewable energy plants to the grid, haveno rotational motion, and therefore, no rotational inertia [7][8]. The lack of iner-tia from these plants lowers the grid’s overall inertia and now, when a disturbanceperturbs the system, the frequency deviations will be greater than on the tradi-tional high-inertial grid. For example, in the South Australian grid, where 48.4% ofgeneration comes from wind and rooftop PV [9], during the 2016 South Australianblackout, there was a frequency drop of 3 Hz, from the 50 Hz grid frequency to below47 Hz [10]. In addition to the grid’s decrease in inertia from the renewable sourcesreplacing the high-inertial generators, certain power grid codes mandate that belowa given frequency, the renewable energy source must be completely disconnectedfrom the grid. For example, in the EDF Island Energy System, the renewable en-ergy sources, such as wind plants and PV arrays, are required to disconnect whenthe frequency falls below a set value, in one case, 49.5 Hz [11]. The disconnection21.2. The Role of AC Synchronization and Power Equilibrium For Power Grid Stabilityof these renewable energy sources further decreases the grid’s inertia, increasingits susceptibility to frequency deviations and decreasing its ability to remain stableduring large disturbances.As the share of renewable energy sources is expected to provide 40% of the world’stotal electricity by 2040 [12], the solution method in the simulation programs usedto model the grid’s behaviour must capture the frequency deviations. However,the traditional transient stability tools that assume the constant 60-Hz frequencyin their models [3][13] are not able to accurately model these fluctuations and thealternative solution programs that can capture these frequency deviations [14][15]are computationally slower than the constant frequency solution methods. There isa need to develop new solution methods for transient stability tools that can cap-ture frequency deviations around the 60-Hz frequency while maintaining the highcomputational speed of the traditional tools.This thesis proposes the Shifted Frequency Analysis (SFA) solution method fortransient stability studies. SFA is able to model these grid frequency fluctuationswithout sacrificing computational time. This thesis is motivated to determine if SFAcan be an alternative solution method to the traditional methods used in currenttransient stability tools.1.2 The Role of AC Synchronization and PowerEquilibrium For Power Grid StabilityThis section provides a brief overview of two factors integral to the AC power grid’sstability: the synchronization of the generators and the network’s power equilibrium.SynchronizationThe AC power grid’s stability depends on the generators rotating in synchroniza-tion with each other. One concern with large frequency deviations is that theycreate a net power imbalance (1.3c). In balanced steady-state operating conditions,Pout = Pin or Pmech = Pelec +Plosses and the generators rotate and produce AC elec-tricity at the grid frequency fo. Sudden disturbances, such as large loads connecting,generators disconnecting, or faults occurring, change the net power balance. Whenthe power is unbalanced, each of the generators releases or intakes more internal31.2. The Role of AC Synchronization and Power Equilibrium For Power Grid Stabilitykinetic energy (1.1) and (1.3c), ↑↓ Tnet → ↑↓ Ek, moving their frequency awayfrom fo and “out of step”[4] with each other. This asynchronicity impacts the grid’sstability. One primary concern of power system stability studies is the grid’s abilityto maintain synchronicity after such sudden disturbances perturb the system [4].StabilityStability is “the quality of the system or part of the system which enables it todevelop restoring forces between these elements equal or greater than the disturbingforces so as to restore a state of equilibrium” [16]. In AC power systems, this stateof equilibrium refers to the synchronous machines remaining in synchronism witheach other [4]. There is a relationship between the maximum power transmittedand the generator’s rotor angle position. A classical equivalent representation of thepower system [4] is given in Figure 1.1 . eP eqX gen relE load0E Figure 1.1: A classical equivalent representation of the power system shows that the electricalpower delivered depends on the voltages, angles and network reactance.In this representation, the generator delivers electrical power Pe (MW) to the loadthrough the network’s equivalent reactance Xeq (Ω). The voltages, Egen and Eload,are the voltages of the generator and load (V) and δrel is the angle between thegenerator’s rotor and stator windings (radians or degrees). The voltage at the loadis taken as the reference bus, with an angle of 0◦. The AC power transmission ofFigure 1.1 is mathematically represented in (1.4)Pe =EgenEloadXeqsin(δrel). (1.4)41.2. The Role of AC Synchronization and Power Equilibrium For Power Grid StabilityMaximum power is transferred when δrel = 90◦, (Pmax = E1E2Xeq ) and this is calledthe system’s “steady-state stability limit” [4] because any extra power transmittedwill cause the rotor to move out of synchronicity (Figure (1.2(a))).After a disturbance, the grid is no longer in steady-state and its stability dependson the imbalance of power in the system. For example, in the case of a fault,the electrical power decreases and Pmech > Pelec. The excess mechanical energyis transferred to rotational kinetic energy and the rotor speeds up, increasing itsrelative angle away from its steady-state position (↑ δrel). If the fault is cleared,but the rotor has accelerated past a critical point in which the generator has gainedmore kinetic energy than it can release, the system is unstable. The equal-area curve,shown in Figure (1.2(b)), depicts this energy balance. In Figure (1.2(b)), Area 1,from δ1 → δ2, is the amount of kinetic energy that the generator has accumulatedas it accelerates, where δ2 is the rotor angle at the time the fault is cleared. Area2, from δ2 → δ3 is the amount of kinetic energy that must be borrowed from therotational masses so that the rotor can decelerate. When the two areas are equal,the system is stable. The transient stability limit occurs when the rotor angle isat a point where the two areas are equal. This angle is called the critical clearingangle, δc. If the fault is cleared after the rotor angle has moved beyond δc then thesystem is unstable [4]. If the fault is cleared before the rotor angle has reached δc,the system is stable.1oArea 1 whereArea 2 where0 eP e mP P2m echP0 3 o180 eP21Area 1 where 0eP Area 2 where e mechP P3 1180 0mechPelec-maxPePePmaxP0max 90 180(a) The steady-state stabilitypower-angle curve shows that Pmaxoccurs when δrel = 90◦. This is thesteady-state stability limit.1oArea 1 whereArea 2 where0 eP e mP P2m echP0 3 o180 eP21Area 1 where 0eP Area 2 where e mechP P3 1180 0mechPelec-maxPePePmaxP0max 90 180(b) The transient stability power-areacurve is used to determine the tran-sient stability limit, which is whenArea 1 = Area 2.Figure 1.2: AC power transmission limits for (a) steady-state conditions and (b) transientconditions.51.3. Transient Solution Methods Used in Transient Stability Studies1.3 Transient Solution Methods Used in TransientStability StudiesTransient stability simulation tools are used to calculate the transient stability lim-its before the system becomes unstable. These programs depend on mathematicalsolution methods which model the power grid. The validity of the given solutionmethod determines how close the model is to the physical network. This sectionpresents an overview of two solution methods used to model the electrical network inpower system dynamic studies: the phasor solution method and the ElectromagneticTransients Program (EMTP), and introduces the proposed time-domain solutionmethod, SFA. Section 1.3.1 also provides a brief introduction to the electromechan-ical equation used to model the electromechanical dynamics, however Chapter 4gives a more thorough discussion of the coupling between the electromechanical andelectrical network dynamics.1.3.1 The Traditional Phasor Solution MethodBackgroundThe birth of calculating transient stability limits came in the early 20th century,when the introduction of automatic voltage regulators led to the increased develop-ment of transmission lines with high impedances. These lines allowed large amountsof power to be transmitted closer to the steady-state stability limit [4][16]. Con-sequently, the system became more susceptible to disturbances and instability fre-quently occurred [4]. The engineers saw a need to develop analytical methods thatwould determine power transfer limits such that the system would remain stable.These limits were called transient stability limits [16].To determine these transient stability limits, the engineers modelled and analysedthe power system using AC network analysers. The analysers would compute thepower and rotor angle for off-line network contingencies at step-by-step time inter-vals [4]. From these results, time-angle curves called “swing-curves” would be drawnfor the different contingencies. As faults were one of the most critical contingencies[4], these curves would determine the critical rotor angle, as explained in Section1.2. The time that the critical angle occurred at is a critical parameter in transientstability studies as it determines the maximum time that the physical circuit break-ers have to open and successfully clear the fault, called the critical clearing time.61.3. Transient Solution Methods Used in Transient Stability StudiesTo simplify the mathematical equations, an approximation was made to the mainelectromechanical equation that governs the machine’s dynamics (1.3a). Becausethe rotor angle and velocity determined in (1.3a)(1.3b) are used in the electricalnetwork solution to calculate electrical power, which couples the electrical networkand electromechanical dynamics; and it is electrical power, not torque, that is theoutput of the electrical network solution, the net torque in (1.3a) is converted tonet power through frequency ω. However, using ω makes the equations non-linearand for simplification, the constant ωo is used in the conversionPnet = Tnet · ωo and Tnet = I dωdt→ Pnet = ωoI dωdt. (1.5)To model the electrical network, the traditional phasor solution method was used.Phasor Solution MethodIn the phasor solution method, a sinusoidal time-domain signal is frequency trans-formed into a frequency-domain signal. The sinusoidal signal is represented with itsEuler representation [17]v(t) = Vˆ cos(ωot+ δv)= R(Vˆ ej(ωot+δv))= R(Vˆ ejωotejδv)(1.6a)= R(V¯ ejωot).i(t) = Iˆ cos(ωot+ δi)= R(Iˆej(ωot+δi))= R(Iˆejωotejδi)(1.6b)= R(I¯ejωot).where V¯ = Vˆ ejδv and I¯ = Iˆejδi are the “phasor-representations” of the sinusoidalsignal. The phasor is independent of time and rotates at a constant frequency ωowith a magnitude and an angle (V¯ = Vˆ ∠δv, I¯ = Iˆ∠δi).Equations (1.6a) and (1.6b) are applied to the time-instantaneous voltages and cur-rents, v(t) and i(t), in the differential equations that govern the network elementsvL(t) = LddtiL(t) (1.7a) iC(t) = CddtvC(t) (1.7b)71.3. Transient Solution Methods Used in Transient Stability Studiesto obtain their frequency-domain representationsvL(t) = LddtiL(t)R(V¯ ejωot)= Lddt(R(I¯ejωot))R(V¯ ejωot)= jωoL(R(I¯ejωot)) (1.8a)V¯ = jωoLI¯ .iC(t) = CddtvC(t)R(I¯ejωot)= Cddt(R(V¯ ejωot))R(I¯ejωot)= jωoC(R(V¯ ejωot)) (1.8b)I¯ = jωoCV¯ .In (1.8a) and (1.8b), the operator ddt becomes jω in the frequency domain becauseddt(ejωt) = jω(ejωt), (1.9)and the inductors and capacitors are represented as constant impedances and ad-mittances rotating at ωoZL = jωoL. (1.10a) YC = jωoC. (1.10b)The transformation from time to frequency is the basis of the phasor solution methodand as seen in (1.10a) and (1.10b), the constant frequency ωo is used in this method.Constant Frequency AssumptionUsing the constant frequency ωo for the electrical network model (1.8a)(1.8b) isan approximation to the actual grid frequency ω. If the frequency of the powernetwork remains close to synchronous frequency, using the constant frequency inthe phasor solution and in the electromechanical equation (1.5) helped to simplifythe mathematical calculations in the early transient stability programs. The phasorsolution method for modelling the electrical network and the approximate-frequencyelectromechanical equation are implemented in the current computer simulationprograms used for transient stability analysis [3][13][18].81.3. Transient Solution Methods Used in Transient Stability Studies1.3.2 The EMTP Solution MethodIn the 1960s, with the advent of digital computers, H. W. Dommel developed theElectromagnetic Transients Program (EMTP) [19]. The EMTP is a discrete, time-domain solution method, where the continuous-time differential equations (1.7a)and (1.7b) are replaced with difference equations through a numerical discretizationmethod. For example, (1.11) derives the discretization of an inductor using thetrapezoidal discretization rule and (1.12) derives the discretization of a capacitorusing the trapezoidal rule [20].vL(t) = LddtiL(t)∫ tt−∆tvL(t)dt = L[iL(t)− iL(t−∆t)]vL(t) + vL(t−∆t)2∆t = L[iL(t)− iL(t−∆t)]vL(t) = (2L∆t)iL(t) + [−vL(t−∆t)− 2L∆tiL(t−∆t)]vL(t) = (2L∆t)iL(t) + ehL(t),(1.11)andiC(t) = CddtvC(t)∫ tt−∆tiC(t)dt = C[vC(t)− vC(t−∆t)]iC(t) + iC(t−∆t)2∆t = C[vC(t)− vC(t−∆t)]vC(t) = (∆t2C)iC(t) + [vC(t−∆t) + ∆t2CiC(t−∆t)]vC(t) = (∆t2C)iC(t) + ehC (t).(1.12)In (1.11) and (1.12), the terms ehL(t) and ehC (t), (where ehL(t) = [−vL(t −∆t) −2L∆t iL(t −∆t)] and ehC (t) = [vC(t −∆t) + ∆t2C iC(t −∆t)]) are called history sourcesbecause they depend on the voltage and current of the inductor and capacitor atthe previous time-step [19].The EMTP solution method is highly accurate because the discretization of thedifferential equations retains the physical relationship between voltage and current91.3. Transient Solution Methods Used in Transient Stability Studiesin the inductors and capacitors, inherently tracking network fluctuations by trac-ing the voltage and current instantaneous waveforms. In addition, the equivalentvoltage across the inductor and the current through the capacitor are dependent onthe previous value of the voltage and current at the time-step before t −∆t. Thismeans that every discrete value of the solution, e.g., v(t) and i(t), remembers thevalue of all the previous solution points, also called the variable’s history. Keepingthe variable’s history is important as it mathematically represents the physical na-ture of an inductor or capacitor, which store the kinetic energy, or memory, in theirmagnetic and electric fields.The EMTP has been used as an alternative to the traditional phasor solution meth-ods in transient stability studies [21][22]; however, due to the discretization of thecontinuous-time signal, its high level of accuracy comes with two limitations:1. There exists a direct relationship when mapping a signal between time andfrequency. A time-window of length Tc, discretized into k = N samples ata time-interval of ∆t, corresponds to a frequency-window of length Fc, alsodiscretized into k = N samples. Frequency is inversely proportional to time:F = 1T , so Fc =1N∆t . The frequency-window is centered around zero andincludes positive frequencies up to half the of the frequency window, Fc2 , andthe complex conjugates of these frequencies. Half of the frequency window,Fc2 , is the maximum frequency that can exist in the signal and it is called theNyquist frequency [23].2. It takes approximately 8-10 samples to sample one period of a sinusoidal wave-form [20].With the above conditions, if N = 10 samples and F = 1T and the trapezoidal ruleis used to discretize the continuous-time signal, to ensure that there is less than3% error in the discretization of the continuous-time signal, the maximum signalsimulated is Fmax =15fNy [20]. Therefore, the time-step used is∆t =110× Fmax . (1.13)or∆t =12× fNy . (1.14)101.3. Transient Solution Methods Used in Transient Stability StudiesBased on (1.13), for transient stability studies around 60 Hz, the EMTP will requirea time-step of approximately 1 ms to solve the electrical network. In addition, theEMTP time-step is limited by the travelling time of the waves in the transmissionlines. If the transmission lines are modelled with the constant-parameter [19] orfrequency-dependent [24] line models, the time-step may be smaller.Although the EMTP solution can capture the detailed network frequency fluctua-tions due to the discretization of the differential equations, the small discretizationtime-step sizes required makes the computer simulation time longer than if the tra-ditional phasor solution, which uses the constant 60-Hz frequency in its model ofthe electrical network solution, was used. This is a disadvantage in large systems;albeit, the phasor solution is a more approximate solution.1.3.3 Time-varying Phasors and the Introduction of the ShiftedFrequency Analysis Solution MethodIn the past decade, there has been growing research interest in developing a newsolution method that can capture the network fluctuations around the grid fre-quency, like the EMTP, whilst maintaining the computational speed of the tradi-tional phasor method (for example, in [25][26][27][28][29][30][31][32][33]). The ideaof a “time-varying phasor” germinated in the early 1990s when V. Venkatasubra-manian noticed that the phasor method gave “dubious results” for voltage stabilitystudies because it could not capture the large dynamics [25]. He began workingon a phasor solution method which would be able to track the system’s dynam-ics during these transients. This phasor was called a “fast time-varying phasor”[34][35]. Around the same time, G. Verghese [26][36] developed a generalized aver-aging method based on time-varying Fourier coefficients for state-space averaging ofpower electronic converters, which he proposed could also be used as a type of time-varying phasor dynamic approach for the electric power network [36]. G. Verghese’smethod and V. Venkatasubramanian’s method branched off into two directions forthe research into dynamic phasors: the use of time-varying Fourier coefficients andthe use of a phasor transformation.G. Verghese and Dynamic PhasorsThe generalized averaging method developed by G. Verghese [26][36] was based onFourier analysis, in which a continuous-time signal is approximated into a sum of its111.3. Transient Solution Methods Used in Transient Stability Studiesexponentials and a time-window of length[t− Tc, t], “slides” across the continuous-time signal and changes the values of the coefficients Xn and θnfo = Xo +X1ej(ω1t+θ1) +X∗1e−j(ω1t+θ1) +X2ej(ω2t+θ2) +X∗2e−j(ω2t+θ2)+ ...+Xnej(ωnt+θn) +X∗ne−j(ωnt+θn) (1.15)orx(t− Tc + s) =k=∞∑k=−∞(〈x〉k(t)ejkωs(t−Tc+s)), (1.16)and ωs =2piTc, (0, Tc] and〈x〉k(t) is the kth time-varying Fourier coefficient [26]〈x〉k(t) =1Tc∫ Tc0x(t− Tc + s)e−jkωs(t−Tc+s)ds. (1.17)The purpose of decomposing the signal into Fourier coefficients was to create anequivalent state-space model where the coefficients〈x〉k(t) are the state variables.In the power system, these variables are the voltage and current sinusoids with time-varying amplitude and phase. Equivalent complex branch equations were derivedfor the inductor and capacitor, where the time-domain quantity in the original dif-ferential equation ( ddt) is replaced by its phasor representation (jωo +ddt) [36]. Thesubsequent works of [27][28][30] further developed this method and it is commonlyreferred to as the dynamic phasor solution method. However, one concern withusing this method is that because the signal is a continuous-time signal, the entirefrequency spectrum is expressed in each new time-window. When approximatingthe signal, frequencies higher than the Nyquist frequency are then truncated, eventhough they still exist. This truncation may lead to frequency distortions if notcarefully taken care of.V. Venkatasubramanian and the Phasor TransformationV. Venkatasubramanian [34][37], on the other hand, recognized that a phasor issimply a “mathematical transformation, where the only time-varying component,frequency, is eliminated” [34]. Therefore, a time-varying phasor is also a mathemat-ical transformation, but one that does not eliminate the frequency component. V.Venkatasubramanian derived the time-varying phasor in the following manner [37].121.3. Transient Solution Methods Used in Transient Stability Studies1. He first defined time-varying phasors for voltage and current as~E(t) = E(t)ejθe(t), (1.18a)~I(t) = I(t)ejθi(t). (1.18b)2. He then created a transformation operator, P such that the instantaneousvoltage and current sinusoids e(t), i(t)e(t) =√2E(t) cos(ωot+ θe(t)), (1.19a)i(t) =√2I(t) cos(ωot+ θi(t)), (1.19b)can relate to the time-varying phasorsP(e(t)) := E(t)ejθe(t). (1.20a)P(i(t)) := I(t)ejθi(t). (1.20b)3. The network parameters were redefined as linear differential equations in thecomplex domain in terms of (1.18a), (1.18b). For example, the voltage andcurrent relationship in an inductor (1.7a) is redefined using the transformationoperator~EL(t) = P(eL(t))= P(L ddt(iL(t)))= Lddt(~IL(t)) + jωoL~IL(t),(1.21)where ~IL(t) is defined in(1.18b).4. The transformation operator that allows the time-differential relationship tobe defined in terms of phasors in the above equation (1.21) is P = ejωot.This transformation was derived from the Blondel transformation matrix usedin power systems for balanced three-phase signals [34]. Two integral prop-erties of P are that: (1) P is a linear transformation and (2) P( ddte(t)) =ddt(~E(t)) + jω ~E(t). These properties are important because it demonstratesthat using the single-frequency exponential can create a transformation froma time representation to a complex phasor representation. This was furtherexpanded in [31].131.3. Transient Solution Methods Used in Transient Stability StudiesThe phasor representations in (1.18a) and (1.18b) can be rewritten in (1.22a),(1.22b)by decomposing the exponential into a sum of its real and imaginary components.The concept of real and imaginary components was again further expanded in [31].~E(t) = E(t)∠δ(t)= E(t) cos(δ(t)) + jE(t) sin(δ(t))= Ed(t) + jEq(t). (1.22a)~I(t) = I(t)∠δ(t)= I(t) cos(δ(t)) + jI(t) sin(δ(t))= Id(t) + jIq(t). (1.22b)These time-varying phasors were then used to decompose the three-phase powersystems into its symmetrical components to model both balanced and unbalancedconditions. The hope was that the time-varying phasors would “unify transientanalysis and steady-state analysis” [37].H. W. Dommel, S. Henshel and modelling Dynamic Phasors using theEMTPThe work done by H. W. Dommel and S. Henshel in [31] further developed theconcept of dynamic time-varying phasors. The objective was to unify the electro-magnetic and electromechanical transients in the power system by applying thenumerical discretization methods used in the EMTP to the new complex-valueddifferential equations created by time-varying phasors. Instead of applying threetransformations to obtain one complex signal as in [37], the dynamic phasor methodin [31] expressed each phase as one complex signal multiplied by e−jωot, which gener-alized the dynamic phasor concept beyond three-phase networks. The general formof the dynamic phasor approach here is described next.Similar to [34], a complex signal is first constructed from a real signal. Definings(t) as the real signal, and S(ω) as a phasor with real and imaginary components,the relationship between s(t) and S(ω) is: s(t) = R(S(ω)). Treating the powersystem transients as bandpass signals around ωo and applying Fourier analysis, thebandpass signal is decomposed into two low-pass signals, sI(t) and sq(t), with twosinusoidal carrier signals surrounding ωo and −ωo (cos(ωot) and sin(ωot)).s(t) = sI(t) cos(ωot)− sQ(t) sin(ωot), (1.23)141.3. Transient Solution Methods Used in Transient Stability Studieswhere sI(t) and sQ(t) contain all the information about the dynamic transientssurrounding ωo, which now surround 0 Hz. Equation (1.24) is called the“complexenvelope of the real signal” or dynamic phasorD[S(t)] = sI(t) + jsQ(t). (1.24)Another way to represent the complex envelope is in (1.25), where H is the Hilberttransformation of the signal s(t)S(t) = (sI(t) + jsQ(t))ejωot = s(t) + jH[s(t)]. (1.25)In (1.25), the complex envelope is represented in terms of the real and imaginary partof the analytical signal. Replacing the arbitrary signal in (1.25) with the voltage andcurrent sinusoids, the differential equations for the network elements (1.7a)(1.7b) arenow complex-valuedD[V ](t) = jωoLD[I](t) + L ddtD[I](t), (1.26)where D denotes the “dynamic phasor” notation.Creating a low-pass signal surrounding 0 Hz from a band-pass signal surrounding ωois conceptually similar to defining a transformation that shifts a signal down from ωoto 0 Hz. The concept of formulating a transformation that can shift the frequencyof a signal was introduced in [32] and was named Shifted Frequency Analysis (SFA).J. R. Mart´ı and Shifted Frequency AnalysisShifted Frequency Analysis (SFA) is a solution method developed by J. R. Mart´ı [32]as a general framework to unify the concept of dynamic phasors with the discrete-time solution method of the EMTP. The SFA framework has its roots in discrete-time signal processing theory [32], where1. The Fourier decomposition of a sinusoidal signal, y(t) = A cos(ωot + θ) hastwo frequency components: ωo and −ωo.2. The signal e−jωt is the only signal with a single-frequency [38] and using itas the input to a linear time-invariant system, allows all the properties of theoriginal system to be retained in the output.151.3. Transient Solution Methods Used in Transient Stability Studies3. In Fourier analysis, the frequency-shifting property states that a time signalmultiplied by a complex exponential at ωo is equivalent to shifting the signalby the frequency ωox(t)ejωot → X (ω − ωo) (1.27)By making use of these concepts, a frequency transformation T = e−jωot is defined[32], where ωo is the signal’s fundamental frequency. The transformation T “shifts”the original signal by −ωo and the dynamics surrounding the original signal nowsurround the lower frequency-shifted signal. In the voltage and current electricalpower signals, the fundamental frequency is the grid frequency (60 Hz). Using theabove concepts 1 and 3 and applying T to the 60 Hz voltage and current sinusoids,the 60-Hz frequency is shifted down by -60 Hz to 0 Hz (DC) and all the dynamicssurrounding the original 60 Hz waveform are preserved in the shifted signal but nowsurrounding 0 Hz (DC). This is important because the SFA solution method, whichis a discrete-time method, can now use large time-steps corresponding to frequencyfluctuations around 0 Hz.The SFA transformation is analogous to a change of variables from one domain (thephysical domain) to another domain (the shifted domain)vSFA(t) = T vph(t), (1.28a)iSFA(t) = T iph(t), (1.28b)where vph(t) = V (t)ej(ωot+θv(t)) and iph(t) = I(t)ej(ωot+θi(t)).By applying T to the time-differential equations (1.7a) and (1.7b), equivalent net-work branches for the inductor and capacitor are created. These branches preservethe physical nature of the elements through the differential term ddt but combinedwith a complex term jωo. The SFA network elements are derived in Chapter 2.The SFA solution method originated in [32], was formally described in [33][39] andwas applied to synchronous and induction machine modelling in [40][41][42]. Ahybrid multi-rate simulation framework that interfaces the EMTP solution with theSFA solution is currently being developed [43]. In this thesis, a generalized SFA-based solution algorithm is programmed and used for transient stability studies.161.4. Research Objectives1.4 Research ObjectivesThe ability for the SFA solution to trace frequency fluctuations around the 60 Hz gridfrequency using large discretization time-steps led to the following research question:Can the SFA solution method be an alternative to the traditional phasorand the EMTP solution methods in transient stability studies due to itsability to trace the network dynamics using large time-steps?The research in this thesis also reformulates how the electromechanical equation(1.3) used in transient stability studies is evaluated. Traditionally, in (1.3), nettorque is converted to net power by using the approximate velocity ωo (1.3b)(1.3c),which simplifies the equations. In the proposed solution algorithm, (1.3) is solved indiscrete-time and the actual velocity ω is solved for at each time-step. This velocityis used in the conversion from net torque to net power. This observation led to thefollowing research question:Does using the actual velocity in the conversion between torque andpower in the electromechanical equation lead to a difference in calculat-ing transient stability limits, such as the critical clearing time?To answer the above questions, the following objectives are set:1. To compare the SFA solution method with the traditional phasor method andthe EMTP in transient stability analysis case studies, using the EMTP as thereference solution method.2. To determine the critical clearing times of different three-phase short-circuitfault contingencies in each of the case studies and to compare the results givenby the three solution methods.3. To determine if using the actual velocity instead of the approximate velocityin the electromechanical equation makes a difference in the solution accuracy.My work on these research objectives led to the following contributions:171.5. Research Contributions1.5 Research ContributionsThe main contributions in this thesis are as follows:1. A generalized SFA-based algorithm is written for a general power system net-work topology and the algorithm is used for transient stability analysis.2. In this work, the SFA solution method is formulated in a generalized manner,such that SFA can be used with any arbitrary source type, not just sinusoidalsources. This is different than in the previous works, where SFA was only usedwith sinusoidal sources [33][39].3. An interpolation procedure is introduced in the SFA algorithm, which allowsthe system to be solved with large time-steps when modelling the system in theshifted domain, while recapturing the physical solution with high resolutionwhen transformed back into the physical 60 Hz domain.4. The electromechanical equation is generalized to account the use of the actualvelocity at every time-step of the simulation instead of using the approximate60-Hz frequency in the equation. This is implemented in the SFA-based al-gorithm to compare the differences between using the known velocity and theapproximate velocity.5. A zero-crossing detection procedure for the fault current is implemented in theSFA algorithm for the transient stability studies. In this procedure, the faultcurrent is detected to cross zero before clearing the fault in the shifted do-main instead of detecting the zero-crossing of the fault current in the physicaldomain, which would have required shifting the solution back to the physi-cal domain at every time-step. By detecting the zero-crossing in the shifteddomain, the procedure saves computational time.1.6 Thesis OutlineThe thesis outline is as follows:Chapter 2 describes the SFA solution framework and solution method. The equiv-alent network branches are derived in the SFA domain and discretized using twonumerical discretization methods. An example of an RL circuit compares the SFA181.6. Thesis Outlinesolution method with the traditional phasor solution and the EMTP.Chapter 3 presents a numerical error analysis of the two discretization rule used todiscretize the SFA solution method in this work: the trapezoidal rule and the back-ward Euler rule, which demonstrates the accuracy of each discretization methodin the SFA domain. In addition to numerical accuracy, the numerical method’sresponse to a discontinuity, such as a step input, is modelled and analysed. Thisanalysis is useful in understanding how the solution method can handle networkdiscontinuities, such as the sudden opening of a circuit breaker. In this chapter, anexample of an RL circuit energized by a step input is modelled, which demonstratesthat the SFA solution method can be used for any arbitrary source type.Chapter 4 summarizes the main concepts of transient stability and derives the dis-cretization of the electromechanical equation using the correct network frequency.The constraints and objectives for the case studies analysed in Chapter 5, Chapter 6and Chapter 7 are also defined here.Chapter 5 applies the SFA solution method to a Single-Generator Infinite-Bus sys-tem, Chapter 6 applies the SFA solution method to a three-bus system, and Chap-ter 7 applies the SFA solution method to the IEEE 39-bus system. For each casestudy, a three-phase fault is applied to the network and the generator’s electrome-chanical variables are compared using the three solution methods, with the EMTPsolution as the reference. The SFA solution is solved using both the trapezoidal andbackward Euler methods to compare the differences between the two discretizationrules and the electromechanical equation is solved using both the approximate andexact network frequency. The critical clearing times and critical clearing angles arecalculated using each solution method and compared with each other, again, us-ing the EMTP solution as the reference. In the three-bus network of Chapter 6, athree-phase fault is applied such that the fault current has the maximum DC offset,to demonstrate that the SFA solution can capture the DC component of the faultcurrent just like the EMTP. In the three-bus network and the 39-bus network, agenerator is disconnected from the network and the solution is analysed to deter-mine how the loss of a generator has a strong effect on the network frequency.Chapter 8 summarizes the main findings of the thesis and presents ideas for futureresearch.19Chapter 2Shifted Frequency Analysis: ADiscrete-Time Solution Methodand Solution FrameworkShifted Frequency Analysis (SFA) is both a discrete-time solution method and a so-lution framework. In SFA, a continuous-time signal is transformed into an equivalentcontinuous-time complex-valued signal through a left-frequency shifting transforma-tion. The transformation T = e−jωot maps the signal from the physical domain to ashifted domain. This mapping preserves all the dynamics that surround the originalsignal.Section 2.1 introduces the SFA framework. Section 2.2 describes the SFA solutionmethod, Section 2.3 derives the equivalent SFA model for the electrical networkcomponents, Section 2.4 discretizes the equivalent SFA network elements using thetrapezoidal and backward Euler rules, Section 2.5 compares the network branchesof the SFA solution method, the EMTP solution method and the traditional phasorsolution method, and Section 2.6 illustrates the similarities and differences amongthe three solution methods in an RL circuit energized by an AC source.2.1 Shifted Frequency Analysis Solution FrameworkThe backbone of the SFA solution framework is built upon two concepts: (1) theuse of a single-exponential frequency transformation, and (2) the time-frequencyrelationship in converting a continuous-time signal to a discrete-time signal.202.1. Shifted Frequency Analysis Solution Framework2.1.1 A Single-Frequency SignalThe SFA framework defines a single-frequency transformation as T = e−jωot and ap-plies it to a signal in the physical domain. Theoretically, this is equivalent to takingthe sinusoidal time-domain signal of the form y(t) = A(t)ejωot+θ(t), which rotatesaround ωo, and left-frequency shift the signal by −ωo, such that the resulting signalrotates around 0 Hz. This is possible because the exponential signal x(t) = Ae±jωthas a single-frequency and when it is multiplied to another signal rotating at fre-quency ω1, the resulting signal will rotate at ω1 ± ω.Once the physical domain signal have been frequency-shifted to the SFA domainthe second property, the time-frequency relationship, is invoked to discretize thecontinuous-time signal into a discrete-time signal.2.1.2 The Time-Frequency RelationshipAs explained in Section 1.3.2, there is a relationship between time and frequency[23],such that any real continuous-time signal x(t) can be segmented into time-windowsof length Tc and discretized into k = N samples, where each sample has a time-width ∆t. Each time-window corresponds to a frequency-window of length Fc,discretized into k = N samples, with each sample having a frequency width ∆f . Therelationship is that Fc =1∆t and ∆t =1Fc, where ∆t is the discretization time-width,Fc is the length of the frequency-window andFc2 is the maximum frequency thatcan exist in the signal, otherwise known as the Nyquist frequency. This relationshipbears the following consequences:1. The signal in each time-window can be expressed as a sum of its exponentialcomponents, with both the positive frequencies and their complex conjugatenegative frequencies (Fourier decomposition).2. Frequencies greater than the Nyquist frequency do not exist.3. The value of ∆t is limited by (1.13).2.1.3 The SFA Solution FrameworkThe SFA framework defines the single-frequency transformation as T = e−jωot andapplies it to a signal in the physical domain. After the signal has been left-frequency212.1. Shifted Frequency Analysis Solution Frameworkshifted into the SFA world, the continuous-time signal is discretized into a discrete-time signal, adhering to the time-frequency relationship and (1.13). Because theshifted signal is of the form Y (t) = A(t)∠θ(t), which is a time-varying complexnumber (phasor), the magnitude of the phasor is the envelope of the physical signaland it captures the signal’s dynamics. The transformation retains the rotationalrelationship between real and imaginary parts of the signal so that in the shifteddomain, both real and imaginary components exist and comply with their perpen-dicular relationship (being 90◦ apart from each other).In summary, the use of the single-frequency exponential and the discretization ofthe continuous-time signal provides the following groundwork for the SFA solutionmethod:1. As the system shifts from the 60 Hz domain to the 0 Hz domain, all the dy-namics embodied in the original 60 Hz signal are preserved in the shifted 0 Hzsignal.2. The continuous-time differential equations are discretized into a discrete-timesolution method and because the signal is shifted to 0 Hz, using (1.13), thetime-steps theoretically can be infinitely large [32].3. Because the SFA solution method is a discrete-time method and the time-frequency relationship of Sec.1.3.2 holds true, the time-step is fixed to themaximum frequency wished to be accurately simulated. Thus, frequenciesoutside the fixed frequency-window do not exist. This is one difference betweenSFA and the other dynamic phasor methods, such as [26][27][30].4. Because the transformation is a frequency shifting rotational transformation,the perpendicularity between the real and imaginary components of a com-plex signal is preserved. Therefore, the voltage or current source used in thephysical domain can be represented as a complex exponential, with both realand imaginary components existing in the shifted domain without need for theHilbert Transform as used in [33] and [39].Section 2.2 formulates the SFA solution method.222.2. The SFA Solution Method2.2 The SFA Solution MethodTo model the electrical network with the SFA solution method, first the frequency-shifting transformation T = e−jωot is applied to the voltage and current sources.In the power system, the voltage and current time-domain signals v(t) and i(t)rotate sinusoidally at the grid frequency, ωo. The frequency decomposition of asinusoidal signal is two frequencies: ωo and −ωo. When the exponential transfor-mation is applied to v(t) and i(t), the frequency components become ω + ωo andω − ωo. Because the frequency of the transformation is chosen as ω = −ωo, thenthe frequency components resulting are ω = 0 and ω = −2ωo. By disregarding thenegative frequency component and multiplying the amplitude by two, the sinusoidcan be represented as a single-frequency exponential rotating around 0 Hz. This isbecause the real part of the exponential is the sinusoidal signal Rejωt = cos(ωt).v(t) = V¯ (t) cos(ωot+ θv(t))ejωtv(t) =V (t)2(ej(ωot+θv(t))ejωt)+V (t)2(e−j(ωot+θv(t))ejωt)if ω=−ωov(t) =V (t)2(ej((−ωo+ωo)t+θv(t)))+V (t)2(e−j((−ωo−ωo)t+θv(t)))i(t) = I(t) cos(ωot+ θi(t))ejωti(t) =I(t)2(ej(ωot+θi(t))ejωt)+I(t)2(e−j(ωot+θi(t))ejωt)if ω=−ωoi(t) =I(t)2(ej((−ωo+ωo)t+θi(t)))+I(t)2(e−j((−ωo−ωo)t+θi(t)))Dropping the negative frequencies and multiplying by 2V¯ (t) = V (t)ejθv(t). (2.1a) I¯(t) = I(t)ejθi(t). (2.1b)232.2. The SFA Solution MethodAlternatively, (2.1) is equivalent to representing the voltage and current sinusoidsusing their Euler representation at the frequency ωo and applying the frequency-shifting exponential with frequency ω = −ωo to create phasorsV¯ (t) = V (t)ej(ωot+θv(t))e−jωot= V (t)ejωote−jωotejθv(t)= V (t)∠θv(t).(2.2)I¯(t) = I(t)ej(ωot+θi(t))e−jωot= I(t)ejωote−jωotejθi(t)= I(t)∠θi(t).(2.3)The resulting signal shifts to the frequency ω = ω−ωo. When ω = ωo, the frequencyis 0 Hz (ω = ωo − ωo) and all dynamics that surrounded the original signal aroundωo surround 0 Hz as is depicted in Figure 2.1.ttAStart with a sinusoidal signal rotating at ocTttTake one period of the signal (from t=0 to t=Tc)AtDo a Fourier decomposition to obtain the frequency components ffof of2ADisregard the negative frequency components and multiply the amplitude by 2 of fo( ( ))( ) ( ) j t tphA t A t e A ( ) ( ) ( )SFAA t A t tfShift the signal to 0Hz using the transformation 0oj teTA(i) (ii)cT(iii) (iv) (v)Figure 2.1: The process to create a complex signal with frequency side-bands around 0 Hzfrom a sinusoidal signal with frequency side-bands around ωo, such that a sinusoidal signalrotating around 60 Hz can now rotate around 0 Hz.The same discrete-time solution techniques employed in the EMTP for the electricalnetwork are applied to the continuous-time signals in the complex SFA world. First,the differential equations governing the network elements, resistors, inductors andcapacitors, are discretized into complex network equations using the trapezoidal orbackward Euler numerical discretization methods. A discretization time-step whichcomplies with the time-frequency relationship described in 2.1.2 is chosen. Thesecomplex network branches form the network admittance matrix [Y ] and throughnodal analysis, [Y ][V ] = [I], the node voltages are solved for as phasors with atime-varying magnitude and angle, where the magnitude is the envelope that traces242.2. The SFA Solution Methodthe dynamic behaviour. Like with the EMTP, the network branches depend on thesolution at each previous time-step or the variable’s history. The history sourcesare included in the matrix [I] when solving the nodal analysis, such that [Itotal] isa matrix of external current sources [I] plus a matrix of history sources [H]. Thediscrete-time solution is then converted back to continuous-time using a linear in-terpolation scheme. This one-step linear interpolation allows the large time-stepsto be used to solve the system in the SFA domain, while bringing back the detailedsolution to the physical domain. The continuous-time complex signal in the SFAdomain is now transformed back to the physical domain through the inverse trans-formation T −1 = e+jωot. Finally, the imaginary portion of the complex-signal isdropped to obtain the real-valued physical solution. The SFA solution method isdescribed in Figure 2.2.Discretize Circuit in SFA Domain Apply reverse transformation to go back to the physical domain and take the real part. Shifting TransformationInput Signal (exponential) CT-Physical World Transform Input to SFA DomainCT-SFA WorldTransformed Input SignalDT-SFA WorldBackward Euler DiscretizatonInterpolate discrete time SFA signal back to continuous timeCT-SFA World( )v ttReverse TransformationCT-Physical System Output Signalo( ( ))( ) ( ) j tphtx t A t e oT tje() ) )( ( SFAx t A ttotjeT ( ) { ( ) }( )phy t R A tt InterpolatedSolved Figure 2.2: The SFA solution method: a real continuous-time signal is shifted to the SFAworld, discretized and solved for in the discrete-time SFA domain, interpolated into acontinuous-time signal in the SFA domain, and reverse transformed back to the physicaldomain. The imaginary portion is dropped and the real part of the signal is kept to createthe physical continuous-time solution.252.2. The SFA Solution Method2.2.1 The Choice of Time-Step in the SFA Solution MethodThe time-step used in the SFA discrete-time solution method can theoretically beinfinitely large as formulated in (1.13): ∆t = 110×0 ms as the signal now rotatesaround 0 Hz. Physically, there will be frequency deviations in the dynamics of thetransient. For example, if the frequency deviates 4 Hz from 60-64 Hz in the network,the time-step used can be ∆t = 110×4 = 25 ms.With the un-shifted 60-Hz signal, the EMTP would require a time-step of∆t = 110×64 = 1.56 ms. There is a time-step ratio between the EMTP and the SFAof 15 times. In reality, the EMTP time-step would be approximately 10 times less,or 100µs and the time-step ratio between the EMTP and the SFA solution will beroughly 160 times. Note, however, that because the SFA solution is complex-valued,the total number of operations in the solution will decrease this ratio.2.2.2 A Summary of the SFA Solution MethodA summary of the SFA solution method is provided below.1. Apply T = e−jωot to the continuous-time signals vph(t) and iph(t) in the sourcesand network elements.2. Form a discrete-time signal from the shifted continuous-time signal using anumerical discretization method, such as the trapezoidal or backward Eulerrules on the sources and network elements.3. Form the admittance matrix [Y ] and solve for node voltages using nodal anal-ysis, where the current source vector [Itotal] = [Iexternal] + [H], where [H] isthe vector of history sources and [Y ][V ] = [Iexternal + H]. The node voltagesare phasors with a time-varying magnitude and angle.4. Do a linear interpolation of the discrete-time SFA solution to convert back toa continuous-time SFA solution.5. Apply the inverse transformation T −1 = ejωot to the SFA solution and dropthe imaginary portion of the complex signal.262.3. The SFA Equivalent Network Elements2.3 The SFA Equivalent Network ElementsThe shifting transformation T is applied to the voltages and currents in the time-differential equations that define the network elements: resistors, inductors and ca-pacitors. The SFA equivalent representation of these network elements are derivedbelow. In the following derivations, the subscript “SFA” is used for the shifteddomain and the subscript “ph” is used for the physical domain. Again, vSFA = T vphor vph = T −1vSFA where T = e−jωot.ResistorsThe relationship between voltage and current in a resistor isvRph(t) = RiRph(t). (2.4)By applying T to both sides of (2.4), one obtainsT vRph(t) = T (RiRph(t))vRSFA(t) = RiRSFA(t),(2.5)or in phasor notation:V¯RSFA(t) = RI¯RSFA(t). (2.6)The equivalent relationship between voltage and current in the shifted domain isthe same as in the physical domain, which matches the physical nature of a linearresistor, which dissipates energy.InductorsThe relationship between voltage and current in an inductor isvLph(t) = Lddt(iLph(t)). (2.7)Applying the change of variables transformation T to both sides of (2.7)(vph = T −1vSFA, iph = T −1iSFA)272.3. The SFA Equivalent Network ElementsT −1vLSFA(t) = Lddt(T −1iLSFA(t))e+jωotvLSFA(t) = Lddt(e+jωotiLSFA(t))e+jωotvLSFA(t) = Lddt(iLSFA(t))e+jωot + jωoLiLSFA(t)e+jωotvLSFA(t) = Lddt(iLSFA(t))+ jωoLiLSFA(t),(2.8)or in phasor notationV¯LSFA(t) = LddtI¯LSFA(t) + jωoLI¯LSFA(t). (2.9)CapacitorsThe relationship between voltage and current in a capacitor isiCph(t) = Cddt(vCph(t)). (2.10)Applying the change of variables transformation T to both sides of (2.10)T −1icSFA(t) = Cddt(T −1vcSFA(t))e+jωoticSFA(t) = Cddt(e+jωotvcSFA(t))e+jωoticSFA(t) = Cddt(vcSFA(t))e+jωot + jωoCvcSFA(t)e+jωoticSFA(t) = Cddt(vcSFA(t))+ jωoCvCSFA(t),(2.11)or in phasor notation:I¯CSFA(t) = CddtV¯CSFA(t) + jωoCV¯CSFA(t). (2.12)The equivalent network branch of either an inductor or a capacitor has both a dif-ferential part ddt and a complex impedance jωoL or jωoC. Therefore, the voltageacross the inductor and the current through the capacitor are complex-valued orphasors.282.4. The Discretization of the SFA Network ElementsSourcesThe sinusoidal source is represented using Euler’s representation of a complex signalvph(t) = V (t)ej(ωot+θv(t)), (2.13)andiph(t) = I(t)ej(ωot+θi(t)). (2.14)Applying the transformation to both sides of (2.13) and (2.14) givesT vph(t) = T V (t)ej(ωot+θv(t))vSFA(t) = e−jωotV (t)ejωotejθv(t)V¯SFA(t) = V (t)ejθv(t).(2.15)T iph(t) = T I(t)ej(ωot+θi(t))iSFA(t) = e−jωotI(t)ejωotejθi(t)I¯SFA(t) = I(t)ejθi(t).(2.16)Equations (2.15) and (2.16) are identical to (2.2) and (2.3).2.4 The Discretization of the SFA Network ElementsA numerical discretization rule creates continuous-time differential equations thatmodel the inductors and capacitors in (2.9) and (2.12) into discrete-time differenceequations. The discretization rule chosen must be both numerically accurate incomparison with the continuous-time solution and must maintain numerical stabil-ity during discontinuities. J. R. Mart´ı and J. Lin [38] validated that in the EMTP,both the trapezoidal and the backward Euler discretization rules maintain both nu-merical accuracy and stability when discretizing a first-order differential equation.Because the EMTP solution method is applied to the signals in the SFA domain,it is hypothesized that the trapezoidal and backward Euler rules hold the same nu-merical accuracy and stability in the SFA domain. This hypothesis is verified inChapter 3 through the error analysis of the two rules in the SFA domain.Sections 2.4.1 and 2.4.2 derive the discrete-time branch equivalents in the SFA do-main using the trapezoidal and the backward Euler rules. For simplicity, the sub-script “SFA” is dropped in the derivations in this section.292.4. The Discretization of the SFA Network Elements2.4.1 Equivalent SFA Network Elements Using the TrapezoidalDiscretization RuleThe trapezoidal rule is applied to the SFA equivalent network elements. The trape-zoidal rule calculates the area under the curve of a trapezoid as shown in Figure 2.3.t t tArea( )v t( )v t( )v t ttFigure 2.3: Numerical discretization with the trapezoidal rule, where the area under thecurve is Area = v(t)+v(t−∆t)2 (∆t).ResistorsIn a resistor, there are no differentials in the relationship between voltage and cur-rent, and the continuous-time equation (2.4) remains the same in discrete-time.V¯R(t) = RI¯R(t). (2.17)InductorsThe trapezoidal rule is applied to the SFA equivalent inductor given in (2.9)t∫t−∆tV¯L (t) dt =t∫t−∆tLddtI¯L(t) +t∫t−∆tjωoLI¯L(t)dtV¯L (t) + V¯L (t−∆t)2∆t = L(I¯L(t)− I¯L(t−∆t)) + jωoL( I¯L (t) + I¯L (t−∆t)2∆t).(2.18)302.4. The Discretization of the SFA Network ElementsRearranging for voltage on the left side and current on the right sideV¯L(t) + V¯L (t−∆t) = 2L∆t(I¯L(t)− I¯L(t−∆t))+ jωoL(I¯L(t) + I¯L(t−∆t))V¯L(t) = −V¯L(t−∆t) + 2L∆t(I¯L(t)− I¯L(t−∆t))+ jωoL(I¯L(t) + I¯L(t−∆t))V¯L(t) = (2L∆t+ jωoL)I¯L(t) + [(−2L∆t+ jωoL)I¯L(t−∆t)− V¯L(t−∆t)]V¯L(t) = (2L∆t+ jωoL)I¯L(t)− [(2L∆t− jωoL)I¯L(t−∆t) + V¯L(t−∆t)].(2.19)The equivalent inductor using the trapezoidal rule is an impedance 2L∆t + jωoL inseries with a history voltage source, because it retains the voltage at the previoustime-step t −∆t [19]. The history source is of the opposite polarity as the voltageacross the inductor V¯L and is a complex vectorE¯h(t) = (2L∆t− jωoL)I¯L(t−∆t) + V¯L(t−∆t). (2.20)CapacitorsSimilarly, the trapezoidal rule is applied to the SFA equivalent capacitor given in(2.12)t∫t−∆tI¯C (t) dt =t∫t−∆tCddtV¯C(t) +t∫t−∆tjωoCV¯C(t)dtI¯C (t) + I¯C (t−∆t)2∆t = C(V¯C(t)− V¯C(t−∆t)) + jωoC( V¯C (t) + V¯C (t−∆t)2∆t).(2.21)Rearranging (2.21), one obtainsI¯C(t) + I¯C (t−∆t) = 2C∆t(V¯C(t)− V¯C(t−∆t))+ jωoC(V¯C(t) + V¯C(t−∆t))I¯C(t) = −I¯C(t−∆t) + 2C∆t(V¯C(t)− V¯C(t−∆t))+ jωoC(V¯C(t) + V¯C(t−∆t))I¯C(t) = (2C∆t+ jωoC)V¯C(t) + [(−2C∆t+ jωoC)V¯C(t−∆t)− I¯C(t−∆t)]I¯C(t) = (2C∆t+ jωoC)V¯C(t)− [(2C∆t− jωoC)V¯C(t−∆t) + I¯C(t−∆t)].(2.22)312.4. The Discretization of the SFA Network ElementsThe equivalent capacitance using the trapezoidal rule is an admittance 2C∆t + jωoCin parallel with a history current source, which is a complex vectorI¯h(t) = (2C∆t− jωoC)V¯C(t−∆t) + I¯C(t−∆t). (2.23)To parallel with the inductance, (2.22) can be rearranged in terms of the voltageacross the capacitanceV¯C(t) =12C∆t + jωoC[I¯C(t) + (2C∆t− jωoC)V¯C(t−∆t) + I¯C(t−∆t)]. (2.24)2.4.2 Equivalent SFA Network Elements Using the BackwardEuler Discretization RuleThe backward Euler rule is applied to the SFA network elements. The backwardEuler rule calculates the area of a rectangle under the curve, shown in Figure 2.4.t t tArea( )v t( )v t( )v t ttFigure 2.4: Numerical discretization with the backward Euler rule, where the area underthe curve is Area = v(t)∆t.ResistorsLike with the trapezoidal rule, the continuous-time equation (2.4) remains the samein discrete-timeV¯R(t) = RI¯R(t). (2.25)322.4. The Discretization of the SFA Network ElementsInductorsThe backward Euler rule is applied to the SFA inductor given in (2.9)t∫t−∆tV¯L (t) dt =t∫t−∆tLddtI¯L(t)dt+t∫t−∆tjωoLI¯L(t)dtV¯L(t)∆t = L(I¯L(t)− I¯L(t−∆t)) + jωoLI¯L(t)∆tV¯L(t) =L∆t(I¯L(t)− I¯L(t−∆t))+ jωoLI¯L(t)V¯L(t) = (L∆t+ jωoL)I¯L(t)− L∆t(I¯L(t−∆t)).(2.26)The equivalent inductor using the backward Euler method is an impedance L∆t+jωoLin series with a history voltage sourceE¯h(t) =L∆t(I¯L(t−∆t)). (2.27)CapacitorsSimilarly, the backward Euler rule is applied to the SFA capacitor given in (2.12)t∫t−∆tI¯C (t) dt =t∫t−∆tCddtV¯C(t)dt+t∫t−∆t(jωoCV¯C(t))dtI¯C(t)∆t = C(V¯C(t)− V¯C(t−∆t)) + jωoCV¯C(t)∆tI¯C(t) =C∆t(V¯C(t)− V¯C(t−∆t))+ jωoCV¯C(t)I¯C(t) = (C∆t+ jωoC)V¯C(t)− C∆t(V¯C(t−∆t)).(2.28)The equivalent capacitance using the backward Euler rule is an admittance C∆t +jωoC in parallel with a history current sourceI¯h(t) =C∆t(V¯C(t−∆t)). (2.29)To parallel with the inductance, (2.28) can be rearranged in terms of the voltageacross the capacitance to obtainV¯C(t) =1C∆t + jωoC[I¯C(t) +C∆tV¯C(t−∆t)]. (2.30)332.5. A Comparison of the Traditional Phasor, the EMTP, and the SFA Network ElementsComparing (2.19) and(2.22) with (2.26) and (2.28), the backward Euler discretiza-tion method is easier to implement than the trapezoidal discretization method.2.5 A Comparison of the Traditional Phasor, theEMTP, and the SFA Network ElementsIn the phasor domain, a change of operators transforms the time-domain differentialequations into a frequency-domain complex number: ddt → s = jωo.In the EMTP domain, the continuous-time differential equations are discretizedto discrete-time equations. The relationship between voltage and current in thedifferential equations are preserved in the discrete-time equivalent representation.In the SFA domain, there is a change of variables from the physical time-domaindifferential equations to the complex time-domain differential equations through thetransformation T as well as a discretization of the continuous-time equations. As inthe EMTP, the relationship between voltage and current in the differential equationsare preserved in the discrete-time equivalent representation, and as in the phasorrepresentation, the equivalent branches are complex-valued.Equation (2.31) gives the relationship between voltage and current in an inductorgiven by the traditional phasor, the EMTP and the SFA solution methods, usingthe backward Euler rule for the EMTP and SFA solutions.vL(t) = LddtiL(t)V¯L(ω) = jωoLI¯L(ω) Traditional PhasorvL(t) =L∆t(iL(t)− (iL(t−∆t)) EMTPV¯L(t) = (L∆t+ jωoL)I¯L(t)− L∆t(I¯L(t−∆t)) SFA(2.31)Tables 2.1 and 2.2 compare the network elements in the traditional phasor domain,the EMTP domain, and the SFA domain using the trapezoidal discretization method(Table 2.1) and the backward Euler discretization method (Table 2.2). Both thetraditional and SFA phasor equivalent branches are complex-valued whereas theEMTP’s equivalent branch is real-valued. It can be seen that the SFA solution is ahybrid between the traditional phasor and the EMTP solution.342.5. A Comparison of the Traditional Phasor, the EMTP, and the SFA Network ElementsTable 2.1The branch equivalents for the traditional phasor, the EMTP, and theSFA methods using the trapezoidal discretization rule for the SFA and theEMTP.ElementTraditionalPhasorEMTP SFAResistor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV 2Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t o2Lj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o12Cj Ct ( )CV t ( )hCE t ( )Ci t 2tC ( )hCe t ( )Cv t o1j C CI CV Inductor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV 2Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t o2Lj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o12Cj Ct ( )CV t ( )hCE t ( )Ci t 2tC ( )hCe t ( )Cv t o1j C CI CV Capacitor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV 2Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t o2Lj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o12Cj Ct ( )CV t ( )hCE t ( )Ci t 2tC ( )hCe t ( )Cv t o1j C CI CV Table 2.2The branch equivalents for the traditional phasor, he EMTP, and theSFA methods using the backward Euler discretization rule for he EMTPand SFA.ElementTraditionalPhasorEMTP SFAResistor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV 2Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t o2Lj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o12Cj Ct ( )CV t ( )hCE t ( )Ci t 2tC ( )hCe t ( )Cv t o1j C CI CV Inductor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t oLj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o1Cj Ct ( )CV t ( )hCE t ( )Ci t tC ( )hCe t ( )Cv t o1j C CI CV Capacitor R RI RV R ( )Ri t R ( )RI t ( )Rv t ( )RV t oj L LI LV Lt ( )hLe t ( )Lv t ( )Li t ( )Cv t ( )Ci t ( )hLE t ( )LI t oLj Lt ( )LV t ( )CI t ( )CV t ( )Cv t ( )CI t o1Cj Ct ( )CV t ( )hCE t ( )Ci t t ( )hCe t ( )Cv t o1j C CI CV 352.6. An RL Circuit Energized with an AC Source Using the Traditional Phasor, the EMTP and SFA2.6 An RL Circuit Energized with an AC SourceUsing the Traditional Phasor, the EMTP and SFAAn RL circuit energized with an AC voltage source, Figure 2.5, is simulated todemonstrate the similarities and differences among the three solution methods. InFigure 2.5, the three solution methods’ equivalent circuits are depicted.Nodal analysis is used for all three methods, [Y ][V ] = [Itotal], where• [Y] is the admittance matrix formed by the network’s branches.• [Itotal] =[I] + [H] is a vector of all the external sources injected into the networkplus the history source derived from the discretized L branch (in the EMTPand SFA solutions).• [V] is the vector of all the unknown bus voltages.In the circuit of Figure 2.5, using the trapezoidal rule for the EMTP and SFAsolutions, and where ωo is 2pi60:• [V ] =[v1(t)v2(t)], where the source value vs(t) = 10 cos(ωot+ 0◦) inFigures 2.5(a) and 2.5(c) and V¯s = 10∠0◦ in Figures 2.5(b) and 2.5(d).• [HEMTP] = 0ehL (t)2L∆t, [HSFA] = 0E¯hL (t)2L∆t+jωoL, [Hphasor]=[00].• [I] =[00]and [Itotal] = [H] + [I].The admittance matrix [Y] is compared among the methods• The admittance matrix for the traditional phasor method has complex terms:Yphasor =[1R1+ 1R2 − 1R2− 1R2 1R2 + 1jωoL]• The admittance matrix for the EMTP has no complex terms:YEMTP =[1R1+ 1R2 − 1R2− 1R2 1R2 + ∆t2L]362.6. An RL Circuit Energized with an AC Source Using the Traditional Phasor, the EMTP and SFA• The admittance matrix for SFA has complex terms:YSFA = 1R1 + 1R2 − 1R2− 1R2 1R2 + 12L∆t+jωoL 1 0.1R 1 0.1R 2 10R 2 10R ( )sv t sV ( )Li t ( )Lv t LI LV 20mHL V1 V2 V1 V2 7.54LZ j (a) RL circuit represented in time. 1 0.1R 1 0.1R 2 10R 2 10R ( )sv t sV ( )Li t ( )Lv t LI LV 20mHL V1 V2 V1 V2 7.54LZ j (b) RL circuit represented in the phasor domain. 1 0.1R 2 10R ( )sv t ( )Lv t ( )Li t 2200Lt ( )hLe t 1 0.1R 2 10R ( )sV t ( )LI t ( )LV t 2200Lt ( )hLE t V1 V2 V1 V2 7.54LZ j (c) RL circuit represented in the EMTP. 1 0.1R 2 10R ( )sv t ( )Lv t ( )Li t 2200Lt ( )hLe t 1 0.1R 2 10R ( )sV t ( )LI t ( )LV t 2200Lt ( )hLE t V1 V2 V1 V2 7.54LZ j (d) RL circuit represented in the SFA domain.Figure 2.5: Four representations of an RL circuit energized by an AC source: (a) the time-domain, (b) the phasor domain, (c) the EMTP, and (d) the SFA domain.It can be seen from Figure 2.6 that the phasor steady-state solution and the envelopein the SFA domain produce identical results and that when the SFA solution is in-terpolated back to continuous-time and transformed back into the physical domain,SFA and the EMTP produce the same results. This demonstrates the ability for theSFA solution to give both the dynamic envelope and the instantaneous-time solution.Due to the time constant in the circuit, (τ = 1.98 ms), the time-step used in theEMTP is ∆t = 110τ or ∆t = 200µs (0.2 ms). To solve for the SFA envelope, thetime-step used is ∆t = 50 ms. For interpolation back to the time-domain, theinterpolation time-step used is ∆t = 200µs (0.2 ms).372.6. An RL Circuit Energized with an AC Source Using the Traditional Phasor, the EMTP and SFA0.2 0.3 0.4−10−50510Voltage (V)Voltage at Node 10.2 0.3 0.4−505Voltage (V)Voltage at Node 20.2 0.3 0.4−0.500.5Time (s)Current (A)Current through the InductorEMTP Δt=0.2msSFA-BE (Time Plot) Δt=0.2msSFA-BE (Envelope) Δt=50msClassic Phasor Δt=50ms0.2 0.3 0.4−10−50510Voltage (V) V1EMTP Δt=0.2msSFA-BE (Time Plot) Δt=0.2msSFA-BE (Envelope) Δt=50msTraditional Phasor Δt=50msFigure 2.6: The node voltages and current through the inductor of the RL circuit in Figure2.5 comparing the traditional phasor, the EMTP, and the SFA solution methods.38Chapter 3The Numerical DiscretizationMethods Used in SFA3.1 The Evaluation of a Numerical DiscretizationMethod Using the Frequency ResponseA discrete-time signal discretized from a physical continuous-time signal is an ap-proximation of the physical signal and therefore, the discrete-time solution is an ap-proximation of the continuous-time solution. How close the approximated discrete-time solution is compared to the continuous-time solution depends on both thediscretization time-step ∆t and the discretization method used [38]. The discretiza-tion time-step is important as it determines the maximum frequency that can besimulated without distortion introduced [38] as discussed in Section 1.3.2. The dis-cretization method used is equally important as it determines how much distortionwill be produced onto the discrete-time solution [2].Different discretization methods use different approaches to approximate the continuous-time signal. For example, the trapezoidal discretization rule approximates the areaof a trapezoid under the given signal (Figure 2.3), whereas the backward Euler ruleapproximates the area of a rectangle (Figure 2.4). Due to their different approaches,the two discretization rules contain different intrinsic errors and it is important toevaluate these errors when choosing when it is best to use one discretization ruleover another. How large the error produced by its approximation to the physicalsolution is called the accuracy of the discretization rule.Accuracy, however, is not the only aspect important in evaluating the behaviourof a discretization rule. The rule’s numerical stability is equally important [19],especially when evaluating how the numerical method will respond to discontinu-393.2. The System’s Frequency Response in the Physical Domain and in the SFA Domainities in the system. Therefore, the chosen discretization rule must maintain boththe fidelity of the continuous-time signal and maintain numerical stability. In theEMTP, the trapezoidal and backward Euler discretization methods are proven to beboth numerically stable and accurate for first-order systems [19][38] and these tworules are used in this work for the discretization of the continuous-time SFA solution.The accuracy of a discretization rule can be analysed in the frequency domain byevaluating the frequency response of a discretized differentiator or an integratorand comparing that response with the frequency response of the continuous-timedifferentiator or integrator [2]. This is discussed in detail in the next section, Section3.2. Sections 3.3 and 3.4 analyse the accuracy of the trapezoidal and backward Eulerrules in the SFA domain through the detailed frequency response of both rules. InSection 3.5, the numerical stability of the two rules is derived in the SFA domainand compared with the numerical stability of the rules in the physical domain,using the EMTP as the reference. The results from this section are illustrated inthe response of an inductor to a step current source, shown in Section 3.6. Thisexample illustrates the behaviour of the rules during discontinuities. Finally, theexample of an RL circuit energized by a DC source in Section 3.7 demonstrates thatthe SFA solution is not limited to systems energized by sinusoidal sources.3.2 The System’s Frequency Response in the PhysicalDomain and in the SFA Domain3.2.1 The Frequency Response in the Physical DomainWhen the input to a linear time-invariant system is the single-frequency exponentialejωt the output, y(t) = H(ω)ejωt will retain the system’s inherent properties in itstransfer function H(ω) as shown in Figure 3.1. The system’s transfer function isunique to each system and so by evaluating the transfer function at every frequencyin a set range of frequencies, also known as the frequency response, how the systembehaves and responds to the input is determined [2]. If the continuous-time sys-tem’s frequency response matches the equivalent discrete-time system’s frequencyresponse, the discrete-time solution is 100% accurate and the further the discrete-time solution’s response is from the continuous-time one, the more error has incurreddue to the discretization rule. Therefore, the rule’s accuracy is determined by com-403.2. The System’s Frequency Response in the Physical Domain and in the SFA DomainInput Signal( ) j tx t e ( ) ( ) j ty t H e Output SignalLTI System( )H Figure 3.1: An LTI system’s response to an input signal that is the single-frequency expo-nential [2].paring the frequency response of the continuous-time solution with the frequencyresponse of the discrete-time solution, which has been discretized by that ruleHerror(ω) =HDT(ω)HCT(ω). (3.1)One way to determine the frequency response of the rule is to evaluate the ruleacting as an integrator or as a differentiator. This is possible to do by evaluatingthe response of an inductor or a capacitor:vL(t) = LddtiL(t) (3.2)andvC(t) =1C∫iC(t)dt. (3.3)When L = 1 or C = 1 and current is the input to the system, then the frequencyresponse of an inductor is the frequency response of a differentiator and the fre-quency response of a capacitance is the frequency response of an integrator [2]. Inthe discrete-time solution of, for example, (1.11)(1.12), the relationship between thevoltage and current input is the discrete-time transfer function and the frequencyresponse is now the response of the discrete-time equivalent integrator or differentia-tor. Figure 3.2 shows the frequency response of an inductor acting as a differentiatorand a capacitor acting as an integrator, where current is the input and voltage is theoutput. Here, current is replaced with ejωt and voltage is replaced with H(ω)ejωt,with H(ω) as the transfer function.The response of the differentiator and integrator are inversely related, and so itis sufficient to derive one or the other. In this thesis, the response of the inductor413.2. The System’s Frequency Response in the Physical Domain and in the SFA DomainFrequency Response of an Integrator( )Ci t ( )Cv t( )Cj ti t e ( ) ( )Cj tv t H e 1( )Hj C1()dtC ( ) ( )Lj tv t H e Frequency Response of a Differentiator( )Li t ( )Lv tL ddt( )Lj ti t e ( )H j L Figure 3.2: The frequency responses of an inductor and a capacitor whose values are L=1and C=1 are the frequency responses of a differentiator and an integrator when current isthe input and voltage is the output.acting as a differentiator is derived for both the trapezoidal and backward Euler rulesand the response of the integrator is the inverse of the differentiator. Alternatively,the response of the integrator could have been derived, with the response for thedifferentiator as its inverse. For each rule, the response of the rule acting as both adifferentiator and an integrator are provided (Figures 3.3 3.4 3.6 and 3.7).3.2.2 The Frequency Response in the SFA DomainIn the SFA domain, the discretization rule is applied to the already shifted signaland so to determine the rule’s accuracy, one takes the ratio of the SFA discrete-time solution’s frequency response to the SFA continuous-time solution’s frequencyresponseHerror(ω) =HDT−SFA(ω)HCT−SFA(ω)(3.4)3.2.3 The Continuous-Time Frequency Response in the SFADomainThe SFA continuous-time (SFA-CT) frequency response is the frequency response ofthe continuous-time transformed inductor in the SFA domain (V¯L(t) = Lddt I¯L(t) +jωoLI¯L(t) (2.9)). This is derived below in (3.5), where I¯L(t) is replaced with ejωt,V¯L(t) is replaced with H(ω)ejωt, and the derivative operator is replaced with thefrequency operator, ddt ⇒ s = jω.423.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA DomainV¯L(t) = LddtI¯L(t) + jωoLI¯L(t)Hct(ω)ejωt = jωLejωt + jωoLejωtHct(ω) =jωLejωt + jωoLejωtejωtHct(ω) = jL(ω + ωo).(3.5)By comparing the discrete-time frequency response in the SFA domain with thecontinuous-time frequency response in the SFA domain, the behaviour of the dis-cretization rule applied in the SFA domain can be understood. Sections 3.3 and 3.4derive the frequency response for the trapezoidal and backward Euler rules.3.3 The Frequency Response of the TrapezoidalDiscretization Rule in the SFA DomainThere are two approaches to obtain the frequency response of a discrete-time signal[2]. The first approach is the one used to obtain the SFA-CT frequency response,in which the system’s input was replaced with the single frequency exponentialx(t) = ejωt , and the system’s output was replaced with y(t) = H(ω)ejωt, whereH(ω) evaluated at every ω is the frequency response. The second approach is toapply the z-operator to the discretized equation, where z → t and z−1 → t − ∆t.When z = ejω∆t, H(z) = H(ω), which is the frequency response [2]. Both approachesare derived below and produce identical results, given in the subsequent equations(3.9) and (3.11).3.3.1 Approach One: Applying the Single-Frequency ExponentialThe equivalent inductor discretized using the trapezoidal rule was derived in (2.19)and is reiterated hereV¯L(t) = (2L∆t+ jωoL)I¯L(t)− [(2L∆t− jωoL)I¯L(t−∆t) + V¯L(t−∆t)].Replacing I¯L(t) with ejωt and V¯L(t) with H(ω)ejωtH(ω)ejωt = (2L∆t+ jωoL)ejωt + (−2L∆t+ jωoL)ejω(t−∆t) −H(ω)ejω(t−∆t). (3.6)433.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA DomainAfter cancelling out ejωt from both sides of the equation and rearranging terms,(3.7) gives the response H(ω) in terms of exponentialsH(ω)ejωt +H(ω)ejω(t−∆t) = (2L∆t+ jωoL)ejωt + (−2L∆t+ jωoL)ejω(t−∆t)H(ω)ejωt +H(ω)ejωte−jω∆t =2L∆tejωt + jωoLejωt − 2L∆tejωte−jω∆t + jωoLejωte−jω∆tH(ω)(1 + e−jω∆t) =2L∆t(1− e−jω∆t) + jωoL(1 + e−jω∆t)H(ω) =(2L∆t + jωoL) + (−2L∆t)e−jω∆t + (jωoL)e−jω∆t(1 + e−jω∆t)H(ω) =(2L∆t + jωoL) + (−2L∆t + jωoL)e−jω∆t(1 + e−jω∆t).(3.7)To simplify the response, the exponentials are rewritten in terms of cosines and sinesand e−jω∆t is rewritten as e−jω∆t2 e−jω∆t2H(ω) =e−jω∆t2e−jω∆t2× (2L∆t + jωoL)ejω∆t2 + (−2L∆t + jωoL)e−jω∆t2ejω∆t2 + e−jω∆t2=(2L∆t + jωoL)ejω∆t2 + (−2L∆t + jωoL)e−jω∆t22 cos(ω∆t2 )=2L∆t(ejω∆t2 − e−jω∆t2 ) + jωoL(ejω∆t2 + e−jω∆t2 )2 cos(ω∆t2 )=2L∆t(j2 sin(ω∆t2 )) + jωoL(2 cos(ω∆t2 ))2 cos(ω∆t2 )=2L∆t(j sin(ω∆t2 )) + jωoL(cos(ω∆t2 ))cos(ω∆t2 ).(3.8)The frequency response of the trapezoidal rule in the SFA domain is given in (3.9),where ω is any frequency in a specified range, L = 1, and ∆t is the chosen time-step.Happrox(ω) = j2L∆ttan(ω∆t2) + jωoL (3.9)3.3.2 Approach Two: Applying the Z-Domain Transfer FunctionTaking the discretized equation from (2.19), by rewriting the discretized equationin terms of z, where t → z and t − ∆t → z−1, the z-domain transfer functionis obtained. When z = ejω∆t, one gets the frequency response of the discretized443.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA Domaininductor, which is equivalent to (3.9). This is derived as follows. For convenience,(2.19) is reiterated asV¯L(t) = (2L∆t+ jωoL)I¯L(t)− [(2L∆t− jωoL)I¯L(t−∆t) + V¯L(t−∆t)].Replacing I¯L(t) with I¯(z) and I¯L(t − ∆t) with I¯(z)z−1 (and similarly V¯L(t) withV¯ (z) and V¯L(t−∆t) with V¯ (z)z−1)V¯ (z) = (2L∆t+ jωoL)I¯(z) + (−2L∆t+ jωoL)I¯(z)z−1 − V¯ (z)z−1V¯ (z) + V¯ (z)z−1 =2L∆t(I¯(z)− I¯(z)z−1) + jωoL(I¯(z) + I¯(z)(z−1))V¯ (z)(1 + z−1) =2L∆tI¯(z)(1− z−1) + jωoL(I¯(z)(1 + z−1))V¯ (z)I¯(z)=2L∆t I¯(z)(1− z−1) + jωoLI¯(z)(1 + z−1)(1 + z−1).(3.10)Defining V¯ (z)I¯(z)= H(z),H(z) =2L∆t(1− z−1)(1 + z−1)+ jωoLH(z) =2L∆t(z − 1)(z + 1)+ jωoL.Replacing z with ejω∆t, and replacing the exponentials with sinusoids as in (3.8),and evaluating H(ejω∆t) at every frequency, results inH(ejω∆t) =2L∆t(ejω∆t − 1)(ejω∆t + 1)+ jωoL=2L∆tejω∆t2 (ejω∆t2 − e−jω∆t2 )ejω∆t2 (ejω∆t2 + e−jω∆t2 )+ jωoL=2L∆t(ejω∆t2 − e−jω∆t2 )(ejω∆t2 + e−jω∆t2 )+ jωoL= j2L∆ttan(ω∆t2) + jωoLHapprox(ω) = j2L∆ttan(ω∆t2) + jωoL.(3.11)Equations (3.11) and (3.9) are identical.453.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA Domain3.3.3 The Trapezoidal Rule’s Accuracy in the SFA DomainThe trapezoidal rule’s accuracy when behaving as a differentiator is the ratio of itsdiscrete-time frequency response in (3.11) to the continuous-time frequency responsein (3.5) as shown in (3.12), where fo is the fundamental frequency, ∆t is the chosentime-step, and f is any frequency in a determined frequency range in the SFA domain(for example, from -60-Hz→60-Hz):Errordiff. =∣∣∣∣Happrox(ω)Hexact(ω)∣∣∣∣ = j 2L∆t tan(ω∆t2 ) + jωoLjL(ω + ωo)=2∆t tan(pif∆t) + 2pifo2pi(f + fo)=1pi∆t tan(pif∆t) + fo(f + fo).(3.12)The accuracy of the rule behaving as an integrator (for example, if voltage is theinput, current is the output and the system is an inductor valued at 1 H) is theinverse of (3.12)Errorint. =11pi∆ttan(pif∆t)+fo(f+fo)=1tan(pif∆t)pi∆t(f+fo)+ fo(f+fo)=1tan(pif∆t)pi∆t((f−fo)+fo) +fo(f−fo)+fo=1tan(pif∆t)pif∆t +fof.(3.13)The Accuracy of the Trapezoidal Rule in the Physical Domain,derived using the EMTP, compared to in the SFA DomainTable 3.1 compares the EMTP and the SFA frequency response accuracy ratio’sboth in per-unit (p.u.) of ∆t, where fp.u. = f∆t. In SFA, it is the termfof whichdetermines the Nyquist frequency, as this term shows that the SFA response dependson the chosen time-step (fop.u. = fo∆t). This differs from in the EMTP, in which thefrequency response is independent of the chosen time-step and the Nyquist frequencyis always 0.5 p.u. when using the trapezoidal rule [2].463.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA DomainTable 3.1A comparison of the EMTP and SFA frequency response accuracy ratio’swhen the trapezoidal rule is behaving as an integratorEMTP SFAErrorint. =1tan(pifp.u.)pifp.u.Errorint. =1tan(pifp.u.)pifp.u.+fop.u.fp.u.3.3.4 The Numerical Accuracy of the Trapezoidal Rule:Magnitude and Phase ResponsesIn the SFA domain, the rotational transformation shifts the physical frequencies fby -60 Hz. To show the response in the physical domain in Figures (3.3) and (3.4),the frequencies f are re-shifted by +60 Hz, and the 0 Hz in SFA is re-shifted back to60 Hz. Because the Nyquist frequency will be different for each time-step (fmax =1∆tand fNy =12fmax) it is easier to examine the response in physical units (Hz). Figure(3.3) presents the magnitude and phase response of the trapezoidal rule acting as adifferentiator, with chosen ∆t’s ranging form 0.1 ms → 32 ms. Figure (3.4) presentsthe magnitude and phase response of the trapezoidal rule acting as an integrator.In the trapezoidal rule, there is a discontinuity when the frequency approachesf = 60±mfNy for m = (2m− 1 : mZ). For example, at ∆t = 10 ms, the Nyquistfrequency is 50 Hz and there are discontinuities at 10 Hz (60-50) and 110 Hz (60+50).Because frequencies greater than the Nyquist frequency do not exist in the time-window, only positive frequencies lower than or equal to the Nyquist frequency haverelevant discontinuities as provided in Table 3.2. When ∆t <8.33 ms, fNy > 60 Hzand the discontinuities are less than 0 Hz.Table 3.2The discontinuities at the Nyquist frequency of different time-steps whenusing the trapezoidal rule.∆t(ms) Nyquist Frequency(Hz) Discontinuities at frequencies lower than the Nyquist frequency (Hz)8.33 60 010 50 1016.67 30 3032 15.625 13.1Regardless of the discontinuities, the trapezoidal rule is quite accurate in the SFAdomain, as the error produced hovers around∣∣∣HSFA−DT(ω)HSFA−CT(ω) ∣∣∣ = 1 as depicted in Figure473.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA Domain(3.5) for all the time-steps and except for at the discontinuities, when the phasejumps from 0 to pi, there is no distortion in the phase response.0 10 20 30 40 50 60 70 80 90 100 110 120020406080100120140160180Frequency (Hz)|H(SFA DT)/H(SFA CT)| - MagnitudeAccuracy of the Trapezoidal Rule as a Differentiator(Magnitude Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.33ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(a),(b),(c),Ideal(d)(e)(f)(i)(g)(h)10 20 30 40 50 60 70 80 90 100 110 12000.511.522.53Frequency (Hz)H(SFADT)/H(SFA CT) - PhaseAccuracy of the Trapezoidal Rule as a Differentiator(Phase Response) (a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.3ms(g) SFA Δt=10ms(h) SFA Δt=16.7ms(i) SFA Δt=32msIdeal(e)(a), (b), Ideal(c)(d) (f) (g) (h) (i)Figure 3.3: The accuracy of the trapezoidal rule acting as a differentiator (magnitude andphase response).483.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA Domain0 10 20 30 40 50 60 70 80 90 100 110 12000.511.522.53Frequency (Hz)|H(SFA-DT)/H(SFA-CT)| - MagnitudeAccuracy of the Trapezoidal Rule as an Integrator(Magnitude Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.33ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(i)(d)(a),(b),(c),Ideal(e) (f)(g)(h)10 20 30 40 50 60 70 80 90 100 110 120−0.500.511.522.533.5Frequency (Hz)H(SFADT)/H(SFA CT) - PhaseAccuracy of the Trapezoidal Rule as an Integrator(Phase Response) (a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.33ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(a), (b), Ideal(c) (d) (e) (f) (g) (h)Figure 3.4: The accuracy of the trapezoidal rule acting as an integrator (magnitude andphase response).493.3. The Frequency Response of the Trapezoidal Discretization Rule in the SFA Domain10 20 30 40 50 60 70 80 90 100 110 1200.60.70.80.911.11.21.31.41.5Frequency (Hz)|H(SFA-DT)/H(SFA-CT)| - MagnitudeFrequency Response (Magnitude) of SFA using the Trapezoidal Rule (As a Differentiator) around 0% error(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.33ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal (a),(b),(c),Ideal(f)(i)(g)(h)(d) (e)Figure 3.5: The accuracy of the trapezoidal rule acting as a differentiator (around 0% error,which occurs when∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1).3.3.5 The Distortion of the Inductor Using the Trapezoidal RuleBecause the discretization rule creates an equivalent inductor that is a discretizedversion of the physical inductor (Tables 2.1 and 2.2), the equivalent inductance isa distorted version of the physical one. (Similarly, the capacitor is also distortedby the discretization rule). One way to define how much distortion is created bythe discretization rule is to equate the discrete-time frequency response Happrox(ω)(3.11) with the continuous-time frequency response Hct(ω) (3.5). By defining thecontinuous-time inductor in (3.5) as Leq and the discrete-time inductance as theinductance L and equating Happrox(ω) to Z(ω), as the response in (3.11) is theresponse of the inductor acting as a differentiator,Z(ω) = j2L∆ttan(ω∆t2) + jωoLj(ω + ωo)Leq = j2L∆ttan(ω∆t2) + jωoLLeq = L( 2∆t tan(ω∆t2 )(ω + ωo)+ωo(ω + ωo))Leq = L( 1pi∆t tan(pif∆t)(f + fo)+fo(f + fo)).(3.14)503.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA DomainThe distortion of the inductor evaluated in (3.14) also defines the boundary condi-tions of the frequency response, which correspond to the results obtained in Figures(3.3) and (3.4).Boundary Conditions1. When f = 0, L→ Lct2. When f = −fo, L→∞3. When f → fNy, L→∞. (fNy = 12∆t)Table 3.3 compares the equivalent inductance in the shifted domain, obtained bythe SFA solution, with the equivalent inductance in the EMTP solution.Table 3.3A comparison of the EMTP’s and SFA’s distortion of the inductor whenusing the trapezoidal rule.EMTP SFALeq = Ltan(pif∆t)pif∆tLeq = L(1pi∆ttan(pif∆t)(f + fo)+fo(f + fo))3.4 The Frequency Response of the Backward EulerDiscretization Rule in the SFA DomainThe frequency response analysis performed with the trapezoidal rule is applied tothe backward Euler rule. Equivalent results are obtained when using either thesingle-frequency exponential approach or the z-domain transfer function approach.For simplicity, only the derivation using the z-domain approach is provided.3.4.1 Applying the Z-Domain Transfer FunctionThe equivalent inductor discretized using the backward Euler rule was derived in(2.26) and is reiterated hereV¯L(t) = (L∆t+ jωoL)I¯L(t)− L∆t(I¯L(t−∆t)).513.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA DomainReplacing I¯L(t) with I¯(z) and I¯L(t−∆t) with I¯(z)z−1 and V¯L(t) with V¯ (z),V¯ (z) = (L∆t+ jωoL)I¯(z)− L∆tI¯(z)z−1 (3.15)Multiplying both sides by I(z) and defining V¯ (z)I¯(z)= H(z),V¯ (z)I¯(z)= (L∆t+ jωoL)− L∆tz−1H(z) = (L∆t+ jωoL)− L∆tz−1H(z) =L∆t(1− z−1) + jωoL(3.16)In terms of frequency, one can replace z−1 with e−jω∆t and set H(z) to be theresponse for all frequencies ω,H(ω) =L∆t(1− e−jω∆t) + jωoL=L∆t(e−jω∆t2 )(ejω∆t2 − e−jω∆t2 ) + jωoL=L∆t(e−jω∆t2 )(j2 sin(ω∆t2)) + jωoL,(3.17)with the final frequency response asHapprox(ω) =2L∆t(e−jω∆t2 )(j sin(ω∆t2)) + jωoL. (3.18)3.4.2 The Backward Euler Rule’s Accuracy in the SFA DomainThe accuracy of the backward Euler rule acting as a differentiator given in (3.19)below is the ratio of the discrete-time frequency response in (3.18) to the continuous-time frequency response in (3.5), where in (3.19), fo is the fundamental frequency, ∆tis the chosen time-step for the response to be analysed at, and f is any frequency ina determined frequency range in the SFA domain (for example from -60 Hz→60 Hz)523.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA DomainErrordiff. =∣∣∣∣Happrox(ω)Hexact(ω)∣∣∣∣ = j 2L∆t(e−jω∆t2 )(sin(ω∆t2 )) + jωoLjL(ω + ωo)=2∆t(e−jω∆t2 )(sin(ω∆t2 )) + ωoω + ωo=1pi∆t(e−jpif∆t)(sin(pif∆t)) + fof + fo.(3.19)The accuracy of the rule behaving as an integrator is the inverse of (3.19)Errorint. =11pi∆t(e−jpif∆t)(sin(pif∆t))+fof+fo. (3.20)The Accuracy of the Backward Euler Rule in the Physical Domain,derived using the EMTP, compared to in the SFA DomainAs with the trapezoidal discretization, the SFA’s frequency response accuracy ratiodepends on the time-step. Table 3.4 compares the EMTP and the SFA frequencyresponse accuracy ratio’s in per-unit (p.u.) of ∆t.Table 3.4A comparison of the EMTP and SFA frequency response accuracy ratio forthe backward Euler rule behaving as an integrator.EMTP SFAErrorint. =∆t2L+∆t2L1j tan(pifp.u.)Errorint. =11pi∆t(e−jpifp.u.)(sin(pifp.u.))+fop.u.fp.u. + fop.u.3.4.3 The Numerical Accuracy of the Backward Euler Rule:Magnitude and Phase ResponsesThe numerical accuracy of the backward Euler rule is given in physical units (Hz)with the frequency response re-shifted to the physical domain through frequencyshifting f by +60 Hz. Figures (3.6) and (3.7) show the magnitude and phase re-sponse of the backward Euler rule acting as a differentiator and an integrator, re-spectively. Again, this shift in frequencies back to the physical domain does notchange the rule’s response at these frequencies; as it is just a frequency-shift. Un-533.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA Domainlike the trapezoidal rule, there is a phase distortion in the accuracy response of thebackward Euler rule, which comes from the fictitious damping resistance introducedin the distortion of the equivalent inductance (derived in Section 3.4.4). However,the error is constrained between 0◦ → 90◦ as the system is minimum-phase. Theaccuracy of the magnitude response comparing of the backward Euler rule around∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1 is shown in Figure (3.8).3.4.4 The Distortion of the Inductor Using the Backward EulerRuleSimilar to the trapezoidal rule, the distortion created by the backward Euler rule onthe inductor is derived by equating the discrete-time frequency response Happrox(ω)(3.18) with the continuous-time frequency response Hct(ω) (3.5). Rewriting (3.18)in terms of cosines and sines, one obtainsZ(ω) = j2L∆t(cos(ω∆t2)− j sin(ω∆t2)) sin(ω∆t2) + j2pifoL= j2L∆t(cos(ω∆t2) sin(ω∆t2)− j2 2L∆tsin(ω∆t2) sin(ω∆t2) + j2pifoL=2L∆tsin(ω∆t2)2 + j2L∆t(cos(ω∆t2) sin(ω∆t2).(3.21)Equation (3.21) has both a real and an imaginary part. The real part correspondsto a fictitious resistance, as seen in the EMTP [38][2]. This resistance allows energyto dissipate, which corresponds to a damping of the oscillations that are prominentwith the trapezoidal rule. In the backward Euler rule, the equivalent resistancedetermines how much the oscillations will be dampened by and the equivalent in-ductance determines how much distortion has incurred. Deriving both the real andthe imaginary parts are as follows. For the real and imaginary parts of the distortedinductor, the real part, Req is (continued on page 57):543.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA Domain0 10 20 30 40 50 60 70 80 90 100 110 1200.511.522.533.5Frequency (Hz)|H(SFA DT)/H(SFA CT)| - MagnitudeAccuracy of the Backward Euler Rule as a Differentiator(Magnitude Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.3ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(a)(b) (c) (d)(e)(f)(g) (h)Ideal(i)10 20 30 40 50 60 70 80 90 100 110 120−90−80−70−60−50−40−30−20−100Frequency (Hz)H(SFA-DT)/H(SFA-CT) - PhaseAccuracy of the Backward Euler Rule as a Differentiator(Phase Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.3ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal (a)(b)(c) (d)(e)(f)(g)(h)IdealFigure 3.6: The accuracy of the backward Euler rule acting as a differentiator (magnitudeand phase response).553.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA Domain0 10 20 30 40 50 60 70 80 90 100 110 1200.20.40.60.811.21.41.61.822.22.4Frequency (Hz)|H(SFA DT)/H(SFA CT)|- MagnitudeAccuracy of the Backward Euler Rule as an Integrator(Magnitude Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.3ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(a) (b) (c) (d) (e)(f)(g)(h)(i)(Ideal)0 10 20 30 40 50 60 70 80 90 100 110 12001020304050607080Frequency (Hz)H(SFA-DT)/H(SFA-CT) - PhaseAccuracy of the Backward Euler Rule as an Integrator(Phase Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.33ms(g) SFA Δt=10ms (h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal (a)(b)(c)(d)(e)(f)(g)(h)(i)IdealFigure 3.7: The accuracy of the backward Euler rule acting as an integrator (magnitudeand phase response).563.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA Domain0 10 20 30 40 50 60 70 80 90 100 110 1200.50.60.70.80.911.11.21.31.41.5Frequency (Hz)|H(SFA DT)/H(SFA CT)| - MagnitudeAccuracy of the Backward Euler Rule as a differentiator around 0% error(Magnitude Response)(a) SFA Δt=0.1ms(b) SFA Δt=0.5ms(c) SFA Δt=1ms(d) SFA Δt=2ms(e) SFA Δt=5ms(f) SFA Δt=8.3ms(g) SFA Δt=10ms(h) SFA Δt=16.67ms(i) SFA Δt=32msIdeal(a) (b) (c)(d)(Ideal)(e)(f) (g)(h) (i)Figure 3.8: The accuracy of the backward Euler rule acting as a differentiator (around 0%error, which occurs when∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1).For the real and imaginary parts of the distorted inductor, the real part, Req isReq =2L∆tsin2(ω∆t2). (3.22)Applying the trigonometric law, sin2(θ) = 12(1− cos(2θ)), (3.22) becomesReq =2L∆t(12(1− cos(2(ω∆t2))))=L∆t(1− cos (2(ω∆t2)))=L∆t(1− cos(ω∆t)).(3.23)Finally,Req =L∆t(1− cos(2pif∆t)). (3.24)573.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA DomainThe imaginary portion is derived as followsj(ω − ωo)Leq = j[2L∆tcos(ω∆t2) sin(ω∆t2) + ωoL]Leq =2Lω∆t cos(ω∆t2 ) sin(ω∆t2 ) + ωoL(ω + ωo)= L2∆t cos(ω∆t2 ) sin(ω∆t2 ) +ωoω(ω + ωo)= L2pi∆t cos(pif∆t) sin(pif∆t) + 2pifo2pi(f + fo).(3.25)Finally,Leq = L1pi∆t cos(pif∆t) sin(pif∆t)f + fo+fof + fo. (3.26)Boundary ConditionsThe boundary conditions in SFA for the backward Euler rule are determined bythe behaviour of the equivalent resistance and inductance. The damping resistancein SFA (3.22) is a function of the frequency (cos(ω∆t)). This contrasts with theEMTP, where the damping resistance is independent of the frequency. The followingdeductions are made with regards to the equivalent resistance:1. When f = fNy(or any multiple of fNy) → R = L∆t .2. When f = fmax → Req = 2L∆t , which is identical to the EMTP.3. When f = 0→ R = 0.For the equivalent inductance (3.25) the following boundary conditions are1. When f → −f0, L→∞.2. When f = 0, or L→ Lct.3. When f → fNy, L = fof+fo .The discontinuities observed using the trapezoidal rule do not occur with backwardEuler, as observed in (3.21), where at the Nyquist frequency (fNy), there exists avalue Zeq → L∆t + fof+fo , which makes the response continuous for backward Euler.Table 3.5 compares the behaviour of the equivalent inductance distorted by thebackward Euler rule in the EMTP and the SFA.583.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA DomainTable 3.5A comparison of the EMTP’s and SFA’s distortion of the inductor whenusing the backward Euler rule.EMTP SFAReq =2L∆tReq =L∆t(1− cos(2pif∆t))Leq = L(tan(ω∆t2 )ω∆t2) Leq = L1pi∆t cos(pif∆t) sin(pif∆t)f + fo+fof + fo3.4.5 A Comparison of The Trapezoidal and Backward Euler RuleAround 0 Hz (DC)Around 0 Hz (DC), the trapezoidal rule is close to 100% (∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1) accuracywhen using a time-step smaller than 5 ms, whereas the backward Euler rule has higherror around 0 Hz, as observed in Figure 3.9. This lends to one possibility in whythe trapezoidal rule is more accurate in detecting the DC component of the faultcurrent, as observed in the case study of Chapter 6. Recalling that the discretizationresponse in the shifted domain retains its discretization properties when re-shiftedback to the physical domain, the response at -60 Hz is the physical 0 Hz (Figure 3.9).0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2−0.100.10.20.30.40.50.60.70.80.911.11.21.31.41.5Frequency (Hz)|H(SFA DT)/H(SFA CT)| - Magnitude Accuracy of the Trapezoidal and Backward Euler Rules Around DC (As an Integrator)(a) SFA-Trap Δt=0.1ms(b) SFA-Trap Δt=0.5ms(c) SFA-Trap Δt=1ms(d) SFA-Trap Δt=2ms(e) SFA-BE Δt=0.1ms(f) SFA-BE Δt=0.5ms(g) SFA-BE Δt=1ms(h) SFA-BE Δt=2msIdeal(a) (b) (c)(d)(f) (g) (h)(Ideal)(e) Trapezoidal at Δt=0.1ms and Δt=0.5ms is close to ideal. Figure 3.9: A comparison of the accuracy of the trapezoidal and backward Euler rules actingas an integrator around 0 Hz.593.4. The Frequency Response of the Backward Euler Discretization Rule in the SFA Domain3.4.6 A Summary of the Accuracy of the Two DiscretizationRulesAlthough the backward Euler rule does not have discontinuities close to the Nyquistfrequency, in general, it is less accurate than the trapezoidal rule. For example, whenusing a ∆t = 10 ms, the trapezoidal rule has 100% accuracy(∣∣∣Happrox(ω)Hexact(ω) ∣∣∣ = 1) from58 Hz to 62.1 Hz, which is a range of 4.1 Hz, whereas the backward Euler rule has100% accuracy from 58.4 Hz to 61.6 Hz, which is a range of 3.2 Hz, approximately1 Hz less than with the trapezoidal.Both trapezoidal and backward Euler are suitable for the transient stability casestudies studied in this thesis, where a time-step of ∆t = 8 ms is used for the SFA so-lution. Figure (3.10) shows the error produced by the rules given a 20 Hz frequencydeviation around the 60-Hz frequency using ∆t = 8.3ms and ∆t = 10ms.Looking at Figure (3.10), one can see that trapezoidal introduces 0.55% less errorthan backward Euler using ∆t = 8.33 ms and 0.78% less error using ∆t = 10 msat the 50 Hz frequency, which is a 10 Hz deviation from 60 Hz. At the 55 Hz fre-quency, which is a 5 Hz deviation, using ∆t = 10 ms, trapezoidal has only an errorof 0.05%, while backward Euler has an error of 0.16% in the magnitude responseand a 0.8◦ error in the phase plot. Regardless, using backward Euler at ∆t = 8 ms,as in the case studies, only introduces an error of 0.10%. In addition, due to theextra damping term in the backward Euler rule, backward Euler is a suitable choicefor network discontinuities, such as switches opening to clear short-circuit faults,whereas trapezoidal introduces numerical oscillations.To summarize, the trapezoidal rule reaches discontinuities when approaching theNyquist frequency; however, it is more accurate than the backward Euler rule forfrequencies about ±0.5-1 Hz from 60 Hz and does not introduce a phase distortion.On the other hand, the backward Euler rule does not have discontinuities at theNyquist frequency, but does have a phase distortion and is, on average, less accuratethan trapezoidal. At frequencies around 0 Hz (DC), trapezoidal is highly accuratewhen using small time-steps, whereas backward Euler introduces large errors.603.5. The Discretization Rule’s Stability in the SFA Domain50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 700.990.99511.0051.011.0151.02Frequency (Hz)|H(SFA DT)/H(SFA CT)| - MagnitudeAccuracy of the Trapezoidal vs Backward Euler Rule around 60 Hzusing a Δt=8.33ms and Δt=10ms (Magnitude Response)(a)BE Δt=8.33ms(b)BE Δt=10ms(c)Trap Δt=8.33ms(d)Trap Δt=10msIdeal(a)1.02% error(b)1.47% error(Ideal) (c)-0.47% error(d)-0.69% error(c)-0.053% error(d)-0.08% error(a)0.101% error(b)0.158% error(50Hz)(55Hz)(50Hz)(55Hz)0% error for all rules at 60HzFigure 3.10: A comparison of the magnitude of the accuracy of the trapezoidal and backwardEuler frequency responses acting as a differentiator around 60 Hz using ∆t = 8 ms and∆t = 10 ms.3.5 The Discretization Rule’s Stability in the SFADomainThe second important property of a discretization rule is the rule’s stability, whichcan also be determined by observing the discrete-time transfer function [2]. Recallingthat the solution to a first-order differential equation has both a transient part anda steady-state part [38] [2], if the transient part dies down to 0 as t → ∞, thenthe system is stable. The transient response of the discrete-time solution is of theform ih[k] = C1pk1 + C2pk2 + ... + Cnpkn for k=t∆t , which as a transfer function isH(z) = H (z+z1)(z+z2)...(z+zn)(z+p1)(z+p2)...(z+pn) , where z and p are the zeros and poles, respectively[2]. If the transfer function’s poles are on or inside the unit circle |p| ≤ 1, thenthe rule is stable, as in the case with the trapezoidal and backward Euler rules. Inaddition to stability, to ensure causality, the transfer function’s zeros must also beon or inside the unit circle, |z| ≤ 1 [2]. In the EMTP, both the trapezoidal andbackward Euler discretization rules are derived as numerically stable [19] [38]. Thestability of the rules in the SFA domain are derived in this section.613.6. The Response to a Step Input in the SFA DomainBy discretizing the inductor using the trapezoidal rule and applying the z-transformon the discretized branch, as derived in (3.10), where H(z)=Z(ω),Z(ω) =2L∆t(z − 1z + 1)+ jωoL. (3.27)By discretizing the inductor using the backward Euler rule and applying the z-transform on the discretized branch, as derived in (3.16), where H(z)=Z(ω)Z(ω) =L∆t[1− e−jω∆t] + jωoL=L∆t[1− z−1] + jωoL=L∆t[1− 1z] + jωoL=L∆t[z − 1z] + jωoL.(3.28)Equation (3.27) is the transfer function for the inductor discretized with trapezoidal,which shows that when the trapezoidal rule acts as a differentiator, there is a poleat -1 and a zero at 1 and (3.28) shows that when the backward Euler rule acts as adifferentiator, there is a pole at 0 and a zero at 1. Therefore, the two discretizationrules in the SFA domain are both stable and causality is maintained. By derivingthese equations, certain properties of the rule, such as its behaviour during disconti-nuities can be explained. The following two sections demonstrate these behaviours.3.6 The Response to a Step Input in the SFA DomainThe solution method’s response to a discontinuity can reveal an idea about thesystem’s behaviour to a discontinuity. This understanding is critical because dis-continuities occur during switching action, which is prominent in power electronicconverters. In the modern grid, where renewable energy sources are coupled to thegrid through power electronic converters, if SFA was used in modelling the powerelectronics in future work, understanding the behaviour of the discretization rules isimportant. As described in [2], discontinuities occur whenever nature switches fromone state to another instantaneously, such as the behaviour of the current throughan inductor. One approach to model discontinuities in a system, therefore, is toenergize an inductor of 1 mH with a DC current source of 1 A [2] as in Figure 3.11.623.6. The Response to a Step Input in the SFA DomainThe relationship between voltage and current through an inductor in the SFA do-main using the trapezoidal and backward Euler rules are given in (2.19) and (2.26),respectively. Tables 3.6 and 3.7 give four time-steps of the voltage across the in-ductor when the inductor is energized with a 1 A step input. These tables comparethe response to a step input in the SFA domain with the response in the physicaldomain (using the EMTP solution).Table 3.6The response of a step input using the trapezoidal rule, comparing the SFAsolution with the EMTP solution.SFA EMTPt I¯(t) V¯ (t) i(t) v(t)0 0 0 0 0∆t 1e−jωot (2L∆t + jωoL)e−jωot 1 2L∆t2∆t 1e−jωot (−2L∆t + jωoL)e−jωot 1 −2L∆t3∆t 1e−jωot (2L∆t + jωoL)e−jωot 1 2L∆t4∆t 1e−jωot (−2L∆t + jωoL)e−jωot 1 −2L∆tTable 3.7The response of a step input using the backward Euler rule, comparing theSFA solution with the EMTP.SFA EMTPt I¯(t) V¯ (t) i(t) v(t)0 0 0 0 0∆t 1e−jωot ( L∆t + jωoL)e−jωot 1 L∆t2∆t 1e−jωot (jωoL)e−jωot 1 03∆t 1e−jωot (jωoL)e−jωot 1 04∆t 1e−jωot (jωoL)e−jωot 1 0633.6. The Response to a Step Input in the SFA DomainThe energization of a source into an inductor is shown in Figure 3.11. The EMTPand SFA responses are given in Figure 3.12 through Figure 3.14. The EMTP’sresponse using the trapezoidal rule is given in Figure 3.12(a) and the EMTP’s re-sponse using the backward Euler rule is given in Figure 3.12(b). Figure 3.13 showsthe SFA response using the trapezoidal rule, with Figure 3.13(a) giving the responsein the shifted domain and Figure 3.13(b) giving the response in the physical domain.Figure 3.14 shows the SFA response using the backward Euler rule, where Figure3.14(a) shows the response in the shifted domain and Figure 3.14(b) shows the re-sponse in the physical domain. ( )Lv t ( )Li t ( ) 1Ai t 0t Figure 3.11: A DC current source of I=1 A energizes an inductor of L=1 mH.In Figures 3.13(a) and 3.13(b), the numerical oscillations that occur with the trape-zoidal rule in SFA match the stability response of Table 3.6. In Table 3.6, we seethat the trapezoidal response has sustained oscillations in both the SFA and theEMTP responses, where in the shifted domain, the sustained oscillations occur at(2L∆t + jωoL)e−jωot and (−2L∆t + jωoL)e−jωot whereas in the physical domain with theEMTP solution, the oscillations occur at 2L∆t and−2L∆t .Using a time-step of ∆t = 0.1 ms, the oscillations occur at ±20 + j0.377 and usinga time-step of ∆t = 1 ms, the oscillations occur at ±2 + j0.377. Because both thereal and imaginary parts have sustained oscillations, these are seen in both the realand imaginary components of the shifted signal (when the real part has a peak of±20, the imaginary portion has a peak of ±0.377). When transformed back into643.6. The Response to a Step Input in the SFA Domainthe physical domain, these oscillations correspond to sustained beats between thereal and imaginary part of the physical continuous-time domain solution peakingat ±20 or ±2 (Figure 3.13(b)). Contrastively, the EMTP trapezoidal response inFigure 3.12(a), has sustained oscillations at 2L∆t .0 0.002 0.004 0.006 0.008 0.01−25−20−15−10−5051015202530Time (s)Voltage across L (V)EMTP Response to a Step FunctionEMTP -Trap Δt=1msEMTP- Trap Δt=0.1ms-20202-2(a) Trapezoidal response in the EMTP0 0.002 0.004 0.006 0.008 0.010246810Time (s)Voltage across L (V) EMTP - BE Δt=1msEMTP - BE Δt=0.1ms101(b) Backward Euler response in the EMTPFigure 3.12: The voltage across an inductor in response to a 1A DC current in the EMTP:(a) the response given when using the trapezoidal rule and (b) the response given whenusing the backward Euler rule.653.6. The Response to a Step Input in the SFA DomainFigures 3.14(a) and 3.14(b) above show the voltage across the inductor in responseto the step current for two ∆t′s: ∆t = 0.1 ms and ∆t = 1 using the backward Eulerrule. Table 3.7 shows that the backward Euler response has sustained oscillations atjωoLe−jωot after the first two time-steps in the SFA domain, whereas in the EMTP,the response is zero after the second time-step.Using the larger time-step of ∆t = 1 ms, in the shifted domain, the response sus-tains a ripple in the real and imaginary part with a peak-peak value of −0.07→ 0.07(an amplitude of 0.14). The sustained ripples match the sustained imaginary com-ponent of the oscillations in Table 3.7. When ∆t = 1 ms and L=1 mH, evaluating(jωoL)e−jωoL gives the amplitude of 0.14. In the inverse transformed physical re-sponse, these ripples are sustained and instead of dropping to zero after the firsttime-step, the response is held at 0.07, due to the error in the backward Euler so-lution at this time-step. However, when the time-step is less than 0.5 ms, such asin the response when using ∆t = 0.1 ms, because the ripples are smaller (peak-peakvalue of −0.007 → 0.007), in the transformed physical response, the response issustained at approximately 0, just as in the EMTP (Figure 3.12(b)).663.6. The Response to a Step Input in the SFA Domain0 0.002 0.004 0.006 0.008 0.01−20−100102030Voltage across L (Real) SFA Response to a Step Function using Trapezoidal(in the shifted domain) Real part - Δt=1msReal part - Δt=0.1ms0 0.002 0.004 0.006 0.008 0.01−0.4−0.200.20.40.6Time (s)Voltage across L (Imag.) Imag. part - Δt=0.1msImag. part - Δt=1ms2020.377(a) Trapezoidal response in the SFA domain.0 0.002 0.004 0.006 0.008 0.01−20−1001020Voltage across L (Real) SFA Response to a Step Function using Trapezoidal(in the physical domain) Real part - Δt=0.1msReal part - Δt=10ms0 0.002 0.004 0.006 0.008 0.01−20−1001020Time (s)Voltage across L (Imag.) Imag part - Δt=0.1msImag part - Δt=10ms20-20-20202-22-2(b) Trapezoidal response shifted back to the physical domain.Figure 3.13: The voltage across an inductor in response to a 1A DC current using thetrapezoidal rule in SFA: (a) the response given in the shifted domain and (b) the responsegiven in the physical domain.673.6. The Response to a Step Input in the SFA Domain0 0.002 0.004 0.006 0.008 0.010246810Voltage across L (Real) SFA Response to a Step Function using Backward Euler(in the shifted domain) Real part - Δt=0.1msReal part - Δt=1ms0 0.002 0.004 0.006 0.008 0.01−0.0500.05Time (s)Voltage across L (Imag.) Imag. part - Δt=0.1msImag. part - Δt=1ms0.07~0110(a) Backward Euler response in the SFA domain.0 0.002 0.004 0.006 0.008 0.010246810Voltage across L (Real) SFA Response to a Step Function using Backward Euler(in the physical domain) Real part - Δt=0.1msReal part - Δt=1ms0 0.002 0.004 0.006 0.008 0.0100.10.20.30.4Time (s)Voltage across L (Imag.) Imag part - Δt=0.1msImag part - Δt=1ms0.3771010.009 00.07 00.377(b) Backward Euler response shifted back to the physical domain.Figure 3.14: The voltage across an inductor in response to a 1A DC current using thebackward Euler rule in SFA: (a) the response given in the shifted domain and (b) theresponse given in the physical domain.683.7. An RL Circuit Energized with a DC Source Using the Traditional Phasor, the EMTP and SFA3.7 An RL Circuit Energized with a DC Source Usingthe Traditional Phasor, the EMTP and SFAIn this scenario, a 10V DC source energizes an RL circuit, as shown in Figures3.15 and 3.16. The response of the SFA and the EMTP solutions for the voltageacross the inductor and current through the inductor are given in Figures 3.17a and3.17(b). Because the time-constant of the circuit is 2 ms, the chosen time-step is0.2 ms, using the rule that τ ≈ 110∆t. Both trapezoidal and backward Euler rules areused for both the SFA and the EMTP solutions, with the EMTP being the referencesolution for the SFA response. As analysed in Section 3.6, the response to a stepfunction using the trapezoidal rule is accurate but contains sustained oscillations,whereas the response to a step function using the backward Euler rule, after it istransformed from the SFA domain back to the physical domain, is accurate onlywith small time-steps (∆t < 0.5 ms). This observation is verified in Figure 3.17 forwhen using the backward Euler rule, there is an error of 2.77%, whereas using thetrapezoidal rule, there is no error. The error in the backward Euler drops to zerowhen the time-step is decreased, for example, if ∆t = 0.01 ms. Figure 3.15: An RL circuit energized with a DC voltage source. 10R eq2LRt ( )he t ( ) 10Vv t 0t 10R eq o2 LZ j Lt ( )hE t ( )LV t ( )LI t ( ) 10VV t ( )Lv t ( )Li t 0t Figure 3.16: The RL circuit discretized with the EMTP (left) and the SFA (right).693.7. An RL Circuit Energized with a DC Source Using the Traditional Phasor, the EMTP and SFA0 0.01 0.02 0.03 0.04 0.05−1012345678910Time (s)Voltage (V)Voltage across an inductor when energized by a DC SourceEMTP (Trap) Δt=0.2msEMTP (BE) Δt=0.2msSFA (Trap) Δt=0.2msSFA (BE)1 Δt=0.2msSFA (BE)2 Δt=0.01msEMTP(Trap&BE), SFA(Trap&BE2), V=0VSFA(BE1), V=0.277V(a) Voltage across the inductor in response to a DC source.0 0.01 0.02 0.03 0.04 0.05−0.100.10.20.30.40.50.60.70.80.911.1Time (s)Current (A)Current through an inductor when energized by a DC SourceEMTP (Trap) Δt=0.2msEMTP (BE) Δt=0.2msSFA (Trap) Δt=0.2msSFA (BE)1 Δt=0.2msSFA (BE)2 Δt=0.01msEMTP(Trap&BE), SFA(Trap&&BE2), I=1ASFA(BE1), I=0.9723A(b) Current through the inductor in response to a DC source.Figure 3.17: A comparison of the trapezoidal and backward Euler rules in response to a DCsource, using the EMTP as the reference solution.70Chapter 4Transient Stability AnalysisTransient stability, or rotor-angle stability, is concerned with the ability for therotor angles of the generators to remain in synchronicity with each other within3-5 seconds after a severe disturbance has occurred [44]. These studies analyze thebehaviour of the system after the disturbance disrupts the power flow balance inthe electrical network. During such disturbances, the generator will accelerate ordecelerate due to the imbalance of power between the input mechanical power andthe output electrical power. As the machine’s velocity changes, the generator’s rotorfield windings will deviate in position with respect to the generator’s stator windings.If the stator is taken as the inertial reference frame, the relative position betweenthe rotor and stator is δrel [20]. δrel is the angle observed during transient stabilitystudies (as show in Figure 4.1) to determine if the generator is “first-swing” unsta-ble. First-swing instability means that δrel will steadily increase without swingingback to its initial position and it is the assumption used when the synchronous ma-chine is modelled using the equivalent voltage behind reactance model as describedin Section 4.1.2 [4]. When there are multiple generators in the same synchronousarea, the stability depends on the synchronicity among all the machines. The gener-ators will accelerate and decelerate at different rates due to several external factors,such as their location to the fault and their internal parameters (such as reactances,constant of inertia, MVA capacity). Each generator’s relative angle δrel will increasewith respect to their own stator windings and with respect to each other. If therelative angle between two generators increases beyond the transient stability limit[4], the synchronicity between these two machines will be lost and the entire sys-tem will be considered unstable. The main objective of transient stability analysisstudies is to ensure synchronism among all network components, in particular therotating masses of the synchronous machines.Due to the imbalance of power that creates the asynchronous behaviour of thegenerators, it is relevant to understand the role of power equilibrium between the714.1. A Physical RepresentationX d 'g e nIgenE Vs'dXgenIgenE where δ isstatorrotorFigure 4.1: The rotor angle of the generator is the relative angle between the stator windingand the rotor’s field winding, taking the stator as the reference frame and using the simplifiedmodel of the synchronous generator (Section 4.1.2).mechanical power from the turbine into the generators and the electrical powerdelivered to the network from the generators, as briefly discussed in Section 1.2.Power equilibrium can be described using both the physical representation of thepower imbalance and the mathematical model that describes the physics.4.1 A Physical RepresentationAs described in Section 1.2, the ability for the power grid to maintain synchronicitydepends on the power balance between the mechanical power into the generatorsfrom the turbine and the electrical power delivered to the network. During steady-state operation, the mechanical power provided by the turbines driving the syn-chronous generators is equal to the electrical power consumed by the network (plusthe mechanical and electrical power losses): Pelec = Pmech + Plosses. However, sud-den disturbances on the grid changes the power balance and Pelec 6= Pmech +Plosses.The input mechanical power supplied by the prime mover (e.g., hydro or thermo),Pturbine is converted into both electrical power that feeds the network, Pelec, andstored kinetic energy in the generators: Pturbine ≈ Pelec + dEkineticdt [20], where in nor-mal steady-state conditions, Pturbine = Pelec anddEkineticdt remains as stored energy.In the first few seconds after a disturbance, the mechanical output power Pmech isslow to react to the changes of the electrical power due to the difference in time-constants (the mechanical time constants are larger than the electrical ones [45]) and724.1. A Physical RepresentationPturbine remains constant for about34 of a second before control action is initiated,such as automatic voltage control [4][46]. There are two conditions that can changethe power balance: (1) If there is an increase in electrical power, for example, due toa large increase in load or a sudden loss of generation and (2) if there is a decreasein electrical power, for example, due to a loss of load or a fault.In the first condition, when there is a large increase in load or loss of generation, thereis a lack of Pturbine to supply Pelec. During this time, the power is borrowed fromthe generator’s kinetic energy Pkinetic, which slows the machine down, or decreasesits velocity away from the 60-Hz frequency. Conversely, when there is a loss of loaddue to, for example, a fault, there is an excess of Pturbine compared with Pelec. Theextra power is transferred to the kinetic energy of the rotor, causing the rotor toaccelerate, which increases the electrical frequency.4.1.1 The Equal Area Criteria and The Critical Clearing TimeThe physical interpretation can be illustrated using the equal area criteria, whichis a graphical representation of the kinetic energy balance in the network. If thenetwork is represented as an equivalent single generator feeding an infinite bus [4],then the equal area criteria can be used to determine the maximum rotor angleposition before instability occurs. As illustrated in the equal-area curve ofFigure 1.2(b), when the fault occurs, the electrical power Pe will decrease and theangle will begin to increase. When the fault is cleared, the electrical power will in-crease to its post-fault value and the angle will continue to increase due to the extrakinetic energy stored. The areas A1 and A2 represent the stored kinetic energy.The system is stable if A1<= A2.The critical clearing angle of the machine is δc and the time at which the rotorangle reaches this position is called the critical clearing time (CCT). The CCT isa fundamental parameter necessary to compute when the disturbance is a fault asthe CCT is the greatest time that the circuit breakers have to open to ensure thatthe system remains stable. When clearing a fault, the total clearing time neededis the total time it takes for the relay to detect the fault, for the relay to signalto the circuit breaker to open, for the circuit breaker to begin opening its contactsand create an arc, which will extinguish the fault, plus a margin of time delay [47].The CCT is expressed in cycles on the basis of the grid’s frequency at 60 cycles per734.2. A Mathematical Representationsecond [47]. For example, the breakers are designed as 5-cycle or 8-cycle breakersif the CCT is computed as 5 cycles or 8 cycles [47]. One reason for this is that thecircuit breaker is set to open when the arc energy is as low as possible, which willoccur after the fault current has crossed zero, any time between two 60-Hz cycles.4.1.2 The Approximated Synchronous Machine ModelAs the main motivation for this work is to model the SFA solution for the electricalnetwork compared with the traditional phasor and the EMTP, the simplified modelfor the synchronous machine is used in this thesis. The simplified generator modelis an approximation of the full-order synchronous machine model and is reasonablyaccurate for the first few cycles of the transient disturbance [4][19]. The generator ismodelled as a voltage source behind a transient reactance, where the rotor angle δrelis the same electrical phase angle of the generator (Figure 4.1). This representationassumes that the rotor flux linkages are constant, the damper winding currents havedied down, and saliency is neglected [19]. Because the time constants of the sub-transient and armature currents are in the order of 0.03 s→0.04 s and 0.1 s→0.3 s[48], whereas the electrical transient phenomenon ranges from 0.5 s→10 s [48], it canbe assumed that only the transient reactance x′d remains.4.2 A Mathematical RepresentationThe core mathematical equation that models the machine’s dynamics is formulatedin Newton’s second law of motion for rotational bodiesTnet = IαTnet = Id2δreldt2,(4.1)where Tnet is the net torque between the mechanical and the electrical network, Iis the machine’s constant of inertia (kg·m2), and α is the rotational acceleration( rads2). This equation, commonly referred to as the “Swing Equation”, is a second-order non-linear equation, which can be expressed more intuitively, as two first-orderequations.744.2. A Mathematical RepresentationAs a second-order equation,Tm(t)− Te(t) = I d2δrel(t)dt2+Kdδrel(t)dt(4.2)As two first-order equations,Tm(t)− Te(t) = I dωrel(t)dt+Kωrel(t), (4.3)ωrel(t) =dδrel(t)dt(4.4)In (4.2),(4.3), and(4.4), δrel is the rotor angle (rad), ωrel is the relative velocity be-tween the rotor velocity and the rotating field created by the stator currents ( rads ),Tm(t) is the torque produced by the prime mover moving the field winding (N·m),Te(t) is the torque produced by the rotating field due to the stator currents (N·m),I is the moment of inertia (kg·m2), and K is the self-damping coefficient of themachine(N·mrads), which ranges from 0.5-2.0 per-unit of the turbine rating for hydroturbines [21]. The damping coefficient allows the rotor’s oscillations to decrease inamplitude and to die out [48]. In stability studies, the inertia constant of the ma-chine H is commonly used, where H = Iω2o2Sbase(s), and Sbase is the MVA rating of thegenerator. Equation (4.4) is used for both the mechanical and the electrical veloc-ity of the generator, where the relationship between the mechanical and electricalvelocity is given as,wmech =2welecp, (4.5)where p is the number of poles in the machine. If the machine is two-poled, themechanical and electrical frequency are equal wmech = welec.Equations (4.2) and (4.3) are non-linear because net power is known in the powernetwork and the frequency, ω is used to convert net torque into net power. In thetraditional transient stability programs, the net torque is converted to net powerthrough ωo under the assumption that the actual frequency ω does not deviate farfrom ωo. This assumption is approximately accurate only under small frequencydeviations where ω ≈ ωo. When the frequency deviations are larger, ω will differfrom ωo. For a more precise mathematical model, the correct ω should be used in(4.2) and (4.3) to convert net torque into net power, as done in this thesis.754.2. A Mathematical Representation4.2.1 The Coupling of the Electromechanical and ElectricalNetwork Equations in the Solution MethodsTo mathematically couple the electromechanical and the electrical networks, thephasor solution, the EMTP, and the SFA solution methods take different approaches.Coupling in the Traditional Phasor MethodThe traditional simulation programs (e.g.,[3][18]) that model the electrical networkwith the 60-Hz frequency phasor solution method use step-by-step time-domainmethods to simultaneously solve the electromechanical and electrical network equations.ddt= f(x, V ),I(x, V ) = YNV,(4.6)where x is the vector of states (rotor angle, frequency), V is the vector of voltages,I is the current vector, and YN is the admittance matrix [3].The differential electromechanical equation (4.2) is solved using continuous-timenumerical integration techniques, such as the Runge-Kutte or the Trapezoidal tech-nique [3] and the algebraic network equations (4.6) use the traditional phasor method,described in Section 1.3.1. In actuality, the phasor method is solved in the frequencydomain through the phasor representation [21].Coupling in the EMTP and SFA MethodsThe EMTP solution is solved completely in the time-domain. The electromechanicalequation is discretized using numerical discretization methods, such as the trape-zoidal or the backward Euler, given in (4.9) and(4.10) or (4.11) and (4.12), andthe electrical network equations are solved through the discrete-time EMTP equiv-alent discretization, as described in Section 1.3.2. Likewise, in the SFA solution,the electromechanical equation is discretized like the EMTP ((4.9), (4.10), (4.11),and (4.12)), and the electrical network equations are solved in the SFA domain asa discrete-time solution, as described in Chapter 2. The coupling of the differentialand algebraic equations is completely in the discrete-time domain when using theEMTP and the SFA solutions due to the discretization of the continuous-time equa-tions, as opposed to the phasor solution, where the equations are coupled betweenthe continuous-time domain and the phasor domain.764.2. A Mathematical Representation4.2.2 The Discretization of the Electromechanical EquationIn this work, the actual velocity ω is used to convert between torque and power.Since the electromechanical equation is discretized using numerical discretizationin the solution method, ω is known at every time-step. If the trapezoidal rule isused in this discretization, (4.3) and (4.4) are formulated as (4.9) and (4.10). If thebackward Euler discretization rule is used in this discretization, (4.3) and (4.4) areformulated as (4.11) and (4.12). Net torque is converted into net power usingTnet(t) =Pnet(t)wrel(t−∆t) , (4.7)wherePnet(t) = Pm(t)− Pe(t). (4.8)The following equations discretize the electromechanical equation using the trape-zoidal and the backward Euler discretization rules.Applying the Trapezoidal Discretization Ruleωrel(t) =1− KI ∆t21 + KI∆t2(ωrel(t−∆t))+1I∆t21 + KI∆t2(Pm(t)− Pe(t)ωrel(t−∆t) ) +1I∆t21 + KI∆t2(Pm(t−∆t)− Pe(t−∆t)ωrel(t−∆t) )(4.9)andδrel(t) =∆t2(ωrel(t) + ωrel(t−∆t)) + δrel(t−∆t). (4.10)Applying the Backward Euler Discretization Ruleωrel(t) =11 + K∆tI(∆tI(Pm(t)− Pe(t)ωrel(t−∆t) )) +11 + K∆tI(ωrel(t−∆t)) (4.11)andδrel(t) = δrel(t−∆t) + ωrel(t)∆t. (4.12)The backward Euler discretization rule is used for the electromechanical equationsin this work.774.3. Implementation of SFA in Three Transient Stability Case Studies4.3 Implementation of SFA in Three TransientStability Case StudiesThe SFA solution method is implemented and tested in three power system tran-sient stability case studies: a Single-Generator Infinite-Bus (SGIB) network [4](Chapter 5), a three-bus network [45] (Chapter 6), and the IEEE 39-bus network [1](Chapter 7). In these case studies, SFA’s capability to computationally and ef-ficiently capture network dynamics is demonstrated during three-phase balancedshort-circuit faults. In the three-bus network, two generators of different inertiaconstants feeding a constant impedance load are modelled with the different meth-ods and with the IEEE 39-bus network, the feasibility of the SFA solution methodis demonstrated in a large-scale network. The three-bus and the 39-bus networksalso demonstrate the behaviour of the solution methods when a generator is discon-nected from the network. For all case studies, the EMTP is used as the referencesolution.In the three case studies, some of the observations considered includes:1. The similarities and differences in CCTs, in both 60-Hz cycles and milliseconds,among the three methods when different three-phase balanced short-circuitfaults occur on various locations in the network.2. The similarities and differences among the three solution methods after a three-phase fault has been cleared from the network.3. The behaviour of the solution methods when one generator is suddenly dis-connected in the three-bus network and the 39-bus network.4. The differences between the trapezoidal and backward Euler discretizationmethods in the SFA solution.5. The difference between using the approximate vs actual velocity in the elec-tromechanical equation for the traditional phasor and the SFA solution meth-ods.The commercial Transient Security Assessment Tool (TSAT) software [3] is usedfor the phasor solution, the commercial MicroTran software [14] is used for theEMTP, and a self-written program TSFA is used for the SFA solution method. A784.3. Implementation of SFA in Three Transient Stability Case Studiesphasor-solution solver, TPhasor, was also written, which allows the comparison ofthe correct network frequency with the approximate network frequency in the elec-tromechanical equation for the traditional phasor solution. Both the TSFA andTPhasor solutions were written in the Python computer language and are general-ized solution algorithms.For the discretization method of the electromechanical equations, the solution fromTSAT uses the trapezoidal numerical integration method, MicroTran uses the trape-zoidal method with critical damping adjustment (CDA), [38] and unless otherwisestated, TSFA uses the backward Euler method for both the electrical and the elec-tromechanical network equations. Backward Euler is used for the TSFA solutioninstead of trapezoidal because the backward Euler rule dampens the numerical os-cillations as explained in Chapter 3. However, as seen in Section 3.4.2, backwardEuler introduces high error around 0 Hz, whereas the trapezoidal behaves more ac-curately at 0 Hz, which is demonstrated when observing the DC component of thefault current.The time-steps chosen for the solution methods were the maximum values that thesolutions would give without too much loss of accuracy in the electromechanical vari-ables. For the EMTP ∆t = 100µs, for TSAT ∆t = 8 ms, and for TSFA ∆t = 8 ms.The time-step for the reconstructed SFA solution (in the continuous-time, physicaldomain) is 100µs to match the EMTP solution. In all the cases, the EMTP is usedas the benchmark solution.Conditions and Necessary AssumptionsThe following conditions were set for the case studies:1. The power system is modelled as a three-phase balanced symmetrical networkand only three-phase faults are applied on the system to maintain the balancedcondition for all three phases. Future work will include unbalanced conditions,such as single line-ground faults.2. The controllers, such as the automatic voltage controller and the turbine’s gov-ernor action, are not considered, and therefore Pm is taken to be constant. Thisassumption is used because the first-swing of the rotor angle occurs within thefirst few seconds after the disturbance [4] and the behaviour of the generators794.3. Implementation of SFA in Three Transient Stability Case Studiesgiven by the solution methods is observed for first-swing stability.3. The simplified generator model is used which consists of a voltage source Egenand a direct-axis transient reactance x′d as described in Section 4.1.2 [4]. Tofacilitate the studies, the generators are assumed to be two-poled (3600 rpm insteady-state) and saliency is ignored. In TSAT, the classical model “CGEN”is used for the generators [3]. In MicroTran, because the time constants of thesub-transient and damper winding currents are small, it can be assumed thatthe it can be assumed that the damper winding currents have died down asthe sub-transient reactance time constants are smaller than the the transientreactance time constants (i.e. T ′′d ≈ 0.035 s, whereas the time constant of thetransient reactance, T′d ≈ 1.5s [48] and so only the transient reactances (x′d )of the synchronous generator are used for the simplified generator model. Tomodel the simplified generator in MicroTran, the machine’s reactances can allbe set to equal x′d ( xd = xq = x′d = x′q = x′′d = x′′q).4. The loads are modelled as constant impedances, where Pload ∝ V 2 andQload ∝ V 2.5. The EMTP solution in MicroTran updates the generator voltage with thefrequency changes in the network; that is, Egen =ω(t)ωoEfield, where Efield is thevoltage induced by the field winding current and is the generator’s voltage Egenand ω is the known grid frequency at that time-step. This is a type of internalcontrol mechanism between the generator’s voltage and the network frequency,which is different from an external control mechanism, such as the automaticvoltage regulator or the governor droop control action. This internal controlmechanism is implemented in the TSFA solver.6. The CCT is observed based on first-swing stability [4][49]; that is, if the firstswing of the rotor angle (or difference in rotor angles between the machinesand chosen reference machine) is unstable after the fault has been clearedthen the system is assumed unstable. The angle stability margin index usedto determine the CCT in TSAT [3] is implemented in TSFA. This margin iscalculated by finding the maximum generator angle separation between anytwo generators and by calculating an index η, which must be between -100and 100 [3]. It is important to note that first-swing stability has relativelygood accuracy only under the assumption that the synchronous machine is804.3. Implementation of SFA in Three Transient Stability Case Studiesmodelled using the constant voltage behind reactance model, as described inSection 4.1.2, and if the critical clearing times are larger than 6 cycles [48].However, the first-swing approximation may not always be reliable, especiallyduring fast clearing times (under 6 cycles)[48]. When the clearing time is veryfast, the transient stability limit may exceed the steady state stability limit,which may not show instability until after the first-swing.[48]. This, however,is not the case in the following case studies, where all clearing times are greaterthan 6 cycles.7. As the circuit breaker must wait until after the fault current has crossed zeroto clear the fault [47], both the TSFA and EMTP methods both use a zero-crossing fault detection algorithm. The fault current zero-crossing detectionmethod implemented in TSFA is described in Section 4.3.1.8. The power flow solution for the network’s initial conditions (initial Q and δ ofthe PV buses and V and δ of the PQ buses) is solved using the commercialproduct, Powerflow & Short-circuit Analysis Tool (PSAT) from Powertechlaboratories[50]. The initial generator voltages Egen and the initial value ofthe rotor angles δo are solved for using Kirchoff’s Voltage LawE¯genn = V¯tnejθn + jx′dnI¯tn , (4.13)whereI¯tn =S¯∗V¯tne−jθnI¯tn =Pn − jQnV¯tne−jθn ,where I¯t is the current out of the generator’s terminals, V¯t is the voltage outof the generator’s terminals, and n is the number of generators in the system.9. The fault voltage is plotted as the peak, line-neutral per-unit voltage. There-fore there is a factor of√2√3when converting from the real values to the per-unitvalues.10. The TSFA and TPhasor solutions refer to the solutions using the velocity ωin the electromechanical equation and TSFA* and TPhasor* use the approxi-mated velocity ωo in the electromechanical equation.814.3. Implementation of SFA in Three Transient Stability Case Studies4.3.1 Fault Current Zero-Crossing Detection Algorithm in TSFAWhen a short-circuit fault occurs on the network, the protective devices, such as thecircuit breaker, must open to isolate the fault. At the time of the fault, the relaydetects the fault and sends a signal for the circuit breaker to open its contacts. Thecontacts are electrodes and as they separate, the medium that the circuit breakeris immersed in, such as oil, vacuum, or SF6, becomes ionized and an arc of energyis formed between the contacts [47]. As long as the arc is sustained in betweenthe contacts, the current cannot be interrupted. The fault current will only be in-terrupted when the energy between the contacts has dissipated and the arc is asthin as possible [20]. This occurs when the current is as low as possible (i.e., whenthe current is zero). In AC systems, the fault current is sinusoidal, changing frompositive to negative every half cycle, and therefore it crosses zero twice in one cycle.The circuit breakers will succeed in clearing the fault when the fault current hascrossed zero.A fault current zero-crossing detection algorithm is implemented in the TSFA algo-rithm to detect the zero-crossing in the shifted domain. Because the shifted domainonly gives the envelope and phase angle of the current, however, detecting whenthe current’s amplitude crosses zero will not be feasible. Therefore, this algorithmdetects when the phase of the fault current crosses 90◦ after the set fault clearingtime. The derivation is shown as follows:i(t) = |I(t)| cos(ω(t) + θ(t))=|I(t)|2(ej(ω(t)+θ(t)) + e−j(ω(t)+θ(t))),(4.14)where ω(t) is the grid frequency, |I(t)| is the magnitude of the fault current, andθ(t) is the phase of the fault current.Current i(t) is zero when cos(ω(t) + θ(t)) = 0◦. This occurs when ω(t) + θ(t) = 90◦.Alternatively, when ω(t)+θ(t) = 90◦, ej(ω(t)+θ(t)) +e−j(ω(t)+θ(t)) = 0, which happenswhen ω(t) and −ω(t) are 180◦ apart from each other.(ω(t) + θ(t))− (−ω(t) + θ(t)) = 180◦ or (ω(t) + θ(t)) = 90◦. (4.15)Therefore, the fault clearing time in the SFA domain occurs the first moment that824.3. Implementation of SFA in Three Transient Stability Case Studiesω(t) + θ(t) = 90◦ after the set clearing time. The algorithm ensures that the angleis phase angle is constrained to 360◦. Without this fault current phase zero-crossingdetection algorithm in the shifted domain, at every time-step the fault current wouldneed to be re-transformed into the physical domain, heavily increasing the compu-tational speed.Considerations Made in the TSFA AlgorithmThe following considerations are made in the SFA solution method for the TSFAalgorithm:1. For the change in network topology, such as the fault occurring and clearing,the history voltage is interpolated between the time-step before the topologicalchange and the present time-step.2. The initial conditions for the SFA solution at time t = 0− are set by solvingthe steady-state phasor solution for the solution at t = 0 instead of assumingzero initial conditions.3. As described above, the fault current must cross zero before the fault is actuallycleared.The SFA solution algorithm is given in Figure 4.2.834.3. Implementation of SFA in Three Transient Stability Case StudiesNoNoIs this the final time-step? YesDetermine angle stability index. Calculate CCT from time that fault was cleared. Interpolation scheme. Interpolate the voltage and current solutions at a smaller time-step.Is the stability index negative? NoDetermine which generator went unstableStore results in output text filesSystem is stable! Define SFA branches and admittance matrix [Y] based on the discretization method (trapezoidal or backward Euler). Solve power-flow solution: Obtain initial conditions for bus voltages (magnitude and angle), real and reactive power for loads and generators. Increase the iteration number by one. NoIs this the first iteration in theGauss-Seidel iteration?Check contingency condition. Reset iteration counterNoYesNoYesDefine a Network Topology File. It contains: Branch Type (resistor, inductor, capacitor, impedance, susceptance)Branch Nodes To and From: node numbers that the branches are attached to, with node 0 as groundBranch values (defined in the case data file)Define the Case Data File . It contains:Initial Conditions: from power-flow solution)Time Conditions: (time-step, total simulation time). Simulation Conditions: discretization technique used (i.e., backward Euler, trapezoidal), Electrical Network: base MVA, base kV, Electrical network parametersElectromechanical network: generators: (machine inertia/MVA, initial generator voltages and angles). Define a Source Data File. It contains: Source Type: voltage or current, AC or DCSource Node To and From: node numbers that the branches are attached to, with node 0 as groundSource initial values: from case data file.Read network topology file. Read case data file.Read source data file.Define steady-state phasor braches and admittance matrix for initial conditonsCreate SFA sources (using transformation: =e oj tTSolve initial conditions: Solve phasor solution for the initial node voltages and current branches atInitialize mechanical variables. Solve for the history current at -t= 0-t=0Increase time-step toUpdate the mechanical variables at t=t+Δtt=t-ΔtYesSolve Electrical Network using the EMTP type solution method [16]: Define matrix of known sources ([VB])Multiply the matrix of known admittances [GAB] with [VB]: [GAB]=[GAB][VB]. Solve for the history terms of the previous time-step at t=(t-Δt)Obtain matrix of history terms [H]. Solve for the node voltage: Solve for the branch currents. 1[ ] [ ] [ ]V GAA HA GAB Has a fault occurred? Calculate and change admittance matrix with faulted condition YesHas the circuit breaker opened? Is this the first time-step for the fault? YesNoYesHas the fault current crossed zero? YesCalculate and change to new admittance matrix for post-fault contingency.Interpolate the histories for the next time-step. , Did a generatordisconnectSolve mechanical network: Obtain voltage and current of generatorsSolve for P, Q of each generator.NoNoCalculate torque: ( )PTtSolve for the frequency and rotor angle at t=ΔtIs convergence error less than tolerance?YesCalculate and change to admittance matrix for the faulted generator condition. Start post-processing dataPre-ProcessingMain Solution LoopPost-ProcessingShift the SFA solution back to the physical domain: Apply the inverse transformationTake the real partFigure 4.2: The TSFA solution algorithm. 84Chapter 5Single-Generator Infinite-BusCase StudyThe Single-Generator Infinite-Bus network is shown in Figure 5.1 [4], where a25 MVA, 60-Hz generator delivers 20 MW through two parallel transmission linesto a large electrical network, represented as an equivalent infinite bus. The networkparameters are provided in Table 5.1, with a system base of 25 MVA, 25 kV. 20MWeP basebase25kV25MVAVS Generator 1 ' 0.3pudX 12 0.2puX 13 0.1puX 23 0.1puX Infinite Bus Figure 5.1: The single-line diagram of the Single-Generator Infinite-Bus (SGIB) network,where the electrical network is represented as an infinite bus. A three-phase balanced faultis applied to the terminals of the generator.Table 5.1Network parameters for the SGIB network (Figure 5.1).GeneratorPowerRating(MVA)VoltageRating(kV)TransientReactancex′d (p.u.1)InertiaConstant (H)(sec)Moment ofInertia (I)(kg·m2)DampingCoefficient(K)(p.u.2)Generator 1 25 25 0.3 2.76 970.99 1.5Generator’sInitialConditionsReal Power(MW)Reactive Power(MVAR)Terminal VoltageVt(per-unit1∠(deg))Eo∠δo(p.u.1∠(deg))Generator 1 20 0.7246 1.032∠4.3◦ 1.067∠16.88◦TransmissionLinesReactance(Ω/km)Reactance(p.u.3)Line 1-2 5 0.2Line 1-3 2.5 0.1Line 2-3 2.5 0.11 per-unit of generator base 2 per-unit of torque base 3 per-unit of system base 855.1. Critical Clearing Times5.1 Critical Clearing TimesThe CCTs given by the traditional phasor solution (TSAT), the traditional phasorsolution with the exact velocity (TPhasor), the EMTP, TSFA, and TSFA* (whichuses the approximate velocity ωo) are provided in Table 5.2 for the three-phase faultat the generator’s terminals. The CCT is in both 60-Hz cycles and in millisecondsand the critical clearing angle δc is in degrees. To be conservative, the 60-Hz cy-cle shown is rounded down to the lowest closest cycle (for example, in the EMTP,the rotor angle becomes unstable at 13.6 cycles, which is rounded down to 13 cycles).Table 5.2The CCTs for the three-phase fault on the generator’s terminals of theSGIB system in Figure 5.1 in 60-Hz cycles and milliseconds, and the criticalclearing angle in degrees.Three PhaseFault on Gen.1EMTP(∆t = 100µs)TSAT(∆t = 8ms)TPhasor(∆t = 8ms)TSFA(∆t = 8ms)TSFA*(∆t = 8ms)tc(cycles) 13∗ 14 14 14 14tc(ms) 216 233 234 248 248δc 80◦ 104.8◦ 106.3◦ 97.5◦ 98.1◦* rounded down from 13.8 to 13 cyclesFrom Table 5.2, the following observations are made:1. There is a half-cycle difference in the CCT between all the solutions(TSAT,TPhasor, TSFA and TSFA*) and the reference EMTP. This is due to thesmaller time-step used in the EMTP solution compared with the other solu-tions. If the EMTP time-step is raised to ∆ = 1 ms, then the CCT for theEMTP solution matches the other solutions at 14 cycles.2. Using the actual ω instead of ωo in the electromechanical equations does notchange the CCT for the SFA solution method; however, the traditional phasorsolution changes by 1 ms.865.2. A Three-Phase Fault at the Generator’s Terminals5.2 A Three-Phase Fault at the Generator’s TerminalsA three-phase balanced fault is applied to the terminals of the generator (Bus 1 inFigure 5.1) at t = 1 s and is cleared after 8 cycles. The generator’s electromechani-cal variables: rotor-angle, frequency and electrical power as well as the fault voltageand fault current are simulated using the phasor, the EMTP and the SFA solutionmethods, with the EMTP as the reference solution.The fault is set to clear at 8 cycles, but due to both the differences in the solutionmethods’ time-steps and the zero-crossing fault current detection employed in theEMTP and TSFA solutions, the actual clearing of the fault is different among themethods, varying between 8-9 cycles. This difference in clearing times means thatextra kinetic energy has accumulated in the network [20], producing both a phasedisplacement and a magnitude difference among the methods.When the fault happens at t = 1 s, the generator’s rotor angle and frequency beginto increase until the fault is cleared. When the fault is cleared, the frequency beginsto decrease, but due to extra kinetic energy stored in the generator, the angle con-tinues to increase for another 4.4 cycles (73.3 ms). Because the fault is cleared beforeits CCT of 14 cycles, the system is considered “first-swing stable”. The infinite busacts like a constant voltage source with zero impedance and infinite inertia [4], fixingthe system’s frequency to 60 Hz around the infinite bus. Figures 5.2 and 5.3 presentthe generator’s electromechanical variables: rotor angle, frequency and electricalpower, where in Figure 5.3, the TSFA solution uses both ∆t = 5 ms and ∆t = 8 msto compare the differences in time-steps with respect to the EMTP reference. InFigures 5.2(c) and 5.3(c), the electrical power oscillations in the EMTP solution aredue to the EMTP calculating the instantaneous power, whereas the phasor and SFAsolution calculate the average power, P = R(V I∗).The fault voltage and fault current are shown in Figures 5.4 and 5.5, where thedifference in clearing times are prominent. In the EMTP, the fault is cleared at8.46 cycles (140.6 ms) after zero-crossing detection and in the SFA solution, thefault is cleared at 8.63 cycles (143.9 ms) after zero-crossing detection. In the TSATsolution, the fault is set to clear at exactly 8 cycles (133 ms), while in the TPhasor*solution, used to simulate the fault current, due to the fixed discrete time-step, thefault current clears in 8.15 cycles (136 ms).875.2. A Three-Phase Fault at the Generator’s Terminals0.5 1 1.5 2 2.5 3−200204060Rotor Angle (Deg)Rotor Angle of Generator 1EMTP Δt=100µsTSFA-BE Δt=8ms TSAT Δt=8msFault cleared after 8 cyclesFault at t=1sSee Fig.5.2(a) EMTPTSFATSAT0.5 1 1.5 2 2.5 35859606162Frequency (Hz)Frequency of Generator 1EMTP Δt=100µsTSFA-BE Δt=8ms TSAT Δt=8msFault at t=1sFault cleared after 8 cyclesSee Fig.5.2(b) EMTPTSFATSAT0.5 1 1.5 2 2.5 3−200204060Time (s)Electrical Power (MW) Electrical Power of Generator 1EMTP Δt=100µsTSFA-BE Δt=8ms TSAT Δt=8msFault cleared after 8 cyclesFault at t=1sSee Fig.5.2(c) EMTPTSFATSATFigure 5.2: A comparison of the TSFA and TSAT solutions compared with the EMTPsolution for the electromechanical variables in the system of Figure 5.1 during a three-phasefault at the terminals of the generator.1.6 1.730405060Time (s)Rotor Angle (Deg)EMTP TSFA-BETSATEMTPTSFATSATFig.5.2 (a) Rotor angle1.5 1.660.56161.5Time (s)Frequency (Hz)EMTPTSFA-BETSATEMTP(61.38 Hz)TSFA(61.31 Hz)TSAT (61.43Hz)Fig.5.2 (b) Frequency1.6 1.7 1.8102030405060Time (s)Electrical Power (MW)EMTPTSFA-BETSATEMTPTSFATSATFig.5.2 (c) Electrical Power885.2. A Three-Phase Fault at the Generator’s TerminalsIn Figure 5.2(b), there is a difference in the peak-value of the generator’s frequencyamong the methods. The frequency error between TSAT and the EMTP is 0.081%,which is smaller than the error between TSFA and the EMTP on the first swing,0.114%. However, at the second-swing and onwards, the TSFA and the EMTP so-lution are consistently closer to each other than the TSAT and the EMTP solution(0.106% between TSFA and the EMTP and 0.146% between TSAT and the EMTPon the second-swing).By using a smaller time-step of ∆t = 5 ms for TSFA, as shown in Figure 5.3, theerror between the TSFA and the EMTP solution is less than when using ∆t = 8 msin Figure 5.2. This is due to two factors. First, when the smaller time-step is used,the fault is cleared at 8.39 s (139 ms) (Figure 5.3), which is a difference of 1 ms be-tween the TSFA and the EMTP, whereas when ∆t = 8 ms, the fault was cleared at8.63 cycles (143 ms) in the TSFA solution (Figure 5.2(b)). This is a 4 ms betweenthe TSFA and the EMTP. When the fault is cleared in the TSFA solution at a closertime to when the fault is cleared in the EMTP solution and the time-step used issmaller, the frequency error between TSFA and the EMTP on the first-swing of thefrequency is smaller (0.016% in Figure 5.3(b)). There is no difference when using asmaller ∆t in TSAT. The TSAT and TSFA solutions are close to each other in thiscase study because the infinite bus does not allow large frequency deviations as itfixes the network frequency close to 60 Hz.A comparison of the error in frequency between the three solution methods is pro-vided in Table 5.3.Table 5.3A comparison of the error in frequency among TSAT (at ∆t =8ms) andTSFA (at ∆t =5ms and ∆t =8ms) with respect to the EMTP.TSAT ( ∆t = 5 ms or ∆t = 8 ms) TSFA (∆t = 8 ms) TSFA (∆t = 5 ms)frequency error(first peak) (%)0.083 0.116 0.016frequency error(second peak) (%)0.146 0.106 0.016clearing time(first peak) (cycles)8.15 8.63 8.39895.2. A Three-Phase Fault at the Generator’s Terminals1 1.5 23040506070Time (s)Rotor Angle (Deg)Rotor Angle of Generator 1EMTP Δt=100µs(a) TSFA-BE Δt=5ms (clear at 8.39 cycles)(b) TSFA-BE Δt=8ms (clear at 8.15 cycles)TSAT Δt=8msTSAT Δt=5ms 95100EMTPTSFA(b)TSFA(a)TSATSee Fig.5.3(a)1.5 2 2.560.860.96161.161.261.361.461.5Time (s)Frequency (Hz)Frequency of Generator 1EMTP Δt=100µs(a) TSFA-BE Δt=5ms (clear at 8.39 cycles)(b) TSFA-BE Δt=8ms (clear at 8.15 cycles) TSAT Δt=8msTSAT Δt=5ms 95100EMTPTSFA(b)TSFA(a) TSATSee Fig.5.3(b)1.2 1.4 1.6 1.8 2 2.2203040506070Time (s)Electrical Power (MW) Electrical Power of Generator 1EMTP Δt=100µs(a) TSFA-BE Δt=5ms (clear at 8.39 cycles)(b) TSFA-BE Δt=8ms (clear at 8.15 cycles) TSAT Δt=8msTSAT Δt=5msEMTPTSFA(b)TSFA(a) TSATSee Fig.5.3(c)Figure 5.3: A comparison of the EMTP, TSFA using ∆t = 5 ms and ∆t = 8 ms and thephasor solution using ∆t = 5 ms and ∆t = 8 ms for the electromechanical variables duringthe three-phase fault in the system of Figure 5.1.1.730405060Time (s)Rotor Angle (Deg)EMTP (a)TSFA-BE Δt=5ms (b)TSFA-BE Δt=8ms TSAT Δt=8msEMTP(b)TSFATSAT (a)TSFAFig.5.3 (a) Gen.1’s Rotor angle.1.5 1.660.56161.5Time (s)Frequency (Hz)EMTP (a) TSFA-BE Δt=5ms (b) TSFA-BE Δt=8ms (clear at 8.15 cycle) TSAT Δt=8msEMTP(61.38 Hz)(b) TSFA(61.31 Hz)TSAT (61.43Hz)(b) TSFA(61.39 Hz)Fig.5.3 (b) Gen.1’s Frequency.1.6 1.7 1.8102030405060Time (s)Electrical Power (MW)EMTP(a) TSFA-BE Δt=5ms(b) TSFA-BE Δt=8msTSAT Δt=8msEMTPTSAT(b)TSFA(a)TSFAFig.5.3 (c) Gen.1’s Power.905.2. A Three-Phase Fault at the Generator’s TerminalsThe fault voltage presented in Figure 5.4 and Figure 5.4(a) shows that the continuous-time TSFA solution matches the EMTP solution and the TSFA envelope is almostidentical to the TSAT solution, except for the first time-step post-fault due to thedifferences in fault clearing times.1 1.1 1.2 1.3 1.4−0.8−0.6−0.4−0.200.20.40.60.8Time (s)Voltage(per-unit of system base)Voltage at the faulted bus (Terminals of Generator 1)EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTSAT Δt=8msTSFA (Envelope)TSAT TSFA (Reconstructed)EMTPDue to different fault clearing times Figure 5.4: The voltage at the faulted bus (terminals of Generator 1) during the 3-phasefault.1.42 1.43 1.44 1.45 1.46 1.47Voltage at the faulted Bus (Pre-Fault and Post-Fault)EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTSAT Δt=8msTSFA (Envelope)TSAT TSFA(Reconstructed)EMTPTime (s)0.8 0.82 0.84 0.86 0.88 0.90.810.820.830.840.850.860.870.88Voltage (per-unit of system base) EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTSAT Δt=8msTSFA (Envelope)TSAT TSFA (Reconstructed)EMTPTime (s)Figure 5.4 (a) A close-up of the voltage at the faulted bus shows that TSAT and the envelope ofTSFA match and the reconstructed TSFA and the EMTP match. (Left: pre-fault, Right: post-fault).915.2. A Three-Phase Fault at the Generator’s TerminalsThe fault current is simulated using both ∆t = 8 ms and ∆t = 1 ms in the TSFAsolution (Figure 5.5), resulting in a 2.05% error difference when compared with theEMTP (Figure 5.5(a)). When ∆t = 8 ms, the fault is cleared at 8.63 cycles (143.9 ms)and when ∆t = 1 ms, the fault is cleared 1ms earlier at 8.57 cycles (142.9 ms).1 1.05 1.1 1.15 1.2−20−15−10−5051015202530Current (per-unit of system base)Fault CurrentEMTPΔt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE(Reconstructed) Δt=0.1msTPhasor* Δt=8msTSFA (Envelope)TPhasor* TSFA (Reconstructed) (19.4 pu)EMTP (17.9 pu)Due to largeΔt in TSFAand TPhasor* Due to difference in fault clearing times 1 1.05 1.1 1.15 1.2−20−15−10−5051015202530Time (s)Current(per-unit of system base) EMTP Δt=100µsTSFA-BE (Envelope) Δt=1msTSFA-BE (Reconstructed) Δt=0.1msTPhasor* Δt=8msTSFA (Envelope)TPhasor* TSFA (Reconstructed) (18.9 pu)Due to difference in fault clearing times Large Δt in TPhasor*EMTP (17.9 pu)Figure 5.5: The fault current during a 3-phase fault at the terminals of Generator 1 (Top:using ∆t = 8 ms for TSFA, Bottom: using ∆t = 1 ms for TSFA).1.02 1.03 1.04 1.05 1.061416182022Time (s)Fault Current (per-unit) EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTPhasor* Δt=8ms1.02 1.03 1.04 1.05 1.061416182022Time (s)Fault CurrentEMTP Δt=100µsTSFA-BE (Envelope) Δt=1msTSFA-BE (Reconstructed) Δt=0.1msTPhasor* Δt=8msTSFA (Envelope)TSFA (Reconstructed)(19.4pu) EMTP(17.99pu)TSFA (Reconstructed)(18.9pu) EMTP(17.99pu)TPhasor* TPhasor*TSFA (Envelope)(a) The fault current during a 3-phase fault at the terminals of Generator 1: (Left: TSFA at∆t = 8 ms and Right: TSFA at ∆t = 1 ms).925.2. A Three-Phase Fault at the Generator’s Terminals5.2.1 A Comparison of Using the Actual Frequency with theApproximate FrequencyThe SFA solution using the approximate frequency TSFA* is compared with TSFAto observe the differences between using the actual velocity ω as opposed to theapproximate velocity ωo in the electromechanical equation. As the EMTP solutionuses the approximate velocity, the TSFA* solution is closer to the EMTP reference.The phasor solution is also simulated using the corrected equation in TPhasor andcompared with using the approximate velocity in TPhasor*. The results are shownin Table 5.4), and Figure 5.6.Table 5.4A comparison of the error in frequency between TSAT, TSAT*, TPhasorand TPhasor* with respect to the EMTP referenceTSFA TSFA* TPhasor TPhasor*frequency error (first peak) (%) 0.114 0.0978 0.114 0.147frequency difference (Hz) 0.07 0.06 0.07 0.091.4 1.6 1.8 2 2.2 2.460.86161.261.4Frequency (Hz)Changing the velocity ω in the electromechanical equationcomparing SFA and the phasor solution with the EMTPEMTP Δt=100µsTSFA*-BE Δt=8msTSFA-BE Δt=8ms1.4 1.6 1.8 2 2.2 2.460.86161.261.4Time (s)Frequency (Hz)EMTP Δt=100µsTPhasor* Δt=8msTPhasor Δt=8msTime (s)TPhasor (61.45Hz)TPhasor*(61.47Hz)EMTP (61.38Hz)EMTP (61.38Hz)TSFA*(61.32Hz)TSFA (61.31Hz)Figure 5.6: A comparison of using the grid frequency, ω, with the approximate frequency, ωo,in the SFA solution and in the phasor solution compared to the EMTP for the three-phasefault on the SGIB.935.2. A Three-Phase Fault at the Generator’s Terminals5.2.2 A Comparison of the Discretization Methods in the TSFASolution During the Three-Phase Fault.The trapezoidal and backward Euler discretization methods were compared for thesame three-phase fault condition in Figure 5.1. Figure 5.7 shows that the errordifference between TSFA-Trap and the EMTP reference is 0.26%, while the errordifference between TSFA-BE and the EMTP reference is 0.11% on the first swingof the frequency, making the backward Euler rule more accurate in this scenario.In the electrical power and fault voltage, the trapezoidal rule has oscillations post-fault. To show the oscillations in the fault voltage, the time-step used for bothtrapezoidal and backward Euler in the fault voltage is ∆t = 20 ms. The oscillationsare numerical as they occur peak-peak ∆t apart, which matches the behaviour ofthe trapezoidal method during discontinuities (Chapter 3), in which sustained nu-merical oscillations occur during a step discontinuity, which, like the fault clearing,is a topological change. If the time-step for the trapezoidal rule was smaller, theoscillations would be smaller.Therefore, if the trapezoidal rule was used in TSFA, two approaches to eliminatethese oscillations could be to either to use a multi-rate technique [51] to change toa smaller time-step during the fault and back to the larger time-step after the fault,or to use a multi-rule technique [51], using backward Euler rule during the faultand one time-step after the fault, and then switch back to using the trapezoidal rulefor the post-fault condition (similar to the critical damping adjustment used in theEMTP [38]).945.2. A Three-Phase Fault at the Generator’s Terminals1 1.2 1.4 1.6 1.8 2−200204060Rotor Angle (Deg)Comparison of the Trapezoidal and Backward Euler RulesDuring a Three-Phase Fault EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8ms 1 1.2 1.4 1.6 1.8 25959.56060.56161.5Time (s)Frequency (Hz)Time (s)EMTPTSFA-BETSFA-TrapEMTP 61.38Hz TSFA-BE 61.31HzTSFA-Trap 61.22Hz1 1.2 1.4 1.6 1.8 2−200204060Electrical Power (MW)EMTP Δt=100µsTSFA-BE Δt=8ms (for power), Δt=20ms (for fault voltage) TSFA-Trap Δt=8ms (for power), Δt=20ms (for voltage)1 1.2 1.4 1.6 1.8 2−1−0.500.51Time (s)Fault Voltage (per-nit)Time (s)EMTPTSFA-BETSFA-Trap(*)(* oscillations in trapezoidal are 40ms peak-peak )Figure 5.7: A comparison of the trapezoidal and backward Euler discretization rules onthe rotor angle, frequency, electrical power, and fault voltage for the three-phase faultcondition. Top left: Generator’s rotor angle, Top right: Generator’s frequency, Bottom left:Generator’s electrical power, Bottom right: Fault voltage.95Chapter 6Three-Bus Network Case StudyThe three-bus network is shown in Figure 6.1. The three-bus network topology is atraditional case study studied in transient stability analysis literature [45]. The net-work consists of two generators, a 25 kV, 300 MVA hydro generator and a 13.8 kV,50 MVA hydro generator, delivering a total of 315 MW over a balanced 3-phase250 kV system to a constant impedance load of 310 MW. The 300 MVA generatorhas a kinetic energy of 6 MJ per MVA (rated inertia of H=6 s) at rated velocityand the 50 MVA generator has a kinetic energy of 2 MJ per MVA (rated inertiaof H=2 s) at rated velocity. The direct-axis transient reactance of both generatorsis x′d=0.30 p.u. on the generator base; and the voltage behind the reactance ofGenerator 1 and Generator 2 are 1.18 and 1.48 p.u. on the generator bases. Thetransmission line is modelled using the pi model. The system is a single-line equiv-alent of a balanced, symmetrical, three-phase network and the system’s parametersare given in Table 6.1 below.The simulations were performed under two conditions: a three-phase fault on theterminals of Transformer 2 and when Generator 1 is disconnected from the network. 270MWeP Gen.1 300MVA 15%tX Ideal 25kV:250kV ' 0.3 pudX Ideal 310MW150MVArLLPQ 12%tX 250kV:13.8kV 45MWeP Gen.2 50MVA ' 0.3 pudX Trans. 1 Trans. 2 B A C D E F G Figure 6.1: The one-line diagram of the three-bus network, where two generators feed aconstant impedance load.96Chapter 6. Three-Bus Network Case StudyTable 6.1Network parameters for the three-bus network (Figure 6.1)GeneratorsPowerRating(MVA)VoltageRating(kV)TransientReactancex′d (p.u.1)InertiaConstant (H)(sec)Moment ofInertia (I)(kg·m2)DampingCoefficient(K)(p.u.2)Generator 1 300 25 0.3 6 25330.3 1.5Generator 2 50 13.8 0.3 2 1407.23 1.0GeneratorInitialConditionsReal Power(MW)Reactive Power(MVAR)Terminal VoltageVt (p.u.1∠◦)Eo∠δo(kV6∠◦)Generator 1 270.1 106.58 1∠0◦ 23.25∠13.72◦Generator 2 45 17.40 0.95∠-7.35◦ 12.36∠7.66◦TransformersRated Power(MVA)Voltage Turns Ratio(kV:kV)Reactance Xt(p.u.4 )Trans. 1 600 25:250 0.15Trans. 2 300 250:13.8 0.12LoadsReal Power(MW)Reactive Power(MVAR)Resistance5(Ω)Reactance5(Ω)Load 1 310 150 164.45 339.86TransmissionLinesLineLength(km)R(Ω/km)X(Ω/km)B8 (S/km)R(p.u.3)X(p.u.3)B8(p.u.3)Line 1-2 1800.04(7.2 Ω)70.4(72 Ω)74.3µS,(774µS)70.0115 0.115 0.4836Line 1-3 1500.0267(4 Ω)70.267(40 Ω)74.3µS,(645µS)70.0064 0.064 0.4032Line 2-3 800.04(3.2 Ω)70.4(32 Ω)74.3µS,(344µS)70.00512 0.0512 0.2151 per-unit of generator base 2 per-unit of torque base; 3 per-unit of system base; 4 per-unitof transformer base; 5 Vload = 225.79 kV6 kV, line-neutral peak value [19]; 7 in total Ω or S8 notB2976.1. Critical Clearing Times6.1 Critical Clearing TimesThree-phase balanced faults were applied at different locations in the network toobserve the different CCTs given by the EMTP, TSAT, the TPhasor method, whichuses the correct velocity in the electromechanical equation, TSFA and TSFA*, whichuses the approximated velocity. The comparisons are presented in Table 6.2. Theresults show that the only time the TSFA solution does not match the EMTP so-lution is when the fault is at the terminals of Generator 1. In this scenario, theTSFA solution is 8 ms smaller than the EMTP solution due to two reasons. First,if the EMTP does not check for zero-crossing, then the EMTP solution gives aCCT=284 ms (17 cycles). Second, when the TSFA solution has a reduced time-step to ∆t = 1 ms, the TSFA solution also gives a CCT=284 ms. When these twoconditions happen, the EMTP and the TSFA solutions will match. Therefore, thediscrepancy in the CCT in this case is not due to the SFA solution method becausewith the other fault cases, the TSFA and the EMTP solutions are either identical oronly 1 ms apart. In the case with the phasor solution, TSAT, there is approximatelya 20 ms difference between TSAT and the EMTP, which corresponds to 1 cycle. Thissignifies that the SFA and the EMTP solutions are closer than the TSAT and EMTPsolution.Observing the differences between the CCTs given in the TSAT/TPhasor methodsand the TSFA/TSFA* methods, there is approximately a consistent 10 ms differencebetween the TPhasor and TSAT solutions. In the case when the fault is placed atthe terminals of Transformer 2, in the middle of Line 1-2, and in the middle ofLine 1-3, the CCT given by the TSFA* solution is TSFA* is 1 cycle lower thanwhen using TSFA. These differences indicate that using the correct frequency inthe electromechanical equation can make a slight difference in the CCT, even in thetraditional phasor solution method. The critical clearing angles is the rotor angle ofGenerator 2 and they are all within 10-15◦ apart from each other. For the EMTPsolution, the CCT is calculated with the zero-crossing of the fault current beingdetected, but for the TSFA solution, the CCT is calculated without the zero-crossingdetection because if the zero-crossing detection was implemented to calculate theclearing times, there would be approximately a 3-cycle difference between the stableand unstable clearing time. This is due to the large time-step used in the solutionmethod. A more accurate method to capture the clearing time in TSFA using the986.1. Critical Clearing Timeszero-crossing detection would be to use multi-rate to change the time-step to asmaller one at the time of the fault clearing such that there would not be this largedifference. This improvement is something to be done in future research.Table 6.2The CCTs for three-phase faults at different locations of the three-bussystem in Figure 6.1 in 60-Hz cycles and the critical clearing angle indegrees.ContingencyEMTP(∆t = 100µs)TSAT(∆t = 8ms)TPhasor(∆t = 8ms)TSFA(∆t = 8ms)TSFA*(∆t = 8ms)Fault on theterminals ofGenerator 1(Location A)17(288 ms)99.5◦16(267 ms)90.7◦16(272 ms)103.8◦17(280 ms)101.0◦17(280 ms)104.5◦Fault on theterminals ofGenerator 2(Location G)15(247 ms)87.5◦13(225 ms)83.2◦14(240 ms)100.6◦15(248 ms)96.7◦15(248 ms)98.9◦Fault at theload (Bus 3)17(296 ms)101.4◦16(277 ms)97.6◦17(288 ms)112.0◦17(296 ms)107.8◦17(296 ms)110.3◦Fault on theterminals ofTransformer 1(Location B)16(272 ms)91.9◦15(254 ms)87.0◦15(264 ms)103.4◦16(272 ms)100.1◦16(272 ms)102.7◦Fault on theterminals ofTransformer 2(Location F)15(255 ms)91.5◦14(233 ms)85.4◦14(240 ms)100.4◦15(256 ms)100.1◦14(248 ms)98.94◦Fault in themiddle ofLine 1-2(Location C)18(313 ms)104.2◦17(287 ms)96.1◦17(296 ms)108.0◦18(312 ms)110.0◦18(304 ms)107.0◦Fault in themiddle ofLine 1-3(Location D)18(304 ms)104.0◦16(274 ms)92.2◦17(288 ms)110.1◦18(304 ms)110.0◦17(296 ms)107.8◦Fault in themiddle ofLine 2-3(Location E)16(280 ms)97.8◦15(260 ms)92.4◦15(272 ms)106.1◦16(280 ms)105.0◦16(280 ms)107.9◦996.2. A Three-Phase Fault at the Terminals of Transformer Two6.2 A Three-Phase Fault at the Terminals ofTransformer TwoA three-phase fault is applied to the terminals of Transformer 2 at t = 1 s and iscleared in 8 cycles by opening circuit breaker B (Figure 6.1). Figures 6.2 and 6.3present the electromechanical variables: rotor angle, frequency and electrical power.Figures 6.4, 6.5 and 6.6 present the fault voltage and fault currents, where in Figure6.6, the fault current is simulated using a ∆t = 1 ms. The TPhasor* solution is usedinstead of TSAT for the fault current in Figure 6.5 and Figure 6.6 because TSATdoes not explicitly give the fault current.The fault occurs at t = 1 s, and both generators begin to accelerate. As Generator 2has a lower inertia, its frequency deviations are larger than Generator 1’s deviations.Generator 2 deviates 2 Hz from the 60-Hz frequency while Generator 1 only deviates0.7 Hz from the 60-Hz frequency. Once the fault is cleared, the frequency begins todecrease, but due to the extra kinetic energy stored in the generators, the anglescontinues to increase for another 2.7 cycles (45.4 ms). As the fault is cleared beforethe CCT of 14 cycles, the rotor angle is considered “first-swing stable.”Both the TSFA and EMTP solutions use zero-crossing detection for the fault currentand so the clearing times differ among the solution methods. The EMTP solutionclears the fault after 7.9 cycles (132 ms) and the TSFA solution clears the faultafter 8.63 cycles (144 ms). The TSAT solution clears the fault at exactly 8 cycles(133 ms). Frequency errors are observed for the first swing in Generator 1 and 2when the results are compared among the benchmark, EMTP, with the TSFA andTSAT solutions. For Generator 1, the frequency error was 0.116% between TSATand EMTP, and 0.033% between TSFA and EMTP (Figure 6.2). For Generator2, the frequency error was 0.355% between TSAT and EMTP and 0.161% betweenTSFA and EMTP (Figure 6.3). This is shown in Table 6.3.Table 6.3A comparison of the error in frequency between TSAT and TSFA withrespect to the EMTP reference solution.TSAT (∆t =8 ms) TSFA (∆t =8 ms)Gen.1 frequency error (%) 0.116 0.033Gen.2 frequency error (%) 0.355 0.1611006.2. A Three-Phase Fault at the Terminals of Transformer Two0.5 1 1.5 2 2.5 3−40−2002040Rotor Angle Difference (Deg) Rotor Angle Difference between Gen.1 and Gen.2EMTP Δt=100µsTSFA-BE Δt=8msTSAT Δt=8msSee Fig.6.2(a)Fault at t=1sFault cleared after 8 cycles0TSFAEMTPTSAT0.5 1 1.5 2 2.5 36060.260.460.660.8Frequency (Hz)Frequency of Gen.1EMTP Δt=100µsTSFA-BE Δt=8msTSAT Δt=8msSee Fig.6.2(b)Fault at t=1sFault cleared after 8 cyclesEMTPTSATTSFA0.5 1 1.5 2 2.5 3200240280320360400Time (s)Electrical Power (MW) Electrical Power of Gen.1 EMTP Δt=100µsTSFA-BE Δt=8msTSAT Δt=8msSee Fig.6.2(c)Fault at t=1sFault cleared after 8 cyclesTSFAEMTPTSATFigure 6.2: A comparison of the EMTP, TSFA, and TSAT solutions for the electromechan-ical variables of Generator 1 in the system of Figure 6.1 during a three-phase fault at theterminals of Transformer 2.1.6 1.7−100102030Time (s)Rotor Angle Difference (Deg)EMTP TSFA-BE TSAT EMTPTSFATSATFig.6.2 (a) Gen.2-Gen.1 angle.1.6 1.7 1.860.560.5560.660.6560.7Time (s)Frequency (Hz)EMTPTSATTSFA-BEEMTP 60.62HzTSFA 60.64HzTSAT 60.69HzFig.6.2 (b) Gen.1’s frequency.1.7 1.8 1.9240280320Time (s)Electrical Power (MW)EMTPTSFA-BETSATEMTPTSFATSATFig.6.2 (c) Gen.1’s power.1016.2. A Three-Phase Fault at the Terminals of Transformer Two0.5 1 1.5 2 2.5 359606162Frequency (Hz)Frequency of Gen.2EMTP Δt=100µsTSFA-BE Δt=8msTSAT Δt=8ms See Fig.6.3(a)Fault at t=1sFault cleared after 8 cycles EMTPTSATTSFA0.5 1 1.5 2 2.5 3−4004080120160Time (s)Electrical Power (MW) Electrical Power of Gen.2 EMTP Δt=100µsTSFA-BE Δt=8msTSAT Δt=8msSee Fig.6.3(b)Fault at t=1sFault cleared after 8 cyclesTSFAEMTPTSATFigure 6.3: A comparison of the EMTP, TSFA, and TSAT solutions for the electromechan-ical variables of Generator 2 in the system of Figure 6.1 during a three-phase fault at theterminals of Transformer 2.1.56161.562Time (s)Frequency (Hz)EMTP TSFA-BETSAT EMTP 61.98 HzTSFA 62.08HzTSAT 62.2HzFig.6.3 (a) Gen.2’s frequency.1.5 1.6 1.7−4004080120Time (s)Electrical Power (MW)EMTPTSFA-BETSATEMTPTSFATSATFig.6.3 (b) Gen.2’s electrical power.1026.2. A Three-Phase Fault at the Terminals of Transformer Two0.9 1 1.1 1.2 1.3 1.4−0.8−0.6−0.4−0.200.20.40.60.8Time (s)Voltage(per-unit of system base)Voltage at the faulted bus (Terminals of Transformer 2)EMTP Δt=100µsTSFA-BE(Envelope) Δt=8msTSFA-BE(Reconstructed) Δt=0.1msTSAT Δt=8msDue to difference in clearing timesEMTPTSFA (Reconstructed)TSATTSFA(Envelope)Figure 6.4: The voltage at the faulted bus during a 3-phase fault at the terminals of Trans-former 2.1 1.05 1.1 1.15 1.2−10−5051015Time (s)Current (per-unit of system base)Fault CurrentEMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTPhasor* Δt=8msEMTP (12.25pu)TSFA (Reconstructed) (11.01pu) TPhasor*TSFA (Envelope)Due to difference in clearing timesDue to large Δt in TSFAand in TPhasor*Figure 6.5: The fault current during a 3-phase fault at the terminals of Transformer 2(using a ∆t = 8 ms for the SFA solution).1 1.05 1.1 1.15 1.2−10−5051015Time (s)Current (per-unit of system base)Fault CurrentEMTP Δt=100µsTSFA-BE (Envelope) Δt=1msTSFA-BE (Reconstructed) Δt=0.1msTPhasor* Δt=8msEMTP (12.25pu)TSFA (Reconstructed) (11.17pu)TPhasor*TSFA (Envelope) Due to difference in clearing timesDue to large Δt in TPhasor*Figure 6.6: The fault current during a 3-phase fault at the terminals of Transformer 2(using a ∆t = 1 ms for the SFA solution).1036.2. A Three-Phase Fault at the Terminals of Transformer TwoAs in the SGIB case study, the fault voltage using the TSFA solution interpolatedback into the physical domain is identical to the EMTP solution and the TSFAenvelope is the same as the phasor solution. For the fault current, using a smallertime-step in TSFA makes a difference as shown in Figures 6.5 and 6.6. In Figure6.6, when the ∆t for the TSFA solution is reduced to 1 ms from 8 ms, the differencein fault current between the EMTP and the TSFA solutions at the second swingpost-fault reduces by 1.31%.6.2.1 A Comparison of Using the Actual Velocity with theApproximate VelocityThe SFA solution using the approximate velocity TSFA* is compared with TSFAto observe the differences between using the actual velocity ω as opposed to theapproximate velocity ωo in the electromechanical equation. This comparison isshown in Figure 6.7. Between the TSFA and TSFA* solution there is only a 0.01 Hzdifference in the frequency of Generator 1, but a 0.04 Hz difference in the frequencyof Generator 2. This difference indicates that using the correct network frequencyhas an effect on the simulation results. These results are tabulated in Table 6.4 .1.2 1.4 1.6 1.8 260.560.5560.660.65Frequency of Gen.1 (Hz)Comparison of the SFA Solution changing ω with respect to the EMTPEMTP Δt=100µsTSFA-BE Δt=8msTSFA*-BE Δt=8ms1 1.2 1.4 1.6 1.8 259.56060.56161.562Time (s)Frequency of Gen.2 (Hz)Time (s)TSFA(60.64Hz)TSFA*(60.65Hz)EMTP (60.62Hz) EMTP (61.9Hz)TSFA(62.06Hz)TSFA*(62.1Hz)Figure 6.7: A comparison of using the grid frequency, ω, with the approximate frequency,ωo, in the SFA solution in the electromechanical equation for the three-phase fault.1046.2. A Three-Phase Fault at the Terminals of Transformer TwoTable 6.4A comparison of the error in frequency between the different TSFAsolutions with respect to the EMTP reference solution.TSFA TSFA*Gen.1 frequency error (%) 0.033 0.049Gen.2 frequency error (%) 0.259 0.323The phasor solution method using the approximate frequency, TPhasor*, is alsocompared with the phasor solution using the correct frequency, TPhasor. In thecase of the phasor solution, to compare with the EMTP benchmark, the internalgenerator voltage (Assumption 5 of Section 4.3) was also implemented in the phasorsolution. The four conditions for the phasor solution are tabulated in Table 6.5 andshown in Figure 6.8. For these comparisons, TPhasor* was used instead of TSATto ensure consistency among the solution methods.Table 6.5The four different conditions for the phasor solution in Figure 6.8.Omega(Constant)Omega(Updated)GeneratorVoltage(Constant)GeneratorVoltage(Updated)TPhasor* 3 7 3 7TPhasor 7 3 3 7TPhasor** 3 7 7 3TPhasor*** 7 3 7 3Similar to the comparison between TSFA and TSFA*, between TPhasor and TPha-sor*, there is a 0.01 Hz difference in the frequency of Generator 1 but a 0.04 Hzdifference in the frequency of Generator 2. On the other hand, the TPhasor**and TPhasor*** are closer to the EMTP on the second swing than the TPhasoror TPhasor* solutions, because correcting the generator voltage allows the phasorsolution to settle. In the case when the generator voltage is internally updated usingthe grid frequency ω, there is a 0.57 Hz difference when using ωo in TPhasor** andω in TPhasor***, which indicates that using the correct grid frequency does make aslight difference in the simulation results. Although there is only a slight differencewhen using the corrected frequency, this does not mean that the correction will notmake a difference in a situation where the frequency of the generator deviates morethan 2 Hz from the 60-Hz frequency. These results are tabulated in Table 6.6.1056.2. A Three-Phase Fault at the Terminals of Transformer Two1.2 1.4 1.6 1.8 260.560.660.760.8Frequency of Gen.1 (Hz)Comparison of the Phasor Solutions changing ω with respect to the EMTPEMTP Δt=8msTPhasor* Δt=8msTPhasor Δt=8msTPhasor*** Δt=8msTPhasor** Δt=8ms1 1.2 1.4 1.6 1.86161.56262.5Time (s)Frequency of Gen.2 (Hz)Time (s)TPhasor* (60.75Hz)TPhasor(60.74Hz)EMTP (62.62Hz)TPhasor***(62.20Hz)TPhasor**(62.24Hz)EMTP (62.0Hz)TPhasor* (62.3Hz)TPhasor(62.26Hz)TPhasor**(62.69Hz)TPhasor***(62.68Hz)Figure 6.8: A comparison of using the grid frequency, ω, with the approximate frequency,ωo, in the phasor solution for the three-phase fault.Table 6.6A comparison of the error in frequency between the different TPhasorsolutions with respect to the EMTP reference solution.TPhasor TPhasor* TPhasor** TPhasor***Gen.1 frequency error (%) 3.002 2.986 0.112 0.096Gen.2 frequency error (%) 0.419 0.484 0.387 0.3236.2.2 A Comparison of the Discretization Rules in the TSFASolution During the Three-Phase Fault.There is negligible difference between using the trapezoidal and backward Eulerdiscretization rules in the rotor angle difference between the two generators andin the generators’ frequencies; however, in the electrical power and post-fault faultvoltage, the solution using the trapezoidal rule has sustained numerical oscillations.These oscillations are different from the EMTP oscillations during the fault. Inthe EMTP, the oscillations are due to the instantaneous power captured instead ofthe average power. When using a smaller time-step for the trapezoidal rule, suchas in the electrical power of Generator 2 in Figure 6.9, the numerical oscillationsbecome smoother and the trapezoidal solution is closer to the EMTP solution thanthe backward Euler solution (by 0.1 Hz on the first-swing of Generator 1).1066.2. A Three-Phase Fault at the Terminals of Transformer Two1 1.2 1.4 1.6 1.8 20102030Rotor Angle Difference (Deg)EMTP Δt=100µsTSFA-BE Δt=8msTSFA-TrapΔt=8msEMTPTSFA-BETSFA-Trap1 1.2 1.4 1.6 1.8−0.500.511.5Time (s)Fault Voltage (per-unit) EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msEMTPTSFA-BETSFA-TrapTime (s)Comparison of the Trapezoidal and Backward Euler Rules During a Three-Phase Fault1.2 1.4 1.6 1.8 260.360.460.560.660.7Frequency of Gen.1 (Hz) EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msEMTPTSFA-BETSFA-Trap1.4 1.6 1.8 2 2.260.56161.56262.5Time (s)Frequency of Gen.2 (Hz) EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msEMTPSFA-BESFA-TrapTime (s)1.1 1.2 1.3 1.4 1.5 1.6 1.7100150200250300350400Electrical Power of Gen 1 (MW)EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msEMTPTSFA-BETSFA-Trap1.6 1.7 1.8050100150Time (s)Electrical Power of Gen.2 (MW)EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap(a) Δt=8msTSFA-Trap(b) Δt=1msEMTPTSFA-BETSFA-Trap(b)Time (s)TSFA-Trap(a)Figure 6.9: A comparison of the trapezoidal and backward Euler discretization rules inSFA on the rotor angle, frequency, electrical power, and fault voltage during the three-phasefault at the terminals of Transformer 2.1076.2. A Three-Phase Fault at the Terminals of Transformer Two6.2.3 The DC Offset Captured in the Fault Current in SFAThe solution of a first-order differential equation consists of a transient solution anda steady-state solution. Because the SFA and the EMTP solutions discretize thecontinuous-time differential equations, they can capture the DC offset in the faultcurrent that happens in the network due to the transient part of the differentialequation. The amount of DC offset in the fault current depends on when the faultoccurs. If the fault occurs when the fault voltage crosses zero, the DC offset is at itsmaximum value. This offset is not captured in the traditional phasor solution. Tocapture the maximum DC offset, the time-step used in SFA must be around 1.67 msdue to the time-frequency relationship (Section 2.1.2). The 0 Hz (DC) frequencyfrom the physical domain is a -60-Hz frequency in the shifted domain and the time-step needed to capture it is ∆t = 110×60 = 1.67 ms. To model the DC offset, athree-phase fault was applied to the middle of Line 1-2 when the voltage crossedzero. Both the trapezoidal and the backward Euler rules were used in SFA andcompared with the reference EMTP solution (Figure 6.10). The DC componentis accurately captured when using the trapezoidal rule, but not when using thebackward Euler rule. The error from using the backward Euler rule can be explainedin the error analysis of Chapter 3, where it is observed that at ∆t = 1 ms, thebackward Euler has a large magnitude error (HDTHCT ≈ 5) and a large phase distortion(≈ 80◦) at frequencies between 0-0.5 Hz (Figure 3.5), whereas at ∆t = 1 ms, thesolution using trapezoidal has a smaller magnitude error (HDTHCT ≈ 0.4→ 2.5) and nophase distortion (Figure 3.2).1.05 1.1 1.15 1.2 1.25 1.3−10−505101520Time (s)Current(per-unit of system base)Fault Current with DC ComponentEMTP Δt=100µsSFA-Trap Δt=1msSFA-BE Δt=1msTPhasor* Δt=8msSFA-Trap (14.17pu)EMTP (14.07pu)TPhasor*(11.34pu)SFA-BE(10.98pu)Figure 6.10: A three-phase fault is cleared when the voltage crosses zero to capture themaximum DC offset. The SFA solution, using both the trapezoidal and backward Eulerrules, and TPhasor* solution are compared to the EMTP.1086.3. Generator One Disconnected From the Network6.3 Generator One Disconnected From the NetworkThe second scenario simulated is the situation when Generator 1 is suddenly dis-connected from the system at t = 1 s. In this scenario, because the load remainsconstant, Generator 2 is not capable of supplying all of the load by itself and becausegovernor control action is not considered, the energy that supplies the load comesfrom the kinetic energy of Generator 2. Therefore, Generator 2’s frequency willdecrease and the network frequency will change drastically from 60 Hz. As shownin Figures 6.11 and6.12 the system is very sensitive to using the actual value of ωinstead of the constant ωo. The TSFA solution with the corrected frequency dropsto 54.93 Hz whereas the EMTP solution settles at 55.4 Hz. In this case, the TSFAsolution is actually more accurate than the EMTP solution because the EMTP doesnot use the actual velocity ω when converting between torque and power in the elec-tromechanical equations; hence, TSFA solution has a 0.47 Hz difference comparedto the EMTP, whereas the TSFA* solution has a 0.2 Hz difference.When the trapezoidal rule used, there are strong numerical oscillations are givenby the trapezoidal rule. It is interesting to note that the trapezoidal in SFA hasoscillations whilst the EMTP solution, which uses the trapezoidal rule with CDA,does not have these oscillations.TPhasor* was used instead of TSAT because TSAT failed after 2 ms of simulationtime. Because the TPhasor* solution does not use the internal voltage controlmechanism, the frequency continues to decrease without stabilizing. Figure 6.13shows the difference in the phasor solution with the different frequency changingcombinations of Table 6.5 above. As predicted, when the generator’s voltage isupdated with the corrected frequency, the TPhasor solution stabilizes; however, thesolution stabilizes to a higher frequency than with the EMTP and SFA.1096.3. Generator One Disconnected From the Network0 1 2 3 4 5 6 7 851525354555657585960Frequency (Hz)Frequency of Gen.2 After Gen.1 Disconnects EMTP Δt=100µsTSFA-BE Δt=8msTSFA*-BE Δt=8msTPhasor* Δt=8msGenerator 1 disconnectedfrom the network See 6.12EMTP(55.4Hz)TSFA*(55.2Hz)TSFA (54.93Hz)TPhasor*(52.1Hz)0 1 2 3 4 5 6 7 851525354555657585960Frequency (Hz)Frequency of Generator 2 after Generator 1 Disconnects EMTP Δt=100µsTSFA- BE Δt=8msTSFA* Δt=8msTSFA - Trap Δt=8msTPhasor* Δt=8msGenerator 1 disconnectedfrom the networkSee 6.12TSFA-TrapEMTP(55.4Hz)TSFA*(55.2Hz)TSFA (54.93Hz)TPhasor*(52.1Hz)Figure 6.11: The frequency of Generator 2 after disconnecting Generator 1 from the three-bus network in Figure 6.1 using the backward Euler rule (top) and both rules (bottom).5.8 5.9 6 6.15253545556Frequency (Hz)EMTP Δt=100µsTSFA-BE Δt=8msTSFA*-BE Δt=8msTPhasor* Δt=8msEMTPTSFA-BETSFA-TrapTPhasor*5.8 5.9 6 6.15253545556EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msTPhasor* Δt=8msTPhasor*TSFA*-BEEMTPTSFA-BE Time (s) Time (s) Figure 6.12: The insert of Figure 6.11. Left: without the trapezoidal rule, Right: with thetrapezoidal rule.0 1 2 3 4 5 6 7 85051525354555657585960Time (s)Frequency (Hz)The phasor solution with the corrected electromechanicalequation when Generator 1 disconnects.EMTP Δt=100µsTPhasor Δt=8msTPhasor* Δt=8msTPhasor** Δt=8msTPhasor***Δt=8msSFA Δt=8msSFA* Δt=8msTPhasor** (56.7Hz)TPhasor*** (56.63Hz)EMTP (55.3 Hz)TPhasor*TPhasorTSFA* (55.1 Hz)TSFA (54.93 Hz)Time (s)Figure 6.13: A comparison of the traditional phasor solution using the corrected electrome-chanical equation.110Chapter 7Thirty-Nine Bus Network CaseStudyThe IEEE 39-bus test system [1] was studied to demonstrate the behaviour of theSFA solution method in a large network. This system is shown in Figure 7.1. Thedata for this system is taken from [1], with the modifications provided in Tables A1A1 A3 and A4 of the Appendix. Figure 7.1: The IEEE 39-bus test system.1117.1. Critical Clearing TimesIn this system, ten generators of varying MVA ratings feed eighteen loads. Gener-ator 5 represents a lower-inertial machine as it is rated at 50 MVA, with an inertiaconstant H = 1 s. Generator 1 is the reference generator, rated at 1100 MVA, withan inertia constant of H = 5 s. To model the power imbalance, two scenarios areapplied to the network. In the first scenario, a three-phase balanced fault occurson Bus 16. In the second scenario, Generator 9 is disconnected from the grid. Inaddition, the CCTs of three-phase faults occurring on different locations along the39-bus network are given in Table 7.1.7.1 Critical Clearing TimesThe CCTs for 3-phase balanced faults at different locations along the 39-bus net-work are given in Table 7.1. The CCTs are in terms of 60-Hz cycles, (rounded tothe lowest cycle) and in milliseconds. In addition, the critical clearing angles ofthe generators that are most affected by the fault are recorded. The results showthat the generators with the lowest moment of inertia, Generators 5, 2, 7, and 8,(Table A1) are also the ones that are the most affected when the fault does notoccur on the terminals of one of the generators.The TSFA and the EMTP solution methods give identical CCTs (in 60-Hz cycles)for almost all the cases that were simulated. On average, the difference in millisec-onds in CCTs between the TSFA and the EMTP range by approximately 10 ms(Table 7.1). The difference in milliseconds between the two solutions could be dueto the difference in time-steps that were used for each solution method. In addition,the EMTP solution checks for the zero-crossing of the fault current to calculate theCCT, whereas the TSFA does not (due to its very large time-steps), an additionalprocedure which could account for the small differences. Likewise, the critical clear-ing angle between the two solution methods are very similar, sometimes with only1◦ difference. Within the SFA solution, the CCTs show that there are differencesbetween the TSFA and TSFA* solutions observed on Bus 14 and Bus 16 (both busesare near the middle of the network), which indicates that using the corrected ve-locity in the electromechanical equation can make a difference. Interestingly, whenthe fault occurs on Bus 7, all the solution methods (TSAT, TSFA, TSFA* and theEMTP) give the identical CCT of 20 cycles.1127.1. Critical Clearing TimesTable 7.1The CCTs for three-phase faults at different locations of the three-bussystem in Figure 7.1 in 60-Hz cycles and the critical clearing angle indegrees.Contingency EMTP TSAT TSFA TSFA*(∆t = 100µs) (∆t = 8ms) (∆t = 8ms) (∆t = 8ms)Fault on Bus 16 17 14 18 17(Gen 5’s angle) (285 ms) (238 ms) (304 ms) (296 ms)171◦ 60.2◦ 107◦ 111◦Fault on Bus 31 15 14 14 14(Gen 2’s angle) (260 ms) (244 ms) (240 ms) (240 ms)104◦ 64◦ 101◦ 111◦Fault on Bus 32 14 12 14 14(Gen 3’s angle) (240 ms) (201 ms) (240 ms) (240 ms)108◦ 86◦ 106◦ 108◦Fault on Bus 33 15 12 15 15(Gen 4’s angle) (258 ms) (210 ms) (248 ms) (248 ms)165◦ 91◦ 106◦ 107◦Fault on Bus 34 11 9 11 11(Gen 5’s angle) (195 ms) (148 ms) (184 ms) (184 ms)126◦ 112◦ 147◦ 150◦Fault on Bus 35 22 18 22 22(Gen 6’s angle) (370 ms) (304 ms) (384 ms) (384 ms)110◦ 99◦ 118.4◦ 119◦Fault on Bus 36 10 9 10 10(Gen 7’s angle) (175 ms) (152 ms) (176 ms) (176 ms)96◦ 165◦ 104◦ 105◦Fault on Bus 37 12 10 12 12(Gen 8’s angle) (205 ms) (179 ms) (208 ms) (208 ms)99◦ 93◦ 113◦ 113◦Fault on Bus 38 12 10 12 12(Gen 9’s angle) (205 ms) (170 ms) (208 ms) (208 ms)85◦ 74◦ 91◦ 91.2◦Continued on next page1137.1. Critical Clearing TimesTable 7.1 – Continued from previous pageContingency EMTP TSAT TSFA TSFA*(∆t = 100µs) (∆t = 8ms) (∆t = 8ms) (∆t = 8ms)Fault on Bus 30 40 30 39 39(Gen 10’s angle) (680 ms) (502 ms) (664 ms) (664 ms)169◦ 152◦ 168◦ 171◦Fault on Bus 4 21 19 21 21(Gen 3’s angle) (365 ms) (330 ms) (360 ms) (360 ms)150◦ 132◦ 149◦ 150◦Fault on Bus 7 20 20 20 20(Gen 2’s angle) (348 ms) (334 ms) (336 ms) (336 ms)140◦ 118◦ 146◦ 146◦Fault on Bus 8 20 19 20 20(Gen 2’s angle) (340 ms) (320 ms) (328 ms) (328 ms)138◦ 126◦ 144◦ 146◦Fault on Bus 10 15 13 15 15(Gen 3’s angle) (257 ms) (225 ms) (256 ms) (256 ms)113◦ 83◦ 119◦ 119◦Fault on Bus 14 22 19 22 22(Gen 3’s angle) (370 ms) (317 ms) (368 ms) (360 ms)153◦ 127◦ 155 ◦ 152◦Fault on Bus 19 11 10 11 11(Gen 5’s angle) (194 ms) (175 ms) (192 ms) (192 ms)152◦ 138◦ 159◦ 159◦Fault on Bus 21 18 15 18 18(Gen 7’s angle) (300 ms) (260 ms) (312 ms) (312 ms)155◦ 135◦ 167◦ 167◦Fault on Bus 25 12 11 12 12(Gen 9’s angle) (212 ms) (191 ms) (216 ms) (216 ms)113◦ 122◦ 73◦ 73◦1147.2. A Three-Phase Fault Applied to Bus 167.2 A Three-Phase Fault Applied to Bus 16In the first scenario, a balanced three-phase fault is applied to Bus 16 at t = 1 sand cleared after 8 cycles. Due to the difference in time-steps and the zero-crossingdetection in the TSFA and the EMTP solutions, the fault is cleared at differenttimes. In TSFA, the fault is cleared at 8.15 cycles (135 ms), in TSAT, the fault iscleared at 8 cycles (133 ms), and in the EMTP, the fault is cleared at 8.39 cycles(139 ms). The electromechanical variables of Generator 5 are analysed as Generator5 deviates the most from 60 Hz (±3 Hz from 60 Hz) due to its low moment of inertia(Figure 7.2). Among the three solution methods (TSFA, TSAT and the EMTP),there is negligible error on the first frequency swing, but on the second frequencyswing (Figure 7.2(b)), there is a 0.16% between the TSFA and the EMTP and a0.25% error between the traditional phasor and the EMTP (Table 7.2).In Figure 7.3 the TSFA solution is compared with the TSFA* solution to observe ifchanging the velocity in the electromechanical equation is noticeable. Because theEMTP uses ωo in the electromechanical equation, the TSFA* is closer to the EMTPsolution on the second swing (Table 7.3), which agrees with the 8ms (half-cycle)difference in CCTs between the TSFA and TSFA* when the fault was on Bus 16(Table 7.1).1157.2. A Three-Phase Fault Applied to Bus 160.5 1 1.5 2 2.5 3−40−20020406080100Rotor Angle Difference (Deg) Rotor Angle Difference between Gen.1 and Gen.5 after a 3-Phase Fault on Bus 16EMTP Δt=100µsTSFA-BE Δt=8ms TSAT Δt=8msFault cleared after 8 cyclesFault at t=1sSee Fig.7.2(a)EMTPTSFATSAT0.5 1 1.5 2 2.5 35758596061626364Frequency (Hz)Frequency of Generator 5 after 3-Phase Fault on Bus 16EMTPΔt=100µsTSFA-BE Δt=8ms TSAT Δt=8msFault at t=1sSee Fig.7.2(b)Fault cleared after 8 cyclesEMTP (63.4Hz)TSFA (63.4Hz)TSAT (63.6Hz)EMTP (57.84Hz)TSFA (57.87Hz)TSAT (57.4Hz)0.5 1 1.5 2 2.5 3−80−4004080120160200Time(s)Electrical Power (MW) Electrical Power of Gen.5 after 3-Phase Fault on Bus 16EMTPΔt=100µsTSFA-BEΔt=8ms TSATΔt=8msSee Fig.7.2(c)Fault cleared after 8 cyclesFault at t=1sEMTPTSFATSATFigure 7.2: A comparison of the EMTP, TSFA, and TSAT solutions for the electromechan-ical variables of Generator 5 in the system of Figure 7.1 during a three-phase fault on Bus16.1.3 1.47580859095100Time (s)Rotor Angle Difference (Deg)EMTPTSFATSAT(a) Rotor Angle of Gen.51.4 1.5575859606162Time (s)Frequency (Hz)EMTP(57.35Hz)TSFA(57.37Hz)TSAT (57.0Hz)(b) Frequency of Gen.51.6 1.7−100−50050100150200Time (s)Electrical Power (MW) EMTPTSFATSAT(c) Electrical Power of Gen.51167.2. A Three-Phase Fault Applied to Bus 161 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 25758596061626364Time (s)Frequency (Hz)Changing the velocity ω in the electromechanical equationcomparing SFA with the EMTP for the Fault on Bus 16EMTP Δt=100µsTSFA-BE Δt=8ms TSFA*-BE Δt=8msEMTP (63.83Hz)TSFA (63.78Hz)TSFA*(63.89Hz) EMTP (63.49Hz)TSFA (63.39Hz)TSFA*(63.48Hz)Figure 7.3: A comparison of the SFA solution using the actual frequency ω instead of theapproximated frequency ωo.Table 7.2A comparison of the error in frequency between TSAT and TSFA withrespect to the EMTP reference solution.TSAT TSFAGen.5 frequency error (first swing, Figure (b)) (%) 0.61 0.035Gen.5 frequency error (second top swing) (%) 0.252 0.16Table 7.3A comparison of the error in frequency between the different TSFAsolutions with respect to the EMTP reference solution.TSFA TSFA*Gen.5 frequency error (first swing) (%) 0.078 0.094Gen.5 frequency error (second swing) (%) 0.158 0.016Figures 7.4 and 7.5 (below) show the fault voltage and fault current during thethree-phase fault at Bus 16. The TSFA solution, interpolated back into physicaltime, matches the EMTP references solution almost identically both in the faultvoltage, pre-fault and post-fault, and in the fault current. In the fault voltage, thereis a difference of ≈ 0.04 per-unit between the TSFA solution and the EMTP pre-fault and a difference of ≈ 0.03 per-unit post-fault. In the fault current, there is adifference of ≈ 4 per-unit between the TSFA and the EMTP.1177.2. A Three-Phase Fault Applied to Bus 160.9 1 1.1 1.2 1.3 1.4 1.5−1−0.500.51Time (s)Voltage(per-unit of system base)Voltage at the faulted bus (Bus 16)EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE(Reconstructed) Δt=0.1msTSAT Δt=8msDue to the difference in fault clearing timesFigure 7.4: The voltage at the faulted bus during a 3-phase fault at Bus 16 in Figure 7.1.2.3 2.35 2.4 2.45 2.5Voltage at the faulted bus (Bus 16)EMTP Δt=100µsTSFA-BE (Envelope) Δt=8msTSFA-BE (Reconstructed) Δt=0.1msTSAT Δt=8ms0.9 0.92 0.94 0.96 0.98 1−1−0.8−0.6−0.4−0.200.20.40.60.81Fault Voltage (per-unit) EMTP (0.82pu) TSFA (Reconst.) (0.78pu) EMTP (0.83pu) TSFA (Reconst.) (0.80pu) Time (s) Time (s) (a) A close-up of the voltage at the faulted bus shows that the reconstructed TSFA and the EMTPmatch. (Left: pre-fault, Right: post-fault).1 1.05 1.1 1.15 1.2−100−80−60−40−20020406080100120Time (s)Fault Current (per-unit of system base ) Fault Current EMTP Δt=100µsTSFA (Reconstructed) Δt=0.1msTSFA (Envelope) Δt=8msTPhasor* Δt=8msEMTP (95.5pu)TSFA (Reconstructed) (100.8pu) TPhasor*TSFA (Envelope)Due to the difference in fault clearing timesDue to largeΔt in TSFAand in TPhasor*Figure 7.5: The fault current during a 3-phase fault at Bus 16 in Figure 7.1.1187.3. Generator Nine Disconnected from the Network7.3 Generator Nine Disconnected from the NetworkWhen Generator 9 is disconnected from the network, the frequency of the othergenerators in the network decreases. Using both trapezoidal and backward Eulerin the TSFA solution, Figure 7.6 shows the frequency of Generator 5. There arestrong oscillations with the trapezoidal rule. Between the TSFA and the EMTPsolutions, there is a 0.1 Hz difference in frequency and when comparing the differencebetween using the corrected velocity (ω) with using the approximate one (ωo), thereis negligible difference. One possible reason for the change in frequency not makinga difference is because the frequency drop in this network is small (only 2 Hz). Theseresults differ from the three-bus network (Chapter 6), where there was a frequencydrop of 5 Hz in the network and using the corrected network frequency had an effecton the solution.0 1 2 3 4 5 6 7 857585960Frequency (Hz)Frequency of Gen.5 after Gen.9 is Disconnected EMTP Δt=100µsTSFA-BE Δt=8msTSFA*-BEΔt=8msTSAT Δt=8msGen.9 disconnectedEMTP (58.1Hz)TSFA (58.2Hz)TSAT (56Hz)TSFA* (58.2Hz)0 1 2 3 4 5 6 7 85657585960Time (s)Frequency (Hz)EMTP Δt=100µsTSFA-BE Δt=8msTSFA-Trap Δt=8msTSAT Δt=8msGen.9 disconnectedTSAT (56Hz)EMTP (58.1Hz)TSFA-BE (58.2Hz)TSFA-Trap (Oscillate)Figure 7.6: The frequency of Generator 5 after losing Generator 9 in the 39-bus Network ofFigure 7.1 using the backward Euler rule (top) and both the backward Euler and trapezoidalrules (bottom). Using the trapezoidal rule gives numerical oscillations.119Chapter 8Conclusion and Future Research8.1 ConclusionThe research in this thesis demonstrates the application of Shifted Frequency Anal-ysis for power system transient stability studies. The results show that through theuse of a shifting transformation, SFA is able to trace the frequency dynamics aroundthe 60-Hz grid frequency using large discretization time-steps. This ability allowsthe SFA solution method to be used as an alternative to the traditional solutionmethods in transient stability studies.The research objectives and contributions defined in Chapter 1 are addressed andimplemented in the case studies of Chapters 5, 6, and 7. A generalized SFA solutionalgorithm TSFA was written for a general power system network and applied in thetransient stability case studies of a single-generator infinite-bus network, a three-bus network, and a 39-bus network. The TSFA algorithm implemented a number ofnovel algorithmic schemes. A zero-crossing detection method for the fault currentwas developed, which allows the zero-crossing of the fault current to be detected inthe SFA domain, which saves computational time. An interpolation scheme was in-troduced, which allows the use of large time-steps in the SFA solution while keepingthe high resolution of the reconstructed physical signal.In addition to applying the SFA solution method to the electrical network, the tran-sient stability electromechanical equation was generalized to allow the use of theknown machine velocity at every time-step of the simulation instead of using theapproximate 60-Hz frequency in the equation. This was implemented in the TSFAalgorithm and the two methods: using the correct velocity and using the approxi-mate 60-Hz velocity for the electromechanical equation, were compared in the TSFAalgorithm. There were noticeable differences between the two SFA solutions, in par-ticular when there were high frequency deviations away from the 60-Hz frequency,1208.1. Conclusionsuch as in the case when a generator was disconnected from the network.The TSFA solution was compared with the traditional phasor solution and theEMTP reference solution for the transient stability cases. These cases demonstratethat the TSFA solution matches the EMTP solution when calculating the criticalclearing times of a number of three-phase balanced faults. The TSFA also matchesthe EMTP solution when capturing the frequency dynamics of the electromechanicalvariables of the generators, rotor angle, frequency, electrical power, and the electri-cal variables, fault voltage and fault current, including capturing the DC offset ofthe fault current during a three-phase fault. In addition, TSFA matches the EMTPsolution when capturing the large frequency deviation from 60-Hz frequency whena generator is disconnected from the network, all while using large discretizationtime-steps. These findings validate SFA’s ability to capture the frequency devia-tions around the 60-Hz frequency, preserving the high fidelity of the EMTP solu-tion, while matching the computational time-efficiency of the fixed 60-Hz frequencyphasor solution. This is important for the modern power grid, where an increasingnumber of renewable energy sources create prominent frequency fluctuations aroundthe 60-Hz. Being able to track these deviations in a computationally time-efficientway is important.In addition to studying transient stability, the work done in this thesis re-formulatesthe SFA solution method in a generalized manner as discussed in Chapters 2 and3. The basis of SFA solution framework is the single-exponential, frequency-shiftingtransformation, which shifts the system from a 60-Hz domain to a 0-Hz domain [32].Because it is a rotational transformation, the shifting of frequencies preserves thedynamics from one domain to another and SFA is able to capture dynamics at alower frequency. Because SFA is a discrete-time solution method that adheres to thetime-frequency relationship, shifting the system down to deviate around 0 Hz allowslarge discretization time-steps to be used. In addition, because the reverse trans-formation is another rotational transformation, the dynamics captured in the lowerfrequency domain are maintained when transferred back to the physical domain. InChapter 3, a detailed analysis of the accuracy and behaviour during discontinuitiesof the trapezoidal and backward Euler discretization methods used in this work todiscretize the continuous-time SFA signals are discussed. The findings from thisanalysis indicate that the trapezoidal method is a more accurate solution method,1218.2. Future Researchespecially around DC, while the backward Euler method is more suited to networkdiscontinuities as its internal damping dampens the numerical oscillations that ap-pear when using the trapezoidal method.One observation from this research is that the SFA solution is not restricted tonetworks energized using sinusoidal sources, but can solve networks with arbitrarysources, such as DC sources, as demonstrated in Section 3.7. This is advantageousfor future work.The beauty of the SFA solution method comes from its ability to derive a complexworld from a frequency-shifting transformation of the physical world. There arean extensive number of future applications that can use SFA to model dynamicbehaviour.8.2 Future ResearchThe research presented in this thesis is one application of research using SFA. Be-cause of the advantages above, SFA can be used in modelling the behaviour of powersystems whenever there are system dynamics around the 60-Hz frequency. This sec-tion describes a few possible research ideas for SFA in both expanding the currentwork done in this thesis and in new topics.1. This thesis only covers three-phase balanced networks with three-phase bal-anced faults. The TSFA algorithm can be expanded for unbalanced networks,such as unbalanced fault conditions and unbalanced loads. By expanding theprogram for unbalanced conditions, SFA can be used in other dynamic stabilitystudies. For example, one prominent issue that occurs on both the transmis-sion and the distribution network is the fault-induced delayed voltage recovery(FIDVR) event. FIDVR events occur when a group of single-phase inductionmotors, such as air-conditioning motors stall, due to a single line-to-groundfault on the feeder, causing the bus voltage to drop and stay depressed fora long time. As in transient stability, the current stability simulators do notcapture the dynamics induced by the single line-to-ground fault accuratelyand the EMTP method is too computationally expensive to be used on thelarge network with many dynamic elements, such as the power electronic con-verters and many small induction motors [52]. Current approaches for this1228.2. Future Researchphenomenon include creating a hybrid simulator that interfaces traditionalphasors with the EMTP [52]. However, it may be possible to use SFA tomodel the induction motors and the electrical network to capture these dy-namics without needing to use a hybrid solution.2. The TSFA program can be extended by coupling the electrical and electrome-chanical network equation using the MATE [22] concept. Currently, the time-step used in TSFA is limited to around 10 ms because that is the highesttime-step needed by the electromechanical equations. If, however, the twonetworks were decoupled, the electrical network could be solved with as largetime-steps as necessary, as long as it complies with the time-frequency criteria,while the electromechanical network would be solved at a smaller time-step.The coupling would increase the solution’s computational speed.3. Using multi-rate techniques,[22],[51], a more accurate zero-crossing detectionmechanism can be done, which would allow the SFA solution to use the zero-crossing detection when calculating the critical clearing times. For example,the time-step would be smaller two time-steps before the fault is scheduled toclear and then once the fault is cleared, will return to the large time-steps.4. Research is currently being done in developing a hybrid simulation platformthat couples the EMTP solution with the SFA solution [43]. By employingthe Multi-Area The´venin Equivalent Multi-Rate (MATE-MR) framework [51],this hybrid solution can be an accurate and efficient platform for systems withmultiple time-scales. Because the SFA and the EMTP solution methods usethe same nodal analysis method and discretization of the continuous-timedifferential equations, the simulator would be able to capture the networkdynamics without approximations[43].5. The modelling of non-linear elements, such as non-linear inductors, can bedone in the SFA domain. In addition, the work here does not include thecontrollers and assumes that the mechanical power is constant and so thiswork can be expanded to include controller action.6. There has been work done in developing a discrete-time state-space modelusing the EMTP solution to follow the network’s trajectory and eigenvalues[53]. A derivation of this work can be done by deriving a SFA-state-space1238.2. Future Researchmodel. The SFA solution already solves for the network’s trajectory throughthe envelope’s dynamics, which identifies the network’s eigenvalues directly.There is potential research to be done in interfacing the SFA with state-spaceequations.7. Finally, the SFA framework is not limited to power system studies. The tech-nique of shifting frequencies down to more easily model dynamic behaviourcan be applied in the field of neuroscience. In the brain, the neural oscillationactivity in the different regions play an important role in how information isprocessed [54]. When information flows from one area to another, sometimesthe frequencies are crossed. This is called cross-frequency coupling (CFC)[54]and is the interaction between oscillations at the different frequency bands andhas been found to be prominent in regions of the brain associated with learningand memory [55]. One form of CFC is called phase-amplitude coupling, wherehigh-frequency oscillations cross with the phases of low frequency oscillations[55]. Modelling and measuring PAC is challenging, as the noise affects thehigh-frequencies of the brain activity in the measuring devices of the EEG andMEGs. One method recently proposed [55] is a dynamic phase-amplitude cou-pling, which uses Fourier decomposition to extract the magnitude and phaseof the signal to see dynamic changes, such as amplitude peaks. The idea ofshifting the frequencies of interest down to a lower frequency, as in SFA, canbe applied in a similar manner to develop a similar algorithm proposed in [55],in a discrete-time shifted frequency manner. Measuring CFC is important inunderstanding both normal and abnormal brain behaviour [54].124Bibliography[1] Benchmark systems for small-signal stability analysis and control. Technical re-port, IEEE PES: Power System Dynamic Performance Committee, Aug. 2015.[Online]. Available: http://resourcecenter.ieee-pes.org/pes/product/technical-publications/PESTR18. Retrieved: 2017-06-26. [Pages x, 78, 111,and 131.][2] J. R. Mart´ı. Lecture notes in Network Analysis and Simulation (EECE 560),1994-2004. The University of British Columbia, Vancouver,Canada. [Pages xii,39, 40, 41, 43, 46, 54, 61, and 62.][3] Powertech Labs Inc., Surrey, British Columbia, Canada. TSAT: TransientSecurity Assessment Tool User Manual. [Pages 1, 3, 8, 76, 78, and 80.][4] E. W. Kimbark. Elements of Stability Calculations. Power System Stability:Volume 1. John Wiley & Sons, Inc, 1948. [Pages 1, 2, 4, 5, 6, 71, 73, 74, 78,79, 80, 85, and 87.][5] NERC Resources Subcommittee. Balancing and frequency control. Techni-cal report, North American Electric Reliability Corporation, January 2011.[Page 2.][6] Australian Energy Market Operator. Final report of the investigation commit-tee on the 28 September 2003 blackout in Italy. Technical report, UCTE, April2004. [Page 2.][7] G. Lalor, A. Mullane, and M. O’Malley. Frequency control and wind turbinetechnologies. IEEE Transactions on Power Systems, 20(4):1905–1913, October2005. [Page 2.][8] B. Kroposki, B. Johnson, Y. Zhang, V. Geovorgian, P. Denholm, B. Hodge,and B. Hannegan. Achieving a 100% renewable grid. IEEE Power and EnergyMag., 15(2):61–73, March/April 2017. [Page 2.]125Bibliography[9] South Australian Advisory Functions. South Australian electricity report. Tech-nical report, Australian Energy Market Operator, November 2017. [Page 2.][10] Australian Energy Market Operator. Black system South Australian 28 septem-ber 2016. Technical report, Australian Energy Market Operator, March 2017.[Page 2.][11] G. Delille, B.Franc¸ois, and G.Malrange. Dynamic frequency control support byenergy storage to reduce the impact of wind and solar generation on isolatedpower system’s inertia. IEEE Transactions on Sustainable Energy, 3(4):931–939, August 2012. [Page 2.][12] International Energy Agency. World Energy Outlook 2017. [Online]. Available:https://www.iea.org/weo2017. Retrieved: 2017-09-01. [Page 3.][13] Power World Corporation, Champaign, Illinois, USA. Power World Simulator,2015. [Pages 3 and 8.][14] Microtran Power System Analysis Corporation, Vancouver, British Columbia,Canada. Microtran (R): Transients Analysis Program for Power and PowerElectronic Circuits. [Pages 3 and 78.][15] Power Systems Computer Aided Design (PSCAD), 211 Commerce Drive, Win-nipeg, Manitoba, Canada. User’s Guide on the Use of PSCAD. [Online]. Avail-able: https://hvdc.ca/uploads/ck/files/reference_material/PSCAD_User_Guide_v4_3_1.pdf. Retrieved: 2017-10-01. [Page 3.][16] AIEEE Committee Report. First report of power system stability. ElectricalEngineering, 56(2):261–282, February 1937. [Pages 4 and 6.][17] C. K. Alexander and M. N.O. Sadiku. Fundamentals of Electric Circuits.McGraw-Hill, 2009. [Page 7.][18] Siemens Power Technologies International. Siemens Corporation, Germany.Power Transmission System Planning Software. [Pages 8 and 76.][19] H. W. Dommel. Electromagnetic Transients Program (EMTP Theory Book).Bonneville Power Administration, 1986. [Pages 9, 11, 31, 39, 40, 61, 74, and 97.][20] J. R. Mart´ı. Lecture notes in Power System Analysis II (ELEC 454), April2017. [Pages 9, 10, 71, 72, 82, and 87.]126Bibliography[21] M. Lukic´. Feasibility of Security Assessment in Power Systems Using FullTime-Domain Solutions in the EMTP. MASC Thesis, The Univeristy of BritishColumbia, December 2000. [Pages 10, 75, and 76.][22] M. Armstrong, J. R. Mart´ı, L. R. Linares, and P. Kundur. Multilevel MATE forefficient simultaneous solution of control systems and nonlinearities in the OVNIsimulator. IEEE Transactions on Power Systems, 21(3):1250–1259, August2006. [Pages 10 and 123.][23] J. R. Mart´ı. Lecture notes in Discrete Time Signals and Systems, 2002. TheUniversity of British Columbia, Vancouver,Canada. [Pages 10 and 21.][24] J. R. Mart´ı, L. Mart´ı, and H. W. Dommel. Transmission line models for steady-state and transients analysis. In Athens Power Tech, 1993. APT 93. Pro-ceedings. Joint International Power Conference, pages 744–750. IEEE/NTUA,1993. [Page 11.][25] J. Zaborsky, H. Schattler, and V. Venkatasubramanian. Error estimation andlimitation of the quasi stationary phasor dynamics. In Power Systems Compu-tational Conference (PSCC 1993), pages 721–729, 1993. [Page 11.][26] S. Sanders, J. Noworolski, X. Liu, and G. Verghese. Generalized averagingmethod for power conversion circuits. IEEE Transactions on Power Electronics,6(2):251–259, April 1991. [Pages 11, 12, and 22.][27] A. M. Stankovic´, B.C. Lesieutre, and T. Aydin. Modeling and analysis of single-phase induction machines with dynamic phasors. IEEE Transactions on PowerSystems, 14(1):9–14, February 1999. [Pages 11, 12, and 22.][28] T. Demiray and G. Andersson. Simulation of power systems dynamics usingdynamic phasor models. In X Symposium of Specialists in Electric Operationaland Expansion Planning, May 2006. [Pages 11 and 12.][29] F. Gao and K. Strunz. Frequency-adaptive power system modeling formultiscale simulation of transients. IEEE Transactions on Power Systems,24(2):561–571, April 2009. [Page 11.][30] Y. Huang, M. Chapariha, S. Ebrahimi, N. Amiri, and J. Jatskevich. InterfacingSFA- and GAM-type dynamic phasors for modelling of integrated AC-DC power127Bibliographysystems. In Power and Energy Society General Meeting (PESGM), 2016, July2016. [Pages 11, 12, and 22.][31] S. Henschel. Analysis of Electromagnetic and Electromechanical Power Sys-tem Transients with Dynamic Phasors. PhD Thesis, The Univeristy of BritishColumbia, February 1999. [Pages 11, 13, and 14.][32] J.R.Mart´ı. Shifted Frequency Analysis (SFA) for EMTP simulation of fun-damental frequency power system dynamics. Internal Report, Power SystemLaboratory, The University of British Columbia, Vancouver, Canada, April2005. [Pages 11, 15, 16, 22, and 121.][33] P. Zhang, J. R. Mart´ı, and H. W. Dommel. Shifted Frequency Analysis forEMTP simulation of power-system dynamics. IEEE Transactions on Circuitsand Systems, 57(9):2564–2574, 2010. [Pages 11, 16, 18, and 22.][34] V. Venkatasubramanian, H. Schattler, and J. Zaborszky. Fast time-varyingphasor analysis in the balanced three-phase large electric power system.IEEE Transactions on Automatic Control, 40(11):1975–1982, November 1995.[Pages 11, 12, 13, and 14.][35] V. Venkatasubramanian, H. Schattler, and J. Zaborsky. A time-delaydifferential-algebraic phasor formulation of the large power system dynamics.In Proceedings of IEEE International Symposium on Circuits and Systems -ISCAS ’94, volume 6, pages 49–52, May 1994. [Page 11.][36] C. DeMarco and G. Verghese. Bringing phasor dynamics into the power systemload flow. In 1993 North American Power Symposium, pages 463–471, October1993. [Pages 11 and 12.][37] V. Venkatasubramanian. Dynamic analysis of the general large power systemusing time-varying phasors. [Pages 12 and 14.][38] J. R. Mart´ı and J. Lin. Suppression of numerical oscillations in the EMTP.IEEE Transactions on Power Systems, 4(2):739–747, May 1989. [Pages 15, 29,39, 40, 54, 61, 79, and 94.][39] J. R. Mart´ı, H. W. Dommel, B. D. Bonatto, and A. F. R. Barreto. Shifted Fre-quency Analysis (SFA) concepts for EMTP modelling and simulation of power128Bibliographysystem dynamics. In 2014 Power Systems Computation Conference (PSCC),pages 1–8, Wroclaw, Poland, Aug. 2014. [Pages 16, 18, and 22.][40] P. Zhang, J. R. Mart´ı, and H. W. Dommel. Induction machine modelingbased on Shifted Frequency Analysis. IEEE Transactions on Power Systems,24(1):157–164, February 2009. [Page 16.][41] P. Zhang, J. R. Mart´ı, and H. W. Dommel. Synchronous machine modelingbased on Shifted Frequency Analysis. IEEE Transactions on Power Systems,22(3):1139–1147, February 2007. [Page 16.][42] Y. Huang, F. Therrien, J. Jatskevich, and L. Dong. State-space voltage-behind-reactance modeling of induction machines based on Shifted Frequency Analysis.In 2015 IEEE Power & Energy Society General Meeting (PESGM), July 2015.[Page 16.][43] J. Tarazona. Hybrid Shifted Frequency Analysis Electromagnetic TransientsMultirate Simulation of Power Systems. Internal Proposal: A proposal forthe PhD Dissertation, The University of British Columbia, Vancouver,Canada,August 2017. [Pages 16 and 123.][44] P. Kundur, J. Paserba, V. Ajjarapu, G. Andersson, A. Bose, C. Canizares,N. Hatziargyriou, D. Hill, A. Stankovic´, C. Taylor, T. Van Cutsem, and V. Vit-tal. Definition and classification of power system stability, IEEE/CIGRE JointTask Force on Stability Terms and Definitions. IEEE Transactions on PowerSystems, 19(2):1387–1401, August 2004. [Page 71.][45] J.Machowski, J. Bialek, and J. Bumby. Power System Dynamics: Stability andControl. John Wiley & Sons, Inc., 2008. [Pages 72, 78, and 96.][46] H.Bevrani, M.Watanabe, and Y.Mitani. Power System Monitoring and Control.John Wiley & Sons, Inc., 2014. [Page 73.][47] E. W. Kimbark. Power Circuit Breakers and Protective Relays. Power SystemStability: Volume 2. John Wiley & Sons, Inc, 1948. [Pages 73, 74, 81, and 82.][48] E. W. Kimbark. Synchronous Machines. Power System Stability: Volume 3.John Wiley & Sons, Inc, 1948. [Pages 74, 75, 80, and 81.][49] Y.Yu. Electric Power System Dynamics. Academic Press, 1983. [Page 80.]129Bibliography[50] Powertech Labs Inc., Surrey, BC, Canada. PSAT: Powerflow & Short-circuitAnalysis Tool User Manual. [Page 81.][51] V. A. Galva´n, J. R. Mart´ı, H. W. Dommel, J. A. Gutie´rrez Robles, and J. L.Naredo. MATE multirate modelling of power electronic converters with mixedintegration rules. In Power Systems Computation Conference (PSCC), 2016,August 2016. [Pages 94 and 123.][52] Q. Huang and V. Vittal. Application of electromagnetic transient-transient sta-bility hybrid simulation to fidvr study. IEEE Transactions on Power Systems,31(4):2634 – 2646, 2016. [Pages 122 and 123.][53] J. A. Hollman and J. R. Mart´ı. Step-by-Step eigenvalue analysis with EMTPdiscrete-time solutions. IEEE Transactions on Power Systems, 25(3):1220–1231, August 2010. [Page 123.][54] P. J. Uhlhaas and W. Singer. High-frequency oscillations and the neurobiol-ogy of schizophrenia. Dialogues in Clinical Neuroscience, 15(3):301–313, 2013.[Page 124.][55] S. S. Samiee and S. Baillet. Time-resolved phase-amplitude coupling in neuraloscillations. NeuroImage, 159:270–279, 2017. [Page 124.]130AppendixTables A1 A2 A3 and A4 are the modified data for the IEEE 39-bus case study [1].Table A1Generator data for the IEEE 39-bus case study (Modified from [1]).GeneratorsPowerRating(MVA)VoltageRating(kV)x′d(per-unit1)InertiaConstant(H)(sec)Moment ofInertia (I)(kg·m2)DampingCoefficient(K)(per-unit2)Generator 1 1100 100 0.3 5 77398.1 1.5Generator 2 350 22 0.3 3 14776.0 1.5Generator 3 750 22 0.3 3 31662.8 1.5Generator 4 650 22 0.3 4 36588.2 1.5Generator 5 50 11 0.3 1 703.6 1.5Generator 6 600 22 0.3 5.5 46438.8 1.5Generator 7 700 22 0.3 2 19701.3 1.5Generator 8 600 22 0.3 2.5 21108.5 1.5Generator 9 900 22 0.3 3.5 44328.0 1.5Generator 10 350 22 0.3 6 29552.0 1.51 per-unit of each generator’s system base 2 per-unit of torque baseTable A2Transformer data for the IEEE 39-bus case study.From Bus To BusR(per-unit1)X(per-unit1)6 31 0 0.02510 32 0 0.02012 11 0.0016 0.043512 13 0.0016 0.043519 33 0.0007 0.014219 20 0.0007 0.013820 34 0.0009 0.018023 36 0.0005 0.027229 38 0.0008 0.015630 2 0 0.018135 22 0 0.014337 25 0.0006 0.02321 Transformers are rated 100kV:22kV and reactances are per-uniton the 100kV, 100 MVA base131AppendixTable A3Power flow results for the IEEE 39-bus case study.BusV(per-unit)Angle(Deg)PGen(MW)QGen(MVAR)PLoad(MW)QLoad(MVAR)1 1.0364 1.16352 1.0185 3.22223 0.9908 -0.2058 322.00 2.404 0.9557 -2.1849 500.00 184.005 0.9555 -1.70486 0.9567 -1.15457 0.9488 -3.2679 233.80 84.008 0.9494 -3.6643 522.00 176.009 1.0088 -1.481410 0.9628 2.077911 0.9593 0.984612 0.9398 1.1603 7.50 88.0013 0.9612 1.453314 0.9617 0.041715 0.9695 0.7541 320.00 153.0016 0.9890 2.7813 329.00 32.3017 0.9926 1.481218 0.9904 0.3742 158.00 30.0019 0.9951 10.038420 1.0021 10.0973 30.00 4.9021 0.9961 5.0020 274.00 115.0022 1.0220 9.353523 1.0201 9.6455 247.50 84.6024 0.9972 2.8785 308.60 -92.2025 1.0276 4.7982 224.00 57.2026 1.0174 3.4692 139.00 17.0027 0.9997 1.3200 281.00 75.5028 1.0189 7.1934 206.00 27.6029 1.0205 10.1108 283.50 26.9030 1.0475 5.6530 250.00 173.0931 0.9820 3.2836 300.00 115.14 9.20 4.6032 0.9831 9.9720 650.00 144.5733 0.9972 15.2223 632.00 11.9834 1.0123 10.4757 40.00 55.6435 1.0493 13.2379 508.00 217.8136 1.0635 18.9686 650.00 210.5037 1.0278 11.8574 560.00 21.0238 1.0265 17.1898 830.00 48.3139 1.0300 0.0000 23.25 -41.01132AppendixTable A4Transmission line data for the IEEE 39-bus case study.From Bus To BusR(per-unit1)X(per-unit1)B(total2)(per-unit1)39 1 0.0010 0.0250 0.750039 9 0.0010 0.0250 1.200028 29 0.0014 0.0151 0.249027 17 0.0013 0.0173 0.321626 27 0.0014 0.0147 0.239626 28 0.0043 0.0474 0.780226 29 0.0057 0.0625 1.029025 26 0.0032 0.0323 0.513024 23 0.0022 0.0350 0.361022 23 0.0006 0.0096 0.184622 21 0.0008 0.0140 0.256521 16 0.0008 0.0135 0.254817 18 0.0007 0.0082 0.131916 24 0.0003 0.0059 0.068016 17 0.0007 0.0089 0.134216 19 0.0016 0.0195 0.304015 16 0.0009 0.0094 0.171014 15 0.0018 0.0217 0.366013 10 0.0004 0.0043 0.072913 14 0.0009 0.0101 0.172311 10 0.0004 0.0043 0.07299 8 0.0023 0.0363 0.38048 7 0.0004 0.0046 0.07807 6 0.0006 0.0092 0.11306 11 0.0007 0.0082 0.13895 6 0.0002 0.0026 0.04345 4 0.0008 0.0128 0.13425 8 0.0008 0.0112 0.14764 14 0.0008 0.0129 0.13823 2 0.0013 0.0151 0.25723 18 0.0011 0.0133 0.21383 4 0.0013 0.0213 0.22142 25 0.0070 0.0086 0.14601 2 0.0035 0.0411 0.69871 per-unit on the 100kV, 100 MVA base 2 total susceptance133
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The application of Shifted Frequency Analysis in power...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The application of Shifted Frequency Analysis in power system transient stability studies Martí, Andrea T.J. 2018
pdf
Page Metadata
Item Metadata
Title | The application of Shifted Frequency Analysis in power system transient stability studies |
Creator |
Martí, Andrea T.J. |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Power system engineers use transient stability computer simulation programs to model the power grid’s behaviour when large disruptions cause the grid to deviate from its 60-Hz operating frequency. These programs must be able to capture these frequency dynamics around 60 Hz while being computationally efficient, as an extensive number of simulations are typically run for a given scenario. In the traditional power grid, the large mechanical inertia of the synchronous generators stabilizes the network during disturbances and maintains the system frequency close to the 60 Hz operating frequency. In the modern grid, however, the increase of renewable energy sources lowers the grid’s inertia and larger frequency deviations can occur. The phasor solution method employed in the traditional programs solves the network assuming a constant 60-Hz frequency. When deviations from 60 Hz are prominent, the Electromagnetic Transients Program (EMTP) is used as an alternative to the phasor solution to capture these fluctuations. The EMTP models the electrical network based on the differential equations of the network components, which allows the tracing of the network waveforms. However, this discretization requires small time-steps, which makes the solution method computationally expensive. The Shifted Frequency Analysis (SFA) method discussed in this work is an alternative to the traditional phasor solution and to the EMTP solution. In this work, a generalized SFA-based program is written and used for transient stability analysis. SFA is a discrete-time solution method, like the EMTP, but uses a frequency-shifting transformation to bring the solution domain down to 0 Hz. Because of this transformation, SFA can capture network dynamics around 60 Hz using large time-steps, making it suitable for transient stability analysis studies. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-02-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
IsShownAt | 10.14288/1.0364058 |
URI | http://hdl.handle.net/2429/64709 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_may_marti_andrea.pdf [ 14.14MB ]
- Metadata
- JSON: 24-1.0364058.json
- JSON-LD: 24-1.0364058-ld.json
- RDF/XML (Pretty): 24-1.0364058-rdf.xml
- RDF/JSON: 24-1.0364058-rdf.json
- Turtle: 24-1.0364058-turtle.txt
- N-Triples: 24-1.0364058-rdf-ntriples.txt
- Original Record: 24-1.0364058-source.json
- Full Text
- 24-1.0364058-fulltext.txt
- Citation
- 24-1.0364058.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0364058/manifest