{'fraction_non_alphanumeric': 0.03317511866100425, 'fraction_numerical': 0.00969273045216088, 'mean_word_length': 4.657433578292821, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of the Autocorrelation Function of the Stochastic Single Machine Infinite Bus System

Critical slowing down (CSD) is the phenomenon in which a system recovers more slowly from small perturbations. CSD, as evidenced by increasing signal variance and autocorrelation, has been observed in many dynamical systems approaching a critical transition, and thus can be a useful signal of proximity to transition. In this paper, we derive autocorrelation functions for the state variables of a stochastic single machine infinite bus system (SMIB). The results show that both autocorrelation and variance increase as this system approaches a saddle-node bifurcation. The autocorrelation functions help to explain why CSD can be used as an indicator of proximity to criticality in power systems revealing, for example, how nonlinearity in the SMIB system causes these signs to appear.


I. INTRODUCTION
Insufficient stability monitoring and situational awareness have been identified as critical contributors to recent large power system failures, such as August 14, 2003 [1] and Sept. 8, 2011 [2]. Recent deployment of phasor measurement units (PMU) produce large quantities of time-series data, which open tremendous opportunities to improve power system stability monitoring and control. However doing so requires algorithms and methods that transform these data into useful information about power system stability.
To address this problem, this paper draws from research on critical transitions in stochastic dynamical systems. The statistics of time-series data change notably in most stochastic dynamical systems as they approach critical transitions; systems are more easily perturbed from equilibrium, and take longer to return after being displaced [3], [4]. Collectively, this phenomenon is known as Critical Slowing Down (CSD), and is most easily observed by testing for autocorrelation and variance in time-series data. Increasing autocorrelation and variance have been shown to indicate proximity to critical transitions in climate models [5], ecosystems [6] and the human brain [7]. In prior work by the authors [8], CSD has been shown to be an early warning sign of critical transition in power systems. A deep understanding of CSD in power systems could lead to new tools for wide area measurement and control.
In many systems with critical transitions, slow-varying continuous parameters cause gradual trends in state variables, which generally move the system closer to, or further from, unstable operating points. Simultaneously, random perturbations cause fast changes. As a result, such systems have two time scales: fast and slow. In power systems, loads have slow predictable trends, such as load ramps in the morning hours, and fast stochastic ones, such as random load switching. Whereas conventional generators tend to vary slowly, intermittent renewable generation, such as PV on a partly cloudy day, can have fast stochastic changes. Because of this, the mathematical framework of fast-slow systems [9] and stochastic differential equations [10] can help to explain phenomena in these systems, such as the manner in which autocorrelation and variance increase with proximity to critical points. In [9], the variance and autocorrelation of state variables for several prototypical fast-slow systems is calculated using the Fokker-Planck approach and numerical simulations.
Estimating the proximity of a power system to a particular critical transition (e.g., voltage collapse) has been the focus of a number of papers in the power systems literature. References [11]- [14] present methods to measure the distance between an operating condition and voltage collapse with respect to slowmoving state variables, such as load. While these methods provide useful information about system stability, they require accurate network models all of which contain some error.
Another approach to estimating the distance to critical transitions is to identify statistical patterns in the response of a system to stochastic forcing, such as fluctuations in load, or production from renewable energy sources. To this end, a growing number of papers study power system stability using stochastic models [15]- [18]. Reference [15] models power systems using Stochastic Differential Equations (SDEs) and solves the SDEs using Itô calculus to develop a measure of voltage security. In [10], the Euler and Milstein methods for numerically solving SDEs are used to assess transient stability in power systems, given fluctuating loads and random faults. Reference [18] uses the time evolution of the probability density function for state variables in a Single Machine Infinite Bus (SMIB) system to show how random load fluctuations affect system stability.
The results above clearly show that power system stability is affected by noise in the system. However, more work is needed to identify useful statistical trends in high sample-rate measurements from power systems. Results from the literature on CSD suggest that the combination of both increased autocorrelation and variance in time-series data are needed in order to gain insight into the proximity to critical transitions from data. Reference [8] provides empirical evidence of increasing autocorrelation and variance for an SMIB model and the 9bus test case. Reference [19] shows that voltage variance at the end of a distribution feeder increases as it approaches voltage collapse. However, these results do not provide insight into autocorrelation. To our knowledge, only [20] derives an approximate analytical autocorrelation function (from which either autocorrelation or variance can be found) for state variables in a power system model. However, the autocorrelation function in [20] is limited to the operating regime very close to the threshold of instability. In this paper, we derive exact autocorrelation functions for the state variables of an SMIB system. The results provide insight into the conditions under which autocorrelation and variance signal proximity to critical transitions in power systems. We use the results to explain why CSD occurs in power systems, and describe conditions under which autocorrelation and variance signal proximity to critical transitions.
The rest of the paper is organized as follows. Section II provides a brief review of SDEs and solution methods. Section III describes our model, analytical and numerical results. Section IV discusses the implications for power system operations, and Sec. V summarizes the results and contributions of this paper.

II. A BRIEF REVIEW OF STOCHASTIC DIFFERENTIAL EQUATIONS
There are two different approaches to modeling and solving Stochastic Differential Equations (SDEs): the Itô and the Stratonovich interpretations. In the Itô interpretation [21], noise is considered to be uncorrelated, whereas in the Stratonovich interpretation [22] noise has finite, albeit very small, correlation time [23]. Itô calculus is often used in discrete systems, such as finance, though a few papers have applied the Itô approach to power systems [15], [10]. On the other hand, the Stratonovich method is often used in continuous physical systems where noise is band-limited [24]. The Stratonovich interpretation also facilitates the use of ordinary calculus, which is not possible under the Itô interpretation.
In this paper we organize our model to use the multivariate Ornstein-Uhlenbeck stochastic differential equation, as described in [23]: where Z (t) is the vector of the state variables, A and B are constant matrices, and W (t) represents an n-dimensional Wiener process. The increments of the Wiener process (dW (t)) are independent, normally distributed random variables, such that: Because B is a constant matrix in this paper, the Itô and Stratonovich interpretations result in the same solution [24].
We will follow the Stratonovich interpretation because it allows for the use of ordinary calculus.

III. STOCHASTIC SINGLE MACHINE INFINITE BUS SYSTEM
Analysis of small power system models can be very helpful for understanding the concepts of power system stability. The single machine infinite bus system has long been used for this purpose. Reference [25] explores the small signal stability of synchronous machines using the SMIB system. In [26], Wang et al. use a novel control technique to improve the transient stability and voltage regulation of a SMIB system. In the recent literature, there is increasing interest in stochastic analysis of power systems, in part due to the increasing integration of variable renewable energy sources. A few of these papers use stochastic SMIB models. In [27], it is suggested that increasing noise in the stochastic SMIB system can make the system unstable and induce chaotic behavior. Reference [18] (mentioned in Sec. I) also studied stability in a stochastic SMIB system.
In this section, we derive autocorrelation functions of the state variables for a stochastic SMIB system. Analysis of these functions provides analytical evidence for, and insight into, CSD in a small power system. Fig. 1 shows an stochastic SMIB system. Equation (3), which combines the mechanical swing equation and the electrical power produced by the generator, fully describes the dynamics of this system:

A. Stochastic SMIB System Model
where (η ∼ N (0, 0.01)) is a white Gaussian random variable added to the voltage magnitude of the infinite bus to account for the noise in the system, M and D are the combined inertia constant and damping coefficient of the generator and turbine, and E a and δ are the transient emf and the rotor angle of the generator. The rotor angle is the angle difference between the rotor position and a synchronously rotating reference axis. The reactance X is the sum of the generator transient reactance (X d ) and the line reactance (X l ), and P m is the input mechanical power. The third term in the left-hand side of (3) is the generator's electrical power (P g ). In order to test the system with varying amounts of stress, we solved the system for different equilibria, considering that generator's mechanical and electrical power are equal at each equilibrium:

B. Autocorrelation and Variance of the Differential Variables
In order to solve (3) analytically, we linearized it around the equilibrium point using the first-order Taylor expansion: )))))X' d ) )))))X l ) V g) /θg) where δ 0 is the value of the rotor angle at the equilibrium, and ∆δ is the deviation of the rotor angle from its mean value. Equation (6) is the standard form of (5), which is known as the damped harmonic oscillator equation with noisy forcing: for which the following equalities hold: Equalities in (7) show that f increases with δ 0 while ω 0 decreases with δ 0 . Rewriting (6) as (1) results in a set of two first order SDEs. Here, Z (t) = ∆δ ∆δ T where ∆δ is the deviation of the generator speed from its mean value, and W (t) is a 2-variable Wiener process. The matrices A and B in (1) are as follows: At δ 0 = π/2, ω 0 is equal to zero, so one of the eigenvalues of matrix A becomes zero. As a result, the system experiences a saddle-node bifurcation.
Following the method in [23], the solution of (6) is as follows: where ω = ω 2 0 − γ 2 is the frequency of the underdamped harmonic oscillator, and φ = arctan(ω /γ). Note that ω and φ decrease with δ 0 . Using (9) and (10), we calculated the stationary variances and autocorrelations of ∆δ, ∆δ. Note that the eigenvalues of A have negative real part before the bifurcation point since γ > 0. The variances are as follows: If γ < ω 0 , which holds until δ 0 ≈ π 2 in our system, the autocorrelation functions for ∆δ and ∆δ are as follows: · sin (ω ∆t − φ) σ 2 η where t and s are two different times such that t > s, and ∆t = t − s.

C. Numerical example
Using (11)- (14), we calculated the variances and autocorrelations of ∆δ, ∆δ at different equilibria; see Figs. 2, 3. Here, the bifurcation parameter is the input mechanical power P m . The parameters are given below: Note that M = 2H/ω s , where H is the inertia constant, and ω s is the rated speed in rad/s. Fig. 4 shows the autocorrelation function of ∆δ for different values of P m . It shows that choosing ∆t close to 1/4 of the smallest period of the function allows one to observe the monotonic increase of the autocorrelation as P m increases. For larger values of ∆t, the increase of the autocorrelation may not be observed since the frequency of the function's oscillations varies with P m . For smaller values of ∆t, the increase of the autocorrelation is less noticeable since the function curves become closer to each other as the time lag approaches zero. We chose ∆t = 0.1s.
In Figs. 2 and 3, the analytical and numerical results are compared with each other. In order to calculate the numerical results, (3) was solved using a fixed-step trapezoidal ordinary differential equation solver. At each time step, η (t) changes according to its normal probability density function. The minimum period of oscillations in this system (T = 2π/ω ) is 0.4 sec. We chose the integration step size to be 0.01 sec which is much shorter than the period of the shortest oscillation. For each equilibrium, we integrated (6) 100 times; the average results are shown in the plots. The numerical results are shown for the range of the bifurcation parameter values for which the numerical solutions were stable. The ratio q 4 /q 1 in Figs. 2 and 3 is equal to the value of the variance or autocorrelation for P m = 4pu divided by the corresponding value for P m = 1pu. Fig. 2 shows that the variances of ∆δ and ∆δ increase with P m , and seem to be good indicators of proximity to the bifurcation. However, the growth rates of the two variances are different. The difference becomes more significant near the bifurcation where the variance of ∆δ increases much faster than the variance of ∆δ. This is caused by the term ω 2 0 in the denominator of the expression for the variance of ∆δ in (11) . In Fig. 3, the autocorrelations of ∆δ and ∆δ increase with
Figs. 5,6 show the variance and autocorrelation of ∆V g , ∆θ g at different equilibria for this system. In Fig. 5, although the variances of both variables increase with P m , the increase of the variance of ∆θ g is much more significant. Also, the variance of ∆V g does not increase with P m until the system gets close to the bifurcation while the variance of ∆θ g increases even if the system is far from the bifurcation. Both autocorrelations of ∆V g and ∆θ g in Fig. 6 increase with P m . However, the ratio q 4 /q 1 is much larger for ∆V g than ∆θ g .
IV. DISCUSSION   Fig. 7 shows the trajectory of the eigenvalues of the state matrix A. The system passes through a saddle-node bifurcation as the mechanical power is increased. Near the bifurcation, the eigenvalues are very sensitive to the change of the bifurcation parameter. As a result, the system is in the overdamped regime (ω 0 < γ) for much less than 0.1% distance in terms of P m to the critical transition. This is consistent with [20], where the autocorrelation function is valid when the system is within 0.1% to the saddle-node bifurcation such that the lowest eigenvalue determines the system dynamics. Accordingly, it can provide a good estimate of the autocorrelation and variance of the system states for a very short range of the bifurcation parameter, but the result of [20] can not be used as an early warning sign.
Figs. 2-6 show that the variance and autocorrelation of all four state variables increase when the system is more stressed. This demonstrates that CSD occurs in this system as it approaches the bifurcation, as suggested both general results [9], and prior work for power systems [8], [20].
In addition to validating these prior results, several new observations can be made. For example, the signs of CSD are more clearly observable in some variables, and not others. While all of the variables show some increases in autocorrelation and variance, they are less clearly observable in ∆V g . The variance of ∆V g decreases with P m slightly until the vicinity of the bifurcation. On the other hand, the variance of ∆θ g always increases with P m . Fig. 5 shows the two terms in the expressions for σ 2 ∆Vg and σ 2 ∆θg , from (25) and (26). The first term in σ 2 ∆Vg is always dominant and growing. On the other hand, the second term in σ 2 ∆Vg is larger for small P m . The fact The two terms comprising the variances in (25) and (26) are also shown. that this term decreases with P m can be observed from C 2 in (18). Accordingly, the decrease of C 2 with P m causes σ 2 ∆Vg to decrease until P m reaches the vicinity of the bifurcation. With respect to σ 2 ∆θg , the second term in (26) increases with P m , and adds to the increase caused by the first term. In conclusion, the variance of ∆θ g is a better indicator of proximity to the bifurcation.
The results also show that the nonlinearity of this system causes CSD to occur within it. In (8), one of the elements of the state matrix (−ω 2 0 ) changes with P m because of the nonlinear relationship between the electrical power (P g ) and the rotor angle, and that element causes the eigenvalues to change with P m . If the relationship between P g and δ were linear, the state matrix would be constant. In [23], it is shown that the stationary time correlation matrix of (1) can be calculated using the following equation: where σ is the covariance matrix of the state variables. From (28), it can be inferred that the normalized autocorrelation matrix only depends on A and the time lag. As a result, if the state matrix is constant, the autocorrelations for an specific time lag will be constant as well. Accordingly, in this system, CSD is caused by the nonlinear relationship between P g and the rotor angle. Further research is required to design methods that can use CSD to accurately determine whether a particular state is approaching instability in larger power systems. However, the following are some ideas for using this method in a real power system. Using the analytical approaches in this paper in combination with data, we expect to be able to estimate baseline (normal) values for the variance and autocorrelation of system variables. Doing so will require some extension of the analytical approaches in this paper. While it is unlikely that this approach can be used to produce explicit autocorrelation functions, the method should be able to produce numerical values for autocorrelations and variances. Also, using the approach outlined here, the ratio of the statistics in critical conditions to those in normal conditions can be calculated from off-line studies. Using the base values and critical ratios, one can calculate thresholds for CSD signs, which help in informing system operators regarding a critical condition. Calculating autocorrelation and variance from data is not computationally expensive, which makes this method wellsuited for real-time application in large power grids.

V. CONCLUSION
In this paper, we derive analytical autocorrelation functions for the state variables in a stochastic single machine infinite bus system. The analytical results were confirmed by numerical simulations. The functions show that this system does indeed express critical slowing down, as evidenced by both increased autocorrelation and variance, before the bifurcation occurs. We also showed that the occurrence of CSD in this system is due to its nonlinearity. Moreover, the results reveal that there are some differences in the growth rate of the CSD signs for different state variables; some variables are better indicators of proximity to the bifurcation than others. These findings suggest ways to better understand CSD in large multimachine power systems. Understanding the statistical behavior of stochastic power systems as they approach instability should allow for the development of new indicators of power system stability based on the statistical properties of PMU data.