\PPLIED MATHEMATICS

COMPUTATION

Stability indices for randomly perturbed power systems q

Humberto Verdejo a'*, Wolfgang Kliemann b, Luis Vargasc

a Department of Electrical Engineering, University of Santiago, Santiago, Chile b Department of Mathematics, Iowa State University, IA, USA c Department of Electrical Engineering, University of Chile, Santiago, Chile

ARTICLE INFO ABSTRACT

This paper considers the stability of moments of stochastic systems, such as stability of the mean or mean-square stability. The exponential growth behavior of moments is compared to almost sure exponential growth via Lyapunov exponents. We develop a series of indices that are useful to describe system performance under random perturbations. The theory is applied to two examples, including an electric power system.

© 2014 The Authors. Published by Elsevier Inc. All rights reserved.

1. Introduction

Stability at operating points is one of the key requirements of engineering systems. As long as the system is given by time-invariant dynamics, linearization at the operating point gives local stability information that can be extended through incorporating some nonlinear features, e.g., via the use of normal forms, see [1]. If the system under consideration has time-varying dynamics, the usual modal approach fails since for these systems eigenvalues do not describe the stability behavior of the linearized system. Therefore, one has to approach (exponential) stability directly via the Lyapunov exponents of the system at the operating point.

An important class of systems with time varying dynamics are those systems that are subject to sustained random perturbations, such as load behavior, environmental effects, or intermittent generation in power systems. The interaction between system dynamics and perturbation falls into two groups: (i) the random noise changes the operating point of the system, or (ii) the equilibrium point persists under all perturbations. We have developed performance indices for case (i) in [2], and analyzed one specific approach in case (ii) in [3] using almost sure Lyapunov exponents. This paper develops several performance indices for case (ii), analyzes their relationships, and compares the results for several examples. The key idea is the look at the sample (exponential) growth rates for trajectories and at the growth rates of moments of the trajectories, such as the stability of the mean, or mean square stability involving the second moment. Both points of view result in potentially useful performance criteria for power systems.

2. Mathematical background

2.1. The system model

We start from a nonlinear differential equation y(t) = f (y(t), n(t, m)) in Rd with sustained random perturbation n(t, m). In order to analyze optimal parameter settings for stability at an operating point, we linearize the system equations at the equilibrium point y*. Linearization (with respect to y) at the equilibrium results in the system

q This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-No Derivative Works License, which permits non-commercial use, distribution, and reproduction in any medium, provided the original author and source are credited. * Corresponding author. E-mail address: humberto.verdejo@usach.cl (H. Verdejo).

Contents lists available at ScienceDirect

Applied Mathematics and Computation

journal homepage: www.elsevier.com/locate/amc

CrossMark

Keywords: Stochastic system Stochastic stability Lyapunov exponents Moment stability Electric power systems

0096-3003/$ - see front matter © 2014 The Authors. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2014.01.024

x(t)=A(n(t, x))x(t) in Rd (1)

where A(n(t, x)) is the Jacobian of f (y(t), n(t, x)) at y*. We denote by u(t, x, x) the trajectories of (1) with initial value U(0, x, x)= x e Rd. We think of a given probability space (X, g, P) under the usual conditions on which the Wiener process in (2) is defined. We use the notation x e X, and all expectations E( ) are with respect to the given probability measure P.

The random perturbation can be considered as white noise, leading to a stochastic differential equation for (1), or as a colored, bounded noise. In this paper we discuss the latter situation since macroscopic perturbations in engineering systems generally are non-white; but a similar theory also holds for the white noise case, see [4] for the basics. We start from a background noise g, given by a stochastic differential equation on a compact C^-manifold M

dg = Xo(g)dt + ^X,-(g)o dWi on M (2)

where the vector fields X0,...,Xr are C1, and "o" denotes the Stratonovic stochastic differential. We assume that (2) has a unique stationary, ergodic solution g*(t, x) which is guaranteed by the condition (compare [5])

dimLA{X1,...,Xr}(9) = dimM for all 9 e M. (3)

Here LA{-} denotes the Lie algebra generated by a set of vector fields. The background noise g*(t, x) is mapped via a sur-jective smooth function f : M ! U c Rm, f (g) = n, into the system perturbation n(t, x). This setup allows great flexibility when modeling the statistics of the system noise.

2.2. Lyapunov exponents

Exponential stability of the system (1) is described by Lyapunov exponents; in [6] we gave an overview of the almost sure theory, with applications to power systems. Here we extend the analysis to moment Lyapunov exponents. The individual Lyapunov exponents of the trajectories u(t, x, x) of (1) are given as

k(x, x) = limsup1 log |u(t,x, x)|, (4)

and for p e R the Lyapunov exponent of the pth moment is given by

g(p,x) = limsup1 log E|u(t,x, x)|p. (5)

This includes for p = 1 the exponential growth behavior of the mean, and for p = 2 the exponential mean square stability of the system. We again need the projection of the linear system onto the sphere Sd-1 in Rd:

s(t) = h(f(t, x),s(t)), h({,s) = (A(f)-q(f,s))s, q(n,s)= sTA(n)s, (6)

where '' T'' denotes the transpose. via identification of s and -s Eq. (6) can be considered on the projective space Pd-1. The Lyapunov exponents of all system states x e Rd \ {0} can be analyzed together if the perturbation affects all states. This is expressed in the condition

dim LA j (Xo, h, , (X1, 0,0),..., (Xr, 0,0) j (9, s, t)= dim M + d (7)

for all (9, s, t)e M x Sd-1 x R. Another approach to condition (7) is as follows: Let I be the ideal in LA{Xo + h,X1,...,Xr} generated by {X1,...,Xr}. Then, by [5], Condition (7) is equivalent to dimI(9, s) = dimM + d - 1. This condition, which is needed for the analysis of moment Lyapunov exponents, is slightly stronger than Condition 7 in [6], but it is generally satisfied for systems that appear in applications, compare, e.g., [5] or [10].

Theorem 2.1. Consider the stochastic system (1) under the conditions (3) and (7). Then

1. the moment Lyapunov exponents exist as a limit and they are independent of x e Rd \ {0}, i.e., g(p) = g(p, x) = limt!11 log E|u(t, x, x)|p for all p e R,

2. the trajectory-wise Lyapunov exponents are a.s. constant and independent of x e Rd \ {0}, i.e., k = k(x, x) = limt^i 11 log |u(t, x, x)|.

The proof of Theorem 2.1 is given in [7], Theorem 1 for the first part, and in [10], Theorem 4.1 for the second part upon noticing that Conditions (3) and (7) together imply Conditions (A) and (C) in [10]. With the results from Theorem 2.1 it was shown by Arnold in [7] that the a.s. Lyapunov exponent is the derivative of the moment Lyapunov exponent function at 0:

Corollary 2.2. Consider the stochastic system (1) under the conditions (3) and (7). Then the function g(p) is analytic on R, convex, and satisfies g(0) = 0 and g'(0) = k.

Remark 2.3. The information contained in Corollary 2.2 shows that

1. If the a.s. Lyapunov exponent of the system (1) is negative, i.e., if the system is almost surely exponentially stable, then moments for small p > 0 are also exponentially stable. And vice versa, if moments for small p > 0 are exponentially stable, then the system is a.s. exponentially stable.

2. The moment exponent function g(p) has at most two zeros. Specifically, if the system is a.s. exponentially stable, i.e., if k < 0, then g(p) has at most one zero besides g(0) = 0, and this occurs for p > 0. Assume that such a second zero exists at a > 0, then all pth moments are stable for 0 < p < a, and unstable for a < p < 1.

In light of Remark 2.3 the key question regarding the exponential stability of the moments of (1) is the existence of a second zero of g(p), i.e., the existence of a > 0 with g(a) = 0. Such a point a does not always exist since g(p) = kp is possible, compare [9], Theorem 2.3, Case 2.1.2(a). A detailed analysis of this question involves the function y : R ! R, given by

y(0)=k,y{p)=S-f for p - 0. (8)

According to Corollary 2.2 this function is analytic on R and increasing. We define

y+ = limy(p) (9)

and note that either y+ = k, in which case g(p) = kp, or k = y(0) < y+, in which case y(p) is strictly increasing. Note that there exists a > 0 with g(a) = 0 if and only iff y+ > 0. For the following results see [9].

Remark 2.4. As it turns out, the quantity y+ does not depend on the statistics of the noise process n(t, m), only on its range U c Rm. More precisely, y+ is given by the exponential growth rate of the spectral radius of the (deterministic) semigroup given by (1), see also [13], Chapter 7. Similarly we have

1. y+ = supA(f)stationary k(A(f)), where the supremum is taken over all stationary processes f with values in U c Rm;

2. y+ = limj log\\u(t,x, m)kM, where || ■ ||1 denotes the essential supremum.

In other words, if we have no information about the perturbation n(t, m) except for its range U c Rm, then y+ < 0 assures exponential stability for all stochastic perturbations with values in U. This property is sometimes called 'universal stability'.

Remark 2.5. If there exists a second zero g(a) = 0 for a > 0, then more can be said about the behavior of individual trajectories of the system (1):

1. If we have k < 0 < y+, then there exists a > 0 with g(a) = 0. In this case there is a constant c p 1 such that for all R > 0 and all x 2 Rd with 0 < |x| < R it holds that

1 (x)°6 p{stPp|u(t,x,m)|p R 6 <?)

2. If we have k 6 y+ < 0, then there is a constant c p 1 such that with P-probability 1 it holds that

sup | u(t, x, m) | 6 c|x |.

In other words, if g(p) does not have a second zero, then the system response is almost surely bounded, while in case there exists a second zero at a > 0 then the probability that the system response will exceed a safety threshold R > 0 behaves like (| x| /R)a, which goes to zero for a !i. For these reasons, the point a > 0 is called the 'stability index' of a linear stochastic system in [8]. In [8,14] the authors investigate the 'robustness' of the stability index under small nonlinear model misspecifications.

3. Stability-based performance indices for stochastic systems

3.1. Performance indices

In this section we discuss performance indices for systems under stochastic perturbations. The goal is to identify system parameters that allow for optimal (exponential) stability behavior of a system at a stable operating point. Since exponential stability can be inferred from the system linearization, we consider systems as in (1), and we assume that the perturbation can be modeled by a function of a Markov-diffusion process as in (2). More specifically, our goal is to guarantee stability of the system under the largest possible perturbation range.

The size of the random perturbation is described in the following way: We consider the noise range U c Rm to be convex, compact with 0 e int U, the interior of U. Introducing the size parameter q p 0 we consider Uq := q ■ U together with the maps fq : M ! Uq,fq(9) = q ■ f (9). In this way we obtain a family {nq(t, x), q p 0} of system perturbations with corresponding dynamics (1)q. For q = 0 this model corresponds to the unperturbed system. To be precise, we analyze the family of systems

xq(t) = A(nq(t, x), b)xq(t), x e Rd, nq(t, x)e Uq c Rm, q p 0 (10)

where b e B c Rk is a vector of parameters that are to be tuned in such a way that the system (10) is stable for q p 0 as large as possible.

The almost sure stability radius

r = inf{q p 0, k(q) > 0} (11)

was introduced in [12] and analyzed in detail in [3]. Here k(q) denotes the a.s. Lyapunov exponent of (10). In a similar way, one can define the pth moment stability radius as

r(p) = inf{q p 0, gq(p) > 0} for each p e R (12)

with gq(p) as defined in (5) for the system (10). This stability radius provides an appropriate performance index if emphasis is placed on stability of specific moments of a system, such as the mean (p = 1) or mean square stability (p = 2). In both cases the design problem can be written as

maxr(b) or maxr(p, b) for a given p > 0.

beB beB

Section 2.2 points at other performance indices that can be useful for the evaluation of stability: Specifically, the second zero a(q) > 0 of the moment Lyapunov exponent function gq(p) of (10) not only describes the moments that are exponentially stable (see Remark 2.3), but also the boundedness behavior of individual trajectories (see Remark 2.5). For a given stochastic perturbation with given range Uq the design problem reads in this case

max a(q, b) for a given q > 0. (13)

Following [8] we call (13) the stability index problem. Note that under our conditions k(q) = 0 is equivalent to g(p, q) > 0 for all p > 0 which in turn is equivalent to the stability index a being zero.

Finally, Remark 2.4 shows that the analysis of the limit c+ of the function c(p) contains important information about the stability of the randomly perturbed systems (10): this number describes the worst case exponential stability behavior for any process with values in Uq. This idea leads to the deterministic stability radius

rdet = inf{q p 0, y+ (q) > 0}.

This radius turns out to be described by the maximal Bohl exponent, or the maximum of the Morse spectrum of a deterministic perturbation (or control) system associated with (10), see [13], Chapter 7 for a detailed discussion of these concepts.

In the following sections, we will analyze the moment stability radius and the stability index problem together with their relationships to the a.s. stability radius.

3.2. Computation of indices

While the computation of a.s. Lyapunov exponents has attracted great interest (see, e.g., [6,11,16], or [17] and the references therein), the computation of moment Lyapunov exponents seems relatively unexplored. One way is to write g(p) as the maximal eigenvalue of a certain second order partial differential operator (see, e.g., [7] for the white noise, and [9] for the general case). This idea has been followed, e.g., in [21], in some examples of Chapter 9 in [20], or in [18], but generalizations to high dimensional systems appear more than cumbersome.

The other approach is to follow the definition (5), i.e., simulate trajectories of the system (1), compute the moments and their exponential growth rate, see e.g., [19] or [20], Chapter 9.2 for an idea in this direction. Our experiences from [6] suggest to simulate solutions directly from the linear differential equation, using renormalization at regular time intervals since the trajectories grow or decay exponentially. This leads to the following approach:

We fix a time interval [0, T], T e N, and a step size h = 1 > 0, k e N, for the simulation of the background process gt in (2), resulting in b time series gt (i), i = 1,..., b of length Tk. We pick a initial conditions xj, = sj, e Sd-1, j = 1,..., a. For each initial condition the linear system (1) is solved on the time interval [0,1], resulting in b time series xjn)(i) with n = 0,..., k. We define x1(i) := x(k) (i) and use sj1(i) := x1(i)/||x1(i)|| e Sd-1 to be the starting point for the time interval [1,2]. Denote by x2(i) the

approximated solution at time t = 2; and renormalize as sj2(i) := xj2(i)/||xj2(i)|| e Sd-1 to obtain the starting point for the time interval [2,3]. We continue in this way over the time interval [0, T] to obtain a ■ b approximate (and renormalized) solutions of (1): Note that for t = 1,..., T the starting points for each time interval [t - 1, t] are the st-1 (i) e Sd-1, and the final points xit(i) are used for the computation of the moment Lyapunov exponents as follows:

Computing the expectation on each time interval [t - 1, t] for t = 1,..., T and j = 1,..., a we obtain

E | u(T, x0, m) |p =

E | u(T, x0, m) |p E | x0 |p

E | u(t, x0, m) |p

t=1 | E |u(t - 1, j m)|p | t-i | e |s»-1 (i) |p |

E | x (i) |p

and hence

1log E|U(T, j m) |p = 1 |>g | E | U U 't 0, m)n ^ 1 Ex (i) |p.

We now note that E |x^t(i) |p = ?= | xt (i) |p by averaging over the realizations of the background process, and finally averaging over all a initial values we obtain

g(p, x) = lim1 log E |u(T, x, m) |p :

1 1 a T /1 b

^ te E lo^1L K (0 |

j=11=1

The numerical approximation of moment Lyapunov exponents given by Formula (16) uses only numerical solutions of (1) on time intervals of length t - (t - 1) = 1 and hence avoids solutions growing or decreasing exponentially over a long period of time. If time intervals of length 1 are too long to give reliable numerics for specific systems, this approach can easily be adapted to smaller intervals. Of course, burn-in intervals and choice of initial values have to be considered carefully, see [6] for a discussion of these issues for a.s. Lyapunov exponents; these considerations apply as well to the computation of moment exponents.

The performance indices introduced in Section 3.1 depend on moment Lyapunov exponents, and they require optimization with respect to the size p of the perturbation range and/or with respect to the tunable parameters b 2 B of the system. Given the general setup we have provided in the previous sections, we do not expect analytical results on these optimization problems. Therefore, optimization is performed numerically over a grid in the parameter spaces.

4. Examples

4.1. Three-dimensional linear oscillator

Consider the linear oscillator x = A(nt)x in dimension 3 given by / 0 1 0 \

x = A(n )x =

V-c(1 + n(t))

x with x e R3.

For the computations we have used the values a = 1, bnom = 2, and c = 1. The stochastic perturbation is np(t) = p ■ sin(g(t)), where g(t, m) is an Ornstein-Uhlenbeck process as in [6] and p p 0 is the size of perturbation.

For the following discussion we consider b to be the tunable parameter with values in B = [0.8bnom, 1.2bnom], and the size of perturbation satisfies p 2 [0,1.4].

Fig. 1 shows the moment Lyapunov exponent curves for b = 0.8bnom and p 2 [0,1.4].

Moment p

Fig. 1. Moment Lyapunov exponents of the linear oscillator, b = 0.8bnom.

Fig. 1 allows to determine the pth moment stability radius as in (12): we have, e.g., r(4; b = 0.8bnom) = 1.2, or r(5; b = 0.8bnom) = 0.95. Note that r(p; b = 0.8bnom) > 1.4 for all p 2 (0,3.5). Furthermore, Fig. 1 shows the stability index proposed in (13) which corresponds to the moment p, for a fixed q, where system will not remain stable. In this case if the size of the perturbation is q = 1.4, then all moments p p 3.5 of the system will be unstable.

To compare the a.s. Lyapunov exponents with the derivative of the moment Lyapunov exponents at p = 0 (compare Corollary 2.2), we also computed the a.s. exponent k(q) using the methodology described in [6]. The results are listed in Tables 1 and 2, showing a very good agreement between the two methods.

Fig. 2 shows the case where b = 1.2bnom and the size of perturbation is q 2 [0; 1.4]. Comparing the results to Fig. 1 we see that increased damping, indicated by increasing b, leads to larger moment stability radii. For example, for b = 1.2bnom we have r(5;b = 1.2bnom) = 1.3.

In the case where b = 1 :2bnom, the stability index defined in (13) corresponds to p = 4:8 when q = 1:4.

4.2. One machine - infinite bus power system

In this section we present the application of the proposed methodology to the example of a power system consisting of one machine connected to a infinite bus; see [15] for more information regarding this model.

Model of linearized system with noise

The state vector for the linearized system x = Ax, is

x = [Am, AS, DWfd, v 1, v2, v3]

Table 1

Almost sure Lyapunov exponents for q = 0.8.

q = 0.8 b = 0.8bnom b = 0.9bnom b = 1.0bnom b = 1.1bnom b = 1.2bnom

k(q) —0.15831 —0.19362 —0.22268 —0.25087 —0.27052

Cq(0) -0.15758 —0.19363 —0.22349 —0.25225 —0.27175

Table 2

Almost sure Lyapunov exponents for q = 1.4.

q=1:4 b = 0.8bnom b = 0.9bnom b = 1.0bnom b = 1.1bnom b = 1.2bnom

k(q) —0.13345 —0.16359 —0.19494 —0.22503 —0.24130

Cq(0) —0.13071 —0.16141 —0.19348 —0.22415 —0.24039

,3i-1-1-1-1-1-1-1-1-1-1

01 23456789 10

Moment p

Fig. 2. Moment Lyapunov exponents of the linear oscillator, b = 1.2b,

where vi, v2, v3 are variables of the power system stabilizer, using expressions defined in [15]. The matrix A has the structure

fan ai2 ai3 0 0 0

a2i 0 0 0 0 0

0 a32 a33 a34 0 a36

0 a42 a43 a44 0 0

asi as2 as3 0 ass 0

\aei a62 a63 0 a6s a66

The model for the field circuit is K 3

DWfd = 1

■(DEjd - K4AS)

and the excitation system is given by

DEfd = -KaA v i,

where v1 is the output of voltage transductor. The perturbation has been introduced as an error in the reference signal. This situation is described by changing the element a34 in the matrix A to

AEfd = -Ka(1 + nt)Avi.

To be precise, we consider the (linearized) one machine - infinite bus system x = A(nt)x with system matrix

0 -o.ii -0.12 0 0 0

377 0 0 0 0 0

0 -0.19 -0.42 a34 0 27.4

0 -7.3 20.8 -50 0 0

0 -1 -1.1 0 -0.71 0

0 -4.8 -5.4 0 26.9 30.3

A(nt) =

where a34 = a34 ■ (1 + nt) is a stochastic perturbation in the excitation component of the system. The stochastic perturbation is nt = p ■ sin(gt) with gt an Ornstein-Uhlenbeck process as in [6] and p is the size of perturbation.

In the similar case of Example 4.1, the key parameter in this system is the gain of the PSS, KPSS, whose nominal value Kpss„om was chosen as in [15]. In order to compute Performance Indices, we used gain values Kw := w ■ Kpss„om, with w 2 [0.8,1.2]. For the range of the random perturbation we used p e [0,0.5], with a step size of 0.1. Fig. 3 shows the moment Lyapunov exponents curves for Kw = 0.8 ■ KPSSnom and p e [0,0.5].

Fig. 3 allows us to determine the pth moment stability radius as in (12): we have, e.g., r(2; Kw = 0.8 ■ Kpss„om) = 0.1, or r(1.5; Kw = 0.8 ■ Kpss„om) = 0.19. Similar to the findings of Example 4.1, Fig. 3 shows the stability index proposed in (13) which corresponds to the moment p, for a fixed p, where system will not remain stable. In this case if the size of the perturbation is p = 0.5, then all moments p p 1.2 of the system will be unstable.

Moment p

Fig. 3. Moment Lyapunov exponents of the one machine system, Kw = 0.8KPSSnom

Fig. 4. Moment Lyapunov exponents of the one machine system, Kw — 1.2Kpssnom

Fig. 4 shows the moment Lyapunov exponents curves for Kw = 1.2 ■ KPSSnom and p 2 [0; 0.5]. Comparing the two cases Kw = 0.8 ■ KPSSnom in Fig. 3 and Kw = 1.2 ■ KPSSnom in Fig. 4 we see that the second zeros of the moment Lyapunov exponents curves are not monotone in Kw: we have a(0.5,0.8 ■ KpSS„om) > a(0.5,1.2 ■ KpSS„om), while a(p, 0.8 ■ KpSS„om) > a(p, 1.2 ■ KpSS„om) for all p 2 [0.1,0.4]. This indicates that optimal parameter tuning relative to exponential moment stability cannot simply be achieved by increasing PSS gains.

5. Conclusions

This paper proposes several performance indices for the stability of operating points in dynamical systems affected by sustained random perturbations. These indices are based on moment Lyapunov exponents and they complement the almost sure Lyapunov exponent and stability radius analyzed in [6,2]. The two sets of indices are related by the fact that the almost sure Lyapunov exponent of a system is the derivative at p — 0 of the pth moment Lyapunov exponent function g(p). This means, in particular, that for small moments p > 0 the pth moment Lyapunov exponents contain the 'same' information as the a.s. exponent, while for large moments p > 0 the pth moment Lyapunov exponents contain the 'same' information as the maximal deterministic (robust) exponent c+, compare Remark 2.4. Hence for design purposes stability indices based on moment Lyapunov exponents can be used to strike a balance between almost sure behavior based on specific random perturbations, and behavior based on the range of the perturbation. Design issues surrounding moment stability indices are the topic of a forthcoming paper.

Note that while we always have for the moment Lyapunov function g(p) that g(0) — 0, for a stable operating point the second zero of g(p), i.e., the point a > 0 with g(a) — 0 determines the moment stability behavior. In realistic systems, such as the one machine - infinite bus power system, this second zero a may not depend in a monotone fashion on the size of the random perturbation, and on the amount of damping in the system. This indicates that optimal parameter tuning relative to exponential moment stability cannot simply be achieved by increasing system damping, such as PSS gains.

A key question then is to what extend system design that uses indices based on almost sure or moment Lyapunov exponents, depends on the specific statistics of the system noise n(t, rn). Of course, if one wants to immunize a system against all specific noise statistics, one will have to use the deterministic (robust) exponent c+: this index depends only on the size of the perturbation, not on its statistics, see the comments above. The robustness of the design indices presented in this paper relative to noise statistics is the topic of ongoing research.

Acknowledgment

This research was financed by Project Fondecyt 11130169 (Comision Nacional de Investigacion Cientifica y Tecnologica Chile) and by University de Santiago de Chile, Project DICYT N 91.

References

[1] J.J. Sanchez-Gasca, V. Vittal, M.J. Gibbard, A.R. Messina, D.J. Vowles, S. Liu, U.D. Annakkage, Inclusion of Higher Order Terms for Small-Signal(Modal) Analysis: Committee Report—Task Force on Assessing the Need to Include Higher Order Terms for Small-Signal (Modal) Analysis, IEEE Trans. Power Syst. 20 (4) (2005) 1886-1904.

[2] H. Verdejo, L. Vargas, W. Kliemann, Improving PSS Performance under Sustained Random Perturbations, IET Gen. Transm. Distrib. 6 (9) (2012) 853862.

[3] H. Verdejo, L. Vargas, W. Kliemann, Stability radii via Lyapunov exponents for large stochastic systems, Procedia IUTAM (2013) 188-193.

[4] L. Arnold, E. Oeljeklaus, E. Pardoux, Almost sure and moment stability for linear Ito equations, Springer Lecture Notes in Mathematics No. 1186 (1986) 129-159.

[5] K. Ichihara, H. Kunita, A classification of the second order degenerate elliptic operators and its probabilistic characterization, Z. Wahrsch. Verw. Gebiete 30 (1974) 235-254.

[6] H. Verdejo, L. Vargas, W. Kliemann, Stability of linear stochastic systems via Lyapunov exponents and applications to power systems, Appl. Math. Comput. 218 (2012) 11021-11032.

[7] L. Arnold, A formula connecting sample and moment stability of linear stochastic systems, SIAM J. Appl. Math. 44 (1984) 793-802.

[8] L. Arnold, R.Z. Khasminskii, Stability index for nonlinear stochastic differential equations, Proc. Symp. Pure Math. 57 (1995) 543-551.

[9] L. Arnold, W. Kliemann, Large deviations of linear stochastic differential equations, in: Springer Lecture Notes in Control and Information Sciences No. 96(1987), 117-151.

[10] L. Arnold, W. Kliemann, E. Oeljeklaus, Lyapunov exponents of linear stochastic systems, Springer Lect. Notes Math. 1186 (1986) 85-125.

[11] F. Carbonell, R. Biscay, J.C. Jimenez, QR-based methods for computing Lyapunov exponents of stochastic differential equations, Int. J. Numer. Anal. Model. B 1 (2010) 147-171.

[12] F. Colonius, W. Kliemann, Stability radii and Lyapunov exponents, in: D. Hinrichsen, B. Martensson (Eds.), Control of Uncertain Systems, Birkhäuser, 1990, pp. 19-56.

[13] F. Colonius, W. Kliemann, The Dynamics of Control, Birkhäuser, 2000.

[14] R.Z. Khasminskii, On robustness of some concepts in stability of stochastic differential equations, Fields Inst. Commun. 9 (1996) 131-137.

[15] P. Kundur, Power System Stability and Control, McGraw-Hill, reprinted edition, 1994.

[16] D. Talay, Approximation of Upper Lyapunov Exponents of Bilinear Stochastic Differential Systems, SIAM J. Numer. Anal. 28 (1991) 1141-1164.

[17] D. Talay, Simulation and numerical analysis of stochastic differential systems: a review, in: P. Kree, W. Wedig (Eds.), Probabilistic Methods in Applied Physics, Lecture Notes in Physics 451, Springer-Verlag, 1995, pp. 54-96.

[18] I.G. Vladimirov, The monomer-dimer problem and moment Lyapunov exponents of homogeneous Gaussian random fields, Discr. Cont. Dyn. Syst. 18 (2013) 575-600.

[19] W.-C. Xie, Monte Carlo simulation of moment Lyapunov exponents, ASME J. Appl. Mech. 72 (2005) 269-275.

[20] W.-C. Xie, Dynamic Stability of Structures, Cambridge University Press, 2006.

[21] W.-C. Xie, R.M.C. So, Numerical determination of moment Lyapunov exponents of two-dimensional systems, ASME J. Appl. Mech. 73 (2006) 120-127.