Reachability analysis based transient stability design in power systems
Author's personal copy
Electrical Power and Energy Systems 32 (2010) 782–787
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Reachability analysis based transient stability design in power systems Licheng Jin *, Ratnesh Kumar, Nicola Elia Department of Electrical and Computer Engineering, Iowa State University, Ames, IA 50010, United States
a r t i c l e
i n f o
Article history: Received 9 September 2007 Received in revised form 6 December 2009 Accepted 28 January 2010
Keywords: Backward reachable set Level set method Power system Stability region Transient stability design
a b s t r a c t This paper provides a systematic framework to determine switching control strategies to stabilize the system after a fault if the stabilization is possible. A method to compute the stability region of a stable equilibrium point with the purpose of power system stability analysis is proposed and the validity of discrete controls in transient stability design is studied. First, a Hamilton–Jacobi–Isaas (HJI) partial differential equation (PDE) is constructed to describe the set of backward reachable states as a function of time starting from a target set of states. The backward reachable set of a stable equilibrium point is computed by numerically solving the HJI PDE backwardly in time using level set methods. This backward reachable set yields the stability region of the equilibrium point. Based on such reachability analysis, a transient stability design method is presented. The validity of a discrete control is determined by examining the stability region of the power system with the said control on. If a post-fault initial state is in the stability region of the system with a control on, the control is valid. A control strategy is provided based on the validity of controls. Finally, this method is illustrated by applying to a single machine infinite bus system with the compensation of shunt and series capacitors. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Power system transient stability [1] is related to the ability to maintain synchronism when subjected to a severe disturbance, such as a short circuit on a bus. The resulting system response involves large excursions of the generator rotor angles and is governed by the nonlinear power–angle relationships. Transient stability assessment essentially determines whether or not the post-fault operating state can reach an acceptable steady-state operating point. The conventional method to determine transient stability is to integrate the system equations to obtain a time domain solution of the system variables, given system operating points and contingencies. An alternative method is to determine stability directly [2]. The stability of a post-fault power system can be determined by checking the fault-on trajectory at clearing time. If it lies inside the stability region of a desired stable equilibrium point of the post-fault system, the system is stable. The stability region is defined as the set of all points starting from which the trajectories eventually converge to the stable equilibrium point (SEP) of a general nonlinear autonomous system as time approaches infinity [3]. In the last three decades numerous efforts have been undertaken to determine the stability region with the goal of power system transient stability analysis. The studies of [4–6] provided the theoretical foundations for the geometric structure of the stability region. The authors in * Corresponding author. Tel.: +1 916 221 0458. E-mail address: [email protected] (L. Jin). 0142-0615/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2010.01.014
[6] proved that under certain conditions the stability boundary of a SEP is the union of the stable manifolds of the type one unstable equilibrium points and proposed a numerical algorithm to determine the stability region. Energy function method has also been applied in transient stability assessment. Anthony et al. [7] derived the energy function for machines based on a center of inertia frame of reference. A procedure for first swing transient stability assessment have been developed using the energy function of individual machines and groups of machines. Several techniques were proposed to determine the critical energy related to system’s stability region. One is called the closest unstable equilibrium point (UEP) method [8]. This corresponds to the smallest energy among all the UEPs. The controlling unstable equilibrium point method [9,10] was carried out to provide a less conservative estimation of stability region. This corresponds to the energy at the UEP that is closest to the post-fault trajectory. PEBS (Potential Energy Boundary Surface) was also proposed to approximate the stability region. The energy at the crossing point of the post-fault trajectory and the PEBS was called the critical energy. Chiang and Wu [11] provided theoretical analysis of this approach. Recently, some algorithms have been developed to approximate the stable manifold of an UEP. For example, in [12,13] the Taylor expansion is used to get a quadratic approximation. In [14,15], the stable manifolds around an UEP are approximated by the normal form technique. In [16], the authors apply the singular perturbation theory to decompose a particular power system into slow and fast subsystems based on the assumption that a power system
Author's personal copy
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
can be perfectly separated in time-scale. Then the stability region of a SEP is obtained by numerical simulations. Also controls are very commonly used to improve the transient stability and damping of a power system. The control of a power system can be continuous (generator exciter, governor control, etc.) as well as discrete (line impedance control through series/ shunt capacitors switch, tripping of generators/loads, etc.). Thus, power system control is hybrid, making a controlled power system a hybrid system. The application of hybrid systems based modeling and control techniques for power systems is a recent activity [17– 21]. The application of discrete switching control to stability of power networks has been reported in articles such as [22–29]. Liu and Vu [30] provided a method to obtain voltage stability regions. Vu and Liu [31] applied the stability region in studying the voltage collapse in power systems. In this paper, we provide a novel method to compute the stability region of a given power system with the purpose of mitigating transient stability problems. We then apply this method for validation of a discrete control. This means, when the said control is switched on, following the clearance of a fault, we verify whether the system will reach a desired stable equilibrium point. We also suggest a sequence of discrete controls that may be exercised to steer the system to a desired SEP. The paper is organized as following. Some fundamental concepts of the reachable set analysis are introduced in Section 2. In Section 3 an algorithm to determine the stability region of a SEP is presented. In Section 4 an algorithm is developed to determine the validity of an existing control in power system transient stability assessment. Section 5 presents an illustrative example. Section 6 provides discussions and the conclusion.
2. Reachable set and its characterization Reachable sets are a way of capturing the safety and stability properties of dynamic systems. There are two basic types of reachable sets depending on whether an initial or a final condition is specified. The forward reachable set is defined as the set of all states that can be reached along trajectories that start in a specified initial set. On the other hand, the backward reachable set is the set of states starting from where trajectories can reach the specified target set. The forward and backward reachable sets are shown in Fig. 1. In Section 3, we make use of backward reachable sets to compute the stability region of a stable equilibrium point of a nonlinear system. This is then used to validate a candidate discrete control. A reachable set is a subset of state space. One way of describing a subset of states is via an implicit surface function representation. Consider a closed set S # Rn . An implicit surface representation of S would define a function / : Rn ! R such that /ðxÞ 6 0 if x 2 S and /ðxÞ > 0 if x R S. On the boundary of S, it holds that /ðxÞ ¼ 0. For this reason, /ðxÞ is also called surface function. In [32], the evolution of a backward reachable set is formulated by a Hamilton–Jacobi–Isaacs (HJI) PDE. The viscosity solution of this PDE is an implicit surface representation of a backward reachable set. This HJI PDE can be solved with the very accurate numerical methods drawn from the level set literature [33]. In addition, Tomlin et al. [34] de-
Fig. 1. Illustration of forward and backward reachable sets.
783
rives Hamilton–Jacobi equations whose solutions describe the boundaries of reachable sets for a hybrid system. Consider an autonomous system described by an ordinary differential equation:
dx ¼ f ðxÞ; dt
ð1Þ
where x 2 Rn is the state vector and f(x) is the vector field. As shown in Fig. 2, suppose /ðx; tÞ is the surface function representing the backward reachability set at time t starting from the level surface /ðx; 0Þ. /ðx; tÞ ¼ 0 is a moving surface in (n+1) dimensional space. It is the boundary of the set of all states x 2 Rn from where a target set can be reached in time t or less. The boundary of ‘‘target set” is represented by /ðx; 0Þ ¼ 0. Note that if /ðx; 0Þ 6 0 represents an e neighborhood of an stable equilibrium point, then /ðx; 1Þ 6 0 represents the region of stability of the equilibrium point. Consider the surface /ðx; tÞ ¼ 0 in the (n + 1) dimensional space. For every (x, t) on this surface, the / value is zero. Consider a small variation along this surface, i.e., moving (x, t) to a neighboring point ðx þ dx; t þ dtÞ on the surface, the variation in / will be zero:
d/ ¼ /ðx þ dx; t þ dtÞ /ðx; tÞ ¼ 0; @/ @/ @/ d/ ¼ dx1 þ þ dxn þ dt: @x1 @xn @t From this it follows that,
dx /Tx þ /t ¼ 0: dt
ð2Þ
@/ @/ where /x ¼ ½@x ; . . . ; @x . Superscript T means transpose. Substituting n 1 (1) into (2), we can get
/Tx f ðxÞ þ /t ¼ 0:
ð3Þ
More generally suppose the system is given by a set of differential algebraic equations (DAE) as,
dx ¼ f ðx; yÞ; dt
gðx; yÞ ¼ 0;
ð4Þ
where y 2 Rm is a set of auxiliary variables. In this setting, we will represent the boundary of the backward reachability set at time t by the surface /ðx; y; tÞ ¼ 0 in ðn þ m þ 1Þdimensional space. Then as before,
d/ ¼ /ðx þ dx; y þ dy; t þ dtÞ /ðx; y; tÞ ¼ /Tx dx þ /Ty dy þ /t dt ¼ 0:
Fig. 2. Illustration of backward and forward reachable sets.
Author's personal copy
784
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
(iv) Propagate in time the boundary of the backward reachable set of the target set by solving the following HJI PDE:
From above it follows that,
/Tx
dx dy þ /Ty þ /t ¼ 0: dt dt dx dt
ð5Þ dy , dt
/Tx f ðx; tÞ þ /t ¼ 0;
ð10Þ
We have that ¼ f ðx; yÞ. Now to compute we consider the algebraic constraint equation, gðx; yÞ ¼ 0. Then using argument similar to above,
with terminal conditions kx x k e ¼ 0. The viscosity solution /ðx; tÞ to (10) is the backward reachable set of the target set at time t,
dg ¼ gðx þ dx; y þ dyÞ gðx; yÞ ¼ g Tx dx þ g Ty dy:
fx 2 Rn j/ðx; tÞ 6 0g:
From here it follows that,
dx dy g Tx þ g Ty ¼ 0: dt dt
ð6Þ
Multiplying both sides of (6) by ðg y g Ty Þ1 g y yields,
ðg y g Ty Þ1 g y g Tx
dx dy þ ¼ 0; dt dt
which yields:
dy dx ¼ ðg y g Ty Þ1 g y g Tx : dt dt
The above backward reachable set of the e ball around the stable equilibrium point is computed using a software tool from [38]. This toolbox provides a collection of routines which implement the basic level set algorithms in Matlab for any number of dimensions. We modified one of the routines to make it applicable in calculating the stability region of a small power system. As t goes to infinity, the backward reachable set approaches the true stability region. If the stability region is bounded, the level set based numerical computation of the backward reachability set eventually converges to the stability region within a finite computation time. Example 1. For the nonlinear system (reverse of Van Der Pol Oscillator),
So (5) can be written as:
/Tx f þ /t ¼ /Ty ðg y g Ty Þ1 g y g Tx f
ð7Þ
Thus we obtain the desired PDE and this PDE describes the propagation of the backward reachability set boundary as a function of time. Level set methods are a collection of numerical algorithms for computing the dynamics of moving curves and surfaces. Given a target set defined by an implicit surface function /ðx; t 0 Þ, we use level set methods to compute the propagation of the target set backwardly. By this means, the backward reachable set is obtained. 3. Computation of stability region Section 2 introduces the concept of reachable sets and their computation using HJI PDE. In our past work, it has been applied to power system transient and voltage stability analysis [35–37]. Here we apply the reachable set analysis for the determination of the stability region for power system transient stability assessment. Given a post-fault stable equilibrium point (SEP), there exists an open neighborhood of the SEP that is contained in the stability region. We pick a sufficiently small e ball around the SEP as the target set and propagate its surface function representation backward in time to compute the region of stability of the equilibrium point. The following algorithm summarizes the procedure to determine the stability region of post-fault power system.
dx1 ¼ x2 ; dt dx2 ¼ x1 þ ðx21 1Þx2 ; dt (0, 0) is the only equilibrium point that is stable. So in this case we cannot employ the manifold-based method [6] to find the stability region of the equilibrium point. Our backward reachability based method is applicable in all cases, and using it we compute the stability region as shown in Fig. 3. In Fig. 3, the region inside the dashed line is the stability region of the stable equilibrium point. We validate our result by drawing the phase portrait from which we can see the validity of our computation. 4. Application to control validation Section 3 presented an algorithm to compute the stability region of a given power system. This can be used to validate effectiveness of any discrete control with regards to transient stability. After a fault occurs and is cleared, we need to examine the initial post-fault state. If the initial post-fault state is inside
(i) Form the state space equations of the post-fault power sys¼ f ðxÞ. tem, dx dt (ii) Find the stable equilibrium point of this autonomous nonlinear system, by solving f ðxÞ ¼ 0 and let x 2 Rn be a SEP. (iii) Specify a e ball centered at the stable equilibrium point with sufficiently small radius e. Define an implicit surface function at time t = 0 as,
/ðx; 0Þ ¼ kx x k e:
ð8Þ
Then the target set is the zero ‘‘sublevel” set of the function /ðx; 0Þ, i.e, it is given by,
fx 2 Rn j/ðx; 0Þ 6 0g ¼ fx 2 Rn jkx x k 6 eg:
ð11Þ
ð9Þ
Note that, a point x is inside the target set if /ðx; 0Þ is negative, outside the target set if /ðx; 0Þ is positive, and on the boundary of the target set if /ðx; 0Þ ¼ 0.
Fig. 3. Stability region computed by level set method.
Author's personal copy
785
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
the stability region of the planned discrete control, the trajectory is stable. There may exist many discrete controls in the power system to improve the stability of the trajectory. When the post-fault trajectory becomes unstable, appropriate controls need to be switched on to stabilize it. We first identify the control mode whose stability region includes the post-fault transient state, and switch the control to that mode. If this is not the most desired mode (as determined by the given objective-functions/performance-criteria), it is possible to steer the system to a more desired equilibrium if it holds that there exists a sequence of control modes i1 ; . . . ; in such that for each 1 6 k 6 n 1, the equilibrium state of mode ik is included in the stability region of mode ikþ1 ; i1 is the starting mode (whose stability region includes the post-fault state), and in is the desired target mode. The following algorithm states this more precisely.
5. Example In this section an example is presented to show how the stability region is computed and the control is validated. 5.1. A single-machine-infinite-bus model The classical single-machine-infinite-bus model of power system is shown in Fig. 4. The system model is given as follows:
( dd dt
dx dt
¼x M
¼ Pm Pe
sin dDx M
ð12Þ
Here, d is the machine rotor angle and x is the relative angular velocity of the rotor. Suppose the inertial constant T EU M ¼ xJ0 ¼ 0:026 s2 =rad; P m ¼ 1:0 per unit;P M e ¼ x ¼ 1:35 per unit, and the damping coefficient D ¼ 0:12 s=rad. 5.1.1. The computation of stability region of normal system From the system Eq. (12) and the chosen parameter values, the point (0.8324, 0) is identified to be a stable equilibrium point of this system. We define the target set as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðd 0:8324Þ þ x 6 0:1. The stability region computed using our algorithm lies inside the solid line drawn in Fig. 5. From this
U∠0
E∠δ
Fig. 4. A single-machine-infinite-bus model.
ω (rad / s)
(i) Compute the stability region of a stable equilibrium point x1 of the normal system with no controls on, referred as the mode 1. (ii) Check if the post-fault initial state x0 is inside the stability region of mode 1, if yes, stop. The trajectory is stable without any controls on. Otherwise, the trajectory is unstable, go to next step. (iii) Compute the stability region of a stable equilibrium point xi of the system with control i on, referred to as mode i. (iv) Check if the post-fault initial state x0 is inside the stability region of control i. If yes, control i can be switched on, in which case, the system will eventually reach the equilibrium point xi . Otherwise, control i cannot prevent the instability of the trajectory. (v) In case control i can stabilize the trajectory, but the SEP xi is not a desired operating point, the trajectory can be ‘‘steered” to another SEP xj , whose stability region includes the mode i equilibrium state xi . To achieve this, control j is switched on after switching on the control i and at a time when the system trajectory enters the intersection of two stability regions. (vi) Step v may be repeated to steer the trajectory further to yet another SEP if desired.
jX
δ (rad )
Fig. 5. Stability region and phase portrait for D ¼ 0:12 s=rad.
figure we conclude that if the post-fault initial condition of the state variables is inside the stability region, the trajectories converge to the stable operating point. If the initial condition is outside the stability region, the trajectories will be unstable. We validate our result by drawing the corresponding phase portrait using time domain simulation of some sample trajectories from which we can see that our method can precisely compute the stability region. When the damping coefficient D is increased, the stable equilibrium point remains the same as (0.8324, 0). For D ¼ 0:15, we compute the stability region as shown in Fig. 6. The figure clearly shows that when D is increased, the size of the stability region also increases. The observation is validated by time domain simulations. Figs. 7 and 8 give the time domain responses of the rotor angle and velocity for an initial condition ðd0 ; x0 Þ ¼ ð5; 15Þ when D equals 0.12 s/rad and 0.15 s/rad, respectively. Fig. 8 shows that the trajectories eventually settle at the post-fault stable operating point. However, when D is 0.12 s/rad, the system loses stability for this initial condition as shown in Fig. 7. This is not unexpected since a large D implies a larger stability region.
5.1.2. Transient stability design Fig. 9 shows a single-machine-infinite-bus system with provision for shunt and series controls. Define the system with no controls on as mode 1, with series control on as mode 2, with shunt control on as mode 3, and with both series control and shunt control on as mode 4. When series capacitors are added, the system is switched from mode 1 to mode 2. The reactance X of the transmission line becomes X 1 þ X 2 X series . Suppose the existing series control compen-
Author's personal copy
786
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
jX 1
ω (rad / s )
E∠δ
jX 2
U∠0
jX series
jX shunt
Fig. 9. System model with shunt and series controls.
sation is 40%, then the maximum power transferred equals PM e ¼ EU=ðX 1 þ X 2 X series Þ ¼ 2:25 per unit. When shunts capacitors are added, the system is switched from mode 1 to mode 3. Assume X 1 ¼ 0:5 p:u. Then according to the shunt
Y—D
network
transformation,
X ¼ X 1 þ X 2 XX 1 X 2 ¼ 0:5 þ 0:5 shunt
0:5 0:25 ¼ 0:875. So the maximum power transferred equals δ (rad )
Fig. 6. Stability region and phase portrait for D ¼ 0:15 s=rad.
δ ω
PM e ¼ EU=X ¼ 1:543 per unit. When both the series and shunt capacitors are added, the system is switched from mode 1 to mode 4. In this mode X ¼ 0:475 and the maximum power transferred equals P M e ¼ EU=X ¼ 2:842 per unit. The equilibrium points of these modes related with different controls are shown in Table 1. The stability regions of these four modes are shown in Fig. 10. The stability region of mode 1 is inside the dotted line, that of mode 2 is inside the dashed line, that of
Table 1 Four control modes and their certain parameters. Mode
Series capacitor
Shunt capacitor
X
PM e value
Equilibrium point
1 2 3
Off On Off
Off Off On
X1 þ X2 X 1 þ X 2 X series
1.35 2.25 1.543
(0.8342, 0) (0.4603, 0) (0.7084, 0)
4
On
On
2.3478
(0.4400, 0)
1 X2 X 1 þ X 2 XXshunt
X 1 þ ðX 2 X series Þ X series Þ X 1 ðXX2shunt
Fig. 7. Time domain simulation when D ¼ 0:12 s=rad.
ω (rad / s )
δ ω
δ (rad )
Fig. 8. Time domain simulation when D ¼ 0:15 s=rad.
Fig. 10. System model with shunt and series controls.
Author's personal copy
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
mode 3 is inside the dashed–dotted line, and that of mode 4 is inside the solid line. Based on the stability region method, we can validate the effectiveness of different controls. When post-fault state is inside the stability region of mode 1, no control is needed because the state will finally reach the equilibrium point. When post-fault state is out of the stability region of mode 1, we need to switch on some controls to ensure that the post-fault state lies inside the stability region of one of the modes. For example, if the initial post-fault state is outside of the stability region of mode 2, we can judge that even if the series capacitors are switched on, the system will not maintain stability. That means in this case the series control is ineffective. From Fig. 10, we can also see that the equilibrium points of different modes are all inside the stability region of mode 1. This means whenever a control is switched on and the system finally stabilizes at the equilibrium point of that mode, we can switch off the specific controls so that the system ultimately stabilizes at the equilibrium point of mode 1. It follows that if a transient-fault causes the system state to deviate, then as long as this state lies in the union of 4 regions of stability, it is possible to switch the series/shunt capacitors ‘‘on” and eventually ‘‘off” to return the system to the equilibrium of normal configuration. 6. Conclusion A novel method for computing the stability region of a nonlinear system, such as a power system, is presented in this paper. We also apply this method for power system transient stability design. The proposed method has the following advantages: (i) It computes the stability region accurately. For large systems, the computation may be stopped after a certain number of iterations to get a sub-region contained in the stability region in order to save the computation time. (ii) It is easy to implement. We only need to form the mathematic model of the post-fault power system and identify the stable equilibrium point. After that, we can use level set methods to compute the stability region as a backward reachable set. A limitation is that the complexity suffers from what is commonly known as ‘‘curse of dimensionality”. This is because the numerical propagation of implicit surface function utilizes a gridding of the state space and the number of grid points grow exponentially in the number of dimensions. As part of future research we plan to explore faster and/or approximate techniques for reachability computation. Acknowledgements The work was supported in part by the National Science Foundation under Grants NSF-ECCS-0424048, NSF-ECCS-0601570, NSFECCS-0801763, NSF-CCF-0811541 and NSF-ECCS-0926029. References [1] Kundur P, Paserba J, Ajjarapu V, Andersson G, Bose A, et al. Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions. IEEE Trans Power Syst 2004;19(3):1387–401. [2] Ribben Pavella M, Evans FJ. Direct methods for studying dynamics of large scale power system – a survey. Automatica 1985;32:1–21. [3] Khalil HK. Nonlinear system. 3rd ed. NJ: Prentice Hall; 2002. [4] Varaiya P, Wu FF, Chen R. Direct methods for transient stability analysis of power systems: recent results. In: Proceedings of the IEEE, vol. 73; 1985. p. 1703–15.
787
[5] Zaborszky J, Huang G, Zheng B, Leung TC. On the phase portraits of a class of large nonlinear dynamic systems such as the power systems. IEEE Trans Autom Control 1988;33(1):4–15. [6] Chiang HD, Hirch MW, Wu FF. Stability regions of nonlinear autonomous dynamical systems. IEEE Trans Autom Control 1988;33(1):16–27. [7] Anthony NM, Fouad AA, Vittal V. Power system transient stability using individual machine energy functions. IEEE Trans Circ Syst 1983;30:266–76. [8] Chiang HD, Thorp JS. The closest unstable equilibrium point method for power system dynamic security assessment. IEEE Trans Circ Syst 1989;36:1187–99. [9] Chiang H, Wu F, Varaiya P. A BCU method for direct analysis of power system transient stability. IEEE Trans Power Syst 1994;9:1194–208. [10] Chiang H, Chu C. Theoretical foundation of the BCU method for direct stability analysis of network-reduction power system models with small transfer conductances. IEEE Trans Power Syst 1995;42:252–62. [11] Chiang HD, Wu FF. Foundations of the potential energy boundary surface method for power system transient stability analysis. IEEE Trans Circ Syst 1988;35(6):712–28. [12] Venkatasubramanian V, Ji W. Numerical approximation of (n 1)-dimensional stable manifolds in large systems such as power system. Automatica 1997;33(10):1877–83. [13] Cheng D, Ma J. Calculation of stability region. In: Proceedings of the 42-nd IEEE conference on decision and control, vol. 6, maui, HI, USA; 2003. p. 5615–20. [14] Saha S, Fouad AA, Kliemamm WH, Vittal V. Stability boundary approximation of a power system using the real normal form of vector fields. IEEE Trans Power Syst 1997;12(2):797–802. [15] Qi R, Cook D, Kliemann W, Vittal V. Visualization of stable manifolds and multidimensional surfaces in the analysis of power system dynamics. J Nonlinear Sci 2000;10:175–95. [16] Jing Z, Jia Z, Gao Y. Research on the stability region in a power system. IEEE Trans Circ Syst 2003;50:298–304. [17] Hiskens IA. Power system modeling for inverse problems. IEEE Trans Circ Syst 2004;51(3):539–51. [18] Hiskens IA, Pai MA. Trajectory sensitivity analysis of hybrid systems. IEEE Trans Circ Syst 2000;47(2):204–20. [19] Grijalva S, Hiskens I. Hybrid system modeling of frequency control and load shedding. In: 32nd North American power symposium, Waterloo, Canada; 2000. [20] Qin S, Song Y. The theory of hybrid control systems and its application perspective in electric power systems. In: 2001 international conferences on Info-tech and Info-net, Beijing, China; 2001. p. 85–94. [21] Fourlas GK, Kyriakopoulos KJ, Vournas CD. Hybrid systems modeling for power system. IEEE Circ Syst Mag 2004;4(3):16–23. [22] Haque MH. Improvement of first swing stability limit by utilizing full benefit of shunt FACTS devices. IEEE Trans Power Syst 2004;19(4):1894–902. [23] Chang J, Chow JH. Time-optimal series capacitor control for damping interarea modes in interconnected power systems. IEEE Trans Power Syst 1997;12(1):215–21. [24] Bernard S, Trudel G, Scott G. A 735 kv shunt reactors automatic switching system for Hydro-Quebec network. IEEE Trans Power Syst 1996;11(4):2024–30. [25] Kosterev DN, Kolodziej WJ. Bang–bang series capacitor transient stability control. IEEE Trans Power Syst 1995;10(2):915–24. [26] Chen S, Glavitsch H. Stabilizing switching. IEEE Trans Power Syst 1993;8(4):1511–7. [27] Report IC. Description of discrete supplemental control for stability. IEEE Trans Power Apparatus Syst 1978;97(1):149–65. [28] Kimbark EW. Improvement of power system stability by changes in the network. IEEE Trans Power Apparatus Syst 1969;88(5):773–81. [29] Kimbark EW. Improvement of system stability by switched series capacitors. IEEE Trans Power Apparatus Syst 1966;85(2):180–8. [30] Liu C, Vu K. Analysis of tap-changer dynamics and construction of voltage stability region. IEEE Trans Circ Syst 1989;36:271–89. [31] Vu K, Liu C. Shrinking stability regions and voltage collapse in power systems. IEEE Trans Circ Syst 1992;39:271–89. [32] Mitchell IM. Application of level set methods to control and reachability problems in continuous and hybrid systems. Ph.D. thesis, Stanford University, Stanford, 2002. [33] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. New York: Springer Verlag, Inc.; 2003. [34] Tomlin CJ, Lygeros J, Sastry S. A game theoretic approach to controller design for hybrid systems. In: Proceedings of the IEEE, vol. 88; 2000. p. 949–70. [35] Jin L, Liu H, Kumar R, Ajjarapu V, McCalley JD, Elia N, et al. An application of reachable set analysis in power system transient stability assessment. In: Proceedings of the IEEE power engineering society general meeting, vol. 2; 2005. p. 1715–9. [36] Liu H, Jin L, Ajjarapu V, Kumar R, McCalley JD, Elia N, et al. Reachability analysis based minimal load shedding determination. In: Proceedings of IEEE power engineering society general meeting, vol. 2, San Francisco; 2005. p. 1775–81. [37] Jin L, Liu H, Kumar R, McCalley JD, Elia N, Ajjarapu V. Power system transient stability design using reachability based stability-region computation. In: Proceedings of the 37th annual North American power symposium, Ames; 2005. p. 338–43. [38] Mitchell IM. A toolbox of level set methods version 1.0. .
Electrical Power and Energy Systems 32 (2010) 782–787
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Reachability analysis based transient stability design in power systems Licheng Jin *, Ratnesh Kumar, Nicola Elia Department of Electrical and Computer Engineering, Iowa State University, Ames, IA 50010, United States
a r t i c l e
i n f o
Article history: Received 9 September 2007 Received in revised form 6 December 2009 Accepted 28 January 2010
Keywords: Backward reachable set Level set method Power system Stability region Transient stability design
a b s t r a c t This paper provides a systematic framework to determine switching control strategies to stabilize the system after a fault if the stabilization is possible. A method to compute the stability region of a stable equilibrium point with the purpose of power system stability analysis is proposed and the validity of discrete controls in transient stability design is studied. First, a Hamilton–Jacobi–Isaas (HJI) partial differential equation (PDE) is constructed to describe the set of backward reachable states as a function of time starting from a target set of states. The backward reachable set of a stable equilibrium point is computed by numerically solving the HJI PDE backwardly in time using level set methods. This backward reachable set yields the stability region of the equilibrium point. Based on such reachability analysis, a transient stability design method is presented. The validity of a discrete control is determined by examining the stability region of the power system with the said control on. If a post-fault initial state is in the stability region of the system with a control on, the control is valid. A control strategy is provided based on the validity of controls. Finally, this method is illustrated by applying to a single machine infinite bus system with the compensation of shunt and series capacitors. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Power system transient stability [1] is related to the ability to maintain synchronism when subjected to a severe disturbance, such as a short circuit on a bus. The resulting system response involves large excursions of the generator rotor angles and is governed by the nonlinear power–angle relationships. Transient stability assessment essentially determines whether or not the post-fault operating state can reach an acceptable steady-state operating point. The conventional method to determine transient stability is to integrate the system equations to obtain a time domain solution of the system variables, given system operating points and contingencies. An alternative method is to determine stability directly [2]. The stability of a post-fault power system can be determined by checking the fault-on trajectory at clearing time. If it lies inside the stability region of a desired stable equilibrium point of the post-fault system, the system is stable. The stability region is defined as the set of all points starting from which the trajectories eventually converge to the stable equilibrium point (SEP) of a general nonlinear autonomous system as time approaches infinity [3]. In the last three decades numerous efforts have been undertaken to determine the stability region with the goal of power system transient stability analysis. The studies of [4–6] provided the theoretical foundations for the geometric structure of the stability region. The authors in * Corresponding author. Tel.: +1 916 221 0458. E-mail address: [email protected] (L. Jin). 0142-0615/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2010.01.014
[6] proved that under certain conditions the stability boundary of a SEP is the union of the stable manifolds of the type one unstable equilibrium points and proposed a numerical algorithm to determine the stability region. Energy function method has also been applied in transient stability assessment. Anthony et al. [7] derived the energy function for machines based on a center of inertia frame of reference. A procedure for first swing transient stability assessment have been developed using the energy function of individual machines and groups of machines. Several techniques were proposed to determine the critical energy related to system’s stability region. One is called the closest unstable equilibrium point (UEP) method [8]. This corresponds to the smallest energy among all the UEPs. The controlling unstable equilibrium point method [9,10] was carried out to provide a less conservative estimation of stability region. This corresponds to the energy at the UEP that is closest to the post-fault trajectory. PEBS (Potential Energy Boundary Surface) was also proposed to approximate the stability region. The energy at the crossing point of the post-fault trajectory and the PEBS was called the critical energy. Chiang and Wu [11] provided theoretical analysis of this approach. Recently, some algorithms have been developed to approximate the stable manifold of an UEP. For example, in [12,13] the Taylor expansion is used to get a quadratic approximation. In [14,15], the stable manifolds around an UEP are approximated by the normal form technique. In [16], the authors apply the singular perturbation theory to decompose a particular power system into slow and fast subsystems based on the assumption that a power system
Author's personal copy
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
can be perfectly separated in time-scale. Then the stability region of a SEP is obtained by numerical simulations. Also controls are very commonly used to improve the transient stability and damping of a power system. The control of a power system can be continuous (generator exciter, governor control, etc.) as well as discrete (line impedance control through series/ shunt capacitors switch, tripping of generators/loads, etc.). Thus, power system control is hybrid, making a controlled power system a hybrid system. The application of hybrid systems based modeling and control techniques for power systems is a recent activity [17– 21]. The application of discrete switching control to stability of power networks has been reported in articles such as [22–29]. Liu and Vu [30] provided a method to obtain voltage stability regions. Vu and Liu [31] applied the stability region in studying the voltage collapse in power systems. In this paper, we provide a novel method to compute the stability region of a given power system with the purpose of mitigating transient stability problems. We then apply this method for validation of a discrete control. This means, when the said control is switched on, following the clearance of a fault, we verify whether the system will reach a desired stable equilibrium point. We also suggest a sequence of discrete controls that may be exercised to steer the system to a desired SEP. The paper is organized as following. Some fundamental concepts of the reachable set analysis are introduced in Section 2. In Section 3 an algorithm to determine the stability region of a SEP is presented. In Section 4 an algorithm is developed to determine the validity of an existing control in power system transient stability assessment. Section 5 presents an illustrative example. Section 6 provides discussions and the conclusion.
2. Reachable set and its characterization Reachable sets are a way of capturing the safety and stability properties of dynamic systems. There are two basic types of reachable sets depending on whether an initial or a final condition is specified. The forward reachable set is defined as the set of all states that can be reached along trajectories that start in a specified initial set. On the other hand, the backward reachable set is the set of states starting from where trajectories can reach the specified target set. The forward and backward reachable sets are shown in Fig. 1. In Section 3, we make use of backward reachable sets to compute the stability region of a stable equilibrium point of a nonlinear system. This is then used to validate a candidate discrete control. A reachable set is a subset of state space. One way of describing a subset of states is via an implicit surface function representation. Consider a closed set S # Rn . An implicit surface representation of S would define a function / : Rn ! R such that /ðxÞ 6 0 if x 2 S and /ðxÞ > 0 if x R S. On the boundary of S, it holds that /ðxÞ ¼ 0. For this reason, /ðxÞ is also called surface function. In [32], the evolution of a backward reachable set is formulated by a Hamilton–Jacobi–Isaacs (HJI) PDE. The viscosity solution of this PDE is an implicit surface representation of a backward reachable set. This HJI PDE can be solved with the very accurate numerical methods drawn from the level set literature [33]. In addition, Tomlin et al. [34] de-
Fig. 1. Illustration of forward and backward reachable sets.
783
rives Hamilton–Jacobi equations whose solutions describe the boundaries of reachable sets for a hybrid system. Consider an autonomous system described by an ordinary differential equation:
dx ¼ f ðxÞ; dt
ð1Þ
where x 2 Rn is the state vector and f(x) is the vector field. As shown in Fig. 2, suppose /ðx; tÞ is the surface function representing the backward reachability set at time t starting from the level surface /ðx; 0Þ. /ðx; tÞ ¼ 0 is a moving surface in (n+1) dimensional space. It is the boundary of the set of all states x 2 Rn from where a target set can be reached in time t or less. The boundary of ‘‘target set” is represented by /ðx; 0Þ ¼ 0. Note that if /ðx; 0Þ 6 0 represents an e neighborhood of an stable equilibrium point, then /ðx; 1Þ 6 0 represents the region of stability of the equilibrium point. Consider the surface /ðx; tÞ ¼ 0 in the (n + 1) dimensional space. For every (x, t) on this surface, the / value is zero. Consider a small variation along this surface, i.e., moving (x, t) to a neighboring point ðx þ dx; t þ dtÞ on the surface, the variation in / will be zero:
d/ ¼ /ðx þ dx; t þ dtÞ /ðx; tÞ ¼ 0; @/ @/ @/ d/ ¼ dx1 þ þ dxn þ dt: @x1 @xn @t From this it follows that,
dx /Tx þ /t ¼ 0: dt
ð2Þ
@/ @/ where /x ¼ ½@x ; . . . ; @x . Superscript T means transpose. Substituting n 1 (1) into (2), we can get
/Tx f ðxÞ þ /t ¼ 0:
ð3Þ
More generally suppose the system is given by a set of differential algebraic equations (DAE) as,
dx ¼ f ðx; yÞ; dt
gðx; yÞ ¼ 0;
ð4Þ
where y 2 Rm is a set of auxiliary variables. In this setting, we will represent the boundary of the backward reachability set at time t by the surface /ðx; y; tÞ ¼ 0 in ðn þ m þ 1Þdimensional space. Then as before,
d/ ¼ /ðx þ dx; y þ dy; t þ dtÞ /ðx; y; tÞ ¼ /Tx dx þ /Ty dy þ /t dt ¼ 0:
Fig. 2. Illustration of backward and forward reachable sets.
Author's personal copy
784
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
(iv) Propagate in time the boundary of the backward reachable set of the target set by solving the following HJI PDE:
From above it follows that,
/Tx
dx dy þ /Ty þ /t ¼ 0: dt dt dx dt
ð5Þ dy , dt
/Tx f ðx; tÞ þ /t ¼ 0;
ð10Þ
We have that ¼ f ðx; yÞ. Now to compute we consider the algebraic constraint equation, gðx; yÞ ¼ 0. Then using argument similar to above,
with terminal conditions kx x k e ¼ 0. The viscosity solution /ðx; tÞ to (10) is the backward reachable set of the target set at time t,
dg ¼ gðx þ dx; y þ dyÞ gðx; yÞ ¼ g Tx dx þ g Ty dy:
fx 2 Rn j/ðx; tÞ 6 0g:
From here it follows that,
dx dy g Tx þ g Ty ¼ 0: dt dt
ð6Þ
Multiplying both sides of (6) by ðg y g Ty Þ1 g y yields,
ðg y g Ty Þ1 g y g Tx
dx dy þ ¼ 0; dt dt
which yields:
dy dx ¼ ðg y g Ty Þ1 g y g Tx : dt dt
The above backward reachable set of the e ball around the stable equilibrium point is computed using a software tool from [38]. This toolbox provides a collection of routines which implement the basic level set algorithms in Matlab for any number of dimensions. We modified one of the routines to make it applicable in calculating the stability region of a small power system. As t goes to infinity, the backward reachable set approaches the true stability region. If the stability region is bounded, the level set based numerical computation of the backward reachability set eventually converges to the stability region within a finite computation time. Example 1. For the nonlinear system (reverse of Van Der Pol Oscillator),
So (5) can be written as:
/Tx f þ /t ¼ /Ty ðg y g Ty Þ1 g y g Tx f
ð7Þ
Thus we obtain the desired PDE and this PDE describes the propagation of the backward reachability set boundary as a function of time. Level set methods are a collection of numerical algorithms for computing the dynamics of moving curves and surfaces. Given a target set defined by an implicit surface function /ðx; t 0 Þ, we use level set methods to compute the propagation of the target set backwardly. By this means, the backward reachable set is obtained. 3. Computation of stability region Section 2 introduces the concept of reachable sets and their computation using HJI PDE. In our past work, it has been applied to power system transient and voltage stability analysis [35–37]. Here we apply the reachable set analysis for the determination of the stability region for power system transient stability assessment. Given a post-fault stable equilibrium point (SEP), there exists an open neighborhood of the SEP that is contained in the stability region. We pick a sufficiently small e ball around the SEP as the target set and propagate its surface function representation backward in time to compute the region of stability of the equilibrium point. The following algorithm summarizes the procedure to determine the stability region of post-fault power system.
dx1 ¼ x2 ; dt dx2 ¼ x1 þ ðx21 1Þx2 ; dt (0, 0) is the only equilibrium point that is stable. So in this case we cannot employ the manifold-based method [6] to find the stability region of the equilibrium point. Our backward reachability based method is applicable in all cases, and using it we compute the stability region as shown in Fig. 3. In Fig. 3, the region inside the dashed line is the stability region of the stable equilibrium point. We validate our result by drawing the phase portrait from which we can see the validity of our computation. 4. Application to control validation Section 3 presented an algorithm to compute the stability region of a given power system. This can be used to validate effectiveness of any discrete control with regards to transient stability. After a fault occurs and is cleared, we need to examine the initial post-fault state. If the initial post-fault state is inside
(i) Form the state space equations of the post-fault power sys¼ f ðxÞ. tem, dx dt (ii) Find the stable equilibrium point of this autonomous nonlinear system, by solving f ðxÞ ¼ 0 and let x 2 Rn be a SEP. (iii) Specify a e ball centered at the stable equilibrium point with sufficiently small radius e. Define an implicit surface function at time t = 0 as,
/ðx; 0Þ ¼ kx x k e:
ð8Þ
Then the target set is the zero ‘‘sublevel” set of the function /ðx; 0Þ, i.e, it is given by,
fx 2 Rn j/ðx; 0Þ 6 0g ¼ fx 2 Rn jkx x k 6 eg:
ð11Þ
ð9Þ
Note that, a point x is inside the target set if /ðx; 0Þ is negative, outside the target set if /ðx; 0Þ is positive, and on the boundary of the target set if /ðx; 0Þ ¼ 0.
Fig. 3. Stability region computed by level set method.
Author's personal copy
785
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
the stability region of the planned discrete control, the trajectory is stable. There may exist many discrete controls in the power system to improve the stability of the trajectory. When the post-fault trajectory becomes unstable, appropriate controls need to be switched on to stabilize it. We first identify the control mode whose stability region includes the post-fault transient state, and switch the control to that mode. If this is not the most desired mode (as determined by the given objective-functions/performance-criteria), it is possible to steer the system to a more desired equilibrium if it holds that there exists a sequence of control modes i1 ; . . . ; in such that for each 1 6 k 6 n 1, the equilibrium state of mode ik is included in the stability region of mode ikþ1 ; i1 is the starting mode (whose stability region includes the post-fault state), and in is the desired target mode. The following algorithm states this more precisely.
5. Example In this section an example is presented to show how the stability region is computed and the control is validated. 5.1. A single-machine-infinite-bus model The classical single-machine-infinite-bus model of power system is shown in Fig. 4. The system model is given as follows:
( dd dt
dx dt
¼x M
¼ Pm Pe
sin dDx M
ð12Þ
Here, d is the machine rotor angle and x is the relative angular velocity of the rotor. Suppose the inertial constant T EU M ¼ xJ0 ¼ 0:026 s2 =rad; P m ¼ 1:0 per unit;P M e ¼ x ¼ 1:35 per unit, and the damping coefficient D ¼ 0:12 s=rad. 5.1.1. The computation of stability region of normal system From the system Eq. (12) and the chosen parameter values, the point (0.8324, 0) is identified to be a stable equilibrium point of this system. We define the target set as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðd 0:8324Þ þ x 6 0:1. The stability region computed using our algorithm lies inside the solid line drawn in Fig. 5. From this
U∠0
E∠δ
Fig. 4. A single-machine-infinite-bus model.
ω (rad / s)
(i) Compute the stability region of a stable equilibrium point x1 of the normal system with no controls on, referred as the mode 1. (ii) Check if the post-fault initial state x0 is inside the stability region of mode 1, if yes, stop. The trajectory is stable without any controls on. Otherwise, the trajectory is unstable, go to next step. (iii) Compute the stability region of a stable equilibrium point xi of the system with control i on, referred to as mode i. (iv) Check if the post-fault initial state x0 is inside the stability region of control i. If yes, control i can be switched on, in which case, the system will eventually reach the equilibrium point xi . Otherwise, control i cannot prevent the instability of the trajectory. (v) In case control i can stabilize the trajectory, but the SEP xi is not a desired operating point, the trajectory can be ‘‘steered” to another SEP xj , whose stability region includes the mode i equilibrium state xi . To achieve this, control j is switched on after switching on the control i and at a time when the system trajectory enters the intersection of two stability regions. (vi) Step v may be repeated to steer the trajectory further to yet another SEP if desired.
jX
δ (rad )
Fig. 5. Stability region and phase portrait for D ¼ 0:12 s=rad.
figure we conclude that if the post-fault initial condition of the state variables is inside the stability region, the trajectories converge to the stable operating point. If the initial condition is outside the stability region, the trajectories will be unstable. We validate our result by drawing the corresponding phase portrait using time domain simulation of some sample trajectories from which we can see that our method can precisely compute the stability region. When the damping coefficient D is increased, the stable equilibrium point remains the same as (0.8324, 0). For D ¼ 0:15, we compute the stability region as shown in Fig. 6. The figure clearly shows that when D is increased, the size of the stability region also increases. The observation is validated by time domain simulations. Figs. 7 and 8 give the time domain responses of the rotor angle and velocity for an initial condition ðd0 ; x0 Þ ¼ ð5; 15Þ when D equals 0.12 s/rad and 0.15 s/rad, respectively. Fig. 8 shows that the trajectories eventually settle at the post-fault stable operating point. However, when D is 0.12 s/rad, the system loses stability for this initial condition as shown in Fig. 7. This is not unexpected since a large D implies a larger stability region.
5.1.2. Transient stability design Fig. 9 shows a single-machine-infinite-bus system with provision for shunt and series controls. Define the system with no controls on as mode 1, with series control on as mode 2, with shunt control on as mode 3, and with both series control and shunt control on as mode 4. When series capacitors are added, the system is switched from mode 1 to mode 2. The reactance X of the transmission line becomes X 1 þ X 2 X series . Suppose the existing series control compen-
Author's personal copy
786
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
jX 1
ω (rad / s )
E∠δ
jX 2
U∠0
jX series
jX shunt
Fig. 9. System model with shunt and series controls.
sation is 40%, then the maximum power transferred equals PM e ¼ EU=ðX 1 þ X 2 X series Þ ¼ 2:25 per unit. When shunts capacitors are added, the system is switched from mode 1 to mode 3. Assume X 1 ¼ 0:5 p:u. Then according to the shunt
Y—D
network
transformation,
X ¼ X 1 þ X 2 XX 1 X 2 ¼ 0:5 þ 0:5 shunt
0:5 0:25 ¼ 0:875. So the maximum power transferred equals δ (rad )
Fig. 6. Stability region and phase portrait for D ¼ 0:15 s=rad.
δ ω
PM e ¼ EU=X ¼ 1:543 per unit. When both the series and shunt capacitors are added, the system is switched from mode 1 to mode 4. In this mode X ¼ 0:475 and the maximum power transferred equals P M e ¼ EU=X ¼ 2:842 per unit. The equilibrium points of these modes related with different controls are shown in Table 1. The stability regions of these four modes are shown in Fig. 10. The stability region of mode 1 is inside the dotted line, that of mode 2 is inside the dashed line, that of
Table 1 Four control modes and their certain parameters. Mode
Series capacitor
Shunt capacitor
X
PM e value
Equilibrium point
1 2 3
Off On Off
Off Off On
X1 þ X2 X 1 þ X 2 X series
1.35 2.25 1.543
(0.8342, 0) (0.4603, 0) (0.7084, 0)
4
On
On
2.3478
(0.4400, 0)
1 X2 X 1 þ X 2 XXshunt
X 1 þ ðX 2 X series Þ X series Þ X 1 ðXX2shunt
Fig. 7. Time domain simulation when D ¼ 0:12 s=rad.
ω (rad / s )
δ ω
δ (rad )
Fig. 8. Time domain simulation when D ¼ 0:15 s=rad.
Fig. 10. System model with shunt and series controls.
Author's personal copy
L. Jin et al. / Electrical Power and Energy Systems 32 (2010) 782–787
mode 3 is inside the dashed–dotted line, and that of mode 4 is inside the solid line. Based on the stability region method, we can validate the effectiveness of different controls. When post-fault state is inside the stability region of mode 1, no control is needed because the state will finally reach the equilibrium point. When post-fault state is out of the stability region of mode 1, we need to switch on some controls to ensure that the post-fault state lies inside the stability region of one of the modes. For example, if the initial post-fault state is outside of the stability region of mode 2, we can judge that even if the series capacitors are switched on, the system will not maintain stability. That means in this case the series control is ineffective. From Fig. 10, we can also see that the equilibrium points of different modes are all inside the stability region of mode 1. This means whenever a control is switched on and the system finally stabilizes at the equilibrium point of that mode, we can switch off the specific controls so that the system ultimately stabilizes at the equilibrium point of mode 1. It follows that if a transient-fault causes the system state to deviate, then as long as this state lies in the union of 4 regions of stability, it is possible to switch the series/shunt capacitors ‘‘on” and eventually ‘‘off” to return the system to the equilibrium of normal configuration. 6. Conclusion A novel method for computing the stability region of a nonlinear system, such as a power system, is presented in this paper. We also apply this method for power system transient stability design. The proposed method has the following advantages: (i) It computes the stability region accurately. For large systems, the computation may be stopped after a certain number of iterations to get a sub-region contained in the stability region in order to save the computation time. (ii) It is easy to implement. We only need to form the mathematic model of the post-fault power system and identify the stable equilibrium point. After that, we can use level set methods to compute the stability region as a backward reachable set. A limitation is that the complexity suffers from what is commonly known as ‘‘curse of dimensionality”. This is because the numerical propagation of implicit surface function utilizes a gridding of the state space and the number of grid points grow exponentially in the number of dimensions. As part of future research we plan to explore faster and/or approximate techniques for reachability computation. Acknowledgements The work was supported in part by the National Science Foundation under Grants NSF-ECCS-0424048, NSF-ECCS-0601570, NSFECCS-0801763, NSF-CCF-0811541 and NSF-ECCS-0926029. References [1] Kundur P, Paserba J, Ajjarapu V, Andersson G, Bose A, et al. Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions. IEEE Trans Power Syst 2004;19(3):1387–401. [2] Ribben Pavella M, Evans FJ. Direct methods for studying dynamics of large scale power system – a survey. Automatica 1985;32:1–21. [3] Khalil HK. Nonlinear system. 3rd ed. NJ: Prentice Hall; 2002. [4] Varaiya P, Wu FF, Chen R. Direct methods for transient stability analysis of power systems: recent results. In: Proceedings of the IEEE, vol. 73; 1985. p. 1703–15.
787
[5] Zaborszky J, Huang G, Zheng B, Leung TC. On the phase portraits of a class of large nonlinear dynamic systems such as the power systems. IEEE Trans Autom Control 1988;33(1):4–15. [6] Chiang HD, Hirch MW, Wu FF. Stability regions of nonlinear autonomous dynamical systems. IEEE Trans Autom Control 1988;33(1):16–27. [7] Anthony NM, Fouad AA, Vittal V. Power system transient stability using individual machine energy functions. IEEE Trans Circ Syst 1983;30:266–76. [8] Chiang HD, Thorp JS. The closest unstable equilibrium point method for power system dynamic security assessment. IEEE Trans Circ Syst 1989;36:1187–99. [9] Chiang H, Wu F, Varaiya P. A BCU method for direct analysis of power system transient stability. IEEE Trans Power Syst 1994;9:1194–208. [10] Chiang H, Chu C. Theoretical foundation of the BCU method for direct stability analysis of network-reduction power system models with small transfer conductances. IEEE Trans Power Syst 1995;42:252–62. [11] Chiang HD, Wu FF. Foundations of the potential energy boundary surface method for power system transient stability analysis. IEEE Trans Circ Syst 1988;35(6):712–28. [12] Venkatasubramanian V, Ji W. Numerical approximation of (n 1)-dimensional stable manifolds in large systems such as power system. Automatica 1997;33(10):1877–83. [13] Cheng D, Ma J. Calculation of stability region. In: Proceedings of the 42-nd IEEE conference on decision and control, vol. 6, maui, HI, USA; 2003. p. 5615–20. [14] Saha S, Fouad AA, Kliemamm WH, Vittal V. Stability boundary approximation of a power system using the real normal form of vector fields. IEEE Trans Power Syst 1997;12(2):797–802. [15] Qi R, Cook D, Kliemann W, Vittal V. Visualization of stable manifolds and multidimensional surfaces in the analysis of power system dynamics. J Nonlinear Sci 2000;10:175–95. [16] Jing Z, Jia Z, Gao Y. Research on the stability region in a power system. IEEE Trans Circ Syst 2003;50:298–304. [17] Hiskens IA. Power system modeling for inverse problems. IEEE Trans Circ Syst 2004;51(3):539–51. [18] Hiskens IA, Pai MA. Trajectory sensitivity analysis of hybrid systems. IEEE Trans Circ Syst 2000;47(2):204–20. [19] Grijalva S, Hiskens I. Hybrid system modeling of frequency control and load shedding. In: 32nd North American power symposium, Waterloo, Canada; 2000. [20] Qin S, Song Y. The theory of hybrid control systems and its application perspective in electric power systems. In: 2001 international conferences on Info-tech and Info-net, Beijing, China; 2001. p. 85–94. [21] Fourlas GK, Kyriakopoulos KJ, Vournas CD. Hybrid systems modeling for power system. IEEE Circ Syst Mag 2004;4(3):16–23. [22] Haque MH. Improvement of first swing stability limit by utilizing full benefit of shunt FACTS devices. IEEE Trans Power Syst 2004;19(4):1894–902. [23] Chang J, Chow JH. Time-optimal series capacitor control for damping interarea modes in interconnected power systems. IEEE Trans Power Syst 1997;12(1):215–21. [24] Bernard S, Trudel G, Scott G. A 735 kv shunt reactors automatic switching system for Hydro-Quebec network. IEEE Trans Power Syst 1996;11(4):2024–30. [25] Kosterev DN, Kolodziej WJ. Bang–bang series capacitor transient stability control. IEEE Trans Power Syst 1995;10(2):915–24. [26] Chen S, Glavitsch H. Stabilizing switching. IEEE Trans Power Syst 1993;8(4):1511–7. [27] Report IC. Description of discrete supplemental control for stability. IEEE Trans Power Apparatus Syst 1978;97(1):149–65. [28] Kimbark EW. Improvement of power system stability by changes in the network. IEEE Trans Power Apparatus Syst 1969;88(5):773–81. [29] Kimbark EW. Improvement of system stability by switched series capacitors. IEEE Trans Power Apparatus Syst 1966;85(2):180–8. [30] Liu C, Vu K. Analysis of tap-changer dynamics and construction of voltage stability region. IEEE Trans Circ Syst 1989;36:271–89. [31] Vu K, Liu C. Shrinking stability regions and voltage collapse in power systems. IEEE Trans Circ Syst 1992;39:271–89. [32] Mitchell IM. Application of level set methods to control and reachability problems in continuous and hybrid systems. Ph.D. thesis, Stanford University, Stanford, 2002. [33] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. New York: Springer Verlag, Inc.; 2003. [34] Tomlin CJ, Lygeros J, Sastry S. A game theoretic approach to controller design for hybrid systems. In: Proceedings of the IEEE, vol. 88; 2000. p. 949–70. [35] Jin L, Liu H, Kumar R, Ajjarapu V, McCalley JD, Elia N, et al. An application of reachable set analysis in power system transient stability assessment. In: Proceedings of the IEEE power engineering society general meeting, vol. 2; 2005. p. 1715–9. [36] Liu H, Jin L, Ajjarapu V, Kumar R, McCalley JD, Elia N, et al. Reachability analysis based minimal load shedding determination. In: Proceedings of IEEE power engineering society general meeting, vol. 2, San Francisco; 2005. p. 1775–81. [37] Jin L, Liu H, Kumar R, McCalley JD, Elia N, Ajjarapu V. Power system transient stability design using reachability based stability-region computation. In: Proceedings of the 37th annual North American power symposium, Ames; 2005. p. 338–43. [38] Mitchell IM. A toolbox of level set methods version 1.0. .
