{'fraction_non_alphanumeric': 0.031129147767963713, 'fraction_numerical': 0.013177369300549057, 'mean_word_length': 4.698041349292709, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Power System Simulation Using the Adomian Decomposition Method

For dynamic security assessment considering uncertainties in grid operations, this paper proposes an approach for time-domain simulation of a power system having stochastic loads. The proposed approach solves a stochastic differential equation model of the power system in a semi-analytical way using the Adomian decomposition method. The approach generates semi-analytical solutions expressing both deterministic and stochastic variables explicitly as symbolic variables so as to embed stochastic processes directly into the solutions for efficient simulation and analysis. The proposed approach is tested on the New England 10-machine 39-bus system with different levels of stochastic loads. The approach is also benchmarked with a traditional stochastic simulation approach based on the Euler-Maruyama method. The results show that the new approach has better time performance and a comparable accuracy.


I. INTRODUCTION
NCERTAINTIES exist in operations of power grids [1] Many factors such as random load consumptions and unanticipated relay protection actions contribute to the randomness of grid operations. It can be foreseen that a future power grid will have more uncertainties and stochastic behaviors in system operations due to the increasing penetrations of responsive loads and intermittent renewable generations. Thus, dynamic security assessment (DSA) of power systems should be conducted in both deterministic and stochastic manners. However, most of today's power system simulation software tools are still based on solvers of deterministic differential-algebraic equations (DAEs) that do not involve stochastic variables to model uncertainties in system operating conditions.
In literature, there are three major approaches for the modeling of a dynamic system having stochastic effects as shown in Fig. 1: the master equation, the Fokker-Planck equation [2] [3] and Gillespie method [1] [5]. The master equation and the Fokker-Plank equation are widely applied in the field of computational biology, which both focus on the evolution of probability distribution; the Gillespie method focuses on individual stochastic trajectories. The first two approaches provide a more comprehensive understanding of stochastic effects with a dynamic system but require solving In recent years, some researchers have contributed to power system simulation in a less-deterministic manner. Reference [7] proposed a systematic method to simulate the system behaviors under the influence of stochastic perturbations on loads, bus voltages and rotor speeds. This approach introduces stochastic differential equations (SDEs) to represent stochastic perturbations and solves the equations by Ito calculus, and then a mean trajectory with the envelope on trajectory variations is yielded by repeating simulations for many times. Papers [8]- [10] utilize a similar approach to study power system stability under random effects. To analyze long term stability of a power system with wind generation, a new SDE model is developed in [11], which also applies the singular perturbation theory to investigate the slow dynamics of the system with stochastic wind generation. However, the time performance of such an approach based on Euler-Maruyama method can hardly meet the requirements for online power system simulation. Especially, when the penetration of distributed energy resources (DERs) reaches a high level, the distribution network behaves in a more stochastic manner as seen from the transmission network, and hence a large number of SDEs need to be included in the power system model, which will significantly influence the simulation speed. Also, the nature of the Gillespie method requires a large number of simulations on the same model to yield the mean trajectory as well as the U envelope on variations. Therefore, adding any extra SDE to the existing set of SDEs will result in multiplying computing time by a factor of hundreds or even thousands. In our previous works [12]- [14], a new semi-analytical approach for power system simulation has been proposed. That approach applies the Adomain decomposition method (ADM) to power system DAEs to derive a semi-analytical solution (SAS) for each state variable as an explicit function of symbolic variables including time, the initial system state and other selected parameters on the system condition; then each function is evaluated by plugging in values of its symbolic variables over consecutive small time windows to make up a desired simulation period so as to obtain the simulated trajectory of each state variable. Since the form of every SAS is a summation of finite terms for approximation, its evaluation can be fast and parallelized among terms. Thus, compared to traditional numerical integration based power system simulation, this semi-analytical approach decomposes the computation into offline derivation and online evaluation of an SAS and is better fit for online power system simulation and a parallel computing environment [14]. In fact, such a semi-analytical approach also suggests a viable, alternative paradigm for fast stochastic simulation. For example, early works by Adomian in the 1970s utilized the ADM to solve nonlinear SDEs [15] by embedding explicitly stochastic processes into the terms of an SAS.
For power system simulation in a stochastic manner, this paper proposes an approach as an extension of the ADM based approach proposed in [14]. Utilizing the semi-analytical nature of an SAS yielded by the ADM, this new approach embeds a stochastic model, e.g. a stochastic load model, into the SAS. Evaluation of an SAS with the stochastic model whose parameters are represented symbolically will not increase many computational burdens compared to evaluation of an SAS for deterministic simulation. Thus, an expected number of simulation runs for one single case are achieved by evaluating one SAS for the same number of times.
The rest of this paper is organized as follows: Section II presents the SDE model of a power system that integrates stochastic loads; Section III gives the ADM-based approach for solving the power system SDEs for stochastic simulation; Section IV uses a single-machine-infinite-bus (SMIB) system to compare the fundamental difference between the ADMbased approach and the Euler-Maruyama approach in mathematics; Section V introduces a criterion for defining the stability of a general stochastic dynamical system, which is also applied to power systems; Section VI validates the proposed approach using the IEEE 10-machine 39-bus system with the stochastic loads and compares the results and time performance with those from the Euler-Maruyama approach; finally, conclusions are drawn in Section VII.

A. Synchronous Generator Modeling
For a power system having K synchronous generators, consider the 4th-order two-axis model (1) to model each generator having saliency ignored [16]. All generators are coupled through nonlinear algebraic equations (2) about the network.
In (1) and (2), ω R is the rated angular frequency; δ k , ω k , H k and D k are respectively the rotor angle, rotor speed, inertia and damping coefficient of the machine k; Y k is the kth row of the reduced admittance matrix Y; E is the column vector of all generator's electromotive forces (EMFs) and E k is the kth element; P mk and P ek are the mechanical and electric powers; E fdk is the internal field voltage; ,

B. Stochastic Load Modeling
A stochastic model can be built based on analysis on real data and assumptions on probabilistic characteristics of the stochastic variables. Traditionally, uncertainties in loads of a power system are ignored in time-domain simulation for the sake of simplicity. However their stochastic behaviors are wellrecognized in [17]. Taking stochastic loads into consideration will enable more realistic power system stability assessment.
This paper uses the Ornstein-Uhlenbeck process in [18] to model the stochastic variations of a load in these SDEs: where W(t) is the white noise vector whose dimension equals the number of load buses, a and b parameters are drifting and diffusion parameters of the SDEs, operator "  " is the Hadamard Product, i.e., element-wise multiplication, and y PL and y QL are the stochastic variations in normal distributions. The stochastic dynamic of the load is therefore modeled by (6) where P L0 and Q L0 are the mean values of the active and reactive loads, respectively.
Periodicities and autocorrelations have been observed in historical data of loads on the daily basis. However, in the time frame of seconds, loads at different substations have much lower autocorrelations. Refer to [7], this paper sets the drifting parameter on the autocorrelations of loads as 0.5 p.u./s.

A. Modeling Stochastic Variables
Consider S stochastic variables 1 ( ), , ( ) S y t y t  , which could be stochastic loads following S different distributions. Each y i (t) can be transformed by function ( ) i g  in (7) from some  i in a normal distribution. For example, if y i (t) is a load represented by a normal distribution with certain mean value, then  i specifies a zero-mean normal distribution as in (9) and ( ) i g  shifts it to around the desired mean value like in (5) and (6).
The Ornstein-Uhlenbeck process is utilized to generate each i  from (10).
B. Solving SDEs Using the ADM Consider a nonlinear system modeled by SDE (10) having M deterministic state variables x 1 , …, x M , such as the state variables of generators, exciters and speed governors, and S stochastic variables y 1 , …, y S .
To solve ( ) t x , the procedure in [14] can be used. First, apply Laplace transformation to (10) Then use (13) and (14) to calculate the Adomian polynomials under the assumption of (12), Recursive formulas (15) and (16) can be derived by matching terms of ( ) t x and ( )  f : The next step is to apply inverse Laplace transform to both sides of (15) and (16) to calculate the N-th order SAS of (10): x y x y (17) In the resulting SAS, stochastic variables in y appear explicitly as symbolic variables.

IV. COMPARISON BETWEEN THE EULER-MARUYAMA APPROACH AND THE PROPOSED APPROACH
This section applies both the Euler-Maruyama approach and the proposed ADM-based approach to the SMIB system with a stochastic load shown in Fig. 2 to illustrate the fundamental difference between the two approaches. The stochastic load is connected to the generator bus and has its resistance R L and reactance X L modeled by stochastic variables. Thus, the whole system is now modeled by DEs (18a), (18b) and SDEs (18c), (18d).
Since R L and X L change stochastically, G L and B L cannot be treated as constants in (19). The variances of R L and X L depend on the values of drifting parameters a 1 and a 2 and diffusion parameters b 1 and b 2 , respectively. To find the SAS of this system, the first step is to apply ADM to DEs (18a) and (18b). Once the SAS of the system's DEs is derived, the SAS of the SDEs can be derived and incorporated into it.
For instance, the 2nd order SAS for rotor speed ω is where, The 2nd order SAS of R L can be derived using ADM as Here, B(t) is the Brownian motion starting at origin and dB(t)=W(t)dt. Similarly, the 2nd order SAS of X L is, To derive the SAS of the entire system considering both the DEs and SDEs, replace the symbolic variables in the DEs' SAS representing the stochastic variables with the SDEs' SAS, i.e., the 2nd order SAS of the system (18) can be derived by replacing the symbolic variables R L and X L in (20) with SAS (24) and (28).
For some forms of SDEs, an analytical solution may exist, which can be incorporated into the DEs' SAS to directly derive the SAS of the entire system. For example, the general expression of the SAS terms of (18c) can be written as Therefore the infinite order SAS of (20c) is Apply Maclaurin expansion of an exponential function and lemma 2.3 in [19] to (33), the solution becomes Then apply the integration by parts formula The close form solution can be found as In this case the symbolic variable R L in (20) can be replaced by (36) instead of (24).
On the other hand, for the Euler-Maruyama approach [20] [21], since the deterministic model described by (18a) and (18b) does not permit a close form solution, the sample trajectories of (18) have to be numerically computed. The numerical scheme for R L is shown in (37) and the same scheme also applies to X L .
In practice the value of ΔW is dependent of the step size Δt for integration.
There are a variety of definitions on the stability of a stochastic dynamical system in literature [23]- [26]. The definition of "asymptotic stability in probability" in [26] can be directly applied to a power system with stochastic variables. That definition is a counterpart of the asymptotic Lyapunov stability of a deterministic system.

Definition 2:
Asymptotic stability in probability An equilibrium point is said to be asymptotic stable in probability if it is stable in probability and for given (0,1) whenever ||x 0 -x eq ||<σ 0 To analyzed the stability of numerical simulation results, this paper modifies (40) to (41) so that the stability can be accessed using the results of finite time period simulations.
where t s is a predefined time instant, ||r 0 || is a small positive number.

VI. CASE STUDIES
The proposed ADM-based approach is tested on the IEEE 10-machine 39-bus New England system as shown in Fig. 3. Selected loads are assumed to change stochastically while all generators are represented by deterministic models. In each case study, the stochastic simulation result by the Euler-Maruyama approach is used as the benchmark, and the 2 nd order SASs (i.e. N=2) are used and evaluated every 0.001 s. The value of each stochastic variable is changed every 0.1 s. For each case, 100 sample trajectories are generated. The fault applied in all cases is a 10-cycle 3-phase fault at bus 3 cleared by tripping line 3-4. All simulations are performed in MATLAB R2016a on a desktop computer with an Intel Core i7-3770 3.40GHz CPU and 8 GB RAM.

A. Stochastic Loads at 5% with Low Variances
In the first case, model the loads at buses 3 and 4 (about 5% of the system load) by Ornstein-Uhlenbeck process. The variances of the loads are 2% of their mean values. The results from the ADM-based approach and the Euler-Maruyama approach are shown in Fig. 4. Among all the generators, generator 1 has the shortest the electrical distance to bus 3 and 4, hence the rotor angle of it is presented in the following results. From the simulation results, the deterministic system response is indicated by the mean value and is asymptotically stable. Use the stochastic system stability definition introduced in Section V. When the loads at buses 3 and 4 have small variances, the system behaves similar to a deterministic system, which is asymptotically stable with a probability of 0.9 (t s =15 s, r 0 =0.05 rad/s).

B. Stochastic Loads at 100% with Low Variances
In the second case, extend stochastic loads to all buses with variances equal to 2% of their mean values. As shown in Fig.  5, the simulation results from two approaches agree with each other, which reveal a less stable post-fault system response due to increased uncertainties. When all the system loads are stochastic, the system is asymptotically stable with a probability of 0.6 (t s =15 s, r 0 =0.05 rad/s). Compared to the first case having only two stochastic loads with the same r 0 value, the probability of the system being asymptotically stable reduces from 0.9 to 0.6. Therefore, when the percentage of stochastic loads increases, even though the load uncertainties are low and the equilibrium point of the system is almost the same as its deterministic model, the asymptotic stability of the system in probability downgrades. That justifies the necessity of using stochastic load models to study the stability of power systems with a high penetration of stochastic loads.   In the third case, all the loads are represented by stochastic loads and the variances of the loads are increased to 4% of the mean values. This case may represent a scenario having DERs widely deployed in distribution networks, which make the aggregated bus load seen from each transmission or subtransmission substation behave more stochastically. The simulation results from the ADM-based approach and Euler-Maruyama approach are shown in Fig. 6.

C. Stochastic Loads at 100% with High Variances
The ADM-based approach agrees with the Euler-Maruyama approach on the simulation results. Both of them show that the system loses its stability when the variance of the loads increases to 4% of their mean values. The instability is due to the cumulative effect of stochastic load variations. The 90% confidence envelope can be utilized as an indicator of the system stability. Unlike Fig. 5, the 90% confidence envelope in Fig. 6 is not bounded any more, indicating a 0.9 probability of the system losing stability. Bus voltages also reflect the impact from high load uncertainties as shown in Fig. 7 about the voltage magnitude of bus 30, denoted by V 30 . With loads of high uncertainties, the system has an increased risk of under-and over-voltage issues because the imbalance between generation and load is magnified by increased load uncertainties. That also indicates the importance of stochastic power system simulation when penetration of DERs becomes high.
From results of stochastic power system simulation, how the probability distribution function (PDF) of a system variable evolves in time during a post-contingency period can be estimated and fit into an anticipated probability distribution for analysis. As an example, if we assume V 30 to follow a normal distribution at each time instant with the mean value and variance varying with time, Fig. 8 shows the evolutions of its PDF using simulation results from both the ADM-based approach and Euler-Maruyama approach for comparison. Fig.  8a basically matches Fig. 8b, indicating the accuracy of the proposed ADM-based approach in reflecting the evaluation of the PDF. From Fig.8, as time elapses, the PDF of the bus voltage not only shifts the mean value but also increases the variance indicated by the increasing width of the shape. Such information is not available from deterministic power system simulation. The longer the system is subjected to the effect of stochastic variables the bigger variance and larger uncertainty the system has in post-contingency dynamics.

D. Variances of State Variables
To compare the accuracy of the numerical results from the ADM-based approach and Euler-Maruyama approach, the mean value and standard deviation of the 100 trajectories are compared. For case A, as shown in Fig. 9 and Fig. 10, the ADMbased approach achieves comparable accuracy as the Euler-Maruyama approach in terms of both mean value and standard deviation value.  As more loads are modeled as stochastic, the variance of state variables grows accordingly. The mean value and standard deviation of the rotor angle of generator 1 for case B are shown in Fig. 11 and Fig. 12. In case B, the standard deviation reaches its largest value 0.25 rad/s during the first swing, which is larger than the largest standard deviation 0.1 rad/s in case A.

E. Comparison on Time Performances
The time performances for cases A, B and C of the ADMbased approach and Euler-Maruyama approach are compared in TABLE I, from which the ADM-based approach takes less than 50% of the time cost of the Euler-Maruyama approach. The advantage of the ADM-based approach in time performance is more prominent when many simulation runs are required. As discussed in [14], the ADM-based approach is inherently suitable for parallel implementation, which could help further improve the time performance if high-performance parallel computers are available.

VII. CONCLUSION
This paper proposes an alternative approach for stochastic simulation of power systems. Using the SAS derived from the ADM, the stochastic effects from load uncertainties can be taken into considerations. The result from the proposed approach is benchmarked with that from the Euler-Maruyama approach. Since the evaluation of SASs is faster than the integration with the Euler-Maruyama approach, the proposed approach has an obviously advantage in time performance. This is critical when a large number of simulation runs need to be performed for simulating stochastic behaviors of a future power grid having a high penetration of DERs. The simulation results on different levels of stochastic loads show that when the level of load uncertainty is low, the deterministic simulation is still trustworthy compared to the mean-value trajectory from stochastic simulation, but, once the level of load uncertainty becomes high, the mean-value trajectory no longer represents the true behavior of the system.