Power System Stability Lecture: Transient Solutions and Load Models

This lecture delves into transient stability solutions and load models in power systems. It covers the simultaneous implicit method for solving algebraic and differential equations, focusing on the Trapezoidal approach for linear systems. The discussion extends to nonlinear cases, particularly using Newton's method for solving with the Trapezoidal method. An example involving an Infinite Bus GENCLS system with a classical model is also presented.

• Power System Stability
• Transient Stability
• Load Models
• Simultaneous Implicit Method
• Newtons Method

annabelle Follow

Uploaded on Aug 03, 2024 | 0 Views

Download Presentation

Please find below an Image/Link to download the presentation.

The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.

Presentation Transcript

1. ECEN 667 Power System Stability Lecture 16: Transient Stability Solutions, Load Models Prof. Tom Overbye Dept. of Electrical and Computer Engineering Texas A&M University overbye@tamu.edu

2. Announcements Read Chapter 7 Homework 4 is due on Tuesday Oct 29 1

3. Simultaneous Implicit The other major solution approach is the simultaneous implicit in which the algebraic and differential equations are solved simultaneously This method has the advantage of being numerically stable 2

4. Simultaneous Implicit Recalling an initial lecture, we covered two common implicit integration approaches for solving Backward Euler ( ) For a linear system we have = x ( ) f x ( ) + = + + x x f x ( ) t ( ) t t t t t 1 + = x A x ( ) ( ) t t t I t t ( ) ( ) + = + + + x x f x f x ( ) ( ) t ( ) t ( ) t t t t Trapezoidal 2 For a linear system we have t 1 + = + x A A x ( ) ( ) t t t I t I 2 We'll just consider trapezoidal, but for nonlinear cases 3

5. Nonlinear Trapezoidal = x ( ) f x We can use Newton's method to solve with the trapezoidal ( ( ( ) ( ) ( t t t t 2 Right now we are just considering the differential equations; we'll introduce the algebraic equations shortly t ) ) ( ) + + + + + = x x f x f x 0 ) ( ) t t We are solving for x(t+ t); x(t) is known The Jacobian matrix is f x t t t 2 f x f x 1 1 1 n ( ) + = J x I ( ) The I comes from differentiating -x(t+ t) f x n 1 1 n 4

6. Nonlinear Trapezoidal using Newton's Method The full solution would be at each time step Set the initial guess for x(t+ t) as x(t), and initialize the iteration counter k = 0 Determine the mismatch at each iteration k as ( ( ) ( ) t t t t + + + h x x ( ) t ) ( ) ( ) + + + ( ) k ( ) k ( ) k x f x f x ( ) t ( ) ( ) t t t 2 Determine the Jacobian matrix Solve ( ) 1 k 1 + + = + + + ( ) ( ) k ( ) k ( ) k x x ( ( J x h x ( ) ( ) ) ( ) t t t t t t t t Iterate until done 5

7. Infinite Bus GENCLS Example Use the previous two bus system with gen 4 again modeled with a classical model with Xd'=0.3, H=3 and D=0 Bus 2 Bus 1 GENCLS Infinite Bus X=0.22 slack 11.59 Deg 1.095 pu 0.00 Deg 1.000 pu In this example Xth = (0.22 + 0.3), with the internal voltage ? 1= 1.281 23.95 giving E'1=1.281 and 1= 23.95 6

8. Infinite Bus GENCLS Implicit Solution Assume a solid three phase fault is applied at the bus 1 generator terminal, reducing PE1 to zero during the fault, and then the fault is self-cleared at time Tclear, resulting in the post-fault system being identical to the pre-fault system During the fault-on time the equations reduce to = d That is, with a solid fault on the terminal of the generator, during the fault PE1 = 0 1 , 1 pu s dt d 1 ( ) , 1 pu = 1 0 dt 2 3 7

9. Infinite Bus GENCLS Implicit Solution The initial conditions are ( ) ( ) ( ) pu 0 0 . 0 418 0 = = x 0 Let t = 0.02 seconds During the fault the Jacobian is . 0 0 1 3 77 0 . 0 02 2 ( ) + = = I s J x ( ) t t 0 1 Set the initial guess for x(0.02) as x(0), and 0 ( ) ( ) 0 = f x . 0 1667 8

10. Infinite Bus GENCLS Implicit Solution Then calculate the initial mismatch ( ) . 0 02 2 ( ) ( ) ( ) + + + ( ) 0 ( ) 0 ( ) 0 h x x x f x f x ( . 0 02 ) ( . 0 02 ) ( ) 0 ( . 0 02 ) ( ) 0 With x(0.02)(0) = x(0) this becomes . . 0 418 0 0 418 0 0 0 0 . 0 02 2 ( ) = + + + = ( ) 0 h x ( . 0 02 ) . . . 0 167 0 167 0 00334 Then 1 . . . 0 418 0 1 3 77 0 0 0 4306 0 00334 = = ( ) 1 x ( . 0 02 ) . . 1 0 00334 9

11. Infinite Bus GENCLS Implicit Solution Repeating for the next iteration ( 0 1667 ) . 1 259 ( ) ( ) 1 = f x . 0 02 . . . . . 0 4306 0 00334 0 418 0 1 259 0 167 0 . 0 02 2 ( ) = + + + ( ) 1 h x ( . 0 02 ) . . 0 167 . . 0 0 0 0 = . 0 4306 0 00334 Hence we have converged with = x ( . 0 02 ) . 10

12. Infinite Bus GENCLS Implicit Solution Iteration continues until t = Tclear, assumed to be 0.1 seconds in this example . ( . ) . 0 0167 0 7321 = x 0 10 At this point, when the fault is self-cleared, the equations change, requiring a re-evaluation of f(x(Tclear)) d dt d 1 1 281 1 dt 6 0 52 = pu s . 6 30 0 1078 ( ) ( ) + = f x . 0 1 . . pu = sin . 11

13. Infinite Bus GENCLS Implicit Solution With the change in f(x) the Jacobian also changes . 0 1 3 77 . 0 02 2 ( ) = = I s ( ) 0 J x ( . 0 12 ) 0 305 0 00305 . . 0 1 Iteration for x(0.12) is as before, except using the new function and the new Jacobian ) ( . ) ( . ) ( . 0 12 0 12 0 01 + x x ( ) . 0 02 2 ( ( ) ( ) + + + ( ) 0 ( ) 0 ( ) 0 h x f x f x ) ( . 0 12 ) ( . 0 10 ) 1 . . . . . . 0 7321 0 0167 1 3 77 0 1257 0 00216 0 848 0 0142 = = ( ) 1 x ( . 0 12 ) 0 00305 . . 1 This also converges quickly, with one or two iterations 12

14. Computational Considerations As presented for a large system most of the computation is associated with updating and factoring the Jacobian. But the Jacobian actually changes little and hence seldom needs to be rebuilt/factored Rather than using x(t) as the initial guess for x(t+ t), prediction can be used when previous values are available ( ( ) ( ) ( ) t t t t + = + x x x ) ( ) 0 x ( ) t t 13

15. Two Bus System Results The below graph shows the generator angle for varying values of t; recall the implicit method is numerically stable 14

16. Adding the Algebraic Constraints Since the classical model can be formulated with all the values on the network reference frame, initially we just need to add the network equations We'll again formulate the network equations using the form ( , ) or = x y YV YV = ( , ) x y 0 As before the complex equations will be expressed using two real equations, with voltages and currents expressed in rectangular coordinates 15

17. Adding the Algebraic Constraints The network equations are as before 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 16

18. Coupling of x and y with the Classical Model In the simultaneous implicit method x and y are determined simultaneously; hence in the Jacobian we need to determine the dependence of the network equations on x, and the state equations on y With the classical model the Norton current depends on x as , Ni i s i d i R jX I I jI E = + = + = + = = + i + E 1 = + = i I G jB i + R jX , , , , s i d i G ( )( ) + + cos sin j jB Ni DNi jE QNi E i i i i i ( ) cos sin E j Di Qi i E B i i I E G Recall with the classical model Ei is constant DNi Di i Qi i I E B E G QNi Di i Qi i 17

19. Coupling of x and y with the Classical Model In the state equations the coupling with y is recognized by noting = + P Di Di E I Qi Qi E I Ei ( )( ( ) ( ) ) + = + + I jI E V j E V G jB Di Qi Di Di Qi Qi i i ( ( ) ) ( ( ) ) = I E V G E V B These are the algebraic equations Di Di Di i Qi Qi i = + I E V B E V G Qi Di Di i Qi Qi i ( ) ( ) ( ) ( ) ( ) ( ) = + + P E E V G E V B E E V B E V G Ei Di Di Di i Qi Qi i Qi Di Di i Qi Qi i ( ) ( ) ( ) = + + 2 Di 2 Qi P E E V G E Qi Qi E V G Di Qi E V E V B Ei Di Di i i Qi Di i Hence we have PEi written in terms of the voltages (y) 18

20. Variables and Mismatch Equations In solving the Newton algorithm the variables now include x and y (recalling that here y is just the vector of the real and imaginary bus voltages The mismatch equations now include the state integration equations ( ( ( ) ( ) ( t t t t 2 ) + = ( ) k h x ( ) t t ( ) t ) ( ) + + + + + + ( ) k ( ) k ( ) k x x f x , ( y f x ( ), ( ) t y ) ) t t t t And the algebraic equations ( ( ) , ( t t + g x ) + ( ) k ( ) k y ) t t 19

21. Jacobian Matrix Since the h(x,y) and g(x,y) are coupled, the Jacobian is ( ( ( ) , ( ) t t t t + + = + + x ) + + ( ) k ( ) k x , ( y ( ) ) J t t t t ) ( ) + + ( ) k ( ) k ( ) k ( ) k h x y h x , ( y y ( ) ) t t t t x ( ) ( ) + + ( ) k ( ) k ( ) k ( ) k g x , ( y g x , ( y y ( ) ) ( ) ) t t t t t t t t With the classical model the coupling is the Norton current at bus i depends on i (i.e., x) and the electrical power (PEi) in the swing equation depends on VDi and VQi (i.e., y) 20

22. Jacobian Matrix Entries The dependence of the Norton current injections on is cos sin cos sin QNi i i i i i i I E B E G I E G E B = + i i = = + I E G E B DNi i i i i i i = sin cos DNi i i i i i I QNi sin cos E B E G i i i i i i i In the Jacobian the sign is flipped because we defined ? ?,? = ? ? ?(?,? 21

23. Jacobian Matrix Entries The dependence of the swing equation on the generator terminal voltage is = . i i pu s ( ) 1 ( ) = P P D , , i pu Mi Ei i i pu 2H i ( ) ( ) ( ) = + + 2 Di 2 Qi P E E V G E E V G E V E V B Ei Di Di i Qi Qi i Di Qi Qi Di i 1 ( ) , i pu = + E G E B Di i Qi i V 2H Di i 1 ( ) , i pu = E G E B Qi i Di i V 2H Qi i 22

24. Two Bus, Two Gen GENCLS Example We'll reconsider the two bus, two generator case from the previous lecture ; fault at Bus 1, cleared after 0.06 seconds Initial conditions and Ybus are as covered in Lecture 16 Bus 2 Bus 1 GENCLS GENCLS X=0.22 slack 11.59 Deg 1.095 pu 0.00 Deg 1.000 pu PowerWorld Case B2_CLS_2Gen 23

25. Two Bus, Two Gen GENCLS Example Initial terminal voltages are . D1 Q1 V jV 1 0726 E 1 281 23 95 1 1709 j0 52 I j0 3 0 9343 j0 2 I j0 2 + = + + = . , . j0 22 V jV 1 0 D2 Q2 = = . . , . . E 0 955 12 08 1 2 + . . = = . . 1 733 j3 903 N1 . . . = = . 1 j4 6714 N2 . 1 0 . . j7 879 j4 545 j4 545 j9 545 . j0 333 = + = Y Y N . . 1 0 . j0 2 24

26. Two Bus, Two Gen Initial Jacobian V V V V 1 1 1 2 2 D1 0 Q1 0 D2 0 0 0 Q2 0 0 0 . 3 77 0 0 0 0 1 . . . 0 0076 0 0 3 90 1 73 0 0 1 0 0029 0 0 0 7 879 0 0065 0 0 7 879 0 4 545 1 . 0 0 0 0 0 0 1 3 77 2 . . . 4 545 0 9 545 0 0 0039 0 0 4 67 1 00 1 0 0008 0 4 545 0 . 9 545 0 0039 2 . . . . I I I I 0 0 0 0 D1 . 0 . Q1 . . 0 . D2 . . 4 545 Q2 25

27. Results Comparison The below graph compares the angle for the generator at bus 1 using t=0.02 between RK2 and the Implicit Trapezoidal; also Implicit with t=0.06 26

28. Four Bus Comparison Bus 1 Bus 2 GENCLS Bus 4 GENCLS X=0.1 slack Bus 3 7.72 Deg 1.0551 pu X=0.1 X=0.2 -2.40 Deg 0.966 pu 2.31 Deg 1.005 pu 0.00 Deg 1.000 pu 100 MW 50 Mvar 27

29. Four Bus Comparison Fault at Bus 3 for 0.12 seconds; self-cleared 28

30. Done with Transient Stability Solutions: On to Load Modeling Load modeling is certainly challenging! For large system models an aggregate load can consist of many thousands of individual devices The load is constantly changing, with key diurnal and temperature variations For example, a higher percentage of lighting load at night, more air conditioner load on hot days Load model behavior can be quite complex during the low voltages that may occur in transient stability Testing aggregate load models for extreme conditions is not feasible we need to wait for disturbances! 29

31. Load Modeling Traditionally load models have been divided into two groups Static: load is a algebraic function of bus voltage and sometimes frequency Dynamic: load is represented with a dynamic model, with induction motor models the most common The simplest load model is a static constant impedance Has been widely used Allowed the Ybus to be reduced, eliminating essentially all non-generator buses Presents no issues as voltage falls to zero No longer commonly used 30

32. Load Modeling References Many papers and reports are available! A classic reference on load modeling is by the IEEE Task Force on Load Representation for Dynamic Performance, "Load Representation for Dynamic Performance Analysis," IEEE Trans. on Power Systems, May 1993, pp. 472-48 "Final Project Report Loading Modeling Transmission Research" from Lawrence Berkeley National Lab, March 2010 NERC 2016, Dynamic Load Modeling ; available at https://www.nerc.com/comm/PC/LoadModelingTaskForceDL/Dynamic%20Load%20Modeling %20Tech%20Ref%202016-11-14%20-%20FINAL.PDF 31

33. ZIP Load Model Another common static load model is the ZIP, in which the load is represented as ( ( ) 2 = + + Load k P BaseLoad k P P V P V P , , , , , z k k i k k p k ) 2 = + + Load k Q BaseLoad k Q Q V Q V Q , , , , , z k k i k k p k Some models allow more general voltage dependence ( ( ) n1 n2 n3 = + + Load k P BaseLoad k P a V a V a V , , , , , 1 k k 2 k k 3 k k ) n4 n5 n6 = + + Load k Q BaseLoad k Q a V a V a V , , , , , 4 k k 5 k k 6 k k The voltage exponent for reactive power is often > 2 32

34. ZIP Model Coefficients An interesting paper on the experimental determination of the ZIP parameters is A. Bokhari, et. al., "Experimental Determination of the ZIP Coefficients for Modern Residential and Commercial Loads, and Industrial Loads," IEEE Trans. Power Delivery, 2014 Presents test results for loads as voltage is varied; also highlights that load behavior changes with newer technologies Below figure (part of fig 4 of paper), compares real and reactive behavior of light ballast 33

35. ZIP Model Coefficients The Z,I,P coefficients sum to zero; note that for some models the absolute values of the parameters are quite large, indicating a difficult fit A portion of Table VII from Bokhari 2014 paper 34

36. Discharge Lighting Models Discharge lighting (such as fluorescent lamps) is a major portion of the load (10-15%) Discharge lighting has been modeled for sufficiently high voltage with a real power as constant current and reactive power with a high voltage dependence Linear reduction for voltage between 0.65 and 0.75 pu Extinguished (i.e., no load) for voltages below ( ( ) May need to change with newer electronic ballasts e.g., reactive power increasing as the voltage drops! = Disch P Base P V arg eLighting k ) . 4 5 = Q Q V arg Disch eLighting Base k 35

37. Static Load Model Frequency Dependence Frequency dependence is sometimes included, to recognize that the load could change with the frequency ( ( , , , , Load k BaseLoad k z k k i k Q Q Q V Q ) ( ) ( ) = + + 1 P + 2 Load k P BaseLoad k P P V P V P f 1 , , , , , , z k k i k k p k f k k ) ( ) ( ) = + + 1 Q + 2 V Q f 1 , , k p k f k k Here fk is the per unit bus frequency, which is calculated as A typical value for T is about 0.02 seconds. Some models just have frequency dependence on the constant power load s f k k 1 sT + Typical values for Pf and Qf are 1 and -1 respectively 36

38. Induction Machines Term induction machine is used to indicate either generator or motor; most uses are as motors Induction machines have two major components A stationary stator, which is supplied with an ac voltage; windings in stator create a rotating magnetic field A rotating rotor, in which an ac current is induced (hence the name) Two basic design types based on rotor design Squirrel-cage: rotor consists of shorted conducting bars laid into magnetic material in a cage structure Wound-rotor: rotor has windings similar to stator, with slip rings used to provide external access to the rotor windings 37

39. Induction Machine Overview Speed of rotating magnetic field (synchronous speed) depends on number of poles 120 f p = where N is the synchronous speed in RPM, f is s N s s s the stator electrical frequency (e.g., 60 or 50Hz) and p is the number of poles Frequency of induced currents in rotor depends on frequency difference between the rotating magnetic field and the rotor p 2 = r s m where and is the stator electrical frequency, is the rotor electrical frequency is mechanical speed, s m r 38

40. Induction Machine Slip Key value is slip, s, defined as N N s N = s act s wherer N is the synchronous speed, and N is the actual speed (in RPM) act s As defined, when operating as a motor an induction machine will have a positive slip, slip is negative when operating as a generator Slip is zero at synchronous speed, a speed at which no rotor current is induced; s=1 at stand still 39

41. Basic Induction Machine Model A basic (single cage) induction machine circuit model is given below Model is derived in an undergraduate machines class ( ) 1 s s R s = + R R r r r Circuit is useful for understanding the static behavior of the machine Effective rotor resistance (Rr/s) models the rotor electrical losses (Rr) and the mechanical power Rr(1-s)/s 40

42. Induction Machine Dynamics Expressing all values in per unit (with the base covered later), the mechanical equation for a machine is ds 1 T T dt 2H ( ) = M E where H is the inertia constant, T is the mechanical torque and T is the electrical torque (to be defined) M E Similar to what was done for a synchronous machine, the induction machine can be modeled as an equivalent voltage behind a stator resistance and transient reactance (later we'll introduce, but not derive, the subtransient model) 41

43. Induction Machine Dynamics Define X X X X = + r + m X X s X r m = + X X s X m where is locked (s=1) and X is the synchronous reactance is the apparent reactance seen when the rotor Also define the open circuit time constant ( ) r m o s r R + X X = T 42

44. Induction Machine Dynamics Electrically the induction machine is modeled similar to the classical generator model, except here we use the "motor convention" in which ID+jIQ is assumed positive into the machine = + = + + = + = All calculations are done on the network reference frame V E R I X I D D s D Q V E R I X I Q Q s Q D dE dt dE dt 1 T ( ) ( ) sE E X X I D s Q D Q o 1 T ( ) ( ) Q sE E X X I s D Q D o 43

45. Induction Machine Dynamics The induction machine electrical torque, TE, and terminal electrical load, PE, are then ( D D Q Q E s P V I V I = + ) E I E I + Recall we are using the motor convention so positive PE represents load = T E D D Q Q Similar to a synchronous machine, once the initial values are determined the differential equations are fairly easy to simulate Key initial value needed is the slip 44

46. Specifying Induction Machine Parameters In transient stability packages induction machine parameters are specified in per unit If unit is modeled as a generator in the power flow (such as CIMTR1 or GENWRI) then use the generator's MVA base (as with synchronous machines) With loads it is more complicated. Sometimes an explicit MVA base is specified. If so, then use this value. But this can be cumbersome since often the same per unit machine values are used for many loads The default is to use the MW value for the load, often scaled by a multiplier (say 1.25) 45

47. Determining the Initial Values To determine the initial values, it is important to recognize that for a fixed terminal voltage there is only one independent value: the slip, s For a fixed slip, the model is just a simple circuit with resistances and reactances The initial slip is chosen to match the power flow real power value. Then to match the reactive power value (for either a load or a generator), the approach is to add a shunt capacitor in parallel with the induction machine We'll first consider torque-speed curves, then return to determining the initial slip 46