Power system stochastic transient stability assessment based on Monte Carlo simulation

: Electrical power system is a complicated dynamic system with plenty of random disturbances while this character has been aggravated attributed to the rapid development of renewable energy technology. Under this circumstance, traditional power system transient stability analysis methods may not always be appropriate. Thus, stochastic transient stability analysis has been established by introducing stochastic differential equations to traditional transient stability analysis. Previous studies on stochastic transient stability assessment just provided a preliminary result but failed to assess its accuracy and reliability, which are essential mathematical properties. Based on the theory of probability and statistics, the Monte Carlo simulation is utilised to further improve the numerical simulation method for transient stability evaluation. Two stability assessment indices are proposed: instability probability and expectation of critical clearing time (ECCT). Then, a corresponding method to calculate the indices and their confidence level is proposed. A two-stage adaptive method is introduced to pre-estimate the sample capacity under a given confidence coefficient, so that the corresponding confidence intervals of the indices can be calculated. The effectiveness of the proposed method is validated on the three-machine, nine-bus system, which demonstrates that it would be suitable for future hybrid AC/DC power systems with fluctuating renewable energy generations.


Introduction
The power system is a complex multidimensional non-linear dynamic system disturbed by many random elements, including the type and location of the fault, the load and the parameters of each component, resulting to its operation state changing all the time. Moreover, the continuous development of renewable energy technology injects more random sources into the power system, so inevitably the power system has to face the challenges brought by increasing randomness. On the other hand, the current model used for stability analysis is deterministic in actual engineering. Usually, it is required to verify the system stability under the most serious failure, but consequently system operation safety is guaranteed at the expense of economic efficiency. Therefore, it is necessary to introduce randomness into stability analysis.
Recently, the research on power system stochastic stability has increased sharply. There have been studies of probabilistic stability since the 1970s in which random variables were concerned [1][2][3], but nowadays more attention is paid to continuous disturbance: stochastic process. The stochastic differential equation (SDE) theory has been introduced to the stability analysis, based on which a new stability analysis system has been established. There are two main methods in the current studies: probability method and stochastic dynamics method. The main idea of the former is to transform problems into deterministic problems and to solve them by statistical methods. Hence, stability analysis theories for deterministic system still can be used: Lyapunov stability method [4] and numerical simulation algorithm [5][6][7][8]. This method is the mainstream, mainly focusing on the research of stochastic source modelling [6,[8][9][10][11][12] and numerical integration method of SDE [5][6][7][8], heavily relying on statistical theories. Among all the statistical calculation techniques, Monte Carlo simulation is most commonly used, which is simple and manoeuvrable, especially suitable for solving large-scale, high-dimensional system problems. On the contrary, the stochastic dynamics method is newly proposed, sophisticated yet efficient and capable of revealing the inherent mechanism of the system, but unfortunately its mathematical theory is not thorough enough. With plenty of difficulties for the high-dimensional, non-linear system stability and control problems to be solved, the usage of this theory applied in power system is very limited.
Considering the fact that the accuracy and reliability of the results is not analysed in current studies caused by the ignorance of statistical basis, a novel power system stochastic transient stability assessment method is proposed in this paper based on Monte Carlo simulation. Two power system stochastic stability indices are proposed, and a sample size calculation method is introduced, which can be used to evaluate the accuracy and reliability. The structure of this paper is as follows: first, the random load is modelled and the mathematical model for the whole system is established. Then, the definition of stochastic stability and the numerical calculation method of the SDE are reviewed. Section 4 introduces the two-stage adaptive Monte Carlo simulation method. Finally, the effectiveness of the algorithm is verified in the threemachine, nine-node system.

Random load model
The load can be divided into two categories: one is regular and predictable, like the load prediction curve built based on previous data and the weather forecast, while the other is random and unpredictable, which is caused by accidental and independent human activities or weather variations. According to the central limit theorem, the sum of all the stochastic parts is a random variable with Gauss distribution. Therefore, the load power at a certain interval is fluctuating around an average value. In such a short time interval, load power can be regarded as a stationary stochastic process. Considering all the factors mentioned, we can present an assumption that within a certain short time interval, the load power is a quasi-stationary stochastic process with constant mean value and variance and each random load is independent of each other. Based on this hypothesis, the load power can be divided into two parts: the deterministic part and the random part. The former is a constant number calculated from the deterministic system power flow but as for the latter considers we introduce the Ornstein-Uhlenbeck process to describe its mathematical characteristics [11], so the load power model can be established as follows: J. Eng where P 0 , Q 0 , v 0 are initial values determined by the deterministic system power flow and b, σ are two coefficients describing load fluctuation and can be extracted from the real statistical data. Now that the model of the random sources has been completed, we can establish the mathematical model for the whole system by adding a set of SDEs on to the deterministic system model: , where x, y, z are state variables, algebraic variables, and disturbance variables, respectively. x 0 , y 0 are initial values obtained from the power flow calculation. The value of z 0 may be power flow initial values or random variables of a certain probability distribution. Actually, in addition to load, disturbance variables can also include other random source: fan speed [7-10], motor parameters etc. As for fault types or other discrete random variables, Poisson process is needed [5], so more generally L t is needed instead of B t (Levy process can cover continuous stochastic processes and discrete stochastic processes simultaneously). Since this paper only discusses Wiener process, B t is suitable for the case.

Stochastic stability
After rebuilt, the system model has become a complex stochastic differential algebraic equation. As a brand-new mathematical field, the system stability has different definition. Since there is more than one definition of convergence in SDE theories, stability also has several kinds of definitions [13]. Given the following SDE, The weakest definition is stochastic stability: the zero solution of the equation is supposed to be stochastic stable if for any arbitrary r > 0 and ε ∈ (0, 1), there is a δ > 0 such that x < r}) Similar to ordinary differential equation (ODE), analytical solution of SDE is hard to obtain, so normally two kinds of methods are used in recent research: the Lyapunov method and numerical time-domain simulation both of which are commonly discussed in ODE problems. In order to acquire higher accuracy, this paper chooses time-domain simulation. More concretely, we choose theta-Euler integration method considering its better numerical stability [14]. Even though either way is extended from ODE problems, there' is an extremely huge divergence that SDE is driven by random variables so that its result is also stochastic with a certain probability distribution. Thus, statistical theories are essential to probability methods. We have to use some statistics to describe the result, such as mean and probability; otherwise, the problem remains unsolvable. It should be noted that according to statistical theories, the result is usually given in the form of confidence interval and confidence level.

Monte Carlo simulation
It has been discussed that statistical computation is crucial for solving stochastic stability problems by probability method, among which the Monte Carlo simulation is commonly used method, also known as statistical test method or random sampling technology. The main idea is to transform the problem into a statistical model so as to obtain the results by repeated tests [15]. Applying this method in power system stochastic stability assessment has the following advantages. First, for given accuracy, the sample size is irrelevant to the system scale, especially suitable for large-scale power system calculation. Second, not only the probability index can be obtained, but also expectation of critical clearing time (ECCT) and other indicators can be acquired, enriching the stability information. Third, the computer programme based on this model is relatively simple and easy to simulate all varieties of random factors such as the load, the fault, and other elements [16]. The basic idea of Monte Carlo simulation has the following four steps: 1. A simple and feasible probabilistic statistical model should be set up for practical problems and make sure the solution is exactly the probability distribution of the model or its statistical characteristics. 2. Select an appropriate sampling method for the model built in step1. Generate random numbers and repeat the calculation so as to obtain adequate sample trajectories. 3. Analyse the result and estimate its accuracy and reliability. 4. If necessary, improve the model to reduce the computational cost and increase the efficiency.

Stability indices
According to the idea of Monte Carlo simulation, we should establish a statistical model first. In this paper, two indices to measure the power system stochastic stability are presented.
1. Expectation of critical clearing time: it expresses the expectation of the critical clearing time of the system with a given fault, and the formula is as follows: 2. Probability of instability (PI): it indicates the probability that system is unstable with a given fault. It is decided by whether the power angle deviation between generators is more than 360 degrees at the end of simulation. The calculation formula is as follows:

Two-stage adaptive Monte Carlo simulation
The problem after is that the sample size as well as the corresponding accuracy and confidence level. However, troubled of the lack of prior information in practice, it is a rather tough question. We have to make use of some probability inequalities, but unfortunately, most of the probability inequalities still need a priori information: variance. Actually in practice, people use the observation function to estimate the variance, giving birth to an adaptive algorithm: choose an initial sample with a sample size n 0 and then estimate the function variance σ 2 with the sample variance v¯n 0 . After that, determine the actual sample size n according to the estimated function variance. Then the statistics of the model can be calculated. It is quite obvious that the most essential part of this method is variance estimation. To solve this problem, another statistic is required: kurtosis, defined as the ratio of the four-order centre moment to square of variance [17]. We present the following inequality without proof: is the kurtosis and g is the sample function. So the inequality above can be written as follows: P{σ up (v¯n, n, α, κ) > σ} ≥ 1 − α, σ up 2 (v¯n, n, α, κ) From the above formula, it is obvious that sample size n has to satisfy this inequality: If a variance magnification factor λ 2 is defined and its upper limit L is appointed: then sample size n has to satisfy this inequality: To sum up, the sample size calculation algorithm has the following procedures: 1. For A function g(X) with four-order moments, appoint a positive constant κ max ensuring kurt(g) ≤ κ max . 2. Appoint an error limit ε, a confidence level α, and a sample variance magnification factor upper limit L. Suppose that Find a positive number n 0 satisfying the following inequality: κ max < κ poss (α 1 , n 0 , L) = n 0 − 3 n 0 − 1 + α 1 n 0 1 − α 1 1 − 1 L 3. Calculate the sample variance v¯n 0 in which n 0 is the sample size solved before. 4. Estimate the maximum function variance σ up 2 = σ up 2 (v¯n, n, α, κ). 5. Determine the actual sample size with the estimated variance: if n μ < n 0 which is possible when the variance is quite low or L is close to 1, then, there is no need to calculate the sample trajectories again and regard the first-stage results as the actual results. 6. With the actual sample size n μ achieved, calculate the sample mean μ and the final result is within the confidence interval [μ − ε, μ + ε] of confidence level 1 − α.
In fact, the kurtosis of our problem is also unknown, so this method does not apply to every case. For the random variable of continuous distribution, such as ECCT, its probability distribution is approximately close to the normal distribution if the sample size is sufficient, whose kurtosis is close to 3 no matter how large the variance is; thus, it is reasonable to suppose the kurtosis is 4 or other number a little greater than 3 in practice. However, for random variable of discrete distribution, the kurtosis is hard to estimate. The index PI is just the case, but considering its particularity that the value of the variable is just 0 and 1, we discover that the kurtosis and the mean satisfy the following relevance (see Fig. 1).
That is to say, the more stable (or unstable) the system is, the higher the kurtosis is. Therefore, a kurtosis limit can be specified in actual calculation. If the actual kurtosis is above that limit, it can be considered to be stable (or unstable) enough and further calculation is unnecessary and resource-wasting. For example, when the instability probability is 0.01, the kurtosis is 98.01. Considering the contradiction between accuracy and computation cost, the upper limit of kurtosis can be chosen as 100. When kurtosis exceeds 100, the PI (or stability) is considered less than 0.01 and exact calculation is omitted. Therefore, the above method is still applicable.

Case study
We tested the validity of the algorithm in the three-machine, ninenode system [18][19][20]. The critical clearing time of the deterministic system is 0.163 s. The active and reactive loads of node 5, 6, and 8 are modelled as stochastic process, the coefficients of which are b = 0.02, σ = 0.02. We calculate ECCT and PI with the algorithm mentioned above, and the results are listed as below.

Calculation of ECCT using two-stage adaptive Monte Carlo simulation
Appoint the calculation precision ε = 0.001, the confidence level α = 0.05, the maximum kurtosis κ max = 5, and the sample variance magnification factor upper limit L = 2. The first-stage calculation requires 617 samples, and the estimated variance is 1.9867 × 10 −4 . According to the algorithm, the second-stage calculation needs 2074 samples and the statistics are as follows.
It should be noted that the sample size needed in the first-stage calculation is just connected to the precision, the confidence level, the maximum kurtosis, and the sample variance magnification upper limit, so for a certain index it is a fixed number; for the second-stage calculation, the sample size is influenced by the variance. As long as the variance is limited, the sample size would not be too large. It is not directly relevant to the system scale. Fig. 2 presents the calculated distribution of ECCT, and the result shows that the distribution is close to the normal distribution, so do the statistics listed in Table 1.

Calculation of PI using traditional Monte Carlo simulation
Considering that PI is a random variable with discrete distribution, its kurtosis is hard to estimate, so we use the traditional Monte Carlo simulation first to test its characteristics. We try to calculate the statistics of PI with different sample sizes, and the results are as follows (see Table 2): Actually, for a SDE, the error limit of the numerical result is difficult to acquire. The error limit above was calculated by the Chebyshev inequality. We have to point out that this error limit may be larger than the real one because this inequality is rather conservative. This test is just a reference compared to the next one.

Calculation of PI using two-stage adaptive Monte Carlo simulation
This test appoints the calculation precision ε = 0.01, the confidence level α = 0.05, and the maximum kurtosis κ max = 100, and the sample variance magnification factor upper limit L is not appointed here. As in this test we have already estimated a rather conservative kurtosis, to limit the sample variance magnification factor will further aggravate the conservation and cause redundant calculation. The sample volume calculated of the first stage is 7622. The first-stage results are as follows (see Table 3).
The second-stage calculation of sample capacity is 5332, less than the first-stage result, so no more calculation is needed as we mentioned before in the algorithm step 5. Actually, it is caused by the fact that the kurtosis we suppose in the first-stage calculation is much larger than the real one. This may be a quite common phenomenon in the calculation of PI. As a result of the rather conservative kurtosis we suppose, adequate calculation is an inevitable consequence for most of the cases, but if a smaller kurtosis is appointed and if the PI is rather large (or small), the precision will be affected. So we have to weigh accuracy against efficiency considering the reality and κ max = 100 is just an example listed here.

Conclusion
This paper introduces a thorough evaluation method for power system stochastic stability based on statistics theories and the calculation method of the confidence interval and confidence level. First, two stability evaluation indices are proposed: the ECCT and the PI. Then, considering the lack of prior information in practical problems, the Monte Carlo simulation method based on kurtosis is introduced to determine the sample size. This method is quite effective for the continuous random variable model (such as ECCT). However, for the discrete random variable model (such as PI), considering the particularity of its distribution and the contradiction between cost and accuracy, the method is still applicable only if the kurtosis is limited to a given upper bound and beyond the bound the system is supposed to be sufficiently stable (or unstable).
The most critical problem in this research is the high computation cost. So the future work will focus on computation load reduction. Quasi-Monte Carlo simulation or the variance reduction technique may be introduced to solve this problem.