1

Controlling Unstable Equilibrium Point Theory

for Stability Assessment of Two-Time Scale
Power System Models
Luis F. C. Alberto, Member, IEEE, and H. D. Chiang, Fellow, IEEE

time tcl, protection actuates and the power system undergoes a

Abstract—This paper presents the foundations of the structural change. Mathematically the vector field
controlling unstable equilibrium point (CUEP) theory for instantaneously jumps from the fault-on vector field (fF,gF) to
stability assessment of two-time-scale power system models. A the post-fault vector field (f,g). The value (xcl,zcl) of the state
conceptual two-time-scale algorithm to calculate the CUEP is
variables along the fault-on trajectory at the clearing time tcl is
provided. Taking into account the time scale features of power
systems in the CUEP theory has several advantages from both a the initial condition for the post-fault system. Assuming the
numerical and analytical point of view. Numerically, the two-time post-fault system has an asymptotically stable equilibrium
scale algorithm to calculate the CUEP is faster and more robust point (SEP) (xs,zs), the stability analysis, in the occurrence of
when compared to the traditional BCU algorithm. Analytically, a large perturbation, can be translated into the following
we gain more insight into power system dynamics and obtain less question: Does the post-fault trajectory, initiating at the post-
conservative estimates of stability region and critical clearing
fault initial condition (xcl,zcl), converge to a stable equilibrium
time.
point (xs,zs) as t→∞? In other words, does the post-fault initial
Index Terms— BCU algorithm, Controlling Unstable condition (xcl,zcl) belong to the stability region A(xs,zs) of the
Equilibrium Point, Power system transient stability, Singularly SEP (xs,zs)?
perturbed systems. Direct methods, in particular the CUEP method, have
made significant progress in the past years. A complete
I. INTRODUCTION characterization of the stability region was developed in [1].

P ower system models have multi-time-scale features. Based

on physical insight, we can decompose power system
variables into fast and slow variables, modeling fault-on and
The foundations of direct methods for power system transient
stability analysis were provided in [2]; including optimal
estimates of the stability region based on the CUEP method.
post-fault systems as singularly perturbed systems: In spite of that, the time-scale features of power systems have
⎧ F not been taken into account in the stability region
⎪ x& = f ( x, z )
( )
ΣεF ⎨
F
for 0 ≤ t < tcl characterization and the CUEP theory. Neglecting time-scale
⎩ε z& = g ( x, z )
⎪
(1) features in the stability analysis of power systems may lead to
⎧ x& = f ( x, z ) (i) conservative estimates of the stability region and critical
( Σε ) ⎨ & for t ≥ tcl clearing time (CCT); (ii) great computational effort to solve
⎩ε z = g ( x, z ) the differential equations that model the power system; (iii)
Parameter ε is a positive small real number, x is a vector of numerical problems in the calculation of CUEP’s.
slow variables while z is a vector of fast variables. Both the In this paper, the theoretical foundation of the CUEP
post-fault and the fault-on systems possess two-time-scale method for the stability analysis of two-time-scale power
features and their slow and fast variables are almost the same. system models is provided. The advantages of taking time-
The dynamical system (ΣεF) represents the fault-on system scale features into account in the CUEP theory are: (i) the
while (Σε) represents the post-fault system. The superscript F stability analysis can be decomposed into the analysis of
is used in this paper to represent elements associated with the simpler systems (the so called slow and fast subsystems) with
fault-on system, while the subscript F is used to indicate that a smaller dimension; (ii) better insight into the location of
certain element is associated with the fast subsystem. CUEP’s is gained leading to algorithms for CUEP calculation
Usually, the power system is operated around an that can be faster and more robust; (iii) a more comprehensive
equilibrium (xo,zo) until a fault occurs at time t=0. At clearing understanding of instability modes is provided; (iv) less
conservative estimates of stability region and CCT.
L. F. C. Alberto is with the School of Electrical Engineering of São Carlos in
the University of São Paulo, Av. Trabalhador São Carlense, 400 – São Carlos – II. TWO-TIME SCALE SYSTEMS
SP – 13566-590 – Brazil (e-mail: luis@sel.eesc.usp.br). Consider the following general form of singularly perturbed
H. D. Chiang is with the School of Electrical and Computer Engineering,
Cornell University, Ithaca, NY 14850 USA (e-mail: chiang@ece.cornell.edu). systems:

©2008 IEEE.
Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 2

⎧ x& = f ( x, z ) between the singularly perturbed system (Σε) and the

( Σε ) ⎨ (2) simplified subsystems (Σo) and (ΠBLS(x*)). In particular, we
⎩ε z& = g ( x, z )
explore a relationship between Aε(xs,zs), Ao(xs,zs) and
where x ∈ R n , z ∈ R m . Functions f and g are of class C1 and ABLS(x*,z*). The topological boundaries of these sets will be
ε is a small nonnegative scalar. Let φε(t,xo,zo) denote the tra- respectively denoted by ∂Aε(xs,zs), ∂Ao(xs,zs) and ∂ABLS(x*,z*).
jectory of (Σε) starting at (xo,zo) and E denote the set of
equilibrium points of (Σε). Suppose (xs,zs) is an asymptotically
stable equilibrium point of (Σε) and let C. The Constraint Manifold
{ n m
Aε ( xs , zs ) := ( x, z ) ∈ R × R : ϕε (t , x, z ) → ( xs , zs ) as t → ∞ } Typically, Γ is a manifold composed by several disjoint
connected components. Let NH ⊂ Γ be the set of
denote the stability region of (xs,zs). Our main interest is to nonhyperbolic points on Γ, that is, the points of Γ where Dzg
study the stability region Aε(xs,zs) when ε→0. In particular, we has at least one eigenvalue on the imaginary axis. As shown in
want to establish a relationship between Aε(xs,zs) and the [4], the set NH is a manifold of dimension n-1 that separates
stability regions of two simplified systems, namely the slow each of the components of Γ in smaller connected
system and the fast system, for sufficiently small ε. components Γi’s such that Γ \ NH = ∪i Γi .
The continuity of Dzg guarantees that the number of
A. The Slow System eigenvalues of Dzg with positive real part is constant in each
The slow system is formally obtained by setting ε=0 in (2): component Γi of Γ. Therefore, it makes sense to define a type
⎧ x& = f ( x, z ) of stability for each component Γi with respect to the family
( Σo ) ⎨ 0 = g ( x, z ) (3) of fast systems (ΠBLS(x*)).
⎩
The algebraic equation 0=g(x,z) constrains the slow
Definition 1: The connected set Γi is a type-k component of Γ
dynamics to the following set if the matrix Dzg, evaluated at every point of Γi, has k
{ }
Γ := ( x, z ) ∈ R n × R m : g ( x, z ) = 0 . Assuming that rank[Dxg eigenvalues lying in the right-half complex plane and m-k in
the left-half plane.
Dzg]=m for every ( x, z ) ∈ Γ , the local form for submersions
[3] assures that Γ is a n-dimensional smooth manifold in Rn+m. If all the eigenvalues of Dzg, calculated at points of Γi, have
Let φo(t,xo,zo) denote the trajectory of (Σo) starting at (xo,zo) negative real part, then we call Γi a stable component of Γ.
and Ao ( xs , zs ) := {( x, z ) ∈ Γ : ϕo (t , x, z ) → ( xs , zs ) as t → ∞} Otherwise, it is called an unstable component. Note that if
be the stability region of (xs,zs) with respect to (Σo). (x*,z*) lies on a type-k component Γi of Γ, then z* is a type-k
hyperbolic equilibrium point of (ΠBLS(x*)).

B. The Fast System

D. Equilibrium Points and the Stability Boundary
In order to explore the time-scale property of (Σε), define the
fast time scale τ=t/ε. In this new time scale, system (Σε) takes Suppose (xs,zs) is an asymptotically stable equilibrium
the form: point of the post-fault system (Σε). If system (Σε) satisfies the
following assumptions:
⎧ dx
⎪⎪ dτ = ε f ( x, z )
( ε) ⎨
Π (4) (A1) All the equilibriums are hyperbolic;
⎪ dz = g ( x, z ) (A2) The stable and unstable manifolds of the equilibrium
⎪⎩ dτ points on the stability boundary ∂Aε(xs,zs) satisfy the
Let Φε(τ,xo,zo) denote the trajectory of (Πε) starting at transversality condition [1];
(xo,zo). It is evident that Φε(τ,xo,zo)= φε(ετ,xo,zo). The fast (A3) There exists an energy function [1] associated with (Σε).
system (Πo) (or boundary-layer-system) is obtained by
formally setting ε=0 in (4). In the fast system, variable x is Then the stability boundary ∂Aε(xs,zs) of the asymptotically
frozen, as a consequence, one can think the fast system (Πo) stable equilibrium point (xs,zs) of (Σε) equals the closure of the
as a family of boundary layer systems: union of the stable manifolds of the type-one equilibrium
dz points on the stability boundary [2]. More precisely, if
( Π BLS ( x) ) ⎧⎨ = g ( x, z ) (5)
(xui , zui ) , i=1,2,…, are the type-one equilibrium points on the
⎩ dτ
where x is frozen and treated as a parameter. stability boundary ∂Aε(xs,zs) of the asymptotically stable
Let Φo(τ,xo,zo) denote the trajectory of (Πo) starting at (xo,zo) equilibrium point (xs,zs), then
and
{ m
ABLS ( x*, z*) := z ∈ R : Φ o (τ , x*, z ) → ( x*, z*) as t → ∞ } be
∂Aε ( xs , zs ) = UW
i
s
( Σε ) ( xui , zui )

the stability region of the asymptotically stable equilibrium Therefore it is important to establish a relationship
point z* of (ΠBLS(x*)). between type-one equilibrium points that lie on the stability
The main aim of this paper is to establish a relationship boundary of (Σε) and the equilibrium points that lie on the

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 3

stability boundaries of (Σo) and (ΠBLS(x*)). The next theorem stability boundary ∂Ao(xs,zs) of (Σo) if and only if (xu,zu) lies
is a generalization of a theorem presented in [5]. Its proof is on the stability boundary ∂Aε(xs,zs) of (Σε).
given in [6].
The previous result says that a type-one unstable equilibrium
point (xu,zu) on Γs lies on the stability boundary ∂Ao(xs,zs) of
the slow system if and only if it lies on the stability boundary
∂Aε(xs,zs) of the singularly perturbed system for sufficiently
small ε. Fig. 1 illustrates Theorem 2. On the right side of fig.
1, one can see the unstable equilibrium point (xu,zu) on the
stability boundary ∂Ao(xs,zs) of the slow system. On the right
hand side of the same figure, we can see, for sufficiently small
ε, that (xu,zu) also belongs to the stability boundary ∂Aε(xs,zs)
of the singularly perturbed system.

Fig. 1 Geometric illustration of Theorem 2. The unstable equilibrium point

(xu,zu) belongs to the stability boundary ∂Ao(xs,zs) of system (Σo) if and only if
(xu,zu) belongs to the stability boundary ∂Aε(xs,zs) of system (Σε), for sufficiently
small ε.

Theorem 1: [6] If a hyperbolic type-j equilibrium point, say

(x*,z*), of (Σo) lies on a type-k component Γi of Γ, then there
exists ε*>0 such that (x*,z*) is a hyperbolic type-(j+k)
equilibrium point of (Σε) for all ε ∈ (0, ε *) .

Consider the following assumption:

(A4) All equilibrium points of (Σo) are hyperbolic.

Assumption (A4) is generically satisfied, that is, it holds

Fig. 2. Assumptions of Theorem 3. The point (xu,z*) belongs to the stability
for almost all dynamical systems in the form of (Σo). This
region of the asymptotically stable equilibrium point (xs,zs) of (Σo) and it is an
assumption and Theorem 1 guarantee that all equilibrium asymptotically stable equilibrium point of the fast system (ΠBLS(xu)). The point
points of (Σε) are hyperbolic for sufficiently small ε. In other (xu,zu) is an unstable type-one equilibrium point of (ΠBLS(xu)) on the stability
words, assumption (A4) guarantees the satisfaction of boundary of (xu,z*) with respect to the fast subsystem (ΠBLS(xu)). Fig. 3 shows
the conclusions of Theorem 3.
assumption (A1) for sufficiently small ε.
A direct consequence of Theorem 1 is that every type-one Theorem 2 does not reveal the complete picture of the
equilibrium point of (Σε) lies either on a stable component Γs situation because, in many practical situations, the unstable
of Γ or on a type-one unstable component Γu of Γ. Moreover, type-one equilibrium points lie on a type-one component Γu of
if the type-one equilibrium point (xu,zu) of (Σε) lies on a type- Γ. The next theorem is a generalization of Theorem 2 that
one component Γu of Γ, then (xu,zu) is necessarily an studies a relationship between the type-one equilibrium points
asymptotically stable equilibrium point of (Σo). This nice of (Σε) on Γu and the equilibrium points of the fast and slow
property will be explored in the next section to propose a two- systems.
time-scale algorithm for CUEP calculation.
Zou et al. [5] studied a relationship between unstable Theorem 3: [6] Consider system (Σε) and the associated
equilibrium points (UEP) of (Σo) on the stability boundary systems (ΠBLS(x*)) and (Σo) satisfying assumption (A4).
lying on a stable component Γs of Γ and those on the stability Suppose that for each ε>0 the system (Σε) has an energy
boundary of (Σε) for sufficiently small ε, proving the function. Suppose (xs,zs) is an asymptotically stable
following result: equilibrium point of (Σo) on the stable component Γs and
(xu,zu) is an asymptotically stable equilibrium point of (Σo) on
Theorem 2: [5] Suppose (xs,zs) and (xu,zu) are respectively the type-one component Γu. Suppose that (xu,z*) lies on
stable and unstable equilibrium points of (Σo) on the stable Ao ( xs , zs ) ⊂ Γ s and (xu,zu) lies on the stability boundary
component Γs. Suppose that for each ε, the associated ∂ABLS(xu,z*) of the fast system (ΠBLS(xu)). Then, there exists
singularly perturbed system (Σε) has an energy function and ε*>0 such that, for all ε ∈ (0, ε *) , (xu,zu) is a type-one
its equilibrium points are isolated. Then, there exists an ε*>0 unstable equilibrium point of (Σε) lying on the stability
such that for all ε ∈ (0, ε *) , the u.e.p. (xu,zu) lies on the boundary ∂Aε(xs,zs) of (Σε).

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 4

manifold W(sΣε ) ( xco , zco ) contains the exit point.

Fig. 2 illustrates the assumptions of Theorem 2 while Fig. 3
illustrates its conclusions. Theorems 2 and 3 are the basis for
the two-time-scale CUEP method. They provide a scheme to Generically, the exit-point (xex,zex) will toile on the stable
check if a type-one unstable equilibrium point belongs to the manifold of a type-one UEP of (Σε). The CUEP method uses
stability boundary of (Σε) through the analysis of the stability the constant energy surface passing through the CUEP to
region of the slow and fast systems. approximate the relevant part of the stability boundary to
which the fault-on trajectory is heading. Actually, Theorem 3
of [8] asserts that the fault-on trajectory ϕεF (t , xo , zo ) , starting
at a point ( xo , zo ) ∈ Aε ( xs , zs ) with V(xo,zo)<V(xco,zco),
crosses the surface with constant energy equal to
Vcr=V(xco,zco) before crossing the stability boundary
∂Aε(xs,zs). In other words, the critical clearing time estimate
obtained via critical energy Vcr will be always on the safe side,
that is, the estimates will be conservative.
Conceptually, the CUEP method for stability assessment is
composed of the following steps:
Step 1: Determine the CUEP (xco,zco) for the fault-on
trajectory ϕεF (t , ( xo , zo )) .
Step 2: Calculate the critical energy Vcr=V(xco,zco).
( )
Step 3: Calculate the energy Vcl = V ϕεF (tcl , xo , zo ) at the
clearing time tcl.
Step 4: If Vcl<Vcr then the post-fault system is stable.
Fig. 3. Conclusions of Theorem 3. The unstable equilibrium point (xu,zu) is a
type-one UEP that belongs to the stability boundary of the asymptotically stable
Otherwise, it may be unstable.
equilibrium point (xs,zs) of (Σε) for sufficiently small ε. The key step in Controlling UEP method is how to
determine the correct CUEP. BCU method was proposed in
III. TWO-TIME SCALE CUEP THEORY [8] to solve this problem. However, the current version of
For the purpose of decomposing the CUEP theory in two- BCU method does not take into account the time-scale
time-scales, consider the fault-on and the post-fault systems features of power systems. Neglecting time-scale features in
respectively modeled by the singularly perturbed systems this step may lead to numerical problems or excessive
(ΣεF) and (Σε) of (1). computational work. Moreover, neglecting time-scale features
Let (xo,zo) be the equilibrium state of the pre-fault system in steps 2 and 3 may result in very conservative estimates of
and (xs,zs) be the SEP of the post-fault system (Σε) in the CCT.
stable component Гs. The assessment of stability due to large For power system models, ε is a fixed number. However
perturbations consists of checking if the fault is cleared before we want to establish a relationship between the CUEP of the
the fault-on trajectory leaves the stability region of the post- full system (Σε) and the CUEPs of the family of simplified
fault system. In other words, the stability assessment is systems (Σo) and (ΠBLS(x*)). In other words, we want to study
reduced to the problem of obtaining a good estimate of the the behavior of the CUEP when ε→0. This requirement leads
stability region of the post-fault SEP (xs,zs). Actually, in to the following definition of uniform CUEP.
power system stability analysis, it is sufficient to obtain a
good estimate of the relevant part of the stability boundary, Definition 3: Let (xcoε,zcoε) be the CUEP of system (Σε)
that is, the part that is in the “direction” of the fault-on with respect to the fault-on trajectory φεF(t,xo,zo). Consider the
trajectory, to obtain a good estimate of the CCT. The CUEP map ε→(xcoε,zcoε). If there exists ε*>0 such that the map is
method can provide this local estimate of the stability region. constant for all ε ∈ (0, ε *) , then (xco,zco)=(xcoε,zcoε) is a
Consider the fault-on trajectory ϕεF (t , xo , zo ) starting from uniform CUEP with respect to the fault-on trajectory
φε (t,xo,zo) for all ε ∈ (0, ε *) .
F
(xo,zo) at t=0. Let (xex,zex) be the exit point, that is, the point
where the fault-on trajectory ϕεF (t , xo , zo ) crosses the
stability boundary ∂Aε(xs,zs). Our aim is to study the relationship between the uniform
CUEP of (Σε) and the CUEPs of the fast and slow subsystems.
Definition 2: The CUEP of system (Σε), for a fixed ε>0, with In order to define the concept of fast CUEP, consider the
family of fast systems (ΠBLS(x)). Suppose (xo,zo) is the pre-
respect to the fault-on trajectory ϕεF (t , xo , zo ) , is the UEP
fault SEP and zs is a post-fault SEP of (ΠBLS(xo)) lying on Гs.
(xco,zco) on the stability boundary ∂Aε(xs,zs) whose stable
Definition 4: The CUEP of the fast system (ΠBLS(xo)) with

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 5

respect to the fast-fault-on trajectory φoF (t , xo , zo ) is the type- Step 3: Check whether the trajectory φo(t,xo,zcoF) along Гu
one UEP (xo,zcoF) on the stability boundary ∂ABLS(xo,zS) whose converges to an asymptotically stable equilibrium point
(xco,zco) or the pseudo-trajectory φoP(t,xo,zcoF) reaches the exit
stable manifold W(sΠ ) ( xo , zcoF ) contains the exit point point (xexS,zexS) of the slow system.
BLS ( xo )

(xo,zex) of the fast system. Step 4: If the equilibrium point (xco,zco) is reached first, then
(xco,zco) is the uniform CUEP of (Σε) for sufficiently small ε.
Since the fault-on system and the post-fault system are Otherwise calculate the CUEP (xcoS,zcoS) of the slow system
structurally different, their stable constraint manifolds will be with respect to the pseudo-trajectory φoP(t,xo,zcoF). In this case
different. This fact poses some difficulties to define a CUEP (xcoS,zcoS) is the uniform CUEP of (Σε) for sufficiently small ε.
for the slow system. The CUEP of the slow system belongs to
the constraint manifold Γs, while the fault-on trajectory of the Scenario 2: The fast fault-on system is stable.
slow-fault-on system belongs, if it exists, to a different
{ }
constraint manifold Γ Fs ⊂ ( x, z ) ∈ R n × R m : 0 = g F ( x, z ) . In Step 1: Calculate the asymptotically stable equilibrium point
(xo,zsF) of the fast fault-on system on ГsF.
order to define the CUEP for the slow system, we follow the Step 2: Calculate the trajectory φoF(t,xo,zsF) of the slow fault-
same ideas of [5] considering a pseudo-fault-on trajectory on on system along ГsF starting at (xo,zsF) and the corresponding
Γs. For each point, say (x,z), that lies inside the attraction set pseudo-trajectory φoPF(t,xo,zsF) on Гs.
Atr(Γs) of Γs, there exists a corresponding target-point (x,z*) Step 3: Check whether the trajectory of the slow system
that lies on Γs. The target point z* is an asymptotically stable reaches a singularity on ГsF or the corresponding pseudo-
equilibrium point of the post-fault fast system Π FBLS ( x) . trajectory reaches the exit point of the slow system.
Step 4: If the singularity (xsing,zsing) is reached in the first
Definition 5: The exit point (xexS,zexS) of the slow system is place, then a new fast phase is initiating. Restart the analysis
the point where the pseudo-slow-fault-on trajectory ϕoP F
(t ) from the first step of Scenario 1 with (xo,zo) replaced by
(xsing,zsing). Otherwise calculate the CUEP (xcoS,zcoS) of the
crosses the stability boundary ∂Ao(xs,zs) of the slow system. slow system with respect to the fault-on pseudo-trajectory
F F
φoP (t,xo,zs ). In this case, (xcoS,zcoS) is the uniform CUEP of
Definition 6: The CUEP of the slow system (Σo) is the (Σε) for sufficiently small ε.
type-one UEP (xcoS,zcoS) lying on the stable component Γs of Γ
whose stable manifold W(sΣo ) ( xcoS , zcoS ) contains the exit- Step 2 of Scenario 1 uses the fact that the type-one CUEP
point of the slow system. of (Σε) is an asymptotically SEP of (Σo) in Гu. As a
consequence, the slow trajectory along Гu naturally converges
Using the fast and slow CUEP definitions, a two-time-scale to the CUEP avoiding the necessity of using the gradient
CUEP method for stability assessment of two-time-scale reduced system of BCU method [8]. From a numerical point
power system models will be proposed. Two scenarios have to of view, convenient step sizes can be chosen for the slow and
be considered to propose a two-time-scale CUEP method. In fast phases to speed up computation. The exit point and
the first scenario, the fault-on fast system is assumed to be CUEP of the fast and slow system are calculated as proposed
unstable, while in the second scenario the fast fault-on system in the BCU method; that is, a good approximation for the exit
is assumed to be stable and possess an asymptotically stable point is the point along the fault-on trajectory where the first
equilibrium point. maximum of the potential energy is reached. Numerical
algorithms and alternatives to compute the pseudo-trajectories
are briefly discussed in [5].
A. Conceptual Two-Time-Scale Algorithm for CUEP
calculation
B. Conceptual CUEP Based Two-Time-Scale Stability
This algorithm is divided in two scenarios and only considers
Assessment Algorithm
the CUEP calculation. The stability assessment will be
discussed in the next subsection. Once the correct CUEP is calculated, an energy function
can be used to obtain a local estimate of the stability boundary
Scenario 1: The fast fault-on system is unstable. ∂Aε(xs,zs) according to the CUEP method. One issue is the
conservativeness of the estimates when the energy function
Step 1: Calculate the exit-point (xo,zex) and the associated fast does not take into account the time-scale features of the
CUEP (xo,zcoF) ∈ Гu with respect to the fast-fault-on trajectory system. The concept of two-time-scale energy functions
F
Φo (τ,xo,zo). proposed in [7] may considerably reduce the conservativeness
Step 2: Using the fast CUEP (xo,zcoF) as an initial condition, of stability region estimates. In this paper, the
calculate the trajectory φo(t,xo,zcoF) of the slow system on the conservativeness of estimates will be reduced by decomposing
type-one component Гu and the corresponding pseudo- the stability assessment into stability assessment of fast and
trajectory φoP(t,xo,zcoF) on Гs. slow system.

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 6

The following algorithm explores the CUEPs of the fast • If Vslowcl< Vslowcr, then the slow system (Σo) is first-swing
and slow system to obtain a local estimate of the stability stable. Otherwise it may be unstable.
boundary of (Σε) via energy functions. For this purpose, we
assume the existence of a slow energy function Vslow for the The first-swing stability of the slow system may not be
slow system (Σo), a family of fast energy functions Vfast for sufficient to guarantee the stability of (Σε). The next step is a
the family of fast systems (ΠBLS(x)), and an energy function refinement that monitors the energy along the post-fault
Vε for (Σε) for every sufficiently small ε. Again this algorithm trajectory to guarantee stability of (Σε).
is divided in two different scenarios:
Step 2: Monitoring the post-fault trajectory
Scenario 1: The fast-fault-on system is unstable • Calculate the critical energy Vcr=Vε(xcoS,zcoS) at the CUEP
of (Σε).
Step 1: Assessing the stability of the fast fault-on system: • Calculate the energy Vcl=Vε(φFε(tcl,xo,zo)) of the post-fault-
• Calculate the critical energy function of the fast system system at the clearing time tcl.
Vfastcr=Vfast(xo,zcoF) at the CUEP of the fast system • If Vcl<Vcr then the system is first-swing stable. Otherwise
( xo , zcoF ) ∈ Γu . continue with the analysis.
• Calculate the fast energy Vfastcl=Vfast(Φo(τcl,xo,zo)) of the fast • Calculate the post-fault trajectory until a time t* when
post-fault system at the clearing time τcl. V(φε(t*,xcl,zcl))< Vcr.
• If Vfastcl< Vfastcr then the fast system is first-swing stable. • If (x*,z*)= φε(t*,xcl,zcl) belongs to the connected
Otherwise it may be unstable. component Sc(Vcr) of the level set
{( x, z) ∈ R n
× R m : V ( x, z ) < Vcr } that contains (x ,z ) then
s s
The first-swing stability of the fast system may not be
the system is first-swing stable. Otherwise, it may be
sufficient to guarantee the first-swing stability of (Σε). In the
unstable.
next step, the energy of the post-fault trajectory is monitored
to refine the stability assessment of the previous step.
IV. EXAMPLE

Step 2: Monitoring the post-fault system trajectory. The following example illustrates concepts that were
• Calculate the critical energy Vcr=Vε(xco,zco) at the CUEP of introduced in the previous sections. Some tests are conducted
(Σε). to evaluate the proposed two-time-scale algorithm for CUEP
• Calculate the energy Vcl=Vε(ΦFε(τcl,xo,zo)) of the post-fault- calculation and stability assessment. Consider the one-
system at the clearing time τcl. machine-infinite-bus system of fig. 4.
• If Vcl<Vcr then the system is first-swing stable. Otherwise
continue with the analysis.
• Calculate the post-fault trajectory until a time τ* when
Vε(Φε(τ*,xcl,zcl))< Vcr.
• If (x*,z*)= Φε(τ*,xcl,zcl) belongs to the connected
component Sc(Vcr) of the level set Fig. 4. One generator is connected to the infinite bus through four parallel

{( x, z) ∈ R n m
× R : Vε ( x, z ) < Vcr } that contains (xs,zs),
transmission lines. The generator is modeled by the one-axis-model with the
following parameters: H=2.77, T’do=6.6, xd=1.15, x’d=x’q=0.3. The terminal
voltage is 1.0 pu at the pre-fault equilibrium point. The infinite bus is used as a
then the system is first-swing stable. Otherwise, it may be reference for phase angle measurement, and its voltage magnitude is V∞=1.0.
unstable.
The reactance of each line will be denoted by xline. The
Step 2 checks whether the clearing state (xcl,zcl) belongs to
differential equations that model this system are:
Sc(Vcr) to guarantee first-swing stability. If it does not, then
the post-fault trajectory is monitored for a short period of time dEq′ 1 ⎡ ( xd + xline ) ( x − x 'd ) ⎤
= ⎢− E 'q + d V∞ cos δ + E fd ⎥
to check if it enters Sc(Vcr). Usually the fast CCT estimation is dτ T 'do ⎣⎢ ( x 'd + xline ) ( x 'd + xline ) ⎥
⎦
sufficient to guarantee first-swing stability of (Σε). dδ
=ω (6)
dτ
Scenario 2: The fast fault-on system is stable V∞ E 'q sin δ
dω 1 ⎡ ⎤
= ⎢ Pm − − Dω ⎥
Step 1: Assessment of the stability of the slow system.
dτ M ⎢⎣ ( x 'd + xline ) ⎥
⎦

• Calculate the critical energy of the slow system Since T’do is a relatively large time constant when
Vslowcr=Vslow(xcoS,zcoS) at the CUEP of the slow system. compared to electromechanical oscillation time constants, the
• Calculate the energy of the post-fault slow system at the variable E’q is a slow variable while δ and ω are fast variables.
clearing time Vslowcl=Vslow(φoP(tcl,xo,zsF)) along the pseudo Fig. 5 illustrates the time-domain behavior of these variables
fault-on trajectory. for contingency scenario #1. The system is written in the fast

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 7

time scale τ given in seconds. We define ε=1/T’do and t=ετ, the fault is self-cleared at τcl=0.277s. Fig. 5 illustrates the
the slow time scale. Rewriting the previous system in the slow time-domain simulation of this scenario. The CUEP belongs
time-scale and setting ε=0, we obtain the constraint set Г: to the unstable component Γu of Γ as illustrated in Fig. 6.
⎧ V∞ E 'q sin δ ⎫
⎪
( )
Γ = ⎨ E 'q , δ , ω ∈ R3 : ω = 0 and Pm −
⎪ x 'd + xline
⎪
= 0⎬
⎪
⎩ ⎭
Set Г contains all the equilibrium points of (1), including
the CUEPs. Suppose (E’qs,δs,0) is the asymptotically SEP of
interest on the stable component Гs of Г. The singularities in Г
occur at the points (E’qssinδs,(2n+1)π/2,0), n=…-2,-1,0,1,2,…
A component of Г and the corresponding singular point is
shown in Fig. 6. The singular point separates the connected
component of Г in two components, a stable component Гs,
that contains the stable equilibrium point (xs,zs) of interest,
and the type-one unstable component Гu.
Three different contingency scenarios will be studied in the
next subsections. The following energy functions will be used
for stability assessment:
Mω2 V∞ Eq' Fig. 6 - Illustration of contingency scenario #1. The fast fault-on system is
V fast ( Eq' , δ , ω ) = − Pmδ − ' cos δ (7)
2 xd + xline unstable. The two-time scale CUEP algorithm uses the fast-fault-on trajectory to
calculate the fast CUEP. The slow trajectory of the slow system starting at the
'2
xd + xline Eq fast CUEP converges on Гu to the CUEP of (Σε). The shaded area is the stability
Vslow ( Eq' , δ ) = − E fd Eq' − region estimate Sc(Vcr) obtained with the traditional CUEP algorithm with the
xd' + xline 2 (8) energy function (9).
xd − xd' The proposed algorithm calculates the correct CUEP.
− V∞ Eq' cos δ − ( xd − xd' ) Pmδ
xd' + xline Since the fast fault-on system is unstable, the algorithm
calculates the fast CUEP (1.03;2.74;0) with respect to the
Mω2 V∞ Eq'
V ( Eq' , δ , ω ) = − Pmδ − ' cos δ + sustained fast fault-on trajectory and then uses the trajectory
2 xd + xline (9) of the slow system along Γu to approach the CUEP
( xd + xline ) Eq' 2 E fd (0.41;1.77;0) of (Σε). Fig. 6 illustrates the calculation steps
+ − Eq'
( xd' + xline )( xd − xd' ) 2 ( xd − xd' ) and the location of CUEP’s. In this case, the step size used to
numerically calculate the slow post-fault system trajectory can
be bigger than the step size used to calculate the fast-fault-on
trajectory. Step sizes of 0.01s for the fast phase and 0.1s for
the slow phase were used to speed up the CUEP calculation.
The two-time-scale stability assessment based on CUEP
guarantees stability for clearing times lesser than 0.21s while
the traditional CUEP method guarantees stability for clearing
times lesser than 0.16s, clearly a significant improvement
(30%) in the CCT estimation. Fig. 6 shows the projected
trajectory for a clearing time of 0.277s. The shaded area in
fig. 6 represents the connected component Sc(Vcr) of the
critical level set that contains the stable equilibrium point
(xs,zs). The post-fault trajectory was monitored until it enters
Sc(Vcr).

Fig. 5 – Time domain simulation of contingency scenario #1. While the fast B. Contingency Scenario #2
variables, the rotor angle and rotor speed, oscillate with a frequency close to In this scenario, xline=1.40, Pm=0.815, Qpre=0.1187, D=0.0588
1Hz, the slow variable E’q takes about 7 seconds to settle down.
and Efd=1.32. Three lines trip out and after 0.4085s one of
A. Contingency Scenario #1 them reconnects automatically. The point (xo,zo) in Fig. 7
In this scenario, four lines with reactance xline=1.50 p.u. are represents the pre-fault operating point. The fault excites the
connecting the generator to the infinite bus. The input fast variables. The fast system is unstable and the fault is
mechanical power is Pm=0.60 pu and D=0.0294. The reactive cleared to guarantee the first-swing stability of the fast
power that is injected at the terminal generator bus at pre-fault system. However, the fault is cleared close to the critical
equilibrium point is Qpre=0.0684p.u. and Efd=1.1821. A three- clearing time and the post-fault trajectory is pulled by the
phase short-circuit occurs at the generator terminal bus and slow dynamics in the neighborhood of the unstable

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 8

component Гu of Г. As a consequence of that, the slow dynamic phase takes place with the fault-on trajectory slowly
variable is excited and when the fast dynamic brings the sliding along the slow manifold in the neighborhood of ΓsF.
trajectory back to the stable component Гs, this component is This slow phase is a consequence of slow flux decay. The
reached outside the stability region of the asymptotically SEP trajectory reaches the neighborhood of a singularity in ΓsF and
of the slow system. The result is the lack of synchronism that an unstable fast dynamic phase is about to initiate when the
in reality was ignited by a slow flux decay of the generator. third line reconnects at t=28s. The post-fault trajectory
The flux decay, on the other hand, was a result of the fast initiates with a fast phase that brings the trajectory back to the
variables that drove the system trajectory to the neighborhood neighborhood of the stable component Гs of Г. The post-fault
of Гu. trajectory “reaches” the stable component Гs inside the
stability region of the SEP (xs,zs) of the slow system. As a
consequence, the slow phase of the post-fault trajectory slides
along the slow manifold of the post-fault system bringing the
system back to a stable equilibrium point.
Since the fast fault-on system is stable, the two-time-scale
algorithm searches the CUEP along the pseudo-fault-on
trajectory along Гs. In this example, the slow fault-on
trajectory reaches a singularity on ΓsF before the pseudo-fault-
on trajectory reaches the exit point of the slow system. The
instant at which the slow fault-on trajectory reaches a
singularity is an approximation of the corresponding CCT. In
this case, the CCT estimated by the proposed two-time scale
algorithm equals 25s while the real CCT is within the interval
(28s,29s). The CCT estimation obtained with the traditional
CUEP method is very conservative and equals 0.51s.
Fig. 7. Illustration of contingency scenario #2. The fast fault-on system is
unstable. The two-time scale CUEP algorithm uses the fast-fault-on trajectory to
calculate the fast CUEP. The slow trajectory of the slow system starting at the
fast does not converges to the CUEP because the CUEP of (Σε) belongs to the
stable component Γs. The pseudo slow post-fault trajectory associated with the
slow post-fault trajectory reaches the exit point of the slow system and then the
slow CUEP is calculated.
In this scenario, the CUEP belongs to the stable component
Γs of Γ. The proposed two-time scale algorithm calculates the
correct CUEP. Since the fast fault-on system is unstable, it
first obtains the fast CUEP (xo,zcoF). Fig. 7 shows the location
of the fast CUEP and the corresponding fast sustained fault-
on trajectory. Starting from (xo,zcoF) ∈ Гu, the slow post-fault
trajectory along Гu and the corresponding pseudo-trajectory
along Гs are calculated. In this case, the exit-point of the slow
system along the pseudo-trajectory is reached before the slow
system trajectory on Гu reaches an equilibrium point or
singularity. Fig. 7 illustrates the location of the CUEP of the Fig. 8. Illustration of contingency scenario #3. The fast fault-
slow system. In this case, the CUEP of the slow system is the on system is stable. The two-time-scale CUEP algorithm
CUEP of (Σε) for sufficiently small ε. searches the CUEP of the slow system in the direction of the
The two-time-scale algorithm for stability assessment pseudo slow fault-on trajectory. However, a singularity is
based on CUEP theory guarantees the stability for clearing reached on ΓsF before the pseudo trajectory reaches the exit-
times lesser than 0.37s while the traditional CUEP method point of the slow system. The instant at which the slow fault-
cannot be employed because the initial condition (xo,zo) is on trajectory reaches the singularity is an estimate of CCT.
outside the critical level set Sc(Vcr) (the shaded area) as shown
V. CONCLUSIONS
in fig. 7.
In this paper, the foundations of the two-time-scale CUEP
C. Contingency Scenario #3 method were provided. By studying a relationship between the
In this scenario, xline=1.50, Pm=0.815, Qpre=0.1276, stability boundary of the real system and the stability
Efd=1.33 and D=0.0588. At time t=0, two lines trip out boundaries of the slow and fast subsystems, a conceptual two-
exciting the fast dynamic. The fast dynamic is stable and the time-scale CUEP algorithm was proposed.
trajectory rapidly approaches the stable component ΓsF of the The proposed algorithm has the following advantages
constraint set of the fault system (see fig. 8). Then a slow when compared to the traditional CUEP method: (i) the

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.
 9

CUEP calculation is divided in fast and slow phases, (ii)

appropriate step sizes for each phase speed up the CUEP
calculation, (iii) CUEP calculation is robust, (iv) critical
clearing time estimates are much less conservative than those
obtained with the traditional CUEP method, especially when
the unstable mode is associated with the slow dynamic, (v) a
more comprehensive understanding of unstable modes is
provided. The proposed algorithm was tested in a simple
power system confirming all the aforementioned advantages.
Further research is still necessary to apply the proposed
conceptual algorithm to larger power system models. One
issue is to consider more comprehensive models. For this
purpose, slow and fast variables must be identified and energy
functions must be found. Another issue is the development of
efficient numerical algorithms to implement the conceptual
algorithm presented in this paper. The calculation of pseudo-
trajectories is an issue that deserves further discussions.

VI. REFERENCES
[1] H. D. Chiang, M. W. Hirsch and F. F. Wu, “Stability region of nonlinear
autonomous dynamical systems”, IEEE Trans. on Automatic Control, vol.
33, n. 1, pp. 16-27, Jan 1988.
[2] H. D. Chiang, F. F. Wu and P. P. Varaya, “Foundations of direct methods
for power system transient stability analysis”, IEEE Trans. on Circuits
and Systems, vol. CAS-34, n.2, Feb 1987.
[3] J. Palis and W. Melo, Geometric Theory of Dynamical Systems: An
Introduction. New York: Springer-Verlag, 1981.
[4] V. Venkatasubramanian, H. Schattler and J. Zaborsky, “Dynamic of large
constrained nonlinear systems - a taxonomy theory”, Proc. of the IEEE,
vol. 83, n. 11, pp. 1530-1560, Nov 1995.
[5] Y. Zou, M. H. Yin, and H. D. Chiang, “Theoretical foundation of the
controlling uep method for direct transient stability analysis of network-
preserving power system models”, IEEE Trans. on Circuits and Systems,
vol. 50, no.10, pp.1324-1336, October 2003.
[6] L. F. C. Alberto, “Characterization and Estimates of the Stability Region
of Nonlinear Dynamical Systems (in portuguese)”, Livre-Docencia Thesis
at EESC-University of Sao Paulo (USP), Brazil, 2006.
[7] L. F. C. Alberto, H. D. Chiang, “Two-Time Scale Energy Functions for
Stability Analysis of Nonlinear Systems”, submitted for possible
publication in the IEEE Trans. on Circuits and Systems I, 2007.
[8] H. D. Chiang, F. F. Wu and P. P. Varaya, “A BCU Method for Direct
Analysis of Power Systen Transient Stability”, IEEE Trans. on Power
Systems, vol. 9, n.3, Aug 1994.
[9] M. A. Pai, Power System Stability. Amsterdam, The Netherlands: North-
Holland, l98l.
[10] K. W. Chang, “Two Problems in Singular Perturbations of Differential
Equations”, Journal of Australian Mathematical Society, pp.33-50, 1969.

Authorized licensed use limited to: UNIVERSIDADE DE SAO PAULO. Downloaded on August 28,2021 at 23:34:30 UTC from IEEE Xplore. Restrictions apply.