Stability Assessment of Stochastic Differential-Algebraic Systems via Lyapunov Exponents with an Application to Power Systems

In this paper we discuss Stochastic Differential-Algebraic Equations (SDAEs) and the asymptotic stability assessment for such systems via Lyapunov exponents (LEs). We focus on index-one SDAEs and their reformulation as ordinary stochastic differential equation (SDE). Via ergodic theory it is then feasible to analyze the LEs via the random dynamical system generated by the underlying SDE. Once the existence of well-defined LEs is guaranteed, we proceed to the use of numerical simulation techniques to determine the LEs numerically. Discrete and continuous QR decomposition-based numerical methods are implemented to compute the fundamental solution matrix and to use it in the computation of the LEs. Important computational features of both methods are illustrated via numerical tests. Finally, the methods are applied to two applications from power systems engineering, including the single-machine infinite-bus (SMIB) power system model.

Authors

• 1 publication
• 1 publication
• 2 publications
• 15 publications
03/29/2020

Solving the inverse problem for an ordinary differential equation using conjugation

We consider the following inverse problem for an ordinary differential e...
11/08/2019

Coupled Systems of Differential-Algebraic and Kinetic Equations with Application to the Mathematical Modelling of Muscle Tissue

We consider a coupled system composed of a differential-algebraic equati...
09/13/2021

Computation of the nearest structured matrix triplet with common null space

We study computational methods for computing the distance to singularity...
02/14/2021

Reflecting stochastic dynamics of active-passive crowds in a queueing theory model

Stochastic differential equation (SDE) models have been extensively used...
07/05/2021

An analytical solution of the isentropic vortex problem in the special relativistic magnetohydrodynamics

The isentropic vortex problem is frequently solved to test the accuracy ...
04/28/2019

Stability conditions of an ODE arising in human motion and its numerical simulation

This paper discusses the stability of an equilibrium point of an ordinar...
10/26/2014

A Novel Statistical Method Based on Dynamic Models for Classification

Realizations of stochastic process are often observed temporal data or f...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Modeling the dynamic behavior of systems employing differential-algebraic equations is a mathematical representation paradigm widely used in many areas of science and engineering. On the other hand, the dynamics of systems perturbed by stochastic processes have been adequately modeled by stochastic differential equations (SDEs). The need of a generalized concept which covers both DAEs and SDEs, and allows the modeling and analysis of constrained systems subjected to stochastic disturbances, has led to the formulation of stochastic differential-algebraic equations (SDAEs). While there has been a broad and fruitful development both in the field of DAEs and SDEs (see, e.g. [6, 7, 23] and [20, 30, 32], respectively), studies of the concepts and the numerical treatment of SDAEs are rather limited, see e.g. [10, 37, 42]. The main reason for this is that the proper treatment of algebraic constraints in SDAEs faces many difficulties, except in the case that all the constraints are explicitly given and can be resolved during the numerical integration process, which is the case that we will discuss.

As tool for the stability analysis we use Lyapunov exponents (LEs) introduced as characteristic exponents in [28]. The theory of LEs experienced a crucial development with the work [33] which, via the Multiplicative Ergodic Theorem (MET), ensures the regularity and existence of the LEs belonging to a linear cocycle over a metric dynamical system, see [3]

for a detailed presentation of the theory. We review and extend the main concepts of this approach for asymptotic stability assessment of differential-algebraic equations driven by Gaussian white noise and apply the technique in the setting of power systems.

Following the ideas from [10, 11, 24, 37, 42], properties such as the existence and uniqueness of solutions are reviewed. Analogously to the DAE case, we define the class of strangeness-free (index-one) SDAEs and show that SDAE systems with this structure can be reduced to a classical SDE system that describes the dynamics of the original SDAE for which the MET can be applied to define the LEs of the system.

Once we have extended the theoretical framework that guaratees the existence of well-defined LEs, we study numerical methods, based on the factorization of the fundamental solution matrix, that allow the numerical computation of spectral values asociated to the Lyapunov spectra. The first technique requires computing the fundamental solution matrix and forming an orthogonal factorization, while the second one involves performing a continuous decomposition of the fundamental solution matrix. Both techniques have been extensively studied in deterministic ODE and DAE systems, see [4, 5, 15, 14, 26, 27]. This paper follows the ideas exposed in [9], where these methods were extended to the stochastic case.

Finally, these concepts and computational techniques are used to assess the asymptotic stability of power systems affected by stochastic fluctuations. We illustrate the results with elementary test cases such as a single-machine infinite-bus system.

The paper is organized as follows: In Section 2 we recall the main theory of strangeness-free SDAEs, their relation with the SDEs, and the existence of LEs generated by such SDEs. Section 3 presents the discrete and continuous -based decomposition methods and their evaluation. Interesting study-cases in the power systems area are presented in Section 4. We finish with some conclusions in Section 5.

2 Review of the theory

2.1 Stochastic Differential-Algebraic Equations

Consider a system of quasi-linear stochastic differential-algebraic equations (SDAEs) of the form

 Edxt=f0(xt)dt+m∑j=1fj(xt)dwjt,t∈I:=[t0,tf], (1)

with a singular matrix of rank . The function (for some ) is known as drift, and are the diffusions, here is an open set. Furthermore, (for ) form an -dimensional Wiener process

defined on the complete probability space

with a filtration . Each -th Wiener process is understood as a process with independent increments such that

, i.e. is a Gaussian random variable for all

, such that

 w0=0,E[wt−ws]=0,E[wt−ws]2=t−s.

Since the Wiener process is nowhere differentiable, it is more convenient to represent equation (1) in its integral form as

 Ext=Ext0+∫tt0f0(xs)ds+m∑j=1∫tt0fj(xs)dwjs,t∈I. (2)

Here, the first integral is a stochastic Riemann-Stieltjes integral, and the second one is a stochastic Itô-type integral, see e.g. [32].

We assume consistent initial values independent of the Wiener processes

and with finite second moments

[32]. A solution of (2) is an

-dimensional vector-valued Markovian stochastic process depending on

and (the parameter is commonly omitted in the notation of ). Such a solution can be defined as strong solution if it fulfills the following conditions, see e.g. [11, 42].

• is adapted to the filtration ,

• almost sure (a.s.), for all ,

• a.s., for all , and ,

• (2) holds for every a.s.

Because of the presence of the algebraic equations associated with the kernel of , the solution components associated with these equations would be directly affected by white noise and not integrated. In order to avoid this, a reasonable restriction is to ensure that the noise sources do not appear in the algebraic constraints. According to [37, 42], this assumption can be accomplished in SDAE systems whose deterministic part

 E˙xt=f0(xt),t∈I, (3)

is a DAE with tractability index one [25, 42], in which the constraints are regularly and globally uniquely solvable for parts of the solution vector. We slightly modify this assumption and consider SDAE systems whose deterministic part (3) is a regular strangeness-free DAE [23] i.e. it has differentiation index one. A system with these characteristics can be transformed into a semi-explicit form by means of an appropriate kinematic equivalence transformation [26, 6]

, i.e., there exist pointwise orthogonal matrix functions

and such that, pre-multiplying (1) by , and changing the variables according to the transformation one obtains a system in semi-explicit form

 d^xDt =^fD0(^xDt,^xAt)dt+m∑j=1^fDj(^xDt,^xAt)dwjt, (4a) 0 =^fA0(^xDt,^xAt)dt+m∑j=1^fAj(^xDt,^xAt)dwjt, (4b)

where and is a separation of the transformed state into differential and algebraic variables, respectively, that is performed in such a way that the Jacobian of the function with respect to the algebraic variables is nonsingular, see [23] for details of the construction. The condition that the noise sources do not appear in the constraints, implies that , so that the algebraic equations in (4b) can be solved as and inserted in the dynamic equations (4a) yielding an ordinary SDE

 d^xDt =^fD0(^xDt,FA(^xDt))dt+m∑j=1^fDj(^xDt,FA(^xDt))dwjt. (5)

This equation is called underlying SDE of the strangeness-free SDAE. It acts in the lower-dimensional subspace , with (where denotes the number of algebraic equations). The SDE system (5) preserves the inherent dynamics of a strangeness-free SDAE system [25]. Note that in this way, the algebraic equations have been removed from the system, but whenever a numerical method is used for the numerical integration, then one has to make sure that the algebraic equations are properly solved at each time step, so that the back-transformation to the original state variables can be performed.

2.2 Random Dynamical Systems generated by SDEs

In the previous section we have discussed the reduction of an autonomous strangeness-free SDAE to its underlying SDE, which preserves the dynamic characteristics of the original system. Using the back-transformation the definitions and properties attributed to the underlying SDE, and the analysis performed on it can be extended to the original SDAE. For simplicity we use the following representation, where the drift and diffusion terms are combined into one term.

 dxt=f0(xt)dt+m∑j=1fj(xt)dwjt=m∑j=0fj(xt)dwjt,t∈I, (6)

where , and for some and . Here is the Banach space of vector fields on with linear growth and bounded derivatives up to order and the -th derivative is -Hölder continuous. In addition, we assume that the differential operator is strong hypoelliptic in the sense that the Lie algebra generated by the vector fields (with ) has dimension for all  [3]. Once again, (for ) is a -dimensional Wiener process, this time with the convention .

For a given initial value, the solution process generates a Markovian stochastic process, and the SDE (6) generates a random dynamical system (RDS) which is an object consisting of a metric dynamical system (MDS) for modeling the random perturbations, and a cocycle over this system. The ergodic MDS is denoted by with the filtration , and defined by the Wiener shift

 θtω(⋅)=ω(t+⋅)−ω(t),t∈I,

which means that a shift transformation given by is measure-preserving and ergodic [8].

Together with the SDE (6), we define the variational system

 dvt=m∑j=0Jj(xt)vt dwjt,t∈I,withJj(xt):=(∂fj∂x)∣∣∣x(t;t0), (7)

obtained after linearizing (6) along a solution. If we denote by the Jacobian of at , then is the unique solution of the variational equation (7), satisfying

 Φ(t,ω,x)=Id+m∑j=0∫t0Jj(φ(s)x)Φ(s,x)dwjs,t∈I, (8)

where

denotes the identity matrix of size

. Therefore, it is a matrix cocycle over . The system (6)-(7) uniquely generates a RDS over . Moreover, the determinant of satisfies the Liouville equation

 detΦ=exp(m∑j=0∫t0tr(Jj(φ(s)x)dwjs),t∈I, (9)

being then a scalar cocycle over (see Theorems 2.3.32 and 2.3.39-40 in [3] for details).

2.3 Lyapunov Exponents of Ergodic RDSs

An important result of the theory of RDSs is the so-called Multiplicative Ergodic Theorem (MET) developed in [33]. This concept allows the definition of LEs for linear cocycles over a ergodic MDS. First, the MET assumes for the linear cocycle that the integrability condition

 log+∥Φ(t,ω,x)∥∈L1,

is satisfied, where denotes the positive part of . This guarantees that the variational equation (7) associated with (6) is well-posed. Additionally, let be an ergodic invariant measure with respect to the cocycle  [3]. Then, the MET assures the existence of an invariant set of full -measure, such that for each there is a measurable decomposition

 Rd=L1(ω)⊕⋯⊕Lp(ω),

of into random linear subspaces , which are invariant under . Here , where denotes the dimension of the subspace (with ), and . This splitting is dynamically characterized by real numbers which quantify the exponential growth rate of the subspaces. These are called Lyapunov exponents, and are defined by

 λi:=limt→∞1tlog∥Φ(t,ω,x)∥\rm wheneverx∈Li(ω)∖{0}.

According to [2, pag. 118], the LEs are independent of and thus they are universal constants of the cocycle generated by (7) under the ergodic invariant probability measure . Finally, the following identity holds for when the system is Lyapunov regular [34, 9]

 d∑i=1λi=limt→∞1tlog|detΦ(t,ω,x)|. (10)

In practice, it is hard (if not impossible) to verify Lyapunov regularity for a particular system [3]. One of the key statements of the MET is that linear RDS (whether these are constant, periodic, quasi-periodic, or almost-periodic) are a.s. Lyapunov regular.

The concept of LEs plays an important role in the asymptotic stability assessment of dynamical systems subjected to stochastic disturbances. Under appropriate regularity assumptions, the negativity of all LEs of the system of variational equations implies the exponential asymptotic stability of both the linear SDE and the original nonlinear SDE system.

3 QR Methods for computing LEs

In this section we derive the numerical techniques to compute the finite-time approximation of the LEs. In the same way as [9], this paper proposes an adaptation of the ideas from [12, 15, 14] to the stochastic case (linear RDSs). The methods take advantage on the existence of a Lyapunov transformation of the linear RDS to an upper-triangular structure, and the feasibility to retrieve a numerical approximation of the LEs from that form. The transformation is performed through an orthogonal change of variables. The approach is made under the assumption of Lyapunov regularity of the system. In order to explain the methods, let us consider again the SDE as an initial value problem of the form

 dxt=m∑j=0fj(xt)dwjt,t∈I,xt0=xt0, (11)

where are sufficiently smooth functions. The corresponding variational equation of (11) along with the solutions , turned into a matrix initial value problem, is given by

 dVt=m∑j=0Jj(xt)Vtdwjt,V0=Id, (12)

with the identity matrix as initial value, where are the Jacobians of the vector functions , and is the fundamental solution matrix, whose columns are linearly independent solutions of the variational equation. A key theoretical tool for determining the LEs is the computation of the continuous factorization of ,

 Vt=QtRt,

where is orthogonal, i.e., , and is upper triangular with positive diagonal elements for . Applying the MET theory presented in Subsection 2.3, and taking into account the norm-preserving property of the orthogonal matrix function , we have

 λi=limt→∞1tlog∥Vtpi∥=limt→∞1tlog∥Rtpi∥, (13)

where is an orthonormal basis associated with the splitting of . Lyapunov regular systems preserve their regularity under kinematic similarity transformations. Then, considering the regularity condition (10), the Liouville equation (9), and performing some algebraic manipulations (see details in [9, pag. 150]), the LEs are given by

 λi=limt→∞1tlog|Riit|a.s.,fori=1,…,d. (14)

The methods require to perform the decomposition of for a long enough time, so that the have started to converge. Depending on whether the decomposition is performed after or before integrating numerically the variational equation, the method is called discrete or continuous method.

3.1 Discrete QR Method

The discrete method is a very popular method for computing LEs in ODEs and DAEs. In this approach, the fundamental solution matrix and its triangular factor are indirectly computed by a reorthogonalized integration of the variational equation (12) through an appropriate decomposition. Thus, given grid points , we can write in terms of the state-transition matrices as

 Vtℓ=Z(tℓ,tℓ−1)Z(tℓ−1,tℓ−2)⋯Z(t2,t1)Z(t1,t0)Vt0. (15)

At , we perform a standard matrix decomposition

 Vt0=Qt0Rt0,

and for , we determine as the numerical solution (via numerical integration) of the matrix initial value problem

 dZ(tℓ,tℓ−1)=m∑j=0Jj(xt)Z(tℓ,tℓ−1)dwjt,Z(tℓ,tℓ−1)=Qtℓ−1,tℓ−1≤t≤tℓ, (16)

and then compute the decomposition

 Z(tℓ,tℓ−1)=QtℓR(tℓ,tℓ−1),

where has positive diagonal elements. From (15), the value of the fundamental matrix is determined via

 Vtℓ=QtℓR(tℓ,tℓ−1)R(tℓ−1,tℓ−2)⋯R(t2,t1)R(t1,t0)Rt0,

which is again a factorization with positive diagonal elements. Since this is unique, for the decomposition , we have

 Rtℓ=R(tℓ,tℓ−1)R(tℓ−1,tℓ−2)⋯R(t2,t1)R(t1,t0)Rt0=ℓ∏κ=0Rκ.

Here we denote as the triangular transition matrices with . From (14), the LEs are thus computed as

 λi=limℓ→∞1tℓlog∣∣ ∣∣ℓ∏κ=0Riiκ∣∣ ∣∣=limℓ→∞1tℓℓ∑κ=0log|Riiκ|,i=1,…,d. (17)

3.2 Continuous QR Method

The implementation of the continuous technique requires to determine a system of SDEs for the factor and the scalar equations for the logarithms of the diagonal elements of the factor elementwise. Then, once the orthogonal matrix is computed by numerical integration, the logarithms of the diagonal elements of can also be obtained.

By differentiating in the Itô sense the decomposition and using the orthogonality , we get

 dVt =(dQt)Rt+Qt(dRt), (18) 0 =(dQTt)Qt+QTt(dQt). (19)

Inserting (18) into the variational equation (12), and multiplying by from the left and by from the right, we obtain

 QTt(dQt)+(dRt)R−1t=m∑j=0QTtJj(xt)Qtdwjt. (20)

Since

is upper triangular, the skew-symmetric matrix

satisfies

 dSilt=⎧⎪ ⎪⎨⎪ ⎪⎩ ∑mj=0(QTtJj(xt)Qt)jldwjt,i>l, 0,i=l, −∑mj=0(QTtJj(xt)Qt)jldwjt,i

This results in an SDE for given by

 dQt=QtdSt=m∑j=0QtTjt(xt,Qt)dwjt, (22)

where the matrices (for ) are defined via

 (Tjt(xt,Qt))il=⎧⎪ ⎪⎨⎪ ⎪⎩ (QTtJj(xt)Qt)jl,i>l, 0,i=l, −(QTtJj(xt)Qt)jl,i

A corresponding SDE for can be obtained from (20) and (21) via

 dRt=m∑j=0(QTtJj(xt)Qt−Tit(xt,Qt))Rtdwjt, (24)

and the equation for the th diagonal element is given by

 dRiit=m∑j=0(QTtJj(xt)Qt)iiRiitdwjt,% for i=1,…,d. (25)

Since the computed LEs can be obtained from (14), we make use of the Itô Lemma to introduce the following SDE for the function from (25),

 dψit=d(logRiit)=m∑j=0(QTtJj(xt)Qt)iidwjt−12[m∑j=0(QTtJj(xt)Qt)iidwjt]2. (26)

If we assume that there are no correlations between the diffusion terms in the SDE system, then we do not have terms (for , and , with ) in the SDE (26).Also, using that , , and for , the SDE (26) is reduced to

 dψit=m∑j=0(QTtJj(xt)Qt)iidwjt−12m∑j=1[(QTtJj(xt)Qt)ii]2dt. (27)

By integrating this SDE, it is possible to obtain the LEs from

 λi=limt→∞1tψit,i=1,…,d. (28)

The difference between the discrete and the continuous method is that for the first one, the orthonormalization is performed numerically at every discrete time step, while the continuous method maintains the orthogonality via solving differential equations that encode the orthogonality continuously.

3.3 Computational Considerations

In this section, we discuss additional aspects of the computational implementation of discrete and continuous methods to calculate LEs. The application of the discrete technique mainly requires the numerical integration of the SDEs (11) and (16). This task is performed by using standard weak Euler-Maruyama and Milstein schemes, which preserve the ergodicity property (see [9, 18, 38]).

On the other hand, the numerical integration of the SDEs (11), (22) and (27) in the computational implementation of the continuous technique, must be performed in such a way that it preserves the orthogonality of the factor in each integration step. This can be achieved via projected orthogonal schemes which consist of a two-step process in which first an approximation is computed via any standard scheme, and then the result is projected into the set of orthogonal matrices [13]. Again we use the Euler-Maruyama and Milstein method as in the discrete case.

We have implemented the two -methods in Matlab. However, to obtain a unique factorization in each step, we have modified the decomposition provided by Matlab to ensure this uniqueness, by forming a diagonal matrix with , for ; and then setting and .

3.4 A Simple Numerical Example

In this subsection we illustrate the described procedures via a simple strangeness-free SDAE system in order to compare the computational efficiency and accuracy of both the discrete and continuous methods using the numerical integration schemes Euler-Maruyama and Milstein. The four numerical methods will be denoted as D-EM, D-Milstein, C-EM, C-Milstein, respectively. The computations are carried out with Matlab Version 9.7.0(R2019b) on a computer with CPU Intel Core i7 composed by 6 cores of 2.20GHz, and 16 GB of RAM.

As a simple example consider the SDAE equation

 [1000]⋅d[x1x2]=[−x2−αx1+arctan(x1)+x2]dt+[(x21+1)120]dwt, (29)

with . The nonlinear functions in both the drift and diffusion part are continuous on , with continuous and bounded derivatives, and is a one-dimensional Wiener process. The underlying SDE of (29) is

 d^xt=[−α^xt+arctan(^xt)]dt+(^x2t+1)12dwt, (30)

whose LE exists and can be explicitly represented as the following integral with respect to the solution of a stationary Fokker–Planck equation (see further details in [3])

 λ=−α+12∫R(^x2−2)(^x2+1)p(x)dx, (31)

where is the stationary density of the unique invariant probability law of . By solving numerically (31) for , we obtain the exact value of the LE associated to (30 and its original SDAE (29, which is . The accuracy of the -based methods will be assessed by comparing with this value as reference.

A large number of simulations have been carried out for stepsizes with ; to obtain computed approximations of the LE truncated at the final time , denoted by . To complete our stochastic numerical analysis of the LE, we have calculated the values of expectation

, and variance

; estimated from

independent realizations. Some results are presented in Tables (A.1) to (A.4), taking and .

Observe that the time scale in Figure 1 has been conveniently adjusted to the range , in order to show the exponential drop of the LE for the different realizations in the four methods, along the time evolution. While in Figure 2 the time scale has been adjusted to the range , to better display the convergence of the mean and variance of the LE.

Based on the analytic expression of the LE, given by equation (31), the LE can be considered as a deterministic quantity. According to the numerical results obtained from the four -based methods, the sequences of random variables reveal a trend towards null variance and convergence to the mean as tends to infinity. Such evolution can be seen in Figure 1, and more obviously in Figure 2. For all the methods, an exponential decay is illustrated in and as tends to infinite. This behavior indicates a mean square (m.s.) convergence of those sequences to a degenerate random variable, based on the implication that if is such that , for all , and , then . This means that the limit of can be interpreted as a deterministic value with probability . This enables us to state that the stochastic approximations converge in m.s. sense to a number (a degenerate random variable), which is expected to represent the LE .

In Figure 3 we compare the relative error of the accuracy of the four numerical methods for different stepsize and time interval . From this graphical representation, we observe that continuous methods obtain better results than discrete ones, as expected. We also observe that the Milstein method has, in general, better accuracy than the Euler-Maruyma scheme, since its convergence order is higher, but requies more computational time. This latter fact is evidenced in Figure 4, where a comparison of CPU time (in seconds) is shown for different values of and . Here, we observe that all the methods are affected to the same extent by incrementing the simulation interval , via a logarithmic increment, and by narrowing the stepsizes , via an exponential increment. A more pronounced difference between the methods should be observed in higher dimensional systems.

4 Application of LEs to Power Systems Stability Analysis

The concept of stability (based on Lyapunov exponents) in power systems is, in essence, the same as that for a general dynamical system. In the literature, power system stability is defined as the stability to regain an equilibrium state after being subjected to physical disturbances [21, 22]. Such equilibrium is characterized through three significant quantities during the power system operation: angles of nodal voltages, frequency, and nodal voltage magnitudes. Based on this triplet, there is an entire classification proposed by the Institute of Electrical and Electronics Engineers (IEEE) and the International Council on Large Electric Systems (CIGRE) in [22], which is illustrated in Figure 5.

The test cases presented in this paper are oriented to evaluate the angle and voltage stability of power systems subjected to small or large disturbances. Studies considering small disturbances are commonly known as small-signal stability assessment (SSSA). Here, linear stability analysis via eigenvalues has been one of the traditional analysis tools to predict the degree of stability of the power system

[21, 36]. However, eigenvalue analysis is limited to linear time-invariant systems or systems close to a stationary solution. When time-varying systems are tested, as is the case of systems subjected to stochastic disturbances, then eigenvalue analysis is no longer applicable. On the other hand, the stability analysis of power systems affected by large disturbances, known as transient stability assessment, is mainly performed with verification strategies based on time-domain integration [22, 29].

Since the concept of LEs is based on the trajectories of the dynamical systems, the method is an interesting measure of dynamic stability for power systems under stochastic disturbances in general. So, testing asymptotic stability of power systems via LEs has become an attractive approach for the two areas mentioned before, namely, the SSSA of rotor angles and voltages, by using the linearized set of SDAEs which model the system [40, 39]; and strategies for the rotor angles via transient analysis using the nonlinear SDAE system and its variational equation [19], [41]. For both cases, asymptotic stability is checked via approximations of the largest Lyapunov exponent (LLE) of the system. In particular, a negative LLE indicates that the dynamics of the system is asymptotically stable. In the next subsections test cases are presented that illustrate for strangenss-free SDAE systems the negativity of the LLE.

4.1 Modeling Power Systems through SDAEs

Under the assumption of deterministic dynamic behavior, power systems are typically modelled via a system of quasi-linear DAEs with partitioned variables, see [21],[36], of the form

 E11dxD1t =fD10(xD1t,xAt)dt, (32a) 0 =fA0(xD1t,xAt), (32b)

where is a diagonal block matrix, , , are the dynamic state variables, and are the algebraic state variables and we set . The DAE system (32) is strangeness-free (of index one).

The dynamic behavior of synchronous machines, system controllers, power converters, transmission lines, or power loads are adequately represented through such a DAE formulation. But in current real world systems, the dynamic behavior of power systems is affected by disturbances of a stochastic nature such as renewable stochastic power generation, rotor vibrations in synchronous machines, stochastic variations of loads, electromagnetic transients, or perturbations originated by the measurement errors of control devices, see [31]. Such disturbances can be modeled through Itô SDEs of the form

 dxD2t=fD20(xD1t,xD2t,xAt)dt+fD21(xD1t,xD2t,xAt)dwt. (33)

Here, is the drift, is the diffusion, are the stochastic variables, and is the Wiener process. By combining (33) and (32), and assuming that perturbs (32a) and (32b), we obtain a strangeness-free SDAE system of the form

 E1dxD1t =fD10(xD1t,xD2t,xAt)dt, (34a) dxD2t =fD20(xD1t,xD2t,xAt)dt+fD21(xD1t,xD2t,xAt)dwt, (34b) 0 =fA0(xD1t,xD2t,xAt), (34c)

or in simplified notation as

 Edxt=f0(xt)dt+f1(xt)dwt,xt0=x0, (35)

with

 E:=[E11000Id20000], xt:=⎡⎢ ⎢⎣xD1txD2txAt⎤⎥ ⎥⎦,

drift , and diffusion , where .

The study-cases presented below are formulated as the form (34). An alternative approach for including the stochastic disturbances is to implement the Wiener process directly in the underlying ODE of the system, turning them into SDEs (see [8, 16] for examples).

4.2 Modeling Stochastic Perturbations

In this subsection, we discuss on the modeling of stochastic variations via SDEs. We employ the well known mean-reverting process termed Ornstein–Uhlenbeck (OU) process [31, 17]. The SDE which defines the OU process has the form

 dηt=α(μ−ηt)dt+βdwt,ηt0=η0,t∈I, (36)

where . The OU process is a stationary autocorrelated Gaussian diffusion process distributed as . Another mean-reverting choice, similar to the OU process would be the Cox–Ingersoll–Ross (CIR) process, whose realizations are always nonnegative, in fact, it as a sum of squared OU process [1].

It is usually recommended to ensure the boundedness of the stochastic variations for the numerical implementations. In this regard, suitable resources are odd trigonometric functions such as a

or to guarantee boundedness. For example, if from (36

) we generate a process with a normal distribution

, for and , this value of variance enables us to generate a mean-reverting stochastic trajectory, whose confidence interval of () is inside the threshold of . Then, through the functions

 ξ(ηt)=sinηt,orχ(ηt)=2πarctanηt, (37)

we obtain a bounded stochastic variation inside the interval , and the OU SDEs, that generate the stochastic variations, are represented by (34b).

To couple the parameters of the system in (34a) and (34c) with a bounded stochastic disturbance, we use

 p(ηt)=p0+ρξ(ηt),

where is a constant parameter, is the stochastic process that describes the variations of the parameter, and is a factor that controls the magnitude of the perturbation.

4.3 Test Cases

In this subsection we present results of our implementation of the -based methods for the calculation of LEs at the hand of several test cases of power systems represented by strangeness-free SDAEs models of so-called single-machine-infinite-bus (SMIB) systems. This simplified model is frequently used in the area of power systems in order to understand the local dynamic behavior of a specific machine connected to a complex power network. The SMIB consists of a synchronous generator connected through a transmission line to a bus with a fixed bus voltage magnitude and angle, called infinite bus, which represents the grid. A diagram of the system is shown in Figure 6.

In each test case, we consider a different type of disturbance. For Case 1, the disturbance is a stochastic load connected to the system. In Case 2, the disturbance is due to noise caused by a measurement error in a transducer of the machine control system. In both cases, the maximum disturbance that the system can admit without loosing stability is analyzed, as well as the effect (positive or negative) of the disturbance for the system in the stable region. The whole SMIB system, i.e., the synchronous machine, system constraints, and stochastic disturbances; are modeled by a strangeness-free SDAE system. The dimension of this system is mainly defined by the type of model used in the synchronous machine; we use a classical model and a flux-decay model, see [21, 29, 35, 36] for detailed descriptions.

4.3.1 Case 1: SMIB with stochastic load

In this test case we make use of the LEs to assess the impact of stochastic disturbances associated with an active power load, over the rotor angle stability of a synchronous generator. Both the machine and load are connected to the same bus, and this bus in turn, is linked to the grid through a transmission line. This kind of SMIB model is typically used to analyze the effects of renewable energy sources, or aggregated random power consumption, see Figure 7.

For this version of SMIB system called classical model, the dynamic behavior of the synchronous machine is represented by the swing equations where the rotor angle and the rotor speed are the state variables, see [29, 36]. The algebraic constraint in the system is given by the active power balance, expressed in terms of the mechanical power, the electrical power, and the constant power consumed by the load. A stochastic process is modeled by an OU SDE. We consider that is the stochastic component of power consumption that perturbs additively the active power balance of the system, where is the size of the disturbance. This leads to the system

 dδt =[ωt−ωs]dt, (38a) 2Hdωt =[Pm−Pe−KD(ωt−ωs)]dt, (38b) ηt =−αηtdt+βdwt, (38c) 0 =E′VXeqcosδt+(PL+ρηt)−Pe. (38d)

By computing the LEs of this SDAE system and checking the LLE, we can determine the maximal perturbation size (via successive increments of ) admitted by the SMIB system before loosing rotor angle stability. The numerical tests are performed for the values ; ; ; ; ; ; ; ; ; . Most of the values are expressed in the per-unit system (pu) [35]. The methods are executed with step size and a simulation time .

Figure 8 displays the computed LLE utilizing the four -based methods for incremental disturbance sizes . As expected, at when the system is not affected by a stochastic disturbance, i.e., the system is deterministic, the computed LEs match closely with the real parts of the eigenvalues obtained from the Jacobian matrix of the linearization of (38). When increasing , all methods reveal the same monotonically increasing behavior of the calculated value of the LLE towards the unstable region. First, there is a slow increase for , and then an abrupt increase of the LLE in the interval . In the interval , even though the LLE has not yet reached the instability region, for this particular case the characteristics such as a low damping coefficient and the presence of the stochastic disturbance, provokes a behavior in the system called pole slipping. This is, in a certain sense, a different kind of instability because the system looses synchronism as it reaches another equilibrium point near another attractor, see [36, sec. 5.8] for further details. The numerical results for this case are presented in Table A.5.

This test case shows the large potential of using LEs as an indicator of instability for nonlinear power systems. These could be used also in multi-machine study cases where, however, the computational complexity has to be reduced, e.g. by model reduction.

4.3.2 Case 2: SMIB with regulator perturbed by noise

In this subsection we consider an SMIB system with a synchronous machine described by a third-order flux-decay model. Here, in addition to the rotor angle and the rotor speed associated to the swing equations, the system