Reachability Analysis of Nonlinear Systems Using Trajectory Piecewise

Linearized Models
Zhi Han and Bruce H. Krogh

Abstract— This paper presents a method to compute over- A number of reachability analysis methods for nonlinear
approximations of the reach set for nonlinear dynamic systems dynamic systems have been reported in the literature. The
using trajectory piecewise linearized (TPWL) models. Based methods presented in [1], [10] partition the state space into
on existing methods for reachability analysis of linear dynamic
systems, the method makes it possible to analyze high-order a finite set of cells and then approximate the nonlinear
nonlinear dynamic systems. The effectiveness of the method dynamics by linear or constant differential inclusions for
is illustrated by transient stability assessments of a three- each cell. For high-order models the number of cells is
machine electric power network and a ten-machine electric prohibitively large; construction of the partition consumes a
power network. large amount of computation time and memory. The method
I. I NTRODUCTION presented in this paper avoids the partition of the whole
state space by constructing the TPWL model on the fly; the
Reachability analysis is an essential problem in formal linearized models are computed as the reach set computation
verification of continuous and hybrid models of embedded is carried out.
control systems. The aim of reachability analysis is to This paper is organized as follows. Section II introduces
compute the set of reachable states (the reach set) in the state basic concepts and notation. Section III presents the main
space for a given model and a set of initial states. Verification approach for reach set computation of nonlinear systems.
results are then obtained from the computed reach set. The Section IV discusses the problem of estimating bounds
application of reachability analysis has been limited by the on the error caused by piecewise linearization. Section V
order of the state space of the model, but recent use of presents the application of the proposed method to assess
reduced-order models [9], [1] and use of zonotopes [6] have the first-swing stability of two electric power systems, one
made it possible to compute reach sets for relatively high- with 3 generators (a 6th -order nonlinear system) and one
order linear dynamic systems. with 10 generators (a 20th -order nonlinear system). The
Reachability analysis of nonlinear systems remains a diffi- concluding section summarizes the contribution of the paper
cult problem for systems of any order. Unlike linear systems, and discusses current research directions.
nonlinear systems can exhibit very complex behaviors for
even low-order models, which makes it difficult to analyze II. P RELIMINARIES
them using current computational methods. In addition, non- This paper considers nonlinear dynamic systems modeled
linear models from practical applications can be high-order, by nonlinear state equations of the form
e.g., the nonlinear models of electric power networks, which
makes the reach set computation even harder, if feasible at S : ẋ = f (x), (1)
all. where the set of initial states is given by X0 = {x|x =
In this paper, we propose the use of trajectory piecewise Φ0 w + γ0 , w ∈ Pw } ⊆ Rn , Φ0 ∈ Rn×p , γ0 ∈ Rn and
linearized (TPWL) models to extend computational methods Pw ⊆ Rp , p ≤ n is a polytope satisfying 0 ∈ Pw and
for reachability analysis of linear systems to nonlinear sys- supw∈Pw kwk ≤ 1. The norm k · k used in this paper is
tems. TPWL models are piecewise linear models obtained the infinity norm. We represent X0 with hΦ0 , γ0 , Pw i, which
by linearizing the nonlinear dynamics of a system at several we call the affine representation. Note that γ0 ∈ X0 . We
points along a trajectory or a set of trajectories [16], [5]. call γ0 the nominal initial state. Also note that any polytope
Our method works as follows. First a nominal initial state is can be converted to the appropriate affine representation by
picked from the initial set. Then a TPWL model is created choosing appropriate γ0 and Φ0 .
along the nominal trajectory and the reach set is computed Given a polytope hΦ, γ, Pw i, if a set X satisfies
using the TPWL model. The error introduced by linearization dist(X , hΦ, γ, Pw i) ≤ δ, where dist(·, ·) is the Hausdorff
is then incorporated by ‘bloating’ the reach set to guarantee distance, then X ⊆ hΦ, γ, Pw i ⊕ Bδ , where ⊕ is the
conservativeness [15], [7], [1]. The recent advances in reach- Minkowski sum and Bδ = {x|kxk ≤ δ, x ∈ Rn } is
ability analysis for linear systems makes it possible to apply a square hyper-box. The size δ of the hyper-box Bδ is
this approach to high-order nonlinear system models, e.g., called the bloating factor. To simplify the notation we
the ten-machine power network in section V. represent hΦ, γ, Pw i ⊕ Bδ with hΦ, γ, Pw , δi, which we call
the approximate affine representation.
Z. Han and B. H. Krogh are with Department of Electrical and Computer
Engineering, Carnegie Mellon University, 5000 Forbes Ave., Pittsburgh, PA We assume that (1) has a unique solution for every initial
15213. Email: {zhih|krogh}@ece.cmu.edu state x(0) ∈ X0 . The reach set at time t of the nonlinear
dynamic system S for a set of initial states X0 is defined as Reach sets for linear systems can be computed effectively
the set Reach(S, X0 , t) = {x|x = χ(t) where χ : [0, t] → using the matrix functions ϕi . For linear system ẋ =
Rn , χ(0) ∈ X0 and ∀s ∈ [0, t] : χ̇(s) = f (χ(s))}. The reach Ax + b where the initial set is hΦ0 , γ0 , Pw i, the reach
set of the system over the time S interval [t1 , t2 ] is defined set is Reach(t) = hϕ0 (A, t)Φ0 , γ0 + tϕ1 (A, t)b, Pw i. The
as Reach(S, X0 , [t1 , t2 ]) = t∈[t1 ,t2 ] Reach(S, X0 , t). We computation of reach set is performed in two steps: first the
simply write Reach(t) and Reach([t1 , t2 ]) if other param- matrix functions are computed, then affine representations
eters are clear from the context. For a given time sequence for the results are computed. With the aid of numerical
t1 , t2 , . . ., tk , the reach set at a time point tk is written as computation methods, the matrix functions can be evaluated
Xk as a short-hand for Reach(tk ). effectively for moderate-size to large-size matrices [12], [17],
which in turn enables efficient computation of reach sets for
A. TPWL Models
large-scale models [8].
For system S in (1), let γ(·) be the nominal trajectory, i.e., The following theorem estimates bloating factors for affine
the trajectory from the nominal initial state γ0 . Choosing n systems.
points on the trajectory γ1 = γ(t1 ), . . ., γn = γ(tN ) where Theorem 1 (Bloating factor for affine systems [8]):
0 = t0 ≤ t1 ≤ . . . ≤ tN , a TPWL model of S along Given an affine system ẋ = Ax + b with the initial
the nominal trajectory γ(·) is the linear time-varying system state set X0 = hΦ0 , γ0 , Pw i and t0 , t1 ∈ R satisfying
defined as 0 = t0 < t1 = h, let P0,1 be a polytope such that
¯
∂f ¯¯ CH(X0 ∪ X1 ) ⊆ P0,1 , where CH denotes the convex hull
Sγl : ẋ(t) = (x(t)−γk−1 )+f (γk−1 ), t ∈ [tk−1 , tk ]
∂x ¯γk−1 of the set argument, then Reach([t0 , t1 ]) ⊆ P0,1 ⊕ Bδ where
(2)
where xk−1 = γk−1 , 0 = t0 < t1 < . . . < tN , k = δ = 2h2 ψ2 (A, A2 Φ0 , h) + 2h2 ψ2 (A, A2 γ0 + Ab, h).
1, . . . , N .
Throughout the paper we neglect the numeric error caused III. R EACH S ET C OMPUTATION U SING TPWL M ODELS
by round-off and focus on the error introduced by lineariza-
This section presents a method for computing over-
tion. In practice the TPWL models are generated on the
approximations of the reach set of a nonlinear dynamic
fly. Initially t = t0 and γkl = γ0 . At each time tk the
system using TPWL models. The approach consists of two
nominal state of the system γkl is chosen and then the
steps. First the reach set for the TPWL model is computed.
TPWL model Sγl l is computed for the point γkl . Sγl l is
k k Then the discrepancy between the TPWL and nonlinear
then used to simulate the trajectory from tk to tk+1 where
model is taken care of, where the reach set is bloated to
hk+1 = tk+1 − tk is the step size. Then, starting from this
incorporate the linearization error.
result, we repeat the computation for the next step. Thus
First we consider the problem of computing reach sets Xk ,
the TPWL model is generated along γ l (·), a trajectory of
k = 1 . . . N for a time sequence 0 = t0 < t1 < . . . < tN .
the TPWL model instead of along the state trajectory γ(·)
Figure 1 shows a procedure to compute the reach set for a
of the nonlinear model. As a result the error caused by
nonlinear dynamic system S and a set of initial states. The
the linearization accumulates over time. It might seem more
reach set Xk is computed using a time-stepping strategy, i.e.,
attractive to compute the TPWL along the nominal trajectory
by computing Xk , k = 1, . . . , N sequentially. In each step
of the nonlinear model to avoid the accumulation of error.
an linear dynamic system model is created by linearizing the
The nominal trajectory cannot be computed exactly, however.
nonlinear dynamics at points on the nominal piecewise linear
Instead it is obtained using numerical integration [4], and it
trajectory γkl . The linearized model is then used to compute
is shown in [11] that by choosing appropriate step sizes hk
the reach set for the next step.
and linearization schemes, the integration method using the
TPWL model achieves the same order of accuracy as most Procedure REACH TPWL
of the existing numerical integration methods. Input: S, X0 = hΦ0 , γ0 , Pw i, tk , k = 1, . . . , N
Output: Over-approximation of Xk
B. Reach Set Computation for Linear Systems Procedure:
t ← 0, k ← 1
Following [17], we introduce several matrix functions WHILE n < N
Let step size hk ← tk − tk−1
to facilitate the discussion. Denote the matrix expo- Compute linearized model at γk−1 l

nential function eAt as ϕR0 (A, t) and define the matrix Compute reach set using TPWL model Xk
t = hΦk , γk , Pw i
functions ϕ1 (A, t) = 1t 0 eA(t−τ ) dτ and ϕ2 (A, t) =
1
R t A(t−τ ) Estimate bloating factor ²k
t2 0 e τ dτ. It can be verified that tϕ1 (A, t)b and k ← k + 1, t ← tk
END WHILE
t2 ϕ2 (A, t)b where b ∈ Rn are the state responses of the RETURN hΦk , γk , Pw , ²k i, k = 1, . . . , N
linear dynamic system ẋ = Ax + bu with zero initial state
to the step input and ramp input, respectively. For matrices Fig. 1. Computing reach sets using TPWL models.
A ∈ Rn×n and B ∈ Rn×m , define the scalar-valued matrix
function ψi (A, B, t) = sups∈[0,t] kϕi (A, s)Bk, i = 0, 1, 2, The computation of the reach set for the TPWL model
which can be computed by maximizing the matrix norm over is similar to the reach set computation of linear systems
[0, t]. l
[6], [8]. Assume Xk−1 = hΦlk−1 , γk−1 , Pw i, then Xk =
hΦlk , γkl , Pw i, where bloating factors ²k are based on estimations of the approxi-
mation errors of each segment in the TPWL models. We first
Φlk = ϕ0 (Ak−1 , hk )Φlk−1 ,
compute a bound on the error incurred in each step of the
γkl l
= γk−1 l
+ hk ϕ1 (Ak−1 , hk )f (γk−1 )
computation, which we call the one-step error (OSE), and
and hk = tk − tk−1 . The difference between the TPWL then derive a bound on the accumulated error (AE) which
model and the linear model is that a new matrix Ak−1 is includes the error propagated from the previous steps.
computed at each time step along with the matrix functions
ϕ0 and ϕ1 . For linear time-invariant systems with an evenly A. OSE for Reach Sets
spaced time sequence, the matrix functions need to be
computed only once. The estimation of the bloating factor The OSE of TPWL approximation is caused by approxi-
²k is based on estimating the error introduced by piecewise mating the nonlinear dynamic model with a linearized model.
linearization. Section IV discusses the estimation of ²k . Consider a segment [tk , tk+1 ] and assume that the nominal
l
Next we consider the problem of computing reach set for initial state is γk−1 . Without loss of generality, we let k = 0,
[0, tf ] where tf ≥ 0. The computation of Reach([0, tf ]) is γkl = γ0 and let h = t1 − t0 .
built upon the results of the REACH TPWL procedure. We The nonlinear system is given by
first choose a time sequence 0 = t0 < t1 < . . . < tN = tf
and use REACH TPWL to compute Xk , k = 0, . . . , N . ẋ = f (x), t ∈ [0, h] (3)
The reach set Reach([0, tf ]) is computed as the union
SN where x(0) ∈ X0 . The TPWL approximation to the system
Reach([0, tf ]) = k=1 Reach([tk−1 , tk ]) where each reach
set segment Reach([tk−1 , tk ]) is computed from Xk−1 and trajectory is
Xk . Figure 2 shows the procedure to compute Reach([0, tf ]) ẋl = A(xl − γ0 ) + f (γ0 ) (4)
from the results of REACH TPWL. ¯
¯
Procedure REACH SEG where xl (0) = x(0) and A = ∂f∂x ¯ . Let eOSE (t) = x(t) −
γ0
Input: S, t0 ,. . .,tN , X0 ,. . ., XN
Output: Over-approximation of Reach([t0 , t1 ]), . . ., Reach([tN −1 , tN ]) xl (t) denote the OSE of TPWL approximation. The dynamic
Procedure: equation for eOSE (t) is obtained from (3) - (4)
FOR k ← 1 TO N
Compute a convex polytope Pk−1,k
such that CH(hΦk−1 , γk−1 , Pw i ∪ hΦk , γk Pw i) ⊂ Pk−1,k ėOSE = −A(xl − γ0 ) + f (x) − f (γ0 ), eOSE (0) = 0 (5)
Estimate segment bloating factor δk−1,k using
Theorem 1
Rk ← Pk−1,k ⊕ Bδk−1,k +²k Let r denote the difference of the nonlinear function f and
END FOR its linearization, i.e., r(x) = f (x)−[A(x−γ0 )+f (γ0 )]. Then
RETURN Rk , k = 1, . . . , N .
we have f (x) = [A(x − γ0 ) + f (γ0 )] + r(x). Substituting
Fig. 2. Computing reach set segments. the nonlinear function f in (5) we have

For each segment, the procedure REACH SEG consists ėOSE = AeOSE + r(x), eOSE (0) = 0 (6)
of three steps. First a polytope over-approximation of the
low-dimensional polytope approximations of the reach set Let R ⊆ Rn be a closed region containing the reach
is computed as the convex hull CH((hΦk−1 , γk−1 , Pw i ∪ set Reach([0, h]), and let d = maxx∈R kf (x) − [A(x −
hΦk , γk Pw i). Second the low-dimensional convex hull is γ0 )+f (γ0 )]k. The following proposition establishes an upper
bloated to enclose the reach set of the TPWL model in bound on keOSE (t)k for all t ∈ [0, h] and x(0) = xl (0) ∈
[tk , tk+1 ]. Third the set is bloated to enclose the reach set of X0 .
the nonlinear model in [tk , tk+1 ]. The first two steps consider Proposition 1 (OSE for reach sets): For eOSE defined by
the reach set segment of the TPWL model while the third (6), ∀t ∈ [0, h] : keOSE (t)k ≤ ψ0 (A, I, h)dh, where I is the
step takes care of the error introduced by linearization. The identity matrix with appropriate dimension.
final result is a set of polytopes bloated by the sum of the two Proof: Noticing that eOSE (·) is the state trajectory of
bloating factors computed in steps two and three such that the response of the linear Rsystem (6) with bounded input.
t
they enclose the reach set segments of the nonlinear model. Therefore keOSE (t)k ≤ 0 keAt kds · kr(x)k. From the
The bloating factor for the reach set segment consists of definition of d we know sup R x∈R kr(x)k ≤R d. The conclusion
t t
two components, the bloating factor δk−1,k for the segment follows from the relation 0 keAt kds ≤ 0 ψ0 (A, I, t)ds =
for the TPWL model and the bloating factor ²k for the tψ0 (A, I, t).
nonlinear model. By Theorem 1, the bloating factor δk−1,k As a corollary, the following proposition asserts that the
for the approximation error can be estimated using the matrix reach set of the TPWL model is the same as the reach set
functions ψ2 . The next section focuses on ²k . of the nonlinear model if the nonlinear equation is linear in
the region R.
IV. E STIMATING B LOATING FACTORS Corollary 1: For the notation introduced above, if ∀x ∈
Bloating is needed to compensate for the discrepancy R : f (x) = A(x−γ0 )+f (γ), then ∀t ∈ [0, h] : keOSE (t)k =
between the TPWL model and the nonlinear model. The 0.
B. AE for Reach Sets Reach(Sγl l , X0 , tk ), then Reach(S, X0 , tk ) ⊆ Pk ⊕ B²k
Suppose the over-approximation of the reach set for Xk−1 where ²k is defined as:
is hΦk−1 , γk−1 , Pw , ²k−1 i. For step k, the nonlinear system ²0 = 0
is given by
²k = ψ0 (Ak−1 , I, hk )(dk hk + ²k−1 ), k ≥ 1
ẋ = f (x), t ∈ [tk−1 , tk ]
where hk = tk − tk−1 .
l
where x(tk−1 ) = γk−1 + Φk−1 w + eAEk−1 , w ∈ Pw and Proof: Inductively invoke proposition 2.
eAEk−1 is the AE in the last step satisfying keAEk−1 k ≤ If the nonlinear function is linear in each region Rk , then
²k−1 . The TPWL model is given by as stated in the following corollary, the reach set at a point
in time Xk computed for the TPWL model is exact.
ẋl = f (xk−1 ) + Ak−1 (xl − xlk−1 ), t ∈ [tk−1 , tk ]
Corollary 2: For the notation introduced above, if ∀k ≤
¯
l ∂f ¯ N : ∀x ∈ Rk : f (x) = Ak (x − γkl ) + f (γkl ), then ∀t ∈
where x(tk−1 ) = γk−1 +Φk−1 w, w ∈ Pw , Ak−1 = ∂x ¯ .
γk−1 [0, tN ] : keAE (t)k = 0.
The difference in the initial condition of the two problems is Proof: Prove by induction on k.
due to the accumulation of linearization error from previous For the reach set segment, the bloating factor δk−1,k from
steps. Let eAE (t) = x(t)−xl (t), t ∈ [tk−1 , tk ] denote the AE theorem 1 is combined with the bloating factor ²k . The
for the current segment. The dynamic equation for eAE (t) following corollary proves the correctness of the procedure
over [tk−1 , tk ] is then given by REACH SEG.
l l l Corollary 3: For the notation used in theorem 2, let
ėAE = f (x)−f (γk−1 )−Ak−1 (x −γk−1 ), keAE (tk−1 )k ≤ ²k−1 .
P be a polytope over-approximation to the reach set
l
Rewriting the function f in terms of r, the dynamic equation Reach(Sγ l , X0 , [tk−1 , tk ]), then Reach(S, X0 , [tk−1 , tk ]) ⊆
of the error term eAE is given by P ⊕ B²k +δk−1,k .
This section presents method to compute a conservative
ėAE = Ak−1 eAE + r(x), t ∈ [tk−1 , tk ] (7) bloating factor ²k using a bound d on the nonlinearity of
the system, which can be determined prior to the reach set
where keAE (tk−1 )k ≤ ²k−1 and r(x) is the linearization computation [5], [15], [7].
error.
Let R ⊆ Rn be a closed region containing the reach set V. C ASE STUDY
Reach([tk−1 , tk ]) and d = supx∈R kr(x)k. The following This section applies the reachability analysis method to
proposition gives a bound on the AE for the reach sets of transient stability assessment (TSA) of electric power net-
step k. works. TSA concerns evaluating the capability of a power
Proposition 2 (Total error for reach set): For eAE de- network to return to a stable operating state without loss of
fined by (7), ∀t ∈ [tk−1 , tk ] : keAE (t)k ≤ ψ0 (A, I, h)(dh + synchronization following a system fault, e.g., a short circuit
²k−1 ), where h = tk − tk−1 and I is the identity matrix with caused by lightning. The most widely used techniques for
appropriate dimension. TSA involve time-domain simulation. In these simulations it
Proof: The response of the system (7) can be decom- is assumed that the pre-fault system is at its steady state.
posed into two parts Then simulation is carried for the faulted network until
Z t the clearing time, tcl , the time when the short circuit is
A(t−tk−1 ) A(t−s) eliminated. Then starting from the state at tcl , the post-
eAE (t) = eAE (tk−1 )e + e ds.
tk−1 fault network is simulated to see if the system is stable.
The second term is the OSE, which is bounded by The system is said to be first-swing stable if the rotor angles
ψ0 (A, I, h)dh. Consider the first term: since exhibit the first swing, i.e., the derivative of the rotor angle
becomes less than zero at some point in time [18]. Instead
eAE (tk−1 )eA(t−tk−1 ) ≤ ²k keA(t−tk−1 ) k ≤ ²k ψ0 (A, I, h), of assuming the system initially at a steady state, we apply
reachability analysis to assess first-swing stability for a set
the proposition follows from the sum of the above two of possible initial states. We choose the pre-fault steady state
bounds. to be the nominal initial state.
Let Rk be the closed set containing the reach set The following model is used for each generator in the
Reach(tk−1 , tk ), and dk = maxx∈Rk kf (x) − [Ak−1 (x − power network:
γk−1 ) + f (γk−1 )]k k = 1, . . . , N . The following theorem
δ̇i = (ωi −“ ω0 )
establishes the bloating factor for the reach set segments for P
ω̇i = M1 Pmi − D(ωi − ω0 ) − n+1 Vi Vj bij sin(δi − δj )
nonlinear dynamic systems, thus proving the conservative- P i j=1,j6
” =i
− n+1j=1,j6=i Vi Vj gij cos(δi − δj )
ness of the procedure REACH TPWL.
Theorem 2 (Bloating factor for nonlinear systems): where Mi is the moment of inertia, Di is the damping
For nonlinear system (1) where the initial set is constant, bij and gij are the transfer conductance of the
X0 = hΦ0 , γ0 , Pw i. Assume the TPWL approximation network-reduced admittance matrix and Pi is the power
is (2) let Pk be a polytope over-approximation of injection at bus i.
 In TSA of power networks, the relative angles and speeds 0.095

are defined with respect to a selected reference angle. In 0.0879

δ (In COI reference)

δ (In COI reference)
principle, any rotor angle may be taken as the reference 0.0878 0.09

angle. A commonly chosen reference is the P center-of-inertia 0.0877

n 0.085

(COI) angle, defined as δCOI = Pn 1 Mi i=1 Mi δi . All 0.0876

2
2
i=1
0
the relative bus angles are defined as δi = δi − δCOI .
0.0875
0.08

0.0874
Power networks in TSA are high-order nonlinear systems. 0 0.05 0.1
t (seconds)
0.15 0.2 1 2 3
t (seconds)
4 5

Due to the limitations of the existing tools, reachability (a) Faulted (b) Post-fault
analysis has been applied to very small power networks, e.g.,
Fig. 3. Reach set of the rotor angle at bus 2.
a one-machine-infinite-bus network modeled as a second-
order system [14]. In this section we apply our method of −3
x 10
−3
x 10

reachability analysis to a 3-machine 3-bus network (a 6th -

5 6

4
order nonlinear system) and a 10-machine 39-bus network 4

ω2 (In COI reference)

ω (in COI reference)

(a 20th -order nonlinear system).
2
3
0
2
−2

2
A. 3-Machine 3-Bus Network 1
−4

0 −6
0 0.05 0.1 0.15 0.2 1 2 3 4 5
A three-bus network from [2] is used in this section to il- t (seconds) t (seconds)

lustrate the TPWL reachability analysis method. The electric (a) Faulted (b) Post-fault
power network consists of three buses. At each bus there is Fig. 4. Reach set of the rotor speed at bus 2.
a generator providing power and a constant load impedance.
The nonlinear dynamic model
£ is sixth-order. The state
¤Tvector
of the system is x = δ1 δ2 δ3 ω1 ω2 ω3 . The B. 10-Machine 39-Bus Network
steady-state of the power network is solved from the equa-
The 10-machine 39-bus IEEE benchmark system is repre-
tions δ̇i = ω̇i = 0, i = 1, . . . , 3. The pre-fault steady state of
sentative of the New England power network and has been
the three-bus example are δ1s = 0.4615, δ2s = 0.2619, δ3s = 0
used in many case studies, e.g., [2], [13]. The parameters of
and ω1s = ω2s = ω3s = 1.
the network and machines can be found in the appendix of
In TSA we consider the case that at t = 0 bus 2 is subject [13].
to a fault which is then cleared at t = 0.2 sec without We use second-order models for all the machines in the
changing the configuration of the network. The initial set network. The nonlinear dynamic model consists of 20 contin-
is chosen to be δ1 (0) ∈ [δ1s − 0.01, δ1s + 0.01], δi (0) = δis , uous state variables, two for each machine. The fault under
i = 2, 3 and ωi (0) = ωis (0), i = 1, 2, 3. consideration is a shortage at bus 35, which is connected
Reach set segments are computed using REACH SEG. directly to generator 6. The fault clearing time is assumed to
In this case study the linearization error dk is approxi- be 0.2 sec as suggested in [2]. The safety property to verify
mated by the maximum over the reach set at tk−1 : d¯k = is that the system is first-swing stable for the fault under the
maxx∈Pk−1 kf (x) − [Ak−1 (x − γk−1 ) + f (γk−1 )]k. k = uncertain initial condition δ6 (0) ∈ [δ6s − 0.8 × 10−3 , δ6s +
1, . . . , N . The maximum is computed numerically at each 0.8 × 10−3 ], ω6 (0) ∈ [ω6s − 0.2 × 10−4 , ω6s + 0.2 × 10−4 ].
time step. The global approximation error ²k is estimated as All other rotor angles and speeds are assumed to be at the
²k ≤ 1.78 × 10−5 . pre-fault steady state.
The reach set computed for the faulted and post-fault The reach set is computed using the TPWL models. The
system are shown in Figures 3 and 4. Figure 3 shows the reach set of the COI-referenced machine angle and speed
reach set for the rotor angle. Figure 4 shows the reach set at the faulted bus are shown in Fig. 5 as the dotted lines.
computed for the rotor speed deviation. All state variables are The response of the system with the nominal initial state
in COI framework of reference. The reach set for the faulted (the pre-fault steady state) is plotted in the same figure as
system at 0.2 sec is the set of initial states for the post-fault solid lines. The sets of ω6 reach below zero at 0.4 second.
system. The solid lines in the two figures are the nominal After examining the reach set of all machines, we claim the
trajectories of the faulted and post-fault system. It can be network for the set of initial state is indeed first-swing stable
seen that the nominal trajectory is contained in the reach for this fault.
set. During the on-fault mode, the rotor speed at the faulted
bus increases, causing the rotor angle to increase. After the C. Comparison with CheckMate
fault is cleared the rotor angle and speed oscillate. Starting The reach set computation procedures REACH TPWL and
from any state in the initial set, the trajectories of ω2 reach REACH SEG are implemented in M ATLAB. Among the
below 0, as observed in Fig. 4. The reach set of the other two existing reachability analysis tools, CheckMate [3] is the
machines are examined in the same way. After examining the most comparable implementation since it handles nonlinear
reach set of all machines, we claim the network for the initial systems and it also uses polyhedral approximation for reach
set is indeed first-swing stable for this fault. sets.
 −3
x 10
0.215 2.5
ACKNOWLEDGMENTS
2
0.21
This research was supported in part by the US Defense

ω (in COI reference)

δ (in COI reference)

1.5
0.205
1
Advance Projects Research Agency (DARPA) contract nos.
0.2
0.5
F33615-00-C-1701 and F33615-02-C-4029, US Army Re-
search Office (ARO) contract no. DAAD19-01- 1-0485, and

6
6

0.195
0

0.19 −0.5
the US National Science Foundation (NSF) contract no.
0 0.1 0.2
t (seconds)
0.3 0.4 0 0.1 0.2
t (seconds)
0.3 0.4
CCR-0121547.
(a) Reach sets of δ6 (b) Reach sets of ω6 We thank Professor Marija D. Ilic for helpful discussions
of power systems and suggestions on the power system
Fig. 5. Reach set of the 10-machine 39-bus network.
modeling and the transient stability assessment of power
networks.
Table I shows a comparison of reachability analysis using R EFERENCES
TPWL approach and the adapted CheckMate method. The [1] E. Asarin, T. Dang, and A. Girard. Reachability analysis of nonlinear
computations are performed on a Pentium4 PC with 1G systems using conservative approximation. In Hybrid systems: com-
RAM running Windows XP and M ATLAB 7.0.1. Note that putation and control (HSCC’2003), pages 20–35, 2003.
[2] T. Athay, R. Podmore, and S. Virmani. A practical method for the
the memory usage is the total memory used by M ATLAB, direct analysis of transient stability. IEEE Transaction on power
which by itself consumes approximately 68 Mbyes to load. apparatus and systems, PAS-98(2):573–584, March/April 1979.
[3] A. Chutinan and B. H. Krogh. Compuational techniques for hy-
It can be seen that the TPWL reach set computation is brid system verification. IEEE Transaction on Automatic Control,
more efficient in both cases for memory use and computation 48(1):64–75, Jan 2003.
time. The main reason is due to the fact that CheckMate only [4] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential
Equations. Prentice-Hall Series in Automatic Computation. Prentice-
allows full-dimensional polytopes and our method uses affine Hall Inc., 1971.
representation which can represent low-dimensional poly- [5] A. Girard. Approximate solutions of odes using piecewise linear
topes. For the 3-machine network the TPWL method is one vector fields. In 5th International Workshop on Computer Algebra
in Scientific Computing, 2002.
order of magnitude faster than CheckMate procedure. For [6] A. Girard. Reachability of uncertain linear systems using zonotopes.
the 10-machine network the CheckMate procedure causes In Hybrid Systems: Computation and Control (HSCC’05), pages 291–
M ATLAB to crash due to the large requirement of memory, 305, March 2005.
[7] M. R. Greenstreet and I. Mitchell. Reachability analysis using
while the TPWL method uses a little more memory than the polygonal projections. In Proceedings of the Second International
3-machine case and the computation is completed in seconds. Workshop on Hybrid Systems: Computation and Control, pages 103–
116, Berg en Dal, The Netherlands, Mar. 1999. Springer. LNCS 1569.
Memory (Mbytes) Computation Time (seconds) [8] Z. Han and B. Krogh. Reachability analysis of affine dynamic
CheckMate TPWL CheckMate TPWL systems using affine representation of polytopes. In Hybrid Systems:
3-machine 78.1 76.8 2.813 0.328 Computation and Control (HSCC’06), in preparation.
[9] Z. Han and B. H. Krogh. Reachability analysis of hybrid systems using
10-machine 644.7 78.2 Out of memory 2.45
reduced-order models. In IEEE Proc of American Control Conference
TABLE I (ACC’04), 2004.
C OMPARISON OF REACH SET COMPUTATION METHODS . [10] T. A. Henzinger, P.-H. Ho, and H. Wong-Toi. Algorithmic analysis
of nonlinear hybrid systems. IEEE transaction on automatic control,
43(4):540–554, 1998.
[11] B. V. Minchev and W. M. Wright. A review of exponential integrators
for first order semi-linear problems. Technical report, Norwegian
University of Science and Technology, 2005.
VI. D ISCUSSION [12] C. Moler and C. Van Loan. Nineteen dubious ways to compute the
exponential of a matrix, twenty-five years later. SIAM Review, 45(1):3–
This paper presents a method to compute over- 49, 2003.
approximations of the reach set for nonlinear dynamic sys- [13] M. A. Pai. Energy Function Analysis for Power System Stability.
Kluwer Academic Publisher, Boston, 1989.
tems using trajectory piecewise linearized (TPWL) models. [14] H.-L. Pei and B. H. Krogh. Stability regions for systems with
The method is illustrated for transient stability assessments mode transitions. In Proceedings of the American Control Conference
of a 6th -order and a 20th -order nonlinear models of electric (ACC’01), pages 4834–4839, 2001.
[15] A. Puri, V. Borkar, and P. Varaiya. ²-approximation of differential
power networks. inclusions. In Conference on Decision and Control, 1995.
The accuracy of the reach set over-approximation depends [16] M. Rewieński and J. White. A trajectory piecewise-linear approach
on how well the system dynamics can be approximated to model order reduction and fast simulation of nonlinear circuits and
micromachined devices. IEEE Transaction on Computer-Aided Design
by the TPWL model within the region that encloses the of Integrated Circuits and Systems, 22(2):155 – 170, February 2003.
reach set. The reach set computation can be exact if the [17] Y. Saad. Analysis of some krylov subspace approximations to the
dynamic system is indeed linear in the region. The method matrix exponential operator. SIAM Journal of Numerical Analysis,
20(1):209–228, 1992.
is useful when nonlinear systems can be well approximated [18] P. W. Sauer and M. A. Pai. Power System dynamics and stability.
by piecewise linear models in the neighborhood of the Pearson Education, Inc., 1998.
nominal trajectories. Empirical results show that such cases [19] S. K. Tiwari and R. A. Rutenbar. Scalable trajectory methods for
on-demand analog macromodel extraction. In Proceedings of Design
are common in some applications [19], [16]. Our current Automation Conference (DAC’1005), pages 403–408, June 2005.
research direction is on performing on-the-fly partition of
the reach set to reduce the approximation error.