NUMERICALLY EFFICIENT MODELLING AND SIMULATION OF INTEGRATED AC-DC POWER SYSTEMS USING DYNAMIC PHASOR TYPE SOLUTIONS by Yingwei Huang B.Eng. Tianjin University, Tianjin, China, 2012 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate and Postdoctoral Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2017 © Yingwei Huang, 2017 ii Abstract Modern power systems, from continent-spanning networks down to isolated microgrids, are experiencing unprecedented technological changes with broader use of direct current (dc) in addition to traditional alternating current (ac). Such integrated ac-dc power systems present notable challenges in all aspects of design, analysis, control, and operation, where extensive computer simulations play the essential and enabling role. Due to the use of diverse types of signal representation and component formulation, state-of-the-art power system simulation tools are limited to their distinct time scales of transient phenomena. This thesis considers the dynamic phasor (DP) type modelling approaches, where two types of DP theories, namely the shifted-frequency analysis (SFA) and the generalized averaging method (GAM), are considered. In DP-type simulations, power systems are modelled using analogous low-pass time-phasor signals, thereby offering flexible selection of time-step sizes and superior combination of numerical accuracy and efficiency. The ultimate goal of this research is to increase the numerical efficiency of DP-type simulations for the integrated ac-dc power systems. This is achieved by proposing several new DP component models with desirable features and improved numerical properties. First, the constant-parameter SFA model of synchronous machines is proposed to avoid numerically-costly recalculations of the time-varying stator-network matrices. This model is then extended to induction machines for modelling in state-variable based (SV-based) simulation tools. Next, we propose a new, highly efficient model of line-commutated rectifiers using a parametric DP formulation, which is demonstrated as valid for various system operating conditions. Moreover, the effect of ac side harmonics is incorporated to improve modelling fidelity. Finally, the interface between SFA- and GAM-type DPs is achieved to interconnect the proposed DP models. Rigorous case studies demonstrate the superior numerical efficiency of the proposed models, and their advantageous accuracy in capturing the desired phenomena of ac-dc power systems. It is envisioned that the proposed models will become highly useful to many researchers and engineers worldwide, and facilitate the development of next-generation power system simulation tools. iii Lay Summary The emergent ac-dc power systems necessitate a massive number of computer simulations to analyze their dynamic behaviour. Several types of simulation tools have been developed throughout the years, each targeted for investigation of specific classes of phenomena. This research focuses on the dynamic phasor (DP)-type modelling approaches, which are demonstrated to provide general-purpose simulations of power system transients with a superior combination of numerical accuracy and efficiency. The objective of this thesis is to increase the numerical efficiency of DP-type simulations for integrated ac-dc power systems. This is achieved by proposing several new DP models of system components, including electric machines and line-commutated rectifiers, with improved numerical properties and more desirable features over prior state-of-the-art models. It is envisioned that the proposed models will become highly useful to many researchers and engineers worldwide, and facilitate the development of next-generation power system simulation tools.iv Preface Many of the research results presented in this dissertation have been published in scientific journals, conference proceedings, and/or will be submitted for peer review. In all publications, I have derived all the mathematical formulations, implemented the proposed models, performed the computer studies, analyzed the simulation results, and was responsible for writing the first draft of each paper. My supervisor, Dr. Juri Jatskevich, has provided guidance and instructive comments throughout the process of conducting research, preparing manuscripts, and improving their revisions. Other co-authors have also contributed to the research by providing their feedback and editorial comments, and by editing and revising the manuscripts. The articles resulting from this doctoral work are listed below. Chapter 3 is based on the following journal article and conference paper: Y. Huang, M. Chapariha, F. Therrien, J. Jatskevich, and J. R. Martí, "A constant-parameter voltage-behind-reactance synchronous machine model based on shifted-frequency analysis," IEEE Trans. Energy Convers., vol. 30, no. 2, pp. 761-771, Jun. 2015. Y. Huang, F. Therrien, J. Jatskevich, and L. Dong, "State-space voltage-behind-reactance modeling of induction machines based on shifted-frequency analysis," in Proc. IEEE Power Energy Soc. Gen. Meet., Denver, CO, USA, Jul. 2015, pp. 1-5. The constant-parameter interfacing technique for VBR machine models presented in Section 3.1.1.2 was proposed by Dr. M. Chapariha. Chapter 4 is based on the following journal articles and conference paper: Y. Huang, L. Dong, S. Ebrahimi, N. Amiri, and J. Jatskevich, "Dynamic-phasor-type modeling of line-commutated rectifiers with harmonics using analytical and parametric approaches," IEEE Trans. Energy Convers., vol. 32, no. 2, pp. 534-547, Jun. 2017. Y. Huang, S. Ebrahimi, N. Amiri, Z. Shan, and J. Jatskevich, "Parametric dynamic phasor modeling of thyristor-controlled rectifier systems including harmonics for various operating modes," IEEE Trans. Energy Convers., vol. 32, no. 4, pp. 1626-1629, Dec. 2017. v Y. Huang, S. Ebrahimi, N. Amiri, H. Atighechi, and J. Jatskevich, "A parametric dynamic phasor model of line-commutated rectifier systems," in Proc. IEEE Canadian Conf. Electri. Comput. Eng., Vancouver, Canada, May. 2016, pp. 1-4. The parametric modelling technique in Section 4.1.3 was proposed by Dr. J. Jatskevich. Chapter 5 is based on the following conference papers: Y. Huang, S. Ebrahimi, N. Amiri, M. Chapariha, and J. Jatskevich, "Interfacing SFA- and GAM-type dynamic phasors for modeling of integrated AC/DC systems," in Proc. IEEE Power Energy Soc. Gen. Meet., Boston, MA, USA, Jul. 2016, pp. 1-5. Y. Huang, H. Chang, S. Rezaee, and J. Jatskevich, "Parametric dynamic phasor modeling of synchronous machine-rectifier systems for integrated AC-DC microgrids," in Proc. IEEE Power Energy Soc. Gen. Meet., Chicago, IL, USA, Jul. 2017, pp. 1-5. Finally, Chapters 1 and 2 of this thesis contain excerpts of the introductory sections of the aforementioned papers. vi Table of Contents Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iii Preface ........................................................................................................................................... iv Table of Contents ......................................................................................................................... vi List of Tables ................................................................................................................................ xi List of Figures ............................................................................................................................. xiii List of Abbreviations ................................................................................................................. xix Nomenclature ............................................................................................................................. xxi Acknowledgments .................................................................................................................... xxiv Dedication ................................................................................................................................. xxvi Introduction ........................................................................................................... 1 Motivation ............................................................................................................................ 1 Literature Review ................................................................................................................. 7 Dynamic Phasor (DP) Type Modelling Approaches ................................................. 7 1.2.1.1 Shifted-Frequency Analysis (SFA) .................................................................. 7 1.2.1.2 Generalized Averaging Method (GAM) .......................................................... 9 Numerically Efficient Modelling of AC-DC Power System Components .............. 10 1.2.2.1 Modelling of Electric Machines .................................................................... 11 1.2.2.2 Modelling of Line-Commutated Rectifiers (LCRs) ...................................... 13 Research Objectives and Anticipated Impact .................................................................... 15 Dynamic Phasor Modelling Approaches for Power System Transient Simulations .......................................................................................................... 19 Conventional Signal Representations in Power System Simulations ................................ 19 vii Time-Domain Instantaneous Signals ....................................................................... 21 2.1.1.1 abc-Phase Coordinates ................................................................................... 21 2.1.1.2 qd0 Coordinates ............................................................................................. 22 Conventional Fundamental-Frequency Phasors ...................................................... 24 DP-Type Modelling Approaches ....................................................................................... 25 Shifted-Frequency Analysis (SFA) .......................................................................... 26 Generalized Averaging Method (GAM) .................................................................. 28 Discussions and Computer Studies .................................................................................... 31 Modelling Basic RLC Components ......................................................................... 32 Modelling System Transients at Different Frequencies ........................................... 33 2.3.2.1 Analytical Modelling and Eigenvalue Analysis ............................................ 34 2.3.2.2 Computer Studies........................................................................................... 39 Modelling Power Systems under Unbalanced Conditions....................................... 46 2.3.3.1 SFA-Type DPs of Three-Phase Power System Signals ................................. 46 2.3.3.2 GAM-Type DPs of Three-Phase Power System Signals ............................... 48 2.3.3.3 Computer Studies........................................................................................... 50 Numerically Efficient Modelling of Electric Machines Based on the Shifted Frequency Analysis ............................................................................................. 57 Constant-Parameter VBR Modelling of Synchronous Machines Based on SFA-Type DPs ...................................................................................................................................... 58 Time-domain VBR Synchronous Machine Models ................................................. 58 3.1.1.1 Variable-Parameter Time-Domain VBR Model (VP-TD) ............................ 58 3.1.1.2 Constant-Parameter Time-Domain VBR Model (CP-TD) ............................ 60 viii VBR Synchronous Machine Models Based on SFA ............................................... 62 3.1.2.1 Variable-Parameter Stator Interface Based on SFA ...................................... 63 3.1.2.2 Constant-Parameter Stator Interface Based on SFA ...................................... 64 3.1.2.3 Rotor Subsystem and Its Interface ................................................................. 65 Computer Studies ..................................................................................................... 68 3.1.3.1 Transition from Fast Transient to Steady State ............................................. 69 3.1.3.2 Slow Transient ............................................................................................... 74 State-Space VBR Modelling of Induction Machines Based on SFA-Type DPs ............... 78 Time-Domain VBR Induction Machine Model (VBR-TD-IM) .............................. 78 VBR Induction Machine Model Based on SFA ...................................................... 80 3.2.2.1 SFA-Type DP Stator Interface....................................................................... 80 3.2.2.2 Rotor Subsystem and Its Interface ................................................................. 81 Computer Studies ..................................................................................................... 84 Modelling Line-Commutated Rectifier Systems Including Harmonics for Various Operating Modes Based on the Generalized Average Method ........ 89 GAM-Type DP Modelling of Diode LCR Systems with Harmonics Using Analytical and Parametric Approaches ................................................................................................ 90 Time-Domain Dynamics of LCR Systems and Switch Functions........................... 91 4.1.1.1 Piecewise-Linear-Approximated Switch Functions (SF-1) ........................... 93 4.1.1.2 Fourier-Series-Approximated Switch Functions (SF-2) ................................ 94 Analytical DP Modelling ......................................................................................... 95 4.1.2.1 Analytical DP Model Using SF-1 (ADP-1) ................................................... 95 4.1.2.2 Analytical DP Model Using SF-2 (ADP-2) ................................................... 96 ix Parametric DP Modelling (PDP) ............................................................................. 97 Constructing Parametric Functions .......................................................................... 99 Implementation ...................................................................................................... 103 Computer Studies ................................................................................................... 105 4.1.6.1 Steady State Analysis in CCM-1 and DCM ................................................ 106 4.1.6.2 Transient Response from CCM-1 to DCM .................................................. 110 4.1.6.3 Predicting THD in Wide Range of Operating Conditions ........................... 110 4.1.6.4 Flexible Step Size Study .............................................................................. 112 4.1.6.5 Discussions of the Ringing of Harmonic DPs ............................................. 116 4.1.6.6 Fixed Step Size Study .................................................................................. 119 Parametric DP Modelling of Thyristor-Controlled LCR Systems Including Harmonics for Various Operating Modes .......................................................................................... 121 Relating the AC-DC DP Dynamics ....................................................................... 122 Constructing Parametric Functions ........................................................................ 123 Computer Studies ................................................................................................... 123 4.2.3.1 Steady-State Response ................................................................................. 124 4.2.3.2 Large-Signal Transient Response ................................................................ 125 4.2.3.3 Frequency-Domain Study ............................................................................ 130 Interfacing SFA- and GAM-Type DPs for Modelling the Synchronous Machine-Rectifier Systems ............................................................................... 132 DP Modelling of Synchronous Machine-Rectifier System Components ........................ 132 Synchronous Machine Model in SFA-Type DPs................................................... 133 Diode LCR Model in GAM-Type DPs .................................................................. 135 x DC Subsystem Model in the Time Domain ........................................................... 136 SFA- and GAM-Type DP Interface and System Implementation ................................... 136 Constructing Parametric Functions .................................................................................. 137 Computer Studies ............................................................................................................. 139 Small Step Size Study ............................................................................................ 139 Flexible Step Size Study ........................................................................................ 142 Conclusions and Future Work ......................................................................... 144 Conclusions and Contributions ........................................................................................ 144 Future Work ..................................................................................................................... 146 References .................................................................................................................................. 149 Appendix A: Parameters for the Case Studies of Section 3.1.3 ............................................ 162 Appendix B: Parameters for the Case Studies of Section 3.1.3 ............................................ 164 Appendix C: Parameters for the Case Studies of Section 3.2.3 ............................................ 165 Appendix D: Parameters for the Case Studies of Sections 4.1.6 and 4.2.3 .......................... 166 xi List of Tables Table 1–1. Properties of state-of-the-art and proposed electric machine models. ........................ 15 Table 1–2. Properties of state-of-the-art and proposed models of LCR systems. ........................ 16 Table 2–1. Properties of different representations of power system signals. ............................... 32 Table 2–2. Dynamic equations governing the basic RLC components using different signal representations. ............................................................................................................................. 33 Table 2–3. Properties of system models for the series RLC circuit using different signal representations. ............................................................................................................................. 38 Table 2–4. Parameters of the series RLC circuit for different resonant frequencies. ................... 41 Table 2–5. Eigenvalues of the subject models for the series RLC circuit with different resonant frequencies. ................................................................................................................................... 43 Table 2–6. Properties of different signal representations for the three-phase power systems under balanced/unbalanced conditions. .................................................................................................. 49 Table 2–7. Parameters of the considered three-phase transmission system. ................................ 51 Table 3–1. Simulation Efficiency and Accuracy for the Fast-transient-to-steady-state Study. .... 74 Table 3–2. Simulation Efficiency and Accuracy for the Slow-transient Study. ........................... 77 Table 3–3. Simulation Efficiency and Accuracy of the Subject Models ...................................... 88 Table 4–1. Models Simulation Efficiency of the Subject Models for the Flexible Step Size Case Study Using a Step Change of Load Resistance. ........................................................................ 116 Table 4–2. Simulation Efficiency of the Subject Models for the Flexible Step Size Case Study Using a Sigmoid Change of Load Resistance. ............................................................................ 119 xii Table 4–3. Total CPU Time (in Sec) Taken by the Subject Models Solved Using Fixed Step Trapezoidal Rule ......................................................................................................................... 120 Table 4–4. Parameters of the close-loop thyristor firing controller for dc voltage regulation. .. 126 Table 5–1. Simulation Efficiency of the Subject Models. .......................................................... 143 xiii List of Figures Figure 1–1. The ac-dc microgrid system presently installed (solid lines) and to be further developed (dashed lines) in the Kaiser Building at the University of British Columbia, Vancouver, Canada [11]. ................................................................................................................ 2 Figure 1–2. Time frame of different power system transient phenomena. ..................................... 4 Figure 2–1. Frequency spectrum of a band-pass power system signal tu represented by the sum of closely spaced sine waves. ........................................................................................................ 28 Figure 2–2. Frequency spectra of the band-pass signal in different representations during the process of SFA: a) instantaneous signal tu ; b) analytic signal tUˆ ; and c) SFA-type DP tU ........................................................................................................................................................ 28 Figure 2–3. Frequency spectra of the power system signal in different representations during the process of GAM: a) instantaneous signal, and b) GAM-type DP. ................................................ 30 Figure 2–4. Illustrated “sliding window” of the signal moving along the time axis. ................... 30 Figure 2–5. The example single-phase series RLC circuit in different signal representations: a) instantaneous signal; b) conventional phasor; c) SFA-type DP; and d) GAM-type DP. .............. 34 Figure 2–6. System response of the RLC circuit with resonant frequency Hz 10nf . ................ 41 Figure 2–7. System response of the RLC circuit with resonant frequency Hz 300nf . ............... 42 Figure 2–8. Transfer function from input Sv to output Li as predicted by the TD model for the series RLC circuit with resonant frequency Hz 10nf . ................................................................ 44 Figure 2–9. Transfer function from input TiSrS VV ] ,[ ,, to output TiSrS II ] ,[ ,, as predicted by the SFA model for the series RLC circuit with resonant frequency Hz 10nf . .......................................... 44 xiv Figure 2–10. Transfer function from input Sv to output Li as predicted by the TD model for the series RLC circuit with resonant frequency Hz 300nf . ............................................................... 45 Figure 2–11. Transfer function from input TiSrS VV ] ,[ ,, to output TiSrS II ] ,[ ,, as predicted by the SFA model for the series RLC circuit with resonant frequency Hz 300nf . ................................ 45 Figure 2–12. Circuit diagram of the considered 230 kV three-phase transmission system. ........ 50 Figure 2–13. System response of csi as predicted by the subject models for the single-phase-ground fault in phase-c voltage source. ........................................................................................ 53 Figure 2–14. Trajectories of system states during the fault transient for the subject models: (a) Phasor; (b) TD-abc; (c) TD-qd0; (d) DP-abc; and (e) DP-qd0. ................................................... 54 Figure 2–15. Trajectories of system states after the fault transient for the subject models: (a) Phasor; (b) TD-abc; (c) TD-qd0; (d) DP-abc; and (e) DP-qd0. .......................................... 54 Figure 2–16. Step size Δt as chosen by the subject models. ........................................................ 55 Figure 3–1. Implementation of the proposed CP-SFA VBR model. ............................................ 67 Figure 3–2. Fast transient and steady-state responses to a three-phase fault as predicted by the subject models. .............................................................................................................................. 72 Figure 3–3. Detailed view of current asi during the three-phase fault: (a) magnified view during fast transient; (b) further magnified view of the portion in part (a). ............................................. 73 Figure 3–4. Detailed view of current asi after removal of the three-phase fault when the system enters steady state. ........................................................................................................................ 73 Figure 3–5. Step size Δt as chosen by the subject models: (a) during the entire fault study; and (b) magnified view from part (a). ................................................................................................. 74 xv Figure 3–6. Slow transient response of the system to a torque change as predicted by the subject models. .......................................................................................................................................... 76 Figure 3–7. Magnified view of current asi for the slow transient study depicted in Figure 3–6. . 76 Figure 3–8. Step size Δt as chosen by the subject models for the slow transient study. .............. 77 Figure 3–9. Block diagram for implementation of the VBR-SFA-IM model. ............................. 83 Figure 3–10. Start-up and load-change transient responses as predicted by the subject models. . 87 Figure 3–11. Magnified view of stator current asi as depicted in Figure 3–10 during: (a) start-up transient; (b) torque-change transient. .......................................................................................... 88 Figure 3–12. Step size Δt as chosen by the subject models. ......................................................... 88 Figure 4–1. Typical configuration of the three-phase six-pulse LCR system. ............................. 90 Figure 4–2. Time-domain waveforms of phase a voltage and current switch functions: (a) vas of SF-1, (b) vas of SF-2, (c) ias of SF-1, and (d) ias of SF-2. ........................................................ 92 Figure 4–3. The large-signal transient response to an exponential increase of load resistance for calculating parametric functions. ................................................................................................ 101 Figure 4–4. Numerically calculated functions )(1 z , )(1 z , and )(1 z as obtained from detailed simulations using (a) the MSS approach; and (b) the LST approach. ........................................ 102 Figure 4–5. Numerically calculated functions )(7,5 zα and )(7,5 zφ as obtained from detailed simulations using (a) the MSS approach; and (b) the LST approach. ........................................ 102 Figure 4–6. Implementation of the DP models: (a) ADP-1 and ADP-2 models, (b) the proposed PDP model. ................................................................................................................................. 105 Figure 4–7. Ac and dc system variables for steady-state in CCM-1 as predicted by the subject models. ........................................................................................................................................ 108 xvi Figure 4–8. Harmonic content of the ac-side phase a current and voltage for the considered operating point in CCM-1. .......................................................................................................... 108 Figure 4–9. Ac and dc system variables for steady-state response in DCM as predicted by the subject models. ............................................................................................................................ 109 Figure 4–10. Harmonic content of the ac-side phase a current and voltage for the considered operating point in DCM. ............................................................................................................. 109 Figure 4–11. Ac and dc system variables for the load-change transient as predicted by the subject models. ............................................................................................................................ 110 Figure 4–12. THD of the phase a voltage and current as predicted by the subject models over a wide range of operating conditions. ............................................................................................ 111 Figure 4–13. Step size Δt as chosen by each model for the flexible step sizes study using a step change of load resistance. ........................................................................................................... 114 Figure 4–14. Current asi as predicted by the subject models using flexible step sizes for steady-state response in CCM-1. ............................................................................................................ 114 Figure 4–15. System transient response to the step change of load resistance as predicted by the subject models using flexible time steps: (a) current asi , (b) 5th order DP 52 asi , (c) 7th order DP 72 asi , and (d) magnified view from subplot (a). ..................................................................... 115 Figure 4–16. Step size Δt as chosen by each model for the flexible step sizes study using a sigmoid change of load resistance. ............................................................................................. 118 xvii Figure 4–17. System transient response to the sigmoid change of load resistance as predicted by the subject models using flexible time steps: (a) current asi , (b) 5th order DP 52 asi , and (c) 7th order DP 72 asi . ........................................................................................................................ 118 Figure 4–18. Current asi transient response to the load step-change as predicted by the subject models using fixed step size Δt = 1ms. ....................................................................................... 120 Figure 4–19. Benchmark configuration of the thyristor-controlled LCR systems. .................... 121 Figure 4–20. Implementation of the PDP model of thyristor-controlled rectifiers..................... 122 Figure 4–21. Parametric functions as numerically calculated from detailed simulations: (a) ,1, dv z ; (b) ,di z ; (c) ,1 dz ; (d) ,5, dv z ; (e) ,5 dz ; (f) ,7, dv z ; and (g) ,7 dz . ..................................................................................................................................... 124 Figure 4–22. Steady-state responses: THD % of asi for different thyristor firing angles. ........ 125 Figure 4–23. System response to the large-signal transients: (a) profile of lR ; (b) profile of ; (c) resulting dc voltage dcv ; and (d) resulting dc current dci . ..................................................... 127 Figure 4–24. Magnified views of ac variables asv and asi response to: (a) step change of lR ; (b) ramp change of lR ; (c) step change of ; (d) ramp change of . ...................................... 128 Figure 4–25. Magnified views of system close-loop dynamic response to: (a) activation of the dc voltage controller; (b) step change of the reference voltage. ................................................. 129 Figure 4–26. The close-loop thyristor firing controller for dc voltage regulation. ................... 129 Figure 4–27. Transfer function from the input firing angle to the output dc voltage as predicted by the subject models. ................................................................................................................. 131 xviii Figure 5–1. Circuit diagram of the synchronous machine–rectifier system for ac-dc backup generation. ................................................................................................................................... 133 Figure 5–2. Block diagram of the DP modelling of synchronous machine-rectifier systems. 137 Figure 5–3. Parametric functions as obtained from detailed simulations: (a)(b)(c) )(1 z , )(1 z , and )(1 z , and (d)(e) )(7,5 zα and )(7,5 zφ . ................................................................................. 138 Figure 5–4. System response to the load change as predicted by the subject models. ............ 140 Figure 5–5. Magnified view of the ac variables in steady states as depicted in Figure 5–4: (a)(b) asv and asi in the original CCM-1 mode; and (c)(d) asv and asi in the new CCM-2 mode. ........ 141 Figure 5–6. Magnified view of the ac and dc variables during the load-change transient as depicted in Figure 5–4. ............................................................................................................... 142 Figure 5–7. Step size Δt as chosen by the subject models. ......................................................... 143 xix List of Abbreviations ac Alternating current AAVM Analytical average-value model ADP Analytical dynamic phasor (model) AVM Average-value model CCM Continuous-conduction mode CP Constant parameter (model) CPU Central processing unit d axis Direct axis DAE Differential-algebraic equation dc Direct current DCM Discontinuous-conduction mode DER Distributed energy resource DP Dynamic phasor emf Electromotive force EMT Electromagnetic transient (program) EMTP Electromagnetic transients program EMTP-type Nodal-analysis-based programs based on the original EMTP FACTS Flexible alternating current transmission systems FAST Frequency-adaptive simulation of transients FDNE Frequency-dependent network equivalent FFT Fast Fourier transform GAM Generalized average modelling HFR High-frequency resonance HVDC High-voltage direct current IM Induction machine LCC Line-commutated converter xx LCR Line-commutated rectifier LST Large-signal-transient (approach) MATE Multi-area Thevenin equivalent MRF Multiple reference frames (model) MSS Multiple-steady-states (approach) ODE Ordinary differential equation PAVM Parametric average-value model PC Personal computer PD Phase-domain (model) PDP Parametric dynamic phasor (model) PMU Phasor measurement unit PWM Pulse width modulation q axis Quadrature axis RLC Resistor-inductor-capacitor (circuit) rms Root mean square SF Switch function SFA Shifted-frequency analysis SM Synchronous machine SMIB Single-machine infinite-bus (system) SSCI Sub-synchronous control interaction SSR Sub-synchronous resonance SV State variable (-based programs) TD Time domain (model) THD Total harmonic distortion TS Transient stability (programs) VBR Voltage-behind-reactance (model) VFD Variable-frequency drive VP Variable parameter (model) xxi Nomenclature In this thesis, scalars are written using italic fonts [e.g., asi ], and vectors and matrices are denoted by bold letters [e.g., abcsi and sR ]. To distinguish between the different representations of electrical quantities (voltages or currents) that are discussed in this thesis, lowercase letters are used to denote time-domain instantaneous signals [e.g., tias ], uppercase letters to denote the SFA-type DPs [e.g., tIas ], hatted uppercase letters to denote analytic signals [e.g., tIasˆ ], and angle brackets “< >” are used to denote the GAM-type DPs [e.g., tikas]. The conventional fundamental-frequency phasors are also used in Chapter 2 but only at the system fundamental frequency, where they are denoted by uppercase letters in italic [e.g., sasI ]. Only basic variables are aggregated in this section; all other variables are defined explicitly throughout the thesis. 0 Null matrix a Phasor rotation operator (by 120 degrees counter clockwise) cos Cosine function xdiag Diagonal matrix containing vector x on the main diagonal ''abcse Subtransient voltage vector ''de d-axis subtransient voltage ''qe q-axis subtransient voltage xfde Scaled field winding voltage (synchronous machine) H Combined machine-load moment of inertia (in s) abcsi Stator/network current vector dci dc current dsi d-axis stator current qsi q-axis stator current xxii J Combined machine-load moment of inertia (in kg∙m2) j Imaginary number sK Park’s transformation matrix DL Inductance of the three-phase branch of the VBR interfacing circuit lfdL Field winding leakage inductance (synchronous machine) .,...,1, NjLlkdj Leakage inductance of the jth d-axis damper winding (synchronous machine) .,...,1, MjLlkqj Leakage inductance of the jth q-axis damper winding (synchronous machine) lrL Rotor winding leakage inductance (induction machine) lsL Stator leakage inductance (synchronous and induction machines) mL Magnetizing inductance (induction machine) mdL d-axis magnetizing inductance (synchronous machine) mqL q-axis magnetizing inductance (synchronous machine) M Transformation matrix of sequence components M Number of q-axis damper windings (synchronous machine) N Number of d-axis damper windings (synchronous machine) P Number of poles Dr Resistance of the three-phase branch of the VBR interfacing circuit fdr Field winding resistance (synchronous machine) .,...,1, Njrkdj Resistance of the jth d-axis damper winding (synchronous machine) .,...,1, Mjrkqj Resistance of the jth q-axis damper winding (synchronous machine) rr Rotor winding resistance (induction machine) sr Stator winding resistance (synchronous and induction machines) s Laplace variable iabcs Current switch function vector vabcs Voltage switch function vector sin Sine function xxiii t Time eT Electromagnetic torque mT Mechanical torque abcv Stator/network voltage vector dcv dc voltage dsv d-axis stator voltage fdv Field winding voltage (synchronous machine) qsv q-axis stator voltage t Simulation step size (unique step size) x 2-norm relative error of x dr d-axis flux linkage of the rotor circuit (induction machine) fd Field winding flux linkage (synchronous machine) .,...,1, Njkdj Flux linkage of the jth d-axis damper winding (synchronous machine) .,...,1, Mjkqj Flux linkage of the jth q-axis damper winding (synchronous machine) qr q-axis flux linkage of the rotor circuit (induction machine) Commutation angle r Electrical rotor position Angular velocity of the arbitrary reference frame (in rad/s) s System fundamental frequency (in rad/s) r Electrical rotor speed (in rad/s) xxiv Acknowledgments I would like to express my sincere gratitude to my research supervisor, Dr. Juri Jatskevich, whose aspiring motivation, exceptional vision and expertise, patient guidance, continuous caring and support were essential to the realization of this research and for my doctoral degree. His unwavering passion for power systems modelling and simulation inspired and motivated me constantly to pursue my research and advance the horizon of engineering science. I could not have imagined having a better supervisor. In addition, I want to acknowledge the financial support from the Natural Science and Engineering Research (NSERC) of Canada under the Discovery Grant entitled “Next Generation Tools for Modelling and Analysis of Evolving Power and Energy Systems” and the Strategic Project Grant entitled “Advanced Integrated AC-DC Systems for Energy Efficient Buildings and Communities in Canada,” both led by Dr. Juri Jatskevich as the principal investigator. I would like to thank the members of my examining committee of the departmental and university exams, Drs. Hermann Dommel, William Dunford, Peter Lehn, Jose Marti, Shahriar Mirabbasi, Farrokh Sassani, and Douw Steyn, for their kind help, valuable discussions and constructive comments, which have contributed to improving this work. I would like to thank all the faculty members, instructors, and professors at UBC whose classes and lectures I attended. Specifically, I would like to thank Drs. Christine Chen, William Dunford, Luis Linares, Ebrahimi Vaahedi, Mr. Nathan Ozog, and my supervisor Dr. Juri Jatskevich for the great teaching experience I gained while working with them as a Graduate Teaching Assistant. I would like to express appreciation to Drs. Hamid Atigechi, Mehrdad Chapariha, Lianghui Dong, Zhenyu Shan, Francis Therrien, and the current graduate students Navid Amiri, Hua Chang, Milad Ebrahimi, and Saeed Rezaee, for their numerous discussions, valuable feedback and xxv comments on my research and the papers written with their co-authorship during my studies at UBC. My special thanks are extended to my friends, fellow colleagues and graduate students from the Power and Energy Systems research group at UBC, for their help and support, and for making my time here more enjoyable: Abdullah Al-Digs, Arash Alimardani, Hamed Ahmadi, Soroush Amini, Patrick Chow, Wonbae Choi, Shuan Dong, Qiang Han, Shreya Iyer, Keven Lei, Mofei Liu, Zhibang Liang, Andrea Marti, Oleksandr Pizniur, Matin Rahmatian, Nana Rong, Melody Ren, Vanness Shen, Arash Tavighi, Javier Tarazona, Yajian Tong, Zhi Qu, Liwei Wang, Xiaotong Wang, Zemeng Wang, Jiayue Xu, Shuai Xu, Tingting Xu, Zejun Yang, Kai Zhang, Qian Zhang, and Sajjad Zadkhast. This list is obviously far from exhaustive. Last but not least, I am deeply indebted to my beloved family and friends for supporting me from the other side of the planet. I am particularly thankful to my parents Xinzhong and Dongmian for raising me in an unconditionally loving environment, and for being the role models of faith, dedication, hard work, and curiosity. Your love is my inspiration and motivation throughout this journey, and for the new adventures ahead. xxvi Dedication To My Wonderful Parents and Beloved Friends 1 Introduction Motivation Modern power systems are evolving quickly, while facing unprecedented technological changes in all aspects of energy generation, transmission, distribution, and consumption. The traditional alternating current (ac) electrical grid, whose aged infrastructures have long been under stressed conditions, is operating close to its physical limits [1]. Direct current (dc) power technology, owing to expanding applications of high-voltage direct current (HVDC) [2] and flexible ac transmission systems (FACTS) [3], is increasingly seen as a viable option for long-distance transmission, interconnection of asynchronous ac grids and renewable generation, effective power-flow control, and power quality support functions. Over the last 20 years, Europe saw the evolution of HVDC systems in voltage levels, converter structures, network topology, control strategies, etc. [4]. More recently, China became the world’s leader in installation of new power infrastructure with over 40 advanced HVDC projects built or launched by 2015 [5]. At the same time, due to economical, ecological, and political factors, renewable energy sources are now supplementing or replacing older fossil-fuel and nuclear power plants at a rapid pace. In 2016, renewables accounted for nearly 62% of the net addition to the global power generating capacity [6]; and as of March 2017, wind and solar, for the first time, exceeded 10% of the total electricity generation in the United States [7]. Unlike traditional power plants, these distributed energy resources (DERs) are usually interfaced with the ac grid using power electronic converters with fast-responding capabilities and extended controllability. Moreover, the advent of smart grid technologies [8] has promoted the proliferation of microgrids [9]-[10] to integrate the increasing 2 Figure 1–1. The ac-dc microgrid system presently installed (solid lines) and to be further developed (dashed lines) in the Kaiser Building at the University of British Columbia, Vancouver, Canada [11]. penetration of DERs by introducing energy storage systems and controllable loads. Figure 1–1 depicts an example of such envisioned integrated ac-dc microgrid power systems [11], which comprises power sources from utility grid and/or a variety of DERs, various types/levels of loads and appliances, some form of energy storages for peak shaving of demands, a backup generation system supporting medium- to long-term interruptions, bi-directional power converter modules, as well as the communication infrastructure that intelligently coordinates and controls all these components. The aforementioned are but a few examples of today’s power system “electronification”, where generation, transmission, distribution and loads are increasingly moved behind the power 3 electronic apparatus. However, reconfiguring the traditional ac grid supplied by reliable generation managed by few operators, into the integrated ac-dc systems laden with intermittent energy sources while involving more participants introduces a high degree of uncertainty and variability [12]. Moreover, numerous challenges arise due to the non-linear properties of power electronic components. For instance, conventional line-commutated rectifier loads are known to contribute into the grid significant amounts of harmonics, which lead to degraded power quality, increased losses, and can cause malfunction of sensitive devices [13]. In addition, these harmonics tend to be aggravated by the oscillatory interactions between the switching operation of power converters, the dynamics of controller systems, and the impedance of the interconnected systems [14]-[16], etc. Unlike large-scale stiff ac grids, power-electronic-based systems may provoke oscillatory interactions over a broad range of frequencies, including the sub-synchronous resonances (SSR) frequently seen in weak or standalone grids [14], the sub-synchronous control interactions (SSCI) between wind farms and HVDC systems [15], and some super-synchronous high-frequency resonances (HFR) due to the inductive/capacitive behaviour of converters [16], etc. To investigate these phenomena, researchers and engineers worldwide are dependent on a massive number of computer simulations. Accurate and efficient modelling of each system component is therefore critical in all stages of design, analysis, monitoring, control, and energy management to ensure the secure, reliable and optimized operation of the integrated ac-dc power systems. The study of power system transients has been an active area of research for many decades, and various analysis tools and simulators have been developed for different objectives, time scales, types of disturbances, etc. [17]-[19]. Figure 1–2 depicts the time frames of transient phenomena frequently seen in power systems, which covers a range from a fraction of a microsecond (lightning and wave propagation) to hours and days (mid- and long-term dynamics). In particular, it is seen in Figure 1–2 that despite the various causes of disturbance, the power system transients can be generally classified into electromagnetic transients and electromechanical transients. This 4 Figure 1–2. Time frame of different power system transient phenomena. classification thus leads to the well-known “time-scale modelling” of power systems [18], which first describes the complete power system dynamics using several sets of ordinary differential equations (ODEs) or differential-algebraic equations (DAEs) with distinct ranges of time constants, and then, by focusing on the (electromagnetic transient or electromechanical transient) phenomena of interest, selectively models the desired components using proper signal representations based on assumed approximations, while neglecting/simplifying the rest of the system. The electromagnetic transients may be induced by a perturbation from external causes (e.g., lightning), and/or a change in the physical network configuration (e.g., action of breaker, operation of power electronic devices, equipment failures, faults, etc.). To study these high-frequency and relatively-fast phenomena, the so-called electromagnetic transient (EMT) programs are used, which can generally be divided into two categories: the nodal-analysis-based tools derived from the original electromagnetic transient program EMTP introduced by H. W. Dommel [20] 5 (hereinafter referred to as EMTP-type [20]−[26]), and the state-variable (SV-based) tools [27]−[32]. In EMTP-type programs, the system circuit/components are first discretized at the branch level (typically using the trapezoidal rule), and then combined into a system nodal equation [20] which is then solved as a difference equation at each time step. In SV-based tools, the system components are represented in state-space formulation using ODEs (or DAEs), which are then combined using state model generation algorithms [33]-[34] into a first-order system of ODEs (or DAEs) and solved at the system level using fixed- or variable-step integration rules [27]−[28]. In both categories of the available state-of-the-art EMT simulators, power system components are modelled in detail and represented using instantaneous time-domain signals. This endows EMT programs with high modelling fidelity, thus making them suitable for studies such as insulation coordination, protection scheme calculation, detailed short-circuit analysis, power electronics controller design, etc. [20]. However, EMT simulators typically require very small step sizes (ranging from a few to hundreds of microseconds) to accommodate for the fast transients, thus yielding significant computational cost/simulation time, and limited size of the systems that can be practically simulated [35], [40]. For studying the slower electromechanical dynamics due to the mismatch of the energy stored in the rotating machine masses and the electrical network, the transient stability (TS) programs [36]−[39] are commonly used. Therein, the stator-network transients are generally neglected [19], and the electric quantities are represented as conventional fundamental-frequency phasors. In TS programs, the complete power system is mathematically modelled by a set of DAEs, where the differential equations represent dynamics of rotating machines (and their excitations, turbine-governor, etc.) and the algebraic equations describe the power network, loads, connected devices, etc. As such, the TS programs are typically designed to handle the (multi-machine) system transients characterized by mid- to long-term dynamics (e.g., low-frequency oscillations from 1 to 3 Hz [17]). Moreover, the TS programs usually use additional modelling approximation such as the positive-sequence single-phase reduced-order representation [18]−[19], which yield less accurate results but permit the use of larger step sizes (up to several milliseconds). It is also known 6 that the TS programs are incapable of accurately representing nonlinear elements in the network (e.g., HVDC/FACTS applications, power electronic apparatuses, etc.), or capturing detailed fast dynamic events (e.g., line energization) [40]. Despite the different time scales of power system transient phenomena, there is not a clear distinction between the electromagnetic or electromechanical transients, as the power system components may change conditions with time and interact with each other. Furthermore, it is noted in Figure 1–2 that state-of-the-art EMT and TS programs, as confined to their distinct time frames, are not well suited for the transient phenomena associated with integrated ac-dc power systems, which may span the range from high-frequency transients (hundreds of hertz to several kilohertz) to sub-synchronous frequencies (a few to tens of hertz). To overcome this limitation, the so-called hybrid simulation approaches have been proposed in the literature to interface the EMT and TS programs [40], by splitting the network into subsystems and then using either an EMT or TS solution for each subsystem, based on the required simulation fidelity. To facilitate the simulation process, efforts were made toward simplifying the system models and/or obtaining subsystems equivalents, such as multi-area Thevenin equivalent (MATE) [41]−[42], frequency-dependent network equivalent (FDNE) [43], or making use of high-performance computing techniques such as parallel processing [44], etc. There are also programs that can switch between EMT and TS simulations by specifying a time or condition, with the use of proper component models at the corresponding stages [39]. Despite various attempts, challenges with using these hybrid simulation approaches [39]-[46] include the presence of numerical discontinuities when switching between simulation methods [45], and the interfacing between the various subsystems and their solutions [40], [46]. Consequently, there is a significant need for accurate and efficient modelling and simulation techniques that will make the simulations as detailed as EMT and as fast as TS simulators. New-generation tools with such capabilities will be particularly useful for general-purpose transient simulations of the integrated ac-dc power systems. 7 Literature Review Dynamic Phasor (DP) Type Modelling Approaches Instead of tracking the detailed but slow-evolving time-domain waveforms (as in EMT programs), or dealing with reduced-order fundamental-frequency phasor representations (as in TS programs), this thesis considers modelling and simulating general ac-dc power systems in hybrid time-phasor signal representations. These conceptually similar approaches are generally referred to in the literature as dynamic phasor (DP) type modelling approaches. Both types of DP theories, namely the shifted-frequency analysis (SFA) [50]-[51] and the generalized averaging method (GAM) [62] are considered in this work. In both types of DP-type simulations, the power systems are modelled using analogous time-phasor signals, which possess the following advantages/properties: 1) good accuracy in capturing both the electromagnetic and electromechanical transients due to the detailed (full-order) modelling of power system components, while avoiding the numerical errors that would be present from interfacing EMT and TS simulators; and 2) flexible selection of time-step sizes (by variable-step solvers) which permits using large time steps in steady-states and during electromechanical trainsets [47]-[48], thus achieving superior efficiency than the time-domain EMT simulations typically requiring small time steps. 1.2.1.1 Shifted-Frequency Analysis (SFA) The first definition of DP, which stems from signal processing theory, was originally applied to power system simulations in [49]. This type of DP represents the envelope of analytic (complex) signals [50], where the frequency spectra of power system signals (around the 50/60 Hz fundamental frequency) can be down-shifted (to around 0 Hz) to enable flexible selection of step 8 sizes during simulations, thus referred to as the shifted-frequency analysis (SFA) [51]-[52]. The concept of SFA-type DPs is first introduced in [50] (therein defined as the envelopes of analytic signals), and then systematically described in [51]-[52], where several numerical aspects are discussed, including signal representation, integration rule, simulation program structure, and models of general system components. However, the generator model in [49] is expressed in reduced-order form (two-axis model, [18]) and only the fast dynamics of simple RLC networks are preserved. In [50], the (constant-parameter) transmission line model is included, and the synchronous machine is modelled in full-order using analytic signals, which yields an indirect interface of qd0-coordinate model to external abc-phase networks that may result in numerical instability. To overcome this limitation, following efforts were made toward direct machine-network-interfaced modelling of rotating machines, by using the so-called voltage-behind-reactance (VBR) formulation [53]-[54] (referred to as VBR-SFA models). Nevertheless, the established synchronous machine model [53] has variable-parameter (rotor-position-dependent) inductance/conductance matrices, which greatly increase computation cost and complicity of implementation. Furthermore, as modelled for EMTP-type solutions, these conductance matrices (after discretization) are formulated as a function of the step size [53]-[54]. Accordingly, varying the step size to reflect the active modes at a given point of simulation (as is particularly useful for DPs) requires re-discretizing all components and reformulating/re-factorizing the conductance matrix, which gives rise to prohibitive cost computations. Other recent work on SFA-type DP simulations include studies validating the SSR transient phenomena [51], the inrush dynamics during energization of transformers [55], and potential application for increasing speed of large-scale power system simulations [26]. A similar concept to the SFA-type DPs, namely the frequency-adaptive simulation of transients (FAST) [56]-[57], has also been developed and shown to effectively simulate multi-scale transients with the carrier frequency shifting techniques. In this method, power system electromechanical transients are simulated by shifting the carrier frequency by system fundamental 9 frequency (i.e., using SFA-type DPs) [56]; and for electromagnetic transients, the carrier shift frequency is set to zero (i.e., using instantaneous time-domain signals). Recent application of the FAST method include modelling of rotating machines [58]-[59], transmission lines [60], and power-electronic-interfaced systems [61]. Nevertheless, it is noted that due to the frequency shift of only one frequency, the SFA-type DPs are mainly applicable to modelling the power system transients where the fundamental frequency is dominant. 1.2.1.2 Generalized Averaging Method (GAM) To account for high-order harmonics that may exist in ac-dc power systems, another type of DP is proposed based on the generalized averaging method (GAM) [62], and defined as the time-varying Fourier coefficients of “sliding-window” waveforms. In this theory, it is observed (and assumed) that the time-domain signals in power systems attain near periodical waveforms, whose Fourier coefficients thus become constant or slow-varying, and can be handled separately with many beneficial numerical features (continuous, time-invariant, etc.) [63]. Moreover, the GAM allows the user to retain the desired frequency components of power system signals (i.e., flexible modelling accuracy), and can be readily augmented to include higher-order harmonics of interest, if so needed. Due to these desirable properties, the GAM-type DPs have been comprehensively applied to power system simulations [63]−[65], component modelling of rotating machines [63]−[68], FACTS devices [69]−[72], HVDC and high-power converter systems [73]−[78], power electronic loads/devices [79]−[80], etc., as well as case studies on unbalanced fault conditions [67]-[68], [81], SSR and SSCI phenomena [70],[82],[83]. Recently, there has also been a growing interest in extending the GAM-type DPs to simulating medium-/large-scale power systems in off-line [21] and real-time [25], [32] industry-grade EMT simulators. However, one major limitation of using GAM-type DPs [73]−[78] is the use of the switch functions that relate the ac- and dc-side dynamics when modelling power converters. In [73]−[76], the switch functions are developed to emulate the switching behaviour of each diode, and are expressed as piecewise-linear-approximated waveforms, thus rendering difficulty to derive 10 compact mathematical expressions for high-order DPs. Accordingly, only the fundamental frequency dynamics of ac systems are preserved in [73]−[74], where higher-order harmonics are truncated. This issue was later addressed by using the Fourier-series-approximated switch functions [77], which requires costly computations per time step. Moreover, all these established analytical DP models [73]−[77] have been analytically derived considering only the common continuous-conduction mode (CCM) of operation (hereafter categorized into ADP-1 [73]−[76] and ADP-2 [77] models based on their use of piecewise-linear- and Fourier-series-approximated switching functions, respectively). This results in switch functions with predetermined fixed-shape waveforms that are not accurate in the other operating modes, thus limiting the model accuracy in predicting actual currents and voltages. Another recent attempt proposes to analytically model the ac-side DP dynamics in qd reference frame (namely ADP-qd model [78] ), which nevertheless considers only the fundamental ac components and the CCM mode of operation limited to diode rectifiers. Numerically Efficient Modelling of AC-DC Power System Components In addition to modelling approaches, the simulation performance of the integrated ac-dc power systems is highly dependent on the proper modelling of each system component. A representative unit of an ac-dc energy conversion system is the backup generation unit shown in Figure 1–1, which includes the rotating electric machines and some kinds of power electronic switching devices, such as the line-commutated rectifiers (LCRs). Such machine-rectifier systems are commonly used in many practical power systems and ac-dc energy conversion applications. To study their dynamic behaviour, thousands of researchers and engineers worldwide run extensive computer simulations everyday, where higher simulation speed and accuracy can make a significant difference and impact. Models of electric machines and power converters are typically the bottleneck in most simulation programs [20]−[32], and therefore are considered to be the focus of this thesis. 11 1.2.2.1 Modelling of Electric Machines Synchronous machines are present in almost every bulk/islanded power system as the main source of electrical energy, as well as being the functional unit for power stability or reactive power support [19], [84]. Induction machines account for 60 to 70% of today’s electrical energy consumption [19], and are also prevalent in wind generation [85]. To accurately represent the dynamics of these rotating machines, full-order general-purpose lumped-parameter models, based on magnetically coupled circuits of machine’s stator and rotor windings [84], are widely used in power system EMT simulations. Depending on the interfacing circuit with the power network, rotating machine models can be generally classified into three types: coupled-circuit phase-domain (PD) models [84]; classical qd models [19], [84]; and hybrid models known as the VBR models [88]-[97]. Despite their algebraic equivalence, these models can possess vastly different numerical properties, thus affecting numerical accuracy, efficiency, and stability during simulation [86]. The PD models [84], as originally derived based on physical variables, use the same set of coordinates (i.e., abc phases) that can be directly interfaced with external power systems. However, the mutual inductances (and self-inductances for salient-pole synchronous machines) of PD models are dependent on the rotor position. These time-varying inductances require costly calculations per time step, and greatly complicate the analysis of machine dynamics due to poorly-scaled eigenvalues [84], [86]. R. H. Park [87] proposed to transform the stator windings of the PD model to the orthogonal quadrature and direct axes fixed on the rotor (thus referred to as qd models [19], [84]). In qd models, the circuit parameters are time-invariant and the state variables (voltages, currents, fluxes etc.) become constant in steady-states, which make these models highly efficient. As such, qd models become the classical default (built-in) models in most state-of-the-art simulators even today. However, the qd model can cause a major issue when interfaced with external inductive 12 abc-phase networks (e.g., transformers, transmission lines, etc.). In SV-based programs, this issue is presented as the incompatibility between the qd models, represented as voltage-input current-output components, and the external inductive systems that require inputs of voltages [86]. A common approach to solve this problem is to use artificial (relatively-large resistive or capacitive) snubber circuits [30], which however add to system stiffness decreasing simulation efficiency (forcing the use of small step sizes), and can even cause continued numerical instability [86]. In EMTP-type simulators, this interfacing issue also exists by requiring prediction of machine state variables (speed voltage, stator current, field voltage, etc.), thus leading to limited accuracy and requiring small time-steps [86]. To achieve direct machine-network interface while attaining good numerical efficiency, the VBR models [88]-[97] have been proposed with a circuit configuration following the widely established machine representation in the power system community − an equivalent voltage source behind a reactance [18]−[19]. The VBR modelling of synchronous machines has been first proposed in [88], where the stator states (currents) are in abc phases, and the rotor currents/fluxes are formulated in qd coordinates. This unique formulation endows the VBR models with better-scaled eigenvalues, and thus, higher efficiency than the PD models. However, due to the dynamic saliency of synchronous machines [89], the interfacing circuit of model [88] still contains rotor-position-dependent variable-parameter inductances (referred to as the VP-VBR model [88]). Consequently, a considerable amount of following research has been done toward eliminating such rotor-position-dependency of stator inductances, by adding a fictitious high-frequency damper winding [89]−[91], or by re-positioning the time-varying part into the subtransient voltage source (followed by numerical approximations) [92]−[93]. In particular, it is shown in [90]−[93] that with a constant-parameter stator-network interface (referred to as CP-VBR model [90]−[93]), the VBR modelling of synchronous machines can achieve a superior combination of numerical efficiency and accuracy. This allows the VBR modelling to be successfully applied to induction machines [94], and extended for EMTP-type solutions [95]-[97]. Owing to the numerical advantages, the VBR formulation of rotating machines has been used for SFA-based DP modelling [53]-[54]. 13 1.2.2.2 Modelling of Line-Commutated Rectifiers (LCRs) The second category of power system components considered in this thesis are the high-power ac-dc converters. While pulse-width-modulated (PWM) converters are commonly found, their applications are usually limited by the power ratings of transistors [98]. Owing to high reliability, low cost, and simple configuration, line-commutated rectifiers (LCRs) or converters (LCCs) are widely used in a myriad of industrial and commercial applications [99], including the input stage of variable frequency drives (VFDs) [98]-[99], DERs [100], electric systems of vehicles, ships and aircrafts [101], and over 70% of today’s HVDC systems [4]-[5], etc. However, due to the automatic line current commutation, the conventional LCR applications, e.g., three-phase (six-pulse) diode rectifiers, may work in various operating modes (switching patterns), including several CCM modes (CCM-1/2/3, [115]) and the discontinuous-conduction mode (DCM) [98], [115]-[116]. These operating modes may contribute considerable ac harmonics into the systems, resulting in distorted voltages, degraded power quality, and even adverse impacts on other equipment [102]. Numerically accurate and efficient models of LCRs are therefore required to predict system-level transients as well as the varying harmonics under different loading conditions. The traditional detailed switch-level models of LCRs can be readily established using the built-in library of industry-grade EMT simulation tools [21]-[32]. In detailed models, the switching of all semiconductor devices is fully included (i.e., high modelling fidelity), which however requires small integration step sizes (ranging from tens to hundreds of microseconds) and can result in prohibitively long simulation time, especially for studies where multiple runs of simulations are required for controller design and tuning [103] . Alternatively, a variety of modelling approaches of ac-dc power converters have been proposed, including exact linear time varying modelling [104]–[105], sampled-data modelling [106], and the average models of various kinds [107]-[108]. In particular, the so-called dynamic average-value models (AVMs) [84], [109]–[118], that neglect/average the fast switching of time-domain dynamic variables (voltages, currents, etc.) over a prototypical switching interval [84], are 14 shown to significantly outperform the detailed counterparts in terms of required simulation time. Moreover, due to the time-invariant formulation, AVMs can be linearized around any desired operating point, and thus are suitable for small-signal analysis (e.g., deriving local transfer functions). The AVMs of LCRs can be mainly classified into two categories: analytical AVMs (AAVMs) [84], [109]-[110], and parametric AVMs (PAVMs) [111] – [118]. In AAVMs, the relationships between the averaged ac and dc variables are derived analytically based on a specific switching pattern of diodes/thyristors, thus are valid for only that specific operating mode (i.e., CCM-1 as defined in [115]). However, in practice the LCR circuits may become complicated including non-idealities, losses and parasitic elements, such that the system dynamics vary frequently with operating modes. Due to these limitations, the PAVMs of diode-rectifier systems [111]–[116] have been proposed to facilitate the construction of AVMs, which numerically relate the rectifier/dc-link dynamics through a set of algebraic parametric functions. The PAVM methodology has been proven highly effective for simulations over a wide range of operating conditions, and has been extended to (machine-fed) thyristor-controlled rectifiers [117]. Nevertheless, the established AVMs of LCRs in [109]–[117] have only considered the fundamental components of ac variables. To overcome this limitation, a recent work proposed to include the main ac harmonics by utilizing multiple reference frames rotating at different speeds/directions (referred to as MRF-PAVM [118]). However, this approach formulates the ac variables in several qd coordinates. As discussed in Section 1.2.2.1, models developed in qd coordinates [109]−[118] also require special consideration for their interfacing with external networks represented in physical abc phase variables [86]. Moreover, since double-fundamental frequencies may exist in qd coordinates under faulty or unbalanced conditions [63], the simulation efficiency with these models can be compromised. Therefore, the GAM methodology, as the “generalization” to time-domain AVMs, has attracted a renewed attention with the goal to develop LCR models [73]−[78] that possess direct interfacing circuit and include effects of ac harmonics, as discussed in Section 1.2.1.2. 15 Research Objectives and Anticipated Impact The properties of state-of-the-art rotating machine models discussed in Section 1.2.1.1 and 1.2.2.1 are summarized in Table 1–1. It is recalled that for modelling of rotating machines, indirect interfacing with external abc-phase networks can add to system stiffness, resulting in degraded simulation efficiency or even numerical instabilities. In addition, the machine-network interfacing circuit, if containing time-varying (rotor-position-dependent) parameters, will greatly increase the computation cost and complicate model implementation. Moreover, to take advantage of flexible time-step sizes (i.e., using variable-step solvers), it is desirable to derive DP models of rotating machines in state space formulation (for SV-based solution). As can be seen in Table 1–1, there is no state-of-the-art DP model of synchronous or induction machines that possess all the following desirable properties: direct abc interface, constant-parameter interfacing circuit, and state-space formulation for the SV-based solution. However, as seen in Table 1–1 (last two rows), the new models developed in this thesis fulfill this gap in research and advance the state of the art. Table 1–1. Properties of state-of-the-art and proposed electric machine models. Models Signal Types Machine Type Direct abc Interface Interfacing Circuit Parameters Solution Type PD [84] Time-Domain Signals SM/IM Direct Variable SV-based qd [19], [84] Time-Domain Signals SM/IM Indirect Constant SV-based VP-VBR [88] Time-Domain Signals SM Direct Variable SV-based CP-VBR [89]−[94] Time-Domain Signals SM/IM Direct Constant SV-based VP-VBR [95] Time-Domain Signals SM Direct Variable EMTP-type CP-VBR [96]−[97] Time-Domain Signals SM/IM Direct Constant EMTP-type qd-SFA* [50] Analytic Signals SM Indirect Constant EMTP-type VBR-SFA [53] SFA-Type DPs SM Direct Variable EMTP-type VBR-SFA [54] SFA-Type DPs IM Direct Constant EMTP-type 16 Models Signal Types Machine Type Direct abc Interface Interfacing Circuit Parameters Solution Type Proposed VBR-SFA [Sect.3.1] SFA-Type DPs SM Direct Constant SV-based Proposed VBR-SFA [Sect.3.2] SFA-Type DPs IM Direct Constant SV-based * The model [50] is derived using analytic signals, while the SFA concept is applied in the integration rule. Similarly, the properties of state-of-the-art models of LCR systems discussed in Section 1.2.1.2 and 1.2.2.2 are summarized in Table 1–2. Therein, in addition to the interfacing issue with external networks, it is recalled that the main challenges for accurate modelling of LCRs are: 1) the representation of various operating modes (i.e., DCM and CCM-1/2/3 [115]), and 2) the inclusion of main ac harmonics. Moreover, it is desirable to include the model of thyristor-controlled LCR systems. Finally, it is recalled in Section 1.2.2.2 that the use of different modelling techniques (detailed, analytical, or parametric) can result in distinct formulations and different numerical properties (accuracy, efficiency, etc.) of the models. As seen in Table 1–2, there is no state-of-the-art model of LCR systems that uses DPs in abc–phase coordinates, covers all operation modes, includes harmonics, and predicts thyristor operation. However, as seen in Table 1–2 (last two rows), the new models developed in this thesis combine the desired properties and advance the state of the art. Table 1–2. Properties of state-of-the-art and proposed models of LCR systems. Models Signal Types Modelling Technique Direct abc Interface Operation Modes Harmonics Inclusion Thyristor Operation Detailed* [21]-[32] Time-Domain Signals Detailed Direct CCM-1/2/3, DCM Yes Yes AAVM [109]-[110] Time-Domain Signals Analytical Indirect CCM-1 No Yes PAVM [111]-[116] Time-Domain Signals Parametric Indirect CCM-1/2/3, DCM No No 17 Models Signal Types Modelling Technique Direct abc Interface Operation Modes Harmonics Inclusion Thyristor Operation PAVM [117] Time-Domain Signals Parametric Indirect CCM-1/2/3, DCM No Yes MRF-PAVM [118] Time-Domain Signals Parametric Indirect CCM-1/2/3, DCM Yes No ADP-1 [73]-[74] GAM-Type DPs Analytical Direct CCM-1 No Yes ADP-1 [75]-[76] GAM-Type DPs Analytical Direct CCM-1 Yes Yes ADP-2 [77] GAM-Type DPs Analytical Direct CCM-1 Yes Yes ADP-qd [78] GAM-Type DPs Analytical Indirect CCM-1 No No Proposed PDP [Sect. 4.1] GAM-Type DPs Parametric Direct CCM-1/2/3, DCM Yes No Proposed PDP [Ch. 4.2] GAM-Type DPs Parametric Direct CCM-1/2/3, DCM Yes Yes * Despite accuracy and relative simplicity, the detailed models are known for being numerically expensive due to switching/commutation of all semiconductor devices. The following objectives are specifically considered to advance the goal of this research: Objective 1: Summarize the fundamentals of state-of-the-art DP-type modelling approaches. Pinpoint the differences and limitations through analytical and simulation studies, and discuss the notable numerical features associated with DP-type simulations. Objective 2: Develop an SFA-type DP synchronous machine model with direct stator-network interface and constant-parameter interfacing circuit that is suitable for SV-based EMT simulators. Objective 3: Extend Objective 2 to developing induction machine models in state-space formulation with improved numerical efficiency. 18 Objective 4: Develop a GAM-type DP model of diode LCR systems that possesses abc–phase interfacing circuit, accounts for harmonic dynamics, and covers wide range operating modes. Objective 5: Extend Objective 4 to numerically efficient modelling of thyristor-controlled rectifier systems. Objective 6: Develop an interface between the SFA- and GAM-type DPs to connect the proposed DP models. The ultimate goal of this thesis is to increase the numerical efficiency of DP-type simulations of integrated ac-dc power systems, which is achieved by proposing several new DP models with improved numerical properties and new desirable features. Rigorous case studies demonstrate the advantageous numerical efficiency (i.e., less required CPU simulation time) and accuracy (i.e., lower relative error of predicted trajectory against the reference solution) of the proposed DP models. The enhanced numerical efficiency of the proposed DP models will greatly accelerate the simulations of transient studies, which is highly desirable for simulating dynamics of mid- to large-scale ac-dc power systems, or even for real-time simulations. Moreover, with the advantageous accuracy in capturing the desired (electromagnetic and/or electromechanical) phenomena, it is expected that the proposed DP models will be applied to studying existing and/or new transient phenomena (as shown in Figure 1–2) occurring in the emerging integrated ac-dc power systems. Finally, as a general-purpose simulation methodology to bridge the gap between the traditional EMT and TS simulators, the proposed DP-type modelling and simulation techniques will become highly useful to many researchers and engineers worldwide, and facilitate the development of next-generation power system simulation tools. 19 Dynamic Phasor Modelling Approaches for Power System Transient Simulations To set the stage for the derivation of the proposed models, this chapter presents the fundamentals of the DP-type modelling approaches. We begin by briefly reviewing the conventional signal representations in power system simulations. Next, the fundamentals of DP-type modelling approaches are set forth, including the shifted frequency analysis (SFA) and the generalized averaging method (GAM) types of DP theories. Moreover, for demonstration and comparison, we present several case studies of DP modelling for basic power system components and phenomena, while discussing the differences and limitations of the approaches. Conventional Signal Representations in Power System Simulations Dynamic behaviour of power systems can be studied using various transient simulation programs developed for different objectives, time scales, types of disturbances, etc. [17]-[19]. Without loss of generality, here the state-space representation of power system models (which is used in SV-based tools [27]−[32]) is briefly discussed to illustrate the fundamentals of general-purpose power system transient simulation. Further details and specific aspects can be found in literature [17]-[19]. The EMTP-type modelling of power systems (for EMT simulations) can also be found in numerous references [20]−[25]. Generally, the dynamics of power systems can be represented by a set of first-order ODEs in the state-space formulation as: 20 tdtd,,uxfx , (2–1) t,,uxgy . (2–2) Here, x is the state vector, u is the input (or control) vector, t is time, and y is the output vector. In some cases, it is also desirable to include algebraic constraints, thus transforming the set of ODEs into a system of DAEs [33]-[34]. For dynamic systems that are linear (and/or piecewise-linear) and time-invariant, (2–1) and (2–2) can be expressed in the well-known matrix-form state-space representation as BuAxx dtd, (2–3) DuCxy , (2–4) where A , B , C , and D are the state-space matrices. The eigenvalues of matrix A represent the system’s dynamic modes [19]. Equations (2–1)–(2–2) [or (2–3)–(2–4)] can be solved using various numerical integration methods such as MATLAB/Simulink’s ODE solvers [27]–[28]. Specifically, depending on the problem type (continuous/discrete states, stiff/non-stiff) and requirements for accuracy/efficiency, etc., the appropriate ODE solvers [27]–[28] may be selected with features such as fixed- or variable-step, explicit or implicit, etc. In general, for transient simulation of integrated ac-dc power systems with dynamic modes spanning a wide range of time frames (i.e., stiff problem), the variable-step implicit solvers (e.g., ode23tb or ode15s [27]–[28]) may offer good numerical accuracy and efficiency. In addition to the simulation framework and numerical integration methods, another critical factor affecting the simulation performance is the selection of proper signal representation for modelling of power system components. 21 Time-Domain Instantaneous Signals In industry-grade EMT simulators (both EMTP-type and SV-based categories), power system components are typically modelled in detail and represented using instantaneous time-domain signals. For illustration, the general ac power system electrical variable (voltage or current) in steady state can be expressed by a sinusoidal signal tu as ttu scosU , (2–5) where s denotes the system fundamental frequency (hereafter assumed to be 3772 ss f rad/s, where 60sf Hz), and U and are the signal magnitude and phase angle, respectively. 2.1.1.1 abc-Phase Coordinates For commonly used three-phase (a-b-c) power systems, the signals in abc-phase coordinates can be expressed as cscbsbasacbaabcttttutututcosUcosUcosUu , (2–6) where the abc-phase signals may attain arbitrary magnitudes and phase offsets. In particular, when the system is balanced and in positive sequence, it satisfies .3232,UUU cbacba (2–7) However, condition (2–7) may not always hold. Under unbalanced conditions, the asymmetrical abc-phase signals in (2–6) can be decomposed into three sets of amplitude-invariant symmetrical components [19], as 22 tttt zabcnabcpabcabc ,,, uuuu , (2–8) where 32cos32coscosU,pspspsppabcttttu , (2–9) 32cos32coscosU,nsnsnsnnabcttttu , (2–10) zszszszzabcttttcoscoscosU,u . (2–11) Here, pU , nU and zU denote the magnitudes of the positive-, negative- and zero-sequence components, respectively; and p , n and z denote the phase angles of the positive-, negative- and zero-sequence components, respectively. 2.1.1.2 qd0 Coordinates The main challenge with using abc-phase representation for power system signals is the presence of rotor-position/speed-dependent (i.e., time-varying) parameters when modelling electric machines [84]. To simplify the analysis of power systems, the well-known Park’s transformation [87] has been widely used, which, in effect, replaces the ac variables associated with abc-phase machine stator windings or stationary circuits, with dc-like signals associated with fictitious windings/circuits rotating at an arbitrary angular velocity, i.e., in the qd0 reference frame. The transformations of power system signals between abc-phase and qd0 coordinates are given by tt abcsqd uKu 0 , (2–12) tt qdsabc 01uKu . (2–13) 23 Here, the transformation matrices are defined as [84] 212121)32sin()32sin(sin)32cos()32cos(cos32sK , (2–14) 1)32sin()32cos(1)32sin()32cos(1sincos1sK , (2–15) where denotes the angular displacement of the arbitrary reference frame. For instance, for transient and dynamic studies of large-scale power systems, the variables of system components (except for synchronous machines) are commonly referred to the reference frame rotating at the synchronous speed [18], thus tsdt00 , (2–16) where 0 is often set as zero for simplicity. Based on (2–16), the time derivatives of the transformation matrices are given by )( sssdtdKJK , (2–17) )( 11 JKK sssdtd , (2–18) where 000001010J . (2–19) Moreover, for unbalanced systems, it is noted that applying (2–14) to (2–8) yields 24 zsznsnsnpppqdttttcos00U02sin2cosU0sincosU0u , (2–20) Conventional Fundamental-Frequency Phasors For studying the electromechanical dynamics of power systems (in TS programs), the conventional fundamental-frequency phasor is widely used based on the assumption of neglecting stator-network transients [19]. For the steady state power system signal tu shown in (2–5), its fundamental-frequency phasor is defined as UsinUcosUU jeU js , (2–21) where sU is assumed to rotate at the constant fundamental frequency s with peak (not rms) magnitude U and phase shift . The time-domain signal is retrieved by ]Re[ tjs seUtu . (2–22) The fundamental-frequency phasor possesses the advantage of using steady-state (i.e., algebraic) relationships for representing the machine stator and the interconnecting transmission network, thus making it convenient for power flow and short circuit studies [18]. It is also noted that for unsymmetrical fault analysis, the unbalanced power systems can be treated as an algebraic combination of sequential networks [19]. In particular, the time-domain equation (2–8) can be transformed into fundamental-frequency phasors as spnzszsnspscsbsasabcUUUUUU MUMU , (2–23) 25 where 3222 ,11111jeaaaaa M , (2–24) and 1111131221aaaaM . (2–25) Other definitions of matrix M are also possible, e.g., with the scaling factor of 31 for power-invariant form [67]-[68], which requires adjustment in the inverse transformation. DP-Type Modelling Approaches Regardless of the signal representations in Section 2.1, it is noted that there exist certain limitations to their application in power system simulations. For time-domain instantaneous signals, the major drawback is the maximum integration step size limited by the signal bandwidth (i.e., maxmax 21 ft , [120]). Specifically, abc-phase signals attain a band-pass bandwidth around the fundamental frequency [as noted in (2–8)]. In qd0 coordinates, under unbalanced conditions, the signal bandwidth is expanded due to the negative-sequence components at double the fundamental frequency [as noted in (2–20)]. In contrast, the conventional fundamental-frequency phasors, by neglecting stator-network transients of power systems (i.e., reduced-order representation [18]−[19]), simply cannot predict the fast dynamic phenomena associated with integrated ac-dc power systems. To overcome these challenges, instead of tracking the detailed but slow-evolving time-domain waveforms, or dealing with reduced-order fundamental-frequency phasors, this thesis considers modelling general ac-dc power systems using hybrid time-phasor signal representations, i.e., using the SFA- and GAM-type DP approaches. 26 Shifted-Frequency Analysis (SFA) The concept of SFA-type DPs stems from the modulation techniques in signal processing theory, by investigating the analytic representation of band-pass signals [50]. Typically, for the power system signal tu with a frequency spectrum as depicted in Figure 2–1, it can be represented as a band-pass signal formed by the sum of closely-spaced sine waves, as iisi tiatu coslim0. (2–26) Equation (2–26) can then be written as ttuttutu sQsI sincos , (2–27) where iiiI tiatu coslim0, (2–28) iiiQ tiatu sinlim0. (2–29) The instantaneous band-pass signal tu can thus be considered as the combination of two low-pass signals tuI and tuQ , modulated by the carrier signals tscos and tssin , respectively. Since the two carrier signals have the same frequency and are out of phase by 90°, tuI and tuQ are often referred to as the in-phase and quadrature components of tu , respectively. The SFA-type DP of the instantaneous signal tu is then defined as [49]-[50] tjututU QI . (2–30) In signal processing, tU is also called the complex envelope of the signal tu , which excludes the fundamental frequency carrier but preserves the neighboring frequencies. 27 Another commonly used representation of the instantaneous signal tu is the analytic signal [50], [119], defined as tjssQIsQsIsQsIsetUtjttjututtuttujttuttutujHtutUsincoscossinsincosˆ (2–31) where H denotes the Hilbert transform defined as [119] 1dttutuH . (2–32) Manipulating (2–31), the relationship between the instantaneous signal tu and its SFA-type DP tU can be obtained as tjtj ss etujHtuetUtU ˆ , (2–33) tj setUtUtu ReˆRe . (2–34) As seen in (2–33) and (2–34), the SFA-type DP tU can be considered as the analytic signal tUˆ after frequency modulation, i.e., shifted by s from the original frequency. The instantaneous signal tu , as retained in the real part of the analytic signal tUˆ , can also be readily retrieved from the DP tU . This process is commonly referred to as the shifted-frequency analysis (SFA) [51]-[52], and can be visualized in Figure 2–2, where it is seen that the SFA-type DP is obtained as an analogous low-pass representation of the original band-pass instantaneous signal. 28 Figure 2–1. Frequency spectrum of a band-pass power system signal tu represented by the sum of closely spaced sine waves. Figure 2–2. Frequency spectra of the band-pass signal in different representations during the process of SFA: a) instantaneous signal tu ; b) analytic signal tUˆ ; and c) SFA-type DP tU . Generalized Averaging Method (GAM) The SFA assumes that power system signals are within a band-pass spectrum, i.e., with frequency content condensed around the fundamental frequency. However, this low-pass signal representation (and its resulting numerical advantages) can be compromised when considering power systems where harmonics may exist, since the frequency shift applied to the analytic signal is limited to only s [as noted in (2–33) and Figure 2–2]. This limitation thus necessitates a new modelling theory that can account for wide-band signals with high-frequency components. In the considered ac-dc power system, the power system signal tu is assumed to be quasi-periodic and possess a wide-band frequency spectrum as depicted in Figure 2–3 (a), which comprises the fundamental frequency and potential high-order harmonic components. The GAM 29 then views the time-domain waveform as u over a “sliding window” with period ssT 2 as depicted in Figure 2–4, which can be expressed using Fourier series [62] as kjkksetuu , (2–35) where tTt s , . Here, since u may not attain periodicity (as shown in Figure 2–4), the kth Fourier coefficient tuk is time-varying, and therefore defined as the kth order GAM-type DP, which is given by deuTtu ssjktTtsk 1. (2–36) Similar to SFA, in frequency domain the Fourier coefficient, i.e., GAM-type DP tuk, can be viewed as the low-pass representation of the instantaneous signal tu , which is visualized in Figure 2–3, except that for GAM-type DPs, the frequency content of each high-frequency component has also been shifted accordingly (by its harmonic frequency). In addition, based on the definition of (2–36), it can be derived that the magnitude of 1st order GAM-type DP (i.e., k =1 for fundamental frequency) is half of the magnitude of the SFA-type DP, as also noted in Figure 2–2 and Figure 2–3. This magnitude definition of GAM-type DPs [(2–35) and (2–36)] will be used consistently throughout this thesis. As seen in (2–35) and (2–36), the GAM enables to select a number of DPs of interest (i.e., selection of k), to construct an adequate approximation of the original time-domain waveform. Specifically, this approximation yields the instantaneous reconstruction [124] as 30 Kktjkksetutu, (2–37) where K={1, 2, 3, …} denotes the set of selected orders of GAM-type DPs. However, it is noted that the inclusion of higher-order DPs also adds to the modelling complexity and computational costs. Moreover, several operation features of GAM-type DPs, which are also used in the following sections, can be obtained as [62] *kkuu , (2–38) iiikkvuuv , (2–39) kskkujkudtddtdu , (2–40) where the superscript * refers to conjugate operation. Figure 2–3. Frequency spectra of the power system signal in different representations during the process of GAM: a) instantaneous signal, and b) GAM-type DP. Figure 2–4. Illustrated “sliding window” of the signal moving along the time axis. 31 Discussions and Computer Studies As shown in Figure 2–2 and Figure 2–3, both the SFA-type DP tU and GAM-type DP tuk are obtained as the analogous low-pass representation of the instantaneous signal tu . From the power system simulation perspective, the use of DPs can be particularly helpful to achieve a significant reduction of the number of necessary time-steps [51]-[52], since many fewer samples are required to accurately represent such low-pass DPs than the original instantaneous signals according to the Nyquist sampling theorem (i.e., maxmax 21 ft , [120]). Several aspects of the aforementioned representations of power system signals are summarized in Table 1–1. Therein, it is noted that as hybrid time-phasor signals, the DPs are complex-value variables, which require real and imaginary parts (i.e., increased system order) in the state-space formulation [as (2–1)–(2–2) or (2–3)–(2–4)] for using regular ODE solvers. Moreover, with the approximated expression (2–37) of GAM-type DPs (i.e., selection of K), one should anticipate a trade-off between the simulation accuracy and modelling complexity. Typically, it is advised that the fundamental component and/or the first few dominant harmonics are considered for common system-level modelling and analysis with GAM-type DPs [62]−[83]. It is worth noting that in the respective literature, the concept of “dynamic phasor” may refer to different methodologies (e.g., for dynamic estimation in phasor measurement unit (PMU) applications [128]), or be named using other terminologies (e.g., “dynamic harmonics” [124], [126]-[127], “time-varying harmonics” [108], etc.). Throughout this thesis, the terminology “dynamic phasor” refers to the SFA- and GAM-type DP modelling approaches. 32 Table 2–1. Properties of different representations of power system signals. Representation Notation Signal Types Assumption Accuracy Instantaneous Signal tu Time-Domain Signal n/a Accurate Conventional Phasor sU Fundamental-Frequency Phasor Neglecting Stator-Network Transients Low SFA-Type DP tU Hybrid Time-Phasor Signals Band-pass Signal Accurate GAM-Type DP tuk Hybrid Time-Phasor Signals Quasi-Periodical (Wideband) Signal Approximated Modelling Basic RLC Components To demonstrate the application of DP-type modelling approaches to power circuits, the examples of basic RLC components are presented here. In general, the procedure of developing power system component models in the DP domain can be summarized as follows: 1) Obtain the dynamic equations in the time-domain using instantaneous signal representation [in state-space formulation as (2–1)–(2–2) or (2–3)–(2–4)] ; 2) Follow the definition of the corresponding type of DPs to transform the time-domain equations. Specifically, for the SFA-type DPs, one should first rewrite the instantaneous signals (voltages, currents, etc.) in the form of (2–27), and then follow the procedure of (2–33) (i.e., Hilbert transform, analytic signal formulation, and frequency shift) [51]-[52]. For the GAM-type DPs, one should first select the set of desired frequencies (i.e., K) based on the properties of non-linear components in the system, and then derive the corresponding DP model based on (2–36) [63]; 3) Manipulate and simplify the DP models derived in 2) [e.g., using the operation features (2–38)-(2–40)] into a compact state-space formulation; 4) Decompose the complex-valued DP dynamic equations into real and imaginary parts to facilitate the numerical integration using regular ODE solvers [27]–[28]. 33 For better illustration of the above procedure, the equations in different signal representations governing the basic RLC components are derived and summarized in Table 2–2. As shown in Table 2–2, both the equations in SFA- and GAM-type DPs closely resemble those formed using the conventional phasors [19], with the difference being the time-derivative terms. In steady state, the SFA-type DP (which is two times the magnitude of the 1st order GAM-type DP) should yield identical solution to those of conventional phasor, thus allowing for potential interfacing between DP models and TS programs [63]. Moreover, for EMT simulations where the sinusoidal signals are simulated, the magnitudes of such DPs essentially represent the envelopes of the time-domain waveforms [51]. Therefore, the DP models may also find potential application in the EMT simulations [75]-[76], by providing straightforward information about the dynamics of the power systems. Table 2–2. Dynamic equations governing the basic RLC components using different signal representations. Representation Resistor R Inductor L Capacitor C Instantaneous Signal tRitv RR tidtdLtv LL tvdtdCti CC Conventional Phasor sRsR RIV sLssL LIjV sCssC CVjI SFA-Type DP tRItV RR tIdtdLtLIjωtV LLsL tVdtdCtCVjωtI CCsC GAM-Type DP tiRtv kRkR tidtdLtiLjktvkLkLskL tvdtdCtvCjktikCkCskC Modelling System Transients at Different Frequencies To demonstrate the effectiveness and limitations of DP-type modelling approaches, here a single-phase series RLC circuit is considered as shown in Figure 2–5 (a). Depending on the system parameters, this circuit is able to exhibit resonant transients from sub- to super-synchronous frequencies. 34 Figure 2–5. The example single-phase series RLC circuit in different signal representations: a) instantaneous signal; b) conventional phasor; c) SFA-type DP; and d) GAM-type DP. 2.3.2.1 Analytical Modelling and Eigenvalue Analysis This single-phase series RLC circuit can be represented by the state-space model in time-domain instantaneous signal representation as tvLtvtiCLLRtvtidtdsCLCL0/1011. (2–41) Neglecting dynamics, the conventional algebraic phasor representation for this circuit is sSsCsLssVLVIjωCLjωLR0/111120 . (2–42) Following the equations derived in Table 2–2, the SFA-type DP model is 35 tVLtVtIjωCLjωLRtVtIdtdSCLssCL0/111, (2–43) and GAM-type DP model is tvLtvtijkωCLjkωLRtvtidtdkSkCkLsskCkL0/111. (2–44) In summary, the RLC circuits in different signal representations as shown in Figure 2–5 can be represented by (2–41)-(2–44). It is noted that the conventional phasor yields pure algebraic relationships (2–42) by neglecting system dynamics. As also noted in Figure 2–5, compared with the time-domain equation (2–41), the complex-valued terms in the diagonal elements of the DP state matrices [i.e., sj in (2–43) and sjk in (2–44)] can be interpreted as the additional voltage source (in series) and current source (in parallel) for inductive and capacitive elements, respectively [e.g., tLIj Ls and tiLjk kLs as shown in Figure 2–5 (c) and (d), respectively]. This circuit representation may find particular use for power system modelling in EMTP-type solutions where the system components are discretized and formulated at branch level [75]-[76]. Moreover, to operate with complex variables using regular ODE solvers [27]–[28], the SFA-type DP model (2–43) and GAM-type DP model (2–44) have to be further decomposed into real and imaginary components, which yield the decomposed SFA- and GAM-type DP models as 36 tVtVLLtVtVtItIωCωCLLRωLωLRtVtVtItIdtdiSrSiCrCiLrLssssiCrCiLrL,,,,,,,,,,0000/100/10100011001, (2–45) ikSrkSikCrkCikLrkLssssikCrkCikLrkLvvLLtvtvtitikωCkωCLLRkωLkωLRtvtvtitidtd,,,,,,,,,,0000/100/10100011001, (2–46) where the additional subscripts “r” and “i” denote the real and imaginary parts, respectively. To investigate the system transients/dynamics at different frequencies, the eigenvalue studies are usually performed in the power system community [19]. The eigenvalues (i.e., system modes) of the time-domain model (2–41) are given by 12 nnTD , (2–47) where n and denote the resonant frequency and damping factor of the series RLC circuit, respectively, as defined by LCRLCn2 ,1 . (2–48) Similarly, the eigenvalues of the decomposed SFA-type and GAM-type DP models (2–45) and (2–46) (which use real and imaginary components of DPs as state variables) are respectively obtained as 37 22_ 1 snnDPSFA jω , (2–49) 22_ 1 snnDPGAM jkω . (2–50) To better observe the resonant transients as predicted by the subject models, the underdamped series RLC circuit can be considered. Specifically, with a small value of damping factor that j12 , the eigenvalues (2–47) and (2–49)-(2–50) are approximated as nnTD j , (2–51) snnsnnDPSFAjj _ , (2–52) snnsnnDPGAMkjkj _ . (2–53) The subject system models and their corresponding properties for a series RLC circuit based on different signal representations are summarized in Table 2–3. Therein, it is noted that the conventional phasor results in an algebraic model (thus no state variables or system eigenvalues), and the SFA- and GAM-type DPs yield dynamic models with increased system orders due to the decomposed real and imaginary representation of complex state variables. Moreover, based on the eigenvalues of the subject models (see last column in Table 2–3), several conclusions can be drawn: 1) Compared with the eigenvalues of the original time-domain model (2–51), the resonance/oscillation modes of the DP models are obtained at complemented frequencies where the original resonant frequency is shifted by the fundamental or harmonic-order frequencies, as 38 shown in (2–52)-(2–53). This is due to the additional elements in the state matrices [i.e., additional off-diagonal elements s in (2–45) and sk in (2–46)]. 2) During simulations, these complemented modes (frequencies) may become excited with the additional ringing in the dynamics of the states (i.e., the real and imaginary components of DPs) after a system input/configuration change. For the considered underdamped series RLC circuit, the additional ringing is introduced at sn and sn k for the decomposed SFA- and GAM-type DP models (2–45) and (2–46), respectively. Moreover, since the selection of integration step sizes in variable-step ODE solvers depends on the signals [27]–[28], these ringing dynamics of DPs may force the ODE solvers to use small-step sizes to track system transients. 3) Despite the additional ringing in DPs, it is noted in (2–51)-(2–53) that the real parts of eigenvalues are identical for different signal representations, i.e., DP modelling approaches do not affect the time constants of the system. Therefore, as the system transients/dynamics die out, the variable-step ODE solver is able to use large time-steps for simulating (quasi-) steady-state conditions. Table 2–3. Properties of system models for the series RLC circuit using different signal representations. Representation System Model State Variables Eigenvalues Instantaneous Signal (2–41) TCL tvti , nn j Conventional Phasor (2–42) N/A N/A SFA-Type DP (2–45) TiCrCiLrL tVtVtItI ,,,, , , , snnsnnjj GAM-Type DP (2–46) TikCrkCikLrkLtvtvtiti,,,, ,... , , snnsnnkjkj 39 2.3.2.2 Computer Studies To facilitate the previous discussions, computer studies are conducted on the series RLC circuit shown in Figure 2–5 with system parameters summarized in Table 2–4. To demonstrate the capability of predicting system transients at different frequencies, the system parameters are tuned [based on (2–48)] to excite resonant transients at 10 Hz (sub-synchronous frequency) and 300 Hz (super-synchronous frequency). The system is assumed to have zero initial conditions. At 1.0ts, the voltage source is turned on (with magnitude V 6.11V s and frequency Hz 60sf ) to excite the dynamic modes of the series RLC circuit. The system response with resonant frequency nf at 10 Hz and at 300 Hz are shown in Figure 2–6 and Figure 2–7, respectively, where the solutions predicted by the time-domain (TD) model (2–41), the conventional steady-state Phasor model (2–42), the SFA-type DP model (2–45), and the GAM-type DP model (2–46) (with K={1}) are shown. For the DP models, the trajectories of the DP magnitudes [calculated using the real and imaginary parts of DPs, and specifically, line d-1) is two times the magnitude of 1st order GAM-type DPs, e.g., 12 Li ] and the reconstructed time-domain values [using (2–34) and (2–37) for SFA- and GAM-type DPs, respectively] are included to compare with the other solutions. As seen in Figure 2–6 and Figure 2–7, for the trajectories of time-domain values, both the SFA and GAM models produce results that accurately match the TD model [see lines a), c-2), and d-2)]. This validates the accuracy of DP models (i.e., full-order modelling) in capturing the resonance transients of the RLC circuit at frequencies ranging from sub-synchronous to high-order harmonic frequencies. 40 It is also observed in Figure 2–6 and Figure 2–7 that the trajectories of DP magnitudes [see lines c-1) and d-1)] represent the dynamic envelopes of the time-domain waveforms. Specifically, it is seen in Figure 2–7 (and in Figure 2–6 provided that the simulation time is long enough) that when the resonance transients die out (in quasi-steady-state conditions), the DP models yield identical trajectories [see lines c-1) and d-1)] to the classical Phasor model solution [see lines b)], which simply neglects the system dynamics. Moreover, it is noted in Figure 2–6 and Figure 2–7 that during the resonance transients, some ringing is presented in the trajectories of the magnitudes of DPs [see lines c-1) and d-1)]. As explained in Section 2.3.2.1, the additional elements in the state matrices [i.e., the additional diagonal elements sj in (2–43) and sjk in (2–44), and the additional off-diagonal elements sin (2–45) and sk in (2–46)] can introduce complemented modes in the dynamics of the states (i.e., the real and imaginary components of DPs), which thus result in additional ringing in the DP magnitudes. Such ringing in the DPs has also been noticed in respective literature [124]-[125], where the authors explain that a direct step change in the DP domain can in fact introduce additional dynamics/oscillations into the associated DPs during transients, while precisely reproducing the correct time-domain results. Basically, when mapping a step change of a time-domain signal into the DP domain, two different approaches can be employed [124], i.e., implementing: 1) a set of continuous time-varying changes of DP inputs; or 2) a direct step change in the DPs, which is simpler and more widely used [124], and has been considered throughout this thesis. However, the direct step change in the DP domain can cause inconsistency between responses of the DP systems before and after the input change [125], and provoke oscillations in the DPs that are damped according to the time constants of the system [as also concluded from (2– 41 51)-(2–53), Section 2.3.2.1]. Finally, it is important to mention that, despite the ringing during transients, both the SFA- and GAM-type DP models reproduce accurate time-domain simulation results that are well matched with the TD model, as validated in Figure 2–6 and Figure 2–7 [lines a), c-2), and d-2)]. This validates the DP models in predicting system responses at different frequencies. Table 2–4. Parameters of the series RLC circuit for different resonant frequencies. System Parameters Resonant Frequency nf Damping Factor 1 =R , H 0.1 =L , F 2533 C 10 Hz 0.0796 10 =R , H 0.1 =L , F 815.2 C 300 Hz 0.0265 Figure 2–6. System response of the RLC circuit with resonant frequency Hz 10nf . 42 Figure 2–7. System response of the RLC circuit with resonant frequency Hz 300nf . The origin of the ringing in DPs can also be verified by investigating system frequency-domain characteristics and eigenvalues of the subject models. For the TD model (2–41), the input admittance of the RLC circuit [i.e., the transfer function from input Sv to output Li ] is depicted in Figure 2–8 and Figure 2–10 for the system with resonant frequencies at 10 Hz and 300 Hz, respectively. For DP models, since with K = {1} the GAM model (2–46) and SFA model (2–45) yield identical formulation, only the results from the SFA model are included. It is also noted that due to the decomposed real and imaginary components of DPs, the input admittance is now represented as a matrix of 2×2 transfer functions from input TiSrS VV ] ,[ ,, to outputTiSrS II ] ,[ ,, . The 43 DP model results are depicted in Figure 2–9 and Figure 2–11, each corresponds to the results shown in Figure 2–8 and Figure 2–10, respectively. As seen in Figure 2–8 and Figure 2–9, compared with the TD model with Hz 10nf (sub-synchronous resonance), the SFA model attains resonant frequencies at 50 Hz and 70 Hz [see the double peaks f admittance magnitudes], which matches the complemented frequencies at sn ( sn ff ) as concluded in Table 2–3. Similarly, for the super-synchronous resonances seen in Figure 2–10 and Figure 2–11, the SFA model is shown to shift the resonant frequency from the original Hz 300nf by the fundamental frequency sf , thus resulting in the two complementary resonances at 240 and 360 Hz [i.e., sn ff ]. These complemented resonances of the SFA model are also verified by the eigenvalues as summarized in Table 2–5, where it is seen that the frequencies of system modes (i.e., the imaginary components of the eigenvalues, in rad/s) match well with the previous frequency-domain analysis. Table 2–5. Eigenvalues of the subject models for the series RLC circuit with different resonant frequencies. Model System Equation Resonant FrequencyHzfn 10 Resonant FrequencyHzfn 300 TD (2–41) 63.625 j 28.18845 j SFA (2–45) 36.314562.4395jj 29.1507528.22615jj 44 Figure 2–8. Transfer function from input Sv to output Li as predicted by the TD model for the series RLC circuit with resonant frequency Hz 10nf . Figure 2–9. Transfer function from input TiSrS VV ] ,[ ,, to output TiSrS II ] ,[ ,, as predicted by the SFA model for the series RLC circuit with resonant frequency Hz 10nf . 45 Figure 2–10. Transfer function from input Sv to output Li as predicted by the TD model for the series RLC circuit with resonant frequency Hz 300nf . Figure 2–11. Transfer function from input TiSrS VV ] ,[ ,, to output TiSrS II ] ,[ ,, as predicted by the SFA model for the series RLC circuit with resonant frequency Hz 300nf . 46 Modelling Power Systems under Unbalanced Conditions As demonstrated in Section 2.3.2, the SFA- and GAM-type DP models are shown to accurately predict the system transients at different frequencies. However, it is reminded that the superior numerical efficiency of DP simulations is achieved through the low-pass representation of power system signals, e.g., the dc-like trajectories of envelopes/states in quasi-steady state as shown in Figure 2–7. As also noted in (2–8) and (2–20), the time-domain waveforms of power system signals may attain different frequencies in the abc-phase and qd0 coordinates, especially under unbalanced conditions. To achieve numerically efficient simulation, it is therefore crucial to select the proper coordinates for DP representation of power system signals. For this purpose, in this section the SFA- and GAM-type DPs for representing the symmetrical components of unbalanced power systems are set forth. 2.3.3.1 SFA-Type DPs of Three-Phase Power System Signals For SFA-type DPs, it is recalled in (2–8) that under unbalanced conditions the abc-phase power system signals attain a band-pass bandwidth for all symmetrical components. Taking phase b as an example, it is recalled in (2–8) that the time-domain signal is zsznsnpspb ttttu cosU32cosU32cosU , (2–54) Applying (2–31) then yields the analytic signal as zsnsps tjztjntjpbbbeeetujHtutU UUUˆ3232, (2–55) and the SFA-type DPs are then obtained after shifting the frequency as 47 znpsjzjnjptjbbeeeetUtUUUUˆ3232 . (2–56) Therefore, following the procedure in Section 2.3.1, the SFA-type DPs of abc-phase power system signals under unbalanced conditions are expressed as )()()()32()32()()32()32()(UUUzzznnnpppjjjzjjjnjjjpabceeeeeeeeetU . (2–57) Moreover, it is convenient to define the SFA-type DPs for symmetrical components tpnzU as )()()(UUUznpjzjnjpznppnzeeetUtUtUtU . (2–58) Therefore, (2–57) can be manipulated into ttUtUtUtUtUtUt pnzznpcbaabc MUMU , (2–59) where M is defined in (2–24). It is noted that (2–59) is identical to the fundamental-frequency phasor expression of symmetrical components (2–23), which infers the potential interfacing between SFA-type DP models and TS programs [63]. However, for the qd0 coordinates, it is noted that under unbalanced conditions, the power system signals are no longer of band-pass due to the zero- and negative-sequence components at fundamental and double-fundamental frequencies, respectively [as noted in (2–20)]. In this case, it becomes impractical to use the SFA-type DPs for transformed qd0 signals, since low-pass representation will no longer be achieved. Specifically, applying (2–31) to (2–20) yields 48 )()22()2(0 00U0U0sincosUˆzsnsnstjztjtjnpppqdeeetU , (2–60) where it is noted that the Hilbert transform of a constant value is zero. Therefore, after shifting the frequency (2–33), the SFA-type DPs of qd0-coordinate signals are obtained as )()2()(0 00U0U0sincosUznsnsssjztjtjntjptjppqdeeeeetU . (2–61) 2.3.3.2 GAM-Type DPs of Three-Phase Power System Signals Similarly, following the procedure in Section 2.3.1, the GAM-type DPs of abc-phase signals can be derived by selecting the fundamental frequency components (i.e., K={1}) and then applying (2–36) to (2–8), which yields )()()()32()32()()32()32()(1 2U2U2Uzzznnnpppjjjzjjjnjjjpabceeeeeeeeetu . (2–62) For the qd0 coordinates under unbalanced conditions, the order of GAM-type DPs should be selected as K = {0, 1, 2} to adequately represent the positive-, zero- and negative-sequence components as dc, fundamental, and double-the-fundamental frequencies, respectively. Therefore, the GAM-type DPs of qd0-coordinate signals are obtained by applying (2–36) to (2–20) as 0sincosU00 pppqd t u , (2–63) 49 Table 2–6. Properties of different signal representations for the three-phase power systems under balanced/unbalanced conditions. Representation Coordinates Equation Balanced ( pU ) Unbalanced ( znp U;U;U ) Signal Freq. No. of State Variables Signal Freq. No. of State Variables Instantaneous Signal abc (2–8) s 3 sss ; ; 3 qd0 (2–20) dc 3 ssdc ;2 ; 3 Conventional Phasor abc (2–23) dc N/A dc N/A SFA-Type DP abc (2–57) dc 3 complex (6 total) dc 3 complex (6 total) qd0* (2–61) s 3 complex (6 total) dc ; ; ss 3 complex (6 total) GAM-Type DP abc (2–62) dc 3 complex (6 total) dc 3 complex (6 total) qd0 (2–63)- (2–65) dc 3 real dc 3 real 6 complex (15 total) * While deriving the SFA-Type DPs for qd0-coordinate signals is possible, it is impractical since the low-pass signal representation vanishes and no numerical benefit has been gained. )(1000U21zjzqdetu , (2–64) 0U21 )2()(20nnjjnqd eetu . (2–65) In summary, the properties of different signal representations for the three-phase power systems under balanced/unbalanced conditions are summarized in Table 2–6. Therein, the signal frequencies are given assuming that the system enters (quasi-) steady state (i.e., not during transients), where pU , nU and zU are constant or slow-varying values. As seen in Table 2–6, the low-pass (dc) representation of power system signals (superior numerical efficiency) can be achieved with the SFA-type DPs of abc-phase signals, and the GAM-type DPs of both abc-phase variables and transformed qd0 signals. Moreover, to operate with complex variables using regular 50 ODE solvers [27]–[28], the SFA- and GAM-type DPs have to be further expanded into real and imaginary components, thus increasing the number of system order (as shown in Table 2–6). For example, to adequately represent the unbalanced/faulty power systems, the order of GAM-type DPs should be selected as K = {0, 1, 2}, which requires a total of 15 variables (i.e., Tiqdrqdiqdrqdqd],,,,[,20,20,10,1000uuuuu ) that significantly adds to model complexity and associated computation costs. Therefore, it can be concluded that the abc-phase coordinates are the most numerically efficient for DP modelling of power system signals. 2.3.3.3 Computer Studies To validate the previous discussions, computer studies are conducted on a simplified 230 kV three-phase transmission system as shown in Figure 2–12. Therein, the power grid is modelled by the Thevenin voltage source with an impedance corresponding to a 2000 MVA short circuit level and X/R = 10. The 50-km transmission line is modelled by a three-phase coupled π section RLC circuit (short line model) [20], the detailed parameters of which can be found in Appendix A. The load (75 MW, 20 Mvar per phase) is modelled by the parallel RL loads. This system is assumed to be three-phase symmetrical with parameters summarized in Table 2–7. Figure 2–12. Circuit diagram of the considered 230 kV three-phase transmission system. 51 Table 2–7. Parameters of the considered three-phase transmission system. System Component Parameters Power Grid peak kV 187.8 =Vphase , Hz 60sf , 2.645 =sR , mH 70.2 =sL Transmission line 5.667 =_ stlR , 3.667 =_ mtlR , mH 73 =_ stlL , mH 30.8 =_ mtlL , F0.2615 =_ stlC , F0.035- =_ mtlC System Load 235.111 =loadR , H 2.339 =loadL To emulate an unbalanced situation, the system (initially assumed in the rated steady state) is subjected to a single-phase-ground fault in the phase c voltage source at t = 0.1s, and the simulation is run until t = 5s. For comparison, the subject models have been implemented in MATLAB/Simulink [27]–[28], including the conventional phasor model [a) Phasor], the time-domain model represented in abc-phase [b) TD - abc] and qd0 coordinates [c) TD - qd0], the GAM-type DP model in abc-phase [d) DP- abc, with K ={1}] and qd0 coordinates [e) DP- qd0, with K ={0, 1, 2}]. Here, only the GAM-type DP models have been implemented, since the SFA- and GAM-type DPs (with K={1}) yield the same model formulation for abc-phase signals (with the only difference being the phasor magnitudes). To verify the numerical accuracy of the subject models, the system is first solved with small time step sizes using the ODE23tb solver [27] and the following settings: relative and absolute error tolerances of 10−3, and maximum and minimum step sizes of 10µs and 0.1 µs, respectively. For the purpose of discussion, only the responses of phase c current csi is depicted in Figure 2–13. As can be seen in Figure 2–13, all of the TD and DP models can accurately predict the time-domain trajectories of current csi [see lines b) to e)], except for the reduced-order Phasor model [see line a)]. In particular, it is observed that the current csi , after the fault, has oscillatory 52 behaviour at three frequencies: 1) the 60 Hz system fundamental frequency, 2) the high-frequency RLC resonant transients that decays rapidly, and 3) the dc offset that decays exponentially in a few seconds. It is noted that such dc offset essentially represents the stator-network dynamics of power systems (Section 3.7.2, [19]), which is generally neglected in classical phasor models and TS programs. Therefore, Figure 2–13 confirms the validity and numerical consistency of the DP models in capturing the details of fast transients similarly to TD models, provided that simulation step sizes are sufficiently small. Next, to investigate the numerical efficiency of the subject models, it is desirable to study the behaviour of the system state variables during and after the fault transient, as shown in Figure 2–14 and Figure 2–15, respectively. As can be seen in Figure 2–14 and Figure 2–15, in the initial steady state, the state variables of the subject models (except for the TD-abc model) are dc signals, which is consistent with Table 2–6. However, as the single-phase-ground fault occurs at t = 0.1s (i.e., unbalanced condition), the trajectories of state variables for the subject models differ greatly. As seen in Figure 2–14, the Phasor algebraic model yields state trajectories as magnitude step change, and the TD-abc model yields time-domain waveforms that are same as those shown in Figure 2–13(b). For the TD-qd0 model, in addition to the rapidly-damped transients (due to RLC resonances), it is seen that the unbalanced fault introduces negative- and zero-sequence components into the system and the waveforms oscillate at s2 and s , respectively [even in the post-transient condition shown in Figure 2–15 (c)]. Nevertheless, this loss of low-pass representation does not occur with the DP models: for both the DP- abc and DP- qd0 models, the system state variables become the dc values (low-pass representation) as the system enters close to the quasi-steady state, as shown in Figure 2–15 and Table 2–6. 53 The numerical advantages of DP models are further demonstrated by solving this case study using the same ODE23tb solver and error tolerances, but increasing the maximum step size to 16.67 ms (i.e., one cycle at 60 Hz). This allows the solver to adaptively select the step sizes in a wide range (even at step sizes comparable with TS programs). Figure 2–16 depicts the step size Δt as chosen by the subject models, where it is seen that in the initial steady state the maximum step sizes are used (except for the TD-abc model). When the fault occurs at t = 0.1s, the system state variables vary accordingly, which results in a drop of Δt to very small values in order to trace the transient response with good accuracy. Moreover, after the transient, it is noted in Figure 2–16 that the step sizes of DP models increase rapidly due to the low-pass representation of system states. In contrast, Δt remains fairly small with the TD model, which validates the superior numerical efficiency of the DP models. Figure 2–13. System response of csi as predicted by the subject models for the single-phase-ground fault in phase-c voltage source. 54 Figure 2–14. Trajectories of system states during the fault transient for the subject models: (a) Phasor; (b) TD-abc; (c) TD-qd0; (d) DP-abc; and (e) DP-qd0. Figure 2–15. Trajectories of system states after the fault transient for the subject models: (a) Phasor; (b) TD-abc; (c) TD-qd0; (d) DP-abc; and (e) DP-qd0. 55 Figure 2–16. Step size Δt as chosen by the subject models. To conclude this chapter, several remarks are made: 1) Both the SFA- and GAM-type DP models are shown to yield time-domain results that match well to the reference solution produced by the TD model. This is particularly favorable for investigating the transient phenomena associated with integrated ac-dc power systems (at sub/super-synchronous frequencies), as the DP solutions can accurately capture both the fast and slow transients/dynamics (provided that the time step is sufficient). 2) As system transients/dynamics die out (in quasi steady states), the DPs (and their real and imaginary components) become constant or slow-changing variables, which allow the variable-step ODE solver to adaptively use flexible (large) time-step sizes sufficient for accurate capturing 56 of the active modes of the systems. It is also noted that the additional ringing with DP models does not affect the system’s natural time constants. 3) For general power system dynamics around the fundamental 50/60 Hz frequency (e.g., transient stability studies), the abc-coordinate SFA- and GAM-type DPs (K={1}) yield equivalent solutions (except for the factor of 1/2 scaling in magnitude for GAM-type DPs). For qd0-coordinates, only the GAM-type DPs can obtain the low-pass representation of signals (i.e., use of large time-step sizes). It is also noted that when simulating unbalanced system conditions, the GAM-type DPs in qd0-coordinates do require higher system orders (K = {0, 1, 2}), which come at additional associated computation costs. 4) For modelling the major/dominant harmonics (not transients) in the system, the SFA-type DP approach gains no additional numerical benefits, since it only shifts the signal by the fundamental frequency. In this case, the GAM-type DPs are preferred, which however increases the modelling complexity by requiring more terms. 57 Numerically Efficient Modelling of Electric Machines Based on the Shifted Frequency Analysis This chapter presents numerically efficient modelling of electric rotating machines based on the SFA-type DPs. Owing to the advantageous numerical properties discussed in Section 1.2.2.1, the VBR formulation of electric machines has been used for the state-of-the-art SFA-based models [53]-[54], and is considered in this chapter. We begin by proposing a new VBR synchronous machine model based on SFA. By re-formulating the machine equations, a constant-parameter DP stator-network interface is achieved with this model, which is simple to implement and numerically-free from the rotor-position-dependent inductances owned by the prior state-of-the-art model [53]. Next, this DP modelling approach is extended to obtain a SFA-based VBR model of induction machines in state-space formulation. Case studies in conjunction with error and efficiency analysis are presented to highlight the superior combination of numerical accuracy and efficiency of the proposed models. 58 Constant-Parameter VBR Modelling of Synchronous Machines Based on SFA-Type DPs Without loss of generality, a three-phase synchronous machine model with a field winding fd and one damper winding kd in the d-axis, and two damper windings kq1 and kq2 in the q-axis, is considered for this section. Motor sign convention is used, and all parameters are referred to the stator side. It is assumed that the machine is magnetically linear, and that the q-axis leads the d-axis by 90° [84]. Time-domain VBR Synchronous Machine Models Before deriving the SFA models, the time-domain equations of VBR synchronous machine models are first reviewed [88], [92]. The mechanical subsystem is represented by a single rigid body with its dynamics represented by the following equations rrdtd , (3–1) mer TTJPdtd2 . (3–2) Here, P is the number of poles and J is the combined rotor-load inertia. The rotor position and angular electrical speed are respectively denoted by r and r . Variables mT and eT denote the mechanical and electromagnetic torque, respectively. 3.1.1.1 Variable-Parameter Time-Domain VBR Model (VP-TD) The stator voltage equation of the variable-parameter time-domain (VP-TD) VBR model can be expressed as [88] 59 '''' ])([ abcsabcsrabcsabcssabcsdtdeiLiRv . (3–3) Here, the stator resistance matrix sR is ssss rrrdiag , ,R , (3–4) and the subtransient inductance matrix rabcs "L is )342()2()322()2()342()322()322()322()2()(''rSrMrMrMrSrMrMrMrSrabcsLLLLLLLLLL, (3–5) where the entries are cosbalss LLLL , cos2baM LLL , (3–6) 3''''mdmqaLLL , 3''''mqmdbLLL . (3–7) The inductances ''mdL and ''mqL are 1" 111lkdlfdmdmdLLLL , 121" 111lkqlkqmqmqLLLL . (3–8) The subtransient voltages ''abcse in (3–3) are Tdqrsabcs ee ]0[][''''1'' Ke , (3–9) where rsK is Park’s transformation matrix referred to the rotor reference frame [84] and 60 2222''2111'''''')()(lkqkqmqkqmqlkqkqmqkqmqdrqLrLLrLe , (3–10) . )()(2''''2''''''lfdfdmdfdmdfdlfdmdlkdkdmdkdmdqrdLrLvLLLrLe (3–11) The subtransient flux linkages are )(2211''''lkqkqlkqkqmqqLLL , )(''''lkdkdlfdfdmddLLL , (3–12) and the magnetizing flux linkages are '''' qqsmqmq iL , ''''ddsmdmd iL . (3–13) The rotor dynamics are represented by the following state equations 2,1);( kqkqjLrdtdmqjljjj , (3–14) kdfdjvLrdtdjmdjljjj ,;)( . (3–15) Here, j denotes rotor flux linkages. Finally, the electromagnetic torque is calculated as )(43dsmqqsmde iiPT . (3–16) 3.1.1.2 Constant-Parameter Time-Domain VBR Model (CP-TD) For general synchronous machines, ''mdL and ''mqL are not equal – a property known as dynamic saliency [89]. In such cases, the subtransient inductance matrix rabcs "L in (3–3) 61 becomes rotor-position-dependent [see (3–5)–(3–8)]. Consequently, numerically expensive recalculation of the state matrices is required at every time step. In order to obtain a constant-parameter stator interface and thus a much more efficient solution, one approach is to add an artificial damper winding, which requires specific tuning of its parameters [89]-[91]. Another approach proposed in [92] is to first derive an implicit constant-parameter VBR formulation, and then use numerical approximations to break the algebraic loop. This approach is simple to implement and has good numerical accuracy [91]-[92], and is therefore considered for this section. The stator voltage equation of the constant-parameter time-domain (CP-TD) VBR model is given as ''abcsabcsabcssabcsdtdeiLiRv , (3–17) where the inductance matrix becomes constant as ''''''''''''''''''''SMMMSMMMSLLLLLLLLLL , (3–18) where 32 ''" mdlsSLLL , 3''" mdMLL . (3–19) The new subtransient voltages in qd coordinates in (3–17) and (3–9) are modified as qssqsqdqlkqkqmqkqmqlkqkqmqkqmqdsdqdrqdqirvLLLLrLLrLiLLLLe ''''''2222''2111'''''''''''''')()(])([, (3–20) 62 fdlfdmdlfdfdmdfdmdlkdkdmdkdmdqsdqrqrdvLLLrLLrLiLLe''2''2'''''''''')()()(, (3–21) where the subtransient inductances "dL and "qL are '''' mdlsd LLL , ''''mqlsq LLL . (3–22) It is noted that in general cases where the machine is connected to an inductive network, the voltages abcsv are unknown (i.e., neither states nor inputs). Therefore, the presence of qsv in (3–20) results in an algebraic loop (i.e., an implicit model). To achieve an explicit formulation, qsv is herein approximated by a first-order low-pass filter as [92] qsqs vpspv00~ , (3–23) where qsv~ indicates the approximated value. It is worth noting that, depending on the required approximation accuracy, the pole(s) of such a low-pass filter can be properly selected following the procedure in [91] by investigating its impact on the machine’s operational impedances. In summary, the CP-TD VBR model is comprised of the stator interfacing equation (3–17), the subtransient voltages (3–9), (3–20), and (3–21), the voltage approximation (3–23), the rotor state model (3–14), (3–15), and the electromagnetic torque (3–16). VBR Synchronous Machine Models Based on SFA The SFA VBR models can now be derived by transforming the electrical signals of the time-domain VBR models in Section 3.1.1 into SFA-type DPs following the approach set forth in Section 2.2.1. 63 3.1.2.1 Variable-Parameter Stator Interface Based on SFA The variable-parameter SFA-based VBR synchronous machine model was proposed in [53]. The corresponding stator interfacing equation using SFA-type DPs is ''*)22(''12''1*)22(''12''1''0)2()2(][abcsabcstjabcsjabcstjsrabcsjsrabcssabcsabcssabcsdtdedtdeejejjdtdsrrsrrEILILILILIILIRV , (3–24) where *abcsI refers to the complex conjugate of abcsI , and the stator inductances are comprised of two parts: alsaaaalsaaaalsLLLLLLLLLLLL222222''0L and 222''11112aaaaaaLbL , 3/2jea . (3–25) To obtain a compact formulation and facilitate complex variable operations using regular SV-based simulation programs, (3–24) is decomposed into real and imaginary parts as '' ,'' , , , , ,'' , ,) ,() ,(iabcsrabcsiabcsrabcsrSFAVPiabcsrabcsrSFAVPiabcsrabcstdtdtEEIIRIILVV , (3–26) where the additional subscripts “r” and “i” denote the real and imaginary parts, respectively. The equivalent inductance and resistance matrices trSFAVP ," L and trSFAVP ,R are rriiiirrrSFAVP t ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1'' ) ,(BABABABAL , (3–27) rriiiirrrSFAVP t ,2 ,2 ,2 ,2 ,2 ,2 ,2 ,2) ,(BABABABAR , (3–28) where the entries are 64 .)2()2()22(''122''12)22(''112''11tjsrjsrsstjjsrrsrrejejjeeLBLLRALBLLA''0''0 (3–29) 3.1.2.2 Constant-Parameter Stator Interface Based on SFA As can be seen from (3–27)–(3–29), in order to evaluate trSFAVP ," L and trSFAVP ,R in (3–26), the rotor-position-dependent complex matrices A1, A2, B1, and B2 must be recalculated, decoupled, and manipulated at every time step, which significantly increases the computational burden. To achieve a constant-parameter stator interface in DP form, the SFA approach is applied to the CP-TD VBR model. First, (3–17) is rewritten in the form of (2–27) as ].sincos[]cossin[]sincos[]sincos[sincostωtωtωtωωtωdtdtωdtdtωtωtωtωs''abcsQs''abcsIsabcsQsabcsIssabcsQsabcsIsabcsQsabcsIssabcsQsabcsIEEiiLiiLiiRvv (3–30) Applying the Hilbert transform to (3-30) then yields ].cossin[]sincos[]cossin[]cossin[cossintωtωtωtωωtωdtdtωdtdtωtωtωtωs''abcsQs''abcsIsabcsQsabcsIssabcsQsabcsIsabcsQsabcsIssabcsQsabcsIEEiiLiiLiiRvv (3–31) Combining (3-30) and (3-31) into the form of (2–33), the stator equation in the SFA-type DP form becomes ''abcsabcsssabcsabcs jdtdEILRILV . (3–32) Finally, (3–32) is decoupled into real and imaginary components, yielding 65 '' ,'' , , , , ,'' , ,iabcsrabcsiabcsrabcsSFACPiabcsrabcsSFACPiabcsrabcsdtdEEIIRIILVV, (3–33) where the constant-parameter equivalent subtransient inductance and resistance matrices " SFACPL and SFACPR are given as LLL00''SFACP, (3–34) ssssSFACPRLLRR. (3–35) 3.1.2.3 Rotor Subsystem and Its Interface In time-domain VBR models, the rotor dynamics are expressed in qd coordinates, and included in the subtransient voltage source ''abcse through the inverse Park’s transformation. For SFA VBR models, it is convenient to keep the rotor subsystem in time-domain qd coordinates as in (3–14)–(3–16), since these rotor variables also change slowly during electromechanical transients. This requires transforming the time-domain q- and d-axis subtransient voltages in the rotor reference frame into abc-coordinate DPs in order to construct the equivalent DP subtransient voltage source "abcsE . In order to obtain such a transformation, the time-domain subtransient voltage "ase is recalled from (3–9) as rdrqas eee sincos'''''' , (3–36) 66 which can be interpreted as the combination of two orthogonal low-pass signals ''qe and ''demodulated by rcos and rsin , respectively. In particular, it is observed that this formulation is identical to (2–27). Applying the Hilbert transform to (3-36) then yields rdrqas eeeH cossin]['''''' . (3–37) Combining (3-36) and (3-37) into the form of (2–33), the subtransient voltage "asE in DP form can be expressed as ''2''''dtjqtjas eeeeEsrsr . (3–38) Using the same approach for "bse and "cse , the subtransient voltage source "abcsE in SFA-type DP form is thus given as TdqabcUqdruabcs ee ]['''',,''KE , (3–39) where )322()32()322()32()2()(,,tjtjtjtjtjtjabcUqdrusrsrsrsrsrsreeeeeeK . (3–40) To compute the rotor equations and the electromagnetic torque, the DP variables from the stator interfacing equation must be transformed into qd0-coordinates. The transformation from abc-phase DPs to instantaneous signals in qd0-coordinates, i.e., from abcsF to sqd 0f ( ivf , ), can be readily obtained by combining (2–34) and Park’s transformation as . ][0tjabcsrssqdseFKf (3–41) 67 Figure 3–1. Implementation of the proposed CP-SFA VBR model. In summary, the proposed constant-parameter model based on SFA (herein referred to as CP-SFA) consists of the stator interfacing equation (3–33), the subtransient voltages (3–20)-(3–21), the voltage approximation (3–23), the transformations between rotor and stator variables (3–39) and (3–41), and the rotor dynamics and electromagnetic torque (3–14)–(3–16). A block diagram depicting the implementation of this CP-SFA VBR model in a SV-based simulation environment is shown in Figure 3–1. Therein, the machine stator and network dynamics [(3–33), etc.] are represented by DPs using a constant-parameter state-space formulation, while the rotor dynamics [(3–14)–(3–16)] are expressed in time-domain instantaneous signals. The stator and rotor variables are interfaced using (3–39) and (3–41), and are solved simultaneously within the overall SV-based program. Moreover, to achieve an explicit formulation, it is noted that qsv , obtained by applying (3–41) to abcsV , is fed into the voltage approximation (3–23). This approximation can be fine tuned (i.e., pole(s) of the low-pass filter) based on the required accuracy 68 [91]. The implementation of the variable-parameter SFA-based VBR model (VP-SFA) is similar. However, the stator equation is replaced by (3–26) [i.e., with variable-parameter state matrices], and the subtransient voltages are replaced by (3–10) and (3–11). Additionally, the VP-SFA model does not require the approximation of qsv . Computer Studies To analyze the aforementioned numerical properties of the models, a single-machine infinite-bus (SMIB) system is considered. This allows focusing the investigation on the numerical accuracy and efficiency of each individual machine model. Synchronous machine parameters are obtained from [84] (pp 220, Steam Turbine Generator) and are shown in Appendix B. To take advantage of existing variable-step solvers, all four state-space VBR models (i.e., the time-domain models VP-TD and CP-TD, and the SFA-based models VP-SFA and CP-SFA) have been implemented in MATLAB/Simulink [27]-[28]. Despite the machine being connected to an infinite bus, the voltage filter (3–23) is used in the CP models as to emulate a general scenario. The reference solutions are obtained by solving the conventional qd model with the explicit Runge-Kutta 4th order solver (ode4 in MATLAB) using a very small integration step size of 1 μs. For evaluating the numerical accuracy of all subject models, the 2-norm cumulative relative error [121] of the predicted solution trajectory is used 100%~)(22 ffff . (3–42) Here, f denotes the solution trajectory obtained from the subject model, and f~ denotes the reference solution, obtained over the duration of the study. In this thesis, the 2-norm [129] is 69 selected to evaluate the numerical accuracy of the subject models (i.e., yielding the cumulative error). Other norms, e.g., the one norm 1f or the infinity norm f , can also be used. However, since all norms are equivalent on finite-dimensional vector spaces [121], [129], they will yield equivalent assessment of numerical accuracy of the subject models. For consistency, all case studies are conducted on a personal computer (PC) with a 3.40GHz Intel i7-2600 CPU. 3.1.3.1 Transition from Fast Transient to Steady State To validate the numerical accuracy of the proposed CP-SFA VBR model while also illustrating its computational efficiency, a case study comprising the transition between fast transient and steady-state conditions is considered first. The system is assumed to initially operate in a steady-state condition with mechanical torque mT = 0.85 pu and field excitation voltage xfdE= 2.48 pu. At t = 1.0 s, a symmetric three-phase fault is applied to the stator terminals. The fault is subsequently cleared at t = 1.3 s, and the simulation is run until t = 8.0 s with the purpose of covering both fast transient and steady-state conditions. Since the selection of step sizes by explicit solvers (such as ode45) is limited by stability requirements [27], here the implicit ode23tb solver is chosen. The relative and absolute error tolerances are set to 10−4, and the maximum and minimum step sizes are set to 0.2 s (i.e., twelve cycles at 60 Hz) and 0.1 µs, respectively. The resulting transients of the stator current asi , the electromagnetic torque eT , and the field current fdi as predicted by the four subject models are depicted in Figure 3–2. Magnified fragments of Figure 3–2 showing asi during fast transient and steady-state conditions are also shown in Figure 3–3 and Figure 3–4, respectively. The SFA models yield stator currents in SFA- 70 type DP form asI (i.e., the envelope of the time-domain waveform). These are shown using dashed-dotted lines [see lines e) and g) in Figure 3–3 and Figure 3–4], and are labeled as VP-SFA (DP) and CP-SFA (DP). To adequately compare the SFA models with the reference solution and the time-domain models, the resulting trajectory of asI has also been converted to time-domain instantaneous values asi using (2–34) at the points where the DP solution exists. These converted solutions are shown using circles and asterisks in Figure 3–3 and Figure 3–4 [see lines d) and f)], and labeled as VP-SFA and CP-SFA. The current asi predicted by the time-domain VP-TD and CP-TD models are shown in dotted and dashed lines [i.e., lines b) and c)]. As can be observed in Figure 3–3, during the fast transient, the models of both types, TD and SFA, all use relatively small step sizes and produce results very close to the reference solution. This observation confirms the validity and numerical consistency of these models. Moreover, it shows that the SFA models are able to capture the details of fast transients similarly to the TD models, provided that the simulation step sizes are sufficiently small. After the fault has been cleared and the fast transient has vanished, i.e., at approximately t = 3.5 s [see Figure 3–2], the system reaches a steady state. From this time on, the SFA and TD models choose very different step sizes as shown in Figure 3–4. On the one hand, despite being in steady state, the TD models still require relatively small step sizes to accurately capture the 60 Hz sinusoidal waveform of the instantaneous signal asi within the specified error tolerance. On the other hand, since the envelope of asi is constant in steady state, the SFA models compute the DP asI with the same specified error tolerance but using significantly larger step sizes. 71 To further investigate the behaviour of the SFA and TD models during the transient simulation, the step size Δt for each model is depicted in Figure 3–5 as a function of time. As can be seen in Figure 3–5 (b), the step sizes of the TD models remain relatively small even when the system enters the steady state. At the same time, as shown in Figure 3–5, when the fast transient and oscillations have essentially damped out (i.e., starting from around t = 3.5 s), the step sizes of the SFA models begin to increase rapidly. In this scenario, near t = 8.0 s, the SFA models use step sizes of approximately 0.2 s, whereas Δt does not grow larger than 0.5ms with the TD models. The number of simulation time steps and the CPU time for the considered fault study are summarized in Table 3–1 for each subject model. As shown in Table 3–1, the SFA models require approximately six times fewer time steps than the TD models (4903 and 4871 steps for the SFA models vs. 29707 and 30412 steps for the TD models). Consequently, the CPU time of the VP-SFA model is about 4.05 times smaller than that of the VP-TD model (547.49 vs. 2217.93 ms). Moreover, the models with constant parameters significantly outperform their variable-parameter counterparts in terms of CPU time. Specifically, due to its constant-parameter stator interface, the proposed CP-SFA model is about twice as fast as the VP-SFA model, and it runs about 4.72 times faster than the CP-TD model (272.68 vs. 1287.99 ms). To compare the numerical accuracy of the considered models, the 2-norm relative error of the stator current asi is summarized in the last column of Table 3–1. Only the accuracy of asi is considered here, since the errors of the other variables are comparatively smaller. For the SFA models, the values of stator current asi are calculated from the corresponding DP asI according to (2–34). It is shown that despite requiring six times fewer time steps than the TD models, the SFA models have acceptable errors (1.13% and 1.39% for the VP-SFA and CP-SFA models, respectively). Moreover, the errors of the CP-TD and CP-SFA models are only slightly increased 72 compared to their variable-parameter counterparts, which is due to the approximation of qsv using the low-pass filter (3–23) with p0 = 1000. Additional information about the accuracy of the filters can be found in [91], [92]. Figure 3–2. Fast transient and steady-state responses to a three-phase fault as predicted by the subject models. 73 Figure 3–3. Detailed view of current asi during the three-phase fault: (a) magnified view during fast transient; (b) further magnified view of the portion in part (a). Figure 3–4. Detailed view of current asi after removal of the three-phase fault when the system enters steady state. 74 Figure 3–5. Step size Δt as chosen by the subject models: (a) during the entire fault study; and (b) magnified view from part (a). Table 3–1. Simulation Efficiency and Accuracy for the Fast-transient-to-steady-state Study. Model No. of Time Steps CPU Time, ms 2-Norm Error of asi , % VP-TD 29,707 2,217.93 0.0601 CP-TD 30,412 1,287.99 1.0549 VP-SFA 4,903 547.49 1.1383 CP-SFA 4,871 272.68 1.3914 3.1.3.2 Slow Transient To further demonstrate the advantageous numerical properties of the proposed CP-SFA model, a slower electromechanical transient on the same system is considered. The system is 75 assumed to initially operate in the same steady-state condition. At t = 0.1 s, the load torque is changed to mT = 1.35 pu and the simulation is run until t = 4.0 s. The same solver, error tolerances, and step size settings as the study in Section 3.1.3.1 are used here. The predicted trajectories of asi , eT , and fdi by the four subject models are shown in Figure 3–6. For a better comparison, a magnified view of asi is also shown in Figure 3–7. As can be observed in Figure 3–6 and Figure 3–7, all models accurately predict the solution of asi , and produce visibly indistinguishable trajectories of eT and fdi . Since the inertia constant of the machine is relatively large (H = 5.6 s), the response of the system to this electromechanical disturbance is fairly slow. The step size Δt as chosen by each model is depicted in Figure 3–8 as a function of time, where it is observed that the SFA models use considerably larger step sizes. The numerical performance of the subject models for this slow-transient study is summarized in Table 3–2. As it can be seen in Table 3–2, the numbers of time steps taken by the SFA and TD models differ by two orders of magnitude (157 and 155 steps vs. 16,781 and 16,801 steps). As a result, the average step size utilized by the new CP-SFA model is 25.8 ms, which is over one hundred times larger than that of the TD models (about 0.24 ms). The CPU times also demonstrate that the SFA models are about 39 to 50 times faster than their corresponding TD models. Moreover, in each case, the model with constant-parameter stator interface is about two times faster than its variable-parameter counterpart. Specifically, the proposed CP-SFA model took 13.09 ms, which demonstrates its computational advantage over the previously established VP-SFA model (which took 31.88 ms) as well as the conventional time-domain models (655.01 and 1251.73 ms). Finally, as seen in Table 3–2, the current asi is predicted with very little error by all models (less than 0.1%). 76 Figure 3–6. Slow transient response of the system to a torque change as predicted by the subject models. Figure 3–7. Magnified view of current asi for the slow transient study depicted in Figure 3–6. 77 Figure 3–8. Step size Δt as chosen by the subject models for the slow transient study. Table 3–2. Simulation Efficiency and Accuracy for the Slow-transient Study. Model No. of Time Steps CPU Time, ms 2-Norm Error of asi , % VP-TD 16,781 1,251.73 0.0193 CP-TD 16,801 655.01 0.0272 VP-SFA 157 31.88 0.0899 CP-SFA 155 13.09 0.0869 78 State-Space VBR Modelling of Induction Machines Based on SFA-Type DPs Similar to synchronous machines, an SFA-based model of induction machines [54] has been proposed for EMTP-type solution using the so-called dimension reduction technique. A major constraint of SFA application in EMTP-type modelling, however, is that the main conductance matrix is directly formulated with the step size. Accordingly, varying the step size to reflect the active modes at a given point of simulation (as is particularly useful with DPs) requires re-discretizing all components and reformulating/re-factorizing the conductance matrix, which gives rise to very costly computations. This section proposes an SFA-based VBR modelling of induction machines using a state-space formulation. This proposed model is a complement to the prior EMTP-type model [54], and can be readily implemented in SV-based programs and solved using built-in variable-step solvers. Consequently, the step size is adjusted automatically “on-the-fly”, reflecting the active modes without additional inputs from the user. Time-Domain VBR Induction Machine Model (VBR-TD-IM) Without loss of generality, a general-purpose lumped-parameter three-phase symmetrical induction machine model is considered here. All rotor variables are referred to the stator side and motor sign convention is used; the dynamics of the mechanical system are similar to (3–1) and (3–2) as used for synchronous machines. The VBR formulation of the induction machine ([94], see model VBR-III) is considered, which assumes ungrounded stator windings and results in an advantageous structure with diagonal stator resistance and inductance matrices. Additional information regarding other VBR formulations considering grounded stator windings can be found in [94]. Herein, the stator voltage equation is represented as 79 ''abcsabcsDabcsDabcsdtdeiLiRv , (3–43) where DDDD rrrdiag , ,R , (3–44) DDDD LLLdiag ,,L , (3–45) with rlrmsD rLLrr22 , mlsD LLL , (3–46) and 111lrmmLLL . (3–47) The subtransient voltages abce in (3–43) are defined as Tdqsabcs ee 01 Ke , (3–48) where sK is Park’s transformation matrix [84] and qrlrqrmlrrmlrdrmrqLLLrLLLe 2 , (3–49) drlrdrmlrrmlrqrmrdLLLrLLLe 2 . (3–50) The rotor state equations are expressed as drrmqqrlrrqrLrdtd , (3–51) qrrmddrlrrdrLrdtd , (3–52) where is the speed of the reference frame and 80 lrqrqsmmqLiL , (3–53) lrdrdsmmdLiL . (3–54) The electromagnetic torque is calculated the same way as (3–16) VBR Induction Machine Model Based on SFA 3.2.2.1 SFA-Type DP Stator Interface The VBR-TD-IM model formulates the stator circuit with abc-phase variables – an advantageous property for direct interfacing with external inductive networks. This special structure also allows the stator-network interface to be formulated in SFA-type DP form. First, the stator variables in (3–43) are rewritten in the form of (2–27) as tωtω sabcsQsabcsIabcs sincos fff , (3–55) where the variable “ f ” may represent stator voltages or currents, or subtransient voltages. This transforms (3–43) into ].sincos[]cossin[]sincos[]sincos[sincostωtωtωtωωtωptωptωtωtωtωs''abcsQs''abcsIsabcsQsabcsIDssabcsQsabcsIDsabcsQsabcsIDsabcsQsabcsIEEiiLiiLiiRvv (3–56) The Hilbert transform is then applied to (3-56), which yields 81 ].cossin[]sincos[]cossin[]cossin[cossintωtωtωtωωtωptωptωtωtωtωs''abcsQs''abcsIsabcsQsabcsIDssabcsQsabcsIDsabcsQsabcsIDsabcsQsabcsIEEiiLiiLiiRvv (3–57) Finally, combining (3-56) and (3-57) into the form of (2–33), the VBR stator equation in DP form can be expressed as ''abcsabcsDabcsDsDabcsdtdj EILILRV ][ . (3–58) To facilitate complex variable operations using regular SV-based simulation programs, (3–58) can be further decomposed into real and imaginary parts as '' ,'' , , , , , , ,iabcsrabcsiabcsrabcsSFAiabcsrabcsSFAiabcsrabcsdtdEEIIRIILVV, (3–59) and the SFA-based equivalent subtransient inductance and resistance matrices SFAL and SFAR are DDSFALLL00, (3–60) DDsDsDSFARLLRR. (3–61) 3.2.2.2 Rotor Subsystem and Its Interface As can be observed from (3–51) and (3–52), the qd-coordinate rotor variables (rotor flux linkages, field currents, etc.) in the VBR-TD-IM model vary slowly during electromechanical transients if the synchronous reference frame is used. It is therefore convenient to keep the rotor 82 subsystem in qd coordinates (using the synchronous reference frame), since no frequency shifting is required. However, the VBR model structure requires transforming the time-domain qd-coordinate rotor variables to abc-phase DPs. To obtain such transformation, the time-domain subtransient voltage ''ase is first recalled from (3-48) as tetee sdsqas sincos'''''' , (3–62) Where sst is the angle between the synchronous reference frame qd0 and the stationary reference frame abc. It is observed that (3-62) is identical to (2–27) as the combination of two orthogonal low-pass signals ''qe and ''de modulated by tscos and tssin , respectively. The Hilbert transform is then applied to (3-62), which yields teteeH sdsqas cossin]['''''' . (3–63) Combining (3-62) and (3-63) into the form of (2–33), the subtransient voltage ''asE in DP form is obtained as '''''' dqas jeeE . (3–64) Applying the same approach to ''bse and ''cse , the subtransient voltage source ''abcsE in DPs can be obtained as TdqabcUqdsuabcs ee ]['''',,''KE , (3–65) where the transformation matrix abcU qdsu,,K is 83 Figure 3–9. Block diagram for implementation of the VBR-SFA-IM model. )322()32()322()32()2()0(,,jjjjjjabcUqdsueeeeeeK . (3–66) It is noted that the synchronous reference frame is considered here, which renders (3-66) to be constant. The same approach, however, can be applied to an arbitrary reference frame, at the expense of making (3-66) time-dependent [e.g., see (3–39) and (3-40)]. Conversely, the transformation for stator current from SFA-type DPs abcsI to instantaneous signals abcsi , can be readily obtained by applying (2–34) as . ][0tjabcsssqdseIKi (3–67) In summary, the proposed SFA-based VBR induction machine model (herein referred to as VBR-SFA-IM) consists of the SFA-type DP form stator-network interface (3–59), the rotor-stator variable transformations (3–65) and (3–67), the subtransient voltages (3–49) and (3–50), and the rotor subsystem equations (3–51) and (3–52). The detailed implementation of the proposed VBR-SFA-IM model is depicted using a block diagram in Figure 3–9. 84 Computer Studies The induction machine is interfaced to an infinite bus to focus the investigation on the numerical properties of the subject models. The machine parameters are summarized in the Appendix C. The subject models (i.e., the VBR-TD-IM and VBR-SFA-IM models) have been implemented in MATLAB/Simulink [27]-[28] and the synchronous reference frame is used for all models. The simulation is executed on a PC with a 3.40GHz Intel i7-2600 CPU and 8 GB of RAM. To emulate a general-purpose electromechanical transient situation which spans the entire speed range (from zero to synchronous speed), a case study comprising the start-up and load-change transients is considered. The induction machine is initially freely accelerated from stall; then, at t = 3.0 s, a load torque mT = 1980 N·m is applied; the simulation is then run until t = 5.0 s to reach a steady-state condition. A reference solution is produced by the conventional qd model solved using the ODE4 solver with a very small integration step of 1 μs. Due to the numerical stability limitations of explicit solvers (such as ode45) for the selection of step size [28], the implicit ode23tb solver is herein used with the following settings: relative and absolute error tolerances of 10−4, and maximum and minimum step sizes of 0.2 s (i.e., twelve cycles at 60 Hz) and 0.1 µs, respectively. Due to space constraints, only the transient responses of the stator current asi , the electromagnetic torque eT , and the rotor speed r by the subject models are depicted in Figure 3–10, where no visual difference is observed for eT and r . Two fragments of Figure 3–10 showing asi during the start-up and torque-change transients are magnified and shown in Figure 3–11. As can be seen in Figure 3–10 and Figure 3–11, the VBR-SFA-IM model yields the DP solution asI [VBR-SFA-IM (DP), see line c)], i.e., the envelope of the time-domain waveform asi . 85 To adequately compare the details, the trajectory of asI has also been converted using (7) to the time-domain values asi as shown using circles in Figure 3–11 [VBR-SFA-IM, see line d)]. As shown in Figure 3–11 (a), during the start-up transient both the VBR-TD-IM and VBR-SFA-IM models use relatively small step sizes and closely match the reference, which validates the numerical consistency of the subject models. The trajectory of asI exhibit 60 Hz oscillations, which result from the decaying DC component in stator current asi . It is verified that when deviating largely from nominal speed, the VBR-SFA-IM model is able to give an accurate solution similar to the VBR-TD-IM model, provided that the simulation step sizes are sufficiently small. After the induction machine has reached close to nominal speed, the subject models choose distinct step sizes during the load-change transient as shown in Figure 3–11 (b). In particular, the VBR-SFA-IM model is able to correctly simulate the DP asI with noticeably large time step sizes, since the envelope of asi varies slowly. This is unlike the VBR-TD-IM model, for which small step sizes are still required to accurately track the fundamental 60 Hz sinusoidal waveforms of asi . To give a more comprehensive insight into the simulation behaviour of the subject models, the step size Δt is depicted in Figure 3–12 as a function of time. As can be seen in Figure 3–12, the step sizes of the VBR-TD-IM model remain relatively small during the whole simulation period. In contrast, the step sizes of the VBR-SFA-IM model begin to increase rapidly when the induction machine has reached close to nominal speed [Figure 3–12, see after t = 2 s]. In particular, it is observed that the VBR-SFA-IM model uses a maximum step size of 0.2 s in steady-state conditions, whereas with the VBR-TD-IM model the step sizes remain fairly small (less than 0.5 ms). 86 A quantitative assessment of the numerical performance of the subject models is summarized in Table 3–3. As can be seen in Table 3–3, the VBR-SFA-IM model requires almost 9 times fewer time steps than the VBR-TD-IM model (2216 vs. 19846 steps) due to the use of DPs. Consequently, the VBR-SFA-IM model outperforms the VBR-TD-IM model in terms of CPU time by about 7.76 times (84.1 vs. 652.9 ms). To better quantify numerical accuracy, the 2-norm relative error [121] of the stator current asi is calculated using (3–42) as 100%~~)(22 asasasasiiii , (3–68) where asi~ denotes the reference solution trajectory. For the VBR-SFA-IM model, the instantaneous values of asi are converted from the corresponding DP asI using (2–34). It is shown that despite requiring noticeably fewer time steps and less CPU time, the VBR-SFA-IM model has a slightly larger error than the VBR-TD-IM model (0.26% vs. 0.12%), which does remain in a highly acceptable range (less than 0.5%). 87 Figure 3–10. Start-up and load-change transient responses as predicted by the subject models. 88 Figure 3–11. Magnified view of stator current asi as depicted in Figure 3–10 during: (a) start-up transient; (b) torque-change transient. Figure 3–12. Step size Δt as chosen by the subject models. Table 3–3. Simulation Efficiency and Accuracy of the Subject Models Model No. of Time Steps CPU Time, ms 2-Norm Error of asi , % VBR-TD-IM 19846 652.9 0.1168 VBR-SFA-IM 2216 84.1 0.2579 89 Modelling Line-Commutated Rectifier Systems Including Harmonics for Various Operating Modes Based on the Generalized Average Method This chapter presents the numerically accurate and efficient modelling of LCR systems based on the GAM-type DPs. It is recalled in Section 1.2.2.2 that the bottleneck for accurate modelling of LCR systems is threefold: 1) the existence of various system operating modes, 2) the representation of harmonic dynamics, and 3) the inclusion of thyristor-controlled operation. Moreover, the use of (detailed, analytical, or parametric) modelling techniques can result in distinct formulations and thus different numerical properties (accuracy, efficiency, etc.). We begin by reviewing the state-of-the-art analytical DP (ADP) models of LCRs, which relate the ac/dc subsystems through complicated switch functions. Then, we propose a new parametric DP (PDP) model of diode LCRs, where the DP dynamics of rectifier/dc-link are represented using a set of explicit algebraic functions that are numerically established. Next, this PDP modelling is extended to thyristor-controlled LCR systems. Computer studies validate the proposed PDP methodology in accurately predicting the steady-state and transient responses of LCR systems under a wide range of operating modes, while highlighting its computational advantages over the conventional detailed model and the established ADP models. 90 GAM-Type DP Modelling of Diode LCR Systems with Harmonics Using Analytical and Parametric Approaches For consistency with previous publications [115]−[116], the same benchmark diode LCR system as depicted in Figure 4–1 is considered in this section. This system has been validated with detailed models, dynamic average-value models, and hardware [109]−[117]. Therein, the three-phase diode rectifier is assumed to be supplied by the ac network as represented by its Thevenin equivalent voltages abcse , series resistance thr and inductance thL , respectively. Since the presence of shunt filters does not significantly affect the modelling of such rectifier systems [115]−[116], only series filters are considered. The optional ac-side series filter is denoted by acr and acL ; and the optional dc filter consists of dcr , dcL and dcC . Similar to [115]−[118], the dc subsystem can be represented by an equivalent resistive load loadR that consumes the needed amount of real power. Finally, the combined equivalent series impedance of the ac network can be further simplified as sr and sL , respectively. Figure 4–1. Typical configuration of the three-phase six-pulse LCR system. 91 Depending on the application, this LCR benchmark system can be readily modified/extended to represent various ac/dc energy conversion systems. For example, if representing a synchronous-machine-rectifier configuration, the benchmark system of Figure 4–1 can be fed from a machine represented in the VBR formulation [86] (with possible coupled/variable inductances in the equivalent ac network) or as a constant-parameter interfacing circuit [117]. In other low-power applications, such as variable frequency drives [97]-[99] or battery charging systems [100], the dc filter is typically simplified to only the dc capacitor for voltage smoothing. On the basis of loading condition (from open circuit to short circuit) and resulting switching patterns, one discontinuous conduction mode (DCM) and three continuous conduction modes (CCM-1, CCM-2, and CCM-3) can be defined for the LCR system [115]. In particular, when the LCR is in the DCM, which can be frequently encountered at light load conditions [97]−[99], the phase current waveforms (see Fig. 2, [118]) can be highly distorted with considerable harmonics (5th, 7th, etc.). Thereafter, for DP modelling of LCR systems, it is desirable to include the dominant 1st, 5th, and 7th order DPs for ac system variables; and for the dc subsystem, the 0th order DP is typically sufficient due to existence of the dc filter. Time-Domain Dynamics of LCR Systems and Switch Functions To develop the GAM-type DP models of power converters, the modulation theory has been traditionally used [73]−[77], which regards the LCR as a modulator of voltages/currents. Using switch functions [122]−[123], the time-domain currents on the ac- and dc-sides are related as iabcdcabcs i si . (4–1) Similarly, the voltages are related as 92 Figure 4–2. Time-domain waveforms of phase a voltage and current switch functions: (a) vas of SF-1, (b) vas of SF-2, (c) ias of SF-1, and (d) ias of SF-2. vabcTabcsvccsvbbsvaasdc svsvsvv sv . (4–2) Here, the vectors abcsv , abcsi denote the voltages and currents of the ac-side; and dcv , dci are the dc-side voltages and currents, respectively, as shown in Figure 4–1. Also, vabcs and iabcs are the time-domain voltage and current switch functions, respectively. Despite the general modulation form (4–1)−(4–2), the switch functions vabcs and iabcs can be developed exploiting different approximation approaches. Two types of switch functions, i.e., the piecewise-linear- and Fourier-series-approximated switch functions (herein referred to as SF- 93 1 and SF-2, respectively) have been investigated in the literature [73]−[77], [122]−[123]. For illustration, the time-domain waveforms of phase a voltage/current switch functions using SF-1 and SF-2 are depicted in Figure 4–2. The switch functions of phase b and c can be obtained by shifting phase a waveforms by 32 . 4.1.1.1 Piecewise-Linear-Approximated Switch Functions (SF-1) Traditionally, the switch functions have been developed to emulate the switching of each diode [73]−[74]. When the commutating inductance is considered, the dc voltage waveform will have steps, while currents can only commutate from one phase to the other gradually [122]−[123], which yields otherwise. ,0)1(332,)1(332 ,21)1(332,)1(3 ,1)1(3,)1(3 ,21)(000000,mmmmmmts mv (4–3) otherwise. ,0)1(332,)1(332,)1(3321 )1(332,)1(3 ,1)1(3,)1(3 ,)1(3)(00000000,mmmtmmmmmtts mi (4–4) 94 Here, mvs , and mis , (m=1, 2, … 6) are the time-domain voltage and current switch functions corresponding to the rectifier diodes S1 ~ S6 as shown in Figure 4–1. The starting instance of the commutation is denoted by 0 [See Figure 4–2, 230 when compared to SF-2]. To include the effect of the commutating inductance, the commutation angle μ (in both SF-1 and SF-2) is calculated as )21(cos 1linedcssViL . (4–5) Also, the patterns in (4–3)−(4–4) are periodical for every switching period sT . Based on the switch functions for each diode, the SF-1 for the LCR are obtained as [73]−[76] 2,5,6,3,4,1,---vvvcvvvbvvvasssssssss;2,5,6,3,4,1,---iiiciiibiiiasssssssss. (4–6) 4.1.1.2 Fourier-Series-Approximated Switch Functions (SF-2) In contrast to SF-1, an alternative type of switch functions can be expressed as the summation of a number of Fourier series as [77], [123] )32()32()()cos()2cos(max1tstststnnAsvcvbvannnva, (4–7) 95 )32()32()()cos(2)2sin(max1tstststnnnAsicibiannnia. (4–8) Here, the highest order of Fourier series is selected as maxn ( maxn ≥ 7) in order to adequately represent the system dynamics, and 6cos2sin4 nnnAn . (4–9) Analytical DP Modelling Next, the analytical DP models of LCRs can be derived by transforming the time-domain ac-dc relationships (4–1)−(4–2) into DPs following the GAM method, which yields Kkiikiabcdckiabcdckabcs ,0ssi , (4–10) , Re2 ***0000Kkkvckcskvbkbskvakasvccsvbbsvaasdcsvsvsvsvsvsvv (4–11) where K={1, 5, 7, …} denotes the order of corresponding DPs. 4.1.2.1 Analytical DP Model Using SF-1 (ADP-1) Nevertheless, due to the piecewise-linear properties of SF-1, developing the respective DPs of (4–3)−(4–4) can be highly complicated. For example, the 1st order DPs of the diode switch functions are obtained as [74] 96 )(21)(4)]1(332[)]1(332[)]1(3[)]1(3[1,)]1(332[)]1(332[)]1(3[)]1(3[1,00000000mjmjmjmjmimjmjmjmjmveeeeseeeejs. (4–12) Accordingly, in [73]−[76] only the 1st order DPs using SF-1 have been derived, and (4–6) is transformed into 12,15,116,13,114,11,1vvvcvvvbvvvasssssssss; 12,15,116,13,114,11,1iiiciiibiiiasssssssss. (4–13) 4.1.2.2 Analytical DP Model Using SF-2 (ADP-2) As seen from (4–12)−(4–13), the ADP-1 model yields complicated expressions even for the fundamental ac components. Moreover, truncating high-order harmonics may not be sufficient for the purpose of system-level transient studies or power quality analysis. Alternatively, the Fourier-series-approximated SF-2 allows to transform (4–7)−(4–8) into DPs as [77] ,,2cos2)32()32()0(KkesessekAskjkvckjkvbkvajkkva (4–14) 97 .,2sin)32()32()0(KkesessekkAskjkickjkibkvajkkia (4–15) In summary, for analytical DP modelling of LCRs, both the ADP-1 and ADP-2 models share the same ac-dc equations (4–10)–(4–11), with the difference being the DP-form switch functions, i.e., (4–12)–(4–13) for the ADP-1 model, and (4–14)–(4–15) for the ADP-2 model, respectively. Parametric DP Modelling (PDP) As shown in Section 4.1.2, in order to analytically capture the ac- and dc-side dynamics of LCRs, complicated expressions are inevitably used in both the ADP-1 and ADP-2 models, which degrade the numerical efficiency of simulations. Moreover, since the switch functions SF-1 and SF-2 were derived upon a single switching pattern, the established ADP models are valid for only one mode of operation, i.e., the CCM-1 mode. A possible solution is to derive several ADP models, and to switch between them as the system operating condition evolves, which however adds to the modelling complexity. To overcome these challenges, this paper extends the parametric approach [112]−[113] to effectively model the ac and dc DP dynamics of the LCR system under various operating conditions. In traditional PAVMs [112]−[118], the average values of ac variables (expressed in the qd synchronous reference frame) and dc variables are related through algebraic parametric functions that are numerically established from detailed simulations over a range of operating conditions. Similarly, assuming a balanced system where the magnitudes of three-phase DPs are identical (i.e., kcskbskasfff , where f = {v, i}), the ac- and dc-side DPs can be related 98 through a set of non-linear algebraic functions. Specifically, this yields the ac- and dc-side DP relations as Kkvv dckkas ,)( 0 , (4–16) 110)( asdc ii . (4–17) Here, kα and 1 denote the parametric functions describing the switching behaviour of the LCR system that are determined by the present operating condition. To simplify the modelling, it is assumed that only the fundamental ac currents affect the predicted dc current waveform. To account for the phase values of DPs, the angular differences between the phase a voltages and its fundamental current are expressed as Kkiv askask ),ang()ang()( 1 . (4–18) Similarly, the DP angles for phase b and c can be obtained by a phase shift of 32 k , respectively. Therefore, combining (4-16) and (4-18), the three-phase ac voltage DPs are obtained as .,)()32()32()()ang(01Kkevevvevvkjkcskjkbskasijdckkaskas (4–19) Finally, to formulate the input of the parametric look-up tables of kα , 1 , and kφ , the operating condition is defined in terms of the dynamic impedance of the LCR as .10asdcivz (4–20) 99 In summary, the proposed parametric DP model (herein referred to as PDP) consists of the algebraic equations (4–17), (4–19), and (4–20) that relate the ac- and dc-side DPs, and the parametric functions of zkα , z1 , and zkφ . Constructing Parametric Functions The parametric functions zkα , z1 , and zkφ can be readily established using the detailed simulation of the LCR system under study, using approaches similar to those in [112]−[113]. Specifically, this is done by examining the desired Fourier series of the system ac/dc variables corresponding to a wide range of operating conditions. Since in steady state the GAM-type DPs become constant values [62], a straightforward approach is to run detailed simulations in multiple steady-state conditions corresponding to the operating points of interest, and then tabulate the fast Fourier transform (FFT) results of the time-domain waveforms of system variables. Typically [115]−[110], a few data points at different operating conditions with subsequent proper curve fitting will suffice for decent accuracy, which however requires multiple runs of detailed simulations. This method of constructing parametric functions is herein referred to as the multi-steady-states (MSS) approach. Alternatively, a faster procedure can be used to generate the desired parametric functions using a single large-signal transient study that spans a wide range of operating conditions [113]. To apply this approach for the purpose of GAM-type DPs, the sliding window FFT of time-domain waveforms is used. To emulate a wide range of loading conditions, it is possible to employ a variable load resistance that increases exponentially as tload rRtR )1()( 0 , (4–21) 100 Here, 0R is the initial load of a small value, and r denotes the exponential growth rate. To better illustrate this, for the considered LCR system with parameters summarized in Appendix C, the large-signal transient response to an exponential increase of load resistance ( 0R = 1 , r = 85) is shown in Figure 4–3. The corresponding DPs are extracted by running the SWFFT of the time-domain waveforms from detailed simulation using MATLAB [27]. This method is herein referred to as the large-signal-transient (LST) approach. For comparison purposes, the parametric functions zkα , z1 , and zkφ (K={1, 5, 7}) extracted using the MSS and the LST approaches are depicted in Figure 4–4 and Figure 4–5. To achieve a reasonable representation of parametric functions, the MSS approach with 32 steady-state solution points has been used. It is noted that the data points are more tightly-spaced in the region where parametric functions are highly nonlinear. These data points are then stored in a look-up table, where interpolation/extrapolation may be needed if necessary. Alternatively, it is shown in Figure 4–4 and Figure 4–5 that the LST approach yields parametric functions that match perfectly with the steady-state solution points, which validates their mathematical equivalence to the MSS-generated parametric functions. Despite featuring similar accuracy in various operation modes, the MSS and LST approaches can differ greatly in terms of computations required to generate the parametric functions. The MSS approach can be implemented by running detailed simulations in a loop to calculate the steady state solutions for the desired number of points, which can take several minutes of total CPU time overall. In contrast, the LST approach only requires a one-run simulation to produce the desired parametric functions, which may take only a few seconds. 101 Figure 4–3. The large-signal transient response to an exponential increase of load resistance for calculating parametric functions. 102 Figure 4–4. Numerically calculated functions )(1 z , )(1 z , and )(1 z as obtained from detailed simulations using (a) the MSS approach; and (b) the LST approach. Figure 4–5. Numerically calculated functions )(7,5 zα and )(7,5 zφ as obtained from detailed simulations using (a) the MSS approach; and (b) the LST approach. 103 Implementation The ADP-1, ADP-2 and the proposed PDP models can be implemented using the block diagrams depicted in Figure 4–6. As shown in Figure 4–6, the overall system consists of three parts: the ac system comprising a number of subsystems for the fundamental frequency and each of the desired harmonics; the LCR model; and the dc subsystem. For ADP models, the ac subsystems can be formulated as an algebraic equation, since the combined network inductance sL has been included into the commutation angle μ [see (4–5)]. For the PDP model, the ac network shown in Figure 4–1 has a convenient voltage-behind-impedance formulation that is governed by the time-domain equation abcsabcsabcssabcss rdtdL evii . (4–22) Therefore, applying the GAM and (5) yields the kth order DP equation as kabcskabcskabcsssskabcssLjkrdtdL evii )( , (4–23) which is then implemented in each of the ac subsystems as shown in Figure 4–6 (b). For the dc subsystem, since the considered 0th order DP is equivalent to averaging time-domain signals [i.e., k = 0], the DP model is directly interfaced with the dc subsystem expressed in the time domain. It is also noted that, compared with ADP models, the PDP model possesses reversed inputs and outputs as shown in Figure 4–6 (b), which is preferable for interfacing considerations [112]. In this case, the dc filter with the PDP model can be modelled using a proper transfer function [112]. Here the control of the LCR system has not been considered since only passive diodes are used. However, many applications of integrated ac/dc systems will include feedback control such as speed and voltage regulation of rotating machines, active and reactive power control in 104 distributed generation, etc. In the DP modelling of control systems, one can choose: 1) to directly model all control signals in DPs; or 2) to selectively model some control signals in the time domain [75], since many control reference signals (magnitudes, angles, etc.) may be close to dc in steady state (e.g., voltage regulation in HVDC systems [75]-[77]). It is also worth noting that for the ADP-1, ADP-2, and the proposed PDP models, only the balanced system configuration and symmetric (healthy) operations of LCRs are considered. For unbalanced and/or fault conditions of the LCR system (e.g., due to unbalanced loads, asymmetric faults on ac lines, internal faults of LCRs, etc.), there will exist some non-characteristic ac harmonics (2nd, 3rd, 4th, etc.) [98], [102], [130], which requires additional order of DPs for accurate modelling. Moreover, since the unbalanced system and/or severe harmonics can also affect switching patterns of LCRs [130], the proposed PDP model, which assumes only symmetric operations, may have limited accuracy in predicting the actual waveforms. It has been shown in the literature that the parametric modelling approach can achieve a decent accuracy even for LCRs with unbalanced ac network [115], [131]. However, to consider severe asymmetric conditions (e.g., internal fault), it would be required to first decompose the system variables into both positive and negative sequences [131], [132], and then properly use switch/parametric functions to relate the ac-dc harmonic dynamics. However, for the scope of research in this thesis, the balanced system conditions are deemed acceptable. 105 Figure 4–6. Implementation of the DP models: (a) ADP-1 and ADP-2 models, (b) the proposed PDP model. Computer Studies To analyze the models discussed in this section, computer studies are carried out on the benchmark LCR system shown in Figure 4–1, with parameters summarized in Appendix D. The DP models (i.e., the ADP-1 [73]−[74], ADP-2 [77], and the proposed PDP model) have been implemented in MATLAB/Simulink [27]−[28]. For the purpose of reference, the detailed switch-level model has been implemented using PLECS Blockset [29] within the Simulink Environment. To demonstrate the proposed approach, only the 1st, 5th, and 7th order DPs of the ac-side are 106 considered in this paper. However, with its simple structure, the proposed PDP model can be readily augmented to include higher-order (11th, 13th, etc.) harmonics of interest, provided that the dynamics of these harmonics are desired. To emulate a situation when the LCR spans various operation modes, a case study comprising the transition from the normal operation in CCM-1 to a light load in DCM is investigated here. The system is assumed to initially operate in CCM-1 with a dc load 150loadR . At t = 2.0 s, the system is switched to a very light load 950loadR ; and the simulation is run until t = 4.0 s. To investigate the numerical accuracy of the subject models in different operation modes, the models are first solved with small step sizes using the ode23tb solver and the following settings: relative and absolute error tolerances of 10−3, and maximum and minimum step sizes of 100µs and 0.1 µs, respectively. For consistency, all case studies are conducted on a PC with a 3.40 GHz Intel i7-2600 CPU. 4.1.6.1 Steady State Analysis in CCM-1 and DCM To evaluate the properties of the subject models in terms of capturing harmonics, it is instructive to consider their performance in steady states first in CCM-1 and DCM. Due to space constraints, only the responses of the phase a and dc voltages and currents, as predicted by the subject models, are depicted in Figure 4–7 and Figure 4–9. The corresponding harmonic contents are shown in Figure 4–8 and Figure 4–10, respectively. Since the simulations of DP models are carried out in the DP domains, to adequately compare with the detailed reference solution, the resulting DP trajectories kasi have been converted back into time-domain instantaneous values asi using (2–37) (the instantaneous reconstruction, [124]). As can be seen in Figure 4–7, in CCM-1 all DP models produce reasonably accurate results on the dc-side compared to the detailed 107 solution. For the ac side, the ADP-1 model only predicts the fundamental frequency response, and the ADP-2 model yields the asi trajectory similar to the SF-2 waveform shown in Figure 4–2. Furthermore, it is observed that the proposed PDP model, by including both 5th and 7th harmonics, is capable of producing results that match well to the detailed model. This is validated in Figure 4–9, which shows the steady-state responses in DCM. As shown in Figure 4–9, in DCM the dc current drops to zero and thus causes significant distortions in the asi waveform (see the double peaks with zero in between). However, due to the analytical derivation based on a single switching pattern, both the ADP-1 and ADP-2 models fail to reproduce the distorted waveforms in DCM [see Figure 4–9, lines b) and c)]. Nevertheless, the PDP model still matches the detailed simulation very well. To provide a comprehensive insight into the accuracy of harmonic prediction by the subject models, the extracted harmonic contents of the phase a voltage and current in CCM-1 (Figure 4–7) and in DCM (Figure 4–9) are shown in Figure 4–8 and Figure 4–10, respectively. As seen in Figure 4–8 and Figure 4–10, the ADP models result in less accurate predictions for the considered harmonics in both CCM-1 and DCM. Specifically, the ADP-1 model captures only the fundamental frequency component; and the ADP-2 model predicts 5th and 7th harmonics of current asi with noticeable error even in CCM-1. This error is due to the switch function ias [see Figure 4–2 (d)] being different from the actual current waveform, which is another limitation of the ADP methods. In contrast, it is observed that harmonics predicted by the PDP model match the detailed simulation very well. 108 Figure 4–7. Ac and dc system variables for steady-state in CCM-1 as predicted by the subject models. Figure 4–8. Harmonic content of the ac-side phase a current and voltage for the considered operating point in CCM-1. 109 Figure 4–9. Ac and dc system variables for steady-state response in DCM as predicted by the subject models. Figure 4–10. Harmonic content of the ac-side phase a current and voltage for the considered operating point in DCM. 110 4.1.6.2 Transient Response from CCM-1 to DCM Next, a fragment of study is shown in Figure 4–11 to demonstrate the predicted trajectories of asv , asi , dcv , and dci by the subject models during the load-change transient. On one hand, it is shown in Figure 4–11 that the ADP-1 and ADP-2 models do not trace the dynamics very well and predict lower dc voltage in the new DCM steady state (which is expected since they have been derived for only CCM-1). On the other hand, the transient response predicted by the PDP model is highly consistent with the detailed solution. This large-signal transient response verifies the superior accuracy of the proposed PDP model over the established ADP models. Figure 4–11. Ac and dc system variables for the load-change transient as predicted by the subject models. 4.1.6.3 Predicting THD in Wide Range of Operating Conditions To investigate the modelling accuracy of LCRs under a wide range of operating conditions, the total harmonic distortion (THD) level is considered [97], which is defined as 111 Figure 4–12. THD of the phase a voltage and current as predicted by the subject models over a wide range of operating conditions. 1,12%FFTHDKkkk , (4–24) where kF denotes the rms value of the k th order harmonic of voltage or current. The THD of phase a voltage and current, as predicted by the subject models for the load ranging from open circuit to short circuit, is calculated and plotted in Figure 4–12. As can be seen in Figure 4–12, high current distortion is observed at light loads due to the discontinuous conduction, and the voltage is more distorted at heavy loads when the current waveform is close to sinusoidal. The ADP-1 model is not included here since it only predicts the fundamental frequency response. The THD results predicted by the ADP-2 model are considerably less accurate due to its analytical derivation and fixed waveform of the switching functions, especially at light loading conditions (which correspond to DCM). Whereas, the PDP model yields THD results very close to the detailed simulation over a wide region of operating conditions. 112 4.1.6.4 Flexible Step Size Study To further demonstrate the advantageous properties of the proposed PDP model, the same 4-second transient study is solved using the same error tolerances, but the maximum step size is increased to 0.1s (six cycles at 60Hz), which allows the solver to adaptively select step sizes “on-the-fly” in a wider range. To study the numerical efficiency, the step size Δt as chosen by each model is depicted in Figure 4–13 as a function of time. As seen in Figure 4–13, the step size of the detailed model changes (due to switching events) but remains fairly small during the entire simulation period. In contrast, the DP models operate at much larger steps Δt (with the maximum step size of 0.1s in steady states). When the load is changed, the system states vary accordingly, which results in the use of small time steps to trace the transient response until the LCR system reaches a new steady state. These observations are verified in Figure 4–14 and Figure 4–15(a), where the responses of asi as predicted by the subject models are shown in steady state CCM-1 and during the load-change transient, respectively. Therein, the trajectory of current asi converted from DP solutions are shown using circles and asterisks [see lines b), d) and f)]. Additionally, two times the magnitude of the 1st order DPs, i.e., 12 asi is shown using the dashed-dotted lines and labeled as 1st order DPs [see lines c), e) and g)]. As can be seen in Figure 4–14 and Figure 4–15(a), the trajectories of 12 asi essentially represent the envelopes of the fundamental frequency waveforms, and thus can be simulated using larger step sizes. It is also noted in Figure 4–13 that after the load-change instance, the step sizes of the ADP models increase more quickly than those of the PDP model. This is due to the different 113 representations of ac systems by the ADP and PDP models. Specifically, for ADP models the ac network inductance Ls is only included into the commutation angle μ [see (4–5)], which simplifies the ac subsystems [as shown in Figure 4–6 (a)] to be purely algebraic (i.e., resulting in a reduced-order representation of the ac network), thus omitting the detailed fast transients. Whereas, the PDP model represents each of the ac DPs as a dynamic state using (21). To better illustrate this point, the trajectories of 52 asi and 72 asi during the load-change transient, i.e., the envelopes of the 5th and 7th harmonics, are plotted in Figure 4–15 (b) and (c), respectively. As can be seen in Figure 4–15 (b) and (c), the ADP-2 model gives the damped dynamics of the 5th and 7th DPs, since it formulates the ac network in algebraic equations (reduced-order representations). In contrast, as observed after the load change, the response predicted by the PDP model has some ringing in the 5th and 7th DPs, which causes the variable-step solver to use small step sizes during the transient. This DP response with ringing, however, reproduces the accurate time-domain simulation results, as verified in Figure 4–15 (d), which depicts a magnified view of current asi during the load-step-change transient. As shown in Figure 4–15 (d), the ADP models yield inaccurate transient results [see lines b) and c)], while the PDP model reproduces time-domain simulation results that precisely match the detailed solutions [see lines a) and d)]. The origin of the ringing in harmonic DPs is discussed in Section 4.1.6.5. Next, the quantitative assessment of simulation efficiency is summarized in Table 4–1, wherein it is seen that all DP models show time steps at least three orders of magnitude less than the detailed model, which results in the drastic reduction of total CPU time. Due to this difference in time steps (modelling accuracy) by the subject models, the CPU time per step is also included in the last column of Table 4–1 to allow a more thorough comparison. Therein, it is seen that the proposed PDP model yields significantly less computational time per steps than the ADP models. 114 Figure 4–13. Step size Δt as chosen by each model for the flexible step sizes study using a step change of load resistance. Figure 4–14. Current asi as predicted by the subject models using flexible step sizes for steady-state response in CCM-1. 115 Figure 4–15. System transient response to the step change of load resistance as predicted by the subject models using flexible time steps: (a) current asi , (b) 5th order DP 52 asi , (c) 7th order DP 72 asi , and (d) magnified view from subplot (a). 116 Table 4–1. Models Simulation Efficiency of the Subject Models for the Flexible Step Size Case Study Using a Step Change of Load Resistance. Model No. of Time Steps Total CPU Time, ms CPU Time per Step, µs Detailed 146925 3521.2 23.97 ADP-1 80 15.4 192.5 ADP-2 80 10.5 131.25 PDP 388 23.9 61.60 4.1.6.5 Discussions of the Ringing of Harmonic DPs As noted in Section 4.1.6.4, some ringing is presented in the harmonic envelops/DPs during the fast load transient. Such oscillations have been noticed and meticulously analyzed in literature [124]-[125]. The authors of [124]-[125] explain that a direct step change in the DP domain can in fact introduce additional dynamics/oscillations into the associated harmonic DPs during transients, while precisely reproducing the correct time-domain simulation results. When mapping a step change of a time-domain signal into the DP domain, two different approaches are employed in the literature, i.e., implementing: 1) a set of continuous time-varying changes of DP inputs; or 2) a direct step change in the harmonic DPs [124], which is simpler and has been considered for the purpose of this thesis. However, the direct step changes in the DP domain can cause inconsistency between responses of the DP systems before and after the input change [125], and can provoke oscillations in the harmonic DPs that are damped according to the system’s time constants. This can also be demonstrated by examining the time-domain and DP equations of the ac subsystems, i.e., (4–22)-(4–23). Specifically, with some simplifications, it can be shown that the eigenvalues of (4–23) for the kth order DPs are related to the modes defined by sss jkLr / , which can become excited as shown in Figure 4–15. Finally, it is noted that despite the transient ringing of 117 harmonic DPs, the PDP model reproduces accurate time-domain simulation results as shown in Figure 4–15 (d) (which is also explained in [124]). To verify the possibility of avoiding excitement of the harmonic DP ringing modes, the same case study is performed with the load step change at t = 2s being replaced by a sigmoid function )(11)(ctaetS , (4–25) where the parameters a = 50 and c = 2s are chosen to implement a smooth transition. The resulting step sizes Δt and the system response as produced by the subject models are shown in Figure 4–16 and Figure 4–17, respectively. As can be seen in Figure 4–17, since the sigmoid function implements a smoother “S” shape change in load, the ringing in the 5th and 7th DPs has been removed. This slower transient also allows both the ADP and PDP models to use large step sizes as shown in Figure 4–16. The quantitative assessment of this case study is also summarized in Table 4–2. Therein, it is noted that all DP models yield a similar number of time steps, which is much smaller compared to the detailed switching model. At the same time, the total CPU time taken by the PDP model is now only 5.4 ms (as opposed to 23.9 ms reported in Table 4–1). The results of Table 4–1 and Table 4–2 also suggest that the advantages of the PDP model should be noticeable in fixed time step solution. 118 Figure 4–16. Step size Δt as chosen by each model for the flexible step sizes study using a sigmoid change of load resistance. Figure 4–17. System transient response to the sigmoid change of load resistance as predicted by the subject models using flexible time steps: (a) current asi , (b) 5th order DP 52 asi , and (c) 7th order DP 72 asi . 119 Table 4–2. Simulation Efficiency of the Subject Models for the Flexible Step Size Case Study Using a Sigmoid Change of Load Resistance. Model No. of Time Steps Total CPU Time, ms CPU Time per Step, µs Detailed 142390 3481.3 24.45 ADP-1 55 11.5 209.09 ADP-2 56 9.5 169.64 PDP 59 5.4 91.53 4.1.6.6 Fixed Step Size Study To fairly compare the efficiency of the considered DP models, the same transient study is performed using the fixed step trapezoidal rule (that is widely used in commercial EMTP-type simulation tools [20]−[25]). To demonstrate a sample of results, Figure 4–18 depicts the fragment of current asi transient response to the load step change as predicted by the subject models using a fixed step size Δt = 1 ms. As can be seen in Figure 4–18, the simulation results of the PDP model match the detailed solution very well, even at large time steps, while the ADP models do not capture the details of this transient. This is expected since the ADP models have less accuracy in predicting the transients of harmonics (see Figure 4–15 and Figure 4–18). In addition, Table 4–3 summarizes the total CPU time (in seconds) taken by the subject models using different step sizes ranging from 10 µs to 1ms. As can be seen in Table 4–3, the computational performance of the PDP model is consistent with its CPU time per step observed in Table 4–1 and Table 4–2, and is noticeably faster than the ADP models. 120 Figure 4–18. Current asi transient response to the load step-change as predicted by the subject models using fixed step size Δt = 1ms. Table 4–3. Total CPU Time (in Sec) Taken by the Subject Models Solved Using Fixed Step Trapezoidal Rule Step size Δt, µs 10 100 1,000 ADP-1 58.3539 5.9496 0.5878 ADP-2 44.5843 4.5444 0.4632 PDP 24.3828 2.4981 0.2554 121 Parametric DP Modelling of Thyristor-Controlled LCR Systems Including Harmonics for Various Operating Modes As shown in Section 4.1, the proposed PDP modelling methodology is effective for accurate modelling of diode LCR systems (including ac harmonics) in various operating modes, while providing significant numerical advantages. Therefore, it is desirable to extend this PDP modelling methodology [see Section 4.1.3] to the thyristor-controlled LCR systems that are utilized in many industrial applications [98]-[101]. For this section, the benchmark thyristor-controlled LCR system is depicted in Figure 4–19 [117], where the rectifier can be fed from two cases of ac source: I) a rotating machine, or II) a distribution feeder/transformer (represented in its Thevenin equivalent circuit). In particular, for Case I, the CP-VBR formulation [89]−[93] of rotating machines is preferred, which possesses a direct machine-network interfacing circuit represented as decoupled constant-parameter RL branches. Such an interfacing circuit can be achieved even for salient pole rotor machines [89]−[93]. Accordingly, the thyristor firing pulses can be generated based on the filtered terminal voltages or the sensed rotor position, respectively [117]. Without loss of generality, the thyristor firing control via the filtered terminal voltages abcsv and a firing angle (i.e., Case II) is considered here [117]. Figure 4–19. Benchmark configuration of the thyristor-controlled LCR systems. 122 Relating the AC-DC DP Dynamics As addressed in Section 4.1, the rectifier switching cell in Figure 4–19 can be represented by an algebraic block depicted in Figure 4–20. Therein, assuming a balance system, the ac-dc voltage and current DPs are related through parametric functions kv,ω , iω , and kφ , where kK, and K={1, 5, 7, …} denotes the order of desired DPs. These parametric functions are defined in (4–26)-(4–28), which are determined by the present system operating mode using two inputs: 1) the dynamic impedance dz defined in (4–20) which indicates the loading condition; and 2) the thyristor firing angle . 0,),( dckasdkv vvz , Kk (4–26) 10),( asdcdi iiz , (4–27) )ang()ang(),(1askasdkivz , Kk (4–28) 10 asdcdivz . (4–29) Figure 4–20. Implementation of the PDP model of thyristor-controlled rectifiers. 123 Constructing Parametric Functions The procedure of finding the algebraic relations between ac and dc DPs of the thyristor rectifier are as follows: 1) Select the set of desired harmonics, i.e., K, based on the properties of non-linear components in the system; 2) Run detailed simulations of the subject system for a range of operating conditions, which can be achieved using a myriad of steady-state points at different dz and α, or several large-signal transients that span various modes [see Section 4.1.4]; 3) Extract desired Fourier coefficients (i.e., DPs) from the time-domain waveforms of the system’s ac/dc variables, and establish the parametric functions at each operating point using (4–26)-(4–28) ; 4) Store these functions in two-dimensional look-up tables, and use proper curve fitting or inter/extrapolation if needed. Computer Studies The system of Figure 4–19 with parameters identical to Appendix D and the ac filter same as Eq.(11) [117], is considered for study, wherein the dominant 1st, 5th, and 7th order DPs for ac variables and the 0th order DP for dc subsystem are considered. The detailed model implemented using PLECS Blockset [29] and the analytical DP (ADP-2) model [77] are used for comparison. For the PDP model, the parametric functions (4–26)-(4–28) are constructed following the steps in Section 4.2.2, and are depicted in Figure 4–21. Therein, 2520 steady-state solution points have been used, which are more tightly-spaced in the region where the parametric functions become nonlinear. 124 Figure 4–21. Parametric functions as numerically calculated from detailed simulations: (a) ,1, dv z ; (b) ,di z ; (c) ,1 dz ; (d) ,5, dv z ; (e) ,5 dz ; (f) ,7, dv z ; and (g) ,7 dz . 4.2.3.1 Steady-State Response To consider a wide range of steady state operating conditions, the computed THD [See (4–24)] of ac current asi for different thyristor firing angles with a dc load resistance 25lR is shown in Figure 4–22. As shown in Figure 4–22, the ADP-2 model gives inaccurate THD results 125 Figure 4–22. Steady-state responses: THD % of asi for different thyristor firing angles. at conditions with 40 , which is due to the fixed-shape switch functions [See Section 4.1.1] that are different from the actual current waveforms. In contrast, the THD results are well predicted by the PDP model for the entire firing angle range, which verifies its accuracy in steady state conditions. 4.2.3.2 Large-Signal Transient Response To emulate various large-signal transients, the system (assumed initially in steady-state with 150lR and 10 ) is subjected to several load and firing angle changes: i) at t=0.05s, step change of lR ; ii) at t=0.15s, ramp change of lR ; iii) at t=0.35s, step change of ; and iv) at t=0.45s, ramp change of . In addition, to illustrate the close-loop dynamic responses of the subject models, two thyristor firing control events are included: v) at t=0.65s, activation of the dc voltage controller; and vi) at t=0.75s, step change of the reference voltage. The details of these large-signal transients are also shown in Figure 4–23(a)(b), and the close-loop thyristor firing controller is depicted in Figure 4–26 with parameters summarized in Table 4–4. 126 The resulting system response of dc variables is depicted in Figure 4–23 (c)(d), which shows noticeable error with the ADP-2 model yet excellent agreement between the PDP and detailed model results. This is further confirmed with the magnified views of ac variables asv and asi shown in Figure 4–24. Specifically, it is noted in Figure 4–24 (a) and (d) that as the system enters DCM due to increase of lR or , the ADP-2 model fails to reproduce the distorted current waveforms (dips between double peaks), while the PDP model yields accurate waveforms as the system operating condition evolves [see lines a) and c)]. It is also seen in Figure 4–24 (a) and (c) that during the step-change transients the ADP-2 model can result in dynamic errors due to the reduced-order representation of the ac network [see Section 4.1.6.5]. These errors, however, do not occur with the PDP model. Moreover, Figure 4–25 shows magnified views of system response to the activation/change of the dc voltage controller, which again demonstrates significant error with the ADP-2 model [see line b)] and verifies the conformity between the PDP and detailed models in close-loop dynamic responses. Table 4–4. Parameters of the close-loop thyristor firing controller for dc voltage regulation. Controller Parameters PI Controller 1.0PK , 10IK Low-Pass Filter 11)(2 bsassH 61045 a , 31014 b 127 Figure 4–23. System response to the large-signal transients: (a) profile of lR ; (b) profile of ; (c) resulting dc voltage dcv ; and (d) resulting dc current dci . 128 Figure 4–24. Magnified views of ac variables asv and asi response to: (a) step change of lR ; (b) ramp change of lR ; (c) step change of ; (d) ramp change of . 129 Figure 4–25. Magnified views of system close-loop dynamic response to: (a) activation of the dc voltage controller; (b) step change of the reference voltage. Figure 4–26. The close-loop thyristor firing controller for dc voltage regulation. 130 4.2.3.3 Frequency-Domain Study Finally, to portray the frequency-domain characteristics of the system, a small-signal analysis is performed around the same initial steady-state operating point as in Section 4.2.3.2. The open-loop transfer function )(sF is considered as )(~)(~)(ssvsF dc . (4–30) where )(~ svdc is the change in the output dc voltage due to the small-signal perturbation in the input firing angle )(~ s . The transfer functions as predicted by the subject models are depicted in Figure 4–27, where a frequency-sweep technique has been used for the detailed model. As shown in Figure 4–27, the ADP-2 model shows inaccurate results for frequencies higher than 20 Hz, while the transfer functions evaluated by the detailed and PDP models are well matched, especially in the range from 1 to 200 Hz. In general, due to the Fourier transformation of system variables [i.e., averaging of frequency components, as shown in (2–36)], the DP models should yield accurate frequency-domain response up to about the system switching frequency, e.g., 360 Hz for the considered six-pulse rectifier system. When the frequency approaches close to or beyond the switching frequency, it would be expectant to observe deviation of results between the detailed and DP models, since the switching patterns of detailed model may be changed and the basic assumption of generalized averaging is no longer valid. However, since the input is the firing angle command, lower frequency dynamics are of more significant importance for studying stability and designing the controller of the system. 131 Figure 4–27. Transfer function from the input firing angle to the output dc voltage as predicted by the subject models. 132 Interfacing SFA- and GAM-Type DPs for Modelling the Synchronous Machine-Rectifier Systems As shown in Figure 1–1 and discussed in Section 1.2.2, the synchronous-machine rectifier system is used as a representative sub-circuit/unit of many integrated ac-dc power systems. This chapter presents a possible interface between the SFA- and GAM-type DPs, by interconnecting the previous DP models of synchronous machines and LCRs proposed in Sections 3.1 and 4.1, respectively. Computer studies validate the accuracy of the proposed DP interface in predicting the ac and dc response of the synchronous machine-rectifier system in all desired operating modes, as well as during large-signal transients. DP Modelling of Synchronous Machine-Rectifier System Components Synchronous machine fed line-commutated converters are commonly found in many ac-dc energy conversion applications, including high-power dc supplies, excitation systems of large electric generators, variable-speed electric drive systems, power systems for aircrafts, vehicles and ships [109]-[110]. In particular, as shown in Figure 1–1 the synchronous machine-rectifier system can be used as the backup generation unit for supporting medium- to long-term interruptions in integrated ac-dc power systems. Without loss of generality, the synchronous machine-rectifier system is depicted by the circuit diagram shown in Figure 5–1. Therein, the synchronous machine is considered in the 133 Figure 5–1. Circuit diagram of the synchronous machine–rectifier system for ac-dc backup generation. CP-VBR formulation (as shown in Section 3.1.1.2) assuming ungrounded stator windings, M damper windings in q-axis, and N damper windings plus a field winding in the d-axis; the three-phase six-pulse diode LCR is connected directly to the synchronous machine; the dc subsystem is represented by filter with fr , fL , and fC , and appropriate dc load. Similar to [109]-[112], the dc subsystem can be represented by an equivalent resistor dcr that consumes the needed amount of real power. This basic configuration (with possible modifications) is used in many industrial applications, and has been validated with different models [109]-[112]. Various operating modes (as discussed in Section 4.1) may exist for such a machine-rectifier system depending on the loading conditions [115]-[116], where high current distortion is observed at light loads (CCM-1), and more distorted voltages at heavy loads when the current is close to sinusoidal waveform (CCM-2). The DCM is typically not observed in this system due to the existence of relatively large machine stator inductances. Synchronous Machine Model in SFA-Type DPs With the use of finely distributed stator/rotor winding and short-/fractional-pitched winding, the synchronous machine can be generally assumed to attain uniformly sinusoidal 134 distributed air gap flux (assuming magnetic linearity), thus inducing sinusoidal emf and can be viewed as a fundamental frequency power source [84]. Moreover, as seen in Figure 5–1, the synchronous machine is directly connected to the diode rectifier in the physical abc-phases. To take advantage of existing DP models of synchronous machines with direct abc-phase stator interface, the CP-SFA model as discussed in Section 3.1.2.2 is considered here, which gives the DP equation of the stator circuit in Figure 5–1 as ''abcsabcsDsDabcsDabcs jdtdEILRILV ][ "" , (5–1) where abcsV and abcsI are the SFA-type DPs of stator voltages and currents, respectively; and DDDD rrrdiag ,,R ; DDDD LLLdiag ,,L , (5–2) where sD rr ; "mdlsD LLL . (5–3) The subtransient voltage sources "abcsE in (5–1) is recalled as TdqabcUqdruabcs ee ]['''',,''KE , (5–4) where abcU qdru,,K is defined in (3–39) and qssqslsmqmdmqMjkqjmqlkqjkqjmqdsmdmqdrlsmqlsmdq irvLLLLLrLiLLLLLLe ''''''12'''''''''''''' )(])([ , (5–5) .)()()(''2''12''''''''''fdlfdmdfdmdlfdfdmdNjkdjmdlkdjkdjmdqsmdmqrqrd vLLLrLLrLiLLe (5–6) In addition, due to the relatively slow electromechanical dynamics, the rotor state equations can be retained in qd coordinates in the time domain as shown in (3–14)–(3–15). 135 Diode LCR Model in GAM-Type DPs Next, it is recalled in Section 4.1 that the diode LCR can be modelled in GAM-type DPs using the analytical or parametric modelling approaches. For comparison, the ADP-2 model [77] in Section 4.1.2.2 is chosen here due to its relatively-higher accuracy, which can be implemented using the ac-dc equations (4–10)–(4–11) and the switch functions (4–14)–(4–15). It is also recalled that the switch functions of fixed-shape waveforms are based on a single mode of operation, i.e., CCM-1. Moreover, since the stator-network inductance has been included in the commutation angle μ [see (4–5)], the ac subsystems of ADP-2 models are formulated as algebraic equations (i.e., reduced-order representation). Alternatively, the PDP modelling approach proposed in Section 4.1.3 is extended to the considered synchronous machine-rectifier system. Specifically, assuming a balanced system and that diode LCR does not contain energy-storing elements, it is recalled that the ac- and dc-side GAM-type DPs can be related through a set of non-linear algebraic functions as Kkevevvevvkjkcskjkbskasijdckkaskas,)()32()32()()ang(01, (5–7) 110)( asdc ii . (5–8) Here, k and 1 are the parametric functions that are determined by the present system operating condition, and the dc current in (5-8) is assumed to be affected by only the fundamental ac currents. Also, the angular differences of phase a GAM-type DPs are recorded as Kkiv askask ),ang()ang()( 1 , (5–9) 136 Finally, the system operating condition, i.e., the input of the parametric functions of kα , 1 , and kφ is defined as the dynamic impedance given by 10asdcivz . (5–10) DC Subsystem Model in the Time Domain It is noted that similar to qd-coordinate rotor variables (in the rotor reference frame), the dynamics of dc subsystems are slow-changing. Considering (2–35) and (2–36), when k = 0, the 0th order GAM-type DP [i.e., )(0tu ], in effect, degrades to the dynamic averaging in the time domain, as is widely used in power electronic converter modelling [107]. Therefore, it is convenient for the dc subsystems to be modelled directly in the time domain. SFA- and GAM-Type DP Interface and System Implementation As discussed in 5.1, the DP modelling of the considered synchronous machine-rectifier system can be summarized into three subsystems: 1) the synchronous machine model represented by SFA-type DPs as a 60Hz VBR source with stator equations (5-1)-(5-4) [and stator-rotor interface (5-5)-(5-6) and rotor dynamics (3–14)–(3–15)]; 2) the GAM-type PDP model of diode LCR comprising the algebraic equations (5-7)-(5-10) and the parametric functions of zkα , z1 , and zkφ ; and 3) the dc subsystem modelled in time-domain signals, which is equivalent to averaging the 0th order GAM-type DPs. 137 To interconnect these subsystems, the synchronous machine-rectifier can be implemented by considering the DP interface as depicted in Figure 5–2. For the considered six-pulse diode LCR, the 5th and 7th harmonics are dominant, and are therefore modelled for this section, i.e., K = {1, 5, 7}. As shown in Figure 5–2, the PDP model can be interfaced to the ac system comprising a number of subsystems for the fundamental frequency and each of the harmonics. In particular, in the 1st order DP subsystem the CP-SFA model of the synchronous machine is viewed as a controlled VBR power source (with adjustment in the magnitude); for 5th and 7th order DPs, it is assumed that a zero voltage source is provided. Finally, the dc subsystem is modelled directly in the time domain [k = 0 in (2–35) and (2–36)]. It is also noted that such a proposed DP interface can be readily augmented to include higher-order harmonics of interest, provided that the dynamics at these frequencies are significant or desired. Figure 5–2. Block diagram of the DP modelling of synchronous machine-rectifier systems. Constructing Parametric Functions For PDP modelling of the considered synchronous machine-rectifier system, the parametric functions of zkα , z1 , and zkφ (K={1, 5, 7}) are established from detailed 138 simulations of the benchmark system as shown in Figure 5–3. Specifically, this is done by running the detailed simulations into multiple steady-states corresponding to a wide range of operating conditions (i.e., the MSS method in Section 4.1.4), and then tabulating the FFT results of the time-domain waveforms of system variables according to (5-7)-(5-10). As seen in Figure 5–3, these parametric functions can become highly nonlinear in the region of heavy loading, where more tightly-spaced data points are typically required. For the considered synchronous machine-rectifier systems [109]-[110], 53 steady-state solution points suffice for modelling accuracy. Figure 5–3. Parametric functions as obtained from detailed simulations: (a)(b)(c) )(1 z , )(1 z , and )(1 z, and (d)(e) )(7,5 zα and )(7,5 zφ . 139 Computer Studies Computer studies are carried out on the benchmark synchronous machine-rectifier system shown in Figure 5–1 with parameters as in [109], [112]. The prior ADP-2 and PDP models have been implemented in MATLAB/Simulink [27]−[28] based on the block diagram shown in Figure 5–2, and the detailed switch-level model is implemented in PLECS [29] for reference. To validate the accuracy of subject models, a case study spanning different operating modes is investigated. The system initially operates in the common CCM-1 mode with a excitation xfdE =19.5 V and an equivalent dc load dcr = 21 ; at t = 0.1 s, a large-signal transient occurs by switching the rectifier to a very heavy load dcr = 0.5, and the simulation is run until t = 0.4s for the system to reach a new steady state in the uncommon CCM-2 mode. Small Step Size Study The system is first solved using small step sizes with the ode23tb solver [27] and the following settings: relative and absolute error tolerances of 10−3, and maximum and minimum step sizes of 100 µs and 0.1 µs, respectively. The system responses of phase a and dc voltages and currents as predicted by the subject models are shown in Figure 5–4. The fragments of Figure 5–4 showing the system responses pre-, post- and during the load change are magnified and shown in Figure 5–5 (in steady states) and Figure 5–6 (in the transient), respectively. From Figure 5–4 to Figure 5–6, the trajectories of asv and asi of DP models are converted from the corresponding DPs kasv and kasi using (2–37). As seen in Figure 5–5 (a)(b), all subject models yield consistent results when the rectifier is in the CCM-1 mode. However, in the new CCM-2 mode, where voltages are highly distorted 140 and currents are close to sinusoidal as shown in Figure 5–5(c)(d), the ADP-2 model fails to reproduce the accurate voltage waveforms [see line b)], while the PDP model results still match the detailed solution very well. Furthermore, this advantageous accuracy of the PDP model can be highlighted in Figure 5–6, which shows the ac and dc responses during the load-change transient. In particular, it is seen in Figure 5–6 that the ADP-2 model renders an inaccurate but quickly converged transient solution (which is due to the inclusion of stator inductance DL into the commutation angle μ [see (4–5) ], thus simplifying the ac system representation). Whereas, the PDP model, by including full-order modelling of harmonic dynamics, can produce large-signal transient results that match well to those predicted by the detailed model [see a) and c)]. Figure 5–4. System response to the load change as predicted by the subject models. 141 Figure 5–5. Magnified view of the ac variables in steady states as depicted in Figure 5–4: (a)(b) asv and asiin the original CCM-1 mode; and (c)(d) asv and asi in the new CCM-2 mode. 142 Figure 5–6. Magnified view of the ac and dc variables during the load-change transient as depicted in Figure 5–4. Flexible Step Size Study Next, the case study in Section 5.4.1 is solved using flexible step sizes [by setting the maximum step size to 0.1s], which yields the step size Δt depicted in Figure 5–7 as a function of time. Therein, very small Δt is used by the detailed model during the entire simulation due to fast switching [see line a)]. In contrast, the Δt curves of DP models start at large values, drop at the load-change instance, and then begin to increase rapidly, which reflects the system operating mode as it starts in steady state, encounters and recovers from a disturbance. The system responses of this case study are similar to Figure 5–4 to Figure 5–6 (except at much larger time steps), and therefore are not included here. 143 Finally, Table 5–1 summarizes the numerical efficiency of the subject models. As can be seen in Table 5–1, the DP models require at least three orders of magnitude fewer time steps than the detailed model and thus significantly outperform the detailed model in terms of CPU time. It is also noted that due to the full-order representation of harmonic dynamics, the PDP model uses more time steps compared to the ADP-2 model, which however does not require much longer CPU time due to its simple structure. Figure 5–7. Step size Δt as chosen by the subject models. Table 5–1. Simulation Efficiency of the Subject Models. Model No. of Time Steps CPU Time, ms Detailed 22763 2234.3 ADP-2 130 65.5 PDP 272 88.7 144 Conclusions and Future Work Conclusions and Contributions Technical challenges resulting from the newly emergent integrated ac-dc power systems necessitate increasing need and reliance on computer simulations to study and understand the system’s dynamic behaviour, and to ensure the system’s operation in an optimal, stable, and secure manner. Several types of simulation programs/tools have been developed throughout the years, each targeted for the investigation of specific classes of transient phenomena. This thesis focuses on the so-called DP-type modelling and simulation techniques, namely the SFA and GAM, which permit general-purpose simulation of both electromagnetic and electromechanical transients. The DP-type modelling approaches represent the power systems using low-pass time-phasor signals, thus offering flexible selection of integration step sizes and a superior combination of numerical accuracy and efficiency. Due to these desirable features, DP-type simulations are particularly suitable for studying the existing and/or new transient phenomena occurring in integrated ac-dc power systems (as highlighted in Figure 1–2). The global objective of this thesis is to increase the numerical efficiency of DP-type simulations for integrated ac-dc power systems, where models of electric machines and power converters are typically the bottleneck in most simulation programs. This thesis addresses these challenges by proposing several new DP component models of electric machines and LCRs with improved numerical properties and more desirable features over the state-of-the-art DP models. Specifically, with respect to the initial objectives of this research, the contributions of this thesis can be summarized as follows: The initial Objective 1 was addressed in Chapter 2, which summarized the fundamentals of DP-type modelling approaches, while pinpointing the limitations of different power system 145 signal representations and the notable numerical features associated with DP-type simulations. The computer studies presented in Section 2.3.2 demonstrated that DP models are numerically accurate in predicting power system transients at different frequencies. Furthermore, it was shown in Section 2.3.3 that the DP models can obtain superior numerical efficiency than the TD models by enabling flexible selection of time step sizes during the simulations. Several important observations on the increased number of state variables, the additional oscillatory modes of DPs, and the signal representation for unbalanced power systems had also been discussed. The initial Objectives 2 and 3 were achieved in Chapter 3, which presented numerically efficient modelling of electric machines based on the SFA. Specifically, a new SFA VBR model of synchronous machines (CP-SFA) was first proposed. This new model achieves a constant-parameter abc-phase stator-network interface, thus avoiding the time-varying terms in the prior variable-parameter SFA model (VP-SFA). The computer studies presented in Section 3.1.3 showed that the CP-SFA model accelerates the simulation speed by at least a factor or two when compared with the VP-SFA model, and reduces the number of required time steps by up to two orders of magnitude compared with classical TD models. This technique was extended to induction machines to achieve a SFA-based VBR model (VBR-SFA-IM) in state-space formulation. The computer studies presented in Section 3.2.3 validated that the VBR-SFA-IM model yields accurate solutions for both fast and slow transients and greatly improved numerical efficiency compared with classical TD models. As highlighted in Table 1–1, this thesis for the first time proposed the SFA-type DP models of synchronous and induction machines that possess direct abc interface, a constant-parameter interfacing circuit, and state-space formulation for SV-based solution. These models represent my contribution and advance the state of the art in this area, as can be seen in Table 1–1 (last two rows). In Chapter 4, the initial Objectives 4 and 5 were achieved by presenting a new highly-efficient parametric DP (PDP) modelling of LCRs based on the GAM. This PDP model for the first time combines the parametric modelling technique with GAM-type DPs in the abc-phase 146 coordinates, and features: 1) simple formulations of ac-dc dynamic subsystems, 2) easy numerical construction of parametric functions, and 3) accurate representations of desired harmonics in all operating modes. Rigorous computer studies presented in Section 4.1.6 demonstrated that the PDP model is capable of accurately capturing the LCR harmonics over a wide range of loading conditions, while providing considerable computational advantages over the detailed model as well as previous ADP models. This PDP methodology was also extended to include the thyristor firing control, which was validated by case studies presented in Section 4.2.3. As highlighted in Table 1–2, this thesis for the first time proposed the GAM-type DP modelling of LCR systems that uses DPs in abc–phase coordinates, covers all operation modes, includes harmonics, and predicts thyristor operation. These models represent my contribution and advance the state of the art in this area, as can be seen in Table 1–2 (last two rows). The initial Objective 6 was achieved in Chapter 5, which for the first time presented the interface between SFA- and GAM-type DPs. By interconnecting the previous DP models proposed in Sections 3.1 and 4.1, the proposed interface obtains a new DP model of the synchronous machine-rectifiers systems, which also contributes to the state-of-the-art DP modelling. Computer studies presented in Section 5.4 validated the accuracy of the proposed DP interface in predicting the ac and dc response in all desired operating modes, as well as during large-signal transients. Future Work To conclude this thesis, several areas for potential extension of the proposed work are discussed below: • Inclusion of Magnetic Saturation into DP Modelling of Rotating Machines Magnetic linearity is assumed for modelling of electric rotating machines in this thesis. However, it is noted that magnetic saturation in rotating machines is a very common phenomenon, 147 which results in changes of effective inductances as the magnetizing flux changes [84]. This requires dedicated consideration and modelling technique to adequately incorporate the effect of magnetic saturation into general-purpose lump-parameter models [133]. A future research task would be to derive SFA-based DP models of rotating machines with the inclusion of magnetic saturations, which will facilitate more accurate prediction of power system transients. • Interfacing DP Models with EMT and/or TS Programs As shown in Chapter 2, the solutions of DP models essentially represent the dynamic envelopes of time-domain signal waveforms, which in steady states can go to constant values as obtained from the conventional phasors. These desirable features infer the potential of the proposed DP models being interfaced with EMT and/or TS programs to provide a straightforward but effective vision into the dynamics of power systems. Several research-based and industry-grade simulation tools [21], [25], [26], [32] have expressed a growing interest in applying DP-type modelling into their simulation for mid- to large-scale ac-dc power systems, and even for real-time simulations [25], [32]. Developing the advanced interfacing techniques between the SFA- and GAM-type DP models and EMT/TS programs is a much-needed research topic that is currently being considered by other students at UBC’s Electric Power and Energy Systems research group. • Combining DP Models and Impedance-Based Stability Criterion for Analyzing Sub-/Super-Synchronous Resonances in Power Electronic Systems The oscillatory interactions between the power electronic component/subsystem and the grid have emerged as a notable practical problem that has gained much attention in recent years [14]-[16]. Unlike the large-scale stiff ac grid, power-electronic-based systems may provoke electrical resonances ranging from the sub- to super-synchronous frequencies, as depicted in Figure 1–2. To analyze and mitigate such resonance problems, the impedance-based stability theory has been proposed [14] based on the Nyquist stability criterion. As shown in Chapter 2, compared with conventional phasor models that neglect dynamics, the DP models yield accuracy 148 simulation results at various frequencies (due to their full-order modelling), which thus permit the investigation of sub-/super-synchronous resonances from the phasor perspective. Therefore, a future research task would be to combine DP models with impedance-based stability criterion for analyzing the sub-/super-synchronous resonances in power electronic based systems. 149 References [1] M. Shahidehpour, "Our aging power systems: infrastructure and life extension issues," IEEE Power Energy Mag., vol. 4, no. 3, pp. 22-76, May. 2006. [2] J. Arrillaga, High Voltage Direct Current Transmission, 2nd ed. London, U.K.: Inst. Electr. Eng., 1998. [3] Y. Song and A. T. Johns, Flexible AC Transmission Systems (FACTS). IET, 1999. [4] M. Barnes, D. Van Hertem, S. P. Teeuwsen, and M. Callavik, "HVDC systems in smart grids," in Proc. IEEE, to appear. [5] J. Cao and J. Y. Cai, “HVDC in China,” presented at EPRI 2013 HVDC & FACTS Conference, USA, 2013. [6] R. E. P. N. for the 21st Century, “Renewables 2017 Global Status Report,” Tech. Rep. 2017, REN21. [7] “Electric Power Monthly with Data for March 2017,” U.S. Energy Information Administration (EIA), U.S. Department of Energy, Washington, DC, Jun. 2017. [Online]. Available: https://www.eia.gov/electricity/monthly/ . (Consulted on Jun. 28, 2017). [8] U.S. Department of Energy (DOE). (2014) Smart Grid System Report. [Online]. Available: https://www.smartgrid.gov/files/2014-Smart-Grid-System-Report.pdf (Consulted on Jun. 28, 2017). [9] N. Hatziargyriou, H. Asano, R. Iravani, and C.Marnay, “Microgrids,” IEEE Power Energy Mag., vol. 5, no. 4, pp. 78–94, Jul./Aug. 2007. [10] D. E. Olivares, et al., "Trends in microgrid control," IEEE Tran. Smart Grid, vol. 5, no. 4, pp. 1905-1919, Jul. 2014. [11] Univ. of British Columbia (UBC). (2017) UBC’s Living Lab Report [Online]. Available: https://sustain.ubc.ca/our-commitment/campus-living-lab (Consulted on Jun. 28, 2017). [12] “Challenges to the Integration of Renewable Resources at High System Penetration,” California Institute for Energy and Environment, Berkeley, CA, May 2014. [Online]. 150 Available: http://www.energy.ca.gov/2014publications/CEC-500-2014-042/CEC-500-2014-042.pdf. (Consulted on Jun. 28, 2017). [13] S. Santoso, Fundamentals of Electric Power Quality. CreateSpace Independent Publishing Platform, 2009. [14] J. Sun, "Impedance-based stability criterion for grid-connected inverters," IEEE Trans. Power Electron., vol. 26, no. 11, pp. 3075-3078, Nov. 2011. [15] M. Amin and M. Molinas, "Understanding the origin of oscillatory phenomena observed between wind farms and HVdc systems," IEEE Trans. Emerg. Sel. Topics Power Electron. , vol. 5, no. 1, pp. 378-392, March 2017. [16] X. Wang, F. Blaabjerg and W. Wu, "Modelling and analysis of harmonic stability in an AC power-electronics-based power system," IEEE Trans. Power Electron., vol. 29, no. 12, pp. 6421-6432, Dec. 2014. [17] J. Arrillaga and N. R. Watson, Computer Modelling of Electrical Power Systems, 2nd ed. Chichester, UK: John Wiley & Sons, Ltd, 2001. [18] P. W. Sauer and M. A. Pai, Power System Dynamics and Stability, Prentice-Hall, 1998. [19] P. Kundur, Power System Stability and Control. New York: McGraw-Hill, 1994. [20] H. W. Dommel, EMTP Theory Book. MicroTran Power System Analysis Corporation, Vancouver, Canada, May 1992. [21] “PSCAD/EMTDC X4 On-Line Help,” Manitoba HVDC Research Centre and RTDS Technologies Inc., Winnipeg, MB, Canada, 2016. [Online]. Available: http://www.pscad.com. (Consulted on Jun. 28, 2017). [22] “Electromagnetic Transient Program,” EMTP-RV, CEA Technologies Inc., 2016. [Online]. Available: http://www.emtp.com. (Consulted on Jun. 28, 2017). [23] “MicroTran Reference Manual,” MicroTran Power System Analysis Corp., Vancouver, BC, Canada, 1997. [Online]. Available: http://www.microtran.com. (Consulted on Jun. 28, 2017). 151 [24] “Alternative Transients Programs,” ATP-EMTP, ATP User Group, 2007. [Online]. Available: http://www.emtp.org. (Consulted on Jun. 28, 2017). [25] RTDS Simulator Software, RTDS Technologies Inc., Winnipeg, Manitoba, Canada, 2016. [Online]. Available: http://www.rtds.com. (Consulted on Jun. 28, 2017). [26] CloudPSS Simulation Platform, Energy Internet Research Institute, Tsinghua University, Beijing, China, 2017. [Online]. Available: http://www.cloudpss.net . (Consulted on Jun. 28, 2017). [27] MATLAB 7 (R2016b), Natick, MA: The MathWorks Inc., Natick, MA, USA, 2016. [28] Simulink – User’s Manual, The MathWorks, Inc., Natick, MA, USA, 2016. [29] PLECS Blockset, Plexim GmbH, Zurich, Switzerland, 2016. [Online]. Available: www.plexim.com. (Consulted on Jun. 28, 2017). [30] Simscape Power Systems, Hydro-Quebec and the MathWorks, Inc., Natick, MA, 2016. [31] “Automated State Model Generator (ASMG) Reference Manual,” P. C. Krause and Associates, Inc., West Lafayette, IN, 2002. [32] RT-LAB, OPAL-RT Technologies, Montréal, Québec, Canada, 2016. [Online]. Available: http://www.opal-rt.com. (Consulted on Jun. 28, 2017). [33] O. Wasynczuk and S. D. Sudhoff, “Automated state model generation algorithm for power circuits and systems,” IEEE Trans. Power Syst., vol. 11, no. 4, pp. 1951–1956, Nov. 1996. [34] J. V. Jatskevich, “A state selection algorithm for the automated state model generator,” Ph.D. dissertation, Purdue University, West Lafayette, IN, USA, 1999. [35] D. Jakominich, R. Krebs, D. Retzmann, and A. Kumar, "Real time digital power system simulator design considerations and relay performance evaluation," IEEE Trans. Power Del., vol. 14, no.3, pp. 773-781, Jul. 1999. [36] TSAT, Powertech Labs Inc., Vancouver, BC, Canada, 2017. [Online]. Available: http://www.powertechlabs.com/areas-of-focus/smart-utility/dsatools-software/transient-security-assessment-tool/ (Consulted on Jun. 28, 2017). 152 [37] PSS/E, Siemens Power Technologies International, Schenectady, NY, USA, 2017. [Online]. Available : http://w3.siemens.com/smartgrid/global/en/products-systems-solutions/software-solutions/planning-data-management-software/planning-simulation/pages/pss-e.aspx. (Consulted on Jun. 28, 2017). [38] Power System Analysis Software Package (PSASP) User Manual. Electric Power Research Institute of China (CEPRI), Beijing, China, 2002. [39] DIgSILENT PowerFactory, DIgSILENT GmbH, Gomaringen, Germany, 2017. [Online]. Available: http://www.digsilent.de/index.php/products-powerfactory.html (Consulted on Jun. 28, 2017). [40] V. Jalili-Marandi, V. Dinavahi, K. Strunz, J. A. Martinez, and A. Ramirez, "Interfacing techniques for transient stability and electromagnetic transient programs IEEE Task Force on interfacing techniques for simulation tools," IEEE Trans. Power Del., vol. 24, no.4, pp.2385-2395, Oct. 2009. [41] J. R. Marti and L. R. Linares, “Real-time EMTP-based transients simulation,” IEEE Trans. Power Syst., vol. 9, no. 3, pp. 1309–1317, Aug. 1994. [42] J. R. Marti, L. R. Linares, J. Calvino, H. W. Dommel, and J. Lin. "OVNI: an object approach to real-time power system simulators," in Proc. IEEE Int. Conf. Power Syst. Technology (Powercon), Beijing, China, Aug. 1998, pp.977-981. [43] A. S. Morched and V. Brandwajn, “Transmission network equivalents for electromagnetic transients studies,” IEEE Trans. Power App. Syst., vol. PAS-102, no. 9, pp. 2984–2994, Sep. 1983. [44] D. M. Falcao, E. Kaszkurewicz, and H. L. S. Almeida, “Application of parallel processing techniques to the simulation of power system electromagnetic transients,” IEEE Trans. Power Syst., vol. 8, no. 1, pp. 90–96, Feb. 1993. [45] H.T. Su, K. W. Chan, and L. A. Snider, "Evaluation study for the integration of electromagnetic transients simulator and transient stability simulator," Elec. Power Syst. Res., vol. 75, no.1, pp. 67-78, Jul. 2005. 153 [46] S. Abhyankar, “Development of an implicitly coupled electromechanical and electromagnetic transients simulator for power systems,” Ph.D. dissertation, Dept. Elec. Comput. Eng., Illinois Inst. Technol., Chicago, USA, 2011. [47] M. Ilić, and J. Zaborszky, Dynamics and Control of Large Electric Power Systems. New York: Wiley, 2000. [48] C. E. Shannon, “Communication in the presence of noise,” Proc. IRE., vol. 37, no.1, pp. 10–21, Jan. 1949. [49] V. Venkatasubramanian, "Tools for dynamic analysis of the general large power system using time-varying phasors," Int. J. Elec. Power Energy Syst., vol. 16, no. 6, pp. 365-376, Dec. 1994. [50] S. Henschel, “Analysis of electromagnetic and electromechanical power system transients with dynamic phasors,” Ph.D. dissertation, Dept. Elec. Comput. Eng., Univ. British Columbia, Vancouver, BC, Canada, 1999. [51] P. Zhang, J. R. Marti, and H. Dommel, “Shifted-frequency analysis for EMTP simulation of power-system dynamics,” IEEE Trans. Circuits Syst. I: Reg. Papers, vol. 57, no. 9, pp. 2564–2574, Sep. 2010. [52] J. R. Martí, H. W. Dommel, B. D. Bonatto, F. R. Barreto, "Shifted frequency analysis (SFA) concepts for EMTP modelling and simulation of power system dynamics," presented at the 18th Power Systems Computation Conference, Wroclaw, Poland, Aug. 18, 2014. [53] P. Zhang, J. R. Marti, and H. W. Dommel, "Synchronous machine modelling based on shifted frequency analysis," IEEE Trans. Power Syst., vol. 22, no.3, pp. 1139-1147, Aug. 2007. [54] P. Zhang, J. R. Marti, and H. W. Dommel, "Induction machine modelling based on shifted frequency analysis," IEEE Trans. Power Syst., vol. 24, no.1, pp. 157-164, Feb. 2009. [55] M.A. Elizondo, F. K. Tuffner, K. P. Schneider, “Simulation of inrush dynamics for unbalanced distribution systems using dynamic-phasor models,” IEEE Trans. Power Syst., vol. 32, pp. 633-642, 2017. 154 [56] K. Strunz, R. Shintaku, and F. Gao, “Frequency-adaptive network modelling for integrative simulation of natural and envelope waveforms in power systems and circuits,” IEEE Trans. Circuits Syst., vol. 53, no.12, pp. 2788–2803, Dec. 2006. [57] F. Gao and K. Strunz, “Multi-scale simulation of multi-machine power systems,” in Proc. 16th Power System Computation Conf. (PSCC), Glasgow, U.K., Jul. 2008. [58] H. Ye, F. Gao, K. Strunz, and Y. Xia, "Multi-scale modelling and simulation of synchronous machine in phase-domain," in Proc. 23rd IEEE PES Innovative Smart Grid Technologies (ISGT Europe), Berlin, Germany, Oct. 2012. [59] H. Ye, Y. Tang and Y. Xia, "DQ-domain modelling for multi-scale transients in a synchronous machine," in Proc. 5th Int. Conf. Electric Utility Deregulation and Restructuring and Power Technologies (DRPT), Changsha, China, 2015, pp. 285-289. [60] H. Ye, and K. Strunz, "Multi-scale and frequency-dependent modelling of electric power transmission lines," IEEE Trans. Power Del., to appear. [61] H. Ye, B. Yue, and X. Li, “Modelling and simulation of multi-scale transients for PMSG-based wind power systems,” Wind Energy, vol.20, pp. 1349-1364, 2017. [62] S. R. Sanders, J. M. Noworolski, X. Z. Liu, and G. C. Verghese, "Generalized averaging method for power conversion circuits," IEEE Trans. Power Electron., vol. 6, pp. 251–259, Apr. 1991 [63] T. H. Demiray, "Simulation of power system dynamics using dynamic phasors models," Ph.D. dissertation, Swiss Fed. Inst. Technol., Zurich, Switzerland, 2008. [64] T. Yang, "Development of dynamic phasors for the modelling of aircraft electrical power systems," Ph.D. dissertation, Univ. Nottingham, Nottingham, UK, 2013. [65] T. Yang, S. Bozhko, J. M. Le-Peuvedic, G. Asher and C. I. Hill, "Dynamic phasor modelling of multi-generator variable frequency electrical power systems," IEEE Trans. Power Syst., vol. 31, no. 1, pp. 563-571, Jan. 2016. [66] A. M. Stankovic, B. Lesieutre, and T. Aydin, “Modelling and analysis of single-phase induction machines with dynamic phasors,” IEEE Trans. Power Syst., vol. 14, no. 1, pp. 9-14, Feb. 1999. 155 [67] A. M. Stankovic, and T. Aydin, "Analysis of asymmetrical faults in power systems using dynamic phasors," IEEE Trans. Power Syst., vol. 15, no. 3, pp. 1062-1068, Aug. 2000. [68] A. M. Stankovic, S. R. Sanders, and T. Aydin, "Dynamic phasors in modelling and analysis of unbalanced polyphase AC machines," IEEE Trans. Energy Convers., vol. 17, pp. 107-113, Mar. 2002. [69] P. Mattavelli, G. C. Verghese, and A. M. Stankovic, "Phasor dynamics of thyristor-controlled series capacitor systems," IEEE Trans. Power Syst., vol. 12, no. 3, pp. 1259-1267, Aug. 1997. [70] P. Mattavelli, A. M. Stankovic, and G. C. Verghese, "SSR analysis with dynamic phasor model of thyristor-controlled series capacitor," IEEE Trans. Power Syst., vol. 14, pp. 200-208, Feb. 1999. [71] P. C. Stefanov and A. M. Stankovic, "Modelling of UPFC operation under unbalanced conditions with dynamic phasors," IEEE Trans. Power Syst., vol. 17, pp. 395-403, May. 2002. [72] Z Zhijun, D.Z. Fang, K.W. Chan, and S.Q. Yuan, “Hybrid simulation of power systems with SVC dynamic phasor model,” Elec. Power Syst. Res., vol. 31, pp. 175-180, Jan. 2009. [73] Q. Qi, S. Chen, Y. Ni, and F. F. Wu, “Application of the dynamic phasors in modelling and simulation of HVDC,” in Int. Conf. Advances in Power System Control, Operation and Management, Hong Kong, Nov. 2003. [74] H. Zhu, Z. Cai, H. Liu, Q. Qi, and Y. Ni, "Hybrid-model transient stability simulation using dynamic phasors based HVDC system model," Elec. Power Syst. Res., vol. 76, pp. 582-591, Apr. 2006. [75] M. Kulasza, “Generalized dynamic phasor-based simulation for power systems,” M.Sc. thesis, Univ. Manitoba, Winnipeg, MB, Canada, 2014. [76] M. Daryabak, S. Filizadeh, J. Jatskevich, A. Davoudi, M. Saeedifard, V.K. Sood, J. A. Martinez, D. Aliprantis, J. Cano, and A. Mehrizi-Sani, “Modelling of LCC-HVDC systems using dynamic phasors,” IEEE Trans. Power Del., vol. 29, no.4, pp. 1989-1998, Aug. 2014. 156 [77] C. Liu, A. Bose, and P. Tian, “Modelling and analysis of HVDC converter by three-phase dynamic phasor,” IEEE Trans. Power Del., vol. 29, no.1, pp. 3-12, Feb. 2014. [78] T. Yang, S. Bozhko, and G. Asher, “Fast functional modelling of diode-bridge rectifier using dynamic phasors,” IET Power Electron., vol. 8, no.6, pp. 947-956, Jun. 2015. [79] A. Emadi, “Modelling and analysis of multiconverter DC power electronic systems using the generalized state-space averaging method,” IEEE Trans. Ind. Electron., vol. 51, pp. 661-668, Jun. 2004. [80] A. Emadi, “Modelling of power electronics loads in AC distribution systems using the generalized state-space averaging method,” IEEE Trans. Ind. Electron., vol. 51, no.3, pp. 992-1000, Oct. 2004. [81] Z. Miao, L. Piyasinghe, J. Khazaei and L. Fan, "Dynamic phasor-based modelling of unbalanced radial distribution systems," IEEE Trans. Power Syst., vol. 30, no. 6, pp. 3102-3109, Nov. 2015. [82] M. C. Chudasama and A. M. Kulkarni, "Dynamic phasor analysis of SSR mitigation schemes based on passive phase imbalance," IEEE Trans. Power Syst., vol. 26, no. 3, pp. 1668-1676, Aug. 2011. [83] S. Chandrasekar and R. Gokaraju, "Dynamic phasor modelling of type 3 DFIG wind generators (including SSCI phenomenon) for short-circuit calculations," IEEE Trans. Power Del., vol. 30, no. 2, pp. 887-897, April 2015. [84] P. C. Krause, O. Wasynczuk, and S. D. Sudhoff, Analysis of Electric Machinery and Drive Systems, 2nd ed. Piscataway, NJ, USA: IEEE Press, 2002. [85] B. Wu, Y. Lang, N. Zargari, S. Kouro, Power Conversion and Control of Wind Energy Systems, Wiley-IEEE Press, Piscataway, NJ, 2011 [86] L. Wang, J. Jatskevich, V. Dinavahi, H. W. Dommel, J. A. Martinez, K. Strunz, M. Rioual, G. W. Chang, and R. Iravani, "Methods of interfacing rotating machine models in transient simulation programs," IEEE Trans. Power Del., vol. 25, no.2, pp. 891-903, Apr. 2010. [87] R. H. Park, "Two-reaction theory of synchronous machines generalized method of analysis: part I," Trans. Amer. Instit. Electr. Eng., vol. 48, no. 3, pp. 716–727, Jul. 1929. 157 [88] S. D. Pekarek, O. Wasynczuk, and H. J. Hegner, "An efficient and accurate model for the simulation and analysis of synchronous machine/converter systems," IEEE Trans. Energy Convers., vol. 13, no. 1, pp. 42–48, Mar. 1998. [89] S. D. Pekarek and E. A. Walters, "An accurate method of neglecting dynamic saliency of synchronous machines in power electronic based systems," IEEE Trans. Energy Convers., vol. 14, no.4, pp. 1177-1183, Dec. 1999. [90] M. Chapariha, L. Wang, J. Jatskevich, H. W. Dommel, and S. D. Pekarek, “Constant-parameter R-L branch equivalent circuit for interfacing AC machine models in state-variable-based simulation packages,” IEEE Trans. Energy Convers., vol. 27, no. 3, pp. 634-645, Sep. 2012. [91] F. Therrien, M. Chapariha, and J. Jatskevich, “Pole selection procedure for explicit constant-parameter synchronous machine models,” IEEE Trans. Energy Convers., vol. 29, no.3, pp. 790-792, Sep. 2014. [92] M. Chapariha, F. Therrien, J. Jatskevich, and H. W. Dommel, "Explicit formulations for constant-parameter voltage-behind-reactance interfacing of synchronous machine models," IEEE Trans. Energy Convers., vol. 28, no.4, pp. 1053-1063, Dec. 2013. [93] M. Chapariha, F. Therrien, J. Jatskevich, and H. W. Dommel, "Constant-parameter circuit-based models of synchronous machines," IEEE Trans. Energy Convers., vol.30, no.2, pp.441-452, Jan. 2015. [94] L. Wang, J. Jatskevich, and S. D. Pekarek, “Modelling of induction machines using a voltage-behind-reactance formulation,” IEEE Trans. Energy Convers., vol. 23, no. 2, pp. 382–392, Jun. 2008. [95] L. Wang and J. Jatskevich, “A voltage-behind-reactance synchronous machine model for the EMTP-type solution,” IEEE Trans. Power Syst., vol. 21, no. 4, pp. 1539–1549, Nov. 2006. [96] L. Wang, J. Jatskevich, C. Wang, and P. Li, “A voltage-behind-reactance induction machine model for the EMTP-type solution,” IEEE Trans. Power Syst., vol. 23, no. 3, pp. 1226–1238, Aug. 2008. 158 [97] U. Karaagac, J. Mahseredjian, I. Kocar, and O. Saad, “An efficient voltage-behind-reactance formulation-based synchronous machine model for electromagnetic transients,” IEEE Trans. Power Del., vol. 28, no. 3, pp. 1788–1795, Jul. 2013. [98] N. Mohan, and T. M. Undeland, Power Electronics: Converters, Applications, and Design. New York, NY, USA: Wiley, 2007. [99] G. Seguir, Power Electronic Converters: AC-DC Conversion. New York, NY, USA: McGraw-Hill, 1986. [100] A. M. De Broe, S. Drouilhet, and V. Gevorgian, “A peak power tracker for small wind turbines in battery charging applications,” IEEE Trans. Energy Convers., vol. 14, no. 4, pp. 1630–1635, Dec. 1999. [101] A. Emadi, M. Ehsani, and J. M. Miller, Vehicular Electric Power Systems: Land, Sea, Air, and Space Vehicles. New York, NY, USA: Marcel Dekker, 2004. [102] J. Arrillaga, Power System Harmonic Analysis, New York, NY, USA: Wiley, 1997. [103] A.M. Gole, et al., “Optimization-enabled electromagnetic transient simulation,” IEEE Trans. Power Del., vol. 20, no. 1, pp. 512–518, Jan. 2005. [104] H. Visser and P. van den Bosch, “Modelling of periodically switched networks,” Proc. IEEE PESC, pp. 67–73, 1991. [105] P. W. Lehn, "Exact modelling of the voltage source converter," IEEE Trans. Power Del., vol. 17, no. 1, pp. 217-222, Jan. 2002. [106] K. L. Lian, B. K. Perkins and P. W. Lehn, "Harmonic Analysis of a Three-Phase Diode Bridge Rectifier Based on Sampled-Data Model," IEEE Trans. Power Del., vol. 23, no. 2, pp. 1088-1096, Apr. 2008. [107] D. Maksimovic, A. M. Stankovic, V. J. Thottuvelil and G. C. Verghese, "Modelling and simulation of power electronic converters," Proc. IEEE, vol. 89, no. 6, pp. 898-912, Jun. 2001. [108] S. Baacha, L. Munteanu, and A. Luliana Bratcu, Power Electronics Converters Modelling and Control: With Case Studies. London: Springer-Verlag, 2014. 159 [109] S. D. Sudhoff and O. Wasynczuk, “Analysis and average-value modelling of line-commutated converter-synchronous machine systems,” IEEE Trans. Energy Convers., vol. 8, no. 1, pp. 92–99, Mar. 1993. [110] S. D. Sudhoff, K. A. Corzine, H. J. Hegner, and D.E. Delisle, “Transient and dynamic average-value modelling of synchronous machine fed load-commutated converters,” IEEE Trans. Energy Convers., vol. 11, no.3, pp. 508-514, Sep. 1996. [111] I. Jadric, D. Borojevic, and M. Jadric, “Modelling and control of a synchronous generator with an active dc load,” IEEE Trans. Power Electron., vol. 15, no. 2, pp. 303–311, Mar. 2000. [112] J. Jatskevich, S. D. Pekarek, and A. Davoudi, “Parametric average-value model of a synchronous machine-rectifier system,” IEEE Trans. Energy Convers., vol. 21, no.1, pp. 9-18, Mar. 2006. [113] J. Jatskevich, S. D. Pekarek, and A. Davoudi, “Fast procedure for constructing an accurate dynamic average-value model of synchronous machine-rectifier system,” IEEE Trans. Energy Convers., vol. 21, no.2, pp. 435-441, Jun. 2006. [114] S. Chiniforoosh, J. Jatskevich, A. Yazdani, V. Sood, V. Dinavahi, J. A. Martinez, and A. Remirez, “Definitions and applications of dynamic average models for analysis of power systems,” IEEE Trans. Power Del., vol. 25, no. 4, pp. 2655–2669, Oct. 2010. [115] S. Chiniforoosh, H. Atighechi, A. Davoudi, J. Jatskevich, A. Yazdani, S. Filizadeh, S. Saeedifard, J. A. Martinez, V. Sood, K. Strunz, J. Mahseredjian, and V. Dinavahi, “Dynamic average modelling of front-end diode rectifier loads considering discontinuous conduction mode and unbalanced operation,” IEEE Trans. Power Del., vol. 27, no.1, pp. 421-429, Jan. 2012. [116] S. Chiniforoosh, H. Atighechi, A. Davoudi, J. Jatskevich, J. A. Martinez, M. Saeedifard, D. C. Aliprantis, and V. K. Sood, “Steady-state and dynamic performance of front-end diode rectifier loads as predicted by dynamic average-value models,” IEEE Trans. Power Del., vol. 28, no.3, pp. 1533-1541, Jul. 2013. 160 [117] H. Atighechi, S. Chiniforoosh, K. Tabarraee and J. Jatskevich, "Average-value modelling of synchronous-machine-fed thyristor-controlled-rectifier systems," IEEE Trans. Energy Convers., vol. 30, no. 2, pp. 487-497, June 2015. [118] H. Atighechi, S. Chiniforoosh, S. Ebrahimi, and J. Jatskevich, “Using multiple reference frame theory for considering harmonics in average-value modelling of diode rectifiers,” IEEE Trans. Energy Convers., vol. 31, no. 3, pp. 872 – 881, Sept. 2016. [119] L. Cohen, Time-Frequency Analysis, New York: Prentice-Hall, 1995. [120] C. E. Shannon, "Communication in the presence of noise,” Proc. IRE., vol. 37, no.1, pp. 10–21, Jan. 1949. [121] W. Gautchi, Numerical Analysis: An Introduction. Boston, MA: Birkhauser, 1997. [122] L. Gyugyi, and B. R. Pelly, Static Power Frequency Changers. – Theory, Performance, and Application. New York, USA: Wiley, 1976. [123] L. Hu, and R. Yacamini, “Harmonic transfer through converters and HVDC links,” IEEE Trans. Power Electron., vol. 7, no.3, pp. 514-525, Jul. 1992. [124] F. Yahyaie, and P. W. Lehn, “On dynamic evaluation of harmonics using generalized averaging techniques,” IEEE Trans. Power Syst., vol.30, no.5, pp. 2216-2224, Sep. 2015. [125] U. Vargas, and A. Ramirez, “Reformulating extended harmonic domain models for accurate representation of harmonics dynamics,” IEEE Trans. Power Del., vol.31, no.6, pp. 2562-2564, May. 2016. [126] J. Rico, M. Madrigal, and E. Acha, “Dynamic harmonic evolution using the extended harmonic domain,” IEEE Trans. Power Del. vol. 18, no. 2, pp. 587-594, Apr. 2003. [127] P. F. Ribeiro, Time-Varying Waveform Distortions in Power Systems. Hoboken, NJ: Wiley-IEEE Press, 2009. [128] T. Bi, H. Liu, Q. Feng, C. Qian and Y. Liu, "Dynamic phasor model-based synchrophasor estimation algorithm for M-Class PMU," IEEE Trans. Power Del., vol. 30, no. 3, pp. 1162-1171, Jun. 2015. 161 [129] R. A. Horn, and C. R. Johnson, Matrix Analysis, Cambridge, England: Cambridge University Press, 1990. [130] M. Bauta, and M. Grotzbach, “Noncharacteristic line harmonics of AC/DC converters with high DC current ripple,” IEEE Trans. Power Del., vol. 15, no. 3, pp. 1060-1066, Jul. 2000. [131] S. M. Ebrahimi, N. Amiri, H. Atighechi, J. Jatskevich, and L. Wang, “Performance verification of parametric average-value model of line-commutated rectifiers under unbalanced conditions,” in Proc. IEEE 16th Workshop on Control and Modeling for Power Electronics (COMPEL), Vancouver, Canada, 2015, pp. 1–6. [132] Z. Bing, K. Karimi, J. Sun, "Input impedance modeling and analysis of line-commutated rectifiers," IEEE Trans. Power Electron., vol.24, no.10, pp.2338-2346, Oct. 2009. [133] F. Therrien, “Numerically efficient modelling of saturable ac machines for power system electromagnetic transient simulation programs,” Ph.D. dissertation, Dept. Elec. Comput. Eng., Univ. British Columbia, Vancouver, BC, Canada, 2015. 162 Appendix A: Transmission Line Parameters for the Case Study of Section 2.3.3.3 Transmission line parameters: Positive Sequence Zero Sequence Resistances km/ 0.01273 = posr km/ 0.3864 = zeror Inductances km/mH 0.9937 =posl km/mH .12644 =zerol Capacitances km/nF 12.74 =posc km/nF 7.751 =zeroc For the 50 km transmission line ( km50linel ), the parameters of the lumped-parameter π section model (short line model) are calculated as (excluding hyperbolic corrections) [20] lineposposlineposposlineposposlcCllLlrR , linezerozerolinezerozerolinezerozerolcCllLlrR. (A-1) The three-phase coupled π section RLC matrices are formulated as [20] stlmtlmtlmtlstlmtlmtlmtlstltlRRRRRRRRR_________R , (A-2) stlmtlmtlmtlstlmtlmtlmtlstltlLLLLLLLLL_________L , (A-3) stlmtlmtlmtlstlmtlmtlmtlstltlCCCCCCCCC_________C . (A-4) where the diagonal and off-diagonal elements are calculated as )2(31_ zeroposstl RRR , )(31_ poszeromtl RRR (A-5) 163 )2(31_ zeroposstl LLL , )(31_ poszeromtl LLL (A-6) )2(31_ zeroposstl CCC , )(31_ poszeromtl CCC . (A-7) 164 Appendix B: Parameters for the Case Study of Section 3.1.3 Synchronous machine parameters [84]: 835 MVA, 26 kV, 0.85 power factor, 2 poles, 3600 r/min, J = 0.0658×106 J·s2, rs = 0.00243Ω, Xls = 0.1538Ω, Xmq = 1.3032Ω, rkq1 = 0.00144Ω, Xlkq1 = 0.6578Ω, rkq2 = 0.0068Ω, Xlkq2 = 0.07602Ω, Xmd = 1.3032Ω, rfd = 0.0075Ω, Xlfd = 0.1145Ω, rkd = 0.0108Ω, Xlkd = 0.06577Ω. 165 Appendix C: Parameters for the Case Study of Section 3.2.3 Induction machine parameters [84]: 500 HP, 2300 V, 4 poles, 1773 rpm, J = 11.06 kg·m2, Tb = 1980 N·m, rs = 0.262Ω, rr = 0.187Ω, Xls = 1.206Ω, Xlr = 1.206Ω, XM = 54.02Ω. 166 Appendix D: Parameters for the Case Study of Sections 4.1.6 and 4.2.3 Testing System Parameters [118]: Vline = 480 V, fs = 60 Hz, rs =0.2Ω, Ls = 10 mH, rdc = 0.5 Ω, Ldc = 1.33 mH, Cdc = 500 µF.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerically efficient modelling and simulation of integrated...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerically efficient modelling and simulation of integrated ac-dc power systems using dynamic phasor… Huang, Yingwei 2017
pdf
Page Metadata
Item Metadata
Title | Numerically efficient modelling and simulation of integrated ac-dc power systems using dynamic phasor type solutions |
Creator |
Huang, Yingwei |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | Modern power systems, from continent-spanning networks down to isolated microgrids, are experiencing unprecedented technological changes with broader use of direct current (dc) in addition to traditional alternating current (ac). Such integrated ac-dc power systems present notable challenges in all aspects of design, analysis, control, and operation, where extensive computer simulations play the essential and enabling role. Due to the use of diverse types of signal representation and component formulation, state-of-the-art power system simulation tools are limited to their distinct time scales of transient phenomena. This thesis considers the dynamic phasor (DP) type modelling approaches, where two types of DP theories, namely the shifted-frequency analysis (SFA) and the generalized averaging method (GAM), are considered. In DP-type simulations, power systems are modelled using analogous low-pass time-phasor signals, thereby offering flexible selection of time-step sizes and superior combination of numerical accuracy and efficiency. The ultimate goal of this research is to increase the numerical efficiency of DP-type simulations for the integrated ac-dc power systems. This is achieved by proposing several new DP component models with desirable features and improved numerical properties. First, the constant-parameter SFA model of synchronous machines is proposed to avoid numerically-costly recalculations of the time-varying stator-network matrices. This model is then extended to induction machines for modelling in state-variable based (SV-based) simulation tools. Next, we propose a new, highly efficient model of line-commutated rectifiers using a parametric DP formulation, which is demonstrated as valid for various system operating conditions. Moreover, the effect of ac side harmonics is incorporated to improve modelling fidelity. Finally, the interface between SFA- and GAM-type DPs is achieved to interconnect the proposed DP models. Rigorous case studies demonstrate the superior numerical efficiency of the proposed models, and their advantageous accuracy in capturing the desired phenomena of ac-dc power systems. It is envisioned that the proposed models will become highly useful to many researchers and engineers worldwide, and facilitate the development of next-generation power system simulation tools. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-01-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0363001 |
URI | http://hdl.handle.net/2429/64322 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_ february_huang_yingwei.pdf [ 7.73MB ]
- Metadata
- JSON: 24-1.0363001.json
- JSON-LD: 24-1.0363001-ld.json
- RDF/XML (Pretty): 24-1.0363001-rdf.xml
- RDF/JSON: 24-1.0363001-rdf.json
- Turtle: 24-1.0363001-turtle.txt
- N-Triples: 24-1.0363001-rdf-ntriples.txt
- Original Record: 24-1.0363001-source.json
- Full Text
- 24-1.0363001-fulltext.txt
- Citation
- 24-1.0363001.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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0363001/manifest