Electrical and Computer Engineering Publications Electrical and Computer Engineering

2018

Sensitivity of Transient Stability Critical Clearing

Time
Shikha Sharma
Iowa State University, ssharma@iastate.edu

Sai Pushpak
Iowa State University, pushpak@iastate.edu

Venkatesh Chinde
Iowa State University

Ian Dobson
Iowa State University, dobson@iastate.edu

Follow this and additional works at: https://lib.dr.iastate.edu/ece_pubs

Part of the Power and Energy Commons
The complete bibliographic information for this item can be found at https://lib.dr.iastate.edu/
ece_pubs/193. For information on how to cite this item, please visit http://lib.dr.iastate.edu/
howtocite.html.

This Article is brought to you for free and open access by the Electrical and Computer Engineering at Iowa State University Digital Repository. It has
been accepted for inclusion in Electrical and Computer Engineering Publications by an authorized administrator of Iowa State University Digital
Repository. For more information, please contact digirep@iastate.edu.
Sensitivity of Transient Stability Critical Clearing Time
Abstract
Once the critical clearing time of a fault leading to transient instability has been computed, it is desirable to
quantify its dependence on system parameters. We derive for a general power system model a new and exact
formula for the first order sensitivity of the critical clearing time with respect to any system parameter. The
formula is evaluated by integrating variational equations forward in time along the base case faulton trajectory
and integrating adjoint variational equations backward in time along the post-fault trajectory. The
computation avoids recomputing the critical clearing time for each parameter change and gives insight into
how parameters influence power system transient stability. The computation of the sensitivity of the critical
clearing time with respect to load impedances and generator inertias is illustrated on a 39-bus system.

Keywords
Power system transient stability, Numerical integration, Nonlinear dynamical systems

Disciplines
Electrical and Computer Engineering | Power and Energy

Comments
This is a preprint of an article published as Sharma, Shikha, Sai Pushpak, Venkatesh Chinde, and Ian Dobson.
"Sensitivity of Transient Stability Critical Clearing Time."

Rights
Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for advertising or promotional
purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any
copyrighted component of this work in other works.

This article is available at Iowa State University Digital Repository: https://lib.dr.iastate.edu/ece_pubs/193

 to appear in IEEE Transactions on Power Systems, accepted June 2018

Sensitivity of Transient Stability Critical Clearing Time

Shikha Sharma, Sai Pushpak, Venkatesh Chinde, Ian Dobson

Abstract—Once the critical clearing time of a fault leading sensitivity computations pioneered by Hiskens [18], [19]. The
to transient instability has been computed, it is desirable to power system model assumed for our calculation is general
quantify its dependence on system parameters. We derive for and widely applicable, only requiring that the fault-on and
a general power system model a new and exact formula for the
first order sensitivity of the critical clearing time with respect to post fault systems be described by smooth, index one, semi-
any system parameter. The formula is evaluated by integrating explicit differential algebraic equations.
variational equations forward in time along the base case fault- Several authors approximately reduce multimachine systems
on trajectory and integrating adjoint variational equations back- to a single machine system to facilitate analysis. Ayasun [20]
ward in time along the post-fault trajectory. The computation expresses the critical clearing time of a one machine infinite
avoids recomputing the critical clearing time for each parameter
change and gives insight into how parameters influence power bus power system model in terms of the load power using the
system transient stability. The computation of the sensitivity of equal area criterion and then linearizes to obtain the first order
the critical clearing time with respect to load impedances and sensitivity of the critical clearing time with respect to system
generator inertias is illustrated on a 39-bus system. load. Ayasun points out the importance of the sensitivity of
Index Terms—Power system transient stability, Numerical the critical clearing time for probabilistic transient stability
integration, Nonlinear dynamical systems assessment, and uses the sensitivity to compute the probability
density function of the critical clearing time. Ayasun assumes
that the multi-machine systems has first been reduced to a one
I. I NTRODUCTION
machine infinite bus power system, whereas our paper directly
To maintain transient stability, a power system fault must computes the critical clearing time sensitivity for the general
be cleared quickly enough so that the fault-on transient multi-machine case. Trajectory sensitivities along the fault-on
remains inside the stability boundary. The critical clearing trajectory are also used by Xu [21] in a one machine infinite
time is the maximum such clearing time, and if the critical bus equivalent of a larger power system to devise preventive
clearing time is exceeded, stability is lost by generators controls to limit angle deviations to stabilize the system.
losing synchronism. An exact computation of critical clearing The functional dependence of the critical clearing time
time requires numerical integration of fault-on and post-fault on parameters can also be approximated from numerically
trajectories and identification of the controlling unstable obtained samples. Chiodo and Lauria [22] use the extended
equilibrium point that determines the relevant portion of the equal area criterion on a 6-bus 3-generator power system
stability boundary. Critical clearing time is a well-established to sample the critical clearing times under variations of the
engineering metric of transient stability and its exact 3 loads. The functional dependence of the logarithm of the
computation by nonlinear analysis and numerical integration critical clearing time on the loads is then obtained by linear
[1]–[3] and its approximate computation by energy methods regression. A multivariate Gaussian model for load power then
[4]–[13] has been extensively studied. leads to a lognormal distribution of critical clearing time to
After critical clearing time has been computed by numerical enable a probabilistic evaluation of transient stability.
integration for a base case, it is desirable to evaluate how In previous work, trajectory sensitivities have been used to
changing the base case parameters affects the critical clearing approximately estimate the critical clearing time. Laufenberg
time. The influential parameters drive the input data require- [23] and Nguyen [10] numerically compute trajectory sensi-
ments, give insight into what affects transient stability, and tivities in the post-fault system and associate the maximum
guide the engineering to increase the critical clearing time if size of the sensitivity trajectory to the proximity to the
it is too short. stability boundary. In particular, Nguyen et al. compute the
One way to approach critical clearing time sensitivity is sensitivities of machine angles and speeds to the clearing
by brute force numerical differencing [14]–[16]. That is, time by computing the trajectory sensitivity forward in time
the critical clearing time is recalculated with the parameter along a fault-on trajectory and further forward in time along
changed to evaluate the change in the critical clearing time. the subsequent post-fault trajectory. They note that these
For example, Khan [17] analyzes the effect of a variety of sensitivities become large during the post-fault trajectory as
parameter changes on the critical clearing time of a single the clearing time approaches the critical clearing time, and
machine infinite bus system by direct simulation. In this paper, therefore use the maximum norm of all the sensitivities as
we avoid this time-consuming recalculation by analyzing the an indicator to estimate the critical clearing time. The high
first-order sensitivity of the critical clearing time to parameters sensitivity is caused by the unstable equilibrium point, but
and taking advantage of the efficient power system trajectory the unstable equilibrium point does not need to be explicitly
located. While Nguyen’s calculation also exploits the trajec-
The authors are with the Department of Electrical and Computer Engineer-
ing, Iowa State University, Ames, IA USA. ssharma@iastate.edu, tory sensitivity techniques of [18] along the fault on and post-
dobson@iastate.edu fault trajectories, it differs from this paper in evaluating the
Preprint 2018
c IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective
works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
 2

critical clearing time, not the sensitivity of the critical clearing left eigenvector of the unstable eigenvalue of the linearization
time, and using the trajectory sensitivity forward in time of the at xu . It will be shown that the normal vector N (xc ) can
post-fault differential equations, whereas we apply trajectory be obtained by starting with the normal vector N (xu ) at the
sensitivities to the adjoint differential equations of the post- unstable equilibrium xu and integrating an adjoint differential
fault system integrated backward in time. Nguyen [10] also equation backward in time along the trajectory ψ.
uses trajectory sensitivities while evaluating the sensitivity If a general parameter is changed, all the equilibria and
of the energy function evaluated at the controlling unstable trajectories and xc change as shown in Fig. 2. The linearized
equilibrium point. Nguyen [10] then extrapolates two samples changes in the equilibria xs and xu are easily obtained.
of these sensitivities at different clearing times to estimate the The linearized change in φ can be obtained by integrating a
critical clearing time. variational fault-on differential equation forward in time from
In this paper we exploit trajectory sensitivities in a novel 0 to tc . The linearized change in the stable manifold W s can
way to give an exact and general formula for first-order be obtained by integrating a variational differential equation
sensitivity of the critical clearing time. The exact formula backward in time along ψ. The change in the stable manifold
and its derivation are new, and include integration backward W s changes its intersection xc with the fault-on trajectory,
in time of an adjoint variational equation along the post- causing a change in the final part of the movement along
fault trajectory. After giving a pictorial overview in Section the fault-on trajectory to the new intersection. The linearized
II, Section III describes the power system model. The new change in clearing time caused by the movement along the
sensitivity formulas are derived in Section IV, and their fault-on trajectory to the new intersection is given by the
numerical application is described in Section V. A 39-bus tangent velocity of the fault-on trajectory at xc (not shown in
example of the calculation is presented in Section VI and Fig. 2), which is given by the fault-on differential equations
Section VII concludes the paper. evaluated at xc . It is plausible, and proved in the following
sections, that suitably combining all these linearized changes
II. P ICTORIAL OVERVIEW yields the sensitivity of the critical clearing time.

N (xc) W s(α1)
W s(α0) xc(α1)
Ws c
x (α0)
xc ψ
φ ψ
u
N (x ) φ
x u xu (α0) xu (α1)
xs xs(α0) xs(α1)

Fig. 1. Visualization of the base case in two dimensions. Fault-on trajectory Fig. 2. Change in the equilibrium points, trajectories, and stability boundary
φ is from the stable equilibrium xs to xc on the stability boundary W s . for a parameter change from the base case value α0 to α1 .
Post-fault trajectory ψ in W s is from xc to the unstable equilibrium xu . N
is the normal vector to W s .
We emphasize that Figs. 1 and 2 require additional di-
We first give an overview of the computation in the two- mensions to be imagined for the intended higher-dimensional
dimensional case that is easiest to depict. (Readers familiar applications. In particular, the stable manifold W s is a hyper-
with the geometry of the nonlinear dynamics of transient sta- surface. For example, Fig. 3 shows W s , its normal vector N ,
bility may skip this overview.) The equilibria and trajectories and the post-fault trajectory ψ in three dimensions.
to be computed in the base case are shown in Fig. 1. The stable
equilibrium of the pre-fault system is xs , and the trajectory φ III. P OWER SYSTEM MODEL
is integrated from xs in the fault-on system until it reaches the This section specifies fault-on and post-fault power system
stability boundary W s at xc at the critical clearing time tc . dynamic models and notation.
The stability boundary W s is in the post-fault system and is
the stable manifold of the unstable equilibrium point xu . The A. Fault-on differential equations
critical post-fault trajectory ψ starts from xc and ends at xu . The fault-on power system differential equations are
In this two dimensional case, the relevant part of the stable
manifold W s coincides with the one dimensional post-fault ẋ = f (x, α) (1)
trajectory ψ. The stable manifold W s can be approximated
where x ∈ Rn is the state, f is smooth, and α is any parameter.
near xc by the dashed tangent line shown in Fig. 1. The
The general solution to (1) with initial condition x0 is
normal vector N (xc ) to W s at xc determines the inclination
of the tangent line. The normal vector N (xu ) of W s at xu is a φ(x0 , t, α) (2)
 3

In particular,
0 = S(xc0 , α0 ) (8)

xc where xc0 c
= x (α0 ).
Ws

ψ D. System critical trajectory

The base case critical trajectory is
φ(xs0 , t, α0 ), 0 ≤ t < tc0

fault-on
(9)
ψ(xc0 , t, α0 ), tc0 ≤ t < ∞ post-fault
xu where tc0 = tc (α0 ). The base case critical trajectory starts at
time zero at xs0 = φ(xs0 , 0, α0 ), passes through xc0 at time
Fig. 3. Post-fault trajectory ψ on the stability boundary W s and the normal tc0 and then tends to xu (α0 ) as time tends to infinity. The
vectors along ψ in three dimensions. critical clearing time tc0 is chosen to make the base case critical
trajectory marginally stable and tend to xu (α0 ).
More generally, accounting for the variation of the critical
xs (α) is the pre-fault stable equilibrium and operating point. trajectory with respect to the parameter α, the critical trajec-
The base case parameter is α0 and the pre-fault base case tory is
stable equilibrium is xs0 = xs (α0 ). The nominal fault-on
φ(xs (α), t, α), 0 ≤ t < tc (α) fault-on

trajectory is φ(xs0 , t, α0 ) for t ≥ 0. (10)
ψ(xc (α), t, α), tc (α) ≤ t < ∞ post-fault
B. Post-fault differential equations The critical trajectory starts at time zero at xs (α), passes
The post-fault power system differential equations are through xc (α) at time tc (α) and tends to xu (α) as time tends
to infinity.
ẋ = F (x, α) (3) When α changes, both xs (α) and xu (α) change and this
where F is smooth. The general solution to (3) with initial affects the fault-on and post-fault trajectories respectively. In
condition x0 is addition, the fault-on and post-fault flows φ and ψ directly
ψ(x0 , t, α) (4) depend on α. These changes in the fault-on and post-fault
trajectories cause the clearing time to change. The sensitivity
The controlling unstable equilibrium is xu (α) and the stable formula derived below quantifies these dependencies.
manifold of xu (α) is W s (xu (α)). W s (xu (α)) is part of
the basin boundary of the post-fault stable operating point IV. S ENSITIVITY FORMULA DERIVATION
xspost (α) [1]. In particular, if the fault-on trajectory has not
This section derives the sensitivity formula using variational
reached W s (xu (α)) when the fault clears, then the system
and dynamical systems methods. A subscripted variable indi-
will restabilize at xspost (α). If the fault-on trajectory crosses
cates (partial) differentiation of that variable with respect to
W s (xu (α)) before the fault clears, then the system is tran-
the subscripted variable and | means “evaluated at”.
siently unstable.
Differentiating (7) with respect to α yields
In the theory derivation in this paper, the power system
model dynamics are expressed as differential equations to 0 = Sx (φx xsα + φt tcα + φα ) + Sα (11)
simplify their expression, similarly to [19]. In practice, differ-
ential equations (1) and (3) are routinely obtained from index Rearranging and using φt = f (xc , α) gives
one semi-explicit differential-algebraic equations. The com- −1
tcα = − (Sx f (xc , α)) (Sx (φx xsα + φα ) + Sα ) (12)
putations can be adapted to apply directly to the differential-
algebraic equations [18], as indicated in the appendix. and evaluating at α0 gives the desired sensitivity formula:

−1
tcα |α0 = − Sx |(xc0 ,α0 ) f (xc0 , α0 ) ×
C. Critical clearing time
s
 
The fault starts at time zero at xs (α). The critical clearing Sx |(xc0 ,α0 ) (φx |(xs0 ,tc0 ,α0 ) xα |α0 + φα |(xs0 ,tc0 ,α0 ) ) + Sα |(xc0 ,α0 )
time tc (α) is the first time that the fault-on trajectory intersects
W s (xu (α)). Write xc (α) for the first intersection of the fault-
(13)
on trajectory with W s (xu (α)): The quantities that need to be computed to evaluate the
xc (α) = φ(xs (α), tc (α), α) (5) sensitivity formula (13) are:
s u 1) Sx |(xc0 ,α0 ) = N (xc0 , α0 ) is a normal vector to the stable
Suppose that W (x (α)) has equation
manifold W s (xu ) at xc0 .
S(x, α) = 0 (6) 2) Sα |(xc0 ,α0 ) is the sensitivity of the stable manifold with
respect to α at xc0 .
near xc (α). A suitable S is defined in section IV-B. Then
3) φx |(xs0 ,tc0 ,α0 ) is the sensitivity of the fault-on trajectory
0 = S(xc (α), α) = S(φ(xs (α), tc (α), α), α) (7) with respect to the initial condition xs (α) at xc0 .
 4

4) φα |(xs0 ,tc0 ,α0 ) is the sensitivity of the fault-on trajectory

with respect to α at xc0 . W s (a ) H F( H )

5) xsα |α0 is the sensitivity of the stable equilibrium with

respect to α and is obtained by solving W s (a 0 )

xT(a)
w(a ) F ( xT (a ))
fx |(xs0 ,α0 ) xsα |α0 = −fα |(xs0 ,α0 ) (14) T
F w(a )
x 0
w0 xu (a )
6) f (xc0 , α0 ) is evaluated directly from the fault-on differ- F ( x0T )
w0 F(xu (a))

ential equation (1). x 0u

F ( x0u )
The following subsections derive the first four quantities in
the list.

Fig. 4. Detail near the controlling unstable equilibrium point xu 0 and its
A. Sensitivity of fault-on trajectory mapping under Φ to locally flatten the stable manifolds for the construction
It follows from differential equation (1) that of the function S. Two stable manifolds are shown by dotted lines, one for
the base case parameter value α0 and the other for the parameter α close to
φt |(xs (α),t,α) = f (φ(xs (α), t, α), α) (15) α0 . The 2 dimensional case is shown.

Differentiating (15) with respect to xs (α) gives

φxt |(xs (α),t,α) = fx |(φ(xs (α),t,α),α) φx |(xs (α),t,α) (16) Now we establish new coordinates near xu (α) with a
transformation Φα in which W s (xu (α)) becomes locally a
Evaluation along the fault-on trajectory gives hyperplane passing through xu (α). Let B ⊂ Rn be a small
φxt |(xs0 ,t,α0 ) = fx |(φ(xs0 ,t,α0 ),α0 ) φx |(xs0 ,t,α0 ) (17) enough ball containing xu0 , and suppose that α is sufficiently
s
close to α0 . Write Wloc (xu (α)) = W s (xu (α)) ∩ B. Let Φα :
Integrating (17) along the fault-on trajectory from time zero to B → Rn be a diffeomorphism for which Φα (Wloc s
(xu (α))) is
tc0 with initial condition the identity matrix yields φx |(xs0 ,tc0 ,α0 ) . a hyperplane E s (α) through xu (α) for each α, and Φα x |xu (α)
Differentiating (15) with respect to α gives is the identity matrix I for each α. This follows from a
φαt |(xs (α),t,α) =fx |(φ(xs (α),t,α),α) (φx |(xs (α),t,α) xsα |α parameterized version of the stable manifold theorem [26]1 .
Then w(α) = w(α)I = w(α)Φα x |xu (α) is normal to both
+ φα |(xs (α),t,α) ) + fα |(φ(xs (α),t,α),α) (18) s
Wloc (xu (α)) and E s (α). In particular, the equation in the
s
Evaluation along the fault-on trajectory gives variable x of Wloc (xu (α)) is
φαt |(xs0 ,t,α0 ) =fx |(φ(xs0 ,t,α0 ),α0 ) (φx |(xs0 ,t,α0 ) xsα |α0 w(α)[Φα (x) − Φα (xu (α))] = 0. (22)
+ φα |(xs0 ,t,α0 ) ) + fα |(φ(xs0 ,t,α0 ),α0 ) (19) Let T be a time such that ψ(xc0 , T, α0 ) is in B. (That
Now φα |(xs0 ,tc0 ,α0 ) is calculated by integrating (19) along the is, integrate along the critical post-fault trajectory until a
fault-on trajectory from time zero to tc0 with initial condition time T when the state is near enough xu0 . We can increase
T later if needed.) We write xT (α) = ψ(xc (α), T, α) and
φα (xs0 , 0, α0 ) = xsα |α0 (20) xT0 = xT (α0 ) = ψ(xc0 , T, α0 ). Let H be a hyperplane through
This fault-on variational trajectory calculation is also done in xT0 transverse to the post-fault trajectory. The local stable
s
Laufenberg [23] for a 17-bus system, and for an equivalent manifold Wloc (xu ) intersects H in a manifold Wloc
s
(xu ) ∩ H
T
single machine system in [21]. of dimension n − 2 near x0 .
In some neighborhood of the post-fault critical trajectory,
we can define τ (x, α) as the time for the trajectory starting at
B. Defining the function S that describes W s
x to first reach H ∩ B. τ generally satisfies
This subsection uses nonlinear dynamical systems methods
to define a function S(x, α) so that the stability boundary and ψ(x, τ (x, α), α) ∈ H ∩ B (23)
stable manifold W s (xu (α)) has equation and, if x is on the base case post-fault trajectory,
S(x, α) = 0 (21) ψ(x, τ (x, α0 ), α0 ) = xT0 .
c
Now we can define
near the critical trajectory, including near x (α). The tools
S(x, α) = w(α) Φα ψ(x, τ (x, α)) − Φα xT (α)



used are standard constructions in nonlinear dynamical sys- (24)
tems, but the derivation is new. Suitable background material
It follows from (22) and (23) that W s (xu (α)) satisfies
for these methods is in [24]–[26].
S(x, α) = 0 near the post-fault critical trajectory. In essence,
We make the generic assumption that the linearized dynam-
S measures the distance of x from the stable manifold
ics at the controlling unstable equilibrium point xu0 has all
W s (xu (α)) by following the trajectory through x until it hits
eigenvalues with negative real parts except for one eigenvalue
the hyperplane H near xu0 and then projecting perpendicular
that is real and positive. It follows that xu (α) is a smooth s
to the local stable manifold Wloc (xu (α)).
function of α sufficiently near α0 , and that a suitably nor-
malized left eigenvector w(α) corresponding to the unstable 1 [26, theorem 6.2] applies to a neighborhood of maps, but this adapts to
eigenvalue is a smooth function of α. We write w0 = w(α0 ). the result for parameterized flows needed here.
 5

Now we discuss the sign of S. There is an ambiguity in and (30) gives

the sign of the left eigenvector w(α) in the formula (24) and 
d

hence an ambiguity in the sign S. Since S(x, α) = 0 describes Sx |ψ(x,t) ψx |(x,t) = −Sx |ψ(x,t) Fx |ψ(x,t) ψx |(x,t) (35)
dt
part of the stable manifold of xu (α), it separates transiently
stable and unstable trajectories. For definiteness we could now Since ψx |(x,t) is invertible,
choose the sign of w(α) so that S(x, α0 ) > 0 for unstable d
trajectories and S(x, α0 ) < 0 for stable trajectories. However, Sx |ψ(x,t) = −Sx |ψ(x,t) Fx |ψ(x,t) (36)
dt
any consistent sign for w(α) and S can be used, since the and evaluating on the base case critical trajectory and restoring
sensitivity formula (13) does not depend on the sign of S. the dependence on α gives
The way that S(x, α) is defined by first moving along the
d
trajectory through x until it meets H also ensures that S is Sx |(ψ(xT0 ,t−T,α0 ),α0 ) =
invariant along trajectories. This can be shown explicitly as dt
− Sx |(ψ(xT0 ,t−T,α0 ),α0 ) Fx |(ψ(xT0 ,t−T,α0 ),α0 ) . (37)
follows:
The initial condition is Sx |(xT0 ,α0 ) = w at t = T and
S(ψ(x, t, α), α)
integrating (37) backward in time from t = T to t = tc0
= w(α) Φα ψ(ψ(x, t, α), τ (ψ(x, t, α), α), α)


yields Sx |(ψ(xT0 ,tc0 −T,α0 ),α0 ) = Sx |(xc0 ,α0 ) Note that (37) is
− Φα xT (α)


(25) the differential equation adjoint to (30) [27].
 α

= w(α) Φ ψ(ψ(x, t, α), τ (x, α) − t, α), α)
D. Computing Sα |(xc0 ,α0 )
− Φα xT (α)


(26)
 α
α T

 The invariance of S (28) on a trajectory through x0 (α) gives
= w(α) Φ ψ(x, τ (x, α), α) − Φ x (α) (27)
S(ψ(x0 (α), t, α), α) = S(x0 (α), α) (38)
= S(x, α) (28)
Differentiating (38) with respect to α,
Equation (25) follows from the definition of S in (24), (26)

follows since the time for an initial point to reach H along its Sx |(ψ(x0 (α),t,α),α) ψx |(x0 (α),t,α) x0α |α + ψα |(x0 (α),t,α) +
trajectory is reduced by t if the initial point is moved for time Sα |(ψ(x0 (α),t,α),α)
t along its trajectory, (27) follows from the basic property of
= Sx |(x0 (α),α) x0α |α + Sα |(x0 (α),α) x0α |α (39)
differential equations that moving along a trajectory for time
t, and then for time τ − t has the same result as moving along Differentiating (39) with respect to t,
a trajectory for time τ , and (28) recalls the definition of S. 
d



Sx |(ψ(x0 (α),t,α),α) ψx |(x0 (α),t,α) x0α |α + ψα |(x0 (α),t,α)
dt
C. Computing the stable manifold normal vector Sx with + Sx |(ψ(x0 (α),t,α),α) ψxt |(x0 (α),t,α) x0α |α + ψαt |(x0 (α),t,α)


adjoint variational equations
d
This subsection computes the stable manifold normal vector + Sα |(ψ(x0 (α),t,α),α) = 0 (40)
dt
Sx |(xc0 ,α0 ) by integrating equations adjoint to the post-fault or, more briefly,
variational equations backward in time.  
It follows from the post-fault differential equation (3) that d d
Sα + Sx (ψx x0α + ψα ) + Sx (ψxt x0α + ψαt ) = 0
dt dt
ψt |(x0 (α),t,α) = F (ψ(x0 (α), t, α), α) (29)
Using (37), (30) and (31),
Differentiating (29) with respect to x0 (α) gives d
Sα −Sx Fx (ψx x0α + ψα ) +
ψxt |(x0 (α),t,α) = Fx |(ψ(x0 (α),t,α),α) ψx |(x0 (α),t,α) (30) dt
Sx (Fx ψx x0α + Fx (ψx x0α + ψα ) + Fα ) = 0 (41)
Differentiating (29) with respect to α gives
so that
ψαt |(x0 (α),t,α) =Fx |(ψ(x0 (α),t,α),α) ψx |(x0 (α),t,α) x0α |α d
Sα + Sx Fx ψx x0α + Sx Fα = 0 (42)


+ ψα |(x0 (α),t,α) + Fα |(ψ(x0 (α),t,α),α) (31) dt
or,
It is convenient to temporarily omit the dependence on α d
from the notation to reduce clutter. The invariance of S along Sα |(ψ(x0 (α),t,α),α) = −Sx |(ψ(x0 (α),t,α),α) Fα |(ψ(x0 (α),t,α),α)
dt
trajectories (28) becomes in the case of the fault-on trajectory − Sx |(ψ(x0 (α),t,α),α) Fx |(ψ(x0 (α),t,α),α) ψx |(x0 (α),t,α) x0α |α
S(ψ(x, t)) = S(x) (32) (43)
Evaluate (43) on the base case post-fault critical trajectory to
Differentiating with respect to x,
get
Sx |ψ(x,t) ψx |(x,t) = Sx |x (33) d
Sα |(ψ(xT0 ,t−T,α0 ),α0 ) =
Differentiating with respect to t gives dt
  − Sx |(ψ(xT0 ,t−T,α0 ),α0 ) Fx |(ψ(xT0 ,t−T,α0 ),α0 ) ψx |(xT0 ,t−T,α0 ) xTα |α0
d
Sx |ψ(x,t) ψx |(x,t) + Sx |ψ(x,t) ψxt |(x,t) = 0 (34) − Sx |(ψ(xT0 ,t−T,α0 ),α0 ) Fα |(ψ(xT0 ,t−T,α0 ),α0 ) (44)
dt
 6

Integrating (44) backward in time from t = T to t = tc0 starting and only outline a simple version of the computations here2 .
from the initial condition Sα |(xT0 ,α0 ) = −w0 xuα derived in the The previous work includes a detailed introduction to finding
next subsection yields Sα |(ψ(xT0 ,tc0 −T,α0 ),α0 ) = Sα |(xc0 ,α0 ) . the controlling unstable equilibrium from both theoretical
and computational viewpoints in [2, chapters 11 and 12],
E. The initial condition Sα |(xT0 ,α0 ) and new continuation [28], optimization [29], and integration
methods [3]. In general terms, to find the unstable equilibrium
The initial condition Sα |(xT0 ,α0 ) is derived and approximated point, one increases the clearing time until one finds the
as follows: Differentiating (24) with respect to α and evaluat- first trajectory diverging from the stable equilibrium, and then
ing at (xT0 , α0 ) gives iterates to find the critical trajectory and clearing time more
Sα |(xT0 ,α0 ) precisely. First, the fault-on critical trajectory φ(xs0 , t, α0 ) is
computed by numerical integration of (15). Then, starting with
= w 0 Φα
 0

x |xT
0
ψt |(xT0 ,0,α0 ) τα |(xT0 ,α) + ψα |(xT0 ,0,α) several points along the fault-on critical trajectory, the post-
α α α
− Φx 0 |xT0 xTα |α0 + Φα0 |ψ(xT0 ,0,α0 ) − Φα0 |xT0


fault critical trajectory is computed by numerical integration
+ wα Φα0 ψ(xT (α0 ), 0, α0 ) − Φα0 xT (α0 )


 of (29) with a shooting method (first bracket the clearing time
 α by finding a transiently stable clearing time and a transiently
= w0 Φx 0 |xT0 F (xT0 , α0 )τα |(xT0 ,α0 ) + ψα |(xT0 ,0,α0 )


unstable clearing time and then shrink the interval containing
α
− Φx 0 |xT0 xTα |α0

(45) the clearing time by an interval-halving algorithm). This yields
the quantities tc0 , xc0 , the post-fault trajectory ψ(xc0 , t, α0 ) for
Since ψ(xT0 , 0, α) = xT0 , ψα |(xT0 ,0,α0 ) = 0. Moreover, as tc0 ≤ t ≤ T and ψ(xc0 , T, α0 ) = xT0 . Then a standard Newton-
α
T → ∞, xT (α0 ) → xu0 , xTα |α0 → xuα |α0 , Φx 0 |xT0 → Raphson algorithm with initial condition xT0 is used to locate
α
Φx 0 |xu0 = identity and w0 F (xT0 , α0 ) → 0. Hence (45) implies the controlling unstable equilibrium xu0 . For the following
that Sα |(xT0 ,α0 ) → −w0 xuα |α0 . Therefore if T is increased as sensitivity computation, the numerical integration of the post-
needed we can use the approximation fault critical trajectory must be accurate enough and T large
enough so that xT0 is close enough to xu0 .
Sα |(xT0 ,α0 ) = −w0 xuα |α0 . (46)
Recalling that S(x, α) measures the distance of x from the
C. Sensitivity computations
stable manifold W s (xu (α)), (46) states that the first order
change in S(x, α) at xT0 due to a change in α is the first order The sensitivity computations are now summarized:
change in xu (α) projected perpendicular to W s (xu0 ). 1) Compute φx |(xs0 ,tc0 ,α0 ) by integrating the fault-on varia-
tional equation (16) from t = 0 to t = tc0 with initial
V. O UTLINE OF COMPUTATIONS condition φx |(xs0 ,0,α0 ) = I.
2) Compute xsα |α0 by solving (14).
This section summarizes the overall computations involved
3) Compute φα |(xs0 ,tc0 ,α0 ) by integrating (19) along the
in evaluating the sensitivity formula (13).
fault-on trajectory from time zero to tc0 with initial
condition φα |(xs0 ,0,α0 ) = xsα |α0 .
A. General requirements 4) Compute the Jacobian Fx |xu0 and the left eigenvector w0
We summarize what is required to apply the sensitiv- corresponding to the unstable eigenvalue. (w0 is normal
ity computation. The sensitivity computation is general and to W s (xu0 ) at xu0 .)
widely applicable. In particular, the sensitivity computation is 5) Use the approximation w = w0 .
applicable if 6) Compute Sx |(xc0 ,α0 ) by integrating the adjoint equation
1) The power system has a smooth, index one, semi-explicit (37) backward in time from t = T to t = tc0 from initial
differential-algebraic model for the fault-on system and condition Sx |(xT0 ,α0 ) = w. This requires the evaluation
for the post-fault system. of the Jacobian Fx |(ψ(xT0 ,t−T,α0 ),α0 ) along the post-fault
2) The critical fault-on and post-fault trajectories, the con- critical trajectory.
trolling unstable equilibrium, and the critical clearing 7) Compute xuα |α0 by solving
time have been determined numerically.
Fx |(xu0 ,α0 ) xuα |α0 = −Fα |(xu0 ,α0 ) (47)
3) The variational methods of Hiskens can be applied
to the critical fault-on and post-fault trajectories. The 8) Evaluate Sα |(xT0 ,α0 ) = −wxuα |α0 (46).
elaboration of usual power system numerical integration 9) Compute ψx |(xT0 ,t−T,α0 ) along the post-fault trajectory
methods to these variational methods is not difficult [18]. by integrating (30) backward from t = T to t = tc0
with initial condition ψx |(xT0 ,0,α0 ) = I. This requires the
B. Preliminary computations
2 Although there are marginal conditions (critical trajectory passing near a
Before performing the sensitivity computations that are the type 2 unstable equilibrium) in which the controlling unstable equilibrium
subject of this paper, it is first necessary to use standard meth- point changes [12], usually the controlling unstable equilibrium is robust.
ods to compute the critical clearing time, critical trajectories, For a robust controlling unstable equilibrium point, the critical trajectory is
sensitive to the exact value of the clearing time, and so should be calculated
and the controlling unstable equilibrium in the base case. For with a robust method such as interval halving to determine the trajectory that
the subtleties of this computation, we refer to previous work, very nearly approaches the controlling unstable equilibrium point.
 7

evaluation of the Jacobian Fx |(ψ(xT0 ,t−T,α0 ),α0 ) along αH is the generator inertia parameter, and it enters (49)
the post-fault trajectory. according to
10) Compute Sα |(xc0 ,α0 ) by integrating (11) backward
in time from t = T to t = tc0 from ini-
tial condition Sα |(xT0 ,α0 ) . This requires the evalua- Hiα = Hi + αiH , i = 30, 31, ..., 39. (54)
tion of Fx |(ψ(xT0 ,t−T,α0 ),α0 ) , Fα |(ψ(xT0 ,t−T,α0 ),α0 ) , and
ψx |(xT0 ,t−T,α0 ) along the post-fault trajectory. Although inertia is a constant for any given generator, here it
11) Compute the sensitivity tcα |α0 using formula (13). can be a parameter since the generator is an equivalent lumped
VI. P OWER SYSTEM EXAMPLE model of group of generators. Indeed, decrease in lumped
inertia is a growing concern as inverter-based generation
This section considers the 39-bus 10-generator power sys- sources displace spinning generators.
tem shown in Fig. 5 to illustrate the computation of the sensi-
tivity of critical clearing time with respect to load impedances The preliminary base case computations are now summa-
and generator inertias. The bus and system parameter values rized. A three phase ground fault is introduced at time zero
are taken from [30]. between bus 17 and bus 18 at 200 km from bus 18 and the
fault-on critical trajectory is computed. The base case critical
clearing time tc0 of 0.34 s, xc , and the corresponding post-fault
trajectory starting at 0.34 s are computed using the shooting
method. The relative rotor angles of the critical trajectory are
shown in Fig. 6. This critical trajectory needs to be computed
to determine the base case critical clearing time. Variations
around this critical trajectory (the fault-on portion forward in
time and the post-fault portion backward in time) are central
to computing the sensitivity of the critical clearing time.

250
Relative rotor angle (degrees)

200

150

100

50

0
Fig. 5. 39 bus example system. A three phase ground fault occurs on the line
between bus 17 and bus 18. -50
0 1 2 3 4
The generator at bus i has fifth order dynamics representing Time (s)
the swing, two axis flux, and field voltage dynamics: Fig. 6. Relative rotor angles of the critical trajectory. This critical trajectory
δ˙i = ωi − ωs , (48) is calculated before the sensitivity calculation that linearizes deviations about
this trajectory is done.
ωs
ω̇i = [Pmi − di (ωi − ωs )
2Hiα
− Vi Iqi cos(θi − δi ) + Vi Idi sin(δi − θi )] (49) To compute the sensitivity of the critical clearing time with
0 0
respect to load impedances αZ and generator inertias αH
Tqi Ėqi = −Eqi + (xdi − x0qi )Idi + Ef di (50) as parameters, we use MATLAB to evaluate all the steps in
0 0
Tdi Ėdi = −Edi + (xqi − x0qi )Iqi (51) Section V and hence evaluate formula (13). All the numerical
TAi Ėf0 di = −Ef0 di + KAi (Vref,i − Vi ) (52) integrations use the MATLAB ode15s solver. The time-varying
matrix variational differential equations (16) have the overall
for i = 30, 31, ..., 39. ωs is the synchronous speed, Pmi is the form Ẋ = M (t)X. They are numerically integrated by
constant mechanical power input and Pei = Re(Vi eθi Ii e−φi ). converting Ẋ = M (t)X to a vector differential equation. The
The current injection at each generator bus is computed from columns of the m × m state matrix X are stacked into a
Ii ∠φi = Idi + jIqi , i = 30, 31, ..., 39. The bus phasor voltage vector x of length m2 , and the vector differential equation is
−1
vector V is computed from the network equations V = Ybus I. ẋ = A(t)x, where the A(t) is the block diagonal m2 × m2
αiZ is the load impedance parameter. αiZ enters the equa- matrix diag{M (t), M (t), ..., M (t)}.
tions via the Ybus matrix entry; e.g., for load at bus 15:
We first focus on the sensitivities to two parameters, the
1 1 1 impedance of load 15 and the inertia of generator 33. The
Y15,15 = + + Z
(53)
Z15,14 Z15,16 Zload15 + α15 critical clearing time first-order sensitivities computed with
 8

formula (13) are

∂tc  0.35
= tcαZ |0 = −5.3816 s/pu (55) Estimated CCT using Sensitivity
∂αZ αZ =0
 15
c  Actual CCT
∂t 
= tcαH |0 = 0.1256 s/pu (56)
∂αH αH =0
 33

These critical clearing time sensitivities can be used in a

linearized model relating the clearing time to the parameter 0.3
change ∆α = α − α0 relative to the base case parameter α0 :
tc (α) = tc0 + tcα |0 ∆α (57)
where tc0 is the base case critical clearing time.
To confirm the sensitivity calculation with the linearization 0.25
(57), we also computed the actual critical clearing time tc as 0 10 20 30 40
a function of the parameter by brute-force re-computing tc as
the parameter varies. The actual critical clearing times and the
Fig. 7. Critical clearing time tc as a function of load impedance parameter
linearized critical clearing times computed from the formula αZ estimated using the sensitivity of the critical clearing time (solid line) and
(13) and the linearization (57) are shown for each parameter the actual critical clearing time (dashed curve).
in Figs. 7 and 8. The tangency of the dashed and solid lines
confirms the correctness of the computation of the sensitivity
of the critical clearing time with formula (13). Figs. 7 and 8
also show the mild nonlinearity of the critical clearing time 1
Estimated CCT using Sensitivity
with respect to load 15 impedance and generator 33 inertia. Actual CCT
To show one way in which the sensitivities can be applied, 0.9

compare the dependence on twenty parameters of the critical

0.8
clearing time based on sensitivities with the actual critical
clearing times in Figs. 9,10 and in Figs. 11,12. It is clear 0.7
that the sensitivities can be used to select the parameters that
0.6
affect the critical clearing time the most, and approximately
quantify this dependence. 0.5
The computation of sensitivity of the critical clearing time
is performed on a 2.4 GHz Intel Core i7 processor in the 0.4

MATLAB R2017b environment. The major effort of the cal-

0.3
culations is integrating variational equations (both forward and 0 2 4 6 8 10 12 14
backward). For our 39-bus simulation example, the overall
computation time including the base case computations is 28 s, Fig. 8. Critical clearing time tc as a function of generator 33 inertia parameter
and the new sensitivity computations by themselves are 9 s αH estimated using the sensitivity of the critical clearing time (solid line) and
per parameter. (The base case computations of 19 s include the actual critical clearing time (dashed curve).
finding the base case critical trajectory and its base case critical
clearing time.) The calculation related to the inverse of Sx is
also of main interest when addressing the computational effort
of the sensitivity of critical clearing time. 0.34

VII. C ONCLUSION 0.32

Given an exact calculation of the transient stability critical
clearing time with its associated critical trajectories, we derive 0.3
a new formula for the first order sensitivity of the critical
clearing time with respect to any power system parameter and
0.28
show how to numerically evaluate the formula using trajectory
sensitivities. The formula and its derivation are novel in power
systems analysis. The new formula is exact but its evaluation 0.26
requires numerical methods. The formula is general and
is widely applicable to power system differential-algebraic 0.24
0 10 20 30 40
models for the fault-on and post-fault systems.
The computations include a conventional variational equa-
tion evaluated along the fault-on critical trajectory and a novel Fig. 9. Critical clearing time tc as a function of eleven load impedance
adjoint variational equation evaluated backward in time along parameters αZ estimated using the sensitivity of the critical clearing time.
the post-fault critical trajectory. Both the normal vector to the
 9

stability boundary hypersurface and the first order variation

of the stability boundary with respect to the parameter are
0.34
propagated backward in time with the adjoint variational
equation, and this is a new method in power systems analysis.
0.32 More generally in computational nonlinear dynamics, it is
challenging to compute with higher dimensional stable mani-
folds because of the complexities of tracking hypersurfaces
0.3 in higher dimensions [31]. Our computation avoids such
difficulties by computing the adjoint variational equation along
the one-dimensional post-fault trajectory that lies in the stable
0.28 manifold. More generally, our computation leverages Hiskens’
efficient trajectory sensitivities calculations [18] and applies
nonlinear dynamics to give a new calculation of the first order
0.26
0 10 20 30 40 sensitivity of a classical metric of transient stability.
Computing the first-order sensitivity of the critical clearing
time avoids the tedious brute-force recomputation of the
Fig. 10. Actual critical clearing time tc as a function of eleven load impedance
parameters αZ .
clearing time with parameters varying from the base case while
quantifying how much various parameters affect the critical
clearing time. Insight into which parameters strongly influence
critical clearing time is basic to increasing critical clearing
time when transient stability is a limiting condition. More-
2
over, there continues to be interest in developing approximate
Gen8
1.8 methods for evaluating transient stability [3], [12], [13]. While
1.6
Gen3 we do not address these approximate methods in this paper,
Gen7
we note that our exact sensitivity calculation can be used to
1.4 Gen9 test and validate any sensitivities that could be obtained via
Gen5
1.2 the approximations. The first-order sensitivity of the critical
Gen2 clearing time is also a useful linearization for probabilistic
1 Gen6
approaches to transient stability [20].
Gen1
0.8
Gen4

0.6 A PPENDIX : ADJUSTMENTS NEEDED FOR

0.4
DIFFERENTIAL - ALGEBRAIC MODELS

0.2 The main text writes the power system model as differential
0 2 4 6 8 10 12 14 equations (1) and (3) for simplicity of expression, whereas
the power system model is often differential-algebraic. This
Fig. 11. Critical clearing time tc as a function of nine generator inertia appendix summarizes the necessary adjustments for the fault-
parameters αH estimated using the sensitivity of the critical clearing time. on model [18].
Suppose that the fault-on power system differential-
algebraic equations are in the semi-implicit form

4
ẋ = g(x, y, α) (58)
0 = h(x, y, α) (59)
3.5 Gen8

where y ∈ Rm is the algebraic state and g and h are smooth

3 Gen3
functions. Write z = (x, y) and write the solution of (58), (59)
2.5 Gen7 with initial condition z0 = (x0 , y0 ) as
Gen9

2 Gen5 φDA (z0 , t, α) = (φD (z0 , t, α), φA (z0 , t, α)) (60)

Gen2
1.5
Gen6
where D indicates differential and A indicates algebraic.
1
Gen1 We assume that we are working in an open set in which we
Gen4
can (in principle and not usually explicitly) solve (59) to obtain
0.5 y = k(x, α), and that hy is nonsingular to ensure index one.
(Also, the MATLAB numerical integration routines that we
0
0 2 4 6 8 10 12 14 use require index one.) Then the fault-on differential equations
equivalent to (58) and (59) are given by
Fig. 12. Actual critical clearing time tc as a function of nine generator inertia ẋ = g(x, k(x, α), α) = f (x, α) (61)
parameters αH .
 10

which is identical to (1). The derivation for the fault-on power [15] Siemens Siguard DSA, “Dynamic Security
system then proceeds exactly as in the main text with the Assessment,” Tech. Rep. [Online]. Available:
http://w3.siemens.com/smartgrid/global/en/products-systems-
differential equations (61). However, to implement numerical solutions/software-solutions/planning-data-management-
methods to compute the results, it is much better to work with software/operation/pages/siguard-dsa.aspx
the original differential-algebraic equations (58), (59). [16] Eurostag, “Dynamic Security Assessment,” Tech. Rep. [Online]. Avail-
able: http://www.eurostag.be/en/products/eurostag/functions/dynamic-
Noting that (15) becomes security/dynamic-security-assessment/
[17] A. Z. Khan, “Effects of power system parameters on critical clear-
s s s
φD D A
t (z (α), t, α) = g(φ (z (α), t, α), φ (z (α), t, α), α) (62) ing time: comprehensive analysis,” Electric Power Systems Research,
vol. 49, no. 1, pp. 37–44, 1999.
0 = h(φD (z s (α), t, α), φA (z s (α), t, α), α), (63) [18] I. A. Hiskens and M. A. Pai, “Trajectory sensitivity analysis of hybrid
systems,” IEEE Trans. Circuits and Systems Part I, vol. 47, no. 2, pp.
the variational equations (17) and (19) along the base case 204–220, 2000.
fault-on trajectory become [19] ——, “Power system applications of trajectory sensitivities,” in IEEE
PES Winter Meeting, vol. 2, 2002, pp. 1200–1205.
[20] S. Ayasun, Y. Liang, and C. O. Nwankpa, “A sensitivity approach for
φD D A
zt |(z0s ,t,α0 ) = gx φz |(z0s ,t,α0 ) + gy φz |(z0s ,t,α0 ) (64) computation of the probability density function of critical clearing time
D A
0 = hx φz |(z0s ,t,α0 ) + hy φz |(z0s ,t,α0 ) (65) and probability of stability in power system transient stability analysis,”
Applied Maths. & Computation, vol. 176, no. 2, pp. 563–576, 2006.
[21] Y. Xu, Z. Y. Dong, J. Zhao, Y. Xue, and D. J. Hill, “Trajectory sensitivity
s analysis on the equivalent one-machine-infinite-bus of multi-machine
φD D D
αt |(z0s ,t,α0 ) = gx (φz |(z0s ,t,α0 ) zα |α0 + φα |(z0s ,t,α0 ) ) systems for preventive transient stability control,” IET Generation,
+ gy (φAz |(z0s ,t,α0 ) zαs |α0 + φAα |(z0s ,t,α0 ) ) + gα (66) Transmission and Distribution, vol. 9, no. 3, pp. 276–286, 2015.
[22] E. Chiodo and D. Lauria, “Transient stability evaluation of multimachine
0= hx (φz |(z0s ,t,α0 ) zαs |α0 + φD
D
α |(z0s ,t,α0 ) ) power systems: a probabilistic approach based upon the extended equal
s area criterion,” IEE Proceedings Gener. Transm. Distrib, vol. 141, no. 6,
+ hy (φz |(z0s ,t,α0 ) zα |α0 + φAα |(z0s ,t,α0 ) )
A
+ gα , (67) pp. 545–553, 1994.
[23] M. J. Laufenberg and M. A. Pai, “A new approach to dynamic secu-
where gx , gy , hx , hy , gα , hα in (64), (65), (66), (67) are eval- rity assessment using trajectory sensitivities,” in IEEE Power Industry
uated at (φDA (z0s , t, α0 ), α0 ). Computer Applications Conf., Columbus OH USA 1997, pp. 272–277.
[24] S. H. Strogatz, Nonlinear dynamics and chaos: with applications to
physics, biology, chemistry, and engineering. Westview press, 2014.
[25] J. Guckenheimer and P. J. Holmes, Nonlinear oscillations, dynamical
R EFERENCES systems, and bifurcations of vector fields. Springer, NY 2013.
[26] J. Palis and W. de Melo, Geometric theory of dynamical systems; an
[1] H.-D. Chiang, F. F. Wu, and P. P. Varaiya, “A BCU method for direct introduction. Springer Verlag, New York., 1982.
analysis of power system transient stability,” IEEE Transactions on [27] J. Hale, “Functional differential equations,” Analytic Theory of Differ-
Power Systems, vol. 9, no. 3, pp. 1194–1208, 1994. ential Equations, pp. 9–22, 1971.
[2] H.-D. Chiang, Direct methods for stability analysis of electric power [28] L. Chen, Y. Min, F. Xu, and K.-P. Wang, “A continuation-based method
systems: theoretical foundation, BCU methodologies, and applications. to compute the relevant unstable equilibrium points for power system
John Wiley & Sons, 2011. transient stability analysis,” IEEE Transactions on Power Systems,
[3] N. Yorino, A. Priyadi, H. Kakui, and M. Takeshita, “A new method for vol. 24, no. 1, pp. 165–172, 2009.
obtaining critical clearing time for transient stability,” IEEE Transactions [29] J. Lee, “An optimization-driven framework for the computation of the
on Power Systems, vol. 25, no. 3, pp. 1620–1626, 2010. controlling UEP in transient stability analysis,” IEEE Transactions on
[4] T. Athay, R. Podmore, and S. Virmani, “A practical method for the direct Automatic Control, vol. 49, no. 1, pp. 115–119, 2004.
analysis of transient stability,” IEEE Transactions on Power Apparatus [30] R. D. Zimmerman, C. E. Murillo-Sánchez, and D. Gan, “Matpower:
and Systems, no. 2, pp. 573–584, 1979. A matlab power system simulation package,” Manual, Power Systems
[5] A. Michel, A. Fouad, and V. Vittal, “Power system transient stability Engineering Research Center, Ithaca NY, vol. 1, 1997.
using individual machine energy functions,” IEEE Transactions on [31] B. Krauskopf, H. M. Osinga, E. J. Doedel, M. E. Henderson, J. Guck-
Circuits and Systems, vol. 30, no. 5, pp. 266–276, 1983. enheimer, A. Vladimirsky, M. Dellnitz, and O. Junge, “A survey of
[6] G. A. Maria, C. Tang, and J. Kim, “Hybrid transient stability analysis methods for computing (un)stable manifolds of vector fields,” Intl.
(power systems),” IEEE Trans. Power Systems, vol. 5, no. 2, pp. 384– Journal Bifurcation and Chaos, vol. 15, no. 03, pp. 763–791, 2005.
393, 1990.
[7] A.-A. Fouad and V. Vittal, Power system transient stability analysis
using the transient energy function method. Pearson Education, 1991.
[8] V. Vittal, E-Z. Zhou, C. Hwang, and A.A. Fouad, “Derivation of stability
limits using analytical sensitivity of the transient energy margin,” IEEE
Trans. Power Systems, vol. 4, no. 4, 1989.
[9] M. Pavella and P. G. Murthy, Transient stability of power systems: theory
and practice. John Wiley and Sons, New York, NY, 1994.
[10] T. B. Nguyen, M. A. Pai, and I. A. Hiskens, “Sensitivity approaches
for direct computation of critical parameters in a power system,”
International Journal of Electrical Power & Energy Systems, vol. 24,
no. 5, pp. 337–343, 2002.
[11] M. A. Pai, Energy function analysis for power system stability. Springer
Science & Business Media, 2012.
[12] L. G. W. Roberts, A. R. Champneys, K. R. W. Bell, and M. di Bernardo,
“Analytical approximations of critical clearing time for parametric
analysis of power system transient stability,” IEEE Journal Emerging
& Selected Topics Circuits & Systems, vol. 5, no. 3, pp. 465–476, 2015.
[13] T. L. Vu, S. M. A. Araifi, M. S. E. Moursi, and K. Turitsyn, “Toward
simulation-free estimation of critical clearing time,” IEEE Transactions
on Power Systems, vol. 31, no. 6, pp. 4722–4731, 2016.
[14] Powertech Labs Inc., “Dynamic Security Assessment Software,”
Tech. Rep. [Online]. Available: http://www.powertechlabs.com/software-
modeling/dynamic-security-assessment-software/