1 Introduction
Modeling the dynamic behavior of systems employing differentialalgebraic 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 differentialalgebraic 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 differentialalgebraic 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 (indexone) 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 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 singlemachine infinitebus 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 based decomposition methods and their evaluation. Interesting studycases 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 DifferentialAlgebraic Equations
Consider a system of quasilinear stochastic differentialalgebraic equations (SDAEs) of the form
(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 thatSince the Wiener process is nowhere differentiable, it is more convenient to represent equation (1) in its integral form as
(2) 
Here, the first integral is a stochastic RiemannStieltjes 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 andimensional vectorvalued 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
(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 strangenessfree DAE [23] i.e. it has differentiation index one. A system with these characteristics can be transformed into a semiexplicit form by means of an appropriate kinematic equivalence transformation [26, 6]
, i.e., there exist pointwise orthogonal matrix functions
and such that, premultiplying (1) by , and changing the variables according to the transformation one obtains a system in semiexplicit form(4a)  
(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
(5) 
This equation is called underlying SDE of the strangenessfree SDAE. It acts in the lowerdimensional subspace , with (where denotes the number of algebraic equations). The SDE system (5) preserves the inherent dynamics of a strangenessfree 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.
2.2 Random Dynamical Systems generated by SDEs
In the previous section we have discussed the reduction of an autonomous strangenessfree SDAE to its underlying SDE, which preserves the dynamic characteristics of the original system. Using the backtransformation 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.
(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
which means that a shift transformation given by is measurepreserving and ergodic [8].
Together with the SDE (6), we define the variational system
(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
(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(9) 
being then a scalar cocycle over (see Theorems 2.3.32 and 2.3.3940 in [3] for details).
2.3 Lyapunov Exponents of Ergodic RDSs
An important result of the theory of RDSs is the socalled 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 denotes the positive part of . This guarantees that the variational equation (7) associated with (6) is wellposed. 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 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
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]
(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, quasiperiodic, or almostperiodic) 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 Methods for computing LEs
In this section we derive the numerical techniques to compute the finitetime 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 uppertriangular 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
(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
(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 ,
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 normpreserving property of the orthogonal matrix function , we have
(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
(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 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 statetransition matrices as
(15) 
At , we perform a standard matrix decomposition
and for , we determine as the numerical solution (via numerical integration) of the matrix initial value problem
(16) 
and then compute the decomposition
where has positive diagonal elements. From (15), the value of the fundamental matrix is determined via
which is again a factorization with positive diagonal elements. Since this is unique, for the decomposition , we have
Here we denote as the triangular transition matrices with . From (14), the LEs are thus computed as
(17) 
3.2 Continuous 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
(18)  
(19) 
Inserting (18) into the variational equation (12), and multiplying by from the left and by from the right, we obtain
(20) 
Since
is upper triangular, the skewsymmetric matrix
satisfies(21) 
This results in an SDE for given by
(22) 
where the matrices (for ) are defined via
(23) 
A corresponding SDE for can be obtained from (20) and (21) via
(24) 
and the equation for the th diagonal element is given by
(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),
(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
(27) 
By integrating this SDE, it is possible to obtain the LEs from
(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 EulerMaruyama 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 twostep 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 EulerMaruyama 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 strangenessfree SDAE system in order to compare the computational efficiency and accuracy of both the discrete and continuous methods using the numerical integration schemes EulerMaruyama and Milstein. The four numerical methods will be denoted as DEM, DMilstein, CEM, CMilstein, 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
(29) 
with . The nonlinear functions in both the drift and diffusion part are continuous on , with continuous and bounded derivatives, and is a onedimensional Wiener process. The underlying SDE of (29) is
(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])
(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 EulerMaruyma 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 smallsignal 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 timeinvariant systems or systems close to a stationary solution. When timevarying 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 timedomain 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 strangenssfree 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 quasilinear DAEs with partitioned variables, see [21],[36], of the form
(32a)  
(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 strangenessfree (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
(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 strangenessfree SDAE system of the form
(34a)  
(34b)  
(34c) 
or in simplified notation as
(35) 
with
drift , and diffusion , where .
4.2 Modeling Stochastic Perturbations
In this subsection, we discuss on the modeling of stochastic variations via SDEs. We employ the well known meanreverting process termed Ornstein–Uhlenbeck (OU) process [31, 17]. The SDE which defines the OU process has the form
(36) 
where . The OU process is a stationary autocorrelated Gaussian diffusion process distributed as . Another meanreverting 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 meanreverting stochastic trajectory, whose confidence interval of () is inside the threshold of . Then, through the functions(37) 
we obtain a bounded stochastic variation inside the interval , and the OU SDEs, that generate the stochastic variations, are represented by (34b).
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 strangenessfree SDAEs models of socalled singlemachineinfinitebus (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 strangenessfree 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 fluxdecay 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
(38a)  
(38b)  
(38c)  
(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 perunit 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 multimachine 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 thirdorder fluxdecay model. Here, in addition to the rotor angle and the rotor speed associated to the swing equations, the system
Comments
There are no comments yet.