Power System Dynamics and Stability: Multimachine Simulation Overview

ECE 576 Power System Dynamics and Stability Lecture 19: Multimachine Simulation Prof. Tom Overbye Dept. of Electrical and Computer Engineering University of Illinois at Urbana-Champaign overbye@illinois.edu 1

Announcements Read Chapter 7 Homework 6 is due on Tuesday April 15 A useful reference is B. Stott, "Power System Dynamic Response Calculations," Proc. IEEE, vol. 67, pp. 219- 241 Another key reference is J.M. Undrill, "Structure in the Computation of Power-System Nonlinear Dynamic Response," IEEE Trans. Power App. and Syst., vol. 88, pp. 1-6, January 1969. 2

Angle Reference The initial angles are given by the angles from the power flow, which are based on the slack bus's angle As presented the transient stability angles are with respect to a synchronous reference frame Sometimes this is fine, such as for either shorter studies, or ones in which there is little speed variation Oftentimes this is not best since the when the frequencies are not nominal, the angles shift from the reference frame Other reference frames can be used, such as with respect to a particular generator's value, which mimics the power flow approach; the selected reference has no impact on the solution 3

Subtransient Models The Norton current injection approach is what is commonly used with subtransient models in industry If subtransient saliency is neglected (as is the case with GENROU and GENSAL in which X"d=X"q) then the current injection is ( q d q Nd Nq s R jX R + ) d + + j + E jE + = = I jI jX s Subtransient saliency can be handled with this approach, but it is more involved (see Arrillaga, Computer Analysis of Power Systems, section 6.6.3) 4

Subtransient Models Note, the values here are on the dq reference frame We can now extend the approach introduced for the classical machine model to subtransient models Initialization is as before, which gives the 's and other state values Each time step is as before, except we use the 's for each generator to transfer values between the network reference frame and each machine's dq reference frame The currents provide the coupling 5

Two Bus Example with Two GENROU Machine Models Use the same system as before, except with we'll model both generators using GENROUs For simplicity we'll make both generators identical except set H1=3, H2=6; other values are Xd=2.1, Xq=0.5, X'd=0.2, X'q=0.5, X"q=X"d=0.18, Xl=0.15, T'do = 7.0, T'qo=0.75, T"do=0.035, T"qo=0.05; no saturation With no saturation the value of the 's are determined (as per Lecture 11) by solving ( s q E V R jX = + + ) I Hence for generator 1 1.0946 11.59 = ( )( ) + = 1.431 30.2 0.5 1.052 18.2 E j 1 1 6

GENROU Block Diagram E'q E'd 7

Two Bus Example with Two GENROU Machine Models Using the approach from Lecture 11 the initial state vector is . . . . . 2q1 0 2446 E 0 0 0 5392 0 E 0 9044 0 8928 0 3594 E 0 Note that this is a salient pole machine with X'q=Xq; hence E'd will always be zero 0 5273 0 0 1 1948 1 1554 1 1 q1 E 1d1 The initial currents in the dq reference frame are Id1=0.7872, Iq1=0.6988, Id2=0.2314, Iq2=-1.0269 d1 = = x ( ) . 2 2 q2 . . Initial values of "q1= -0.2236, and "d1 = 1.179 1d2 . 2q2 d2 8

Solving with Euler's We'll again solve with Euler's, except with t set now to 0.01 seconds (because now we have a subtransient model with faster dynamics) We'll also clear the fault at t=0.05 seconds For the more accurate subtransient models the swing equation is written in terms of the torques d dt 2H d 2H d T T dt dt T i i = Other equations are solved based upon the block diagram = = i i s i ( ) = = i i i i D Mi Ei i i s s with , , Ei d i qi q i di 9

Norton Equivalent Current Injections The initial Norton equivalent current injections on the dq base for each machine are ( 1 1 1 1 1 6.55 1.242 2.222 6.286 ND NQ I jI j I jI j + = + + = ) ( ) q d + j 0.2236 + 1.179 (1.0) 0.18 j 1 + = = I jI Nd Nq jX j j = = + + Recall the dq values are on the machine's reference frame and the DQ values are on the system reference frame 1 1 4.999 1.826 2 2 Nd Nq 1 5.227 I jI j 2 2 ND NQ 10

Moving between DQ and dq Recall I I I I sin cos cos sin di Di = qi Qi The currents provide the key coupling between the two reference frames And I I I I sin cos cos sin Di di = Qi qi 11

Bus Admittance Matrix The bus admittance matrix is as from before for the classical models, except the diagonal elements are augmented using 1 Y R jX + = i , , s i d i 1 0 . . j10 101 j4 545 j4 545 j10 101 . j0 18 = + = Y Y N . . 1 0 . j0 18 12

Algebraic Solution Verification To check the values solve (in the network reference frame) 1 . . . . j10 101 j4 545 j4 545 j10 101 2 222 j6 286 j5 227 = V . . . 1 + 1 072 = . . j0 22 . 1 0 13

Results The below graph shows the results for four seconds of simulation, using Euler's with t=0.01 seconds PowerWorld case is B2_GENROU_2GEN_EULER 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20 -25 -30 -35 -40 0 0.5 1 1.5 2 2.5 3 3.5 4 Rotor Angle_Gen Bus 1 #1 Rotor Angle_Gen Bus 2 #1 g f e d c b g f e d c b 14

Results for Longer Time Simulating out 10 seconds indicates an unstable solution, both using Euler's and RK2 with t=0.005, so it is really unstable! 26,000 24,000 22,000 20,000 18,000 16,000 14,000 12,000 10,000 8,000 6,000 4,000 2,000 35,000 30,000 25,000 20,000 15,000 10,000 5,000 0 -2,000 -4,000 -6,000 -8,000 -10,000 -12,000 -14,000 -16,000 0 -5,000 -10,000 -15,000 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 Rotor Angle_Gen Bus 1 #1 Rotor Angle_Gen Bus 2 #1 Rotor Angle_Gen Bus 1 #1 Rotor Angle_Gen Bus 2 #1 g f e d c b g f e d c b g f e d c b g f e d c b RK2 with t=0.005 Euler's with t=0.01 15

Adding More Models In this situation the case is unstable because we have not modeled exciters To each generator add an EXST1 with TR=0, TC=TB=0, Kf=0, KA=100, TA=0.1 This just adds one differential equation per generator ( A REF t A dt T dE 1 ( ) ) = K V V E FD FD 16

Two Bus, Two Gen With Exciters Below are the initial values for this case from PowerWorld Because of the zero values the other differential equations for the exciters are included but treated as ignored Case is B2_GENROU_2GEN_EXCITER 17

Viewing the States PowerWorld allows one to single-step through a solution, showing the f(x) and the K1 values This is mostly used for education or model debugging Derivatives shown are evaluated at the end of the time step 18

Two Bus Results with Exciters Below graph shows the angles with t=0.01 and a fault clearing at t=0.05 using Euler's With the addition of the exciters case is now stable 45 40 35 30 25 20 15 10 5 0 -5 -10 -15 -20 -25 -30 -35 0 1 2 3 4 5 6 7 8 9 10 Rotor Angle_Gen Bus 1 #1 Rotor Angle_Gen Bus 2 #1 g f e d c b g f e d c b 19

Constant Impedance Loads The simplest approach for modeling the loads is to treat them as constant impedances, embedding them in the bus admittance matrix Only impact the Ybus diagonals The admittances are set based upon their power flow values, scaled by the inverse of the square of the power flow bus voltage ( , , , load i i load i i load i load i S V I V G jB S G jB V In PowerWorld the default load model is specified on Transient Stability, Options, Power System Model ) 2 = = * , = , load i , , load i load i 2 i Note the positive sign comes from the sign convention on I 20 load,i

Example 7.4 Case (WSCC 9 Bus) PowerWorld Case Example_7_4 duplicates the example 7.4 case from the book, with the exception of using different generator models = Bus 5 Example: Without the load = + . 2 553 j17 339 - . Y 55 . . and V =0.996 j0 5 S 1 25 , 5 . load 5 ( ) . 1 25 j0 5 + = Y = . - . . . 2 553 j17 579 3 813 j17 843 55 2 0.996 21

Nonlinear Network Equations With constant impedance loads the network equations can usually be written with I independent of V, then they can be solved directly (as we've been doing) = 1 V Y x ( ) In general this is not the case, with constant power loads one common example Hence a nonlinear solution with Newton's method is used We'll generalize the dependence on the algebraic variables, replacing V by y since they may include other values beyond just the bus voltages 22

Nonlinear Network Equations Just like in the power flow, the complex equations are rewritten, here as a real current and a reactive current YV I(x,y) = 0 The values for bus i are ( ( ) Di ik Dk ik QK k 1 = This is a rectangular formulation; we also could have written the equations in polar form n ) = = x,y g G V B V I 0 NDi n ( ) = + = x,y ( ) g G V B V I 0 Qi ik Qk ik DK NQi k 1 = For each bus we add two new variables and two new equations If an infinite bus is modeled then its variables and equations are omitted since its voltage is fixed 23

Nonlinear Network Equations The network variables and equations are then n ( ) = x,y ( ) 0 G V k QK B V I 1 1 1 k Dk ND = 1 k n ( ) + = V V x,y ( ) 0 ik Qk G V ik DK B V I 1 D 1 NQ = 1 k n 1 Q ( ) = x,y ( ) 0 G V B V I V 2 2 2 k Dk k QK ND 2 D = = y ( , ) g x y = 1 k V V Dn n ( ) = x,y ( ) 0 nk Dk G V nk QK B V I Qn NDn = 1 k n ( ) + = x,y ( ) 0 nk Qk G V nk DK B V I NQn = 1 k 24

Nonlinear Network Equation Newton Solution The network equations are solved using a similar procedure to that of the Netwon-Raphson power flow ( ) v = y y Set 0; make an initial guess of , v ( ) v g y While ( ) Do + ( 1) ( ) v ( ) v 1 ( ) v v = = y y J y g y ( ) ( ) + 1 v v End While 25

Network Equation Jacobian Matrix The most computationally intensive part of the algorithm is determining and factoring the Jacobian matrix, J(y) ( , ) D g V = x y ( , ) x y ( , ) x y g g 1 1 V 1 V D D 1 1 D Q Qn ( , ) x y ( , ) x y ( , ) x y g g g 1 V 1 V 1 V Q Q Q ( ) J y 1 1 D Q Qn x y x y x y ( , ) V ( , ) V ( , ) V g g g Qn Qn Qn 1 1 D Q Qn 26

Network Jacobian Matrix The Jacobian matrix can be stored and computed using a 2 by 2 block matrix structure The portion of the 2 by 2 entries just from the Ybus are ( , ) ( , ) Di Di g g V V G g g B V V The "hat" was added to the g functions to indicate it is just the portion from the Ybus x y x y = B Dj Qj ij ij ( , ) x y ( , ) x y G Qi Qi ij ij Dj Qj The major source of the current vector voltage sensitivity comes from non-constant impedance loads; also dc transmission lines 27

Example: Constant Current and Constant Power Load As an example, assume the load at bus k is represented with a ZIP model ( ( , , , Load k BaseLoad k z k k Q Q Q V The base load values are set from the power flow ) = + + 2 Load k P BaseLoad k P P V P V P , , , , , z k k i k k p k ) = + + 2 Q V Q , , i k k p k The constant impedance portion is embedded in the Ybus ( ( , , , Load k BaseLoad k i k k Q Q Q V ) ( ) ( , p k ) = + = + Load k P BaseLoad k P P V P BL i k P V BL p k P , , , , , , , , i k k p k k ) = + = + Q BL i k Q V BL p k Q , , , , k Usually solved in per unit on network MVA base 28

Example: Constant Current and Constant Power Load The current is then P I I jI = + = + + = * + jQ , , Load k Load k , , , , , Load k D Load k Q Load k V k ) ) ( ( + + 2 2 2 2 BL i k P V V BL p k P j Q V V BL p k Q , , , , , , , , DK QK BL i k DK QK V jV Dk Qk Multiply the numerator and denominator by VDK+jVQK to write as the real current and the reactive current 29

Example: Constant Current and Constant Power Load + + V P V Q + V P V Q , , V , , , , V , , Dk BL p k QK 2 QK V BL p k Dk BL i k QK V BL i k = + I , , D Load k 2 + 2 2 DK DK QK V P V Q + V P V Q , , V , , , , V , , Qk BL p k DK 2 QK V BL p k Qk BL i k DK V BL i k = + I , , Q Load k 2 + 2 2 DK DK QK The Jacobian entries are then found by differentiating with respect to VDK and VQK Only affect the 2 by 2 block diagonal values Usually constant current and constant power models are replaced by a constant impedance model if the voltage goes too low, like during a fault 30

Example: 7.4.ZIP Case Example 7.4 is modified so the loads are represented by a model with 30% constant power, 30% constant current and 40% constant impedance In PowerWorld load models can be entered in a number of different ways; a tedious but simple approach is to specify a model for each individual load Right click on the load symbol to display the Load Options dialog, select Stability, and select WSCC to enter a ZIP model, in which p1&q1 are the normalized about of constant impedance load, p2&q2 the amount of constant current load, and p3&q3 the amount of constant power load Case is Example_7_4_ZIP 31

Example 7.4.ZIP One-line Bus 7 Bus 8 Bus 9 Bus 3 Bus 2 1.016 pu 163 MW 7 Mvar 85 MW -11 Mvar 1.025 pu 1.026 pu 1.032 pu 1.025 pu 100 MW 35 Mvar Bus 5 0.996 pu Bus 6 1.013 pu 125 MW 50 Mvar 90 MW 30 Mvar Bus 4 1.026 pu Bus1 1.040 pu 72 MW 27 Mvar slack 32

Example 7.4.ZIP Bus 8 Load Values As an example the values for bus 8 are given (per unit, 100 MVA base) ( , . . . . BaseLoad 8 P 0 983 = = + ) = 0 4 1 016 + 0 3 1 016 + 2 . . . 1 00 BaseLoad 8 P 0 3 , ( ) 0 3 1 016 + 2 . . . . . . 0 35 BaseLoad 8 Q 0 4 1 016 0 3 , 0 344 = . BaseLoad 8 Q , * + . 1 j0 35 j0 0129 + + = = . . I jI 0 9887 j0 332 , , , , D Load 8 Q Load 8 . . 1 0158 33

Example: 7.4.ZIP Case For this case the 2 by 2 block between buses 8 and 7 is 1.155 9.784 9.784 1.155 These entries are easily checked with the Ybus And between 8 and 9 is 1.617 13.698 13.698 1.617 The 2 by 2 block for the bus 8 diagonal is 2.876 23.632 23.352 3.745 The check here is left for the student 34

Additional Comments When coding Jacobian values, a good way to check that the entries are correct is to make sure that for a small perturbation about the solution the Newton's method has quadratic convergence When running the simulation the Jacobian is actually seldom rebuilt and refactored If the Jacobian is not too bad it will still converge To converge Newton's method needs a good initial guess, which is usually the last time step solution Convergence can be an issue following large system disturbances, such as a fault 35