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.


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 strangenessfree (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 welldefined LEs, we study numerical methods, based on the QR 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 QR decomposition of the fundamental solution matrix. Both techniques have been extensively studied in deterministic ODE and DAE systems, see [4,5,14,15,26,27]. This paper follows the ideas exposed in [9], where these QR 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 strangenessfree SDAEs, their relation with the SDEs, and the existence of LEs generated by such SDEs. Section 3 presents the discrete and continuous QR-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 with a singular matrix E ∈ R n×n of rank d < n. The function f 0 ∈ C k (D x , R n ) (for some k ≥ 1) is known as drift, and f 1 , . . . , f m ∈ C k+1 (D x , R n ) are the diffusions, here D x ⊆ R n is an open set. Furthermore, w j t (for j = 1, . . . , m) form an m-dimensional Wiener process defined on the complete probability space (Ω, F, P) with a filtration (F t ) t≥t 0 . Each j-th Wiener process is understood as a process with independent increments such that (w t −w s ) ∼ N (0, t − s), i.e. is a Gaussian random variable for all 0 ≤ s < t, such that Since the Wiener process is nowhere differentiable, it is more convenient to represent equation (1) in its integral form as 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 x t 0 = x 0 independent of the Wiener processes w j t and with finite second moments [32]. A solution x t = x(t, ω) of (2) is an n-dimensional vectorvalued Markovian stochastic process depending on t ∈ I and ω ∈ Ω (the parameter ω is commonly omitted in the notation of x). Such a solution can be defined as strong solution if it fulfills the following conditions, see e.g. [11,42].
Because of the presence of the algebraic equations associated with the kernel of E, 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 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 [6,26], i.e., there exist pointwise orthogonal matrix functions P and Q such that, pre-multiplying (1) by P, and changing the variables x t according to the transformation x t = Qx t one obtains a system in semi-explicit form wherex D t andx A t 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 functionf A 0 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 m j=1f A j ≡ 0, so that the algebraic equations in (4b) can be solved asx A t = F A (x D t ) and inserted in the dynamic equations (4a) yielding an ordinary SDE This equation is called underlying SDE of the strangeness-free SDAE. It acts in the lowerdimensional subspace R d , with d = n−a (where a 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 backtransformation to the original state variables can be performed.

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. where for some k ≥ 1 and is the Banach space of C k vector fields on R d with linear growth and bounded derivatives up to order k and the k-th derivative is δ-Hölder continuous. In addition, we assume that the differential operator L := f 0 + 1 2 m j=1 (f j ) 2 is strong hypoelliptic in the sense that the Lie algebra L(f 0 , f 1 , . . . , f m ) generated by the vector fields f j (with j = 0, . . . , m) has dimension d for all x t ∈ R d [3]. Once again, w j t (for j = 0, . . . , m) is a m-dimensional Wiener process, this time with the convention dw 0 t ≡ dt.
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 ϕ : R + × Ω × R d → R d over this system. The ergodic MDS is denoted by θ ≡ (Ω, F, P, (θ t ) t∈R ) with the filtration (F t ) t≥t 0 , 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 obtained after linearizing (6) along a solution. If we denote by Φ(t, ω, x) the Jacobian of ϕ(t, ω) at x t , then Φ is the unique solution of the variational equation (7), satisfying where I d denotes the identity matrix of size d. Therefore, Φ it is a matrix cocycle over Θ.

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 is satisfied, where log + denotes the positive part of log. 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 of R d into random linear subspaces L i (ω), which are invariant under Θ. Here p ≤ d, where d i ∈ N denotes the dimension of the subspace L i (ω) (with 1 ≤ i ≤ p), and p i=1 d i = d. This splitting is dynamically characterized by real numbers λ 1 > . . . > λ p which quantify the exponential growth rate of the subspaces. These are called Lyapunov exponents, and are defined by According to [2, pag. 118], the LEs λ i are independent of (ω, x) 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 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.

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,14,15] 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 where f j are sufficiently smooth functions. The corresponding variational equation of (11) along with the solutions x t (x 0 ), turned into a matrix initial value problem, is given by with the identity matrix I d ∈ R d×d as initial value, where J j (x t ) := ∂f j ∂x are the Jacobians of the vector functions f j (x t ), and V ∈ C 1 (I × R d×d ) 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 QR factorization of V t , where Q t is orthogonal, i.e., Q T t Q t = I d , and R t is upper triangular with positive diagonal elements R ii t for i = 1, . . . , d. Applying the MET theory presented in Subsection 2.3, and taking into account the norm-preserving property of the orthogonal matrix function Q t , we have where {p i } is an orthonormal basis associated with the splitting of R d . 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 The QR methods require to perform the QR decomposition of V t for a long enough time, so that the R ii t 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 QR method.

Discrete QR Method
The discrete QR method is a very popular method for computing LEs in ODEs and DAEs. In this approach, the fundamental solution matrix V t and its triangular factor R t are indirectly computed by a reorthogonalized integration of the variational equation (12) through an appropriate QR decomposition. Thus, given grid points 0 = t 0 < t 1 < . . . < t N −1 < t N = T , we can write V t in terms of the state-transition matrices as At t 0 = 0, we perform a standard matrix QR decomposition and for = 1, 2, . . . , N , we determine Z (t ,t −1 ) as the numerical solution (via numerical integration) of the matrix initial value problem and then compute the QR decomposition where R (t ,t −1 ) has positive diagonal elements. From (15), the value of the fundamental matrix V t is determined via which is again a QR factorization with positive diagonal elements. Since this is unique, for the QR decomposition V t = Q t R t , we have Here we denote as R κ the triangular transition matrices R (t ,t −1 ) with κ = 0, 1, . . . , . From (14), the LEs are thus computed as

Continuous QR Method
The implementation of the continuous QR technique requires to determine a system of SDEs for the Q factor and the scalar equations for the logarithms of the diagonal elements of the R factor elementwise. Then, once the orthogonal matrix Q is computed by numerical integration, the logarithms of the diagonal elements of R can also be obtained.
By differentiating in the Itô sense the decomposition V t = Q t R t and using the orthogo- Inserting (18) into the variational equation (12), and multiplying by Q T t from the left and by R −1 t from the right, we obtain Since (dR t )R −1 t is upper triangular, the skew-symmetric matrix dS t := Q T t (dQ t ) satisfies This results in an SDE for Q t given by where the matrices T j t (x t , Q t ) (for j = 0, . . . , m) are defined via A corresponding SDE for R t can be obtained from (20) and (21) via and the equation for the ith diagonal element R ii t is given by Since the computed LEs can be obtained from (14), we make use of the Itô Lemma to introduce the following SDE for the function ψ i t = log R ii t from (25), If we assume that there are no correlations between the diffusion terms in the SDE system, then we do not have terms dw k t dw t (for 1 < k < m, and 1 < < m, with k = ) in the SDE (26).Also, using that dt dt ≡ 0, dt dw k t ≡ 0, and dw j t dw j t ≡ dt for 1 < k < m, the SDE (26) is reduced to By integrating this SDE, it is possible to obtain the LEs λ i from The difference between the discrete and the continuous QR method is that for the first one, the orthonormalization is performed numerically at every discrete time step, while the continuous QR method maintains the orthogonality via solving differential equations that encode the orthogonality continuously.

Computational Considerations
In this section, we discuss additional aspects of the computational implementation of discrete and continuous QR methods to calculate LEs. The application of the discrete QR 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 QR technique, must be performed in such a way that it preserves the orthogonality of the factor Q 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 QR-methods in Matlab. However, to obtain a unique QR factorization in each step, we have modified the QR decomposition provided by Matlab to ensure this uniqueness, by forming a diagonal matrix I with I i,i = sign(R i,i ), for i = 1, . . . , d; and then setting Q := QI and R := IR.

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 QR 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 with α ∈ R + . The nonlinear functions in both the drift and diffusion part are continuous on R + , with continuous and bounded derivatives, and w t is a one-dimensional Wiener process. The underlying SDE of (29) is 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]) where p(x) is the stationary density of the unique invariant probability law ofx t . By solving numerically (31) for α = 2, we obtain the exact value of the LE associated to (30 and its original SDAE (29, which is λ = −1.3385. The accuracy of the QR-based methods will be assessed by comparing with this value as reference. 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 QR-based methods, the sequences of random variables λ t 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 E[λ t ] and V[λ t ] 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 λ t is such that E[λ t ] = µ λ , for all , and    This means that the limit of λ t can be interpreted as a deterministic value with probability 1. This enables us to state that the stochastic approximations λ t converge in m.s. sense to a number (a degenerate random variable), which is expected to represent the LE λ.
In te rv a l T  In Figure 3 we compare the relative error of the accuracy of the four numerical methods for different stepsize h and time interval [0, T ]. 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 h and T . Here, we observe that all the methods are affected to the same extent by incrementing the simulation interval T , via a logarithmic increment, and by narrowing the stepsizes h, via an exponential increment. A more pronounced difference between the methods should be observed in higher dimensional systems.

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 [39,40]; 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.

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 where are the dynamic state variables, and x A t ∈ R a are the algebraic state variables and we set n 1 = d 1 + a. 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 Here, f D 2 ∈ R d 2 are the stochastic variables, and w t is the Wiener process. By combining (33) and (32), and assuming that x D 2 t perturbs (32a) and (32b), we obtain a strangeness-free SDAE system of the form or in simplified notation as with E := 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).

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 [17,31]. The SDE which defines the OU process has the form where α, µ, β ∈ R + . The OU process is a stationary autocorrelated Gaussian diffusion process distributed as N (µ, β 2 /2α). 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 sin or arctan to guarantee boundedness. For example, if from (36) we generate a process with a normal distribution N (µ, σ 2 ), for µ = 0 and σ 2 = 0.16, this value of variance enables us to generate a mean-reverting stochastic trajectory, whose confidence interval of 95% (±2σ) is inside the threshold of ±1. Then, through the functions or we obtain a bounded stochastic variation inside the interval [−1, 1], 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 where p 0 is a constant parameter, η t is the stochastic process that describes the variations of the parameter, and ρ ∈ R + is a factor that controls the magnitude of the perturbation.

Test Cases
In this subsection we present results of our implementation of the QR-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.

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 δ t and the rotor speed ω t are the state variables, see [29,36]. The algebraic constraint in the system is given by the active power balance, expressed in terms of P m the mechanical power, P e the electrical power, and P L the constant power consumed by the load. A stochastic process η t is modeled by an OU SDE. We consider that ρη t 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 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 computed LLE utilizing the four QR-based methods for incremental disturbance sizes ρ = 0.00, 0.05, . . . , 2.00. As expected, at ρ = 0.0 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 0.00 < ρ < 0.60, and then an abrupt increase of the LLE in the interval 0.60 < ρ < 0.75. In the interval 0.75 < ρ < 1.20, 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.

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 δ t and the rotor speed ω t associated to the swing equations, the system includes the effect of the field flux ψ f d described by the field circuit dynamic equations and constraints. In this model the machine is equipped with an automatic voltage regulator (AVR) to keep the generator output voltage magnitude in a desirable range, and a power system stabilizer (PSS) to damp out lowfrequency oscillations, see Figure 9. The AVR and PSS add to the system three more state variables v 1 , v 2 , and v s ; together with their corresponding DAEs, which describe the dynamic behavior and constraints of the controllers into the SMIB system. The resulting model is a nonlinear system of strangeness-free DAEs. We use the LEs to analyze the system stability at a specific operation point in the state-space when it is subjected to small-disturbances. Using the small-signal stability assessment (SSSA), the set of DAEs that describes the dynamics of the power system is linearized around the desired operating point. The final result is a linear DAE system. A comprehensive explanation of this model, its linearization, and reduction to an underlying ODE system can be found in [21, ch. 12]. We consider a disturbance of stochastic nature entering in the exciter block of the AVR as an error of the reference signal [21,40], by adding the stochastic variable η to v 1 in equation (39c). Resolving the algebraic constraints leads to the linearized system of SDEs where ∆δ, ∆ω, ∆ψ f d , ∆v 1 , ∆v 2 , ∆v s , η are the state variables of the linear underlying SDE system (for simplicity, the subscript t has been omitted in this formulation). Once again, the stochastic perturbation is generated via an OU SDE, and the size of the perturbation is controlled by the parameter ρ. increasing perturbation size ρ clearly mark four defined intervals. In the leftmost interval with 0.00 < ρ < 0.40, the calculated LLE is practically constant and equal to the real part of the right-most eigenvalue from the deterministic system. In this region, there is no impact of the disturbance on the system stability. In the interval 0.40 < ρ < 1.30 a curious situation occurs, as the size of the disturbance increases, the distance from the LLE to the positive region increases, in other words, the noise improves the stability of the system. In the interval 1.30 < ρ < 2.60, the situation changes completely, and the LLE converges to zero. Finally, from ρ ≈ 2.60 onwards, the system is unstable. Finally, we have evaluated the computing-times for this 7-dimensional test case. The results are shown in Figure 11. Although the computational cost for all method is similar for the different methods as a factor of the step sizes h and time interval [0, T ], the computational costs strongly increase.

Conclusions
We have revisited the theory of strangeness-free SDAE systems, as well as the concepts of LEs associated with the RDEs generated via such SDAEs. We have adapted and implemented stochastic versions of continuous and discrete QR-based methods to calculate approximations of the LEs, and assessed them by using Euler-Maruyama and Milstein schemes over the corresponding underlying SDE. The results obtained from our numerical experiments illustrate the approximations of the corresponding LE converge to degenerate random variables, i.e., the LE can be interpreted as a deterministic value, since in the limit the variance of the approximations tends to zero. Both QR-based method provide reliable results, but in general, continuous methods provide better accuracy than the discrete counterpart at the expenses of higher computational cost and higher memory requirement. We have illustrated the QR-based methods for SMIB power system problems and shown the usefulness of the LEs as a stability indicator for the rotor angle and voltage stability analysis of power systems affected by bounded stochastic disturbances.
As future work, we suggest the use of discretization schemes for SDAEs in order to directly apply the numerical integration to the SDAE system. Furthermore, methods for computing the LEs based on Singular Value Decompositions, a combination with model reduction, and a careful comparison with QR-based methods would be of interest. Concerning to the applications to power systems and dynamical network systems in general, stability assessment of large-scale cases are remarkable works to be performed.

Appendices
In this Appendix we present the numerical values for several different simulations.