Fast Monte Carlo Simulation of Dynamic Power Systems Under Continuous Random Disturbances

Yiwei Qiu^*, Jin Lin^*, Xiaoshuang Chen^*, Feng Liu^*, and Yonghua Song^†* ^*State Key Laboratory of Control and Simulation of Power Systems and Generation Equipment,
Department of Electrical Engineering, Tsinghua University, Beijing 100084, China
^†Department of Electrical and Computer Engineering, The University of Macau, Macau 999078, China
ywqiu@mail.tsinghua.edu.cn

Abstract

Continuous-time random disturbances from the renewable generation pose a significant impact on power system dynamic behavior. In evaluating this impact, the disturbances must be considered as continuous-time random processes instead of random variables that do not vary with time to ensure accuracy. Monte Carlo simulation (MCs) is a nonintrusive method to evaluate such impact that can be performed on commercial power system simulation software and is easy for power utilities to use, but is computationally cumbersome. Fast samplings methods such as Latin hypercube sampling (LHS) have been introduced to speed up sampling random variables, but yet cannot be applied to sample continuous disturbances. To overcome this limitation, this paper proposes a fast MCs method that enables the LHS to speed up sampling continuous disturbances, which is based on the Itô process model of the disturbances and the approximation of the Itô process by functions of independent normal random variables. A case study of the IEEE 39-Bus System shows that the proposed method is 47.6 and 6.7 times faster to converge compared to the traditional MCs in evaluating the expectation and variance of the system dynamic response.

Index Terms:

Continuous random disturbance, Itô process, Karhunen-Loève expansion, Latin hypercube sampling, Monte Carlo simulation, stochastic differential equations

I Introduction

Continuous-time random disturbances due to the volatile renewable generation, e.g., wind and solar power, have a significant impact on power system dynamics. To evaluate this impact in system operation, the disturbances must be modeled with random processes instead of random variables that do not vary with time to ensure accuracy [1, 2, 3]. However, despite many studies on power system uncertainty evaluation considering random variables, e.g., probabilistic load flow (PLF), dynamic uncertainty evaluation methods for handling continuous disturbances are still inadequate.

Generally speaking, two types of methods, i.e., intrusive methods [4, 5, 6, 7, 8, 3], and nonintrusive methods [1, 2, 9] have been proposed to evaluate the impact of continuous random disturbances on dynamic power systems.

Intrusive methods use numerical computation program derived from highly specialized mathematical knowledge. In this sense, commercial power system simulation software such as PSS/E cannot be employed. For example, in [8, 3] statistical information on power system dynamics under random disturbances is characterized by partial differential equations (PDEs) based on stochastic calculus. Clearly these PDEs cannot be solved by power system simulation software, making them impractical to be used by power utility companies.

In contrast, nonintrusive methods are based on commercial simulation software, which makes them easy for power utilities to use. However, to the authors’ knowledge, currently Monte Carlo simulation (MCs) is the only available nonintrusive method for power system dynamic uncertainty evaluation in the presence of continuous random disturbances [1, 2, 9], but the large number of sampling restricts its application.

To speed up MCs in dealing with random variables, for example in probabilistic load flow (PLF), Latin hypercube sampling (LHS) has been introduced, achieving much faster convergence than simple random sampling [10, 11]. Unfortunately, LHS cannot be directly used to sample random processes, therefore cannot be used to speed up the MCs of dynamic power systems under continuous random disturbances.

In this paper, to overcome this limitation, a fast MCs method for dynamic power systems under continuous random disturbances, which first enables LHS to deal with continuous random processes, is proposed.

Detailedly, the continuous random disturbances are approximated by functions of independent standard normal random variables based on the Itô process model of the disturbances and the Karhunen-Loève expansion (KLE) of Wiener processes that drive the Itô process. Then, LHS is applied to sample the normal random variables, and disturbance signals are followingly reconstructed for simulation in commercial software. A case study in the IEEE 39-Bus System shows that the proposed method is $47.6$ times faster than the traditional MCs in evaluating expectation and variance, respectively.

II Problem Description

II-A Uncertainty Assessment of Dynamics Power Systems Under Continuous Random Disturbances

In evaluating the impact of continuous-time random disturbances from renewable generation on power system dynamics, the disturbances need to be modeled with random processes instead of static random variables that do not vary with time to guarantee accuracy [1, 2, 3]. Then, the dynamic response or performance index of the power system is determined by a function of the paths of the random processes.

Specifically, a power system under continuous random disturbances can be depicted as differential-algebraic equations with the disturbances $\bm{\xi}_{t}$ as the parameters:

	$\displaystyle d{\bm{x}_{t}}$	$\displaystyle=\bm{f}(\bm{x}_{t},\bm{y}_{t};\bm{\xi}_{t})dt,$		(1)
	$\displaystyle\bm{0}$	$\displaystyle=\bm{g}(\bm{x}_{t},\bm{y}_{t};\bm{\xi}_{t}),$		(2)

where $\bm{x}_{t}$ and $\bm{y}_{t}$ are the states and algebraic variables; $t$ is time; and $\bm{f}(\cdot)$ and $\bm{g}(\cdot)$ are the state and algebraic equations.

Fixing certain initial condition, as the solution to (1)–(2), system dynamic response such as the trajectory of rotor angle, or performance index such as CPS1 and CPS2 in automatic generation control (AGC) [3] can be described by an implicit function of the path of the excitation $\bm{\xi}_{t}$ , denoted as

\displaystyle\omega=\omega\big{(}\left\{\bm{\xi}_{\tau}\right\}_{\tau\in[0,t]}\big{)},

(3)

where $\omega(\cdot)$ is referred to as the random response function (RRF) in this paper; and $\omega$ represents the value of $\omega(\cdot)$ .

Clearly, as the paths of the disturbances randomly vary, the RRF varies accordingly. To provide assistance to system operation, statistical information on the RRF, e.g. expectation and variance, needs to be evaluated accurately and efficiently. This is the target of the uncertainty assessment.

II-B Motivation of the Proposed Method

The renowned Monte Carlo simulation (MCs) can be used to evaluate the statistical information on the RRF [1, 9]. In MCs, paths of the disturbances are repeatedly sampled and put into power system dynamic simulation software such as PSS/E to evaluate the RRF, and then statistical information is extracted. However, MCs can be computationally burdensome for practical application, since a large number of samplings are needed to achieve an acceptable precision.

To speed up sampling static random variables in MCs, for example in probabilistic load flow (PLF) problems, Latin hypercube sampling (LHS) is used [11, 9], which achieves a much faster convergence rate. However, currently LHS cannot be directly used to sample continuous-time random processes, which hinders it to be used to speed up MCs in the presence of continuous random disturbances.

To overcome this problem, this paper uses the Itô process to enable LHS to handle continuous random processes, and a fast MCs method for dynamic power systems under continuous random disturbances is proposed. The brief framework of the proposed method is shown in Fig. 1, and the detailed method are given in Sections III and IV.

III Itô Process Model of Continuous Random Disturbances

III-A Modeling Random Disturbances With the Itô Process

Refer to caption — Figure 1: Brief framework of the proposed fast MCs method for power systems under continuous-time random disturbances.

Our previous work [3] has shown that the volatile renewable generations can be modeled with an Itô process, represented by the following stochastic differential equations (SDEs):

\displaystyle d\bm{\xi}_{t}

\displaystyle=\bm{\mu}\big{(}\bm{\xi}_{t},t)dt+\bm{\sigma}\big{(}\bm{\xi}_{t},t)d{\bm{W}}_{t},

(4)

where $\bm{\xi}_{t}$ is the $m$ -dimensional random disturbances; $\bm{W}_{t}=\left[W_{1,t},\ldots,W_{n,t}\right]^{\mathrm{T}}$ is the $n$ -dimensional independent standard Wiener processes (or Brownian motions); $\bm{\mu}(\cdot):\mathbb{R}^{m}\times\mathbb{R}_{+}\rightarrow\mathbb{R}^{m}$ and $\bm{\sigma}(\cdot):\mathbb{R}^{m}\times\mathbb{R}_{+}\rightarrow\mathbb{R}^{m\times n}$ are called the drift and diffusion terms, respectively.

Remark 1.

By selecting different drift and diffusion terms, Gaussian or non-Gaussian random process $\bm{\xi}_{t}$ can be exactly characterized.

For example, for one-dimensional Itô processes subject to several typical probability distributions, the drift and diffusion terms are listed in Table I. For other distributions, the corresponding Itô process can also be constructed; for multidimensional random processes, their correlation is characterized by non-diagonal entries in the diffusion term; see [3] for details.

In addition, we present a method for identifying the Itô process model from empirical data; see Appendix.

TABLE I: Drift and Diffusion Terms of One-Dimensional Itô Processes with Some Typical Probability Distributions

Type

Probability Density

\mu\left(\xi_{t}\right)

\sigma^{2}\left(\xi_{t}\right)

Gaussian

e^{{(\xi_{t}-a)^{2}}/2b}/\sqrt{2\pi b}

-\left(\xi_{t}-a\right)

2b

Beta

\dfrac{\xi_{t}^{a-1}(1-\xi_{t})^{b-1}}{B(a,b)}

-\left(\xi_{t}-\dfrac{a}{a+b}\right)

\dfrac{2\xi_{t}\left(1-\xi_{t}\right)}{a+b}

Gamma

\dfrac{b^{a}}{\Gamma(a)}\xi_{t}^{a-1}e^{-b\xi_{t}}

-\left(\xi_{t}-a/b\right)

2\xi_{t}/b

Laplace

e^{-|\xi_{t}-a|/b}/(2b)

-\left(\xi_{t}-a\right)

2b|\xi_{t}-a|+2b^{2}

III-B Simulation of Continuous Random Disturbances

In Monte Carlo simulation, paths of the random disturbances are created via a stochastic numerical integration scheme, for example the Maryuama-Euler (EM) scheme [12]:

\displaystyle\bm{\xi}_{t+h}

\displaystyle=h\bm{\mu}\left(\bm{\xi}_{t};\bm{q}\right)+\bm{\sigma}\big{(}\bm{\xi}_{t};\bm{q})\sqrt{h}\bm{\zeta},

(5)

where $h$ is the step length; $\bm{\zeta}\sim\mathcal{N}\left(\bm{0},\bm{I}\right)$ is sampled as a vector of independent normal random variables in each step.

The stochastic integration scheme (5) also implies that time-domain discretization of the random disturbances is impractical for uncertainty evaluation since the number of the resulting random variables is proportional to the number of discretization steps, which causes the curse of dimensionality. However, discretizing the disturbances spectrally instead of in time domain can be a feasible solution. This is the basic idea of our proposed method illustrated later.

IV The Proposed Fast Monte Carlo Simulation

IV-A Spectral Representation of the Random Disturbances

According to the Karhunen-Loève theorem [12], a standard Wiener process $W_{i,t}$ can be spectrally decomposed into a series of independent standard normal random variables $\{\zeta_{i,j}\}_{j=1}^{\infty}$ , i.e., the Karhunen-Loève expansion (KLE), as

\displaystyle W_{i,T}=\int_{0}^{T}dW_{i,t}=\sum_{j=1}^{\infty}\zeta_{i,j}\int_{0}^{T}m_{j}(t)dt

(6)

where $\left\{m_{j}(t)\right\}_{j=1}^{\infty}$ are functions defined on interval $t\in[0,T]$ :

\displaystyle{m_{j}(t)=}\begin{cases}\sqrt{{1}/{T}}&,\ j=1\\ \sqrt{{2}/{T}}\cos\left[\dfrac{(j-1)\pi t}{T}\right]&,\ j\geq 2\end{cases}

(7)

For computation, the infinite series (6) is truncated at a given order $K$ . Then, taking the derivative of both sides of (6) yields

\displaystyle dW_{i,t}\approx\sum_{j=1}^{K}\zeta_{i,j}m_{j}(t).

(8)

Substituting (8) into the Itô process (4) yields the following ordinary differential equation with random coefficients:

\displaystyle\frac{d\bm{\xi}^{*}_{t}}{dt}(\bm{\zeta})

\displaystyle=\bm{\mu}\big{(}\bm{\xi}^{*}_{t},t\big{)}+\bm{\sigma}\big{(}\bm{\xi}^{*}_{t},t\big{)}\sum_{j=1}^{K}\bm{\zeta}_{j}m_{j}(t)

(9)

where $\bm{\zeta}_{j}$ is a vector of independent normal random variables with dimension $n$ ; and $\bm{\zeta}\triangleq\{\zeta_{i,j}\}_{1\leq i\leq n,1\leq j\leq K}$ is the vector of all independent normal random variables. For convenience, in the rest of this paper we reindex the entries of $\bm{\zeta}$ as $\{\zeta_{i}\}_{1\leq i\leq M}$ without ambiguity. $M=nK$ is the size of $\bm{\zeta}$ .

By (9), the continuous random disturbances are characterized by independent normal random variables. Then substitute (9) into (3), the system dynamic response or performance index defined by the RRF $\omega(\cdot)$ can be approximated by a function (denoted as $\omega^{*}(\cdot)$ ) of the normal random variables:

\displaystyle\omega\approx\omega^{*}(\bm{\zeta})=\omega\big{(}\left\{\bm{\xi}^{*}_{\tau}(\bm{\zeta})\right\}_{\tau\in[0,t]}\big{)}.

(10)

This is ensured by the following proposition.

Proposition 1.

The approximate RRF $\omega^{*}(\bm{\zeta})$ in (10) converges to the actual RRF $\omega(\bm{\zeta})^{*}\rightarrow\omega\big{(}\{\bm{\xi}_{\tau}\}_{\tau\in[0,t]}\big{)}$ , as $K\rightarrow\infty$ .

Rigorous proof of this proposition can be obtained by combining the result of Section 15.5.3 of [12] and the boundedness of the RRF (3). To save space we omit the detailed proof here.

This far, we have made the RRF be characterized by normal random variables as (10). This allows LHS to be employed to speed up finding the statistical information on the RRF.

IV-B Latin Hypercube Sampling of Random Variables in the KLE

The basic idea of Latin hypercube sampling (LHS) is to make the sampling point distribution close to the probability density function (PDF). Thus, fewer samplings are needed to achieve an acceptable precision in Monte Carlo simulation compared to simple sampling [10, 11]. LHS has two steps:

Step 1: Random Sampling. In our problem, the random variables are the coefficients $\{\zeta_{i}\}_{1\leq i\leq M}$ in the KLE (6). Given sampling size $N$ , for each random coefficient $\zeta_{i}$ , evenly partition the range of the cumulative density function (CDF) into $N$ regions, and pick a random sample uniformly in each region. Then, sampling points are obtained by the inverse CDF of a standard normal distribution. This procedure is illustrated in Fig. 2. The $N$ samplings of each random coefficient constitutes a row of an $M$ -row primary sampling matrix.

Step 2: Permutation. Random permutation on the primary sampling matrix is then performed. Permutation algorithms such as Cholesky decomposition [10] can be employed to minimize the correlation between different columns. Due to space limit, here will not exposit the permutation algorithm. Interested readers are referred to the literature.

Once the two steps are finished, the $N$ samplings vectors are obtained as the $N$ columns of the permutated sampling matrix, denoted as $\{\hat{\bm{\zeta}}_{k}\}_{k=1}^{N}$ . Then, paths of the random disturbances can be obtained using (9) by replacing the random coefficients $\bm{\zeta}$ with the sampling vectors.

Remark 2.

Because the random variables $\bm{\zeta}$ in the KLE (6) are independent, no additional transformation such as Nataf transformation is needed to make them independent. This makes our LHS much easier than that in correlated cases, for example in [9].

IV-C Overall Computation Procedure

The overall procedure of the proposed method is summarized as Algorithm 1.

As clearly can be seen, the proposed method can be easily realized on commercial simulation software, making it very easy for power utilities to use in practical operation.

Algorithm 1 Procedure of the Proposed Fast MCs Method

0: Itô process model (4) of the random disturbances, truncation order

K

of the KLE (8), sampling size

N

1: create the Latin hypercube sampling set

\{\hat{\bm{\zeta}_{k}}\}_{k=1}^{N}

following the procedure described in Section IV-B.

2: for each sampling point

\hat{\bm{\zeta}}_{k}

, use (9) to calculate the paths of corresponding disturbances

3: use power system simulation software such as PSS/E to find the RRF with respect to the paths obtained in Step 2

4: extract statistical information on the RRF such as expectation and variance from the results obtained in Step 3

V Case Study

V-A Case Settings

We compare the proposed fast dynamic uncertainty assessment method and the traditional Monte Carlo simulation method using the IEEE 39-bus system [13]. The proposed method is coded in Python with dynamic simulation performed on PSS/E. Detailed models of GENROU generators, IEEET1 exciters, and TGOV1 governors are included.

A wind farm of rated power $3,000$ MW is connected to bus 15, which acts as a continuous random disturbance imposed on the system. The following Itô process is used to model the per unit wind power, formulated as

	$\displaystyle dP_{t}=$	$\displaystyle\left[0.0535-0.0899P_{t}+0.0349P_{t}^{2}\right]dt$
		$\displaystyle+\left[-0.410+0.919P_{t}-0.505P_{t}^{2}\right]dW_{t},$		(11)

which is identified from the recorded data of an offshore wind farm [14], as shown in Fig. 3.

To validate the Itô process model (11), the probability distribution and autocorrelation of the recorded wind power data and simulated paths of the Itô process model are compared in Fig. 4. Clearly, the result shows that non-Gaussian probability distribution and temporal correlation of wind power are precisely characterized by the identified Itô process model.

Next, we investigate the dynamic performance of frequency control under continuous wind power volatility. At $t=0$ s, the system is assumed at equilibrium, at $t=1.0$ s the generator at bus $30$ is tripped to simulate a fault. The RRF to be evaluated is defined as the root mean square (RMS) value of system frequency deviation within $60$ s.

V-B Comparison Between the Proposed Method and the Traditional Monte Carlo Simulation

In this section, the proposed method is compared to the traditional MCs to exhibit the improved efficiency. The two methods are respectively used to find statistical information, i.e., expectation and variance, on the RMS system frequency deviation. In MCs, paths of the Itô process (11) are sampled using the integration scheme (5). In the proposed method, the order of KLE is set as $K=6$ , and the paths are created using (9) with the random coefficients sampled via Latin hypercube sampling. For visualization, three sampled paths used in MCs are shown in Fig. 5(a), and five paths used in the proposed method are shown in Fig. 5(b).

With different sampling size $N$ , the expectation and variance of the RMS frequency deviation obtained by the traditional MCs are presented in Fig. 6, and those obtained by the proposed method are shown in Fig. 7. Obviously, the proposed method achieves a much faster convergence than the traditional MCs either in evaluating the expectation or variance.

Quantitative comparison of the computational efficiency of the two methods is followingly performed. The maximal difference of the results with five successive sampling size is used to quantify the degree of convergence. With $N$ varies from $1$ to $1,000$ , we find that in evaluating the expectation, the proposed method only need a sampling size of $N=21$ to reach the degree of convergence of the traditional MCs with $N=1,000$ . This means the proposed method is $47.6$ times faster than the traditional MCs. As of variance, the proposed method also only needs $N=150$ to reach the convergence degree of the traditional MCs with $N=1,000$ , about 6.7 times of efficiency compared to the traditional MCs. The related values of degree of convergence are labeled in Figs. 6 and 7 respectively, and summarized in Table II as well.

The above results effectively verify the improved efficiency of the proposed method compared to the existing MCs.

TABLE II: Size of Sampling of the Proposed Method and the Traditional MCs for Convergence in Evaluating Expectation and Variance

	Method	Sampling Size $N$	Degree of Convergence
Expectation	Traditional MCs	$1000$	$2.52\times 10^{-3}$
Expectation	The Proposed	$21$	$2.23\times 10^{-3}$
Variance	Traditional MCs	$1000$	$9.83\times 10^{-5}$
Variance	The Proposed	$150$	$7.01\times 10^{-5}$

VI Conclusions

An uncertainty assessment method for dynamic power systems under continuous disturbances is proposed. The disturbances are approximated by functions of independent normal random variables, enabling Latin hypercube sampling to be applied to speed up the Monte Carlo simulation.

Applying more advanced probabilistic analysis method on the proposed spectral approximation of continuous disturbance to achieve a faster and more precise assessment of variance and high-order moments is a promising future work.

Acknowledgement

Financial supports from National Key R&D Program of China (2018YFB0905200), National Natural Science Foundation of China (51907099, 51577096, 51677100, 51761135015), and China Postdoctoral Science Foundation are acknowledged.

[Method for Identifying the Itô Process Model]

Suppose a set of recorded data with sampling interval $h$ , denoted by $\{\tilde{\bm{\xi}}_{0},\tilde{\bm{\xi}}_{h},\tilde{\bm{\xi}}_{2h},\ldots,\tilde{\bm{\xi}}_{T}\}$ . Construct the drift and diffusion terms $\bm{\mu}\left(\bm{\xi}_{t};\bm{q}\right)$ and $\bm{\sigma}\left(\bm{\xi}_{t};\bm{q}\right)$ in (4) as simple functions of $\bm{\xi}_{t}$ such as polynomials with parameters $\bm{q}$ to be identified, such that the likelihood of the following logarithmic conditional probability is maximized:

\displaystyle\max_{\bm{q}}L=\log\Pr\left[\tilde{\bm{\xi}}_{h},\tilde{\bm{\xi}}_{2h},\ldots,\tilde{\bm{\xi}}_{T}|\tilde{\bm{\xi}}_{0}\right].

(Acknowledgement1)

By the independent incremental property of the Itô process [15], the conditional probability in (1) can be reformed as

\displaystyle L=-\sum_{j=1}^{T/h}\log\Pr\left[\tilde{\bm{\xi}}_{jh}|\tilde{\bm{\xi}}_{(j-1)h}\right].

(Acknowledgement2)

Considering the Itô process (4) in its discrete form (5), and considering that the sampling interval $h$ is short, we obtain

\displaystyle\bm{\xi}_{t+h}\sim\mathcal{N}\left(\bm{\xi}_{t}+h\bm{\mu}\left(\bm{\xi}_{t};\bm{q}\right),h\bm{\sigma}^{\mathrm{T}}\left(\bm{\xi}_{t};\bm{q}\right)\bm{\sigma}\left(\bm{\xi}_{t};\bm{q}\right)\right).

(Acknowledgement3)

Therefore, the conditional probability in (2) is

	$\displaystyle\Pr$	$\displaystyle\Big{[}\tilde{\bm{\xi}}_{t+h}\|\tilde{\bm{\xi}}_{t}\Big{]}=\dfrac{1}{\sqrt{\left(2\pi h\right)^{m}\det\left(\tilde{\bm{\sigma}}_{t}\tilde{\bm{\sigma}}_{t}^{\mathrm{T}}\right)}}$		(Acknowledgement4)
	$\displaystyle\times$	$\displaystyle\exp\left\{-\frac{1}{2}\left[\Delta\tilde{\bm{\xi}}_{t}-h\tilde{\bm{\mu}}_{t}\right]^{\mathrm{T}}\left(\tilde{\bm{\sigma}}_{t}\tilde{\bm{\sigma}}_{t}^{\mathrm{T}}\right)^{-1}\left[\Delta\tilde{\bm{\xi}}_{t}-h\tilde{\bm{\mu}}_{t}\right]\right\},$

where $\Delta\tilde{\bm{\xi}}_{t}=\tilde{\bm{\xi}}_{t+h}-\tilde{\bm{\xi}}_{t}$ represents the change in the recorded data over a sampling interval; $\tilde{\bm{\mu}}_{t}$ and $\tilde{\bm{\sigma}}_{t}$ represent $\bm{\mu}(\tilde{\bm{\xi}}_{t};\bm{q})$ and $\bm{\sigma}(\tilde{\bm{\xi}}_{t};\bm{q})$ , respectively.

Substituting (4) into (2), letting $\tilde{\bm{D}}_{jh}=h\tilde{\bm{\sigma}}_{jh}\tilde{\bm{\sigma}}_{jh}^{\mathrm{T}}/2$ , and neglecting the constant terms in the logarithmic function yields

	$\displaystyle\min_{\bm{q}}L^{\prime}=$	$\displaystyle\frac{1}{4}\sum_{j}\left[\Delta\tilde{\bm{\xi}}_{jh}-h\tilde{\bm{\mu}}_{jh}\right]^{\mathrm{T}}\tilde{\bm{D}}_{jh}^{-1}\left[\Delta\tilde{\bm{\xi}}_{jh}-h\tilde{\bm{\mu}}_{jh}\right]$
	$\displaystyle+$	$\displaystyle\frac{1}{2}\sum_{j}\log\det\left(\tilde{\bm{D}}_{jh}\right)$		(Acknowledgement5)

Since (5) is an unconstrained programming problem, common methods, such as gradient descent, can be used to find the optimal parameters for the Itô process model.

References

[1] Z. Y. Dong, J. H. Zhao, and D. J. Hill, “Numerical simulation for stochastic transient stability assessment,” IEEE Trans. Power Syst., vol. 27, no. 4, pp. 1741–1749, Nov. 2012.
[2] F. Milano and R. Zárate-Miñano, “A systematic method to model power systems as stochastic differential algebraic equations,” IEEE Trans. Power Syst., vol. 28, no. 4, pp. 4537–4544, Nov. 2013.
[3] X. Chen, J. Lin, F. Liu, and Y. Song, “Stochastic assessment of AGC systems under non-Gaussian uncertainty,” IEEE Trans. Power Syst., vol. 34, no. 1, pp. 705–717, Jan. 2019.
[4] D. Apostolopoulou, A. D. Domnguez-Garc, and P. W. Sauer, “An assessment of the impact of uncertainty on automatic generation control systems,” IEEE Trans. Power Syst., vol. 31, no. 4, pp. 2657–2665, 2016.
[5] P. Ju, H. Li, C. Gan, Y. Liu, Y. Yu, and Y. Liu, “Analytical assessment for transient stability under stochastic continuous disturbances,” IEEE Trans. Power Syst., vol. 33, no. 2, pp. 2004–2014, Mar. 2018.
[6] Q. Shi, Y. Xu, Y. Sun, W. Feng, F. Li, and K. Sun, “Analytical approach to estimating the probability of transient stability under stochastic disturbances,” in 2018 IEEE PESGM, Portland, USA, 2018, pp. 1–5.
[7] H. Li, P. Ju, C. Gan, S. You, F. Wu, and Y. Liu, “Analytic analysis for dynamic system frequency in power systems under uncertain variability,” IEEE Trans. Power Syst., vol. 34, no. 2, pp. 982–993, Mar. 2019.
[8] K. Wang and M. L. Crow, “The Fokker-Planck equation for power system stability probability density function evolution,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2994–3001, Aug. 2013.
[9] W. Wu, K. Wang, G. Li, and Y. Hu, “A stochastic model for power system transient stability with wind power,” in 2014 IEEE PESGM Conf. & Exposition, National Harbor, USA, Jul. 2014, pp. 1–5.
[10] H. Yu, C. Chung, K. Wong, H. Lee, and J. Zhang, “Probabilistic load flow evaluation with hybrid latin hypercube sampling and cholesky decomposition,” IEEE Trans. Power Syst., vol. 24, no. 2, pp. 661–667, May 2009.
[11] Y. Chen, J. Wen, and S. Cheng, “Probabilistic load flow method based on nataf transformation and latin hypercube sampling,” IEEE Trans. Sustainable Energy, vol. 4, no. 2, pp. 294–301, Apr. 2012.
[12] P. K. Friz and N. B. Victoir, Multidimensional stochastic processes as rough paths: theory and applications. Cambridge University Press, 2010.
[13] Illinois Center for a Smarter Electric Grid (ICSEG). (2019, Oct.) IEEE 39-bus system. [Online]. Available: http://icseg.iti.illinois.edu/ieee-39-bus-system/
[14] J. Lin, Y. Sun, L. Cheng, and W. Gao, “Assessment of the power reduction of wind farms under extreme wind condition by a high resolution simulation model,” Appl. Energy, vol. 96, pp. 21–32, 2012.
[15] E. Pardoux and A. Rǎşcanu, Stochastic Differential Equations, Backward SDEs, Partial Differential Equations. Springer, 2014.