MINIMUM!IHHIIWNWHHI{INN“HIHIHHUNHHI 570 0853 r This is to certify that the dissertation entitled BIFURCATION SUBSYSTEM METHOD AND ITS APPLICATION TO DIAGNOSIS OF POWER SYSTEM BIFURCATIONS PRODUCED BY DISCONTINUITIES presented by Khadija Ben-Kilani has been accepted towards fulfillment of the requirements for Ph.D degree in Electrical Eng GEM AMM ta Major professor Date é/EO/j7 MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771 ‘w~‘ —" A.—.. “A v, 5 i ; LIBRARY Michigan State University PLACE It RETURN 30X to remove thie checkout from your record. TO AVOID FINES return on or betore date due. DATE DUE DATE DUE DATE DUE MSU ie An Affirmative ActioNEquel Opportunity Inetituion mane-m BIR'RC DIAG BIFURCATION SUBSYSTEM METHOD AND ITS APPLICATION TO DIAGNOSIS OF POWER SYSTEM BIFURCATIONS PRODUCED BY DISCONTINUITIES by Khadija Ben-Kilani A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering 1 997 81R M l 7"?" 59 mt“ bi. was ’V'.“ ‘ ilk ABSTRACT BIFURCATION SUBSYSTEM METHOD AND ITS APPLICATION TO DIAGNOSIS OF POWER SYSTEM BIFURCATIONS PRODUCED BY DISCONTINUITIES by Khadija Ben-Kilani Recurrent problems in diagnosing the cause and location of a stability problem in a power system are that of (a) the size and complexity of the model, (b) the various kinds and classes of bifurcations, and (c) the effects of discontinuous hard limit or equipment outage-induced transitions. Problem diagnosis in power systems is in need of stability assessment methods that (i) take into account both the existence of discontinuities and their effect, and (ii) identifies which physical elements or subsystems are associated with a specific bifurcation observed on the full model. First, diagnosis is initiated by performing a classification of the modeling complexity and of the various bifurcations, based on their kinds and classes. The dynamic behavior of a simulated stability problem depends on the modeling complexity used. Second, an Epoch-Based Trajectory Stability Assessment Methodology that is aimed at assessing asymptotic stability within time intervals free of hard limits discontinuities is presented. In this method, the stability assessment requires testing of quasi equilibria, limit cycles and trajectories within each epoch, rather than attempting to study stability of the entire trajectory. The third major step is to identify the smallest subsystem, a subset of the equations of the full system model that experiences the same bifurcation in the full system model. The bifurcation subsystem must satisfy both the linear and the transversality conditions “1,. Jitter“ .,,., 1... :1: ‘ V W ~13 tn J ‘ . h..-‘ 5.1 5. .. .e. ..Au 5' a». , .l.‘.. s, i . . :Eli \A- A Th: ~31" tier: Li 'ba K; W n, 4...; ”i ‘ »4:_ l.~‘J\J‘V L for the bifurcation that is being experienced in the full system model. The portion of the full system external to this bifurcation subsystem has no effect on the linear conditions for bifurcation to be satisfied in the full system model, and yet the equilibrium point is dependent on the variations in the full system equations which results in bifurcation in both the full system and the bifurcation subsystem. The bifurcation subsystem method utilizes the geometry associated with the various submatn'ces of the differential-algebraic Jacobian J and with the eigenvectors associated with the bifurcating eigenvalue to establish conditions for existence of such subsystems. Two examples for validating and using the bifurcation subsystem method to identify bifurcation subsystems in a model are presented, and the results are compared with right eigenvector and participation information. The bifurcation subsystem method can be in disagreement with both eigenvector and participation information. The bifurcation subsystem and the model dependency of test matrices T, KS and KD can identify the subsystem that experiences the bifurcation, the subsystem that produces the bifurcation, as well as possibly the location, cause and cure for the bifurcation, including loading changes or hard limit discontinuities. Copyright by Khadija Ben Kilani 1997 To my mother and father Fatma and Azzouz Ben Kilani 5.... ACKNOWLEGMENTS I thank Allah the most powerful for giving me the power and the patience to accomplish this dissertation. I am in debt to my parents who although far away, were with me all steps of the way in every way possible. I would like to express my sincere gratitude to my dissertation advisor Dr. Robert Schlueter for his continued academic advising, moral support, and constructive guidance. 1’ am also grateful to my committee members, Dr. Gerald Park, Dr. Hassan Khalil, and Dr. Jerry Schuur for giving me of their encouragements and their professional opinions. The financial support from the Detroit Edison company and the MSU Computing and Technology, as well the Graduate School made this work possible. vi TABLE OF CONTENT LIST OF TABLES ........................................................................................ xii LIST OF FIGURES................ ....................................................................... xiii LIST OF SYMBOLS ..................................................................................... xvii CHAPTER l INTRODUTION ......................................................................... l 1.0 Introduction to and Justification for Bifurcation Subsystems .......................................................... 1 1.1 Description of the Static Bifurcation Phenomena ..................... 9 1.2 Dynamic Bifurcation on Power Systems: A Brief Review ..................................................................... 10 1.3 Objective ............... 13 CHAPTER 2 POWER SYSTEM MODELING EFFECTS ON STABILITY PROBLEMS .................................. 19 2.1 Description of a Power System ............................................. 20 2.1.1 Synchronous Generator ............................................... 20 2.1.2 Electric Load ............................................................ 22 2.2 General Power System Model ................................................ 24 2.2.1 Model Equations ......................................................... 25 2.3 Modeling Complexity ........................................................... 26 2.3.1 Type 1 Model ........................................................... 27 2.3.2 Type 2 Model ........................................................... 27 2.3.3 Type 3 Model ........................................................... 28 2.3.4 Type 4 Model ........................................................... 29 2.4 Classes and Kinds of Bifurcations in a Type 2 Model ............................................................. 31 2.4.1 Kinds of Bifurcation .................................................. 34 2.4.2 Classes of Bifurcation ............................................. 38 2.4.2.1 Classes of Loss of Causality .......................... 39 2.4.2.2 Classes of Static Bifurcation .......................... 40 2.4.2.3 Classes of Hopf Bifurcation ........................... 41 vii CHAPTER 3 POWER SYSTEM STABILITY ASSESSMENT A TAXONOMICAL APPROACH .............................................. 42 3.1 introduction ........................................................................... 42 3.2 A Taxonomy of Stability Assessment in Power Systems ....................................................................... 45 3.3 Development of a Taxonomy for Power System Models with Hard Limits ....................................................... 48 3.4 Stability Security Assessment for Models with Hard Limits ................................................................... 51 3.4.1 Epoch—Based Trajectory Stability Assessment............... 55. 3.4.2 Bounded Stable Trajectory Feasibility Method .............. 62 CHAPTER 4 BIFURCATION SUBSYSTEMS IN A POWER SYSTEM DIFFERENTIAL-ALGEBRAIC MODEL .................................. 65 4.0 introduction ........................................................................... 65 4.0.1 Bifurcation Subsystems for Epochs .............................. 69 4.0.2 Existence of Bifurcation Subsystems ............................ 70 4.0.3 Introduction to Bifurcation Subsystem Theory .............. 73 4.0.4 Chapter Objective and Outline ...... y .............................. 75 4.1 Differential-algebraic Power System Model ............................. 78 4.2 Conditions for Bifurcation Subsystems Experiencing SN and LF Bifurcations ..................................... 82 4.2.1 Equivalent System Models ........................................... 82 4.2.1.1 Equivalent Algebraic System Model ................ 84 4.2.1.2 Equivalent Dynamic System Model ................. 85 4.2.2 Definition of Bifurcation Subsystems ........................... 87 4.2.3 Decoupling Conditions for Algebraic Bifurcation Subsystem and Differential Bifurcation Subsystems ................................................ 94 4.2.4 Geometric Conditions .................................................. 97 4.2.5 Algebraic Bifurcation Subsystem Experiencing LF Bifurcation ........................................ 105 4.3 Bifurcation Subsystems Experiencing SN and LF Bifurcation .......................................................... 106 4.3.1 Differential Bifurcation Subsystems for SN Bifurcation ...................................................... 107 4.3.2 Algebraic Bifurcation Sub-subsystem for LP Bifurcation ...................................................... 109 4.3.3 Differential-Algebraic Sub-subsystems .......................... 110 4.4 Differential Bifurcation Sub-systems for Hopf Bifurcation .............................................................. 112 viii CHAPTER 5 5.3 5.4 CHAPTER 6 6.0 4.4.1 Differential Bifurcation Sub-subsystems for Hopf Bifurcation ................................................... 114 Application of Algebraic Bifurcation Subsystem Experiencing LF Bifurcation in Power System Stability Assessment .................................................. 115 4.5.1 Description of a Load Flow Model .............................. 117 4.5.2 Structural Conditions for Utilizing J), ......................... 119 4.5.3 Comparison of Algebraic Equivalent Model, Algebraic Bifurcation Subsystem, and Load Flow Model ................................................. 120 4.5.4 Bifurcation in LP or in Generator Dynamics ................ 126 4.5.5 Conclusion ................................................................. 127 Effect of Hard Limit Discontinuities on Stability Analysis .............................................................................. 129 Introduction ......................................................................... 129 Linearized Muiti-machine Power System Model ...................... 132 Single Machine Infinite Bus Model with Local Load .............. 134 5.2.1 Model Representation .................................................. 134 5.2.2 Model Implications ...................................................... 136 5.2.3 Model Equations .......................................................... 139 5.2.4 Simplified Linear SMIB Model ................................... 140 5.2.4.1 The E’q Equation ........................................... 140 5.2.4.2 The electrical torque Equation ........................ 142 5.2.4.3 Terminal Voltage Equation .............................. 143 Application Examples of the Bifurcation Subsystem Method on a SMIB model ...................................................... 146 5.3.1 Differential Bifurcation Subsystem for SN (Unregulated Model) .................................................... 146 5.3.2 Differential Bifurcation Subsystem for Hopf (Regulated Model) ....................................................... 154 Effect of Hard Limit Discontinuities on Stability Analysis .............................................................................. 158 5.4.1 Synchronizing and damping torques ............................. 159 5.4.1.1 Unregulated System (Ge(s)=0) ......................... 160 5.4.1.2 Thyristor-Type Excitation System .................... 161 5.4.1.3 IEEE dc Type 1 Excitation System ................. 164 5.4.2 Bifurcation Tests for the SMIB ......... ~ ......................... l 67 5.4.3 Conclusions ................................................................. 170 'A Diagnosis for SN and Hopf Bifurcations in a SMIB Model .......................................................................... 173 Introduction ......................................................................... 173 HENDi) uPENDix that \J I‘ 6.1 6.2 6.3 CHAPTER 7 7.1 7.2 APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX El APPENDIX E2 Diagnostic Models Formulation .............................................. 174 6.1.1 Power Transfer Model (Rsh= X“: infinity) .................................................... 174 6.1.2 Local Load Model ....................................................... 175 6.1.3 Combination Model ..................................................... 176 Diagnostic Study .................................................................... 176 6.2.1 Power Transfer Model ................................................. 179 6.2.1.1 Effect of PM and OM .................................... 179 6.2.1.2 Effect of Transmission Network Inductance ..... 181 6.2.2 Power Local Loading Model ...................................... 183 6.2.2.1 Effect of real-inductive local loading ...................................... 183 6.2.2.2 Effect of real-capacitive local loading ..................................... 184 6.2.2.3 Effect of the transmission network .......................................... 184 6.2.3 Conclusion of the Power Transfer and Local Load Models ..................... . .......................................... l 85 6.2.4 Combination Local Load-Power Transfer Model ............ 187 6.2.4.1 Effect of Power Transfer ................................. 187 6.2.4.2 Effect of Local Loading ................................. 188 6.2.4.3 Effect of the Transmission Network ................ 189 6.2.5 Discussion of the Combination Model Results .............. 192 Diagnostic Summary .......................................................... . ..... 194 6.3.1 Classes of Hopf Bifurcation in a SMIB Model ............. 194 6.3.2 Classes of Static Bifurcation in a SMIB Model ........... 195 SUMMARY AND FUTURE RESEARCH ................................ 197 Thesis Summary ..................................................................... 197 Future Research ..................................................................... 202 POWER SYSTEM DYNAMIC MODEL .................................... 205 COMMON DEFINITIONS IN POWER SYSTEM STABILITY ............................................................. 217 DEFINITION OF HOPF AND SADDLE NODE BIFURCATIONS WITH TRANSVERSALITY CONDITIONS ........................................................................ 221 INPUT DATA FILES FOR THE BIFURCATION SUBSYSTEM EXAMPLES ..................................................... 223 FIGURES FOR DIAGNOSTIC ANALYSIS POWER TRANSFER MODEL ................................................. 227 FIGURES FOR DIAGNOSTIC ANALYSIS LOCAL LOADING MODEL .................................................... 231 i??E\DI.‘ ... wan ht Q.- ~. R fa. APPENDIX E3 FIGURES FOR DIAGNOSTIC ANALYSIS COMBINATION MODEL ........................................................ 235 REFERENCES ............................................................................................... 249 xi Table 5.1 Table 5.2 Table 5.3 Table 5.4 Table 6.1 (a) Table 6.1(b) Table 6.2(a) Table 6.2(b) Table 6.3 (a) Table 6.3 (b) Table 6.3 (c) LIST OF TABLES Computational results for identifying the bifurcation subsystem for SN bifurcation in the SMIB model - - ............................... 149 Participation factor results for the SN bifurcation in the SMIB model ....................... 153 Computational results for identifying the bifurcation subsystem for Hopf bifurcation in the SMIB model: Inertial mode -- 157 Participation factor results for the inertial Hopf bifurcation in the SMIB modeL--- - - -- ......................... 158 Results on the power transfer model bifurcation diagnosis ........................................ 181 Power transfer model non-bifurcation results .............................. 182 Results on the local loading model bifurcation diagnosis .......................................... 185 Non-Bifurcation Results on the local loading model Diagnosis ............................. 185 Results on the combination model bifurcation diagnosis ................. 190 Combination model: Effect of Local loading ............................................................. 190 Combination model: Non-bifurcation cases ............. - 191 xii LIST OF FIGURES Figure 1.1 Identification of the Critical Area in 230 kV buses under area System Stress ......................................................... 5 Figure 2.1 Functional Block Diagram of a synchronous generator excitation control system ........................................... 23 Figure 2.2 IEEE Type 1 Excitation System ................................................ 23 Figure 3.1 2-D Example portraying the effects of the i’th discontinuous control action. (xbyj) is the i’th quasi equilibria ......................... 45 Figure 3.2 Response of a system after a severe disturbance: EI-IV Line fault & line trip which resulted in a large generator off-line .............................................. _. ............. 61 Figure 3.3 Active power flow between central and south-east regions of the Swedish test system after the tripping of a line and generator in the North ................................................... 61 Figure 5.1 Block Diagram of a Power System Model .............................................. 134 Figure 5.2a General Configuration of a single machine connected . to a large system transmission line ............................................. 137 Figure 5.2b Equivalent system where the transmission network is reduced to its Thevenin’s equivalent ....................................... 137 Figure 5.2c Circuit Model of (b) where R5 is neglected ................................. 137 Figure 5.3 Simplified single-machine infinite bus mode. (a) Equivalent Thevenin model where rm xe and Vi are given by equations (5.1), (5.2) and (5.3) respectively. (b) Phasor diagram for the ‘ synchronous machine connected to an infinite bus shown in (a) 138 Figure 5.4 Block diagram of a simplified linear model of a synchronous machine connected to an infinite fa, Figure 5.5 Figure A1 Figure A2 Figure A3 Figure Bla Figure Blb Figure E1.1 Figure El.2 Figure E1.3 Figure 151.4 Figure E2.l Figure E2.2 Figure E23 synchronous machine connected to an infinite bus through an external impedance ............................................. 145 Block Diagram of a linearized machine model showing Bifurcation Subsystems and the test matrix T ............................. 153 Power system stabilizer model ................................................... 214 Speed-governing-turbine system model ....................................... 215 Speed-governing-turbine system model ....................................... 216 Response of a four-machine System during a Transient-, stable System ................................................................ . ........... 219 Response of a four-machine System during a Transient- unstable System ........................................................................ 219 Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter K5 . (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions ....... Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter [(4. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions Power Transfer Model: Effect of real and reactive power uansfer to the infinite bus on the parameter K1 -K2K5/K6. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter K, -K2K3K4. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions Local Loading Model: Effect of the shunt local load on the ............ 227 228 229 _ 230 parameters K4, (a) Inductive loading, (b) Capacitive loading ...................................... 231 Local Loading Model: Effect of the shunt local load on the parameter K5, (a) Inductive loading, (b) Capacitive loading ............. Local Loading Model: Effect of the local load on the parameters xiv 232 .e 5 Iii-3 p p . J. n... . p w r .1 ¥ \ h F: , . 3. nfl —s —.. F: ”in: u s e P.- . PM.” Figure 52.4 Figure 53.1 Figure E3.2 Figure E3.3 Figure 53.4 Figure E3.5 Figtn’e E3.6 Figure E3.7 Figure E3.8 Figure E3.9 Figure E3.10 K1 -K2K3K4, (a) Inductive loading, (b) Capacitive loading ........................... Local Loading Model: Effect of the local load on the parameters K. -K21(5/K6, (a) Inductive loading, (b) Capacitive loading“ - Combination Model: Effect of power transfer on the parameter K5 when the local load is real-inductive. (a) Power flow in the same direction, (b) Real and reactive power flow in opposite directions Combination Model: Effect of power transfer on the parameter K5 when the local load IS real-capacitive. (a) Power flow in the same direction (b) Real and reactive power flow in opposite directions - Combination Model: Effect of power transfer on the parameter K4 when the local load is real-inductive. (a) Power flow in the same direction, (b) Real and reactive power flow in opposite directions--- Combination Model: Effect of power transfer on the parameter K4 when the local load is real-capacitive. (a) Power flow in the sane direction. (b) Real and reactive power flow in opposite directions ....... Combination Model: Effect of power transfers to the infinite bus on the synchronizing term Kl «21(ng (regulated) when the local load is real-inductive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions - ....... Combination Model: Effect of power transfers to the infinite bus on the synchronizing term Kl {21(ng (regulated) when the local load is real-capacitive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. ............................... Combination Model: Effect of power transfers to the infinite bus on the synchronizing term K r 491ng (unregulated) when the local load is real-inductive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions -- Combination Model: Effect of power transfers to the infinite bus on the synchronizing term Kl -K2K3K4 (unregulated) when the local load is real-capacitive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions - - ................ Combination Model: Effect of local loading on the parameter K5 (regulated) when the uansfer to the infinite bus is at PmFOS. Qinf=05 (a) real-inductive loading; (b) real-capacitive loading ------ - - - Combination Model: Effect of local loading on the parameter XV ............. 233 234 ............ 235 236 237 .. 238 ............ 239 ............ 240 -_. ....... 241 ............ 242 ............ 243 Figs: E311 : F1 1 ‘v‘ij-\ r-.J~L—.-l- .‘all: E}. II -11 v iv". 1:53. ‘ Figure 53.11 Figure E3. 12 Figure £3.13 Figure £3.14 K4 (unregulated) when the transfer to the infinite bus is at PmFOB. Q-mf=0.5. (a) real-inductive loading; (6) real-capacitive loading. ......................................... Combination Model: Effect of local loading on the parameter K, -l(2l(3l(4 (unregulated) when the transfer to the infinite bus is at ..... 244 Pint-=08. Olaf-0.5. (a) real-inductive loading; (b) real-capacitive loading .................. 245 Combination Model: Effect of local loading on the parameter K , «21(ng (regulated) when the transfer to the infinite bus is at P-mFOB. Qinf=05 (a) real-inductive loading; (b) real-capacitive loading ................. 246 Combination Model: Effect inductive-real local load on the system synchronizing and damping torque coefficients when PmFOB. Obj-=05 (a) synchronizing torque coefficient KS. (b) damping torque coefficient KD ......... Combination Model: Damping torque coefficient Kd and parameter K5 for different levels of real local loading ............. xvi 247 248 LIST OF SYMBOLS Synchronous machine E’q D internal bus voltage proportional to d-axis flux linkage of field current damping constant due to damper winding effects (pu / (rad/sec) ) angular speed of the machine (rad/sec) nominal angular speed (rad/sec) terminal bus voltage mechanical power input to the rotor (pu) direct-axis (d-axis) current quadrature-axis (q-axis) current 1,, + j Iq angle difference between internal bus (rotor) and terminal bus q-axis synchronous reactance d-axis synchronous reactance d-axis transient reactance d-axis open-circuit transient time constant (sec) resistance of the stator Excitation system Vc output of load (or line drop) compensator xvii output of power system stabilizer (PSS) reference voltage output of automatic voltage regulator (AVR) amplifier output voltage measured voltage stabilizer output voltage dc gain of AVR Time constant of AVR (sec) dc gain of exciter time constant of exciter washout gain of stabilizing feedback of excitation system washout time constant of stabilizing feedback of excitation system output of the exciter (field voltage of synchronous machine) compensator adjustable inductive reactance xviii (I; we I Introduction 1.0 Introduction to and Justification for Bifurcation Subsystems An important aspect of stability analysis of a power system is to determine which physical elements or subsystems are associated with a specific instability phenomena observed on the full system. The recurrent problem in diagnosing a specific power system stability problem is that of size and complexity of the model since a power system may be composed of 10,000 interconnected buses and 1000 generators. However, although the models required to study the dynamic effects of interest in a power system seem to initially involve a large set of variables associated with the many different components of the system, only a subset of these state variables may actually be critical to problem diagnosis for a specific stability problem associated with a specific eigenvalue or mode. The challenge in then not to obtain a reduced-order model that incorporates all of the components in the system, but to obtain a subsystem of the full model which experiences and causes the same bifurcation (a discontinuous change in system response as a parameter changes slowly) as the full model, i.e. the smallest subsystem that preserves the qualitative information of the system instability involving the specific bifurcating eigenvalue. I‘d System modal reduction methods on a power system have mainly been based on modal analysis (or eigenanalysis) of the linearized power system model, a technique that has been found to be promising and of fruitful applications [1,2]. Sparsity, singular perturbation theory, coherency reduction as well as Kalman minimal realization theory [3] have also been applied in an attempt to simplify the large and complex dynamical system models that arise in various power system problems. These traditional methods for identifying subsystems are generally based on apriori knowledge of (a) the system or its control design; (b) time scale information on what system state variables are fast and which are slow variables; (0) eigenvectors and participation factor information from eigenvectors; or (d) decoupling within the state space model. While such methods are often useful, they are not always effective in identifying the smallest subsystem or subsystems that experience and cause the bifurcation in the full system because the criterion for deciding their effectiveness is their usefulness in model stabilization, designing controls, and developing operation schedules for a system. Obtaining a bifurcation subsystem model that experiences and causes the bifurcation in the full model requires a quite different set of methods. The singular perturbation approach for problem diagnosis in power system [4] can not generally be applied for problem diagnosis of a particular stability problem that occurs after a power system disturbance since the states most involved in an unstable mode in a power system cannot generally be classified into either the slow or the fast dynamics groups. When a power system is subject to undamped inter-area oscillations or local generator oscillations, this technique seems to be inadequate to establish the bifurcation subsystem for a single mode since the frequency of all inter-area oscillation modes lie in a i\ E narrow band and the frequency of all local generator oscillation modes also lie in a narrow band. A singular perturbation based subsystem model contains all the modes in the narrow band and not just the bifurcating mode. In modal analysis theory, the location and involvement of specific components of the system in the system instability may be obtained from the right and left eigenvectors associated with the bifurcating eigenvalue. While the magnitude of the eigenvalue of the normalized modal-transformed system is only a relative measure to proximity to bifurcation (collapse) [2], the components of the left eigenvector can be interpreted as indicating a direction normal to the operational boundary of the system [5-6], and right eigenvectors indicate the degree to which given variables are involved in a given mode. The dimensionless participation matrix, obtained from right and left eigenvector information associated with the bifurcating eigenvalue, has widely been used as an indicator to the relative involvement of system components in the unstable system mode or as a tool for identifying the collapse regions in the system. Although Eigenvalue/modal analysis has provided promising results, it seems to be insufficient to locate the voltage collapse regions or system components causing the collapse as can be demonstrated by the following example. Example Modal analysis has been used in [7] on the 1993 BC. Hydro system model to find the suitable locations where reinforcement and remedial actions ought to be taken. The relevant portion of the system in question as well as the three critical regions are shown in Figure 1.1. The system participation factor matrix for the critical mode of the normalized 5(1'1“ V {lu' u'lL 5 , . ilt‘sit’fl Lfitfifti [titers SIC). \lt‘Z'ldi imirlm mimic: [Anifx a system has been computed when the entire BC. Hydro system was stressed to determine the location where reactive power injection would benefit the system the most. Mth a uniform MVA load increase, modal analysis near the point of collapse has resulted in the participation factor plot shown in Figure 1.2. The modal information in the plot indicates that the British Columbia Lower Mainland was the critical area, with North Vancouver having the next larger participation factor elements. The Vancouver Island region in Western British Columbia had relatively small participation and thus was rated relatively ineffective for reactive compensation [7]. A study on the same BC. Hydro system has been conducted by Michigan State University where the reactive compensation device was a synchronous var compensator (SVC). The study based on MSU Voltage Stability Security Assessment and Planning Methodology Programs [8] concluded that the collapse region was caused by reactive imbalance in the Vancouver Island region. It has been found [9] that in this system, loss of excitation voltage control on generators in one generating station on Vancouver Island causes a voltage collapse in Vancouver Island which uncontrollably spreads and brings about voltage collapse in the entire BC. Hydro system. This result was later confirmed because ABB [10] installed one SVC and it was located in the Vancouver Island region at the site suggested by [9]. 5; Participation Factor ('0.001) ‘ N d 0.5 Figure 1.1 Buses Ordered in the Sequence of P Factor Identification of the critical areas in 230 kV buses under area system-stress The al‘ n01 suifidrn fail sixifm- illts huiei rues ulth‘ :13 mi be 1613 8:03-36 Oiliti bum submtdcl a ital tiiupi Baillie pm thpoutr x; more of chargers it susepante 33311113 Will Constraints, lifercntialz dilliO iis’u‘iilti -.-. ‘ "In? . , j‘QNIfntl .al-a Stir. il‘ \ iii Will's-iii; The above example demonstrates that the participation matrix can be misleading or not sufficient for indicating the portion of the system that is causing the instability in the full system, since (i) the voltage collapse can occur due to discontinuities in regions where states have quite small participation for the mode closest to bifurcation or (ii) out of many states with “high” participation, only few may be crucial for remedial actions. This result may not be surprising since for most realistic power systems (a) the Jacobian matrix is very sparse where a node in a 10,000 bus network may be connected to only three to five other buses in the system; (b) weak coupling exists between the real power-angle submodel and the reactive power-voltage submodel in the algebraic network model; (c) weak coupling of coherent bus groups exists in both the real power-angle and in the reactive power-voltage network submodels; (d) decoupling in the different dynamics of the power system model exists and (e) most bifurcations in a power system occur due to a sequence of discontinuities such as generators reaching field current limits, under load tap changers reaching tap position limits, and switchable shunt capacitors reaching susceptance limits. All of these factors suggest that loss of control in some coherent bus groups results in decoupling between some generator dynamics and system algebraic constraints, and thus will produce a bifurcation for the algebraic, differential, or differential-algebraic model associated with that coherent group. This has been demonstrated in [11] for a classical load flow model and is expected to be true in a differential-algebraic model. Selective Modal Analysis is a model reduction approach that recognized that the participation matrix is not sufficient for obtaining the smallest Subsystem that preserves the qualitative information of the system near collapse. In the modal analysis algorithm 0 re I g- 9' '9“ 1.1L .‘JI‘ I L'ttic‘i if. '. . a 1. mill: “'42 ab "*7 litr [12], state variables with high participation are divided into “relevant” and “less relevant” states where the relevant states are preserved in the reduced order model, while the effect of the less relevant states is incorporated in the final reduced model [12]. The relevant states are obtained either by intuition or from participation factor information. This approach demonstrates that only a few out of the many states with high participation factors in a model may be of crucial importance in a particular unstable mode. Center manifold dynamics may be a subsystem of the actual system that produces and causes the bifurcation experienced in the full model. The center manifold dynamics are obtained through a nonlinear u'ansformation of the original state space model. The center manifold dynamics are not considered a bifurcation subsystem despite the fact that it manifests the bifurcation similar to a bifurcation subsystem, because (a) a bifurcation subsystem is a subsystem of the state space model written in terms of physically meaningful variables and not a nonlinear transformation of these dynamics; (b) the nonlinear transformation on very large power system models generally makes it impossible to analyze how operating condition changes or model parameter changes will cause the bifurcation; (c) the needed nonlinear transformation is not quickly and easily computed, which is not the case for finding a bifurcation subsystem. (d) a bifurcation subsystem is not a reduced model that can necessarily be simulated or analyzed apart from the full system model. A bifurcation subsystem is a subset of the equations of the full system model that experiences, produces and causes the full system bifurcation by satisfying the same L71, - '- C‘s.» conditions for bifurcation as satisfied by the full system model. The subsystem that experiences bifurcation may or may not produce the bifurcation in the full system since the structure or parameters external to the subsystem may help produce the bifurcation in both the full system and the subsystem. This is certainly true when geometric coupling of the bifurcation subsystem and external system approaches zero as bifurcation occurs rather than being zero. If the geometric coupling between the subsystem experiencing bifurcation and the external system is zero as bifurcation occurs, which is called a structural geometric decoupling, then the bifurcation subsystem both experiences and produces the full system bifurcation. Thus the subsystem that produces the bifurcation experienced by the bifurcation subsystem and the full system may be of larger dimension than the subsystem experiencing bifurcation. The subsystem causing bifurcation may be still larger than the subsystem experiencing or producing the bifurcation in the full system and the bifurcation subsystem. The system that causes bifurcation must include all model parameters, controls, and discontinuous changes that cause the bifurcation to develop. The discontinuous changes include under load tap changer tap position changes, switchable shunt capacitor switching changes, equipment outages, over-excitation limit relays, etc. that cause instability to occur. The controls include maximum excitation limiter controls, under load tap changer controls, switchable shunt capacitor controls, and excitation system controls. Drift“. 3m stern Patter 5 m a his ncx Una}? .e If “a. 1.1 Description of the Static Bifurcation Phenomena In a power system, the parameter space breaks up into open connected regions, each characterized by a particular behavior of the system in the state space. As parameters move across the boundary from one open region to another, the topology of the state space changes. This corresponds to the bifurcation phenomena. Although the bifurcation can take various forms, a great concentration of the study of bifurcation in power systems has been devoted to the study of saddle node bifurcation. It is only recently that the power system researchers have recognized the significance of other forms of bifurcations in power system stability analysis above and beyond the static bifurcations. The creation of limit cycles in the state space through Hopf bifurcation, although not fully investigated, has nevertheless been a familiar although undesirable phenomena in power systems. Saddle node bifurcation is a kind of static bifurcation that occurs at the steady state operating point of a dynamic system. Significant research on static voltage stability [App. B] has resulted in (a) establishing the modeling requirements for static voltage stability, (b) classification of the various kinds of voltage instability such as algebraic bifurcation, loss of causality, and saddle node bifurcation. Each of these different kinds of bifurcations can affect different algebraic and/or differential states of a model and each is called a class of that particular kind of bifurcation. For example, a saddle node bifurcation can occur in generator flux decay dynamics or generator inertial dynamics. Each of these is a different class of saddle node bifurcation. A loss of causality is associated with singularity of the Jacobian of the algebraic equations at an equilibrium. The Jacobian can be singular due to row dependence of the real power balance equations and is referred to as angle instability, or due to row dependence of the reactive power balance equations and is termed a voltage ' r;- TI.» 1‘ Iii 5.11., L r ‘ \ggnil 0min '. 3 I" t 3\ ell-D i l eu.uo¥mu 7‘" vv -r;'~ ~huAeuL4 .- . 3 ‘1‘)[55 AI «I: ' fi.‘ “fix/311.13 ..',,. mniirid ‘ "‘V r. “‘CJLUIC P‘i’i CI 8}: fill-”frig- lilacs ll Willem 'm. Kim-3d ‘1' 10 instability; both are classes of loss of causality. Although there is a theory for establishing which kind of these bifurcations occurs, there is no general method for defining subsystems associated with classes of bifurcation for each of these kinds of bifurcation. There is no actual test that identifies whether a subsystem actually causes the bifurcation. A method will be developed that (i) identifies subsystems that experience and cause saddle node, load flow, and Hopf bifurcations and the classes of each; (ii) provides tests for whether that specific subsystem causes the static bifurcation observed in the system, which in turn may identify the kind and class of static bifurcation which occurs; and (iii) establishes a diagnostic procedure for determining the causes in terms of the specific parameters in parameter space and the model structure discontinuities that produce this subsystem bifurcations. The bifurcation subsystem method will also be developed for Hopf bifurcations in power systems where less literature exists. A brief review of the literature on dynamic bifurcation in differential power system models is now given. 1.2 Dynamic Bifurcation on Power Systems: A Brief Review As a continuing effort to the understanding of voltage stability [App B] mechanism in a power system, recent attention has been focused on dynamic bifurcation where periodic solutions suddenly appear. It is now evident that poorly damped or unstable low frequency oscillations are a major threat to the security. of power systems for many utilities. The modes of oscillations observed are becoming more numerous and complex in mOdern interconnected systems. Three specific classes of Hopf bifurcation have been identified [25]: The local plant mode is associated with the mechanical dynamics of units 11 at a generating station swinging with respect to the mechanical generator dynamics of the rest of the system and the frequency of such modes is typically in the range of 0.8 to 2.0 Hz. The inter-area mode is associated with the swing of the mechanical dynamics of one group of machines in one part of the system against the mechanical dynamics of groups of machines in other parts of the system. The frequency of the inter-area mode is typically in the range 0.2 to 0.8 Hz [13]. The voltage mode appears in any machine’s flux decay dynamics and exciter dynamics. Its frequency is higher than the local plant and inter-area modes. In 1984, Abed and Varayia investigated generating units dynamics including both generator electro-mechanical and flux dynamics [14]. Their investigation revealed for the first time that generating unit dynamics may experience Hopf bifurcation. In 1989, Sauer and Pai [15] investigated the stability of a two-axis generator model connected to an infinite bus as load level at the generator terminals increased. The study showed that near the critical or maximum possible load level, a pair of complex eigenvalues associated with the excitation system experienced Hopf bifurcation, and that the onset periodic solutions were locally unstable. From 1990-1992, similar investigations have been carried on by others [16-19], where system stability was based on eigenvalue analysis at different loading conditions: The oscillatory instability in all of this literature was due to an apparent subcritical Hopf bifurcation in the generator flux decay and exciter dynamics. During 1992-1993, Ajjarapu, Abed and Varayia [20-23] conducted a thorough testing on a similar single machine—infinite bus voltage collapse model developed by Dobson et al [24] in 1990. The results revealed that Hopf bifurcation and saddle node bifurcation occur on this model at different system loading levels. The oscillatory instability was characterized in will 3\ )l LDC in; “(Lil PCSsl ul“ 12 by the interaction of the generator angle dynamics and the induction motor dynamics. A recent study on bifurcations on system models of increasing complexity has been completed by Ajjarapu [49]. Two models were of interest: one model included a detailed machine model with IEEE type DC-l exciter, and network with a constant power load model, and a more complex model where a generic dynamic load model has replaced the constant power load model. Eigenvalue analysis on both models revealed that before the system collapses due to saddle node bifurcation (which occurs at the tip of the PV curve), the system experienced Hopf bifurcation, node focus bifurcation and singularity induced bifurcation. The sequence of Hopf, node focus, singularity induced and finally saddle node bifurcations occurred on both a 3-generator 9-bus system and a IO-generator 39—bus New England system model. Issues that need additional research include (1) identification of generic sequences of bifurcations in a power system; (2) identification of a bifurcation subsystem for each possible class of saddle node, Hopf, node focus, or load flow bifurcation in the system; (3) development of tests for establishing whether a specific subsystem causes a specific bifurcation observed in the system; and (4) a diagnostic for establishing the specific cause of that subsystem bifurcation. Such a subsystem whose bifurcation causes bifurcation in the full system is called a bifurcation subsystem, and may be composed of a subset of solely differential, algebraic, or a combination of differential and algebraic equations. The method presented in this thesis establishes conditions for which a specific bifurcation subsystem exists. This subsystem method is first developed in Chapter 4 for (1) a general differential-algebraic model for identifying both differential and algebraic bifurcation subsystems; then (2) for an algebraic model for identifying algebraic bifurcation sub- (1) l3 subsystems; (3) to a differential model for identifying differential bifurcation sub- subsystems; and (4) to a differential-algebraic model for identifying differential-algebraic bifurcation subsystems. 1.3 Objectives 1. Initiate a diagnostic classification of power system stability problems that is based on bifurcation kinds and classes for the different model types that have been used for bifurcation studies. The necessity for simplifying assumptions used in power system stability analysis has resulted in four types of power system models of increasing complexity [25], and which are commonly used to study the different forms of power system instabilities. Kinds of bifurcation refer to saddle node, Hopf, LF bifurcation, etc., thus refers to the bifurcation behavior observed, while a class of bifurcation refers to the specific dynamical or algebraic states experiencing the bifurcation. It has been shown in [25] that the type of model used for assessing stability has an effect on the kind and class of bifurcation observed. This classification based on model types as well as kinds and classes of power system stability problems seem to be necessary to diagnose the problem and to obtain information on how system stability can be improved most effectively. This diagnostic classification together with a description of a power system differential-algebraic model is given in Chapter 2. Briefly review in Chapter 3 the recent development on power system stability which extends the classical theory of power systems as dynamical smooth nonlinear systems. to the ongoing development of a method that includes the effects of hard limits and equipment outages which discontinuously modify the model order and system 14 behavior. The theoretical framework for study of stability of a power system model, called a “taxonomy”, has been developed for smooth systems by Zaborsky [26]. The taxonomy addressed three regions that lie in state and parameter space: (i) stability region composed of the region of attraction defined for every parameter vector in parameter space; (ii) a feasibility region in state and parameter space composed of all connected stable equilibria and the parameter vectors that produce them; and (iii) a viability region composed of a subset of the equilibria and associated parameter vectors belonging to the feasibility region that lie within the acceptable range for each current, voltage, real power and reactive power variables. Since hard limits associated with equipment outages, actions of operator control, protective devices such as relays and field current limit controllers, actions of under load tap changers and actions of switchable shunt capacitors cause most of the bifurcations experienced in a power system, the taxonomy is being extended to include such models. The taxonomy that includes the effects of hard limits [27] defines and describes similar stability, feasibility and viability regions as the taxonomy for smooth systems. However, the taxonomic assessment is only valid if one can predict and compute the sequence of quasi equilibria or limit cycles as well as the sequence of hard limits encountered along every system trajectory without simulation. This requirement and the assumptions made in [27] are not true in a power system with hard limits, where there is different time delay between the time that a specified hard limit surface is encountered and when the action associated with the hard limit is implemented. A different approach to stability assessment methodology for systems with hard 15 limit discontinuities is then proposed in Chapter 3. This approach considers an episodal trajectory stability assessment rather than attempting to study the stability of all trajectories and stability of all quasi equilibria and limit cycles encountered in some subregion of state and parameter space. This new approach is aimed at assessing asymptotic stability of the transient during (5.1;) and the stability of quasi equilibria or limit cycles in each epoch, where an epoch is a period of time (if, if” ) where no hard limit transitions occur. This procedure can be applied to a broader class of models with hard limits that need not obey the assumptions made in [27]. The Epoch Based Trajectory Stability Assessment Method is based on a stability simulation method implemented by Van Cutsem in 1992. [28] which does not require time simulation in (135“) , but only computes the quasi equilibria (El-,9.) for i = 1,...,N in a manner similar to how a classical load flow simulates the effects of hard limits. Specific rules and guidelines decide the sequence of discontinuities that occur using a procedure that avoids full time simulation that in turn determines the sequence of quasi equilibria (2.3.) for i = 1..... N which can then be used to test for stability, feasibility and viability via trajectory segments in epochs (if, rg+ l) , i = 1... N. A Bounded Trajectory Feasibility Region is then defined in Chaptm 3, which classifies trajectories produced by continuous parameter changes or by discontinuous parameter or model structure changes that are induced by continuous operating condition or parameter changes. Tire classification of trajectories is based on (i) whether each quasi equilibria (21.3.) or encircling orbit is stable and (ii) whether each trajectory segment in (if, 1.; 1) belongs to the Region of Stability of the associated quasi stable equilibrium point or limit cycle. The subset of trajectories, where (a) all l6 equilibria or encircling orbits are stable and (b) all trajectory segments (or only the last trajectory segment) lie in the region of stability of a stable orbit or equilibrium, belongs to a Bounded Trajectory Feasibility Region. Since the episodal stability assessment method is applied after each hard limit is encountered, the power system model in each interval may be different and yet analytic so that existing dynamical system theory can be applied. A natural desire is to find a reduced—order model or the smallest subsystem within each epoch that experiences the same bifurcation as the full model, which motivates the bifurcation subsystem method investigated in the rest of the thesis. If the bifurcation subsystem is identical in each epoch as discontinuities bring on the bifurcation, it is a bifurcation subsystem for the bifurcation that ultimately occurs and changes the system dynamical behavior. A power system differential-algebraic model is investigated for the purpose of identifying the possible bifurcation subsystems for each kind of bifurcation (SN, Hopf, LF bifurcation) that may occur on a Type 2 power system model where load dynamics and controls are ignored. A Type 2 model with a constant power load model, tap changers, switchable shunt capacitors in the transmission, subtransmission and higher voltage distribution system will be assumed in this thesis because it generally contains the dynamics that experience bifurcation and cause voltage collapse in a large region [25]. The bifurcation subsystem method allows one to find a specific subsystem that experiences, produces and causes the same bifurcation that occurs in the full model. The bifurcation subsystem always experiences the bifurcation that occurs in the full system model. The geometric coupling of the bifurcation subsystem may not be null l7 and may produce the bifurcation in the full system and the bifurcation subsystem as it becomes zero. The system that causes the bifurcation in the bifurcation subsystem contains all the discontinuities, operating changes, equipment outages and model parameter changes that cause the bifurcation to develop in both the full system and the bifurcation subsystem. Thus, this subsystem may be larger than the bifurcation subsystem that experiences the bifurcation of the full system model. or the subsystem that includes the effects of coupling to the external system when this coupling is not null but approaches zero as bifurcation occurs. The bifurcation subsystem method is an extension of the diagnostics employed in producing the taxonomy for feasibility, stability and viability for an equilibrium or closed orbit developed in Chapter 3, the classification of kinds and classes of bifurcation given in Chapter 2, and the episodal epoch by epoch stability assessment given in Chapter 3, because it finds the smallest subsystem that experienced, produced and caused the loss of stability, feasibility and viability for that equilibrium for bifurcations produced by encountering hard limit induced discontinuities. The bifurcation subsystem method allows one to better diagnose the causes of a subsystem bifurcation and separate causes and effects of this bifurcation for bifurcations brought on by hard limit induced discontinuities. Since the sequence of discontinuities are the major cause for many bifurcations, a bifurcation subsystem method based on epoch by epoch diagnosis appears to be the only way to understand the causes and the bifurcation subsystem that produced the bifurcation. This is needed in a power system model when there are hundreds of generators, tens of thousands of dynamic loads, and a network that couples the load and generator dynamics. The bifurcation subsystem l8 method should identify bifurcations that occur in (a) algebraic subsystems composed of generator bus groups, groups of load buses, and groups of generators and load buses; (b) differential subsystem dynamics in one or more groups of generators; (c) differential-algebraic subsystems in one or more groups of generator and in one or more groups of generator and load buses. The bifurcation subsystem method, developed in Chapter 4 would allow a very complex behavior to be seen as being very simple so that one can potentially explain virtually everything about feasibility, viability and stability for each and every subsystem that can produce loss of feasibility, loss of viability or loss of stability. The bifurcation subsystem model is produced to study the causes of a specific bifurcation and may have no value except for studying that bifurcation. It may not be useful for control design or even simulating the contingencies which produce that bifurcation, because the subsystem that produces the bifurcation may be quite small compared to the subsystem needed to simulate a specific contingency or design a control system to prevent that bifurcation. In chapter 5, two examples for validating and using the bifurcation subsystem method in a model are presented, and the results are compared with the participation factor information. In Chapter 6 a diagnostic study for the occurrence of SN and Hopf bifurcations in a single-machine-infinite bus model is performed based on test matrices derived in Chapter 5. Chapter 7 summarizes the dissertation and suggests topics of future investigation. II Power System Modeling Effects on Stability Problems In this chapter, a diagnostic classification of power system stability problems that is based on bifurcation kinds and classes is developed, and a description of the model types that have been used for bifurcation study is initiated. A specific model type (Type 2) is selected for further study in this thesis. A further classification of kinds and classes of bifurcations for this Type 2 model is then given. Kinds of bifurcations include saddle node, H0pf, load fiow bifurcation, thus refer to the kind of qualitative change in dynamical or algebraic network behavior encountered at the point of bifurcation. Classes of bifurcations are the specific subsystem dynamics or algebraic network variables that experience a particular kind of bifurcation for a given model type. For the model Type to be used (Type 2 model), there are two classes of saddle node bifurcation, three classes of Hopf bifurcation, and three classes of loss of causality bifurcation. The effort to classify model types as well as kinds and classes of bifurcation is necessary since later, a bifurcation subsystem method will be developed and applied to saddle node, load flow and Hopf bifurcations and the classes for each. This bifurcation subsystem method will attempt to diagnose the geometric subsystem structure that bifurcates. what locations in a system are vulnerable to this particular kind and class of 19 l-J if] ngr .‘k. 20 bifurcation, what causes it in terms of model parameter changes, equipment outages, or Operating condition changes. This systematic diagnosis needs a classification of kinds and classes of bifurcation in a given model type. Details of the classification to be given in this Chapter comes from a paper published by the author [25] in 1994. As a natural start to this diagnostic classification, a general power system model description including models for generators, transmission network and loads are given in Sections 2.1 and 2.2. Section 2.3 presents the four model types suggested by the author in [25], and finally kinds and classes of bifurcations in a Type 2 model are presented in Sections 2.4. 2.1 Description of a Power System Electric power systems are comprised of generation, transmission and distribution/ load. Electric power is produced at generating stations and transmitted to consumers through a complex network of individual components including transmission lines, transformers, and switching devices. Power systems models, although they vary in size and complexity, all have three basic structural components: synchronous generators, a transmission network, and electric load. 2.1.1 Synchronous generator Prime movers convert the primary sources of energy (fossil, nuclear and hydraulic) to mechanical energy that is in turn converted to electrical energy by synchronous generators. The synchronous generator output voltage is usually regulated to maintain a constant value using an excitation system. A functional diagram of a synchronous generator with an excitation control system is shown in Figure 2.1 [29]. The subsystems 21 inside the dashed-line box in Figure 2.1 constitute the excitation system, of which a typical model widely used in literature is the IEEE Type DCl excitation system model shown in Figure 2.2. For consistency and task simplification, this specific model will be used in this work. First, a brief description of the indicated functional blocks is given. (a). (b) (C) (d) (e) Exciter: provides dc voltage to the synchronous machine field winding based on input from the regulator control; Regulator. Processes and amplifies input control signals (as shown in Figure 2.1) to a level and form appropriate for the control of the exciter. This includes both regulating and stabilizing functions (rate feedback or lead-lag compensation); Terminal voltage transducer and load compensator: senses generator terminal voltage, rectifies and filters it to do quantity, and compares it to a reference which represents the desired terminal voltage. In addition, load compensation may be provided if it is desired to hold constant voltage at some point electrically close to or distant from the generator terminal bus; Power system stabilizer: provides an additional input signal to the regulator to damp oscillations in generator mechanical dynamics. Some common input signals are rotor speed deviation, accelerating power, and frequency deviation; Limiters and protective circuits. They include a wide array of control and protective functions whose role is to insure that the capability limits of the exciter and synchronous generator are not exceeded. A recent textbook [29] provides a detailed description of the various limiters and protective devices commonly used on power systems. I») 1‘») 2.1.2 Electric load The diversity of electric loads present on a power system, makes it virtually impossible to represent every individual component in a power system model. The common practice is to make considerable simplifications of the composite load characteristics as seen from bulk power delivery points. Load models have been broadly classified into two categories, static loads and dynamic loads. The static load model expresses the load as algebraic equations, classically in the exponential form - b P P0(V) Q = 9007)” where P and Q are active and reactive components of the load and the voltage magnitude is V. The subscript “0” refers to the values at the initial operating condition. If the parameters a and b are 0, 1 and 2, the models are respectively constant power, constant current and constant impedance load. The dynamic attributes to motors are usually the most significant aspects of dynamic characteristics of system loads. Typically motors consume 60-70% of the total energy produced by a power system. Other devices that require dynamic consideration are load shedding, operation of protective relays, thermostatic control of loads, response of under load tap changers (ULTC). etc... 7?” ~_—. I. ' t ”I ‘1 Limiters and ‘. Protective circuits .............................................. Term' a1 volta e tran ucers a load compensator - . To wer —i Excrter Generator -u—> sysigm —— v Ref Regulator v C Figure 2.1 Functional block diagram of a synchronous generator excitation control system I < V REF Amplifier KA/(1+sTA)-4| 7L . i—E-“l 0(5) l—q ' Exciter inertial and ‘ flux decay iliz MW“ sKF/(l+sTF V I l/(l-i-STR) fi—L Vc= lV-i-jxcl lF Measurement Line p Compensator Figure 2.2 IEEE Type DC 1 Excitation system model 7’ Gent llic g! xuhtttto d§irnu.x.l tie 31;:ch bUsin the grncrriors ‘ a PILL. I!!! j (7) Val; j«‘mli Par. SM 24 2.2 General Power System Model The general power system model consists of a set of first order differential equations, subject to algebraic constraints. The differential equations represent the mechanical dynamics, fiux decay dynamics, and excitation system dynamics of a generator whereas the algebraic equations represent the real and reactive power balance equations for each bus in the network. A multi-machine system consisting of m network buses and n generators with voltage regulation can therefore be written in the parameter dependent differential-algebraic model where x(t) W) 110) XxY I R I R i = f(x(t).y(t).u(t)) Mt) e P 2.1(a) 0 = 800). YO). 9(1)) 210’) n + m + p -> R" , an n-vector of real smooth nonlinear functions. '1 + m + p —; Rm, an m-vector of real smooth nonlinear functions. e X cR" , a state vector of dynamic state variables. typically time dependent values of generator voltages, rotor phases and induction motor variables. : 6 Y c Rm , a state vector of instantaneous variables, typically time dependent values of voltage and angle of terminal buses, high side transformer buses and load buses. : e P c R” , a vector of slow varying operating parameters. It can represent structural parameters like system components, reactances, under-load tap changers and switchable shunt capacitors, or operational parameters like reactive power load or generation, voltage set points on voltage control devices, as well as generation dispatch. As it varies over the parameter space F, it is assumed that there must be at least one solution to the above differential-algebraic system of equations. : state space P : parameter space. 2.2.1 Model Equations The differential-algebraic model consists of a) Differential equations d5i . W = (”i—(Os [=1 ..... n dcol. Mi—d-t— = Pmi—Di((i)i-O)S)—Pei r=l,...,n d5". Tdoi dt = “Eqi—(Xdi-Xdi)ld;+Efdi i=1.....n dE . J '_ _ -_ TE. m _ (KEi+sEi(Efdi))Efdi+vRi r-l.....n Tfl- —V +V +'(lx) i-l n Rid: " 11' In JTicil — ..... T dV3i _ v KFIVRi KFiEfdi(SEi(Efdi)+KEi) _1 Pi dt - _ 3i..- T - T l- ,. n 51' El T {-l—V-B-i v +K v v v '-1 M d: ‘ Ri Ai( REFi’ ii’ 31') l- ~~~~~ n for the generator and exciter dynamics where the variables and constants are defined in the list of symbols. b) Stator algebraic equations Vicosei + Rsi(ld‘.sin8i + Iqicosfil.) — X'diuqisinbi ~1dicosfii) _[E'disinfii+(X'qi-X'di)lqisin8i+E' icosbi] = 0 q —[E'qisin6i—(X'qi-X'di)lqicosbi+E'dicosbi] = 0 t: I, ..... n that define the lqi and 14,- at the generator terminal buses. lll‘ll Cl: _. 1 (I: .5” 3 1.11 crab: Where SUSCCI 2.3 dimer, em in: t. ' WP: ('73 26 c) Network (Power Balance) Equations Terminal bus equations ,- = Villdisinwi - 9i) + lqicos(5i — 91.)] + PLi(Vi) i: 1..... n Q.- = Viildicosw. — 0,) + Iqisinrsi - 991+ QLi(V,-) i = 1..... n that define voltage and angle at generator terminal buses where PLi(Vi) and Q ”(Vi) are the loads at the terminal buses. Load and high side transformer bus equations The net real and reactive powers at each of the network load and transformer buses can be expressed in terms of the voltage magnitudes, conductances and susceptance as, n pl. = —PU(V'.) = kg, lViv/((Gikcostiik+Biksin0ik) 1=n+l,...,n+m n Qi = -QLl-(Vi) = 2 VivkmikSineik—Bikcoseik) i=n+l,...,n+m k = l where PLi(Vi) and QLi(Vi) are the loads at bus i and G ik and Bik are the conductance and susceprance of the ik’th element of the Y-bus. 2.3 Modeling Complexity In dynamical system theory, the power system is generally a very high-order, high dimensional multi-variable process, that is operating in a constantly changing environment. Simplifying assumptions are therefore a virtual necessity for developing and using power system analytical methods. A classification of four types of power system models of increasing complexity, that are used in the existing literature to study or simulate. the various forms of power system instability is presented. A model type has been defined in [25] by the author as the simplest model that captures observability, €0th de. el all}. [if 2.3) ‘- 27 controllability, feasibility region, region of attraction and modal behavior associated with at least one class of stability problem. This paper [25] discusses the known classes of stability problems for each model type. 2.3.1 Type 1 model Type 1 model reduces a complex differential-algebraic model that includes dynamic models of both generator and loads, to an algebraic model of the generator, network and load, because the dynamic states are assumed at equilibrium. However this model takes into account the time sequence of tap changing transformer, switchable shunt capacitor and low voltage relaying actions as well as secondary voltage control, AGC and over- excitation limiter (OEL) control actions. A classical load flow model approximates a Type 1 model at steady state ( t=infinity) and often incorrectly estimates the sequence of control actions of under load tap changers, switchable shunt capacitors, and field current limit controllers after a contingency or due to slow operating condition changes. The recently developed simulation tool of Van Cutsem [28] is an example of a Type 1 model if relaying actions are represented. 2.3.2 Type 2 model Type 2 model is perhaps the most common power system dynamic model used to simulate and study instability. It currently includes generator inertial, flux decay, exciter, power system stabilizer, and boiler turbine governor dynamics. However this model excludes all load dynamics, for they are assumed to be decoupled from the generator dynamics, and thus do not contribute to the stability problems observed on this model. The 28 general form of the Type 2 model is: x = f(x.ld.lq.V.O) 2.2(a) O = ld—gd(x,V.O) 2.2(b) O = lq—gq(x,V,O) 2.2(e) 0 = PL—gP(x,V,O) 2.2(d) 0 = QL-gQ(x,V,O) 2.2(e) where x : generator dynamics P L, Q L : real and reactive power load v, 6 : voltage and angle of generator terminal, high side transformers and load buses. 2.3.3 Type 3 model In this model, the generator dynamics are modeled by constant power sources, while the load modeling includes all of induction motor and under load tap changers (ULTCs) controller dynamics, as well as dynamic models of switchable shunt capacitors, and thermostatically controlled loads. An extensive use of this model is for studying radial voltage instability, since by definition, this instability problem occurs at points electrically decoupled from the generator dynamics [30]. The vector x representing the generator dynamical states is either a constant vector or not included in the state space model equations. Therefore the complexity of this model depends on the complexity of the dynamic model used to represent the induction motor dynamics. A Type 3 power system model has the following form (2.3a-e), where the induction motor model used is a 3rd order model recently derived by Sauer and Pai [31] as an implicit load dynamic model, although other dynamic load models may alternatively be used. it here SID-d)’ s trough 23.4 decay, the 103 M} an. 29 o = PL—gp(x, v, e) 2.3(a) o = QL—gq(x, v,e) _ 2.3(b) o = h,(PL, QL, v, m,‘%, %,%) 2.3(c) o = he(m, 559,-) 2.3(d) o = hm(PL, QL, v, 03.9%, ‘2—‘f,‘%, ‘13-?) 2.3(e) where V, q and to are N dimensional state variables. A type 3 model ought to be used to study saddle node bifurcation in load dynamics and saddle node and Hopf bifurcations brought on by under load tap changers and switchable shunt capacitors. 2.3.4 Type 4 model This model combines the complexity in generator modeling of a type 2 model and in load modeling of a type 3 model. The generator model includes generator inertial, flux decay, exciter, power system stabilizer, and boiler-turbine-generator dynamics, whereas the load model may be a separate dynamic model for real and reactive power, or a composite model of real and reactive power dynamics. A Type 4 model has the form 2.2( a-c) and 2.3(a—e) if a nonlinear implicit load model is used, although the load model equations 2.3(c-e) may differ. Type 1-3 models will later be viewed as specific subsystem models for the Type 4 power system model. This will be made clear when differential subsystems are discussed in Chapter 6. Thus it will be clear that a Type 4 model is quite general and that Type 1-3 models are only used when specific types and classes of bifurcations are to be studied. A Type 2 model may be able to capture the generator dynamics that ultimately experience saddle node bifurcation, Hopf, or node focus bifurcafion, and which may be {a 51 s n exp U011] 01 t’.’ iriild tier: dynun couch litlie 3() affected by load flow bifurcation. algebraic bifurcation and singularity induced bifurcation. A Type 2 model with a constant power load model. tap changers, switchable shunt capacitors in the transmission, subtransmission and higher voltage distribution system will be assumed in this thesis because it generally contains the dynamics that experience bifurcation and cause voltage collapse in a large region, a utility, or even a county. However, a Type 2 model can not accurately capture the load buildup and overshoot but just approximates the peak load produced by action of distribution system load control dynamics and the voltage control dynamics. A Type 4 model is the only model that can hopefully simulate actual voltage collapse events. For Hopf bifurcation studies, a Type 4 model must be used because both the load dynamics and generator dynamics must be considered Load dynamics may strongly couple to generator dynamics when subcritical (unstable), or low frequency supercritical (stable) Hopf bifurcations occur. In fact, when a pair of stable conjugate eigenvalues approaches the jun-axis, the real part of the eigenvalues also decreases, which causes the orbit size of the associated limit cycles to increase. This causes more control actions of under load tap changers, switchable shunt capacitors and thermostatically controlled loads to take place to prevent the system equipment from operating beyond the maximum tolerable range. The impact of the load dynamics is therefore present and a full Type 4 dynamic load model is necessary to represent them. Although Type 4 models are needed, they have generally been developed only for small regions within a utility or are only being developed at the present time for larger regions. Since a Type 2 model with voltage dependent algebraic load model, with under load tap changers, switchable shunt capacitors, control in the transmission, subtransmission, and f \Olic Mum A'CL‘enr PCP-iii! T3330; "m p435 End 31 and higher voltage rating distribution system is the best model available at present to study voltage collapse, they will be used in this thesis with the following assumptions 1. The load dynamics are electrically far from the generator such that they are decoupled from the generator dynamics; 2. There is reasonably good voltage control at the load buses in the system so that the load variation with voltage is small; 3. If the magnitude of stable orbits associated with generator dynamics remains large enough and if the frequency of the associated oscillations remains small enough, then significant control actions of tap changers, switchable shunt capacitors and thermostatically controlled loads should take place in response to these oscillations. These assumptions indicate a constant power load model would be satisfactory to determine the equilibrium point set by the actions of tap changers, switchable shunt capacitor and thermostatically controlled loads, and that a perturbation load component that is voltage dependent and responsive to orbital dynamics need not exist using this assumption. Indeed, a perturbation load component model has been validated and used in recently developed simplified dynamic load models, that include both steady state and perturbation components. The perturbation load component model can be viewed as the Taylor series approximation of the system load model around equilibrium [31]. 2.4 Classes and Kinds of Bifurcations in a Type 2 Model A classification of stability problems that are observed and/or may occur on a Type 2 power system model is established in terms of kind and class: Kind - Kind (SN, Hopf, Load Flow bifurcation...) of local bifurcations that may 32 occur on the specific model type. Class - The subset of generators, controls, and/or load dynamics which exhibits the specific kind of bifurcation in this model type. The bifurcation phenomena in the nonlinear power system model (2.1) refers to characterizing the qualitative change of the fixed points of the differential-algebraic model in (2.2) for a smooth continuous change in parameter value it over a specified range. Although all system bifurcations usually result in a qualitative change in system behavior, they however have different effects on the dynamics of the system. Kinds of bifurcation simply refers to saddle node, Hopf etc., while a class is a further division of each bifurcation kind into the specific system dynamics experiencing the bifurcation. Kinds of bifurcations observed in a Type 2 power system model have been reported in [25] and include saddle-node, Hopf, load flow bifurcation and period doubling bifurcations. Loss of causality bifurcation may be split into a voltage, angle, or both voltage and angle bifurcation classes. Similar classes may be found for the other kinds of bifurcations. A class of bifurcation is necessary to (a) identify the subsystem dynamics and specific parameter u or parameters a. of a power system model, that cause a specific kind and class of bifurcation to occur in the full dynamic model; (b) to develop methods that can identify all potential bifurcations that exist in a particular model; and (c) to find the underlying continuous or discontinuous parametric, operafing schedule, or structural model control changes that produce such bifurcations. In order to present the kinds and classes observed on a Type 2 power system model, rewrite the differential-algebraic model (2.1) as x = f(x.0T. VT. u) 2.4(a) 3 3 t) = (}P(x,9T,VT,1.t) - 2.4(b) GQ(X, 67" VT’ ll) 2.4(C) C II where x = [6, 5",. y] : generator dynamic state variables: 6 is the angle at synchronous machine internal bus excluding swing bus, E", is the voltage behind transient reactance at generator internal bus and y presents the states of exciter dynamics of synchronous generator. VT , GT : voltage and angle at generator terminal, generator high side transformer and load buses, not including synchronous machine internal buses and the swing bus. u : parameter of interest that has an impact on the stability property of the equilibrium points of (2.4), such as real and reactive loads and generator real generation. Local stability of a steady state equilibrium point of interest for some u e R Y. = (ram). 91001). V1000) may be investigated using a linearized version of (2.4) evaluated at the equilibrium, obtained by solving the system of equations 0 = f(x.GT. V1.11) 2.5(a) 0 = GP(X, 9T, VT, ‘1) 2.5(b) O = 00(x, 9T, VT’ u) 2.5(C) The Jacobian matrix of (2.4) at the equilibrium point Y0 = (xo(u), 91 (u), VT (11)) may be obtained by linearizing the differential and algebraic equations and then combining all the linearized components as shown below ax Ax 0 = 100' A91 2.6(a) o AVT 34 where "er. 2:. 2:2 ax 39, av, DFO—(O) 30,, BGP 30,, ’(“)-_ax - ‘37 Te; 5v; aGQ Bag 809 L 8x BOT 3V1}? A1 31 C1 A3 33 C3 X = (1(9).9T(i~l). V7~(il)) F(X) = [f(X).Gp(X).GQ(X)] An equilibrium point X ‘0 = (x3, 9%,, Va.) of (2.4) is called a bifurcation point and u‘ a bifurcation value if the complete Jacobian matrix J(u") is singular, in which case existence or uniqueness of solutions is lost. Since the power system model as given in (2.4) is composed of both algebraic and differential equations, it is clear that an equivalent differential or equivalent algebraic equation model ought to be obtained to derive the kinds of bifurcation that may occur in the model, as will be shown in the following section. 2.4.1 Kinds of Bifurcation From the implicit function theorem [33], a unique solution to 2.4(b) and 2.4(c) exists locally if and only if the matrix J (u) ___ 32 c2 5 BGP/BBT aGP/BVT ‘ B3 C. aGQ/BBT aGQ/BVT go is nonsingular. If Jc(|.l.) is nonsingular for all x, 8T, VT and u in a neighborhood of )‘to then 2.4(b-c) may be solved uniquely for OT(x(|.1)),V7-(x(u)), which results in the equivalent differential equation model Ax = f(x. 91001)). V1061)» 35 which linearization gives the linear equivalent differential model x = rxru) . Ax 2.7(a) -1 Jx(ll)=Ar-[Bl c1][:2:2] [2] 2.70») 3 3 3 Similarly, a unique solution to 2.5(a) exists locally if and only if the matrix 41005 SE l? x 0 is nonsingular. If A101) is nonsingular for all x, 9T, VT and u in a neighborhood of if" then 2.5(a) may be solved uniquely for x(6T(u), VT(u)) . which results in the equivalent algebraic equation model 0 = g(x(97~(u). V100). err»). VT(u)) 2.8(a) or upon linearization gives the linear equivalent algebraic model 0 = Jy(u) - Ay 2.8(b) _ 32 C2 A1 -1 Jy(|.l) .. [33 c3] - [Ali/4.1 [3, cl] 2.8(c) The power system defined by (2.4) with equivalent differential equation model 2.7(a) is said to be asymptotically stable iff all the eigenvalues of the equivalent differential system Jacobian matrix J,(l1) have negative real parts. Hence a critical state is a state when one or a pair of eigenvalues become non-hyperbolic and have zero real part. Kinds of bifurcation on a Type 2 model are summarized in the following: (0 WW Occurs at u= u“ when Jc(u‘) is nonsingular and the system (2.7) has one real eigenvalue that crosses the imaginary axis of the complex plane and where certain transversality conditions hold [App C]. In this case the system trajectories are attracted by the center manifold, and the system behavior becomes monotonic and (ii) 36 may experience monotonic instability. i.e. collapse type of instability. II [E .E . Ill: Occurs at u= u‘ when 1001‘) is nonsingular and the equivalent differential Jacobian matrix 1.01.) has one pair of imaginary eigenvalues with zero real parts. Additional transversality conditions must also be satisfied [App C]. In the neighborhood of this point, periodic solutions always exist and the stability of the system is determined by the stability of the periodic solutions (limit cycles). A Hopf bifurcation can be stable (supercritical) or unstable (subcritical). Both supercritical and subcritical Hopf bifurcations have been observed on a single-machine-infinite-bus model [25]. (iii) WW) Occurs in the system equivalent algebraic model 2.8(a), at u: u‘ when Al(u") is nonsingular and the equivalent algebraic Jacobian Jy(u‘) in 28(c) has one simple zero eigenvalue. This bifurcation is called Load Flow (LF) bifurcation because it indicates a system at equilibrium is unable to satisfy the system power flow requirements imposed by the given generation and demand requirements.The system in 2.8(a) is called the Exact Load Flow Model since these equations represent the power balance equations at all buses in the system, while the ‘classical’ load flow model will be shown in Chapter 4 to be an approximation of this model. (M LasaQLCausallnr Occurs when Jc(u) becomes singular as it varies. Simulation cannot be performed on the system (2.4) since 2.4(b-c) have no solution or have multiple solutions 37 91",(X(ll)) Vri(x(l1)) i= 1,2... in the neighborhood of an equilibrium 550 = (x0, 9w VTo) . If there is no solution, the system trajectory terminates abruptly. If multiple solutions exist, random discontinuous trajectories can occur as the solution (en-(x00), VTi(x(u))) changes rapidly from time point to time point. Two types of loss of causality may occur. (a) Algebraicfiifirtcaficnflil Occurs when the submatrix of 1(a) . A a c 1' 1m”: 2 2 2 = f" A3 33 C3 JQL becomes row dependent at u: u‘. Hence it requires that the differential- algebraic bifurcation occurs in a subset of the equations. Note that algebraic bifurcation may suggest that there are no solutions of the form (x7, 9T, VT) or multiple solutions (113971,» VTr)’ i=1,2,.. (b) 5' l . I I lE°fi . :51: Occurs when the full system Jacobian matrix 1.01) is nonsingular with eigenvalues crossing from the right half plane to the left half plane or vice versa by going through infinity, rather than across the imaginary axis (jar-axis). At this point, one real eigenvalue of the system is at infinity and the rest of them are bounded which makes simulations using the differential-algebraic model not possible. (V) W Occurs in the case of multiple eigenvalues when a pair of complex conjugate 38 eigenvalues become two equal real eigenvalues. A node focus bifurcation (NF) may be stable or unstable depending on whether the real part of the eigenvalue pair is negative or positive, respectively. (vi) Wee: Occurs in the case of multiple eigenvalues when two pairs of complex conjugate eigenvalues become two equal complex eigenvalues. At (SN), (H) and (SI) bifurcation points, local system steady state stability around the equilibrium point changes. At 1:1 resonance point, the system may experience non- linear oscillatory behavior due to order 1 resonance. In this case, participation factors, eigenvectors and residues can be misleading because eigenvectors are not unique and they lose their physical meaning. 2.4.2 Classes of Bifurcation For each kind of bifurcation, the bifurcation subsystem method will attempt to identify subsystems of smaller and smaller dimension to determine the smallest bifurcation subsystem that produces and causes the bifurcation. The bifurcation subsystem implies the class of bifurcation of the particular kind which will be investigated in Chapters 4-6. Classes of bifurcation in a Type 2 power system model may be (a) the power system differential equations; (b) the power balance equations at load buses; (c) a generator’s inertial, flux and exciter dynamics; (d) a generator’s electrical dynamics; . (e) the control systems dynamics ( PSS, Govemor..) on a generator; (f) a combination of the above. 39 An investigation to find the smallest bifurcation subsystem that experiences and causes the bifurcation in the full system is initiated theoretically in Chapter 4 and continued to the end of the thesis. 2.4.2.1 Classes of Loss of Causality in a Type 2 Model Although the two kinds of loss of causality system bifurcations (algebraic bifurcation, and SI bifurcation) all produce truncation of trajectory, they have quite different effects on generator dynamics. Loss of causality may occur in real power-angle variables, in reactive power-voltage variables, or in both angle and voltage variables. (i) mislelrcsscmausaliut Occurs when [32 C2] = [Bap/89, aoP/avT] is row dependent, thus since [32] is singular, this loss of causality is due to real power-angle instability. (ii) W Occurs when [83 C3] = [609/89, aoQ/aVT] is row dependent, thus since the square matrix [C3] it singular, the loss of causality is due to reactive power-voltage instability. (iii) WWW Occurs when Jc(u) is row dependent but neither [82 C2] nor [B3 C3] is row dependent resulting in a voltage and angle loss of causality. 40 2.4.2.2 Classes of Static Bifurcation in a Type 2 Model Static bifurcation in a Type 2 model occurs when the matrix 1.01) is nonsingular and -l B C A W = [3‘ €11le ll 3 is singular as m changes. Saddle node bifurcation is a special but generic case of static instability which occurs when a single eigenvalue of 1,.(11) is zero and certain transversality conditions hold [App C]. Two classes of static bifurcation may occur on a Type 2 power system model: (i) S'El .. °'l| . Static bifurcation in generator inertial dynamics may be associated with the loss of transient [App. B] stability [36] such as the voltage instability that occurred in the Czechoslovakian system after tripping lines in a major interface [30]. This class of static bifurcation is likely to occur when the disturbances are very large and the transmission interfaces or boundaries are weak, so that loss of stability occurs soon after the contingency before the field current level/duration limits are reached. (ii) W This bifurcation has been called classic voltage instability and is due to insufficient reactive supply in the EHV system, an example of which is the blackout that occurred in the French system [37]. This class of stability has been analyzed with the use of a sensitivity matrix between reactive generation and voltage at generator internal buses 5965' The matrix SQGE is diagonally dominant for inductive networks (networks that absorb reactive power), and not diagonally dominant for capacitive 41 network (network with long high voltage transmission lines). Classic voltage instability has been described as occurring when generators reach their field current limits and the field current limit controllers disable the exciters. An important implication to this result is that studying this saddle node bifurcation must include the field current limit controllers action in the model. 2.4.2.3 Classes of Dynamic Bifurcation in a Type 2 Model Three classes of dynamic bifurcation may exist and have been observed in a Type 2 power system model (i) (ii) IIE' "lffllf°l ': Hopf bifurcation in generator inertial, flux decay and exciter dynamics and is associated with low frequency inter-area oscillations in a power system. This mode of oscillation is usually referred to as a swing mode and appears as stable limit cycles (supercritical) [15]. Wanna: Hopf bifurcation in the flux decay/exciter control dynamics. This mode of oscillation is usually referred to as voltage swing modes and appear as unstable limit cycles (subcritical) in literature [15]. (iii) Ii E' . °lffl l l . : Hopf bifurcation in generator inertial/flux decay dynamics and occurs only when the exciter in completely disabled. This mode of oscillation is called the voltage collapse mode and has been shown to occur in the single machine infinite bus model [15]. III POWER SYSTEM STABILITY ASSESSMENT A TAXONOMICAL APPROACH 3.1 Introduction Power system stability can be broadly defined as that property of a power system which enables it to remain in a state of operating equilibrium under normal operating conditions. and to regain, in a finite time, an acceptable state of equilibrium after being subject to a disturbance [29]. Stability analysis of post contingency equilibrium point or limit cycle using dynamical system theory for smooth nonlinear systems, although necessary, is not generally sufficient for assessing stability of equilibria or limit cycles of the power system dynamical model. Stability analysis for smooth nonlinear systems must also be performed for points in parameter space where parameter change is continuous. Assuring power system security based on describing local stability of an equilibrium, determining local connected regions in state and parameter space where local stability is retained (feasibility region) and region of stability for trajectories when parameters are fixed, has been investigated in [26]. This work is being extended [27] to include effects of hard limits and equipment outages which discontinuously modify the model order, the structure, equilibria and dynamic behaviors 42 43 A realistic power system model is characterized by nonsmooth elements, the most important of which are the hard limits encountered during system operation. In power system terminology, “saturation” is usually reserved for denoting smooth effects associated with saturation phenomena such as field saturation; and “hard limits” are usually associated with nonsmooth effects such as actuation limits, tap changer limits, relay effects, etc. From a conceptual point of view, there exists three types of hard limits: (a) windup limits or actuation limits [27]; (b) nonwindup limits or state limits [27]; and (c) relay type limits or switching limits [27]. Depending on how a system with saturation limits comes back to the limit, windup and non-windup limits are distinguished. Whereas windup limits can come off the constraint at any time, this is not possible for nonwindup limits before the windup variable has been released completely. Hard limits are generally associated with actions of operator control; actions of protective devices such as relays or field current limit controllers; actions of under load tap changers and actions of switchable shunt capacitors. An example that qualitatively portrays the effect of consecutive discontinuities on a state trajectory is shown in Figure 3.1 as well as Figure 3.2 [40] and Figure 3.3 [41] for simulations on actual detailed power system models. The standard bifurcation theory for smooth systems [26], being based on the functions f and g in Equation (2.1) being smooth, is not applicable to assess the stability of the system experiencing such discontinuous transitions. Since voltage collapse time scenarios are characterized by all of the discontinuities described above, it is necessary to extend the stability analysis for models with discontinuities (hard limits and contingencies). The stability analysis of a power system is complicated by the lack of smoothness 44 needed to characterize not only the bifurcations existing on smooth systems, but also because new bifurcations that are directly connected with the hard limits must be described.'Such bifurcations include hard limit induced static bifurcations and hard limit induced dynamic bifurcations [27,42-43]. Although the purpose of this thesis is not to develop an extension of stability theory for hard limits, a brief description of the current development on this subject is necessary for consistency and a better understanding of the theoretical framework under which the bifurcation subsystem method is being developed. In this chapter, a rather brief description of the recent development of the taxonomy of the large differential algebraic power system model is given. Zaborsky’s taxonomical work of 1991 [26] which summarizes a method for applying dynamical system theory to a smooth power system model is reviewed in Section 3.2 and a discussion of the status of the theoretical development on the extensions of this taxonomy to include models with hard limits discontinuities [27] is briefly reviewed in Section 3.3. A parametric episodal approach to security assessment is then proposed that a. analyzes the stability region, feasibility region and viability region for the continuous analytical model in intervals (if, rg+ 1) between instants t,- and ti“ where discontinuities brought on by equipment outages or hard limit induced model changes occur; and uses this analysis to decide stability, feasibility and viability of the entire trajectory; b. enumerates all equipment outages and operating changes for which the trajectories and attracting equilibria and limit cycles of the system are expected to be stable and viable and then establishes for which contingencies and operating changes the trajectories are indeed unstable or nonviable. Viability for equilibria, limit cycles and 45 transients is defined to be those equilibria, limit cycles and transients which (i) do not cause equipment damage; (ii) do not cause relaying actions for protecting equipment of the system; and (iii) have voltage, current, frequency and power variations that are within acceptable design and operating ranges. )’(1i+l') r I ‘i ti+1 (time) Figure 3.1 2-D example portraying the effects of the ith discontinuous control action. (yi) is the ith quasi equilibria 3.2 A Taxonomy of Stability Assessment in Power System A taxonomy theory refers to a general and structural theoretical analysis of the dynamic behavior of very large nonlinear systems. Assessing local and transient stability [App B] of a power system as well as maintaining its security for parameter changes in an operating range is certainly a taxonomical task. Security has been defined in Zaborszky [26] in parameter and state space as operating conditions which belong to the feasibility region, stability region (region of attraction) and viability region. The state space includes 46 both the dynamic state vector x and algebraic state vector y, whereas the parameter space includes structural parameters such as transmission network capacitances and topology, and operational parameters that could include a. the level and variation of real and reactive load as a function of voltage and frequency at each bus; b. the real generation at every generator and its change with frequency, load level, and area control error as well as its upper and lower limits; c. the voltage setpoints on generators exciters, or voltage setpoint upper and lower voltage limits on switchable shunt capacitors and the upper limit on their susceptance, and voltage limits and tap setting limits on under load tap changing transformers; d. the transfer or wheeling actions taken. There are numerous parameters above and beyond those itemized above but these are certainly important parameters in deciding when voltage stability is retained or lost [11]. Assessing stability of the large power system for a specific equilibria or limit cycle requires the assessment of four regions: region of stability, feasibility region, viability region and security region, which definitions are based on Zaborszky [26], the Zaborszky’s discussion [47] of a paper by Schlueter and Benkilani [25], and their reply [48]. The definitions are not new but generalize the definitions found in [26] as suggested by the discussion and reply. 3.2.1 The region in state space where every trajectory originating inside, will finally con- verge to the specific attractor of a particular type (equilibrium point or closed orbit) for a 47 specific set of operating parameters p = p0. The region of attraction in state space is dif- ferent for each specific point pa in parameter space that belongs to a feasibility region of that specific attractor of a particular type. 322 E .1]. E . [5 'E E .H. 1..: I The region in parameter and state space, within which a specific attractor of a partic- ular type (operating equilibrium point or stable closed orbit) can be shifted smoothly from any point in the region to any other point by continuous parameter change while retaining local stability. Structuralstability holds [34,50] within the feasibility region for each spe- cific attractor of a particular type. 323 W The region in parameter and state space belonging to a feasibility region of a specific attractor of a particular type where the values of the state variables, i.e. voltages, current, frequency, active or reactive power for the specific attractor of a particular type remain within prescribed limits. The bounds can be based on physical device limitations on parameters or states; safe operating limits on parameters or states; hard limits based on device limitations or relay actions for protection of the device or system. 32.4 W A trajectory that is stable and belongs to a region of attraction of a specific simple attractor of a particular type that belongs to the feasibility region for that specific attractor of a particular type is viable if the current, voltage, frequency or other variable trajectory does not cause equipment damage, does not exceed design limits, and does not trigger relay actions that protect equipment or system operation. 48 3.3 Development of a Taxonomy for Power System Models with Hard Limits The taxonomy developed in [26,27] is certainly the first rigorous effort to describe stability in a differential-algebraic power system model, and yet the taxonomy is still under development. One major area of development is the inclusion of hard limit-induced discontinuities, which Occur during operation when 1. the generator field current reaches its magnitude and duration limits so the field cur- rent limit controller either causes a relay to disable the generator exciter or reduces the exciter voltage setpoint and thus field current down to a sustainable level that will not produce thermal damage to the generator rotor; the controlled bus voltage remains outside the' prescribed upper or lower voltage lim- its, and under load tap changers discontinuously and abruptly change the turns ratio on a transformer to bring the voltage within limits. The action is possible as long as the tap settings remain within tap position limits. If under load tap changers reach upper or lower tap position limits, a second hard limit induced discontinuity is encountered; the voltage at the controlled bus exceeds lower or upper voltage limits, a switchable shunt capacitor switches one or more banks of capacitors in or out to bring the con- trolled voltage within limits. Voltage collapse in the distribution system can be prevented by actions of under load tap changers and switchable shunt capacitors as long as the under load tap changers remain within the tap setting limits and switchable shunt capacitors remain within the capacitive susceptance limits. However when a sufficient number of generators reach field 49 current limits, tap changers reach tap position limits. and switchable shunt capacitors reach capacitive susceptance limits, voltage collapse will occur due to loss of voltage con- trollability [11]. The discontinuities associated with generator field current limit controllers, under load tap changer controllers, and switchable shunt capacitor controllers can be modeled using actuation (windup); state (nonwindup); and relay hard limits. Due to the discontinu- ous nature of the actuation, state and relay hard limit models [27], the taxonomy described for smooth dynamical systems can not be used to describe the feasibility region and stabil- ity region for voltage collapse which develop when generator field current limit control- lers, under load tap changer controllers and switchable shunt capacitor controllers hit hard limits in response to equipment outages or continuous operating changes. Zaborszky [27] has started the long process of extending the taxonomy so that it can include actuation, state and relay limits. The status of taxonomy development and its structure is now described. (1) The model used [27] is a differential equation model rather than the differential- algebraic model used in the taxonomy without hard limits [26]. When hard limits are ignored in the model, the boundary of the feasibility region has segments due to Hopf bifurcation, saddle node bifurcation and singularity induced bifurcation, whereas singularity induced bifurcation can not be part of a feasibility region bound- ary in the taxonomy being developed for a differential equation model with hard lim- its. In the taxonomy without hard limits [27], the boundary of the region of stability allows (i) stable manifolds of order 1 equilibria; (2) (3) 50 (ii) stable manifold of periodic orbits; (iii) stable manifolds of pseudo saddles, semi saddles and bad anchors; (iv) singular boundary pieces; among others where the stable manifolds of pseudo saddles, semi saddles, bad sad- dles and singular pieces could not belong to the stability boundary for a taxonomy being developed for a differential equation model with hard limits [27]; the feasibility region for a differential equation model with actuation and state hard limits have been derived but not for relay limits. The feasibility boundary for actua- tion (state) limits is composed of segments due to saddle node bifurcation, Hopf bifurcation, actuation (state) induced dynamic bifurcation, and actuation (state) induced static bifurcation. The actuation induced static bifurcation and state induced static bifurcation behaviors near the feasibility boundary are similar to saddle node bifurcation except that the transversality and bifurcation conditions are somewhat different. The actuation induced dynamic bifurcation behavior is similar to Hopf bifurcation except that the transversality conditions at the bifurcation point are dif- ferent; When windup limits are modeled in the analysis, equilibria and their stability prop- erties change and the system transient trajectories are characterized by nonsmooth- ness (“kinks”). The unstable equilibria which anchor the stability boundary may correspond to points where hard limits are encountered, hence results from these equilibria differ from those of the smooth systems. Although the boundary of the stability region is composed of stable manifolds of unstable equilibrium and periodic orbits, these manifolds have comers or kinks (actuation) or have different dimension if state or relay limit are reached. The analysis of the region of stability and charac- terization of its boundary is only complete for actuation limits. This task is far more challenging for state and relay hard limits since the state space for state and relay limits consists of several smooth systems of varying dimensions. The most important conclusion [27] is that there is indeed a stability region bound- ary when hard limits are included in the model and where one or more hard limits are reached either on the trajectories on the stability boundaries or on trajectories that are ini- tiated at some point within the region of stability and converge to some equilibrium. Another conclusion is that within the epoch (2;, (5+ 1 ) , the trajectory is smooth and the fea- sibility and stability regions satisfy all the properties of a model where hard limits are ignored. A final conclusion is that the nonsmoothness of the stable manifolds means that the construction of Lyapunov functions or energy functions is more difficult and requires different techniques than the typical u.e.p based methods used in models without hard lim- its. 3.4 Stability Security Assessment for Model with Hard Limits The desire in developing a taxonomy for studying stability in a power system model is to be able to assess retention or loss of stability without time simulation of the system, or quite possibly with almost no simulation of the contingency or operating change. The stability region, feasibility region, and viability region and their boundaries are known if such a taxonomy is complete. Even if the taxonomy were known, its success of avoiding time simulations appears questionable, since one is unable to predict apriori the post con- tingency or operating change steady state equilibrium point or limit cycle without time 52 simulation. This is because (1) the set of hard limits and the sequence through time at which they are encountered is impossible to predict apriori and (2) knowing the set of hard limits and their sequence over time is necessary to determines the final post contingency equilibrium point or limit cycle. The fact that the set and time sequence of hard limits encountered not only determine the equilibria or limit cycle trajectories converge to but whether they are stable or unstable, is documented in [41]. The reasons that it is impossi- ble to predict the equilibria or limit cycles and their stability are (i) that the hard limit encountered in a power system model are never enabled instantaneously when the limit is encountered, (ii) the delay for each hard limit is different and can depend on the magni- tude of the violation and duration of the hard limit, and (iii) the delay may be large com- pared to the transient response of the dynamics, (iv) the sequence of future hard limit control actions are dependent on previous ones, (v) the system is highly nonlinear so dif- ferent sequence of hard limit control actions give different ultimate results. An approach is now proposed that studies stability of the equilibrium or limit cycle as well as transient stability epoch by epoch where an epoch is a period of time I: < t < t; + 1 (see Figure 3.1) where no hard limit transitions occur. This method would utilize the tax- onomy to study stability of equilibria, limit cycles and trajectories within epochs, rather than attempting to study stability of the entire trajectory as well as the steady state equilib- rium or limit cycle without simulation. The justification for this approach is that the taxon- omy being developed with hard limits [27] indicates that the taxonomy developed for models without hard limits [27] applies to time intervals (12.1}, 1) as well as for subinter- vals in state space where no hard limit transition surfaces are encountered. An Epoch Based Trajectory Stability Assessment Method is discussed in Section 3.4.1. 53 A second change in stability assessment methodology to be discussed is to consider an episodal trajectory stability assessment rather than attempting to use a taxonomy which might determine stability without time simulations for all trajectories resulting from (a) disturbance initiated changes in initial state but no change in model; (b) continuous slow operating changes in parameters of that model and (0) power system models with no delay hard limit so that the set and sequence of hard limits encountered can be known. Utilities currently utilize an episodal stability assessment rather than a taxonomy based method in dynamic security assessment and contingency selection for fault contingencies [52]. The motivation for this taxonomy is based upon the earlier work on fault contingencies. The use of an episodal rather than a taxonomic method is due to the fact that most realistic fault contingencies are followed by a number of discontinuities which could not be pre- dicted to occur without time simulation. Thus taxonomy based methods and episodal methods are competing in an application where taxonomy based methods held dominance for years due to their promise in providing accurate stability assessment without the need for time simulation. An episodal trajectory stability assessment is thus proposed for the equipment outages and operating changes that cause voltage instability. Utilities currently specify the contingencies and the range of operating condition changes they expect to sur- vive without loss of voltage stability. If one can quickly assess stability of trajectories pro- duced by each of the set of contingencies for the ranges of operating conditions specified epoch by epoch to find those that are stable and those that are unstable, then a contingency selection or stability assessment procedure, in the language of power system engineering, could to be produced. The set of trajectories that are viable and transient stable [App B] to attracting sets (equilibria and limit cycles) that are feasible and viable for every epoch are 54 said to belong to a Bounded Trajectory Feasibility Set. A discussion of the contingency selection and Bounded Stable Trajectory Feasibility Methodology is given in Section 3.4.2. 3.4.1 Epoch Based Trajectory Stability Assessment The Epoch Based Trajectory Stability Assessment is proposed because each hard limit induced model change may produce (a) a model xi = fi(xiryi»l*i) ‘:<’<‘i+r (3.1) O = 810pr Hi) where an equilibrium (mi) may or may not exist; (b) a stable or unstable equilibrium point; (c) a stable limit cycle encircling the unstable equilibrium point; ((1) a region of attraction for that stable equilibrium point or limit cycle that is suffi- ciently large to include (x(z;), y(r;)) . All four questions (a-d) must and will be assessed for each time epoch (If, t' i+ r) to be able to diagnose when and why the trajectory does not converge to a final attracting set as time progresses. The Epoch Based Trajectory Stability Assessment is based on a trajec- tory simulation method implemented by Van Cutsem in 1992 [28] that does not require time simulation in (1321;, l) , but only computes the quasi equilibria (ivy!) for i = 1,...,N since it has been found that (a) the transition to the quasi equilibria (5:9,) is stable because the region of attraction is generally large; (b) that xi(t‘-’+ l) = ii and yi(t‘?+ 1) = 5:1. since the time T= lf+ r —ti+ between hard limit 55 model changes is large enough to cause the trajectory to reach steady state at each hard limit model change; (c) The lack of a stable equilibrium without existence of a stable limit cycle may be observed as a diverging set of quasi equilibria. If one has an unstable sequence of quasi equilibria (fl-,5") but a stable encircling limit cycle around each quasi equilibria, this rapid simulation method [28] would not indicate retention or loss of stability. Furthermore it has been found that the short term dynamics and long term dynamics interact and thus such a rapid simulation method is not always adequate [41]. The Epoch Based Stability Assessment Method utilizes the rapid simula- tion method to evaluate quasi equilibria (ivy!) but takes into account not only stability and feasibility of quasi equilibrium (in?) but also stability of trajectory during each epoch. In order to illustrate the argument that the rapid simulation method is needed to sim— ulate the voltage stability trajectories with a large number of discontinuities but is inade- quate to assess the interactions of short and long term dynamics that cause instability without Epoch Based Trajectory Stability Assessment, we bring out two examples from recent literature. Figure 3.2 shows the time simulation of the active power load of a large utility system after a severe fault where the over-excitation limiter controller (OEL) dis- abled the exciter on one machine but all other machine protection devices are enabled [41]. The disturbance was an EHV line fault and line trip which resulted in tripping a large generator unit off-line. The model includes excitation systems, turbine/govemor systems, PSSs, generator protection functions, AGC...etc. The transition dynamics which occur due to the action of voltage protective devices (LTCs, switchable shunt caps.) clearly manifest a nonsmooth behavior. A number of discontinuities occurred before the system was driven 56 to voltage collapse and thus a rapid time simulation method [28] is needed. This example would not need to simulate the trajectories in epochs but only needs to determine quasi equilibria provided by the Epoch Based Trajectory Stability Assessment because the quasi equilibria capture the response of the system and no short term dynamic behavior is observed. Feasibility and viability of equilibria could be utilized to assess why the system experienced voltage collapse. Figure 3.3 is intended to demonstrate interaction of short term and long term dynamic behavior via simulation of the active power flow between two regions following tripping of a line and a generator in a Swedish test system [41]. We notice that due to the severe fault, the dynamic behavior was initiated by large oscillations, which after the acti- vation of the OEL and the tap changers reduced to smaller limit cycles. As the voltage keeps dropping, the system was finally driven to voltage collapse. Three crucial phenom- ena are demonstrated by this figure: (1) Viability of limit cycles is necessary in each epoch; (2) quasi stable/unstable limit cycles rather than quasi equilibria may exist; and (3) both long term and short term dynamics must be considered to assess stability and security of the entire trajectory. Viability, feasibility and region of stability of trajectory segments and the feasibility and viability of equilibria and limit cycles would be needed to assess why voltage collapse occurred in this system. In order to address the stability of the full trajectory (x(t), y(t)) as well as the interac- tion of short term and long term dynamics, time simulation or taxonomy based stability assessment within epochs I: < z < :;+ r is necessary. The Epoch Based Trajectory Stability Assessment that avoids complete time simulation within 031;”) that utilizes a taxon- omy stability assessment in each epoch to computes 57 (1) X(t;) = J‘c,._I and y(:;) = 5L1 by solving for (ii-l’yi—l) from fi-1(xi-l'yi—1)= 0 32(8) 85-1(xi-1,)’1_1) = 0 3.2(b) (2) y(t;‘) by assuming x(t;') = x(t;) = 31.1 using g;(i.-_1.y(t.-*)) = 0 (3.3) This assumption is based on the fact that whereas the algebraic state vector y(t) may change instantly, this is not possible for the dynamic state vector x(t). Infinite power is needed to cause an instantaneous change in x( t). The Epoch Based Trajectory Stability Assessment then checks (i) if solutions (ii—1,534) and (xi_1,y(ti+)) exist; (ii) if (2,; 1, 5L1) is stable and viable; (iii) if (x(r,.+_,), y(z;'_ 1)) lies within the region of attraction of (ii-1’ 5L1) if step (ii) is true; (iv) establishes if there are attracting sets (closed orbits) surrounding the unstable equi- librium and whether they are stable if step (ii) is not true; (v) establishes whether these attracting sets are stable and viable if step (iv) is true; (vi) establishes whether (x(t;‘_l), y((‘-+_1)) lies in the region of attraction of the attracting set if steps (iv) and (v) are true. Stability and viability are assumed to be retained if the answer to step (i) and either steps (ii) and (iii) or steps (vi), (v), and (vi) are true. Such a procedure would eliminate the need to perform time simulation and would only require solving y(t;') using equation (3.3) and checking (i-vi) A slight expansion of the above (i-v) rules for testing trajectories via trajectory seg- 58 ments in epochs ($117+ 1) for the diagnostics of knowing feasibility, viability and stabil- ity is given below (a) (b) (C) (d) (8) Every post contingency steady state operating point or stable limit cycle, and every quasi equilibrium point (ii, 5".) or limit cycle which occurs at the PM discontinuous control action, must exist, be stable and viable and thus reside in the (non-empty) intersection of its viability region and its feasibility region. (feasibility, viability) Assuming (x(t‘-+), y(t;)) accurately reflects the unmodeled transition from (x(t;), y(t;)) to (my), y(t;‘)) for each i, the stability for the transition from (x(t;), y(t;)) to the stable steady state or quasi equilibrium point (it-,1.) or its stable encircling limit cycle can be decided by checking whether (x(ti+), y(r;)) lies within its region of attraction for each i = I , 2 N. Shrinking unstable limit cycles around the stable steady state or quasi equilibrium point constrain the size of its region of attraction. A lower limit on the size of the unstable limit cycles should be prescribed. (viability of unstable limit cycles) Every trajectory initiating in the region of attraction of each steady state or quasi sta- ble equilibria does not stall or trip motors, does not cause equipment damage, relay action to protect equipment , and operator emergency action to trip generators off the system. (viability of transient trajectory) Stable limit cycles around unstable steady state or quasi equilibrium points that are viable must also satisfy the additional condition to be small enough to prevent dam- age of equipment. An upper limit on the size of the limit cycles is prescribed so that (e) equipment damage could not occur or (b) so that the region of attraction could not disappear (as the limit cycles reach the boundary of the region of attraction). 59 (viability of stable limit cycles) Performing steps (ii - vi) or (a-e) is computationally difficult if not impossible at present due to the very large dimension of models. Research is underway to make possible rapid evaluation of steps (ii - vi). This thesis on development of a bifurcation subsystem method would allow one to quickly assess whether the equilibrium is stable or unstable and whether the attracting set (.5 one ext-as) is stable or not with very little computation. The bifurcation subsystem method identifies the small dimension subsystem that experi- ences the same bifurcation as the entire system and for the same value of the bifurcation parameter. The natural sparse structure, weak coupling and geometric structure of the model cause bifurcation subsystems to exist. Furthermore, the bifurcation subsystem method would indicate not only which parameters of a linearized model cause a specific bifurcation to occur, but also what operating changes cause that bifurcation to occur. The bifurcation subsystem method may also identify the subsystem dynamics that determine the stable manifold of a particular unstable equilibrium that belongs to the boundary of the region of attraction. Although investigation of how the dynamics of the bifurcation subsystem defines a stable manifold of an unstable equilibrium, is beyond the scope of this thesis, it is an important subject for the future. Bifurcation subsystems appear to be omnipresent in a power system model. The bifurcation subsystem can not generally be identified using eigenvectors or participation factors because the eigenvectors indicate the states affected by the bifurcation eigenvalue and not the subsystem that specifically experiences the same bifurcation as that experi- enced in the full system. The bifurcation subsystem of a differential-algebraic model may be differential, algebraic, or differential-algebraic depending on the kind and class of 6(1) bifurcation experienced. Each bifurcation subsystem needs to be checked for structural stability and retention of stability at each quasi equilibrium (it-.53) to assume stability dur- ing each epoch. If such assessments are rapid, the epoch based trajectory stability assess- ment is a natural procedure for performing a contingency selection and stability assessment that does not require time simulation in each epoch. 60000 $0000. cocoa soooo 20000 Figure 3.2 61 time (seconds) System active power load (MW) of a system [40] after a severe disturbance: EHV line Fault & line trip which resulted in a large generator ofi-line 200.. 100.. 20 I I l 40. 60 80 lm time (seconds) Figure 3.3 Active power flow between Central and Southwest regions of the Swedish test system [41] after the tripping ofa line & generator in the North 62 3.4.2 Bounded Stable Trajectory Feasibility Method A Bounded Stable Trajectory Feasibility Method is proposed because the trajectories initiated by contingencies or operating changes which experience voltage collapse or other types of instability are marked by a number of discontinuities such as the contingency itself or a hard limit-induced discontinuity such as described in Section 3.3. Since the number of such discontinuous control actions is typically large and quite different, an episodal classification based on trajectory segments and thus trajectories and their attracting sets being stable appeared to be logical and consistent with electric utility practice. Contingency selection methods have been widely used to a. determine contingencies where the currents on every branch (element) of the Hansnussion network lie within thermal limits and those contingencies where the currents on one or more branches exceed thermal limits. b. determine contingencies where voltages at every bus (node) lie within a secure range ( 0.95puSViS 1.05pu) and those trajectories where one or more bus voltages lie outside the range; .6 determine contingencies where voltage collapse does not occur in any subregion and determine the contingencies that cause voltage collapse in some identified subregion. Each of these subregions can be proven to be a bifurcation subsystem where the saddle node bifurcation makes it impossible to obtain solutions or may cause multiple solutions to exist. All of these contingency selection and security assessment procedures are based on load flow solution which attempts to find ilgn~(ii, 5".) where N is the number of discontinuities experienced along the trajectory. A Bounded Trajectory Feasibility Method 63 appears far more reasonable because it uses the Epoch Based Trajectory Stability Assessment that actually simulates (i) the actual sequence of quasi equilibrium points (3,, 53-) that would occur after each discontinuity; (ii) the transitions from ii = x(t;); )7, = y(t;) to x(t;') = x(r;); y(t;') = y(tg); also evaluates feasibility and viability of each equilibrium and limit cycle, and assess stability of trajectories (x(t), y(t)) over epochs (131‘; l) . The Bounded Trajectory Feasibility Method would utilize an Epoch Based Trajectory Stability Assessment rather that a load flow because the load flow often does not calculate a steady state equilibrium ileN(ri, 5".) that actually occurs since it does not determine the correct set or sequence of discontinuities that occur, and thus the correct equilibrium as explained earlier. Therefore the conclusions concerning its stability, viability and whether the initial conditions after a contingency lie within its region of stability are irrelevant Even if load flow calculated the correct steady state equilibrium [lgnNup ii) the hard limit based taxonomy would be difficult to apply to assess stability and viability of trajectories within each epoch. The set of equipment outage and operating condition combinations evaluated by the Bounded Trajectory Feasibility Method is clearly specified by a utility’s planning criterion. For each such combination the Bounded Trajectory Feasibility Method would determine stability and viability based on (a) existence and uniqueness of a solution or equilibrium point (39.5!) for each trajectory 1 segment; (b) stability of the equilibrium point or stable closed orbit for each trajectory segment; (c) asymptotic stability of the trajectory segment towards a stable closed orbit or equilibrium for the epoch (tf, I; + r ) ; (d) viability of the equilibrium, stable closed orbit, and trajectory segment transient for each trajectory segment; (e) classification of the trajectory as being or not being stable, and bounded trajectory feasible based on (a-d) for each segment. The Bounded Trajectory Feasibility Method would then attempts to assess why the bifurcation subsystem was affected by the equipment outage and operating change combination, why it was unstable or nonviable, and what could be done to prevent the problem, based on information that can be determined once the bifurcation. subsystem is known. This diagnostic step of the Bounded Trajectory Feasibility Method is similar to the diagnostic that can be performed for the voltage stability security assessment procedure in [11]. IV Bifurcations Subsystems in a Power System Differential-Algebraic Model 4.0 Introduction The focus of this chapter is to initiate investigation on the bifurcation phenomena in a differential-algebraic power system model for the purpose of identifying the various bifurcation subsystems for each kind of bifurcation (saddle node, Hopf and load flow bifurcation) observed on a Type 2 power system model. The bifurcation subsystem method will attempt to identify subsystems of very small dimension, called bifurcation subsystems, that are associated with each kind and class of bifurcation. Diagnostic classification of power system instability problems that is based on bifurcation kinds and classes, as well as model types used for bifurcation studies has been presented in Chapter 2. A Type 2 power system model was selected for further study in this thesis. The kinds of bifurcations that have been observed on this model include saddle node, Hopf, load flow, and each exhibits a different qualitative change in dynamical or algebraic network behavior encountered at the point of bifurcation. Classes of bifurcations are associated with the specific subsystem dynamics that experience a particular kind of bifurcation for a given model type. Possible bifurcation subsystems of a Type 2 power 65 66 system model are (a) the power system differential equations, (b) the power balance equations at load buses, (c) a generator’s inertial, flux and exciter dynamics, (d) a generator’s electrical dynamics, (e) the control systems dynamics ( PSS, Govemor..) on a generator, (f) a combination of the above One known class of Hopf bifurcation occurs in generator flux decay and exciter dynamics and another occurs in generator flux decay and inertial dynamics. Two other classes of saddle node bifurcation are also known to exist; one in generator flux decay dynamics and the other in generator inertial dynamics. Similiar classes may be found for other kinds of bifurcations. Since a power system trajectory after a disturbance is generally characterized by a sequence of hard limit induced discontinuities, during each epoch (23.15”) , the bifurcation subsystem method will be applied to a differential-algebraic Type 2 power system model to (i) (ii) determine the subsystems that experience the same bifurcation as that observed in the full power system model during each time epoch (93:5,, 1) . These bifurcations capture the effects of the previous discontinuous or continuous parameter changes that had occurred in the past time epochs (13‘, (5+ 1) for k = 1.....i-1. The subsystem that produces the bifurcation in the subsystem and in the full system may at times be larger than the subsystem experiencing the bifurcation of the full system. find the cause and the remedial action for preventing the bifurcation without actually finding the center manifold dynamics. The subsystem containing all the causes 67 including model parameters, controls, and model discontinuities may extend outside the subsystem that experiences or produces the bifurcation. Even though the research ultimate objective is to find the subsystem experiencing, producing, and experiencing bifurcation, the present work is solely aimed at finding the subsystem experiencing bifurcation as the full system experiences bifurcation. The diagnostics on the bifurcation subsystem may be more easily performed because retaining how specific continuous model parameter changes, how continuous operating, or scheduling changes affect the center manifold dynamics is not easily retained on very large power system models. The ability to analyze how the center manifold change along the whole trajectory when discontinuities are encountered has been investigated by Zaborszky [27] and is not an easy task. The traditional methods for identifying subsystems are generally based on apriori knowledge of (a) the system or its control design, (b) time scale‘information on real eigenvalues or imaginary components of complex eigenvalues, (c) eigenvectors and participation factor information from eigenvectors, or (d) decoupling within the state space model. While such methods are often useful in model reduction that preserves time or modal behavior at one operating condition or developing and redesigning particular excitation systems, system stabilizers or FACTS controls, they are not always effective in identifying the smallest subsystem or subsystems that manifests the same bifurcation in the full system. These traditional methods (a-d) for determining subsystems are not effective because obtaining a bifurcation subsystem model which experiences the bifurcation in the full model has a different objective than preserving modal or coherent behavior, designing controls or developing operation schedules for a system. The objective 68 of a bifurcation subsystem is to capture the nonlinear behavior that produces and causes a bifurcation in a specific eigenvalue or a particular model or operating condition changes. Time scales are often not effective in defining a bifurcation subsystems for a specific bifurcating eigenvalue in a power system model, because the frequencies of all inter-area oscillations lie in a narrow band and because the frequencies of all local generator oscillations also lie in a narrow band. A singular perturbation method would preserve all the modes in a band and not just the bifurcating modes for a particular bifurcation parameter. Subsystems used for control design must be satisfactory nonlinear model representations of the full system and the associated linearized model must be a satisfactory representation for several eigenvalues of the specific device or process being controlled and for a range of operating conditions. Despite the fact that eigenvector methods indicate the system state variables most severely affected by each Hopf bifurcation or the controls that are controlling that mode, this information does not always describe the location or cause of the bifurcation as noted in the power system example presented in Chapter 1. In that case, the modal information (participation factor) showed a rather large region near the city of Vancouver was the most severely impacted, but did not capture the fact that the disabling of voltage control on generators on Vancouver Island could actually cause the voltage collapse. Decoupling of subsystem dynamics by inspection of the structure of the system Jacobian matrix could also be used for model reduction and for identification of system states with high participation in a particular bifurcation. 69 4.0.1 Bifurcation Subsystems for Epochs The great majority of the literature on bifurcation in power systems assumes that the differential-algebraic power system model is continuous and differentiable. The discussion in Chapter 3 shows that virtually all bifurcations occur after (1) generators reach field current limits, (2) switchable shunt capacitors act, (3) tap changers reach upper or lower tap position limits, and (4) switchable shunt capacitors reach upper susceptance limits. All of these actions are discontinuous and violate the smoothness assumptions needed to apply bifurcation theory. Thus an extended bifurcation theory is being developed in [27] to properly describe the bifurcation that occurs in power system models that experience actuation, state, and relay hard limits [27]. Since inclusion of these discontinuities is absolutely essential for system stability analysis of power systems, the framework for the bifurcation subsystem method should accommodate and include the effects of such discontinuities. If the effects of discontinuities were neglected, then the bifurcation subsystem method would have limited and possibly little value, since there are relatively few if any bifurcations that occur in an actual power system without one and likely several discontinuities. Therefore, in this chapter we adopt the Epoch Based Trajectory Stability Assessment method, where system stability analysis is performed epoch by epoch where an epoch is a period of time (1:45”) free of hard limit discontinuous transitions. In this method presented in Chapter 3, stability analysis for smooth systems may be performed within epochs (1;, 117+ l) for i =1,... N, where 6+ 1 — z: is very long compared to the time constants of the linearized dynamics, rather than attempting to study stability of the entire trajectory. In each time epoch 031;”), bifurcation subsystems need to be determined and checked 70 for retention and structural stability at each quasi equilibrium (Sq-.5"). Since the progress of discontinuities during system operation is accompanied by loss of control and loss of subsystem coupling, bifurcations may be of different kinds, and bifurcation subsystems may be different for each epoch. Consequently, the bifurcation subsystem method makes the Epoch Based Trajectory Stability Assessment Method a rapid and natural procedure for stability assessment and contingency selection in a power system. 4.0.2 Existence of Bifurcation Subsystems Stability analysis and testing for bifurcations in a power system had always put special focus on its subsystem structure and on finding the critical subsystem dynamics or subsystems which exhibit the bifurcation. For example, it is quite clear that the generator flux decay dynamics and the generator excitation system are a subsystem based on the fact that the exciter controls both the flux and the voltage produced by a generator, and that a subcritical Hopf bifurcation is known to occur in these dynamics [25]. It is also known that a supercritical Hopf bifurcation can be observed in the generator mechanical, flux decay, and excitation system dynamics [25]. Since two Hopf bifurcations affect the exciter and flux decay dynamics, it might be difficult to define the proper subsystem for the subcritical and the supercritical Hopf bifurcations. These classes of bifurcation motivate the investigation of whether a subsystem composed of those dynamics produces and experiences the same bifurcation as the full system, what produces bifurcation in both and causes of the specific bifurcation, thus the name bifurcation subsystem. A nonlinear transformation can reduce a very complex high dimensional nonlinear system into a simple subsystem that describes the system behavior along the center 71 manifold. If one is able to retain how the parameters and equilibrium state change with the changes in operating conditions in the system, one could fully diagnose causes of the bifurcation and how it might be effectively prevented through direct control actions or through better adjustment or scheduling of the operating conditions. Although the subsystem that describes behavior at or near the center manifold is typically of very low dimension, the dimension of the participation factor subvector (obtained from normalized right and left eigenvectors of a specific mode) with large (> 0.1) elements can be quite large. From the example shown in Chapter 1, the right and left eigenvectors do not necessarily indicate the cause or the geographical location of the bifurcation, nor the states in the model where the bifurcation has an effect. In this system, loss of excitation voltage control on generators in one generating station on Vancouver Island causes a voltage collapse on Vancouver Island that uncontrollably spreads and brings about voltage collapse in the entire BC Hydro system. The Vancouver Island subsystem experiencing bifurcation may not be even retained based on the magnitude of the participation factor elements of the bifurcating eigenvalue before this loss of excitation control occurs which is shown in Figure 1.1. Once the loss of excitation occurs, there may be no load flow solution that can be used to calculate participation factor information, or eigenvectors. This result indicates that when discontinuities in the model exist and cause the occurrence of bifurcation [11] in a power system, the magnitude of the eigenvectors and participation factors may be insufficient or may give misleading information on the location, cause, or remedial action needed to prevent bifurcation. The result in [11] indicate that there are small coherent bus groups, called voltage control areas, with independent voltage collapse problems. Each of these coherent bus 72 groups in the electric transmission network is protected from experiencing a voltage collapse bifurcation by retaining reactive supply and thus voltage control on one or more of the generators in the unique subset of generators that protect that voltage control area from voltage collapse. This result indicates there are indeed small subsystems of a power system model, called bifurcation subsystems, that produce and cause the bifurcation that is experienced in the entire power system. (1) (2) (3) The reasons that bifurcation subsystems exist on a power system model are Bifurcations in a power system algebraic network model are typically brought on by a series of hard limit or equipment outage induced discontinuities. While hard limit actions are result of disablement of system (voltage) protective devices and can result in loss of voltage control, equipment outages such as the trip of a tie line result in decoupling or separation of two regions in the system. Therefore in each time epoch (If, 5+ 1) , the effect of the previous discontinuities which occurred at times (k for k=1,..,i-1 may fully be captured through the disappearance of certain dynamics, isolation of some subset of dynamics, increase in matrix sparsity when dynamic equations are substituted by algebraic equations,.. etc. The resulting subsystems within each epoch (If, 1,; I) may experience different kinds of bifurcations. Bifurcation subsystems with different kinds of bifurcations seem to be a natural and essential way for tracking the subsystem experiencing, producing, and causing bifurcation, as well as the actual cause for the final system instability; the extreme sparsity of the transmission network where every bus or node in a 10,000 bus network may only be connected to three to five other buses; the weak coupling between real power and voltage magnitude as well as the weak S! 260 in it PM deCO lIdtnt 73 coupling of reactive power—voltage angle. This decoupling implies the existence of a real power-voltage angle subsystem and a reactive power-voltage magnitude subsystem. Both subsystems have the extreme sparsity discussed above; (4) effective decoupling of coherent bus groups within the real power-voltage angle subsystem and within the reactive power-voltage magnitude subsystem. The decoupling of the two subsystems, the decoupling of coherent groups within each subsystem, and the extreme sparsity make a power system vulnerable to loss of control induced bifurcations; The bifurcation subsystem model is just a truncated portion of the actual power system model where the impact of parameter and operating changes simulated on a full system model can be easily observed and analyzed on the subsystem model. 4.0.3 Introduction to the Bifurcation Subsystem Method The bifurcation subsystem method utilizes the geometry associated with the various submatrices of the differential-algebraic Jacobian matrix I and with the eigenvectors associated with the bifurcating eigenvalue to establish conditions for existence of a bifurcation subsystem. Conditions for existence of a bifurcation subsystem require geometric decoupling that encompasses time scales, matrix decoupling and eigenanalysis in identifying the smallest subsystem that not only experiences the bifurcation but also produces and causes the same bifurcation observed in the full model. Geometric decoupling will be shown to be the most effective as well as a generalized method for identifying bifurcation subsystems. A bifurcation subsystem can therefore be loosely defined as a truncated portion of the .11 \E \lf bli‘dl ll) 1" .1 t... :3” SW: [lime )Iifl/ filer at” Osfic 74 the actual power system model that experiences and causes the same bifurcation in the full system model. Special characterictics of a bifurcation subsystem come from the fact that a bifurcation subsystem (1) is a reduced model developed for a specific kind of bifurcation since the geometric decoupling conditions for a bifurcation subsystem are only on the linearized model and only on the subspace associated with the left and right eigenvectors of the bifurcating eigenvalue as the bifurcation parameter p. approaches the bifurcation value 14*. The bifurcation subsystem should not only experience the bifurcation of the full system model, but may also be useful in capturing what produces and causes the bifurcation in any epochs (1;,1;+1) at the quasi equilibrium (39,511.), but also in tracking the model transitions which uncover the persisting dynamics that finally produced the collapse. (2) is not necessarily an equivalent model to the full differential-algebraic model in the sense that it preserves all the characteristics of the full model, since it can not be used to observe and investigate all other bifurcations of the full model. (3) geometric decoupling can break down at the point of bifurcation so that the entire system produces the bifurcation rather than just the bifurcation subsystem. The fact that the geometric decoupling at a quasi equilibrium for an epoch (1331;+ l) breaks down before the actual bifurcation occurs and does not always solely produce the system bifurcation may discredit application of the bifurcation subsystem method for these cases if one views the bifurcation subsystem in the strict sense of its definition. However as an engineering tool, the bifurcation subsystem method still provides diagnostic information on how much the geometric decoupling and how much the full shit" from future bli’dlt‘ him since J Also. i bifurca no pay pattern 4.0.4 1 ll. idt‘mifyj an equit CXlSlCnCe bifurtan‘g CKamplt‘s example i. bifumadOn These this km)“ k 75 system behavior after the breakdown contribute to the bifurcation in the full system model. It should be noted that geometric decoupling does not ever prevent the entire system from experiencing the bifurcation experienced in the subsystem, but hopefully can in the future help identify the subsystem that produces and causes the bifurcation. The fact that bifurcation subsystems experience and could produce most bifurcations even though the bifurcation affects all dynamic and algebraic states of the system should not be a surprise since almost all system bifurcations are understood in terms of affecting some subsystem. Also, it should be noted that one should not and can not assume all bifurcations develop in bifurcation subsystems, because the bifurcation could be produced by the full system, in no particular identifiable subsystem or subsystems, and in very strange unexplainable patterns of state variables. 4.0.4 Chapter Objective and Outline The three main objectives of this investigation are therefore to (a) develop tests for identifying subsystem dynamics that produce and cause a specific bifurcation to occur in an equivalent dynamic or equivalent algebraic system model; (b) utilize the test for existence of bifurcation subsystems to develop methods that can help identify all potential bifurcations that may exist in a particular model; and (c) introduce several application examples of bifurcation subsystems in stability assessment of power systems. One such example in this thesis is the justification of the classical load flow model as an algebraic bifurcation subsystem for load flow bifurcation of the differential-algebraic model. These objectives are quite ambitious Since the study of bifurcations usually assumes this knowledge while it is often unknown and is rather essential for excitation system red: km \01 mt die dc di CL fit Ct 76 redesign and FACTS control that may prevent the bifurcation from occurring. This knowledge may also be useful in developing unit commitment, generation dispatch and voltage setpoint schedules that could prevent each specific kind and class of bifurcation. The development of this bifurcation subsystem method will be first developed for a general differential-algebraic model and then applied to a Type 2 power system model. A Type 2 power system model is chosen rather than any other power system model type because (1) the bulk of bifurcation studies in power systems has been based on this model and (2) because the kinds and classes of bifurcations are best understood for this model. The differential equations in this Type 2 power system model represent generator mechanical and flux decay dynamics, generator excitation system dynamics, and governor turbine system dynamics. The algebraic equations represent the transmission and distribution networks as well as the load. The Type 2 power system model is described in details in Chapter 2 and structurally as a differential-algebraic model in Section 4.1. The first step in this investigation is to assume that it is possible to obtain a differential (dynamic) system model and/or an algebraic system ’ models, that are equivalent to the full differential-algebraic power system model. Such assumption is generally made because most bifurcations do not occur solely in the differential equations or solely in the algebraic equations of the differential-algebraic model. Formulation of equivalent differential and algebraic system models will be based on Schuur’s theorem as will be shown in Section 4.3.1. The conditions for existence of an algebraic bifurcation subsystem and a differential bifurcation subsystem which experience bifurcation in the equivalent algebraic system model and the equivalent differential system model, respectively, are then derived in Section 4.3.2. 77 The conditions for a differential bifurcation subsystem to exist are further refined to provide conditions for a subsystem of the differential equations to produce bifurcation in the equivalent differential equation model. When the bifurcation is experienced solely in a subsystem of the differential equations, then a bifurcation sub-subsystem is said to exist Conditions for a differential sub-subsystem are given in Section 4.4.1. The conditions for an algebraic bifurcation subsystem to exist are further refined to provide conditions for a subsystem of the algebraic equations to experience bifurcation in the equivalent algebraic system model. When the bifurcation is experienced solely in a subset of the algebraic equations, then an algebraic bifurcation sub-subsystem is said to exist. Conditions for algebraic sub-subsystems to exist are given in Section 4.4.2. The bifurcation produced in the full set of differential and algebraic equations can be experienced solely in a subsystem composed of a subset of the differential equations and a subset of the differential equations. When this occurs, a differential-algebraic sub- subsystem is said to exist. Conditions for a differential-algebraic subsystem to exist are given in section 4.4.3 Under normal conditions, a power system should never experience bifurcation solely in the differential equations, if the generators, governor and exciation system are properly designed. However, during voltage collapse, generator exciters cause the generator to reach field current limits in an effort to produce the reactive power needed by the network, and theexciter is completely disabled on excitation systems built prior to 1970 [51]. This can produce bifurcations that are solely in the generator dynamics. A condition for an algebraic bifurcation subsystem is also derived in section 4.3 where bifurcation in the equivalent algebraic system model is captured totally in the algebraic equations. Voltage 78 collapse is often studied using an algebraic model and thus its use is theoretically justified based on the conditions derived for an algebraic bifurcation subsystem to exist. This understanding motivates the development of Sections 4.3, 4.4 and 4.5. Finally, an application example of bifurcation subsystems in stability assessment of power systems, is given in Section 4.6. A classical load flow model is justified for the first time in terms of being both a bifurcation subsystem for saddle node bifurcation in the differential-algebraic model or equivalent differential model and a bifurcation subsystem for algebraic bifurcation in the algebraic equations of the network far from generator terminal and internal buses. In Chapters 6 and 7, the concept of bifurcation subsystems will be used to study several different kinds and classes of bifurcations in a Type 2 power system model. 4.1 Differential-algebraic Power System Model From power system stability assessment practice it is currently assumed that the underlying cause of voltage collapse is typically a bifurcation associated with the nonlinear power system dynamics, and can be invoked by a small change in system parameters. First, recall the nonlinear power system differential-algebraic model presented in Chapter 2, where the differential and algebraic nonlinear equations describe the generator and control dynamics, the network and load, respectively. The form of the nonlinear model in each time epoch (1‘1“,1;+ I) after the i’th discontinuity occurs in the system, 1 gig N — 1 is given by: Difierential equations: r,- = f1(x1(’))s)’1(’))sflr(‘)) te (1;:1;+ l) 4.1(a) Algebraic equations: 0 = gi(x,-(1), 311(1), 111(1)) 115 (1311:”) 4.1(b) Mn 1} IQ l where )3- 81' xi(t) : Yr“) 11,-(1) XxY 2 R I R N n+m+ n - - ' p —> R . an n-vector of real analytic nonlinear functrons + . . . " + m p —> Rm , an m-vector of real analytic nonlrnear functrons e x cR" , a state vector of dynamic state variables, typically time dependent values of generator and excitation system voltages, rotor phases and induction motor variables for (E (1;',1;+ l) . : e 1’ cRm, a state vector of instantaneous variables, typically time dependent values of voltage and angle at generator terminal buses, high side transformer buses and load buses for 1e (1f,1;+l) . e P cRp , a vector of slow varying operating parameters for the time epoch 1e (1f,1;+1) . It can represent a scaling of any of the parameters discussed in Chapter 3 as belonging to parameter space. The parameter “1“) is likely to scale parameters like reactive power load or generation, generator exciter voltage set points, or active generation or load. As Ill-(t) varies over the parameter space P, there must be at least one solution to the above differential-algebraic system of equations. : state space P : parameter space Note (Notation): The subscript ‘i’ used in the vectors xi(1) , yi(t), f i(x,-(1), yi(1), ui(1)),.. indicates that the model in valid after the i’th discontinuity, during the epoch (131:+ l) . For example, I xi(t) = “‘1’ x12, ""x‘n) e R" is the n-dimentional state vector during (1;,1;+1) . ><| 'tt euurt A fin: ht 0b series The J Uncut tompt there and [nthe tjrj1 0CL‘Ur \ and thfi SIand“ 80 For an arbitrary fixed parameter u‘. = H.» a steady state quasi equilibrium point X .= (time), 53115)), is a point of quasi equilibrium defined by locally solving the system of l equafions 0 = firm.» £111.). 11.) (49-3) 0 = 342411.). Time). it.) (421’) A linearized model of (4.1) about an isolated equilibrium point of interest (2., £11,) may be obtained using Taylor’s theorem where every function is approximated by its Taylor series expansion of the first degree, in the neighborhood of the isolated equilibrium point. The Jacobian matrix of (4.1) at an equilibrium point (Ivy—1,11,) may be obtained by linearizing the differential and algebraic equations and then combining all the linearized [A151] = Ji- [AL] 4.3(a) O Ayi J _ afi/axi afi/ayi = f‘x“ fly“ ‘ egg/3% agi/ayi (:1" )7? He) gixi 3in components, as shown below. where 4.3(b) and Axi = rim-ft,- Ayi = )2“) " )7, In the rest of the thesis, sufficient smoothness of the functions f i(xi(1)), yi(t)), ui(t)) and gi(x,-(1)), yi(1)), lit-(0) is assumed to be preserved since no discontinuous system changes occur within the time epoch (1;',1;+ l) . Therefore, for simplicity, we omit the subscript i and the argument ‘t’ in the analysis within the epoch. That is, unless otherwise noted, g stands for gi(xi(t)),y,-(1)), ui(1)), fx stands for fi(xi(t)), yi(1)), [rt-(1)) ,.. etc. A linearized differential-algebraic power system model has been derived in [52] for \lller 81 a general Type 2 power system model, where each synchronous machine and its control systems are described. The control systems include (a) Excitation control with Load (Line-Drop) Compensator, (b) Power System Stabilizer (P88) [52], and (c) Speed- Govemor-Turbine System. The full linearized system model [52] is given by: TXXAXX Axx AXE Axe 0 Axe Axv FAX; TEEAXE AEX A55 0 Ass A59 AEV AXE TGGAXG Aax 0 A06 0 0 0 AXG = +BOAU0 (4.4) TSSAX'S Asx 0 0 Ass 0 0 AXS -APC APX 0 0 0 Ape APV A9 Where AXX : states of mechanical and flux decay dynamics; AEF : states of flux decay dynamics AXE : states of excitation system AX G : states of speed-governing-turbine systems; AX S : states of power system stabilizer; A9 : angle variables at network buses; AV : voltage variables at network buses; APC : coefficients of non-voltage-dependent active power load demand model; AQC : coefficients of non-voltage-dependent reactive power load demand model; T x x : diagonal matrix composed of inertia constants of synchronous machines; identity matrix, and time constants of flux decay dynamics; TEE : time constants of the excitation systems dynamics; Too : time constants of speed-governing-turbine systems; 82 T55 : time constants of power system stabilizers; 4.2 Conditions for Bifurcation Subsystems Experiencing SN and LF Bifurcations 4.2.1 Equivalent System Models Although testing for singularity of the full Jacobian matrix J in (4.3b) is a sufficient indication of a system bifurcation, this test is not sufficient to tell (a) whether either saddle node bifurcation (when gy is nonsingular), load flow bifurcation (when fx is nonsingular), or algebraic bifurcation ([gx gy ] is row dependent) occurred; (b) whether the bifurcation is associated with strictly an equivalent dynamic system model or an equivalent algebraic system model; or (c) what system dynamics were associated with the bifurcation. Schur’s theorem [53] is used to establish conditions under which equivalent systems to (4.1) comprised solely of differential equations or solely of algebraic equations can be constructed with different eigenvalues than the differential equation model which better reflect the effects of the eigenvalue qualitative changes of the full system model. Existence of equivalent differential or algebraic system models to (4.1) depends on the row dependence of the rows associated with the differential equations and/or the row dependence of the rows associated with the algebraic equations in the linearized full system (4.3). This will be clear from Schur’s theorem which reduces the computation of a determinant of a matrix of order n to the computation of a determinant of a matrix of order r 11* . Under these conditions, singularity of the reduced order matrix fx becomes crucial for analyzing and studying the occurrence of saddle node bifurcation in the equivalent dynamic model of (4.3). The reduced order system is called a bifurcation subsystem of first order for (4.1), since its change in dynamic behavior experiences the change in dynamic behavior in the full system. Recall conditions (i) and (ii) of Proposition 4.1 v,’ - tf,g;‘g,1 = o tf,g;‘g,1 - u, = 0 as u-w" (4.9) These equalities are seldom perfectly satisfied in a real power system Type 2 model. However taking into account the extreme sparsity, decoupling of reactive power and voltage from real power and angle, and the decoupling of coherent bus groups in both the (It In. mal 101 real power and angle as well as reactive power and voltage model, (i) and (ii) are likely to be satisfied to an 0(8) as discontinuities occur that cause J to be singular when fx is nearly singular. Hence, the subsystem Ax = f, - Ax , may effectively be used to assess the causal factors of the bifurcation in the full model. Equations in (4.9) are approximately satisfied for the following five conditions: (a v 1ng n is 0(a) and llg;l||"fy|| - "u,“ is 0(1) near 11* (b) Ilfy n is 0(a) and ||g;‘||"g," - "u,“ is 0(1) near 11* (c) up u is 0(a) and “hung," . "a," is 0(1) near to (d) None of the above but near it“ , the vector v, satisfies vf- [fyg;‘gx] z 0 (e) None of the above but near 11* , a vector u, e N, satisfies [fyg;lgx] - u, = 0 where H . II is any suitable matrix norm. Proposition 4.2 Let gy(tt*) in 4.3(b) be singular for some 11* e ([1,, u?) , fx nonsingular near 11* and let N, = Null(gy(lt*)) and N2 = Null(gy(p.*)T). Then, = gyul) - Ay is an algebraic bifurcation subsystem of first order for (4 .l ) experiencing LF bifurcation if and only if near 111* , condition (i) or (ii) is satisfied (1') [gxf;lfy] - u2 = 0 for some 11215 N, (ii) v§.[gxf;1fy] = 0 forsome v26 N2 Implication of Proposition 4.2 In contrast to (4.9), conditions (i) and (ii) of Proposition 4.2 may be satisfied when the matrix fx associated with the machine dynamics, exciter and governor controls is diagonally dominant, a condition which is generally satisfied in a real power system Type 102 2 model. It can be shown that conditions (i) and (ii) of Proposition 4.2 can be satisfied to an O(e) only in Section 4.6.5 if the bifurcation occurs near or in the distribution system and electrically far from generator terminal and internal buses. The resulting bifurcating model 0 = gym) - Ay may then effectively be used to assess the causal factors of the bifurcation in the full model. Structural and operating conditions and constraints for the validity of using the algebraic bifurcation subsystem will be discussed in Section (4.3.5). Discussion of Propositions 4.1 and 4.2 The geometric conditions that near 11* , [gxfjfy] - u2 = 0, v; - [gxfjfy] = 0, v? - [fyg;lgx] = 0, and [fyg;lgx] ~ u, = O generalize the case of full decoupling of the differential equations and the algebraic equations stated in Corollaries 4.1 and 4.2, while Corollaries 4.1 and 4.2 give conditions for the existence of equivalent reduced-order systems to the full differential-algebraic model and would never occur. The more general geometric conditions given in Propositions 4.1 and 4.2 are for the existence of bifurcation subsystems and can indeed occur. (a) Sparsely connected subsystems Although these conditions seem to be a strong requirement, they are indeed a quite easily satisfied requirement since the submatrices of J and especially g, have few nonzero elements, so that there is a natural column or row dependence in the submatrices of J, because so many elements of each column and row are identical or are zero. In a 10,000 node network, each node may be connected to at most five other nodes making gy virtually null except for five nonzero elements in the block (b) (C) 103 submatrices corresponding to the Jacobians dF/de , dP/uv , dQ/de and dQ/dv . The submatrices of fx in (4.4) are diagonal which makes f, very sparse as well. Weakly connected subsystem Weekly connected may be expressed by elements in off-diagonal matrices being multiples of a small scalar e > 0. If some of the few nonzero elements in different rows and columns are 0(8) and can be approximated by zero, then the combination of sparsity and weak coupling can produce row or column dependence. The natural decoupling of real power-angle submodel and reactive power-voltage submodel as well as the decoupling within bus groups of both these models satisfies the weak coupling assumptions. Discontinuous system changes Time epochs (1,.+,1,r+ ,) are defined in Chapter 3 as the time interval after the i’th discontinuity and before the i+l hard limit or equipment outage discontinuities. These nonsmooth system transitions described in Chaptre 3 usually result in loss of voltage control of certain buses in the system and/or decoupling and islanding of certain areas or group of buses of the system. Loss of control implies that the generator dynamics and controls or other voltage protective devices of the system are disabled which results in a decoupling of the generator dynamics and the system algebraic operational constraints. Loss of control induced discontinuities occur when generators reach field current limits that cause the exciter to effectively lose control of voltage, tap position limits to cause tap changer controls to lose control of voltage 01 If \"C de 104 and switchable shunt capacitor susceptance limits that also cause loss of control of voltage discontinuously affect the structure and in some cases change the dimension of f, and g, . On the other hand an equipment outage such as the trip of a tie line, i.e a line that connects two major areas in a‘power system results in decoupling between the two regions or even an isolation or islanding of a major part of a system. Therefore in each time epoch (1;,1;,,) , the effect of the previous discontinuities which occurred at times 1,, for k=l,..,i-l may fully be captured through the disappearance of certain dynamics, decoupling or isolation of some subset of dynamics, increase in matrix sparsity when dynamic equations are substituted by algebraic equations,.. etc. For example, the effects of the discontinuities in epoch (1 ,t, 1;, , ) , in addition to the extreme sparsity in the network create the weak coupling of the reactive power and voltage and real power and angles as submodels in g,,, as well as the weak coupling of coherent bus groups in both models. Therefore, the resulting model within each epoch (1;,1;, ,) may experience separate bifurcations in the real power inertial (angle) dynamics, flux decay, exciter, voltage and reactive power dynamics within coherent groups. The evolution of the bifurcation subsystems with different kinds of bifurcations for each discontinuity is a natural and essential way for tracking the actual cause for the final system instability. The relative isolation of the real power and angle coherent groups and reactive power and voltage coherent bus groups in submodels of g), allows Proposition 4.1(i) and (ii) to be satisfied for eigenvalues associated with the inertial dynamics and for the flux decay exciter dynamics in coherent groups of buses. U34 105 (d) Geometric Structural Decoupling The decoupling may only be geometric and thus occurring in the direction of the bifurcating eigenvalue. The geometric decoupling may occur as u —+ (1* , or may be such that it holds for all u as u—>u*. This geometric decoupling, asymptotic as u—iut' or structural for all u, is quite easy to achieve given factors (a-c). It is often impossible to determine all of the bifurcation subsystems without comprehensive stress tests at every bus in the network and possibly every subsystem of the generator model since the coordination of controls in different bifurcation subsystems mask the existence of such subsystems. Comprehensive stress test cause the discontinuities that disable the controls to allow one to observe the subsystem or sub- subsystem that experience bifurcation, which result when conditions of Proposition 4.1 and 4.2 are satisfied. 4.2.5 Algebraic Bifurcation Subsystem Experiencing LF Bifurcation Until the last decade, the classical method for analyzing local power system stability has been limited to the use of a load flow model. This model commonly called classical or standard load flow model, will be shown in Section 4.5 to be an approximation of a linearization of the Exact Load Flow model 0 = g(x(y).y(u).u) Under the Conditions of Proposition 4.2, the linearized Exact Load Flow model is shown to be a bifurcation subsystem of first order for the differential-algebraic model. The P-V or Q-V curves computed using the classical load flow model have been the principal methods used to study voltage stability [App B] and proximity to voltage collapse. Studies [54] 106 have shown that voltage collapse occurs at a bifurcation in the classical load flow model and results either in no solution or multiple solutions in the model equations. The results in [l 1] show that each coherent bus group where the Q-V curve have nearly identical minima have unique voltage collapse problems. These voltage collapse problems occur whenever the generators that exhaust reactive reserves and lose control of voltage in computing the Q-V curve minima for each coherent group exhaust reserves for equipment outages or operating changes. The fact the classical load flow is an approximation of an algebraic bifurcation subsystem of the differential-algebraic model supports the discussion of how bifurcations occur due to extreme sparsity, decoupling of real power-angle and reactive power-voltage submatrices, weak coupling of coherent groups in each submodel and loss of voltage control discontinuities actually cause bifurcation in the differential- algebraic model as claimed in Section 4.2.4. 4.3 Bifurcation Sub-subsystems Experiencing SN and LF Bifurcation The results presented in 4.6.5 indicate the classical load flow model is a bifurcation subsystem when generator exciters have high gain and are not disabled and results presented in Chapter 5 indicate that the reduced differential equation model Ax = J,- Ax can be a bifurcation subsystem of first order. when generator exciters are disabled. The fact that the classical load flow is a bifurcation subsystem that produces and cause bifurcation in the differential-algebraic model does not end the application or investigation of bifurcation subsystems. The bifurcation subsystem method will attempt to identify algebraic sub-subsystems of smaller and smaller dimension to determine the smallest l()7 bifurcation sub-subsystem (a coherent bus group) that produces and causes the bifurcation in the complete system as will be shown in Chapter 5. Similarly, differential bifurcation sub-subsystems containing (a) generator mechanical dynamics, (b) electrical dynamics, (c) generator flux decay dynamics, ((1) the power balance equations at terminal buses, (e) the control systems dynamics or (f) a combination of the above would be investigated to find the smallest bifurcation subsystem once it is confirmed that a differential subsystem bifurcation has occurred. If a differential or algebraic bifurcation subsystem does not exist when conditions of Propositions 4.1 and 4.2 do not hold, a differential-algebraic bifurcation subsystem containing both differential and algebraic equations may exist For the following results on finding the smallest algebraic subsystem, the smallest differential bifurcation subsystem and the smallest differential—algebraic subsystems, the Jacobian matrix J must be partitioned as follows x11 ler12 fyrr fy12 f f J fxzr x22 yzr fyzz g g,” gle Yu gym ngI g‘22 gyzr gyzz 4.3.1 Differential Bifurcation Sub-subsystem for Saddle-Node Bifurcation Not all system dynamical states are generally involved in a differential subsystem bifurcation. The differential subsytem Jacobian matrix fx includes the machine electrical and flux decay dynamics, excitation system dynamics, governor control dynamics and power system stabilizer dynamics. The matrices associated with each set of dynamics are 108 shown in the partitioned fx matrix. The matrix fx“(p.) in 1 (WI (u) 1]=1,..sx= ‘11 ’12 [Ail] (4.12) 2 13,2101) fxnm) Ax, [Ax Ax’ = , can represent any subsystem of the the differential equation in (4. 12), and the conditions for differential bifurcation sub-subsystems of first order are given in the following two corollaries. Corollary 4.3.1(a) Given that (4.1) has a differential bifurcation subsystem for SN bifurcation because conditions of Proposition 4.] hold. Let fx”(u*) in (4.12) be singular for some 11* 5 (111.112) where transversality conditions for SN hold at if" , fx,, nonsingular near u* and let N, = Null(fx“(,u*)) and N, = Null(fx”(tt*)T). Then Ax, = fx“(u)-Ax, is a difierential bifurcation sub-subsystem of first order for (4 .1 ) experiencing SN bifurcation if and only if near 11*, condition (i) or (ii) is satisfied (1') [3,12an fle] - u, = Oforsome 11,6 N, .. T -1 . Corollary 4.3.1(b) Given that (4.3) has a drfiferential bifurcation subsystem for SN bifurcation because conditions of Proposition 4.1 hold. Let fx22(u*) in (4.12) be singular for some 11* e (u,, u,) where transversality conditions for SN hold at X“ , fx” nonsingular near 11* and let N, = Null(f([.t*)) and N, = Null(f(l.t*)T). Then Ax, = f Jr220.1) - Ax, is a difl'erential bifurcation sub-subsystem of first order for (4.1 ) experiencing SN bifurcation if and only if [(19 near 11* , condition (i) or (ii) is satisfied - —1 (z) [“21fo [X12] - 11, Oforsome u, e N, .. -1 (ll) V{'[fx2,fx,, [X12] = Ofor some v, 6 N2 4.3.2 Algebraic Bifurcation Sub-subsystem for Load-Flow Bifurcation Let the algebraic bifurcation subsystem 4.3(b) be partitioned as 0 = gym = [g’118’12J.[Ay1] (4.13) g’zr gyzz Ay, - where gm , gym , g),21 , g,22 are arbitrary submatrices of 3),. Depending on the singularity of g,“ and gy22 , algebraic bifurcation sub-subsystems experiencing LF bifurcation may be obtained by applying Propositions 4.1 and 4.2 to the system in (4.13) as given in the two corollaries below. Corollary 4.3.2(a) Given that (4 .3) has an algebraic bifurcation subsystem for LF bifurcation because conditions of Proposition 4.2 hold. Let gy”(|.t*) in (4.13) be singular for some we (u,,u,), 11,2201) in (4.13) nonsingular and let N, = Null(gy“(tt*)) and N2 = Null(gy”(|.t‘)T). Then 0 = gynut) . Ay, is an algebraic bifurcation sub-subsystem of first order for (4.1) experiencing LF bifurcation if and only if condition (i) or (ii) is satisfied . —1 (l) [gylzgyzzgyu] - u1 = Ofor some ul 6 N, -- T -1 (u) v [ ]=0 orsomeveN 1 smgyzzgyu f 1 2 110 Corollary 4.3.2(b) Given that (4 .3 ) has an algebraic bifurcation subsystem for LP bifurcation because conditions of Proposition 42 hold. Let gy22(|.t*) in (4.13) be singular for some at e (u,,11,), gin nonsingular near 111* and let N, = Null(gy22(tt*)) and N, = Null(gy22(tt*)T). Then 0 = gynm) - Ay,2 is an algebraic bifurcation sub-subsystem of first order for (4 .1 ) experiencing LF bifurcation if and only if, near 11* , conditions (i) or (ii) is satisfied 0 forsome u, e N, - —1 (I) [gyzrgyugylz] “1 (ii) vT-[g 3“,; ] er yrl'ylz Ofor some v, 6 N2 The thesis will show that the disablement of generator excitation control produces voltage collapse in a load flow model, which is an algebraic bifurcation subsystem and will also cause bifurcations in the differential subsystem composed of the generator inertial and flux decay dynamics. These bifurcations in this differential bifurcation subsystem can explain the dynamical behavior observed as a system experiences voltage collapse. 4.3.3 Differential-algebraic Bifurcation Sub-subsystem Differential-a1 gebraic subsystems may be the most common type of bifurcation subsystem that produces and causes bifurcation in the full differential algebriac model. It has been reported that the transfer of real power across long transmission network is the reason behind the low-frequency oscillations, which may affect one or many generators in the system, or may affect more than one coherent generator groups. Differential-algebraic subsystems may therefore exist to represent not only the generator and control dynamics involved in the bifurcation, but also the contribution of the power balance constraints in lll causing the failure. These subsystems may be derived by first permuting the order of the state variables as for instance, 1x ylr= [x115 Y1 by —>[x1y1x2 ill Then the full system (4.3) now has the form r - ~ f f f f r . r - ,7 x x y Ax, g 11 gym 8 12 g 12 Ax, .. . Ax, 0 = .x11 )’11 x12 in , Ay, = fx fy , Ay, (414) . t t sz fx fy fx fy M2 8x g), M2 0 21 21 22 22 Ay A)’ “ ‘ g g g s L 2‘ L 2‘ _ ‘21 hr ‘22 222 Differential-algebraic sub-subsystems may be obtained by applying Propositions 4.1 and 4.2 on the Jacobian matrix as partitioned in Equation (4.14), where f;, f; , g; and g; are given by: pr1 f)"; 1 fX" fr” gx = -1 ..l ’ 8y «2 ~~ 4.15(a) £3112, g)’31_ ng22 gy22 - f - f e x y " x y fx = 11 11 ’ fy 12 12 4.15(b) ‘11 8)”, x,2 gy12 Corollary 4.4.3(a) Let f;(u*) in (4.14) be singular for some 11* 6 (111,112). g;(u) in (4.15) nonsingular near 11* and let N, = Nuu(f;(11*)) and N, = Nu11(f;(u*)7) . Then Ax, e AX, = f (11) ~ [ 0 i x [Mr] is a differential-algebraic bifitrcation subsystem of first order for (4 .1 ) if and only if near 11* , condition (i) or (ii) is satisfied ‘ ‘ fi-l . - (I) lfyg, gxl-u, = Oforsome u, e N, 112 (ii) 1" 1 - [figf‘sfél = 0 for some V, 6 N2 Corollary 4.4.3(b) Let g;(u*) in (4.14) be singular for some 11* e (11,, 11,). f;(u) in (4.15 ) nonsingular near 11* and let N, = Null(g;(tt*)) and N, = Nu11(g;(11*)T). Then . . A 0 Ah is a diflerential-algebraic bifurcation subsystem of first order for (4 .1 ) if and only if near 11* , condition (i) or (ii) is satisfied (i) 1g;1;"1;1 . u, 0 for some it, e N, —1e -- T e ‘ (II) VI ~111fo fyl Oforsome v, e N, 4.4 Differential Bifurcation Subsystems for Hopf Bifurcation The analysis and results carried out so far in this chapter have been based on Propositions 4.1 and 4.2 which focused on the system’s change in behavior when a simple real eigenvalue of the differential-algebraic Jacobian matrix crosses the jw—axis. It has been shown that under certain system operating conditions, two bifurcation subsystems of 4.3(a) are of interest: (a) An algebraic bifurcation subsystem for studying load flow bifurcation (also called loss of causality) in the equivalent algebraic system model (4.5) is given by O=Jy-Ay = [gy-gxf;lf,l-Ay = gy-Ay (b) A differential bifurcation subsystem for studying saddle node (SN) bifurcation in the equivalent differential system model (4.6) is given by Ax’=Jx-Ax=[fx-fy8;18x1'Ax = f, Ax Such results are desired and may easily be extended to the case where the system behavior 113 is characterized by the onset of unstable (subcritical) or stable (supercritical) oscillations. With no further analysis Proposition 4.1, Corollaries 4.4.1(b) and 4.4.1(b) may be generalized to cover the case when the equivalent dynamic system Jacobian matrix J, has a pair of pure imaginary eigenvalues is, , = iiwo. If (00:0, the equivalent dynamic system model experiences SN bifurcation, otherwise the system experiences undamped oscillatory behavior whose frequency depends on too. Also it is evident that if the system Ax = 1,01!) - Ax has a pair of pure imaginary eigenvalues "1.2 = :in with associated conjugate (complex) eigenvectors 11,, , e C" then we have l,-u, = two-u, = Jx(u*)-u, [Jx(|-‘*)"“°o’nxn] ' “r = 0 which says that when the matrix I, has a pair of pure imaginary eigenvalues 71,, = :in , the matrix [1,(u* )riwolm] has a zero eigenvalue, i.e. singular. Therefore, with the following definitions, ffflu) = f,(u)-iw,l,,,,,, 4.16(a) Jf‘out) a 1410-1111051”, 4.16(b) _| . = f,(u)-f,g, griffin/"rm, = flaw-1,59, similar proposition and corollaries to Proposition 4.1 and Corollaries 4.4.1(a-b) may be stated for the existence of differential bifurcation subsystems of 4.3(a) for Hopf bifurcation. Proposition 4.3 Let ffflu) in (4.16) be singular for some 11* e (11,, 11,) where transversality conditions for Hopf hold at 2*, gym) nonsingular near 11* and let N, = Null[f:)°(|.t*)]ch and N, = Null[ffO(u*)T] c C". Then, Ax = fx(u) - Ax is a difi’erential bifurcation subsystem of 114 first order for (4.1) experiencing Hopf if and only if near 11*, condition (1‘) or (ii) is satisfied (1') {1,518,} ' u, = 0 for some u, e N, .. T . —1 (U) v r-[fygy 12,] = Oforsome v, 6 N2 4.4.1 Differential Bifurcation Sub-subsystems for Hopf Bifurcation The Hopf bifurcation in a differential bifurcation model for Hopf, may exist in only a subset of the generator inertial dynamics, flux decay dynamics, exciter system control dynamics, governor control dynamics...etc. Condition for bifurcation of Ax = fxur) . Ax or bifurcation sub-subsystems of (4.1) for Hopf are given below. Let the matrix ffom) be partitioned as (D f,"(11) f f?"(u)-Ax= ” “2 [1111,] (4.18) . m0 Ax, 1,2, 52,111) then conditions for differential bifurcation sub-subsystems for Hopf bifurcation are given in the following two corollaries. Corollary 4.4.1(a) Given that (4 .3) has a dijferential bifurcation subsystemof first order for Hopf bifurcation exists because conditions of Proposition 4.3 hold. let 1'2", (11*) in (4.18) be singular for some 11* e (11,, 11,) where transversality conditions for Hopf hold at 2*, figm) nonsingular near 11*, and let N, = Null|:f:)"(u*)] and N, = Null[(f:)0(u*))1]. Then 11 11 Ax, = fxuut) - Ax, is a differential bifurcation sub-subsystem of first order for (4.1) for Hopf if and only if near 11* , condition (i) or (ii) is satisfied 115 (i) folzfégfll) 1,2,] - u, = 0 for some u, e N, .. 1' to -1 (ll) v, -[fxlzf,2"2(u) [X21] = Ofor some v, e N, Corollary 4.4.1(b) Given that (4 3 ) has a dtfierential bifurcation subsystem for Hopf bifitrcation because conditions of Proposition 4.3 hold. Let f2; (11*) in (4.18) be singular for some 11* e (u,, 11,) where transversality conditions for Hopf hold at 71*, fzolm) nonsingular near 11* and let N, = Null[f:;02 (11*)] and N, = Null[(f:)2°2(|.t‘))r]. Then Ax, = fx22(t1)-Ax,, is a dijferential bifurcation sub-subsystem of first order for (4.1) for Hopf bifurcation if and only if near 11* , condition (i) or (ii) is satisfied (i) 1%,,fo (1104fo - u, = Ofor some 11, e N, l .. T 020 —l _ (II) V, -[fx2,fx”(u) fx12]- Oforsome v, e N, 4.5 Application of Algebraic Bifurcation Subsystem Experiencing LF Bifurcation in Power System Stability Assessment When the linearized system model (4.4) satisfies conditions of Proposition 4.2, the algebraic bifurcation subsystem of first order 0 = g, . Ay or consisting only of the network equations 0 = g(x, y, p.) is a bifurcation subsystem of the Exact Load Flow model 0 = g(x(y), y, u) that is used to study load flow bifurcation occurring in the full [if] 1:1 - Considering the similarity between g, and the classically used approximate load flow differential-algebraic model model with Jacobian matrix 1,; [52], such a result would justify why experience has shown that a classical load flow model is such a useful tool in assessing stability that 116 occurs in a differential-algebraic model. The classical load flow model was the first and for some time the only tool for studying voltage collapse or voltage stability [App B]. The P- V or Q-V curves computed using a classical load flow model have been the principal methods used to study voltage instability and proximity to voltage collapse. Studies [54] have shown that voltage collapse occurs as a bifurcation in the classical load flow model and results either in no solution or multiple solutions in the model equations. Over the past ten years far more accurate Type 2-4 differential-algebraic models have been developed. However the classical load flow model is still the model and computer program being used in power system voltage stability assessment programs despite the existence of far more accurate models and simulation tools. The objective is therefore to use the present knowledge of bifurcation subsystems, to assess when and why a classical load flow model is a bifurcation subsystem for the full differential-algebraic model. The discussion will start with a brief description of a classical load flow model, then take hierarchical steps by first looking into the preliminary system condition (f, nonsingular) for utilizing J, as a test matrix that is equivalent to J for testing when a static or algebraic bifurcation occurs in the differential-algebraic model. The second step is to compare the Exact Load Flow model 0 = J, - Ay ( 0 = g(x(y), y(tt), u) ) with the algebraic bifurcation subsystem 0 = g), - Ay (0 = g(x(11), y(u), 11) ), and the classical load flow model 0 = JLF ~ Ay (0 = g’(y, 11)). Finally, based on the bifurcation subsystem conditions, some operational and geometric grounds are established for the use of the classical load flow power system model when algebraic bifurcation subsystem experiences the same algebraic bifurcation as the Exact Load Flow model. 117 4.5.1 Description of 3 Classical Load Flow Model The classical load flow model describes the power system steady state condition of the transmission network and is a set of algebraic equations composed of the active and reactive power balance equations of the transmission network at load buses and only the active power balance equations at generator terminal PV buses. The reactive power balance equations at generator buses are neglected since the excitation systems of these generators are assumed to have infinite gain and thus hold the generator terminal voltage to set point value. The classical load flow model can also be viewed as a simplified model of a power system differential-algebraic model under the following assumptions: (a) the generator armattu‘e 1esistaz..'e is zero; (b) fx is invertible (c) the excitation gain of each generator is infinite; and (d) all time constants in the differential equations are set to zero since the time interval over which the model is to be valid is infinite. Buses in a classical load flow model are conventionally divided into three types: (1) a swing bus is a generator bus with infinite active and reactive power reserves so that its voltage and angle are fixed while its active and reactive power generations are dependent variables that provide the power mismatch between generation and load plus losses in the rest of the system; (2) a generator PV—bus represents a generator bus operating under over- excitation limits (OEL), and thus its reactive power generation and angle are dependent variables while its voltage and real power generation are independent variables. (3) load PQ-bus is any bus with no reactive or real power generation, and thus instead active and reactive power injections are specified at this bus and are independent variables in the model. The voltage and angle at PQ buses are dependent variables that must be solved for using the model. 118 The classical load flow model consists of active power balance equations at all PV- buses and PQ-buses, and reactive power balance equations at all PQ-buses. The reactive power generation at each PV-bus and the active and reactive power generations at the swing bus are output variables (dependent variables) of the load flow model. In addition, if a generator PV-bus exciter is disabled by field current limit controllers, the generator PV- bus becomes a load PQ bus due to action of field current limit controller that disables the exciter and controls the reactive power to be the continuous rating limit. Thus, one more reactive power balance equation at that generator terminal bus needs to be added, since the voltage magnitude at its terminals is no longer held fixed. The general structure of the classical load flow Jacobian has been derived in [52] and presented in Appendix A as: ’1”. Ar 31H 0 C19 011-! 0 3P 8P 1”” A21! BZHH ‘3an C519 021m 02111. JLF = E W ___ J“ = 0 BZLH 321.1. 0 DZLH DZLL (4.18) fig 9!: 1’"? A59 3179 0 C13? 0139 0 393V 19” A B B CPfiD D JQL 4H 4HH 41.11 4 41m 4111. L - 0 B B 0 D 0 Where _ 4LH 4LL 41.14 41.1., I 9 = 191‘ 9”! 9L1] : angle variables at network buses; t v = [vat V”! VL‘] : voltage variables at network buses; t P = [P1, PH, [31“] : voltage-dependent active power load demand; I Q = [QfQ' Q”: M2131 voltage-dependent reactive power load demand; The real and reactive power balance equations in ( 4.18) have been arranged according to the types of bus variables in the transmission network, i.e generator terminal (T) buses, high side (H) transformer buses and load (L) buses. The matrices in the fourth row block and column blocks are solely associated with the generator PV buses which became PQ 119 buses, thus the superscript notation PQ. 4.5.2 Structural Conditions for Utilizing Jy The structural condition (f, nonsingular) imposed by Proposition 4.2 for formulating the algebraic bifurcation subsystem 0 = g, . Ay may be understood by observing the matrix f, as presented in (4.4). The derivation of the model used is given in [52] and summarized in Appendix A. The matrices fx and ACC are given by r r Axx Axe Axo 0 ~ “ A15x ABE 0 ”‘53 OX 66 Asx O 0 Ass ABS 0 Ass 0 0 ASS Matrix fx is composed of diagonal block submatrices where each diagonal element is associated with one machine. The elements of each submatrix of f, is defined based on the single machine infinite bus matrix, since it represents the operating condition when the terminal bus of each machine is assumed to be an infinite (swing) bus with constant terminal voltage. Note that each diagonal block submatrix of fx can be viewed as an element of a single machine model, since the dynamics of each machine are uncoupled in fx. The control systems form an upper triangular block matrix Ace and its eigenvalues are those of each control system: AEE , AGO , and Ass- Based on IEEE general models of speed-goveming-turbine system and power system stabilizers [52], both Ace and ASS are lower triangular matrices with real and negative eigenvalues. Regarding the excitation 120 system matrix AEE, when the excitation system is “properly tuned” and is able to be manually operated, all the eigenvalues of Aer—2 are nonzero negative real values. Moreover, under the assumption that each synchronous machine is operated below its steady state stability limit, the matrix f, can be shown to be nonsingular, with stable eigenvalues. It will be shown in Chapter 6 that fx can be singular if the mechanical dynamics or flux decay dynamics experience saddle node bifurcation. Operationally neither saddle bifurcations should occur if the machine is operated within its real power generation limits and the exciter is properly designed. Therefore, under the above assumptions, the equivalent algebraic Jacobian matrix Jy = gy — gxfjfy can be analytically defined. 4.5.3 Comparison of Exact Load Flow model, Algebraic Bifurcation Subsystem, and Classical Load Flow Model The Exact Load Flow model 0 = g(x(y), y(u), 11) , algebraic bifurcation subsystem 0 = g(x(u), y(u), 11), and the classical load flow model 0 = g’(y,u) are all quite similar because they all have real and reactive power balance equations at generator high side transformer buses. The associated linearized models with J acobians J), g, and J”: are also quite similar. The diflerences in J, and g), is investigated first to establish (i) when and if load flow bifurcation occurs due to loss of causality so that 0 = g(x(u),y(1t), u) is an algebraic bifurcation subsystem of 0 = g(x( y), y(u), u), and (ii) when the classical load flow model is a bifurcation subsystem model of the Exact Load Flow model. The geometric conditions of existence of an algebraic bifurcation subsystem of first order for LF bifurcation. —1 [exfx f,l-u,=0 when gy-u,=Jy-u,=0 121 1"2-[gxf;'fy] = 0 when v, -gy = v, J, = 0 where u, and v, are the right and the left eigenvectors associated with the zero eigenvalue. may be approximately satisfied in different ways such as (Cl) 1ng II is 0(8) and "ff““fy” - “a," is O(1) near 11* (C2) llfy II is 0(8) and ||f;‘||||g,|| - “a," is O(1) near 11* (C3) riff" is 0(a) and "f,““gxu - "u,“ isO(1) near 11* (G1) None of the above but near 11* , the vector v, is such that vT2- [,g,,f;l fy] = 0. (G2) None of the above but near 11* , the vector u, is such that [gxfffy] - u, z 0. where 11 . II is an appropriate matrix norm. Conditions C1 and C2 pertain to the case when the coupling matrices g, and f, are negligible so that the interaction between generator internal and terminal buses is insignificant to the bifurcation study. Likewise, condition C3 implies that the single machine matrix f, must be diagonally dominant. The eigenvalues of j; are in fact proportional to the eigenvalues of the control systems AGC , A S S and A 125 , and the system K -matrices (will be derived in Chapter 6). These eigenvalues may be high when the system is operating below its steady state stability limits and with high excitation gain. Despite the implication of these conditions, it is clear that the geometric conditions G1 and G2 not only encompass conditions C1-C3, but also are more likely to occur due to the sparsity and weak coupling property of the Jacobians as noted earlier. G1 and G2 say that near 11* , [gxf;lfy] need not be zero, but only that gy-Jy be nearly singular and that its null space (approximately) overlap the null space of gym") . Conditions G1 and G2 require examination of the structure and characteristics of the matrix [gy—Jy ] = gxfjfy near 11* and its geometric interaction with the right and left eigenvectors u, and v, of the null space of gy(|.l*) 122 In order to demonstrate and characterize the validity of the use of the algebraic bifurcation subsystem and subsequently the classical load flow model, analytical expressions of g, and [11,, —Jy ] = gxfjfy derived and presented in [52] are ng FAIK BIH 0 CIK Dru 0 8 £5” A211 32111! 321.11 CZH 021m Dzm. _ 145’” 3 0 321.11 321.1. 0 DZLH 021.1. (420) 35” A3K 33H 0 C3K 03H 0 QJQJQJ 9.5110 CD11: on we) Q) <|w <14 3,9” A4” 34111-1 341.}! C411 D411” D4HL £321. ,0 341.11 341.1. 0 041.11 D4LL_ Kgooxgod oooooo [g,f;'f,l=g,-J,= O 00 0 00 (4.21) K3100C3x00 o 00 0 00 [0 00 o 0Q Where t 6 = [A6,' A9”! ABL'] : angle variables at network buses; l v = AVrt AVH‘ AVL‘] : voltage variables at network buses; l P = [Mort APH' APL‘] : voltage-dependent active power load demand: Q [AQr‘ AQn‘ A941! :voltage-dependent reactive power load demand; Note that the submatrices in the fourth row and fourth column blocks of g, model the reactive power Jacobians associated with the generator terminal buses, while in the classical load flow model (4.18) the fourth row and fourth column blocks are associated with unregulated generator PV-buses which became PQ-buses. Therefore, the difference between the two matrices resides mainly in the fourth row and fourth columns, since the row dimension of submatrices in the fourth row of JLF associated with Q’T’Q and the column dimension of the submatrices in the fourth column of J U, associated with Vf‘l equal to the number of generators with PV to PO bus changes, but those in gy are proportional to the number of generators in the model. These are also differences in the terms Am, Am, C”, and C3,, of g), and A, , AfQ, CW and C5? of JLF .The relations between these terms as well as their relation to the equivalent algebraic Jacobian matrix J y is now discussed when the dimension of uPQ and the number of rows associated with V§Q and Q’T’Q in (4.18) are zero. Note that the submatrices A”; and A3,, of gy are given by AIK = Kipqr+41 ; A311 =K31+A3 where A, and A3 do not include the effects of the K131: aPT/ae, and K3,: aQT/ae, coupling between generator terminal and internal buses. Thus A, and A3 are identical to the matrices A, , A§Q of a classical load flow model Jacobian (4.18) since this model includes terminal buses but not internal buses. Also note that ClK = Kits-1C, ? C311 = K317+C3 where C, and C3 do not include the effects of the K #17: ap,/av, and K547: aQT/av, coupling between generator terminal and internal buses. Therefore C, and C3 are identical except in dimension to the matrices CfQ and C5]? of a classical load flow model Jacobian (4.18) since this model includes terminal buses but not internal buses. Expressions for K51, , K51, , K fat, and K3, are easily derived [52] and can be found in Appendix A. When adding [1, -g,, ] of (4.21) and g, of (4.20) to obtain Jy , it turns out that A, K and A3,, of g), are replaced [52] by A, and Ag’Q of the classical load flow Jacobian J L F which implies that the exciter and flux decay dynamics have no effect on aPT/ae, and aQT/ae, terms of J y . Furthermore, in J, the submatrices C , K and C3,; of 124 g, are replaced by Cir = Cir—C111 Cir = CsK‘Csx C3+K37 ‘Crx Matrices C 1x and C3,, are derived in [52] and are dependent on the DC gain of the excitation system KEE which makes diagonal matrices CfK and c), have very large negative elements when all exciters are active. However if the armature resistance of the generator R0 = O , then [52] C”, = Kg, thus ch = C,. This result implies that when the armature resistance of the generator is neglected, the control systems do not have any effect on the active power-angle Jacobian matrix CfK = C,. Since R, is generally small so that CL, = C, and 1:1,, is assumed to be large when exciters are operating at high gain at every generator, the conditions for static bifurcation of J y 0 = 1y. )4, 4.23(a) where n,” = (air. 11;”. uh. 115,. 115),. 115,) 4.23(b) Equation (4.20), (4.21) and (4.23) imply that A3uPT+B3uPH+C§KuQT+D3HuQH z C§KuQT z 0 4.24(a) since or, is very large due to its being linearly dependent to KEE. From equation 4.24(a) and the fact that cat, is diagonally dominant when all generators have excitation control, 4.24(a) suggests that em = 0 4.24(b) Since 4.24(b) is associated with the fourth row of J y not in J L r when all generators have excitation control and uQT = 0 eliminates the fourth column of J, , the remaining 125 equations of 4.23(a), (which do not include 4.24(a), require that J T T T T T T :3 0 LF' “pr “PH “PL “Po “QH “or. where uPQ has dimension equal to zero. Thus Propositions 4.2 implies that the classical load flow model is a bifurcation subsystem of and u = , u 0 Ay for LP bifurcation for the case where all exciters are active and operating with high gain, Ra on every generator is assumed to be zero, and own is invertible. An option has been added to some load flow algorithms to model the generator flux and exciter loop. In this case, the exciter gain on every generator does not have to be large for the load flow to be a bifurcation subsystem for the differential-algebraic power system model, since in this case J L r = J y regardless of whether the exciter is present or is disabled by field current limit controllers. The only difference in J u, = J y when exciters are present or disabled is in the terms C, x and C3,, that depend on whether excitation system gain K or = 115,511,,“ is zero or nonzero, where ng is the inverse diagonal matrix of the inverse DC gain of the exciter, K A is a diagonal matrix of the amplifier gain of each excitation system, and K o is the diagonal matrix of the sensor gain for each generator exciter. 126 4.5.4 Bifurcation in LF Model or in Generator Dynamics Recall that the conditions for g, to be a bifurcation subsystem of J, and J require that [1y —gy] - u2 = O (4.25) when J,- u, = o (4.26) From (4.21), condition (4.25) can be written as 0 4.27(a) K1141 '“PT+C1X'“QT O 4.27(b) K31 '“PT + Csx ' “QT where u, is the eigenvector associated with the zero eigenvalue, T .. T T T T T “2 " (“PT’ “PH’ “er “QT’ “on: “91.) The J y condition in (4.26) for terminal buses can be written A,uPT+B,uPH+(C,+K;¥7-C,X)uQT+D3HuQH = 0 4.28(3) A3“PT + B3uPH + (C3 + [(8, - CBX)“QT + D3HuQH = 0 4.28(b) There is very little chance that J y and J will be singular simultaneously with g, since the C3,, term in equations 4.27(b) and 4.28(b), which is dominant, occurs with different signs. On the other hand, if 11,, , up” , um and 110,, were zero for a specific bifurcating eigenvalue, then from (4.26) and (4.27) we can see that the bifurcation occurs strictly in the load buses far from the generators in the system, which indicates that J and I, can reach bifurcation when gy reaches bifurcation [52]. If 11", “PH , um. and 110,, are nonzero then g, is never singular when J and J), are singular, and singularity in J always occurs due to singularity in 1,. This bifurcation in generator dynamics that occurs whenever uPT, up” , um. and 11,,” are nonzero, will occur when load flow bifurcation occurs if all generator exciters have infinite gain, all generators have zero armature resistance, and f, is nonsingular, or if the load flow modeling option has been selected that makes J, = J ,F. 4.5.5 Conclusion (1) (2) (3) (4) The bifurcation subsystem method suggests that the classical load flow model is a bifurcation subsystem for the differential-algebraic model for all static bifurcations when all exciters are active and high gain at the point of bifurcation; R, is zero on every generator, and f, is invertible. the classical load flow model with complete DC modeling of generator exciter and flux decay effects is a bifurcation subsystem for LP bifurcations for both J -Ay = o and 11 = ,_ u 0 Ay the set of algebraic equations for a differential-algebraic model is not likely to be a bifurcation subsystem of all the differential-algebraic model whether exciters are active or disabled, unless an , u H, , uQT and uQH are zero. This result implies that bifurcation of these algebraic equations can only occur if it occurs solely in the algebraic equations of load buses in the system, which produce an algebraic bifurcation subsystem. as a result of (3), static bifurcations of the differential-algebraic model are either static bifurcations of the reduced differential model with Jacobian JJr when an , up” , uQT and 11,, H are non zero, or an algebraic bifurcation in the real and reactive power balance equations at load buses when 11,-, , up” , “QT and 11,,” are null. The algebraic bifurcation is referred to as clogging voltage instability in [l l] and is due to inability to ship reactive power to load buses. The clogging occurs due to reactive losses that consume the reactive supply. The static bifurcation of the reduced differential model is a saddle node bifurcation because the generic static bifurcation 128 is a saddle node bifurcation. The static bifurcation is most often a loss of control voltage instability that occurs due to exhaustion of reactive reserves and voltage control at generator buses since the generator terminal bus is affected by the bifurcation. There can be clogging voltage collapse problems [11] that are due to static bifurcation in the reduced differential equation model since clogging voltage collapse can exhaust reactive reserves and cause loss of voltage control, but this is by definition not the cause of the instability. The exhaustion of reactive reserves of these generators could only occur if the generator terminal buses were affected and uQTrtO. V Effect of Hard Limit Discontinuities on Stability Analysis 5.0 Introduction The Bifurcation Subsystem Method developed in Chapter 4 and the Epoch Based Trajec- tory Stability Assessment Method developed in Chapter 3 provide a framework for initiat- ing a structural stability analysis on a single—machine—infinite bus (SMIB) model that encompasses local shunt loads, with emphasis on the effects of one type of hard limit dis- continuous system transition, the disablement of excitation control. This discontinuous control action is shown to have an effect on the kind, class and the tests for bifurcation occurring on the SMIB model. The goals of this chapter are to: (a) Demonstrate two examples of using the bifurcation subsystem method to identify bifurcation subsystems in a model. The generator angle-speed dynamics is the bifur- cation subsystem for Hopf bifurcation in the regulated SMIB model, and the angle- speed-flux decay dynamics is the bifurcation subsystems for SN bifurcation, in the unregulated SMIB model, respectively. (b) Analyze the effects of hard limit discontinuities by assessing the effect of excitation control or lack of excitation control on the (i) types, (ii) classes of bifurcations pro- duced, and (iii) the stability test needed for detecting the occurrence of bifurcation. 129 130 (c) Derive expressions for the synchronizing and damping torques in a linearized power system model with no exciter, a thyristor exciter, and type 1 exciter models in terms of the model parameter K.-K,. The damping torques are derived for f<> l. (d) Derive conditions on linearized model parameters [(4 and K, for the occurrence of Hopf and on parameter K,, K, K3, K... K, and K, for Saddle Node bifurcations, and preservation of stability of the equilibrium. Such conditions are derived for the unreg- ulated model and for the regulated model (thyristor and type 1 exciter) SMIB model. The conditions on model parameters K, and [(5 set the stage for establishing diagnostic tests for specific kind and class of Hopf bifurcation in terms of active load, reactive load, shunt capacitive susceptance, and real power transfer levels in Chapter 6. The stability analysis in this chapter is based on the principle of synchronizing and damp- ing torques or a bifurcation sub-subsystem MAS + 0116‘ = Arm — AT e AT, = [(5135 + KDAS 5 ATS + ATD where AT m : incremental change in mechanical power input to the generator AT e : incremental change in the electrical power output of the generator AT s : incremental change in synchronizing torque AT D : incremental change in damping torques of the machine The use of damping and synchronizing torque model for the model shown in Figure 5.3 was developed in [55] where the generator inertial and exciter dynamics are assumed to produce equivalent damping and synchronizing torques on the inertial dynamics. This type of analysis performed in [55] was later extended in [56], in both classical papers on 131 power system dynamics. Matrices KS and K,, were derived in [52] for multi-machine power systems given that the subsystem external to the mechanical dynamics is invertible. The Ks and KD tests for saddle node and Hopf bifurcations respectively, do not necessarily imply that the bifurcation occurs in the generator mechanical dynamics, nor that the inertial dynamics are a bifurcation subsystem, and neither the test matrices nor the bifurcation subsystem indicates anything about the center manifold for the bifurcation in a non-classical power system model. The system is stable if and only if both the synchronizing and damping torques are positive. If a classical generator model is used, where no flux decay or exciter dynamics are modeled, then the damping and synchronizing torques are namely ATS = K, , the machine real power coefficient, and ATD = D, the natural damping of the machine. In the case of more complex generator models, ATS and ATD depend on the single machine infinite bus model constants Kl'KG. The form of ATS and ATD in terms of K,-K6 have been derived for nonclassical generator model with different exciter models and for the unregulated model [56]. The system Kl-K4 constants (unregulated) and Kl-K6 (regulated) are derived as part of this thesis and are expressed as functions of the system operating point, the network parameters, and the infinite bus voltage. Although the parameters or matrices ATS and ATD do not imply that the inertial dynamics contain the center manifold or even a bifurcation subsystem, they are tests for saddle node and Hopf bifurcations in the full system model utilized as effective measures of stability and security, and are for that reason used in this dissertation. 132 5.1 Linearized Multi-machine Power System Model The power system model used in this chapter is based on linearizing the differential and algebraic equations of a multi-machine power system model. The modeled system consists of n synchronous generators represented by a two-axis models, m network buses, and an infinite bus for each area in the model. The model takes the multivariable block diagram form shown in Figure 5.1, where the state variables are the electrical angle and rotor speed deviations of the n generators A5 and A11) , and the q-axis transient EMF deviation AE', . In this model, T M is a diagonal matrix of inertial time constants TM; = 2H ,, T’d is a diagonal matrix of field open circuit time constants T'dai, and D a diagonal matrix of damping coefficients. The inputs to the multi-machine model are the mechanical power deviations AP", and the excitation voltage deviations AE, of the n generators. The nxn linearization coefficients matrices K,, K2, K3, K4, K5 and K6 depend upon the network parameters, the system quiescent operating point, the infinite bus voltage. They change as the system conditions change, and therefore can decouple angle and voltage dynamics or can strongly couple voltage and angle dynamics for other operating condi- tions. Specific system operating characteristics, such as heavy local loading conditions, large transfer of active power, or a capacitive network can dramatically effect K1'K6 and the kind and class of the bifurcation produced. As the significance of these constants will be revealed throughout this chapter, it is important to explain their definitions. In the next section a linearized SMIB model will be derived, where these constants Kl-K6 are scalar. (a) K, = BPG/ (8801p : Change in electric power for a small change in rotor 411 = 15,411" angle with constant flux linkage in the direct axis °. The component of E’q,| = Eq, torque (or power) K,A51‘ is in phase with A5 and hence represents a synchronizing (b) (C) (d) (e) (f) 133 power component. Indeed K, is the synchronizing power coefficient in a classical generator model where no exciter or field circuit is present. K, = arc/(a Eq, )| , _ : Change in electric power for a small change in the d- 5: = 51° axis flux linkage at constant rotor angle 81' = 81°. K 3 = [1 +((26,,,-—x',,;)(xq;+11,,,.))/(A1')]"l :K3 is an impedance factor that takes into account the effect of the external impedance that relates the generator transient closed circuit time constant T’ (,0 to the transient open circuit time constant T’ d,. 1 K4 = E313 E'qi)/(88i)1efd : K, represents the natural coupling between i: constant AlE’,,-| and A81. It is related to the demagnetizing effect of a change in the rotor angle. K5 = (alv,,| )/(a81')llE' , 0: K5 is the change in terminal voltage V,, for a small qt = 15' qi change in rotor angle with constant flux linkage in the direct axis E' ° qi=E' qi K, = (alv,,|)/(a : Change in terminal voltage V,,- for a small change in Eq, ’1 i=5i° the d-axis flux linkage at constant rotor angle 81’ = 81°. 134 AlV,l + <—-——4 K, . Ks A8 + + _, AIE'ql _> (Ta) 1 Figure 5.1 Block Diagram of a Power System Model 5.2 Single Machine Infinite Bus Model With Local Load 5.2.1 Model Representation Considering how large and complex a general power system is, understanding the basic effects and concepts of power system stability phenomena using analytical techniques requires a simple configuration and a reasonably low-order model. A single machine con- nected to a large system through a transmission network and supplying a local constant impedance load (Z,,,) is shown in Figure 5.2 (a). This system may be simplified by using a Thevenin’s equivalent of the transmission network external to the machine and its trans- former as shown in Figure 5.2(b). Furthermore, because of the relative megawatt (MW) generator capacity and inertia of the machine that represents a generating station compared to the generator MW capacity and inertia of the rest of the large interconnected system, dynamics associated with the single machine will cause virtually no change in the voltage and frequency of a Thevenin equivalent source (V0,) used to represent the large intercon- nected system. Such a source of invariable frequency and voltage has been called an infinite 135 bus [29]. The described single machine infinite bus model has widely been used in the lit- erature as a simple yet powerful tool for understanding of power system stability of far more complex systems under small disturbances [55]. Accordingly, in order to develop an appreciation for the dynamical behavior of the system as it undergoes bifurcations, their causes and effects, it is virtually necessary to investigate such a low-order system before dealing with very large complex systems. The Thevenin electric circuit of a complex trans- mission network model is shown in Figure 5.2(a) with simplification shown in Figure 5.2(c). This single machine-infinite bus model includes (1) A generator and its transformer represented by a voltage source 5,128 behind a transient reactance x',, + x, where x7 is the transformer reactance. This generator and transformer are radially connected to a large interconnected system, via a long transmission corridor with equivalent reactance X E that is large. The large interconnected system represented by the infinite bus (i) may be a major bus or group of buses that are stiffly interconnected so that they act as a single bus and (ii) contains very large generation capacity and inertia compared to the single machine connected to this bus group via the transmission corridor reactance X E; (2) A load center at the generator terminal is represented by a constant shunt resistance R,,, in parallel with a shunt capacitance or inductance X 3,, that can include line charging of the transmission line, local shunt capacitive compensation, or inductive load. (3) X E is the equivalent transmission network inductance shown in Figure 5.2(b); (4) The equivalent transmission network resistance R E shown in Figure 5.2(b) that is neglected in Figure 5.2(c); 136 (5) The infinite bus is modeled as a voltage source with voltage 1V°°1 40°. 5.2.2 Model Implications Despite the simplicity of its configuration, the single-machine infinite bus model can be used to represent (a) A generator connected to a large system, via a long equivalent transmission corridor with equivalent reactance x); that is large. The large system may be a major bus or group of buses of a power system of very large capacity compared to the rating of the machine under consideration, or a large load center. Real power is shipped from the generator to the large load center, whereas reactive power may be shipped in either direction across the transmission network; (b) A generator serving a heavy local constant impedance load that can include line charging of the transmission line, local shunt capacitive compensation, or inductance load. With x5 large, the dynamics of the generator serving a local load center can be investigated. Effects of transfer to and from the large remote load center can also be included. 137 Loca—l load __ T (a) Vt le|40° /G\ l e V =RE+1XE E Zsh Locjload 1 (b) V, Ia 1V°°1 40° Xd'+XT ———’ ‘ XE I? (11111 ' 1 011111 E [154’ 28 a Rsh g: xsh E— —'L“— —=" Local load (0 Figure 5.2 Single machine with local load connected to an infinite bus through a transmission network: (a) General configuration of a single machine connected to a large system through transmission lines; (b) Equivalent system where the transmission network is reduced to its Thevenin’s equivalent; (c) Circuit model of (b) where the RE is neglected 138 V, I, v; lifl—M—tlxil___i lEq’28 (a) q-axis E d-axis E ' . (Xd-Xqfld qa ‘ Eq, ‘. , ‘. 1(xq+xe)lq I. 8 Ref , ‘ , (X ’d‘ Xe) Id ///// Id a rela (b) Figure 5.3. Simplified single-machine infinite bus model. (a) Equivalent Thevenin model where r, x, and v; are given by equations (5.1). (5.2) and (5.3) respectively, (b) Phasor diagram for the synchronous machine connected to an infinite bus shown in (a) 139 5.2.3 Model Equations The circuit model in Figure 5.3 (a) can obtained from Figure 5.2 (c) by first transforming the Thevenin’s equivalent XE and 1V°°1 40° to a Norton’s equivalent circuit, aggregating XE, X,h and R5,. and then transforming the resulting Norton’s equivalent to a Thevenin’s equiv- alent. This produces the circuit model shown in Figure 5.3 (a) where re and xe and V; are the equivalent thevenin’s resistance, inductance and infinite bus voltage. The parameter re, x8 and V; are expressed in terms of the original line inductance XE, and the shunt loads R,., and X,,, in equations (5.1), (5.2) and (5.3) respectively. Rsh rc = 2 (5.1) 1+ii,,,(1/x,,,+1/XE)2 Z R3,,(1/X3h-1- 1/x,.) ,, = 2 (5.2) 1+R (1/x +1/x )2 sh :11 E e = 1le a Ive|z_B (5 3) °° l+XE/Xs,,+jXE/RS,, °° ' The differential equations for the single synchronous machine are M8+08 = Pm-PG(8, IE'ql) (5.4) T‘d0(d|E'q|/dt)+|Eq| = of, (5.5) Where PG has the form 5' v; vg,2 P615454) = 1 , '1 lsin8+| I ( 1 — , 1 )sin2(8—B) (5.6) xd+x, 2 14"": xd-t-x, when r, = 0. An equation of PG is derived when r, at 0. The algebraic equations represent- ing the network shown in Figure 5.3 are 5,, = V;+rela+j(xd+x,)ld+j(xq+xe)lq (5.7) E'q = v; +rcla+j(x'd+x,,)ld+j(xq+x,)lq ld<0 (5.8) l 40 5.2.4 Simplified Linear SMIB Model A simplified linear model for the synchronous machine connected to an infinite bus through a transmission line with reactance x, shown in Figure 5.3(a) is derived, based on the phasor diagram in Figure 5.3(b). In the linearized model, the following assumptions have been made: (i) The amortisseur effects are neglected; (ii) Stator winding resistance is neglected; (iii) Balanced conditions are assumed and saturation effects are neglected; (iv) The synchronous machine can be represented (electrically) by a constant volt- age source behind its direct-axis transient reactance. Under these assumptions, the linearization of the E ', equation (5.2), the electric torque/ power equation and the terminal voltage equation will be performed in the next subsection. It should be noted that the equations describing the system are expressed in p.u. 5.2.4.1 The E’q Equation From the system algebraic equations (5.7) and (5.8), E q and E ’q are related by E“, = Eq+ j(xd"x.d)ld (5.9) Using the phasor diagram representation of (5.7)-(5.9), we express 1,, (<0) and Iq as func- tions of 15.41 and lv‘ml. Projecting E 'q and V; on the q-axis and the d-axis and keeping in mind that Id is negative yields -(x’d+x,)ld+rclq = [E'ql-lvglcosw-B) (5.10) (xq+x,)lq+reld = lVf,|sin(5-B) (5.11) |qu = [E'ql-(xd-x'dfld (5.12) Assume steady state operating condition at an operating point such that s = 6" 1, =1; 15,1 =15'.1° 1v1=1v1° 141 B = B ’d = ’4° 514 = 514° Pa = (”13° = Po(8°’ IE'41°) A180 assume a small perturbation around the operating point, such that 5=8+A6 lq B=BO+AB Id lq°+Alq |V,| = |V,|°+A|V,| lo° + No IE'ql = IE'ql° + A1511 Substituting these into equations (5.10) and (5.11), and using Taylor’s expansion method around the nominal solution as the steady operating point yields the following linear equa- tions (5.10’)-(5.12’) —(x',, + x,)AI,, + rem, = AIE’ql + |v;| sin(8° — 15°) - A8 (5.10’) reAld+(xq+x,)Alq = lV,|cos(8°—B°)~A8 (5.11’) Solving for A1,, and A1,, and writing the result in matrix form as _,, -(xq+xe) relVeoolcos(5°-B°)—(xq+xe)|Veoolsin(5°-B°) AlEiq, (5.13) re ’91Ve°°151"(6°‘13°) + (x'd+xe)|Vem|cos(5°—fi°) A“ = r} + (x, + x,)(x',, + x,) Replacing Alq from equation (5.13) into equation (5. 12’), an expression of the incremental quadrature voltageAIE (11 as solely function of AIE’ 41 and A5 is obtained Alqu = [1+(xd—x'd)(xq+x,)A-'] -A|E' «l + [(114 - .11',,)A’l IVfol ((xq + x,)sin(6° - 13°) - r,cos(6° — B°))] - A8 which can be written as l A15111 = [(3 ~A|E'q| +K4-A6 (5.14) Now, substituting into an incremental version of equation (5. 12) expressed in the s-domain, the following transfer function is found K3K, . — —3_. _____. A1541 - 1+K3T’dos AIEfdl 1+K3T’dos A5 (5.15) where 142 K, = [1+ A-1(x,,_x',,)(x,+x,)]" (5.16) K4 = A—l(xd-X'd) V" '[(xq+x,,)sin(5°-15°)—recos(8°-B°)] (5.17) 5.2.4.2 Electrical Torque Equation The electric torque T e in numerically equal to the three-phase generated electric power PG when expressed in per unit, thus T, = PG = Re(v,-1,,*) = led+ v41, p.u (5.18) where V, and V, are the projections of the generator terminal voltage V, on the q-axis and the d-axis, respectively, i.e v, = Vq-i-de (5.19) Vq and Vd can be expressed in terms of I,I and 1,, by projecting IV.) and I, on the q-axis and the d-axis as shown in the phasor diagram in Figure 5.3, we obtain V, = v; + (r, +jx,)l,, ... |vg| cos(8° - 13°) - j|v;| sin(8° -13°)+ (r, + jx,)(1,, + jld) = |vg|cos(8° — 11°) + 51,, —x,1,,—j[|v:,|sin(8°— 13°) -x,1,,-r,1,] Which may be expressed using 5.10 and 5.11 simply as v, = lay + {41,- jqu,, (5.20) Setting real and imaginary components in these two representations of V, in equation 5.19 and 5.20 equal, the following is obtained Vq = 15.41 +x'dld (5.21) Thus, = [lqu "' (xq‘x'd),d]lq (5.22) Linearizing PG as expressed in equation (5.22), we get I43 All)G = A115.4114o " (x, " x'dlA’d’qo + Al4115.111o ‘ (x, ‘ x'dlAlq’do = A|E',,|I,,° + |Eqa|°Alq — (x, -—x',,)AI,,I,,° where we used the q-axis voltage 13,, defined in Figure 5.3 (b) as with initial conditionvalues lqul° = 1qu° + (“"d’i‘qfldo Substituting for 15410 from Equation 5.12, one obtains 15 o =1Eq’ = 15410 " (x, ‘ x-d),do qal °—(xd—x'd)ld°+(xd—xq)ld° Substituting expressions from equation (5.13) into APO to finally compute APG = [1,,°+A-IIE,,,|°r,+n-l(x,-x',,)(x,,+x,)(l,,°)] . AIE'ql + A-1|v;| [|E,,|°(r,sin(6° — 3°) + (x'd + x,)cos(8° — 3°)) _(xq-x',,)lq°(r,cos(8° - 13°) - (x, + x,)sin(8° - B°))] . A8 a K, - A8 + K, - A|E,,’1 (5.23) where K1 = 4"lVilllqul°(resin(5° - 6°) + (x'd+ x.)cos(5° - B°))] + A-1|v;|1,,°[(x, + x,)(x,, - x'd)sin(8° — 3°) — r,(x,, - x',,)cos(8° — 3°)1 (5.24) K, = rc|Eqa|°A‘1+lq°[l+ A"(x,, - x'd)(x,, + xe)] (5.25) 5.2.4.3 Terminal Voltage Equation The square of the norm of V, as expressed in equation (5.19) is given by the equation 2 2 2 V, = Vq+vd (5.26) Linearizing equation (5.25) yields V v" AV =—1;-AV +é-Av (5.27) , l " IV) 1 MI .. AV d and AV ,, can be obtained from an incremental version of equation (5.21) as AV = q AIE'ql + x'dAlq 144 Finally, substituting AV d , AV 4 , and then Ala, and A1,, from (5.13) with re = 1) into equation (5.27) we obtain V,° o v ° A|V,| = “61°.qu +Xe|Vfolcos(8 —3°-) A8+,—9—V V,° -'A|E ,| v ° ' v ° ' " x" -A|E' |-—"—-—,x—--|v; |sin(8°— 3°) A8 lV,|°.x',,+x,, IV ,|° x'd+x,, 1V;1.( x" Vd°cos(8°-B°)- .xd V,°Sifl(5°-B°1)'A5 IV,|° x,,+x,, xd+xe V° x, +_‘l_. __ |V,l° x'd +x A1541 = K5 - A8+K6 - AIE' 41 (5.28) where K = _|Vf°|.-(—x‘1——V °cos(8°—B°)- x'd V°sin(8°-B°)) (529) 5 IV,|° xq+xe d x'd-t-xe q " O K, _ 1’3. x. 15.30) +|V,|°'x'd+xe Equations (5.30), (5.31) and (5.32) below constitute the basic equations for the simplified lin- ear model, which is implemented in the block diagram shown in Figure 5.4 where the dynamic characteristics of the system are expressed in terms of the K constants K.-K6 MA8+DA8 = ATm-AT (5.30) . K 3K4 APG =K,-A8+K,-A|E,,’| (5.32) AIV,| = K5-A8+K6-A|E',,| (5.33) 145 41 Ks A1V,1 { 3F K - , , ’ + + AP“, Avref 5x31321211 + K3 , NE q K2 _ . l '1' _ 1+K3T do 5 + AI, M82 + DS ux Decay Dynamics APO Rotor angle Dynamics K4 Figure 5.4 Block diagram of the simplified linear model of a synchronous machine connected to an infinite bus. through an external impedance. 146 5.3 Application Examples of the Bifurcation Subsystems Method on a SMIB Model Assume that under some system stress, no bifurcation in the model algebraic equa- tions occurs in the SMIB model shown in Figure 5.4. Then, the resulting system bifurca- tion, if one exist in some subsystem must reside in the model differential equations representing the generator angle, speed, and flux decay dynamics. Using a power system toolbox for Matlab [59] will allow us to demonstrate the existence of bifurcation sub- systems in the model differential equations. Example 1 shows that when the system in unregulated, the bifurcation subsystem for SN is the machine angle-speed-flux decay dynamics (A8, A10. AEq’). 1n Example 2, the bifurcation subsystem for Hopf is the generator angle-speed dynamics (A8, Aw), called the inertial mode [25,31]. These examples are now presented. 5.3.1 Differential Bifurcation Subsystem for SN (Unregulated Model) A linearized differential equation model of the unregulated SMIB model is obtained under heavy local loading system stress, i.e. when the shunt inductance Gs, is increases as can be seen in Appendix C. The state vector Ax , the A matrix as well as the critical eigen- value N" and its associated eigenvector are given by Ax = [A8 Ate AEq' Awkd AEd' Awkq]T - '1 0 376.9911 0 0 0 0 —0.1980 —0.0000 40.2289 -00000 -00091 0.0000 A = —0.4740 0.0000 —0.5519 0.0000 0.0000 —0.0000 —2.1235 —0.0000 30.6816 —32.2581 0.0000 0.0000 -0.1333 -00000 .00000 —0.0000 -4.0476 0.0000 _-0.0422 0.0000 0.0000 00000155923 46.3934 147 __ 0.2767 1 j8.6336_ 60.6281 _4.0431 151000,,23 “SW = .00030 2* = —0.0030 as = 0- —32 2581 0.5572 -l6.3934 0.0207 1110213- The eigenvalues of A indicate the existence of static bifurcation (SN), provided that trans- versality and degeneracy conditions for this SN bifurcation [34] hold at this point Since at this operating point the system differential equations experience bifurcation, the Bifurca- tion Subsystem Method is used to possibly find a smaller differential bifurcation sub-sub- system for SN. In order to test whether we have a bifurcation subsystem, we need to partition the square matrix A as A = [A“ A,,] where A,, is ixi, A,, is ixj, A,, is jxi, and A21 A22 A,, is jxj, i+j = 6, i = l, ..,5. The Bifurcation Subsystem Method (Proposition 4.1) says that when the full matrix A is near singularity and for all possible subsystems of order j=2,3..,n, A22 is nonsingular, the term [A ,, - A,,l - A,,] - u,* z 0 as 11 -> 11* for u,* e Null(A,,) , then the dynamics confined in Ax, = A ,, - Ax are a bifurcation subsystem for SN. In a prac- tical sense [A , ,] - u,* is never null as u approaches the bifurcation value 11* since we never find the exact value of the bifurcation parameter 11* , and thus ”511141111531 should be a measure of how close the subsystem is experiencing saddle node bifurcation. The coupling term expressed by C E 111A12A221A21] ' “* 111 is also never null at 11 close but not equal to 11* , and thus is a good indication of how little the system external to [A,,] contributes to producing the saddle node bifurcation in [A]. This is indicated by 148 R E 11[A12A21A21] ' “I 111 111A”) ‘ “1‘11 When saddle node bifurcation is close to occur in [A,,] as measured by N, then if R << 1, the saddle node bifurcation in [A,,] is producing the saddle node bifurcation in [A] with little or no coupling or help from the system external to [A,,] . Bifurcation subsystems exist if N is small and approximately equals 1* as u —) 11* , and C is very small and much less than 71* for all u —~> 11* . This second condition requires R to be much less than 1.0. The subsystem used for testing are decided by the ranking of the participation factor or right eigenvector elements except that the dynamics 8 and or are considered together since they are part of the inertial dynamics. Note that both N and C are equal to or less than 1* = —0.0030 for all subsystems of order greater or equal to third order, which implies that the bifurcation subsystem is (A8, Aw, AEq‘). Note that N = 71* and R = C/N is less than one for any subsystem of order equal to or greater than three, indicating that the third order model is the smallest bifurcation subsystem, but that the fourth and fifth order models are all bifurcation subsystems. This third order model is thus the bifurcation subsystem since it is weakly coupled (R <<1) to the dynamics of AE,‘ and A111,“, from the results in rows 4 and 5 in Table 5.1. 149 subsystem -1 A,, 141242242111,- [Arr] ' "i N C R 145] 4.8397e + 12 inf 4.8397 6+12 0 0 0 0.1237 0.9944 A5 1 0019 0.1244 [Am] 10-1237 10.12431 ' As 0 0.0019 0 -3 . 1.88731: - 04 .0754 A2,“ —0.1887 10 0:0002 000” 1 4,3000 -0.0016 :05] 0 0.0019 —0.1887 -3 0.0002 10 0.003 _ AEq 0,0000 -0.(X)16 1.8873e 04 0.0626 A111,“, 0.0000 .0.0017 - A8 ’ ‘ ' 1 Am $0043 0001: 0.003 ' -13 _ _ AE‘q 0.0218 ,0 -0.0016 4.47e 14 1.49e 14 AW)“, —0.2474 .0.0017 , A54 1-0-3719. 1-0-0001. Table 5.1 Computational results for identifying the bifurcation subsystem for SN bifurcation in the SMIB model. 150 This bifurcation subsystem is quite obvious from Figure 5.5 since three integral blocks have infinite gain when A = s = 0 which implies that the states A8, A0), AE,‘ are con- nected by infinite gain blocks in a loop when saddle node bifurcation occurs. This explains why the flux decay and mechanical dynamics have both been argued to be associated with saddle node bifurcation when in fact both are part of the bifurcation subsystem. It should be noted that the saddle node bifurcation could be due to inertial dynamics or the flux decay dynamics experiencing bifurcation and yet the bifurcation subsystem (A8, Am. AEq‘) would still be from Figure 5.5 because there is a loop with infinite gains relating (A8, Am, AEq’). In order to see whether the information obtained from the participation matrix indi- cates similar conclusions, the participation factor have been computed for this case and summarized in Table 5.2. This information indicates the only state with significant partici- pation in this SN bifurcation is the AE,’, which is one of the states in the bifurcation sub- system. The eigenvector u* is also a poor prediction of the bifurcation subsystem since (A8, A0). AEq‘) are predicted to be associated with the dynamics experiencing bifurcation. It should be noted that the bifurcation subsystem does not necessarily contain the center manifold dynamics since no formal theory has been presented that suggests so for this saddle node bifurcation. The fact that R<<1 as u -> 11* so that the coupling of the exter- nal system is not affecting the conditions for bifurcation [A,, —A,, - A,, -A,,] - u,* = [A,,] . u,* long before bifurcation occurs in the full model matrix [A I and in the bifurcation subsystem matrix [A ,,] at the value of 11* , suggests but does not prove that (i) the center manifold of the full system is contained in the bifurcation subsystem, or (ii) that the bifurcation subsystem necessarily contains any portion of the center manifold dynamics. The proof is beyond the scope of this thesis. 151 Several factors suggest the bifurcation subsystem contains the center manifold beyond the discussion above are: (i) (ii) (iii) (W) (V) (Vi) a root locus of the unregulated SMIB model shows that the eigenvalue associated with AEq’ is the one approaching zero; the participation factor indicates AEq‘ is by ar the largest element; the test matrix [52] Ks for saddle node bifurcation in the full system to be derived for the unregulated model is K, - K,K,K4 < 0 and thus depends on the parameters of the bifurcating subsystem; . a test matrix T developed by [51] for an equivalent flux decay model AE’, = TAE’, is T = (—1/T’,,,,)[1/K3 — K,K,/K,] for the SMIB model shown in Figure 5.5, is singular when the full system is experiencing saddle model bifurcation. This test matrix also depends on all the parameters of the bifurcating subsystem. This test matrix if it exists assures that loss of causality has not occurred and that saddle node bifurcation has not occurred [51] in mechanical dynamics K, at 0; the test matrices T and KS for saddle node bifurcation in the full system, do not theo- retically guarantee that bifurcation occurs in flux decay or inertial dynamics respec- tively as their derivations might suggest, but solely that a saddle node bifurcation in the full system exists. This confirms the above result that as long as K, at 0 , both T and Ks and are singular at the same point, i.e. the point of bifurcation. This would suggest that the bifurcation subsystem (A8, Am. AEq’) containing all the parameters in both T and Ks and no others, contain the center manifold of the full system model; the fact that the flux decay eigenvalue in a root locus is the one approaching zero [59] as the term K,K4 increases and that the participation factor so heavily weighs AEq' 152 might suggest the center manifold is just the AE,‘ dynamics, but the fact that both T and K3 are both singular at saddle node bifurcation, and both depend on all the parameters of the bifurcation subsystem and no others, the center manifold would appear to lie within the bifurcation subsystem. Without determining the center manifold, it is impossible to say whether it is tangent to the eigenvector of AE,’ in the (A8, Am, AEq’) bifurcation subsystem subspace, or even whether it lies outside the bifurcation subsystem. It is clear that the test matrices T and K5 are diagnostic in saying the center manifold lies solely in the inertial or flux decay dynam- ics since both are zero at the bifurcation. The above results (i-vi) suggest that the bifurca- tion subsystem may well contain the center manifold but no proof exists that it does. AIE’ql 153 A8 V (a) Figure 5.5 Block diagram of the a linearized machine model showing the bifurcation subsystem and the test matrix T. State A8 A11) 45,, A‘de A54, Participation 0.(X)00 0.0000 0.2088 0.0000 0.0000 factor Table 5.2 Participation factor results for the SN bifurcation in the SMIB model 154 5.3.2 Differential Bifurcation subsystem for Hopf (Regulated Model) A linearized differential model of a regulated SMIB model is obtained when the transmission network inductance x, or the power transfer to the infinite bus has been grad- ually increased. The input parameters for this bifurcations are given in Appendix E. The state vector Ax , the A matrix as well as its eigenvalues and the critical one 71* , are given by T Ax = [A6 A0.) A5," Awkd AEd' Aqu AVA: AVR] ' 0 376.9911 0 0 0 0 0 6 4.1395 0 -0.0868 -0.0781 —0.0072 —0.0318 0 0 4.0542 0.0000 -2.0818 2.0781 —0.0240 -0.1056 0.0000 1.9373 A = -2.9066 0.0000 30.7985 --33.5716 0.0000 0.0000 0.0000 0.0000 4.0934 0.0000 0.0000 0.0000 -8.3982 6.5047 0.0000 0.0000 —0.9893 0.0000 0.0000 0.0000 15.6517 49.6571 0.0000 0.0000 3.2771 0.0000 4.5117 -8.5604 3.0798 1.5531 4.1499 0.0000 47.0239 0.0000 -163.5591 486.1236 56.3153 247.8250 14.9613 —l72.986Q 471.59 42.99 eig(A) = ‘2559 N" = 0.09: j7.16 0.09 :1: j7.16 -394 __ 1.46 :1; j0.4l_ The eigenvalues of A indicate the existence of a swing mode of frequency to = 7.16 rad/s Since at this operating point the system differential equations are at bifurcation, it is possi- ble to find a difi’erential bifurcation sub-subsystem for Hopf. Proposition 4.3 implies for- mulating the matrix 155 — 171.59 - 17.16 0.3285 + )08108 — 32.99 _ j7.16 — 0.0153 + j0.0064 — 25.59 — )7. 16 0.0424 + 10.0114 6,3,3, = 0.09 ,. = 0,09 ,,., = - 0.0022 - )0.0591 0.09 _ )1432 — 0.0349 + )00082 - 3.94 - j7.16 - 0.0499 - j0.0160 — 1.46 + j6.75 0.3991 + j0.0264 _ — 1.46 - )6.57_ _ 0.3991 + j0.0264_ In order to test whether a bifurcation subsystem exists, the B matrix is partitioned as 3 = [8“ 8‘2] where 3,, is ixi, 3,, is ixj, 3,, iiji, and 3,, iijj, i+j= 8, i = 1, ..,7, and 821 322 form column vectors a, = [u,...u,-], i=1....7. Note that N = "[B,,]-u,*||, ; Cs||{B,,B,,B,,].u*,|| and R = C/N are defined as for saddle node bifurcation for all subsystems where 3,, is nonsingular. N must be less than 2* and C must be << N“ as u —> 11* for a bifurcation to exist. The values of N and C are less than 21* = 0.09 and R<O {fandonlyif K4K2>0 . . K2K3K4 1+0) T’d 5.4.1.2 Thyristor-Type Excitation system This is a simplified thyristor type excitation system that includes only the elements that are considered necessary for representing a specific excitation system, namely the terminal voltage transducer time constant TR with gain K R , the exciter with gain K E and time con- stant TE. The transfer function is given by G (S) ___ IAEfdl = "KEKR t’ A|Vt| (1+sTE)(l +sTR) (5.42) For typical values for the terminal voltage transducer time constant TR and for the gains and time constants in fast exciters, acceptable bounds for TR, the exciter gain TE and exciter time constant TR satisfy KE»1; TE<<1 and TR<<1 Substituting G,(s) in (5.42) into the general transfer function in equation (5.38), and col- lecting the powers of s in the numerator and the denominator, we obtain AT, 00 + ass + azs2 __ = K1 _ (5.43) A5 b0~1»bls+bzsz44:353 Where a1 = K2K3K4(TE+TR) 02 = K2K3K4TETR b0 = 1+K3K6KEKR 5K3K6KEKR b1: T’d+TE+TR ()2 = TETR+T’d[TE+TR] b3 = T’dTETR Evaluating ATS/A6 at s = j0), we get AT 0 -a (oz-1'0 0) _55 =Kl— "22, 1‘ 2 (5.44) A ,,]-a, bo—bza) +1(b10)—b30)) Therefore ATS and ATD are given by: [an - azwz] [b0 — b21021 + alwlblw - b3w3] AT = K ~A5- S ‘ [ha-b2w2]2+[blm-b3w3]2 A8 (5.45a) [aO—aZwZHbl - 17302] —al[b0— 12202] AT = D [ho-b20212 + [blm-b3m312 . (ma) (5.451» 163 (a) Frequency Range (0)« 1) Setting the of terms to zero in the numerator and the denominator of K o and K 5 yields, ATD aobl -alb0 D: — _ wA503=0 bzo (K2K3K4 + K EKRK2K3K5)(T’d + T5 + TR) — K2K3K 4(15 + TR)( 1 + K3K6KEK (1 + K3K6KEKR)2 R) or approximately (T’d+ TE+ TR) K (TE+ TR) KDEK5K2 — 4 2 (5.4621) . K3K125KEKR K6KEKR by assummg [(5 >> 1 Similarly, the synchronizing torque coefficient is K _ K aobo _ K2K3K4 + KEKRK2K3K5 s - 1 — - 1 — b2" 1+K3K6KEKR K 2 K 2 or, K s s K 1 K (5.46b) —K .____ ._ 4 K6KEKR 5 [(6 Therefore positive damping and synchronizing torques are obtained for the following con- ditions: T£+TR K >0 ifwehave K >KKK -,——— D 5 4 3 6 Td+TE+TR K K K K >0 ifwehave K >K ---3+—“.-—2 5 1 5 K6 KEKR [(6 (b) Frequency Range (0)»1) In this range, the terms involving the highest power of 0) dominate. Hence, keeping only these terms from equations S.45(a—b), approximate expressions for K 1) and KS may be obtained as, - azb3 1 K0 - W =KZK4‘m (5.473) Similarly, the synchronizing torque coefficient may obtained by retaining only the highest powers of 0) in the numerator and the denominator 164 s - 1‘ [232036 " 1 13ng K K K or, K5 5 Kl——2——3—5§ (5.4713) “32(7’40) 5.4.1.3 IEEE dc Typel Excitation System The gain of the voltage transducer block is usually close to unity, it time constant is small, in the 0.2 to 0.4 range. The amplifier gain KA is typically on the range 25 to 400, its time constant TA is typically in the range 0.02 to 0.4 sec. The power system stabilizer parameter values are KP in the range 0.02 to 0.1, and T; in the range 0.35 to 2.2 sec. The block dia- gram of this system has been shown in Figure 2.2, from which the transfer function is given by 25995400 0.0251As0.1 0585513 —l.0$KESI.0 0.551‘5310 550.05 The transfer function for this excitation system is given by -KAKR(1+ sTF) 09(5) = [(1+sTA)(1+sTF)(55+KE+STE)+‘KAKF]'(l+STR) (5.50) Substituting Ge(s) into the transfer function in equation (5.38) and collecting terms of pow- ers of “s” in the numerator and denominator we obtain AT! r10+ass-1-0252+03.33+a4s4 = - 5.51) A5 'bbb2b3b4b5 ( 0+ l5+2s-1-3.9-1-454-55 Where a0 = MK2K3K4+K2K3K5KAKR 5K2K3K5KAKR since KAEM [65 a4 = K2K3K4TATETFTR a T04 b0 = M + K3K6KAKR s K3K6KAKR since K3K6KAKR » M 5K3K6KAKRTF+Tbl b = T 4 b4 b5 = T’dTATETFTR = Tb5 M = SE+KE = K3T’do '5 KAKF+TE+M(TA+TF+TR) T T T = KAKFTR + (TE-1» MTR)(TA + TF) + MTATF + TETR T TATETF+TETR(TA+ TF)+MTATFTR T - KAKF+K3K6KAKRTF+M[T’d+TA+TF+TR]+TE @- I T13, = TE[T+TR]+M[TATF+T’4TR]+[TE+MT+MTR][TA+TF] Tb? = [T’dTE+MT’dTR+TETR][TA+TF] +TATF[TE+MT’d+MTR] + T’dTRTE Now, evaluating AT/AS at s = j0), we get ATc A0 ~1-jAl = K -..—___... (5.52) A8 ”1.0) ‘ 30"131 Where 2 4 A0 = ao-azw +040) Bo = bo- b20)2 + b41104 B1: b10)-b30)3+b50)5 The synchronizing and the damping torques can finally be deduced. 148 +A£i ATS = [Kl -°—‘2%]-A5 (5.5323) BO+B1 A03, - .4130 ATD = —-2——2—— - (0)A8) (5.5313) BO+Bl 166 (a) Low" Frequency Range ( 0) « 1) In this range of frequency, we can neglect the terms in ATS and ATD of higher exponents. “obt — (11 b b 0 0 (5.54) Ill K D Evaluating the numerator yields to Replacing in equation (5.54) gives an approximate expression for K D, for 0)<<1. K K T T Tb‘ K2 Ta‘ 555 D ‘ 5 K F’ R+KAKRK3K6 ‘ KGKAKR ("20 Similarly, 5‘ 1 132., ‘ 1 K3K6KAKR K or, KS :—_-K,—K5 - 1?? (5.5513) 6 Therefore, Ta K >0 ifwehave K5>K4K1K6- ‘ D - T +K K K K T —T 131 A R 3 6( F R) KS>0 ifwehave K1>K5-,—(3 6 (b) Frequency Range (0) »1) TaTb Ta KKK _ 4 5 _ 4 _ 2 3 4 K0... 2 2 ..T 2 .. 2T, (5.56a) 7’wa b5“) (1) d Similarly, the synchronizing coefficient may obtained by retaining only the highest powers of 0) in the numerator and the denominator 4 (azbz—alb3)m _ K (azbZ-alb3) ~ K2K3K4 1‘ 2 6 " 1' 2 2 = 1' 2 , 2 b30) b30) 0) Ta, KS=K (5.56b) 167 5.4.2 Bifurcation Tests for the SMIB Given that the system K—constants K2, K3 and K6 are positive, bifurcation tests for SN and Hopf Bifurcations in both the regulated and unregulated SMIB models have been derived in the previous section and are stated in the form of two propositions and two corollaries, followed by their proofs. Proposition 5.1 (Unregulated system) (i) A sufficient (and necessary) condition for a Hopf bifurcation to occur in the unregu- lated SMIB model is that [(4 < 0. (ii) A sufficient (and necessary condition) for a saddle node bifurcation to occur prior to Hopf bifurcation in the unregulated SMIB model is that K 1 - K2K3K4 < 0 when [(4 > 0 Proof: (i) Since all the terms in the expression of the damping coefiicient for the unregulated SMIB given in equation (5 .4 I a) are positive except K4, then its sign depends on the sign of the parameter K.,, i .e. KD<0 when K, < 0. (ii) Similarly, the expression for the synchronizing torque coefficient K s in equation (5 .4 I b) indicates that it is negative if and only if K 1 - K2K3K4 < 0 whenmw=0 as is the case for saddle node bifurcation. The condition K4 > 0 assures that Hopf bifurcation has not occurred. Corollary 5.1 (Unregulated System) The unregulated system is stable if (and only if) the constant [(4 lies in the range 0 0,whichcanberewrittenas 1+0) T’d K 1 2 2 K < (1+0)T’) 4 K2K3 d Since stability is assured only when both synchronizing and damping torque coefficients are positive, this condition combined with K4>0 prove the claim of the corollary. Proposition 5.2 (Regulated system, Inter-area mode) When f<<1 and K4 >0 under voltage regulation, (i) A necessary condition foraHopfbifur- cation to occur in the SMIB model is that K5 < 0 and (ii) A sufiicient condition for a saddle node bifurcation to occur in the SMIB model is that K ,, K 2, K 5 and KO satisfy K 5" 2 l - K6 K 5" 2 K 4" 2 1-.-KT — K 3" 6K R < 0 for Type I exciter < 0 for Thyristor exciter Proof: (i) Equation 5 .46(a) and 5 .55 (a) indicated that sufficient damping forces exist in the SMIB systemif (T’d+TE+TR) (TE+TR) _ KDEKSKZ 2 - 4 2m— 0 for Thynstor type K3K6KEKR 6 E R K K K2[T T Tb‘ ] K2 Ta‘ > 0 f T I 't s -— - +—— — -—— or ype excz er 0 5 K6 F R KAKRK3K6 4 KéKAKR which can be rewritten as -1 K , and K > K - 1 T -T +——-‘—— 5 (Td+TE+TR) 5 4 KAKR( F R KAKRK3K6] for a thyristor type exciter and type 1 exciter, respectively. By observing the terms on the right hand side of both inequalities of K 5, we can conclude that they are generally small and either positive or negative depending on the sign of K4. Thus, K 5 is my??- ciently positive, no Hopf bifurcation occur; and when K 5 is sufiiciently negative Hopf 169 bifurcation occurs. (ii) Similarly. equation 5 .46(b) and 5 .55 (b) indicated that sufficient synchronizing forces are positive in the SMIB system if : _ . ___— - - — > 5 _. . .— KS_Kl K4 K6KEKR K5 K6 0 and KS K1 K5 K6 > 0 for a thyristor type exciter and type I exciter, respectively. If K 5 is sufi'iciently positive, large and increasing, saddle node bifurcation occurs. Corollary 5.2 (Regulated System) The SMIB model is stable for f< 0 8— 1 4 KGKEKR 5 K6 These two conditions may be rewritten respectively as TE+TR K1K6 K4 6'T",4~T,:.+TR K 2 _ K 5" R (ii) For a Type I exciter, equations 5 .55 (a) and (b) indicate that sufi‘icient damping and K5 > K4K3K and K 5 synchronizing forces exist in the system if K :: K ...?[7' ..7‘ +._n,l__]-g.153i >0 0 5 K6 F R KAKRK3K6 4 K6KAKR K, KsK—K-—" >0 S 15K6 170 These conditions may be rewritten respectively as a < +KAKRK3K6(TF- TR) 5 K2 K5 > K4K3K6 - Tb] Since stability is assured only when both synchronizing and damping torque coefii- cients are positive, conditions for positive damping are combined with conditions for positive synchronizing for the thyristor type exciter and for the Type I exciter, which leads to the desired result of the corollary. If K 5 is below the lower limit or K 5 exceeds the upper limit Hopf bifurcation occurs. 5.4.3 Conclusions It is our desire to hopefully find a subsystem of the full model which experiences, produces, and causes the same bifurcation in the full system model. Determining the center manifold dynamics requires finding a nonlinear transformation that is often given in terms of a Tay- ' lor series and may be infinitesimally linked to all the states of the model even if it virtually lies in a subsystem of the full model which experiences, produces, and causes the bifurca- tion. Determining such a subsystem if one exists is essential to diagnosing the location, class, cause, and cure for a bifurcation, and that is important from an engineering sense. Determining the center manifold is thus not what is desired, but a subsystem of the full model which experiences, produces, and causes the bifurcation, and which is normally of larger dimension than that of the center manifold. One may be able to diagnose the loca- tion, class, cause and: cure without resorting to finding the center manifold transformation, or its approximation. The subsystem which experiences bifurcation may or may not produce bifurcation in the full system since the structure or parameters external to the subsystem may help pro- 171 duce the bifurcation both in the full system and subsystem. This is especially true if the geometric coupling of the bifurcating subsystem and the external system in the direction of the bifurcating eigenvector is not zero but approaches zero. Thus the subsystem that pro- duces the bifurcation experienced by the subsystem and the full system may be of larger dimension than the subsystem experiencing bifurcation. The subsystem causing bifurca- tion may be still larger than the subsystem experiencing or producing the bifurcation in the full system model and the subsystem. The subsystem that causes bifurcation must include all model parameters, controls, and discontinuous changes that cause the bifurcation to develop. The discontinuous changes include under-load tap changes, switchable shunt capacitors switching changes, equipment outages, over-excitation limit relays,...etc. The controls include maximum excitation limiter controls, under-load tap changer controls, switchable shunt capacitor controls, and excitation system controls. The results of this chapter along with chapter 4 imply the bifurcation subsystem method may provide such a tool for identifying the subsystem that experiences, produces, and causes a specific kind and class of bifurcation, since (a) the theoretical conditions for existence of bifurcation subsystem was applied to show that they exist and can be determined. The bifurcation subsystem showed some prom- ising results in identifying the subsystem that experiences the bifurcation in the full model; (b) the sensitivity matrices K5, KD and T were derived in two recent Ph.D dissertations with the purpose of identifying the subsystem that causes and produces the bifurca- tions, and which may also contain the center manifold. It has been shown that the test matrices Ks, and T are not necessarily associated with the subsystem that experi- (C) (d) 172 ences, produces, or causes the bifurcation, but the parameters of the system they depend on can help identifying the bifurcation subsystem; The results on Hopf bifurcation indicate that the bifurcation subsystem obtained (Example 2) which experiences the bifurcation of the full system model may be too small to track what produces and even less capable of tracking what causes and cures that bifurcation, since the test matrices suggest thath depends on K, and K5 which lie outside the bifurcation subsystem (A5, A03). This result may be true since the (A8, A00) dynamics are experiencing bifurcation for all u , and this may imply that the bifurcation subsystem producing and causing the bifurcation should be identified as u —> m , and not at the point of bifurcation, so that the subsystem captures the causes and cure for that bifurcation; the determination of what [(5, KD and T depend on for a generator with exciter, a thy- ristor based exciter, and a dc] model exciter, gives implicit information for determin- ing the subsystem which experiences the bifurcation in each case. VI A Diagnosis for SN and Hopf Bifurcations in a Discontinuous SMIB Model 6.0 Introduction The goal of this chapter is to perform a diagnostic study to assess the structural causes, operational stresses, and stabilizing requirements for the occurrence of SN and inter-area modes of oscillations that occur in the generator inertial and flux decay dynamics of a sin- gle-machine-infinite-bus model (S MIB). The tool for this study is based on the bifurcation tests derived and concluded in propositions 5.1-2 and corollaries 5.1-2 of Chapter 5. The diagnostic study is performed by formulating three study models: (1) a power transfer model where the effects of both the network inductance xe and real and reactive power transfer are investigated: (2) a local loading model where the eflect of local loading on system stability may be investigated; (3) a power transfer-local loading model where the effect of local load model and transfer can be investigated. The SMIB local load may be capacitive, inductive or some combination. The generator exciter may be enabled (regu- lated case) or disabled (unregulated case). The objective of the analysis in this chapter is to establish that (i) increased local load causes quite different bifurcation phenomena in the regulated and unregulated load center models and (ii) that increased power transfer causes quite different bifurcations in 173 174 the regulated and unregulated power transfer model: When the exciter is enabled, the regulated model is vulnerable to Hopf bifurcation for increased power transfer, whereas when the exciter is disabled, the system is vulnerable to saddle node bifurcation in the flux decay dynamics as power transfers increase. The analysis also shows that when the local loading in light-real, heavy-capacitive, the regulated local load center model is vulnerable to Hopf bifurcation, and the unregulated local load center model is vulnerable to saddle node bifurcation. The combination local load center/transfer model is then used to validate, extend, or possibly disagree with the diagnostic conclusions obtained from the local load center or the transfer models. 6.1 Diagnostic Models Formulation Based on operational and structural parameter variations of the SMIB, three operationally and structurally different models may be obtained. This kind of diagnostic modeling is per- missible by the SMIB model since this model represents a generator connected to a large system, via a long equivalent transmission corridor with equivalent reactance x, The large system may be (a) a major bus or group of buses of a power system of very large capacity compared to the rating of the machine under consideration, or (b) a large load center. Such a structure allows formulating three study models as follows. 6.1.] Power Transfer Model (Ra: X“: infinity) In a power transfer model, no local load or equivalent line shunts are modeled, so that the generated real and reactive power are solely used to satisfy the power transfer requirements to the ‘large’ system at the side of the infinite bus. Real power is shipped from the generator to the infinite bus, whereas reactive power may be shipped in either direction across the 175 transmission network. The effects of the power transfer as well as the network inductance may then be investigated. In the power transfer model shown in Figure 6.1, the system steady state conditions are computed for given values of P” , Q“, and V,° = 1.010 p.u. V,, = 1.040 9 Vt XT+Xd I Xe mo 1 Tom—fl 15'qu _~__, Q... Figure 6.1 Power transfer model 6.1.2 Local Load Model A local loading model is obtained when a shunt constant load is connected to the generator terminals while the power transfer requirements to the infinite bus are negligible. The local constant impedance load can include line charging of the transmission line, local shunt capacitive compensation, or inductance load. “fith xe large, the dynamics of the generator serving a local load center can be investigated. This results in a model where stability problems are solely associated with the local load center. Note that since the power absorbed by the local load in proportional to VZ/Zsh, the local load is capacitive when X,h<0 and the reactive absorbed is proportional to l/Xm, it is inductive when X,h>0 and the load is proportional to l/Xsh. Similarly, the real shunt load is proportional to 1/R,,,. The local load is shown in Figure 6.2 176 V“, = LOZO V XT+X‘d t 15000 ‘ 015615 E E IqZS Rsh xsh Xe large LoEai-‘Ead Figure 6.2 Local loading model 6.1.3 Combination Local Load/Power Transfer Model Under normal operation, power system models serve both local loading and power transfer requirement patterns, and therefore a diagnostic study of this model is necessary. The more realistic SMIB model consists of real local loading in combination with power transfer model. 6.2 Diagnostic Study From Propositions 5.1 and 5.2, it is clear that since the system constants K2, K3 and K6 are positive for a generator mode of operation of the machine, the occurrence of Hopf bifurca- tion in the unregulated and regulated SMIB model depends on the parameters K, and K5 which approximate the sign change behavior of the test matrix KD, respectively. The oper- ating condit: ‘218 where K5<0 cause undamped oscillations in the generator inertial dynam- ics when the exciter is active. The operating conditions where K4 O and Q” 2 O or Q” s 0) (b) Local loading: inductive or capacitive shunt X3). > 0 or Xsh < 0 in parallel with a real load (R,.,) (c) Network inductance (x,). In this diagnostic study, the real power transfer P” is varied between 0.1 and 1.0 p.u, (Pa = 0.1, 0.2, 1 p.u.) , and the reactive power transfer to the infinite bus Qm is either pos- itive (Q,, = 0.0.1, 0.2, 1 p.u.) indicating reactive power absorption by the infinite bus, or negative (Q,o = —1.0, —0.9, —0.8, 0 p.u.) , indicating reactive power supply by the infinite bus. Power transfer levels of $1.0 indicate maximum system rating since for [le = 1 and Poo 1.0 x s1n(zV‘ AV”) — 2.S|Vt|srn(zV‘ AV”) —|Vi|2+|Vm||VAcos(zV,-ZV”) xe 1.0 = (Q 8 ll = 2,5(_|v°°|2 + IVNHVJ cos(zv, — eve,» 178 so that the system is operating at [Vt—AV” = 15.1° lVll = 1.45 p.u. This level of terminal voltage indicates heavy reactive supply to the infinite bus with large reactive losses sufficient to bring on instability in the network equations let alone the SMIB dynamics. Similarly, when P0:, = 1.0 and Qm = —1.0 the system operates at: 4v,_zv,, = 33.5° [th = 0.72 p.u. and this level of angle difference indicates heavy real and reactive power loads at the infi- nite bus and very large 12X losses sufficient to bring on instability in the SMIB dynamics. Thus P” = 1.0 and Q” = :10 indicate the system is heavily stressed and it is operating at marginal values of voltages and angle differences. Similarly, The chosen ranges of the local load R5,, Xsh, and x, are such that 1 s l/Rsh s 6 , and 1 s l/Xsh s 6 result in acceptable system operating conditions and viable system voltages and angles required for normal operation: When X,.,= 0.17, the resulting voltage at the terminal bus shown in Figure 6.1 is approximately V t = 0.45 , since x’d+xT = 0.2 p.u.. Operating: below this voltage level, is considered outside the system capabilities [59]. However parameter studies in literature [55,56] usually consider the entire range of RM, X3}: = 0.1.1.0 p.u. i.e, ls l/Rshs 10, and Is l/Xshs 10, as it is the approach in this study. For this kind of study, when solvability was not possible, the cases were skipped as will be observed in the tabulated results and the associated figures. The results of this diagnostic study are summarized in tables for operating conditions resulting in bifurcations, and tables for operating conditions where no bifurcations occurred. The tables indicate the test parameter, the kind of bifurcation tested, the network 179 and the power transfer operating conditions used. The arrows in the tables point to the range where the test parameter was more sensitive, i.e was more negative. 6.2.1 Power 'ITansfer Model 6.2.1.1 Effect of P,,,, and Q,,., The effect of P” and Q,° on the parameters [(4, K, as well as K, -K2K3K4 and K, —K2K5/K6 is simulated in Figures El.1-E1.4 in Appendix El, while the bifurcation results and the associated figures are summarized in Table 5.1. In these figures, the param- eters K4, K5 as well as K, — Ksz/K,5 and K, -K2K3K4 are plotted versus sequential 0.1 p.u step increases in real power transfer P, , and for increases in levels of reactive power trans- fer Q“, = O, 0.1, 0.2, ..., 1 p.u. or Q” = —l, —O.9, —0.8, 0 p.u.. The simulation results on bifurcations of this study are summarized in Table 6.1(a) for cases resulting in bifurcations, and in Table 6.1(b) for non-bifurcation cases. A careful review of these tables indicates ranges of P” , Q” 2 0 , Q” s 0 when oscillatory and non-oscillatory instability occurs, the points (pointed by arrows) where the measure of instability of a particular kind (K4, K5, K, -K2K5/K6, K, —K2K3K4) is most negative over the range of P,° and Q”, the kind of bifurcation (SN or Hopf) that occurs, the values of (P9,, , Q”) where the bifurcation occurs (in general at points opposite to where the arrows point to) and whether the instability occurs in the regulated or unregulated model. The tables are so complete so that the discus- sion that follows may be superfluous. The results indicate that Hopf bifurcation occurs in the regulated model but not in the unregulated model, and saddle node occurs in the unreg- ulated model but not in the regulated model. Namely, when the model is regulated, Figure El.1(a) shows that (i) as the real power 180 transfer P“ increases and the reactive power is such that 0 s Qm<0.3 , K5 is positive and decreasing; (ii) for 3 3 Q0, < 0.7 , K 5 = a|V,|/E)8 decreases and becomes negative as P“, increases; and (iii) for higher levels of reactive power transfer, i.e when Q” > 0.7 pu , K5<0. This implies that the regulated SMIB model is subject to undamped oscillations in the iner- tial dynamics when the generator terminal voltage begins to decrease for increases in 8 and occurs for the range of values (Po, , Q”) shown in Table 6.1. The fact that Hopf bifurcation occurs when K 5 = Blvd/35 is negative for a regulated system due to large P” and Q” indi- cates that Hopf bifurcation in the inertial dynamics is due to voltage if not voltage instabil- ity. . Diagnosis for the causes of undamped oscillations in the system is also performed when the reactive power is withdrawn from the infinite bus, i.e. real and reactive power transfer are in opposite directions. From Figure El.1(b), we see that (i) when Q” = —0.3, 0 p.u. K5>0 and (ii) when Qm = —1.0, ...,—0.5 p.u. K5<0. An oscillatory mode in the regulated SMIB model may then be brought on by negative reactive power transfer, i.e generator reactive absorption. Power transfer also affected the synchronizing term K, — K2K3K4 for the unregulated machine, where it became negative for heavy-real and heavy-reactive power transfer to the infinite bus, as indicated by Table 6.1(a). Note that the non-bifurcation cases are summarized in Table 6.1(b). Real and reactive power transfer to the infinite bus has no effect on the parameter K4, since K, is positive and increasing for the whole range of (Po, ,Q” ). If no local load is present at the generator side, K4-associated oscillatory stability problems in the unregulated model do not seem to appear. The synchronizing torques for the regulated model are not affected by the transfer 181 of real/reactive power, whereas when the system is heavily stressed by real/reactive power supply to the infinite bus, the synchronizing coefficient in the unregulated model is driven to zero, indicating the occurrence of saddle bifurcation. The bifurcation subsystem of (A5. A0). AEq‘) explains why heavy real and reactive power transfer cause saddle node bifurcation because K, decreases with real power transfer and K2 and K4 coupling between inertial (real power) and flux decay (reactive power) increases with real and reactive power transfer causing K, - K2K3K4 to approach zero. The diagnosis study will now focus on the effect of the transmission network inductance. 6.2.1.2 Effect of transmission network inductance x,, The effect of the transmission network inductance is simulated when P” = 0.8 p.u. and Q” = 0.5 p.u. in Figure 6.3. As the network inductance increases, the constant K, decreases, and finally changes sign at xe=0.4 p.u. The system damping coefficient also decreases and changes sign, indicating that the system damping forces are insufficient to damp the ongoing oscillatory modes. This implies that long transmission corridors, cause undamped oscillation in the regulated SMIB model. The inter-area oscillations often occur in a multi-machine model between two groups of machines connected by a long line carry- ing significant transfer. Excitation Bifurcation Network Power Transfer Figure System Reg/Un Kind Test x. 11R... 0x... P... G... t R H K5<0 0.4 0 0 .—l‘—. ‘ El.1(a) 0.1 1 0.3 1 R H K,<0 0.4 0 0 ,.L‘___. ,1— El.1(b) (r1 1 .1 .05 U SN K. -K,K,K,0 0.4 0 0 57—1 .'——, o 1 131.2 U SN K| “K2K3K4>0 0.4 0 O (ST—‘1 .1 0 E1.4 R SN K. -K2K,/K,>0 0.4 0 0 6.1—___) fro—1 51.3 Table 6.1(b) Power transfer model non-bifurcation results 04 1 . T r 1 1 I ! . : : : : Pint-0.8 f z 3 =. 2 z taint-0.5 2 3 005 ........ .. . ... ......... —0.05 -O.‘l —O.15 .1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 xe (DU) Figure 6.3 Power Transfer Model: Effect of transmission inductance on the parameter K5 and the damping torque Kd. 183 6.2.2 Local Load Model The effect of the local constant impedance load parameters X,,, and RS, on the occur- rence of system bifurcations is simulated in Figures E2. l—E2.4. The bifurcation results, the associated operating conditions and figures are summarized in Tables 6.2(a-b). In the local load model, Rsh is a positive parameter between 0.1 3 R3,, 51 and Xsh may be negative or positive, depending on whether it is capacitive or inductive, respectively. Note that x3 = 0.9 which effectively isolates the local load bus from the infinite bus. In this case, no transfer of power to or from the infinite bus is permitted. Therefore, three cases are considered: the effect of (i) real-inductive load and (ii) real-capacitive load. The simulation results on bifurcations of this study are summarized in Table 6.2(a) for cases resulting in bifurcations, and in Table 6.2(b) for non-bifurcation cases. A careful review of these tables indicated ranges of (l/RM , l/Xsh ). when oscillatory and non-oscillatory instability occurs, the points (pointed by arrows) where the measure of instability of a particular kind (K4, K,, K, —K2K5/K6 , K, -K2K3K4) is most negative over the range of ( “Ru. , l/Xsh ), the kind of bifurcation (SN or Hopf) that occur, the values of (l/Rsh , UK“) where the bifurcation occurs (in general at points opposite to where the arrows point to) and whether the instabil— ity occurs in the regulated or unregulated model. The tables are so complete so that the dis- cussion that follows may be superfluous. 6.2.2.1 Effect of real-inductive local loading When the shunt local load was real-inductive, the parameter K, was negative for the entire R,,,-X,,, range. Furthermore, K4 was most negative at l/Rsh = 2.5 and decreased as 1/ X5,, was decreasing. This implies that the unregulated SMIB model was subject to undamped oscillations in the inertial dynamics, as summarized in Table 6.2(a). 184 6.2.2.2 Effect of real-capacitive local loading When the shunt local load is real-capacitive load, Table 6.2(b) shows that the param- eter K, is negative and decreasing for the entire R,,,-X,,I range from Figures E2.1(a-b) results. Furthermore, K4 is decreasing when the inductance l/Xsh is negative and ll/X,,,l increasing. This implies that the unregulated SMIB model is subject to undamped oscilla- tions in the inertial dynamics. Table 6.2(a) also indicates that the synchronizing term K, —- K 2K 5/K 6 became negative for light-real, heavy capacitive local loading, implying that the regulated model is vulnerable to SN bifurcation. The capacitive shunt is larger than x’d+.tT = 0.2 p.u. when Ixshl = 1 and becomes the same size as x'd-l-xT as [X34 becomes smaller, causing heavier capacitive loading. This appears to cause the network to induce the saddle node bifurcation. 6.2.2.3 Effect of the transmission network We notice that K,, is negative if no power transfer to the infinite bus occurs, regardless of the local real power load (R,,,), or the shunt capacitance/reactance (X,,,) at the generator side, which indicates that undamped oscillations in the unregulated SMIB model are asso- ciated with (i) heavily loaded load centers; (ii) local shunt capacitive system compensation; (iii) transmission lines with high capacitive/inductive charging; and (iv) when power trans- fer from or to the infinite bus is not possible. On the other hand, under the real-inductive and real-capacitive local loading operating con- ditions, Table 6.2(b) summarizes the non-bifurcations cases. Both K5 and the synchroniz- ing term K, - K2K3K4 are positive for the local load variations. This diagnostic study implies that in a local model, it is unlikely to experience (i) HOpf in the regulated model; or (ii) SN in the unregulated model. 185 Excitation Bifurcation Network Poster Transfer Home System Reg/Un Kind Test x. 1m... 11x... PM 0..., s u H K4<0 0.9 1‘- 10 [4.1—m 0 o EZ.l(a) u H K4<0 0.9 7,—“10 145—.1 0 0 152103) R SN K. -K2K,/K.<0 0.9 [_A_’___w .fi.’__,.6 0 0 E2.4(b) Table 6.2(a) Results on the local loading model bifurcation diagnosis Excitation Bifurcation Network Power Tiansfer Figure System Reg/Un Kind Test x. 1/R,.. 1/X... P... 0.... # R H K,> o 0.9 1 ,0 l——'10 o o 152.2(a) R H K5> 0 0.9 {___—'10 ’17—: 1 0 0 52.2(b) U SN K. -K2K3K.>0 0.9 5——-,,, {—10 o o 132.3(a) U SN K. -K,K,K,>0 0.9 0'__'“10 INT—'1, 0 0 152.303) R SN K, -K,1<,/1<,>o 0.9 (3—“10 ,,,,'——_, o o 132.4(a) Table 6.2(b) Non-bifurcation results on the local loading model diagnosis 6.2.3 Conclusion of the Power Transfer and Local Load Models From the diagnostic stress tests results of the power transfer and local loading models, we deduce the following three conclusions: (i) undamped oscillations in the regulated SMIB model are highly associated with the transfer of active and reactive power (supply or absorption) along transmission corridors: (ii) undamped oscillations in the unregulated SMIB model are associated with heavy capacitive local loading; and (iii) power transfer causes different kinds of bifurcations in the SMIB depending on whether the generator is 186 regulated or unregulated: the effect of relay-type control action. The regulated SMIB model is vulnerable to Hopf bifurcation when considerable reac- tive power transfer is absorbed by the generator —1 s Q” « 0. The unregulated model is vul- nerable to Hopf bifurcation in the generator inertial dynamics under local capacitive and local real loading operating conditions. Therefore, shunt capacitive load (generator absorp- tion) causes Hopf in the unregulated model, and transfer of reactive power to the generator (generator absorption) causes Hopf in the regulated model. The same conclusions hold in the regulated model for (Q,° > 0) and in the unregulated model for (l/Xsh > O). Hopf bifurcation occurs both in the regulated and unregulated models for heavy absorption of reactive power. This is not surprising due to armature reaction that reduces the rotor field flux for reactive absorption and is shown to cause instability due to insuffi- cient synchronizing torque between rotating rotor stator flux waves. The unregulated SMIB model is vulnerable to SN bifurcation when considerable real and reactive power is supplied by the generator 0.4 SPgs 1, 0 « gas 1. The regulated model is vulnerable to SN bifurcation in the generator inertial dynamics under heavy-local capacitive, and heavy-real local loading. Therefore, transfer of reactive power to the infi- nite bus causes SN in the regulated model, and shunt capacitive load (generator absorption) causes SN in the unregulated model. Finally, no saddle node bifurcation occurs in the regulated power transfer model, and no saddle node in the unregulated local loading model. In correlation, no Hopf bifurcation occurs in the unregulated power transfer model, and no Hopf in the regulated local loading model (l/R,,,, l/Xsh). 187 6.2.4 Combination Local Load-Power Transfer Model The diagnostic study now focuses on the more realistic model consisting of a local load in combination with power transfer requirement patterns. The variation of the machine damping and synchronizing torques under the effect of the various local loading patterns, and power transfer patterns is investigated for both the case where the generator is regulated (K5,K, — K 2K S/K,,) and the case where the generator is unregulated (K4,K, -K,K3K4). Various transfer and loading patterns are considered: (i) the effect of power transfer for a local loading pattern, (ii) the effect of local loading for a power transfer pattern, and (iii) the effect of the transmission network. Simulation for the combination model bifurcation diagnosis are shown in the figures in Appendix E3, and the bifurcation results are summarized in Tables 6.3(a-c). 6.2.4.1 Effect of Power Transfer First, the local shunt load and the transmission network inductance are fixed respectively at R5,, = 0.65 p.u. , X,,, = $0.17 p.u. and x,2 = 0.4 p.u. , while the power transfer to the infinite bus (Pm , Q”) is varied. Note that x, = 0.4 p.u. rather than 0.9 p.u. as in the local load center model. The bifurcation results for the effect of P” and Q0 on the parameters K4 and K5 as well as K, -K,K5/K6 and K, --K2K3K4 are summarized in Table 6.3(a). The occurrence of bifurcation as shown in the table indicates that (i) when the load Xsh is inductive, K5<0 for both the case when Q” 2 0 and the case when Q” s 0; (ii) when X,,, is capacitive, K5<0 for both when Q” 2 0 and Q“ s 0; and (iii) the parameter K, >0 for the power transfer stress tests except when the local load is capacitive and reactive power is being supplied by the infinite bus (X3). = —0.17 ,Qm s 0 ). Note that while (i) and (ii) confirm the power transfer model diagnosis results, (iii) occurred only in the combination model, 188 and thus extends the diagnostic conclusions. This result (iii) may be due to the effect of the large shunt capacitive load which caused Hopf bifurcation in the unregulated model for shunt capacitive load and no power transfer. Power transfer stress on the system seem to also ultimately cause (i) SN bifurcation in the regulated model when the local load is real—capacitive X3, = -0.17 and P” increasing, Q0, decreasing; and (ii) SN in the unregulated model when the local load is inductive X3, = 0.17 and Po, and 0,20 are increasing. These results in (11) above extend the SN bifurcation results obtained from the power transfer model. 6.2.4.2 Effect of Local Loading The power transfer to the infinite bus, the local shunt load and the transmission network inductance are fixed respectively at P,° = 0.8 p.u. , Qco = 0.5 p.u. and xe = 0.4 p.u. , while the local shunt load at Rs]: and X,,, are varied. The bifurcation results for the effect of R3,, and X“, on the parameters K, and K, as well as K, -K2K3K4 and K, — Ksz/K,5 , are summarized in Tables 6.3(b-c). The occurrence of bifurcation as shown in the table indicate that when the model is regulated, (i) the parameter K, changed sign for the light- real, heavy-inductive local loading conditions (approximately l/Rshs 1.67, 1/X3h22.0), and is decreasing for lighter l/Rsh loading conditions as can be seen from Figure E3.9(a); and (ii) the parameter K, is positive with no change in sign when the local load is real-capacitive from Figure 3.9(b) results. Given the system is unregulated, (i) when the shunt load is real-capacitive, the parameter K, changes sign under the heavy-capacitive loading conditions (l/Xsh s -1.67 , generator reactive absorption); and (ii) when the shunt load is real-inductive, the parameter K, is positive which implies that shunt real-inductive real loading patterns are unlikely to 189 cause undamped oscillatory behavior in the unregulated machine dynamics in the SMIB. Figures E3.ll and E3.12 in Appendix E3, show the variation of. the machine synchronizing torques under the effect of the machines local loading patterns, for both when the system is regulated (K ,-K2K5/K6) and unregulated (K,—K2K3K4). For the regulated system, the summarized figures in Table 6.3(b) indicate that (i) the synchronizing torque coefficient is positive for real-inductive loading conditions; and (ii) negative for light-real, heavy-capacitive loading (l/Rsh small, |1/X,h| 2 1.67 ). When the generator exciter is disabled, the synchronizing torque is positive for the entire R,,,-X,ll range. This , implies that the system is unlikely to experience static kind of instability when the generator is unregulated. 6.2.4.3 Effect of transmission network The effect of the system’s transmission network inductance x, on the system synchronizing and damping forces is simulated in Figure 6.4. The system operating conditions were at P, = 0.8 p.u., Q, = 0.5 p.u., X3, = 0.4 p.u., and two real loading levels were tested, R5,, = 0.17 and R3,, = 0.87. For this case, the damping torques are plotted versus x, We see that (i) both K, and K(, both decrease as x° increases, and (ii) K0 is very much more negative for heavy loading conditions, but K5 is slightly more negative for light loading conditions. This behavior confirms the results obtained from the power transfer model. 190 Excitation Bifurcation Network Power Transfer Figure System Reg/Uh Kind Test x, R... x... P... G... # R H K5<0 0.4 0.65 +0.17 6?)? 65—1 53101) R H K5<0 0.4 0.65 +0.17 (ST—:51 1’ M 53103) R K5<0 0.4 0.65 0.17 (3?)? 0 ‘1 53.20.) R H K,0 0.4 0.65 +0.17 0, , 6—‘ E3.3(a) U H K,>() 0.4 0.65 +0.17 5.1—1 :5— E3 .3(b) U H K,>0 0.4 0.65 .017 (3.1—1 6?"— E3.4(a) U SN K. -K,K,1<,>o 0.4 0.65 +0. 17 5,—1 .173“ 53.7(13) 0 SN K. -1<,K,x,>o 0.4 0.65 -017 57—, .T_"‘- 53803) R SN K. -K.K,/K,>0 0.4 0.65 +0.17 (ST—1 F E3.5(a) R SN K. -K2K,/K.,>0 0.4 0.65 +0. 17 67““, (fl—“*7 E3.5(b) R H K5>0 0.4 {—10 127—11 0.8 0.5 E3.9(b) R H K,>0 0.4 1 ,0 f ,0 0.8 05 5310(3) R SN K.-1<,1<,/1<,>0 0.4 1 m i 10 0.8 0.5 53.1203) U SN K.-K2K,K,>0 0.4 1 '10 1 ,0 0.8 0.5 E3.1l(a) 0 SN K.-K2K,K,>0 0.4 r 810 2,, T. 0.8 0.5 53.1103) Table 6.3 (c) Combination model: Non-bifurcation cases 0.1 1 _04l- ............................................. _. ---— K5. lash-0.17 _0-5-H,"n"Kd,'.RSh*o.l1?HU-'j .......................... .1 '-"K5._RSh-0.87 . _OIG- -.—.——-.Kd..RshsO.87.. . -f ..................... .... 1' i i 0 0.2 0.4 0.6 0.8 X903“) Figure 6.4 Combination Model: Damping torque coefficient Kd and parameter K5 for different levels of real local loading 192 6.2.5 Discussion of the Combination Model Results How do the power transfer and local loading models predict, disagree, or possibly extend the results obtained from the combination model?. The combination model diagnostic study shows that combining both models may extend or confirm the results from each individual model, depending on Operating conditions or stress is dominating. The following conclusions are drawn: a) b) Although local capacitive or inductive local loading (XS.,<0 or Xsh>0,P°° = Q” = 0) did not cause Hopf bifurcation in the regulated model in the local load diagnostic Study for reactive supply (Q2 0) to the infinite bus, local inductive or capacitive loading in combination with reactive power supply to the infinite bus (X31. = $0.17 , Q0, 2 0) appears to drive the system to oscillatory instability (Figures E3.1(a) and E3.2(a)). This extends the local loading model results and allows the diagnostic conclusion that the regulated model Hopf bifurcation can occur when the reactive power is supplied to the infinite bus. The problem is aggravated when real power supply P” increases and real loading (UR...) decreases. This result may be due to the transfer of power to the infinite bus over a weak but extremely weak line. Although local capacitive loading (X,.,<0 or X,.,>0,P,, = Q... = 0) did not cause Hopf in the regulated model, for reactive power supply from the infinite bus (Q, s 0) local capacitive or inductive loading (Figure E3.1(b), E3.2(b)) in combination with reactive power supply from the infinite bus (X3, = $0.17 , Q“, s 0) seem to drive the system to oscillatory instability. This extends the local loading model results and allows the diagnostic conclusion that the regulated model Hopf can occur when reactive power is supplied from the infinite bus. The problem is aggravated when real C) (d) 193 power supply P” increases and real loading (l/R,.,) decreases. This result may be due to the large active power transfer to the infinite bus, and reactive power transfer from the infinite bus. For saddle node bifurcation, although power transfer stress seems to affect the unregulated model synchronizing torque, and not the regulated model, in the combination model the following bifurcation occurred: 0) (ii) SN bifurcation in the regulated model occurs when the local load is real- capacitive(xs,, = -0.17 ), active power transfer to the infinite bus is increased and reactive power from the infinite bus is increased (Figures E3.6(a-b)). This bifurcation was not observed in the power transfer model diagnosis, but appears to be caused by the presence of the heavy capacitive local load ll/Xshl == 6 , as diagnosed by the local loading model. SN in the unregulated model for power transfer (P002 0,Q” 2 0) to the infinite bus and either X,..>0, or X,.,<0). These results (Figures E3.7(a), E3.8(a)) are consistent with the SN bifurcation results obtained from the power transfer model. The effect of local loading stress seems to be less severe in the presence of real and reactive power transfer to the infinite bus. The local loading model diagnosis Showed that K 4 s 0 for the entire range of R,.,-X,.,. In the combination model, Figure 133.10 (a) shows that the parameter K4>0 for a Significant portion of the operating range- of Rsh-Xs... 194 6.3 Diagnostic Summary The diagnostic study demonstrated that the cause and effeCt of a specific kind of bifurcation in the SMIB model are dependent on (1) the effect of the relay-type control action, i.e whether the system is regulated or unregulated; (2) local loading-power transfer operating requirement patterns; and (3) transmission network inductance. The comprehensive power transfer and local loading pattern variations in the three study models, has been the mean for drawing (if possible) conclusions for diagnosing the causes of SN and Hopf bifurcations in the SMIB. Furthermore, the results establish the classes of saddle node and Hopf bifurcations that were observed in the SMIB model. 6.3.1 Classes of Hopf Bifurcation in a SMIB Model Ill; Power Swing Mode (Regulated) This oscillation mode occurs when K5 is negative in the single machine infinite bus model, when the voltage regulator is still operating within its fields current limits. This oscillation is observed in the generator angle/speed/fiux decay dynamics, and is associated with low frequency inter-area oscillations. It is a supercritical (stable) Hopf bifurcation since these oscillations are stable inter-area oscillations. It is often referred to as a power swing mode, Since it is brought on by considerable power transfer between two areas. The power swing mode can be associated with a single generator, two generators or with groups of generators. There are two subclasses: l. Reactive Power Supply Swing Mode (regulated) 2. Reactive Power Absorption Swing Model (regulated) 195 fig; Voltage Collapse Mode (Unregulated) Occurs in the inertial and fins decay dynamics of the generator in a Single machine infinite bus model, when the exciter is completely disabled, and the single machine parameter K, is negative. These classes of I-Iopf bifurcations are likely to be the limit cycles that develop as a precursor to voltage collapse and thus are called voltage collapse modes. It is not clear whether voltage collapse oscillatory modes are supercritical or subcritical Hopf bifurcations. The diagnostic study showed that these undamped oscillations in the unregulated SMIB model are associated-with (i) heavy loaded local load centers; (ii) local shunt capacitive system compensation; (iii) transmission lines with high capacitive line charging; and (iv) when power transfer from or to the infinite bus is not possible. There are two subclasses: l- Reactive Power Supply Voltage Collapse Mode (unregulated) 2- Reactive absorption Voltage Collapse Mode (unregulated) It should be noted that there is a Voltage Mode that is subcritical Hopf bifurcation that depends solely on the exciter and flux decay dynamics. 6.3.2 Classes of Static Bifurcation in a SMIB Model (SN) Similarly, two classes of Static bifurcation may occur on a Type 2 power system model: SNJ Static Bifurcation in generator inertial dynamics (Regulated) Static bifurcation in generator inertial dynamics may be associated with the loss of transient Stability [36] such as the voltage instability that occurred in the Czechoslovakian system after tripping lines in a major interface [30]. This class of static bifurcation is likely to occur when the disturbances are very large and the l 96 transmission interfaces or boundaries are weak, so that loss of stability occurs soon after the contingency before the field current level/duration limits are reached. It would appear from the results in this chapter that the Static bifurcation in generator inertial dynamics occurs the local network is capacitive near a local generator, real power is Shipped away from that generator, and reactive power is shipped to that generator. 5N1 Static Bifurcation in generator inertial dynamics(Unregulated) This bifurcation has been called classic voltage instability and is due to insufficient reactive supply in the EHV system, an example of which is the blackout that occurred in the French system [37]. This class of Stability has been analyzed with the use of a sensitivity matrix T between reactive generation and voltage at generator internal buses [51] that is diagonally dominant for inductive networks (networks that absorb reactive power), and not diagonally dominant for capacitive network (network with long high voltage transmission lines). Classic voltage instability has been described as occurring when generators reach their field current limits and the field current limit controllers disable the exciters. An important implication to this result is that study- in g this saddle node bifurcation must include the field current limit controllers action in the model. VII Thesis Summary and Future Research 7.1 Dissertation summary This thesis addressed three major difficulties in diagnosing the cause of a specific Stability problem in a power system. These difficulties are (a) the size and complexity of the model, (b) the various kinds and classes of bifurcation that can occur in the model, and (c) the lack of smoothness in the dynamic behavior after a disturbance or during operation, (d) the inability to determine a subsystem of the full model that experiences the same bifurcation as the full system for each kind and class of bifurcation. These different problems in a power system have been and Still are in need of better methods for Stability analysis, developing application software, and problem diagnosis in a large power system. This thesis will address all four problems and hopefully provide the basis for methods that can be utilized to establish where, why, and what can be done to prevent loss of Stability for three kinds (kinds) of bifurcation: Hopf, saddle node, and Load Flow bifurcations observed on a Specific power system model as well as all of the specific subsystems of the model affected by each particular kind of bifurcation (classes). The contributions of this work can be summarized in the following four sections: 197 198 A) Diagnostic Classification First, the diagnostic study started by classifying (a) the model types used in power system studies, (b) the kind and (c) class of bifurcations. This classification is performed in Chapter 2 of the thesis. The kind and class of a bifurcation observed depends on the model used. Four model types are classified: Typel model is the simplest model where both generator and loads are modeled by algebraic equations and under load tap changers (ULTCS) and switchable shunt capacitors controls are modeled. Type 4 model represents both generators and loads by differential equations and yet all the ULTCS controls dynamics and switching capacitor inductors, etc. are included. Type2 and Type 3 models are simplified Type 4 models that eliminate the load dynamics and generator dynamics respectively. The complete Type 4 dynamic model has been developed only for small test systems. B) Epoch Bases Trajectory Stability Assessment Method Virtually every local bifurcation is brought on by a sequence of equipment outages and hard-limit discontinuities enabling or disabling various controls. Although dynamic system theory is fairly well developed for smooth systems, it is currently being developed for nonsmooth systems. These hard limit discontinuities do not cause immediate enablement or disablement of a particular control. The dynamical system theory for nonsmooth systems [27] is not sufficiently well developed to be able to determine the sequence of hard limit discontinuities when each discontinuity has a different length time delay. Current dynamical system theory for smooth system can be applied to intervals (1315”) called epochs where discontinuities occur only at [1,] , the discontinuous change 199 at [1,] can be computed; and the time T, = 1,:+1 4; is long compared to the settling time of the dynamics. Therefore, a method aimed at assessing asymptotic stability of the transient during 0,7,1?) and the Stability of quasi equilibria or limit cycles in each epoch is proposed. This method called the Epoch Based Trajectory Stability Assessment, is based on a trajectory simulation method implemented by [28] that does not require time simulation in (131;, ,) , but only computes the quasi equilibria (22,53) for each epoch for 1' = 1.... N. This method allows the assessment of the stability of individual trajectories for a power system by assessing stability for every epoch through (a) assessing stability of the sequence of quasi equilibria or the limit cycles, and (b) the asymptotic stability to each quasi equilibrium or limit cycle. C) Bifurcation Subsystem Method: A third a major contribution in this diagnostic framework is the establishment of the Bifurcation Subsystem Method, This method attempts to identify the smallest subsystem, a subset of the equations of the full system model that experiences and causes the same bifurcation that occurs in the full system model. A bifurcation subsystem can loosely be defined as a truncated portion of the actual power system model that experiences the same bifurcation in the full system model. The Bifurcation Subsystem Method utilizes the geometry associated with the various submatrices of the full system differential-algebraic Jacobian J = [f‘ f5] since x 8 y Det(J) = Det(gy) - Det(fx—fyg;lgx) when g), is nonsingular; De1(1) = Det(fx)-Det(gy—gxf;lfy) when f, is nonsingular; and with the eigenvectors associated with the bifurcating eigenvalue to establish 2110 conditions for existence of such systems. The fact that the Jacobian of the truncated portion of the system (gy) must be Singular somewhat restricts the set of possible subsystems that one is using to find the bifurcation subsystem. The bifurcation subsystem must satisfy both the linear and the transversality test conditions for bifurcation that is used to test for bifurcation in the full system model. The null space of the effects of the truncated system (fyg;lgx) must be common with the null space of the bifurcation subsystem (15,) in which the right or left eigenvector of the bifurcating eigenvalue resides. The portion of the full system external to this bifurcation subsystem has no effect on the linear conditions for bifurcation to be satisfied in the full system model only at bifurcation (asymptotic), or as bifurcation develops for 11¢ 11" (structural), and yet the equilibrium point is dependent on the variations in the full system equations which results in bifurcation in both the full system and the bifurcation subsystem. Propositions for the existence of bifurcation subsystems for a particular kind of bifurcation were derived and proved. Possible bifurcation subsystems in a power system Type 2 model are: (a) Differential Bifurcation Subsystem Experiencing SN; (b) Differential Bifurcation Subsystem Experiencing Hopf; (c) Algebraic Bifurcation Subsystem for LF; (0) Differential Bifurcation Sub-subsystem Experiencing SN; (e) Differential Bifurcation Sub-subsystem Experiencing Hopf; (g) Differential-algebraic Bifurcation Subsystem. It should be noted that the bifurcation subsystem is not the center manifold of a specific bifurcation since the same bifurcation subsystem could exist for saddle node bifurcation in the flux decay dynamics, and for saddle node bifurcation in the inertial dynamics. The 201 sensitivity matrices Ks, KD and T also do not necessarily identify the center manifold dynamics, and thus do not necessarily identify the exact class and location of a specific bifurcation. The research has shown that what is desired is not to determine the center manifold dynamics, but rather a subsystem of the original model that experiences the same bifurcation as the full system. The bifurcation subsystem may in most cases be an excellent approximation of the possibly larger subsystem that causes and produces the bifurcation in the full model, as it is the desired result. The root locus of the bifurcating eigenvalue and the participation factor information provides additional diagnostic information on such a subsystem. Determining sensitivity matrices Ks, KD or T as a function of the model parameters K.-K., and then studying how K3, KB or T vary with changes in operating conditions and hard limit discontinuities also can assist in indicating the subsystem that produces, causes and cures a specific bifurcation, in special cases. D) Stability Conditions and Diagnosis: Lastly, the Epoch Based Trajectory Stability Assessment Method and the Bifurcation Subsystem Method are the framework for initiating a diagnostic study for Hopf and saddle node bifurcations on a single-machine-infinite bus (SMIB) model to (a) computationally demonstrate and validate the use of the bifurcation subsystem method to identify bifurcation subsystems in the SMIB; (b) identify the classes of saddle node and Hopf bifurcations that may occur in the SMIB model: the generator angle/Speed dynamics and angle/speed/flux decay dynamics are bifurcation subsystems for Hopf and SN bifurcations in both the regulated and unregulated SMIB models, respectively; (0) gain better 202 understanding on the effects of hard limits by assessing the effect of excitation control on the types and classes of bifurcations produced; and (d) establish an understanding to the causes, operational and structural stresses, and stabilizing requirements for SN, local and inter-area modes of oscillations that occur in generator inertial, flux decay and exciter dynamics of the SMIB model. This study was performed by formulating three study model: (1) a power transfer model; (2) a local loading model; and (3) a combination power-transfer/local loading model. The diagnostic study demonstrated that the cause and effect of a specific kind of bifurcation in the SMIB model are dependent on (1) the effect of the relay-type control action, i.e whether the system is regulated or unregulated; (2) Local loading or power transfer operating patterns; and (3) transmission network inductance across interfaces. 7.2 Future Research The following topics are subjects for further investigation: 1. Develop a load flow program based on the Bifurcation Subsystem Method to obtain load flow bifurcation subsystems that experience the same bifurcation as the full sys- tem, but also include what produces and causes the bifurcation. Evaluate how load flow subsystems change with (a) decoupling of P-V and Q-O Jacobian submatrices; (b) generator reaching field current limits; (0) under-load tap changers reaching tap limits; ((1) switchable shunt capacitors reaching shunt susceptance limits; and (e) the effects of multiple contingencies. This knowledge allows the prediction of the most critical bifurcation subsystem in the system for a load flow model, given any equip- ment outage or operating condition combination using proximity measures and diag- 203 nostic tools. A heuristic method has been developed for identifying bifurcation subsystems that experience the bifurcation of the full system, given certain assump- tions are true. The work would establish the bifurcation subsystems when the assumptions are not made and when they are true. Develop a dynamic stability program based on the Bifurcation Subsystem Method to obtain dynamic bifurcation subsystems that not only experience the bifurcation of the full model, but also produce and cause the bifurcation. Evaluate these subsystems change with (a) decoupling the generator/load dynamics of the system; (b) hard limit discontinuities in generator and network controls; and with (c) the effects of multiple contingencies. This knowledge allows the prediction of the most critical bifurcation subsystem in the system for a differential-algebraic model for any equipment outage or operating condition combination using proximity measures and diagnostic tools. Develop proximity measures for each class of bifurcation subsystem found in a load flow or a differential algebraic model. These proximity measures can be (a) reactive reserves on a particular group of generators; (b) measures of properties of sensitivity matrices like K5, KB or T; (c) measures of properties of the system matrices K.-K.,; or ((1) measures of properties of the J acobian of the bifurcation subsystem that produced the bifurcation. Develop diagnostic procedures for identifying the kind, class, location, cause, and remedial action needed for each bifurcation subsystem or class of bifurcation identi- 204 fied in a load flow or differential algebraic model. The diagnosis Should be based on identifying how the system matrix A; sensitivity matrices like K5, KB or T; sensitivity matrices [(1-K6; or Jacobian matrix J are affected by load, transfer, reactive reserve exhaustion, etc. that produce each bifurcation class. APPENDICES APPENDIX A POWER SYSTEM DYNAMIC MODEL Appendix A Power System Dynamic Model A.1 Nonlinear Model The dynamic model of the full power system [52] is described in this section. Without loss of generality, we assume that the connection between each generator and its high-side bus is equivalently one-to-one, and that the power system has 11 buses total buses including m machines with m high—side buses and (n-m) other (load) buses. The subscripts used in the various variables used in the model equations, are described here: Subscript t: represents the generator terminal PV—buses (including switchable shunt capacitor buses. Eg: P0,, P0, V. etc, t = 2, m. Subscript H: represents the generator-transformer high-side PQ-buses. Eg: P0,, PC, Subscript L: represents the rest of the P-Q buses. Eg: PC, Subscript G: represents the active or reactive power generations Subscript C: represents the non-voltage (i.e constant) active or reactive power load or injection A.1.1 Swing Equation For each synchronous machine, 8 = (1)—0),, M0) = —D(0)—0)(,)+Pm—Pe Where 5 : angular rotor position, measures in electrical radians of the relative to a syn- chronous rotating frame 0) : angular speed 205 (no : nominal angular speed D : approximate damping constant due to damper winding effects( pu/(rad/s) ) M = 214/(1)0 : moment of inertia (sec/(rad/sec)) l—l : moment of inertia constant (MW - sec/MVA) Pm : mechanical input power to the machine P6 : electric power output of each generator A.1.2 Flux Decay Equatlon For each synchronous machine, K53 = l-(xd-x'd)-Bq > 1 U4, = Gasin(5—9,)-chos(5—9,) R -x Ga = 7—“—— ; B = _T—L— Ra-l-x’dxq Ra+x’dxq Where T’do : d-axis open-circuit transient time constant Efd : field winding voltage Eq’ : induced stator voltage due to the flux linkage of the filed winding A.1.3 Excitation system A typical excitation system in shown in Figure Al, is responsible for controlling the terminal voltage. Using the following notation: VC: output of load ( or line drop) compensator V3: output of power system stabilizer (PSS) Vref: reference (set point) voltage VD: voltage detector output KD: dc (direct current) gain of voltage detector Efdi l/KSE: 207 time constant of voltage detector transient gain reduction lagging time constant of TGR leading time constant of TGR output of automatic voltage regulator (AVR) dc gain of AVR time constant of AVR output of exciter (field voltage of synchronous machine) dc gain of exciter. Note that KSE is function of the saturation effect of the exciter. TE/KSE:time constant of exciter KF: TF: washout gain of stabilizing feedback of excitation system washout time constant of stabilizing feedback of excitation system The state space model of the [BEE-DC] excitation system is given by .%DV; roF rAv, rev, P -KA L O :2 TB 7'. _(1__E TB H A.1.4 Power System Stabilizer A power system stabilizer improves the dynamic stability of a synchronous machine and/ or a power system, using other regulator input signals in addition to terminal voltage through its excitation system. Some common stabilizer input signals are accelerating power, machine speed, frequency and terminal voltage. A state space model for the PSS model shown in Figure A2 may be written using V50, V52 and V54 shown in the block diagram. '- '1 KDVC 0 TC TC (1 - fi)(an + VS) L 0 d 208 P ..1 o 0) ‘ —KS TSVSU l—Z‘fl —1 0 v30. (1 TS')K Tszvsz = Tsz Vsz + 752 5 VS! T. V54 T33 TSl T53 _1 V54 T33 T51 54 “FT—1’?— 1__ —KS L S4 52 54 L T54 T57 T53751 753751 753 s = EEKS 51+‘fsj1'T—SZ so 754 A.1.5 Speed-Governing-Turblne System Model The speed—goveming system controls the mechanical input power, using a turbine system and rotor speed as an input signal. Figure A3 present the general model of speed governing system, for steam-turbine and hydro-turbine. In the figure, PGV is the power at the gate or valve and represents the input of the turbine system. The state space model of the speed-goveming-turbine system may be obtained by combining both models. a. Speed Governing-Turbine .- ' W P — T1 q TP~ —1 o 0 0 0 0 P ' T1 3 f” -1—1 0 0 o 0 6V T, TCHPVHP = 0 1 -l 0 O 0 PVHP + _KG(1'F‘)(°”(°0)+P0 3 3 3" ° 3 ° 1 -1 T P PIP 0 RH?” _o o o o 1-1, I, 0 _TCOPLP_ - LP 0 PM = FVHPPVHP+FHPPHP+FIPPIP+FLPPLP 209 b. Hydraulic Governing-Turbine T . K (iv—Warm) TIP! -1 0 0 P: 0 T1 " TPGV= “1‘1 0 Pav“ 3 _ 0 a ~-a —KG(l——3)(w-wo)+l’o TWPW W0 WW PVHP l _ 0 _ . “13021 . “13021 PM=FGVPGV+PW 9 “WC: 0211 ’ “WV/'51— , FGv-aza‘ a” A.1.7 Load Compensator Model The load compensator uses the terminal voltage and current as feedback measurement of the excitation system model. A state space representation of this model is given by: ‘ K4254 + Kd7V: < C a. l Where, K42 = -(XcGa+Rch) K42 = Rcaa-XCBq K47 = (l +XCB’d-RCGa)sin(8—9,)—Kd2cos(6-9,) Kq-, = (1 -Kq2)cos(8-9,) -(XCGa +RCB’d)sin(8-9,) A.1.6 Algebraic Load Flow equation (1') At Internal Bus: for each synchronous machine e ’2 I 2 P = E 40515-5 thTPE+VtGEt ,2 , 2 Qc E 4355+‘5quUQE—VzBEx A.2 210 (ii) At terminal Bus: for t = n+1. n+m; h = 1.,m 2 , I 2 ,, ‘PDt Vt (th ‘ VtE thq + V! “’Cth + Gm) ' Viv/1TH: 2 , 2 "Q0: _Vtth+-VtE qth—Vt(BCth+Bth)_ vrvhUth (iii) At high-side bus:for h = l, m; t = n+h, we have 2 w 2 n w I! J: I: _QDh 1:] 1: (iv) At load bus: for 1': m+1, ..., n, we have 2 n ‘ n ‘PDi = Vi .21((’Cij+ Gij)‘ .ZlVngng /= 1= ’3 ' n 2 ‘Qoi = _Vi .ZI(BCij+Bij)" _leiVjUij j: I: “‘1 ll i]- Gijcos(9,- — 91-) + Bijsin(9i- 9}.) Q ll Q ll Q ll Linearized Power System Model 2 2 n n "Vh(ch+Bm) —VthUht'Vh 2- (HChj+th)' 2 VthUtj . . 1 1:} i¢j i¢j i¢j i¢j i¢j Taylor series expansion is used to obtain a linear model about a nominal isolated equilibrium point. Linearization is necessary for the load demand model, power balance equations, swing equation, flux decay equation and load compensator. The jacobian matrix full linear model is expressed in terms of submatrices that are explicitly given below. F . 7 — '1 TXXAXX AXX AXE Axr; 0 Are Axv AXX TEEAXE AEX AEE 0 A55 A59 AEV AXE = +BOAU0 TSSAXS ASX 0 0 ASX 0 0 AXS —APC APX 0 0 0 Are APV A9 L-AQc_ LAQX 0 0 0 AQa AQV‘LAV‘ Where, r n-1c w '- —l- Txx o 0 0 AXXAXEAXG 0 Txx 0 0 oi AMA”) f.= 0 TEE 0 0 AEXAEE 0 A55 :0: o TEE 0 o AseAev ' 0 O TGG 0 AGX 0 AUG 0 O 0 TGG O 0 0 _0 o 0 rs, LASX o o ASX‘ L0 o o TSS‘ Lo 0_ t I AXX = [Amt A5‘ AE'F] : states of mechanical and flux decay dynarrucs; L. r I AEF = LAB"; AE’d‘ Atyld' Atqu’] :states of flux decay dynamics F t - t t t t t . ‘ ' AXE _ LAVD AVF AVA AVB ABM] . states Of CXCltathl'l system I AX“ [AFC “’th APVHPt APHP‘ APIP‘ APLP’] :states of speed-governing-turbine systems; I = I t l : . . ; AXS [Avso Avsz (”54] states of power system stabilizer t A9 = [AGT' A9”! 491.1] : angle variables at network buses; l AV = [AVT' AVH' AVL‘] :voltage variables at network buses; I Me = [APCT‘ APCH‘ APCL‘] : coefficients of non-voltage-dependent active power load demand model; 1 I AQC [AQCT AQCHI AQCLI] : coefficients of non-voltage-dependent reactive power load demand model; diug[M I T’do] : diagonal matrix composed of inertia constants of synchronous machines, identity matrix, and time constants of flux decay dynamics; TEE = diag [TD Tr T A TB TE] : time constants of the excitation systems dynamics; Tar; “mg [Tr Ts TCH TRHl TRHi Tea] : time constants of speed-governing-turbine systems; T55 = diug[T5 T52 T54] : time constants of power system stabilizers; M = diag[Ml M2, Mm] :inertia constants for m synchronous machines 1 .° m x m identity matrix Tdo : m x m diagonal matrix of time constants of flux decay dynamics —D _Kj’.‘l «if, 0 0 0 0 0 0 0 Fvnp Fur Fm FLP =1 0 o Ax15=00000 Axa=oo 0 o o 0 0454-19,, 00001 00 o 0 o 0 “KS 00 -l 0 0 [($100 ‘ 1211‘s 00 Ass-'- ’21 ‘I O Axe: o ()0 O—Kfi'l—Kg’z 04(31ng Am BIH 0 0 0 0 AQX" 0' o o APB: Azu 3211:1321”. 0 0 0 O 0 0 0 32w BZLL Cur 01H 0 Ant 33H 0 CM 03H 0 C2}! 02m: 02m. A06 = A411 34m: 34m. AQV = Cm D4HH 04m. 0 02LH Out. 0 34w B4LL 0 04”! 041.1. EX GX Axv _ q -I 0 0 0 0 0 0 0 —1 -l 0 _ 0 0 0 0 0 TFKF -1 KATBCT43T21 KATBCT43 KATBC AEE = —KATBC -KATBC 4 KA —KATBCTF KF I .T T l .T l . 3c 43 21 BL 43 3t -1 "BC "BC 0 " “’BCTF KP L 0 0 0 . L 0 0 I 0 4‘55 . P 0 K K- K K.“ ~ 1 - , D c1 DLLZ —KDKC100 -KDKC7OO 0 0 0 0 00 0 00 AAw 0 0 Axe“ 0 00 AEV‘ 0 00 A 0 0 0 00 0 00 Bw L L0 0 0 . L 0 00 0 00 I T—IT K 00 ' r - - r - [‘12]G 400000 0 0 00 ..1 -I-I 0 0 O 0 0 o 01 ‘T1T2KG 00 A = 01—1000 3 = KATBCO B = 00 0 00 00 001—100 50 IT 0 GO 00 0 00 0 0 0 I .1 0 ’ BC 00 0 00 _0 0 0 0 ’-. L 0 0_ _00_ . 0 00. qt -KP100 0 00 K5700 tq . _ tq . __ tq . _ '4 AlK=KPl+Al . AIK-KP7+C1 . A3K‘KQ1+A3 . C3K-KQ7+C3 . . ‘1 . ‘1 . ‘1 ’43 =1+T43 .12,= ”721 ~ T43 = T S4Tss . T21= T SZTSI ’ TBC = TB Tc ’Bc = "THC ; AM = KATBCT43T21KS ? BBw = [BCT43T21KS K54 = UPE _(xd _ x’d) ' thvl = klUdt-Gauq-x’dflq, = kZUdt , 2 , E qlePE + V! (xq -x d)(Tdth(- Udtht) , 2 , -E qv,TQE+ 2v,(xq—x (1)0414, -(xd -x’d) - U4, led, — Ga(xq —x’d)Uq, kZTdt 4‘ I Km 2‘5 qGEE‘ VITPE 45’qu — mugs —2V,BE,—E’qUQE KD/( l+s$ ) KD/( l+s'§ ) PSS VAmax + T ' Ef + 1+5 L KA/(1+STA) -——~q 1/(Ksr~:+S"E) d > l+sTB TOR V_/ Exciter Amin SKF/(l-l'STF) (a) PSS '1‘ i‘zc VAmax + V5 + VA Efd - KA/(l-l-STA) —q l/(KSgI-S'E) HSTC + Exciter 1+sT3 VB VAnu'n 5F Tr VF sKFJ(1+sTF]' ‘ (b) Figure Al. [EBB-DC] Excitation System Model V 215 Vsr ——u SKSTS l + STSI 1 + STS3 VS 1 + 31‘s 1 + 5T3; l + sTS4 ' T51 1'51 . T52 T54 + l'TSi/r ,2 . 1°Tsfl54 l + sTsz V82 1 + 8T4 v84 Figure A2. Power system stabilizer model (00 + 1)0 PU a PllSMAX P 0) +‘ KG(1 + s13) 1__ 1 av ’ l + 8T1 - T3 3 pDOWN PMIN + P PGV Fvap + A 1 Pvnp l + STCH m l _ . 1 1+ Slum Pup 1+ srm , Prp 1+ s'rCO Pu. fie?) ........................................................................................... Hydrp mo P° PMAX m_.$‘__ KG(1 + 5T2) (g r ’ Pov (1+sT1) (1 + 5T3) J P a23 '3 (313321 ' “1323)T pM GV 1 + sauTw Figure A3. Speed-Governing-Turbine System Model APPENDIX B COMMON DEFINITIONS FOR POWER SYSTEM STABILITY APPENDIX B COMMON DEFINITIONS FOR POWER SYSTEM STABILITY In order to facilitate analysis of stability problems, identify the essential associated factors and to formulate methods of assessing and improving stable operation, stability problems have been classified into appropriate categories. A classification recently presented in [29], summarizes the most widely used stability definitions, the associated concepts, and the related terms. The classification is based on the following considerations: (a) the physical nature of the resulting instability, (b) The size of the disturbance considered, (c) the devices, processes and time span required to determine stability. 8.1 Rotor angle stability Any unbalance between the generation and load results in a net accelerating (or decelerating) torque exerted on the rotors. From the swing equation, we see that this transient causes the rotors of the synchronous machines to “swing”. If the net torque is sufficiently large to cause the rotor to swing far enough so that one or more machines “slip a pole”, synchronism is lost. However an appreciable increase in the machine’s rotor angular velocities does not necessarily mean that synchronism will be lostThe important factor is the angle difference between machines, where the rotor angle is measured with respect to a synchronously rotating reference. Rotor angle stability is the ability of interconnected synchronous machines of a power system to remain in synchronism after 217 218 being subject to a transient disturbance. 3.2 Small signal ( or small-disturbance) stability It is the ability of the power system to maintain synchronism when subject to small disturbances. The disturbances are considered small enough to cause a small change from the operating equilibrium point, so that linearization of system equations is permissible. The powerful tool of linear system theory may therefore be explored for purpose of analysis. 8.3 Transient stability Transient stability is the ability of a system to maintain synchronism when subjected to a severe disturbance. The system response involves large excursion of the generator rotor angle. Stability depends on both the initial state and the severity of the disturbance.Usually the system is altered so that the post-disturbance steady state operating point differs from that prior to the disturbance. 219 f, Rotor angle, rad time (s) > (a) A 1i A 2 DO 5 B 5 o CZ. C D time (s) ’ (b) Figure B]. Response of a four-machine system during a transient [57] (a) stable system; (b) unstable system 8.4 Voltage stability and voltage collapse A power system at a given operating state is voltage stable, if an injection of reactive power causes the voltage magnitude at this node increases and all other voltages also increase or at least do not decrease. A system enters a state of voltage instability if a 220 disturbance. increase in load demand, or a change in system condition causes a progressive and uncontrollable drop in voltage. APPENDIX C DEFINITION OF HOPF AND SADDLE NODE BIFURCATIONS WITH TRANSVERSALITY CONDITIONS Appendix C Definition of Hopf and Saddle Node, Bifurcations with 'Ii'ansversality Conditions Hopf and Saddle node bifurcations have been recognized as the two generic bifurcations in variety of power system models, and are generally detected by monitoring the eigenvalues at a certain operating point, as a particular parameter of interest is varied slowly. Saddle node involves the coalescence and disappearance of two neighboring equilibrium points, while Hopf involves the emergence of a small amplitude periodic orbit from an equilibrium as the parameter is varied (a1) Saddle Node Bifurcation Let x = f(x,u) be a system of diflerential equations in R" depending on the single parameter p. At 11 = u‘ , assume that there is an equilibrium point X: (2(tt"), u‘) for which the following hypotheses are satisfied (SNI) fx = af/axam‘), if) has one simple eigenvalue Mir“) = 0 with right eigenvector v and left eigenvector w. fit all other (n-l) eigenvalues with negative real part (5N2) w- [Bf/3u(i(u‘). u')] #0 (SN3) w- [(af/allb’cw'). ll'))(v. V)] #0 Then the system undergoes a saddle node bifurcation from Y: (2(u"), u“) 221 37') (a2) Hopi-Andronov Bifurcation Let x = f (x, u) be a system of dzfierential equations in R" depending on the single parameter it. At ii = u‘, assume that there is an equilibrium point Y= (201‘), u‘) for which the following hypotheses are satisfied (HI) f; = af/ax(i(p.‘),u‘) has a simple zero eigenvalue Mu") = O and no other eigenvalues with zero real part. (H2) £1(Re7t(u))|u=p,¢0 (Transversality condition) Then a small amplitude non-constant periodic orbit of x = f(x,l,t) emerges from 7= (301‘). 11‘) at 11 = 11‘ APPENDIX D INPUT DATA FILES FOR THE BIFURCATION SUBSYSTEM EXAMPLES APPENDIX D INPUT DATA FOR BIFURCATION SUBSYSTEM EXAMPLES EXAMPLE 1 (SN) 09 The single machine infinite bus system for power system % stabilizer design % datasmib.m % bus data format % bus: number, voltage(pu), angle(degree), p_gen(pu), % q_gen(pu), p_load(pu), q_load(pu), gshunt, bshunt, % bus_type % bus_type - 1, swing bus % - 2, generator bus (PV bus) % - 3, load bus (PQ bus) % the voltage on bus 2 is adjusted for zero Q at machine 1 bus = [ l 1.05 0.0 0.7 0.4 0.0 0.0 0.00 0.300 2; 2 1.08103 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1; 3 1.0 O 0 0 0 0 0.0 0.0 3]; % line data format % line: from bus, to bus, resistance(pu), reactance(pu), % line charging(pu), tap ratio, phase shifter angle line = [l 3 0.0 0.61 0. l. 0.; 2 3 0.0 0.1000 0. 1. 0.; 2 3 0.0 0.1000 0. l. 0.]; % Machine data format leakage reactance x_l(pu) r_a(pU) sychronous reactance x_d(pu) transient reactance x'_d(pu) subtransient reactance x'_d(pu) open-circuit time constant T'_do(sec) open-circuit subtransient time % machine: 1. machine number % 2. bus number % 3. base mva % 4. % 5. resistance % 6. d—axis % 7. d-axis % 8. d-axis % 9. d-axis % 10. d-axis % constant T‘_do(sec) % 11. q-axis sychronous reactance x_q(pu) 223 224 % 12. q-axis transient reactance X'_q(pU) % 13. q-axis subtransient reactance x"_q(pU) % l4. q-axis open-circuit time constant T'_qo(sec) % lS. q-axis open circuit subtransient time % constant T'_qo(sec) % 16. inertia constant H(sec) % 17. damping coefficient d_o(pu) % l8. dampling coefficient d_1(pu) % 19. bus number % 20. S(1.0) - saturation factor % 21. S(l.2) - saturation factor % note: machine 1 uses mac_sub model, machine 2 uses mac_em % model mac_con = [ l 1 991 0.15 0 2.0 0.2 0.2 5.0 0.031 1.91 0.2 0.2 0.66 0.061 2.8756 0.0 0 1 0 0, 2 100000 0.00 0 0. 0.01 0 0 0 0 O 0 0 3.0 2.0 0 2 0 O], Exciter data format exciter: 1. exciter type - 3 for ST3 2. machine number input filter time constant T_R voltage regulator gain K_A voltage regulator time constant T_A voltage regulator time constant T_B voltage regulator time constant T_C . maximum voltage regulator output V_Rmax 9. minimum voltage regulator output V_Rmin 10. maximum internal signal V;Imax 11. minimum internal signal V_Imin 12. first stage regulator gain K_J 13. potential circuit gain coefficient K_p 14. potential circuit phase angle theta_p 15. current circuit gain coefficient K_I 16. potential source reactance X_L 17. rectifier loading factor K_C 18. maximum field voltage E_fdmax l9. inner loop feedback constant K_G 20. maximum inner loop voltage feedback V_Gmax exc_con = [ 3 1 O 7.04 0.4 0.00 1.0 7.57 0 0.2 -0.2 200 4.365 20 4.83 0.091 1.096 6.53 1 6.53]; (DQOWU'Iobw deddedePo‘Po‘PoPdeOdePdePdePdePo‘PoPo‘PdOdPON EXAMPLE 2 (Hep!) % The singl machine infinite bus system for power system stabilizer design datasmib.m W09 bus data format bus: number, voltage(pu), angleidegree), p_gen(pu), q_gen(pu), p_load(pu), q_load(pu), gshunt, bshunt, bus_type 090909019 % bus_type - 1, swing bus % - 2, generator bus (PV bus) % - 3, load bus (PQ bus) % the voltage on bus 2 is90djusted for zero Q at machine 1 % # V ang Pg Qg Pl Ql Gsh Bsh type bus = [ l 1.05 0.0 0.90 0.80 0.0 0.0 0.0 0.0 2; 2 1.08103 0.0 0.0 0.0 0.0 0.0 0.0 0.0 l; 3 1.0 0 0 0 0.7 0.4 0.0 0.0 3], % line data format % line: from bus, to bus, resistance(pu), reactance(pu), % line charging(pu), tap ratio, phase shifter angle % fr to r x chg tap pshft line = [1 3 0.0 0.6 0 1 0.; 2 3 0.0 0 01 0 1 0.; 2 3 0.0 0 01 0 1 0.]; % Machine data format % machine: 1. machine number % bus number % base mva leakage reactance x_l(pu) resistance r_a(pu) d-axis sychronous reactance x_d(pu) . d-axis transient reactance x’_d(pu) d-axis subtransient reactance x'_d(pu) d-axis open-circuit time constant T'_do(sec) d-axis open—circuit subtransient time constant T'_do(sec) 11. q-axis sychronous reactance x_q(pu) 12. q-axis transient reactance x'_q(pu) 13. q-axis subtransient reactance x'_q(pu) l4. q-axis open—circuit time constant T'_qo(sec) 15. q-axis open circuit subtransient time constant T‘_qo(sec) l6. inertia constant H(sec) 17. damping coefficient d_o(pu) 18. dampling coefficient d_l(pu) 19. bus number H OKDmQO‘lUlewN oPo‘PdOoWdPoWoPdfldPonPdPoPdePdePdPoPo‘OWdP 20. S(1.0) - saturation factor 21. S(1.2) - saturation factor note: machine 1 uses mac_sub model, machine 2 uses mac_em model mac_con = [ l l 991 0.15 0 2.0 0.245 0.2 5.0 0.031 1.91 0.42 0.2 0.66 0.061 2.8756 0.0 0 l 0 0; 2 2 100000 0.00 0 0. 0.01 0 0 0 0 0 0 O 0 3.0 2.0 0 2 0 0]; % Exciter data format % exciter: l. exciter type - 3 for ST3 % 2. machine number % 3. input filter time constant T_R o‘PoOdPanf’o'PopoPoOo'PoPoP o‘Po‘PdePo‘P exc_con = [ 004mm» 9. 10. ll. 12. 13. 14. 15. l6. 17. 18. 19. 20. 226 voltage regulator gain K_A voltage regulator time constant T_A voltage regulator time constant T_B voltage regulator time constant T_C maximum voltage regulator output V_Rmax minimum voltage regulator output VflRmin maximum internal signal V_Imax minimum internal signal V_Imin first stage regulator gain K_J potential circuit gain coefficient K_p potential circuit phase angle theta_p current circuit gain coefficient K_I potential source reactance x_L rectifier loading factor K_C maximum field voltage E_fdmax inner loop feedback constant K_G maximum inner loop voltage feedback V_Gmax 3 1 0.00 1.9 0.19 3.7 0.08 7.57 0 0.2 -0.2 400 4.365 23 4.83 1.0 1.096 6.53 2.949 6.53]; APPENDIX El FIGURES FOR DIAGNOSTIC ANALYSIS POWER TRANSFER MODEL APPENDIX E1 FIGURES FOR DIAGNOSTIC ANALYSIS POWER TRANSFER MODEL T V T V I Fowor Transfor Model glocnl load-O arc-0.4 It) . x I gains-0.0 o .— § hint—1.0 "0'10 0.2 o 4 0.6 0.3 1 Pinf (pu) (a) 0-4 '. r I .6 :Powor Transfdr Model _ 0.3 ,. .............. ............ ........... no“. '0..d.-o .. . . . . . 0.2 - » “x7 o _ . _0'1 _ .............. .................................................................... _0‘2 .. ........... ................................................. ............... -O.3 .. ............. I .......................................................... —o.4 l 1 z 1 a O 0.2 0.4 0.6 0.8 1 Pint (pu) (b) Figure 81.1 Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter Ks . (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. 227 K4 K4 228 r r I l I 2.5._..Qinfa°_.7 ................ ;""W"T ......... l A.” ............... ... 15_..........;, ........3. ............. ................ ; ............. _. . flooalload-Oé gxesOA QSL .............. ............... ................ ............... . Qinf=0.0 ; 2'5 T. l E ! fi 15.. .. . .. .. ............ ant-'51JO'H" Power Transfer Model : Elocal load-O i 0.5 i— .................... . .......... 3.x..e'.’.o...4. ............................ . ......... - o 1 l l J A O 0.2 0.4 0.8 1 0.6 Pint (pu) (b) Figure E1.2 Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter K4 . (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. 3 f I I T ' pint-0.7 2 5 _ _ Power Transfer Model .- local load-O 3 Eire-0.4 ' 2 ~- ............................ _ <0 ‘2‘. g 1.5L .......... . ’f SE 3 QIDI80.0 ‘ ._ ........................................... . ............... q 0.5 - . . ...... .. 0 l l 4 l l O 0.2 O 4 0.6 0.8 1 Pint (nu) (a) l :6 l i r l I 1 4 Power Transfer Model Glut-0.0 ' local load-O ' ‘ ire-0.4 1 .................................................... / -.. .............. .. g S , Oink—1 .0 $0.8 ............ ............... ... | . SE ; 0.6 . . . — , . . . ............................................... .. 0‘4 . ............. § ................................................ .. 0.2 .......................... . . .............. ............................................... q o L 1 1 l l 0 O 2 0.4 0.6 O 8 1 Pint (w) (b) Figure E13 Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter K1 -K2K5/K5. (a) Real and reactive power flow in the same direction. (b) Real and reactive power flow in opposite directions. 2M) I I I T I V . ............... Fdwé}. Tr'a'ns‘fe'r'M'o’d-el ...... . . . . ......... '1 ~ : f rocaI load-0 f ' \\ ; Exec-0.4 ' : : : §oinr.0.o 0.5.... .m g ............... .L .............. _. K1-K2K3K4 . . . . ._ .............. _. ................ .-.. . . . a ; . . . . ’ ‘ Oil ' -10 , n . _ .5 i i i i i O 0 0.2 0.4 0.6 0.8 1 Pint (pu) (a) ‘ r f T I I 0,9f- ....... . ...... Poww.Tr.Mr.M°dd ....... ., ................ g.nt-..-.1:.o.....( Flocat load-O ; i i e O 0.2 0.4 0.6 0.8 Pint (pu) (b) Figure E1.4 Power Transfer Model: Effect of real and reactive power transfer to the infinite bus on the parameter K1 ~K2K3K4. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. APPENDIX E2 FIGURES FOR DIAGNOSTIC ANALYSIS LOCAL LOADING MODEL APPENDIX E2 FIGURES FOR DIAGNOSTIC ANALYSIS LOCAL LOADING MODEL K4 f r : 3 _O.1’- ............ ... ..... u ..:... .... .:. ....1... 2 .. .. ...d —O.15‘- ............... .............. ;. ..... _3 ..... LocalLoadModel : Pinte-Qint-O ‘ are—OB 1 1 1 4 1 “0'20 0.2 0.4 0.6 0.8 1 Fish (pu) (a) -o.5)-- .. .. . .. .. 1 5 ... .............................. . ........... . .............................................. .... - . . . . . K4 fiPinr-Qinr-of : . I _235- ............... "-09 ............... ............... ........... 4 _2_ ............... Leca.|.l.-9e.d.Mec19! ....... .............. ........... -3p--~ -- .3. XSh--—O.6 _3_5_.. , ...... _ . .. . ..2. .L _4 1 1 1 1 1 O 0.2 0.4 0.6 0.8 1 Fish (pu) (b) Figure E2.1 Local Loading Model: Effect of the shunt local load on the parameters K... (a) Inductive loading, (b) Capacitive loading 231 K5 232 I I I I l : ; LocalLoadModel 3 oo25......... .............. ... . ore-0.9 . 002... ..... ......... ........... .............. ............. .4 0.015? ............. : . ......... : ....... .....E..... ........é .......... ............... —1 0.01,. .. . . ..... .. 0.005 stem-0.1 —o.oos 1 1 1 1 1 o 0.2 0.4 0.6 3 0.8 1 0-5 T. r r. 1 ! 0.46~-~~ .............. ........... §.x.h._o,a..... 0.4» .......... ,...,......Locglteadqual ..... ...... i Pint-Qiht-O : § ; 0_35 .. .. ............. E ........ x..0_9§- ................ ............... .3 ............... ............... _ 0.3L. .............. -............. .. ............. . ....o: ................ ............... 1 gozs- ......... .............. ............... _ o_2L. ............... ............... ...... ........... ............... —1 Q1. .............. 5 ............... 3. . . ................ Xeh-—10_ Figure E2.2 Local Loading Model: Effect of the shunt local load on the parameter K5, (a) Inductive loading, (b) Capacitive loading 233 0.12 l I I iLooalLoadMiodel 3 . 0'1? ............... P'nmn'fl. ......... .......... ......1 19-0.9 I 0.04,. ...................... . ............................................................... 0.02? ....................................... 3................................: ............... ___...——— i EXsh-o1 0... ................................................ ................................. C' .............. . 0 02 0.4 . 03 1 o_5_....... ................ ..... ...........; ....... .1 ............... .1 LocalLoadModel ; g / ; 0.4l—. Wxa-D.9M'. ............... ............... r' ............. ............... DAL. ........... ............. ................ ............... ... (b) Figure E23 Local Loading Model: Effect of the local load on the parameters K1 -K2K3K4 , (a) Inductive loading, (b) Capacitive loading Kl—K2K5/K6 Kl-KZWKG 0.1 L- . ................ - PinIf-Qinf-O XO-O.9 0.05 " .. ./ 3 z 2 2 1 1 1 1 1 O 0.2 0.4 0.6 0.8 1 95" (W) (a) " XShI—I .0 —4 —5 l») b.) A I I I .Lesell-eedMede'f ........... ....... .......... I I ................................ .............................................. 3. . . . . . .. Xsh--0.1 3 l i (b) Figure E2.4 Local Loading Model: Effect of real-inductive loading on the parameters K1-K2K5/K5 (a) Inductive loading, (b) Capacitive loading APPENDIX E3 FIGURES FOR DIAGNOSTIC ANALYSIS COMBINATION MODEL APPENDIX E3 FIGURES FOR DIAGNOSTIC ANALYSISLOCAL LOADING MODEL r ' T l 7 :Combination Model _ " \, 1 § Xsh-O.1 7 3 ’ . ‘ . 5"“0-4 . : _0.04L .............. ............... . . .. - ...-......é ................ 2: ............... - K5 E 1 1 "II-1.0 _o.061.... ............ ........ . . .. . Q . .......................... —0.08 I ................................................. . . I h I I I 0.12... ............. : .......... ................ ; ........................... _ ............... _. _0.14L..... ....... . ........ ........... ...;.QIDIQ0.0...._ 1 '1 i i 1 O 0.2 0.4 0.6 0.8 1 Pinf(pu) (a) r I I I I .f ........... ....sqwoereww .......... .............. - , §x:n-o.17 _ _o_05...... ........ . “-0.4, ................ .......... _, .............. .. ; x QIHI-0.0 - d _o.25._ ............... .. .......... . ..... 2................3................S ............... _( 0.6 Pint (pu) (b) Figure E3.1 Combination Model: Effect of power transfer on the parameter K5 when the local load is real-inductive. (a) Power flow in the same direction. (b) Real and reactive power flow in opposite directions. 235 2 1 ! I 1 I Combination Model ..Rshe-Q..65. . i ............. _ ‘ Xshs—OJ 7 (a) 2 r 1 r r I 1 . Combination Model: 1- ..;.. --~~-Fish=-0-.65 ........... .............. .. ' ' Xsh-—0.1 7 0 _ XO-O.4 . . _ _1 l' < (b) Figure E32 Combination Model: Effect of power transfer on the parameter K5 when the local load is real-capacitive. (a) Powerflow in the same direction. (b) Real and reactive power flow in opposite directions. K4 1 237 .5 r f F .' ' :Cornblnatlon Model fish-0.65 ' §Xsh-o.1 7 30-0.4 1 ’- .1 0.5-... .....E...,t. ....§.-.. : .: ............... d ant-t .0 o ; ; g 1 ; 0 0.2 0.4 0.6 0.8 1 Pint (pu) (a) 1 -5 I r 1 l f . 1 g E OWN-0.0 1 .. ..... HH'HEHQAih'f-H—‘OIBI q :5 . . Combination Model chh-O.65 . 3 0.5....”H. XSh-O.17. ................ .............. 4 gxo-OA . O O l I l 1 1 O 0.2 0.4 0.6 0.8 1 Pint (pu) (b) Figure 63.3 Combination Model: Effect of power transfer on the parameter K4 when the local load is real-inductive. (a) Power flow in the same direction, (b) Real and reactive power flow in opposite directions. t0 '9.) OC 30. . . . . ., f. . . . .. ,,§.Qinr-o.a.... .................. 1 ............ f-Qinf-tt'.0""‘ ECombinatlon Model S :Rsh-0,65....g ............... .. szh--o.t7 ' 4 _ p E Site-0.4 o l L i l l O 0.2 0.4 0.6 0.8 1 Plnf (pu) (a) 0 1. T ' 3 ! Combination Model Rah-0.65 Xsh.-_0.17 . , _1o_............. . ................ ........... ...- - anfI—LO x _15.. ......... q _20. ..... q "flu-0.1 _25 ; t L 1 i 0.2 0.4 0.6 0.8 1 Pinf(pu) Flgure E3.4 Combination Model: Effect of power transfer on the parameter K4 (b) when the local load is real-capacitive. (a) Power flow in the same direction, (b) Real and reactive power flow in opposite directions. Kl-K2K5/K6 Kl-KZKSIKB 239 ‘ form-0.0 9h- --------------- I .............. l ................ t.................:. .............. j ............... q 8"""" ....... Co.mbin.ationModel ......... ....... ....... ........ ........ ,......— 7_ .......... ....fFish.-Q-.6§....§ ................ 3. ............ ........... _ :XSh-0J7 ‘ ' 6,. . fire-0.4 a Sr... ................ _l 4i. ........................................................ Oi "181......0... 3... ....................................... _://./. ............................... fl 2.. ................................. _,/ . .............................. .. ‘LPW.."_______:—:—::"”r" ..E ....... . ........ i................é ................ i ............... Cl 0 i i L J j 0 0.2 0.4 0.6 0.8 1 Pinf(pu) (a) 30 r r r r r ‘ EON--03 Combination Model _ . § EX.h-O.17 ' Ezra-0.4 20,. .:.. . ...... .. 15'- -4 10w » « 5'- r o L L O 0.2 0.4 0.6 0.8 1 Pint (pu) (b) Figure 53.5 Combination Model: Effect of power transfers to the infinite bus on the synchronizing term K, -K2K5/K6 (regulated) when the local load is real-inductive. (a) Real and reactive power flow in the same direction. (b) Real and reactive power flow in opposite directions. Kl-K2K5/K6 240 I I l 1' I 35 iComblnatlon Model . :Rsh-O-GS ‘ ._ . a . l . 3_ Xsh-—017 , - PXB-OA’ . 2_5 ,. : .,. . “4 0mm .0 g Q Q J. x ngnt-O.3 ..1 i . 0 0.2 0.4 0.6 0.8 1 Plnf (pu) (a) 1 F .T I T ! :Combinatlon Model gRsh:O.65 _ Z 05? ....... Xsh-v-Otz ............... ............... _ xe-O.4 _ ‘ faint-4 .o _15 L ; ; ; t 0 0.2 0.4 0.6 0.8 1 Pint (pu) (b) Flgure E3.6 Combination Model: Effect of power transfers to the infinite bus on the synchronizing term K, -K2K5/K6 (regulated) when the local load is real~capacitive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. 241 1.5 r i I I ‘ 1_ _ _ j .. x 'g .:.Qinf-0.0 a Kl -K2K3K4 Combination Model _o_5-._...,RshuO‘;65 .......... _ ,. ... _. f, ; XSh-O.17 OWN-1.0 XB-O.4. ‘ ‘ ‘ _1 a 1 z i i 0 0.2 0.4 0.6 0.8 1 Pinf (pu) (a) T 3 f r ‘. geomblnation Model 1 RSh-O.65 .Xsh-O.17 ,xe-OA , . r 3 o’er ............... g .............. ............ ........ .......... ......§.Q'.'.‘_".'.‘.9-§.._ 0.6 L- : i ........ - 0.4... .................................... - ‘4 .4 f . : : ; 0.2... ..... W ................ ................ .” .............. _, /” j ‘ ‘ Qtnf=0 0 OF. ................ .................................................................. a O 0.2 0.4 0.6 0.8 1 Pint (pu) (b) Figure E3.7 Combination Model: Effect of power transfers to the infinite bus on the synchronizing term K, -K2K3K4 (unregulated) when the local load is real-inductive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. K1-K2K3K4 242 ‘ ' V T l 1 '8 - ............ - ... .......... . e . . . i ........... - . s . E ........ 'o - v ... o . e ’M. ............... Eu - . u , .......... d ; 3 ;Combination odel ; ° 0 0.2 0.4 0.6 0.3 1 Pinf (pu) (a) F Y t ', F 1.2i— ............... ;,...............é................E................é ................ é. .............. _1 Contain-tion Model : 5 Q,fl,__,_o {Fish-0.65 i ' 1 1n- .............. :xghro.’.17. ............................................................ .- :2 : g g '1 0.6» s ..... .. a: ? : pink—0.1 o_4.. ...... . A ..., i ”‘/ /' é ; s 2 o 4 ; t i i O 0.2 0.4 0.6 0.8 1 Pinf (pu) (b) Figure E33 Corrbination Model: Effect of power transfers to the infinite bus on the synchronizing term K, -K2K3K4 (unregulated) when the local load is real-capacitive. (a) Real and reactive power flow in the same direction, (b) Real and reactive power flow in opposite directions. K5 0.15 0.1 0.05 -0.05 -O.1 0.8 243 I V I l I . fiPlnt-oa Quit-0.5 _ ............... , ................ . ................ f-XSh-O.~1----- . 1 l l _L O 0.2 0.4 0.8 1 0.6 Rsh (pu) (a) 7 7 Y f Y : 3 Combination Model _ 0.7- ............... .............. iPinf-O.8..Oirr:f-0.5 .......... g ............... _ - 5 bro-0.4 ' 3 QGPH. . ............. .......... . ..... ..... . ......... _ ; ; ; : XSh'-°-‘ 0.5.. ............... .............. ................ ................ ............... _ 0'4. ........... _ 0.3)....”W .............,:.... ..; ........... ......... ............... .4 0.2,. ....................................... ................................ E ............... _, fiXSh-—1.0 0.1 l l I l l. O 0.2 0.4 0.6 0.8 1 Rsh(pu) (b) Figure E3.9 Combination Model: Effect of local loading on the parameter K5 (a) (regulated) when the transfer to the infinite bus is at P,,,,=O.8, Q...=0.5. real-inductive loading; (b) real-capacitive loading. i ' ! F ! 1 .4 _. ............... ........ . ..... : ........ . ....... Combination Mbdfil .......... ............... .4 ' gPinr-oe cameos {XS-0.4 1.2” : " ' ‘ SXSh-fl .0 ............... 0'20 o 2 o 4 0.6 o e 1 Fish (pu) 3 i I f r l ‘x’ 0.. .................................................................................. ._( _1L.... ..... CombinauonModei ............... _ Pintfi-Ofl aim-0.5 ' ’ ' xe-0.4 S : t i ' ' i 3 X8h--O.4 -2- ............. ..... . ....... .............. 3 ............... _3 1 1 l l l O 02 0.4 06 08 1 Figure £3.10 Combination Model: Effect of local loading on the parameter K4 (unregulated) when the transfer to the infinite bus is at Pim=0.8, Q...=0.5. (a) real-inductive loading; (b) real-capacitive loading. 245 I , . :XSh-LO 0.95 ._ ............ ....... V ...... ; ............... $ ,,,,,,,,,,,,,, ~ ............ _( 0.85” ............. ...... _ .. _ ................ ............... .4 3 / '/ i i :Xsh-Q1 .0 O l i K1—K2K3K4 .0 \l on / goombination Model / ' " ' fipinreo‘e amigos“ Q " "" 2 ; EKG-0.4 2 2 0.65). ............. ................ ..... . .......... x ........... r. .; ............... —: .0 V 1 °~55F ............. ................ g ................ ................ ............... _ 05 i i i L 1' 0.6 Fish (pu) (a) 4 1 V I U I 3.5 ~ Ram—0:4 . . . A combination Model 2 3 b .............. germane. . QQfHO.-5. .......... ............... g .............. g ............... 4 xe-O.4 I f 2.5 p ........... _ 1 . ................ g ................ .............. j ................ § .......... , . . d K1-K2K3K4 Xsh-—1.0 1 .. ........................................................................ I ............... .1 0.5). ........... . ................ ................ Z ............... -1 (b) Figure E3.1 1 Combination Model: Effect of local loading on the parameter K1 -K2K3K4 (unregulated) when the transfer to the infinite bus is at PiM=0.8. Q...=O.5. (a) real-inductive loading; (b) real-capacitive loading. K1-K2K5/K6 10 246 I T l I I ,........ ..,...b°mb|naflgnMode| ..... .. ................ ............. ...... . ........ .4 Pint-0.8 emf-0.5 ' ' q E :Xsh-OJ .......................................................................................... O 0.2 0.4 0.6 0.8 1 Fish (pu) (a) Kl-K2K5/K6 / / _2_ ............. Piaf-0.8..Qini.o.s. .......... .............. g ............... a Eire-0.4 ' ' ' ' 3 : : Xsh-—04 _3,_ ............... ........ . ...... . ............... ............... -4 _4 i ; i 4' ; Figure E3.12 Combination Model: Effect of local loading on the parameter K1 -K2K5/K5 (regulated) when the transfer to the infinite bus is at Pw=0.8, Q...=O.5. (a) real-inductive loading; (b) real-capacitive loading. 247 Ks 1.2P—"" .... ..E .. ..... . , ..: ....... . ....... ,.....W ... .4 ............... x4 Sh...1.'Q-I 1i— ............................................. _ 0.8- . ct 0.6? ................. ........................................................... ... 0'4 ' gramme ‘ KIM-0.5 i : ate-0.4 . 02..., ................. . ....... ,... ....... _( 4 i L L ' O 02 0.4 0.6 08 1 RSh(PU) (a) 0.1 f . r . i 0.08-. ...... .. . ..é... . ....... Pin-f;o-I‘8 ........ I ................. ....... flint-0.5 0.06 0.04 —0.02 -0.04._. - .' ............. . .............. Xs‘fiiO”.1'-l _0-06 1 1 l 1 1 O 0.2 0.4 0.6 0.8 1 Fish (pu) (D) Figure £3.13 Corrbination Model: Effect inductive-real local load on the system synchronizing and damping torque coefficients when Pm=0.8, 05“,:05 (a) synchronizingforque coefficient KS, (b) damping torque coefficient KD 248 Figure £3.14 Combination Model: Damping torque coefficienth and parameter K5 for different levels of real local loading BIBLIOGRAPHY [l] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] REFERENCES P. Kundur. G. J. Rogers. D. Y. Wong. L. Wang. and MG. Lauby. “A Comprehensive Com- puter Program Package for Small Signal Stability Analysis of Power Systems”. IEEE Transactions on Power Systems. Vol. 5. No.4, pp. 1076-1083. Nov 1990 B. Gao. G. K. Morison, and P. Kundur, “Voltage Stability Evaluation Using Modal Analy- sis", IEEE Transactions on Power Systems. Vol. 7. No.4, pp. 1529-1542, Nov 1992 Bruce C. Moore, “Principle Component Analysis in Linear Systems: Controllability, Observability. and Model Reduction”, IEEE Transactions on Automatic Control, Vol. AC- 26, No. 1. Feb 1981 Kokotovic, P.V.. Allemong, J. J., Winkelman, J. R., and Chow, J. H., “Singular Perturba- tion and Iterative Seperation of Time Scales”, Automatica, Vol. 16: pp. 23—33. Jan 1980 I. Dobson. “Observations on the Geometry of Saddle Node Bifurcation and Voltage Col- lapse in Electric Power Systems”, IEEE Transactions on Circuits and Systems, Part I, 39(3):pp. 240-243. March 1992 Chhewang Rinzin. Fernando L. Alvarado, and Rambanu Adapa. “Formulas for Geo- graphic Differentiation of Security Limits”, 24th American Power Symposium. Oct 1992 Yakout Mansour. Wilsun Xu, Femando Alvarado, Chhewang Rinzin. “SVC Placement Using Critical Modes of Voltage Instability”, IEEE Transactions on Power Systems, Vol. 9, No. 2. May 1994 “MSU Voltage Stability Security Assessment and Planning Methodology Programs”. Department of Computing and T ed“ aology, ‘rlichigan State Univ" East Lansing, 1994 R. A. Schlueter, T. Y. Guo. T. Lie, F. Li. J. Ashley. “Voltage Collapse Security Assessment on the B. C. Hydro System”. Division of Engineering Research, Report PSL-BCI to B. C. Hydro, March 1991 “SVC for Dynamic Voltage Stabilization of 132 kV System in Western Canada”, ABB Advertising Application Note A02-0144 E, ABB Power Systems. 1994 R. A. Schlueter, K. B. Kilani. S. Liu, “A Voltage Stability Security Assessment Method", submitted to IEEE Working Group on Voltage Collapse for IEEE Tutorial, Jan 1996 Perez-Arriaga, G.C Verghese, F.C. Schweppe, “Selective Modal Analysis With Applica- tions to Power Systems, Part I: Heuristic Introduction”, IEEE PES I 982 Winter Meeting, 249 [13] [14] [15] [16] [17] [18] [19] [20] [21} [22] [23] [24] [25] 250 New York Vittal. V.. N. Bhatia. and A. A. Fouad. “Analysis of the Inter-area Mode Phenomena in Power Systems Following Large Disturbance”. IEEE Transactions on Power Systems. (6):pp. 1515-1521. 1991 E. H. Abed and P. P. Varaiya. “Nonlinear Oscillations in Power Systems”, Electric Power and Energy Systems, Vol. 6:pp. 37-43. No. 1. January 1984 C. Rajagopaian. P.W. Sauer and MA. Pai. “ Analysis of Voltage Control Systems exhibit- ing Hopf Bifurcations.” Proceedings 28th Decision and C0ntrol Confi, Tampa, FL. pp. 332-335. December 1989. J.H. Chow. A. Gebreselassie. “Dynamic Voltage Stability Analysis of a Single Machine Constraint Power Load System", Proceedings of the 29th Conference on Decision and Control, Honolulu. Hawaii. Dec 1990 Venkatasubramanian, V.. H. Schattler and J. Zaborsky. “A taxonomy of the Dynamics of a Large Power System with Emphasis on its Voltage Stability”, Proc. Bulk Voltage Stability and Security, Deep Creek Lake, Aug, 1991 C. Rajagopalan. B. Lesieutre. P. W. Sauer. M. A. Pai, “Dynamic Aspects of Voltage / Power Characteristics", IEEE Transactions on Power Systems. Vol. 7. No. l3zpp. 990- 1000. Aug 19992 F. Alvarado. M.Gu. Y. Hu. “Direct Detection of Instability Points Using an Augmented Formulation”. Proceedings of the 24th North American Power Symposium. Reno, Nevada. Oct. 1992 Lee. B. and V. Ajjarapu. “Observation of Chaos in an Electrical Power System Via Period Doubling”, IEE Proceedings of Parc C. Vol. 140. No. 6. pp. 490-496. Nov, 1993 Ajjarapu, V. and B. Lee, “The application of Bifurcation Theory to Study the Nonlinear Dynamical Phenomena in an Electric Power System”, IEEE Transactions on Power Sys- tems. Vol. 7. No. 1. pp. 424432. Feb 1992 Abed. E. H., J. C. Alexander, H. Wang, A. M. A. Hamdan. and H. C. Lee. “Dynamic Bifurcation in a Power System Model Exhibiting Voltage Collapse". Technical Research Report, System Research Center. U. Maryland, College Park. MD, U.S.A.. Feb 1992 Chiang, H. D., C. W. Liu. P. P. Varaiya. F. F. Wu, and M. G. Lauby, “Chaos in a Simple Power System". IEEE PES Winter Meeting Paper 92 WM 1 5 1 -I PWRS Chiang, H-D, 1. Dobson. RJ. Thomas. J.S. Thorp. F.A. Lazhar, “On the Voltage Collapse in Power Systems”. IEEE Trans. on Power Systems. 5:pp. 601-611, May 1990 Robert A. Schlueter, Khadija B. Kilani, and Uhi Ahn. “Impact of Modeling Accuracy on Type, Kind, and Class of Stability Problems in a Power System Model”, Proceedings of the ECC & NSF International Workshop on Bulk Power System Voltage Stability, Security and Control, Phenomena — Ill. pp. 117-156. Aug 1994 [26] [27] [28] [29] [30] [31] [33] [34] [35] [36] [37] [38] [40] [41] 251 V. Venkatasubramanian, X. Jiang. H. Schttler, and J. Zaborszky. “A Taxonomy of the Dynamics of the Large Electric Power System with Emphasis on Its Voltage Stability”. Proceedings of the NSF International Workshop on 8qu Power System Voltage Phenom- ena - II. pp. 9-52, 1991 V. Venkatasubramanian. X. Jiang. H. Schttler, and J. Zaborszky. “Cunent Status of the Taxonomy Theory of Large Power System Dynamics. DAE Systems with Hard Limits” Proceedings of the ECC & NSF International Workshop on Bulk Power System Voltage Stability, Security and Control, Phenomena - Ill. pp. 15-103. Aug. 1994 T. Van Cutsem. Y. Jaquemart. J.-N. Marquet. P. Pntvot.. “A Comprehensive Analysis of Mid-term Voltage Stability”, IEEE Transactions on Power Systems. Vol. 10, No. 3:pp. 1173-1182. August 1995 P. Kundur. Power System Stability and Control. Vol 1 of the EPRJ Power System Engi- neering Series, McGraw Hill. New York 1994 W.R. Lashs. “Control of voltage stability on EHV power systems”. Ph.D. Dissertation. University of South Wales. Australia, May 1992 BC. Lesieutre, P.W. Sauer, M.A. Pal, “Development and comparative study of P, Q load models”. IEEE Winter Power Meeting, Paper 94 WM 166-9 PWRS. Hassan K. Khalil, “Nonlinear Systems”, Macmillan Publishing Company. New York. 1992 J. Guckenheimer and P. Holmes. “ Nonlinear Oscillations, Dynamical Systems and Bifur- cations of Vector Fields,” Applied Mathematical Sciences. Springer. New York, 1986 R.A Schlueter, “Unification and Classification of Algebraic Tests for Voltage Stability," Electric Mach. Power Systems.. 12. pp. 557-590, 1993 G. Taylor. Voltage Stability, Part 1. Survey of Voltage Collapse Phenomena. NERC Rep.. National Electric Reliability Council, Princeton, NJ , Aug, 1991 IEEE Systems Dynamic Performance Report. “Voltage Stability of Power Systems: Con- cepts, analytical tools. and industry Experience”, 1555 Document No. 90 TH 0358-3 PWR, 1990 [BEE Systems Dynamic Performance Report. “Suggested Techniques for voltage stability analysis”. IEEE Document No. 93 TH 0620-5PWR, 1994 NW. Miller. W.W.Price. “Influence of Load Modeling Assumptions on the Dynamics of the Large System Voltage Collapse Simulations”, Proceedings of the ECC & NSF Interna- tional Workshop on 8qu Power System Voltage Stability, Security and Control, Phenom- ena - II], pp. 635-641, Aug 1994 M. Amorouayeche. M. Brown. R. Craven, J.L. Jardim. A. Johannesen. P. Kundur. R. Mar- conato. K. Maslo, F. Mc Namara, B. Meyer. J.V. Mitsche. Y. Nakanishi. W.W. Price. J.L. [42] [43] [44] [45] [47] [48] [49] [50] [51] [52] [53] [54] [55] Sancha. K. Waive, D. Xia. R. Yokoyama. “Tools for Simulating Long Term Dynamics”. Electra No. 163, pp. 151-165. December 1995 R. A. Schlueter, I-P. Hu and T-Y. Huo. “Dynamic/Static Voltage Stability Security Crite- ria”, Proceedings of the NSF International Workshop on Bulk Power System Voltage Phe- nomena - Voltage Stability and Security. pp. 265-303, 1991 Modeling of Voltage Collapse Including Dynamic Phenomena. Final Report. CIGRE Task Force 38-02-10. December 1992 M. Amorouayeche. M. Brown, R. Craven, J.L. Jardim. A. Johannesen. P. Kundur, R. Mar- conato, K. Maslo, F. Mc Namara. B. Meyer, J.V. Mitsche. Y. Nakanishi, W.W. Price. J.L. Sancha, K. Waive, D. Xia, R. Yokoyama, “Tools for Simulating Long Term Dynamics”, Electra No. 163. pp. 151-165. December 1995 I. Dobson. L. Lu, “Immediate Change in Stability and Voltage Collapse When Generator Power Limits are Reached”. Proceedings of the NSF International Workshop on 8qu Power System Voltage Phenomena - Voltage Stability and Security, pp. 65-73, 1991. V. Venkatasubramanian, X. Jiang. H. Schattler, and J. Zaborszky, “Discussion of Impact of Modeling Accuracy on Type. Kind, and Class of Stability Problems in a Power System Model”, Proceedings of the EC C & NSF International Workshop on Bulk Power System Voltage Stability, Security and C ontrol, Phenomena - III, pp. 171-172, Aug 1994 R. A. Schlueter, “Closure of Discussion of Impact of Modeling Accuracy on Type, Kind, and Class of Stability Problems in a Power System Model”, Proceedings of the ECC & NSF International Workshop on Bulk Power System Voltage Stability, Security and Con- trol, Phenomena - III, pp. 172-173. Aug 1994 Lee. B. and V. Ajjarapu, “A Unified Frame Work to Study Voltage Stability of Structure Preserving Power System Model”. Submitted for Publication in IEEE Trans. on Power Systems S. Wiggins, “Introduction of Applied Nonlinear Dynamical Systems and Chaos”. Springer-Verlag, New York, 1990 1. Hu. “Voltage Collapse Bifurcation on a Power System Transient Stability Model”, Ph.D. Dissertation, Michigan State University, July 1990 T.Y. Guo. “Structural bifurcation tests for voltage collapse and low frequency oscillations in power systems”, Ph.D. Dissertation, Michigan State University, March 1992 FR. Gantmatcher “ The Theory of Matrices”, vol 1(2), Chelsea Publishing Company, New York, 1977. Y. Tamura. H. Mori. and S. Iwamoto, “Relashionship between Voltage Instability and Multiple Load Flow Solutions in Electric Power Systems”, IEEE Transactions on Power Apparatus and Systems. Vol. PAS-102. pp. 1115-1123, May 1983 Francisco P. Demelio and Charles Concordia. “Concept of Synchronous Machine Stability [56] [57] [58] [59] [59] 253 as Affected by Excitation Control”, IEEE Transactions on Power Apparatus and S ystems. V01. PAS-88. No. 4. pp. 316-329. April 1969 F. P. deMello, T. F. Laskowsky, “Concept of Power System Dynamic Stability”, IEEE Transactions on Power Apparatus and Systems. Vol. PAS-94. No. 3. pp. 827-833. May- June 1975 AA. Fouad and P. M. Anderson, “Power System Control and Stability”, Iowa State Uni- versity Press. Ames. IA, 1977 Sotomayor, J., “Generic Bifurcations of Dynamical Systems”, Edited by PM. Peixoto, Academic Press, NY, 1973 “Power System Toolbox”. Cherry Tree Scientific Software, RR#5 Calborne, Ontario KOK ISO, Canada. Arthur R. Burgen. “Power System Analysis”, Prentice Hall, lnc., Englewood Cltfl’s, New Jersey 07632, 1986 MICHIGAN smrE UNIV. LIBRQRIES llllllllllllllllllllWilllllHHllllllllllllllllllllllllll 31293015700853