Nonparametric Stochastic Analysis of Dynamic Frequency in Power Systems: A Generalized Itô Process Model

Renewable energy generation is regarded as one of the most important stochastic resources in modern power systems, of which the uncertainty can be accurately quantified by nonparametric probabilistic forecasting [18]. Without loss of generality, renewable energy discussed in this paper mainly focuses on wind power. As an important form of nonparametric probabilistic forecasting, quantiles are used to describe the predictive uncertainty of renewable generation without any probability distribution assumption. The cumulative probability distribution function (CDF) of stochastic resources is defined as $F$, and the corresponding quantile $q^{\alpha}_{t}$ is defined as

where $\rm{Pr(\cdot)}$ represents the probability operator,$x_{t}$ is stochastic resource value at time $t$, $\alpha$ is the nominal proportion of the quantile $q^{\alpha}_{t}$. The series of predictive quantiles for stochastic resources can be obtained by direct quantile regression approach [3], expressed as

where $\hat{q}_{t}^{\alpha_{r}}$ represents the estimation of actual quantile $q^{\alpha}_{t}$, and $\hat{F}_{t}$ is the predictive quantile series with nominal proportion $\alpha_{r}$ need to be estimated.

Gaussian mixture model (GMM) is a typical nonparametric model to describe probability distribution with the linear combination of several Gaussian distribution functions [19]. Theoretically, GMM can fit any type of distribution. The probability density function (PDF) of the one-dimensional GMM $f_{\rm{GMM}}\left(x_{t}\mid\theta\right)$is expressed as

where $\omega_{i}$, $\mu_{i}$ and $\sigma_{i}^{2}$ represent the weight, expectation and variance of the i-th sub-Gaussian component, respectively, and $N$ is the number of sub-Gaussian components in GMM. For each GMM, the set of parameters $\theta$ is unknown, expressed as:

Given the sub-Gaussian components number of GMM, the expected maximization (EM) algorithm is used to estimate probability model parameters with hidden variables [20], which consists of two steps in each iteration:

1. (1)

   Step 1 (E-step): Based on the current parameters, calculate the probability $\gamma_{ij}$ that each data $x_{t,j}$ comes from the i-th Gaussian component, expressed as:

   $$\gamma_{ij}^{s}=\frac{\omega_{i}^{s}\frac{1}{\sqrt{2\pi\left(\sigma_{i}^{2}\right)^{s}}}\exp\left[\left(x_{t,j}-\mu_{i}^{s}\right)^{2}/2\left(\sigma_{i}^{2}\right)^{s}\right]}{\sum_{i=1}^{N}\omega_{i}^{s}\frac{1}{\sqrt{2\pi\left(\sigma_{i}^{2}\right)^{s}}}\exp\left[\left(x_{t,j}-\mu_{i}^{s}\right)^{2}/2\left(\sigma_{i}^{2}\right)^{s}\right]}$$

   (6)

   where $s$ is the number of iterations.

2. (2)

   Step 2 (M-step): Assuming that the result of (6) is true, the estimated value of the parameter to be evaluated is calculated according to the maximum likelihood method. The formulas are given as follows:

   $$\mu_{i}^{s+1}=\sum_{j=1}^{M}\left(\gamma_{ij}^{s}x_{t,j}\right)/\sum_{j=1}^{M}\gamma_{ij}^{s}$$

   (7)

   $$\left(\sigma_{i}^{2}\right)^{s+1}=\sum_{j=1}^{M}\gamma_{ij}^{s}\left(x_{t,j}-\mu_{i}^{s+1}\right)^{2}/\sum_{j=1}^{M}\gamma_{ij}^{s}$$

   (8)

   $$\omega_{i}^{s+1}=\sum_{j=1}^{M}\gamma_{ij}^{s}/M$$

   (9)

where $M$ is the number of training sample $x_{t,j}$.

Repeat the calculation of the two steps until the result of the parameter to be solved converges, and the maximum likelihood solution of the GMM parameter is obtained.

The original EM algorithm relies on the selection of initial values, which has a slow iteration speed. In this paper, the initial dataset is divided into $N$ different classes $\Omega_{N}$ by the k-means clustering algorithm, expressed as

The expectation and variance of each class $D^{(i)}$ are used as the initial expectation and variance of the EM algorithm, and the data proportion $M_{i}/M$ in each class $D^{(i)}$ is used as the initial weight. This can reduce the sensitivity of the EM algorithm to initial values and the possibility of falling into local optima, while improving the iteration speed.

The classic Itô process has been used to express wind power with an arbitrary parametric probability distribution [17], such as Gaussian, Beta, Weibull, etc. The SDE expression in Itô process is consistent with the ordinary differential equation (ODE) of the SFR dynamic model, which would benefit for unifying description of stochastic SFR. In this paper, without loss of generality, wind power is considered as a stochastic resource, here let $P_{\rm{w}}$ represent the wind power, the following SDE can be obtained as

where $m(\cdot)$ is the drift function driving $P_{\rm{w}}$ to a set point, $\tau(\cdot)$ is the diffusion function describing the stochastic characteristics, and $W_{t}$ is a standard Wiener stochastic process [21]. It can be found that the Itô process is actually an integral with respect to the standard Wiener stochastic process.

Given a certain parametric probability distribution, the corresponding Itô process can be constructed via the method in [21]. As the functions $m(\cdot)$ and $\tau(\cdot)$ are not unique, let the drift function $m(\cdot)$ be a linear function of the stochastic resource for easy calculation, expressed as

where $\lambda_{\mathrm{w}}$ is an optional positive real number, let it equal to $1$ in this paper, and $c$ is a constant related to parameters of the corresponding probability distribution, especially for Gaussian distribution $c$ is the corresponding expectation.

Then the diffusion function can be calculated by using the following formula, expressed as

where $p_{(}\cdot)$ is a given probability density function (PDF), and $z$ is an auxiliary variable [21].

The above classic Itô process models can only describe specific parametric probability distributions by constructing drift and diffusion function of corresponding probability density functions, while the uncertainty of actual wind power cannot be accurately approximated by a specific parametric distribution. Therefore, GMM $f_{\mathrm{w}}\left(P_{\mathrm{w}}\mid\theta_{\mathrm{w}}\right)$ is utilized to describe the nonparametric probability distribution of wind power $P_{\rm{w}}$, which can be represented as a linear combination of several Gaussian distributions, expressed as

where $n_{\mathrm{w}}$ is the number of sub-Gaussian components, $\theta_{\mathrm{w}}$ represents the corresponding GMM parameter set, and the i-th sub-Gaussian component is expressed as

where $\mu_{\mathrm{w},i}$ is the corresponding expectation, and $\sigma_{\mathrm{w},i}^{2}$ is the corresponding variance.

According to (12)-(14), Itô process model corresponding to sub-Gaussian component of wind power $P_{\rm{w}}$ can be obtained, which are described as follow

By decomposing the probability distribution of wind power $P_{\rm{w}}$ into $n_{\mathrm{w}}$ sub-Gaussian components, a generalized Itô process, which can describe any probability distribution unlike classic Itô processes, can be represented by a linear combination of $n_{\mathrm{w}}$ sub-Gaussian Itô processes, expressed as

where each sub-Gaussian Itô process has a corresponding sub-Gaussian distribution, the weight of the i-th sub-Gaussian Itô process is the same as the i-th sub-Gaussian component of wind power $P_{\rm{w}}$, which is equal to $\omega_{\mathrm{w},i}$.

The frequency response model of power systems, which mainly includes generator, turbine, governor and load, is a closed-loop control system. In this paper, it aims to study the overall system frequency response characteristics, so the dispersion of frequency and angle stability are not considered [10]. The secondary frequency regulation and detailed nonlinear links such as complex governor control, amplitude limit and dead band are also neglected [9]. The SFR model is widely used in system dynamic frequency analysis, where a complex power system dynamic model is simplified to a low-order system model, and all synchronous generators can be simplified to the closed-loop control model [22].

The transfer function of SFR model is shown in Fig. 1, where $H$ is the inertia time constant of synchronous generators, $D$ is the damping coefficient, $a$ is the turbine characteristic coefficient, $T$ is the turbine time constant, $R$ is the governor regulation coefficient, and $\Delta P_{m}$ is the mechanical power increment.

For the power system with high penetration of renewables, it is necessary to consider the frequency control of renewable generation. As an important way for renewable energy to participate in frequency regulation, virtual synchronous generator (VSG) can design control algorithm of grid-connected converter by simulating the external characteristics of the synchronous generator, such as inertia, damping and active power frequency regulation [23].

The inertial support power of traditional synchronous generators can be expressed as

where $\Delta P_{e}$ is the inertia support power of the synchronous generator, $P_{N}$ is the rated power of the synchronous generator and $f_{0}$ is the reference frequency of the power system.

With applying the virtual inertia control to the inverter [23], wind turbine VSG realizes inertial support by (20), and its corresponding expression is represented as follow

where $\Delta P_{H,\mathrm{w}}$ is the inertia support power of the wind turbine VSG, $H_{\mathrm{w}}$ is the equivalent virtual inertia time constant of wind turbine VSG and $P_{\mathrm{w}}^{\max}$ is the rated power of the wind turbine VSG.

The wind turbine VSG can participate in primary frequency regulation by reserving part of spare power. The expression of primary frequency regulation power can be simplified as a linear function of the system dynamic frequency, expressed as

where $\Delta P_{k,\mathrm{w}}$ is the support power of the primary frequency regulation, $\delta_{\mathrm{w}}$ is the equivalent regulation coefficient of the wind turbine VSG, and $\Delta f$ is the system dynamic frequency denoting the system frequency deviation between the system frequency $f$ and its reference value $f_{0}$.

To suppress the influence of unbalanced torque on the turbine, a first-order inertia link with a time constant of $T_{w}$ is usually added after (21) and (22) [24], which is shown in Fig. 2.

Due to the time constants of VSGs being approximately 2-3 orders of magnitude lower than that of synchronous generators, it can be assumed that $T_{w}\approx 0$ [25]. Therefore, a new VSG-SFR model can be proposed by considering the participation of wind turbine VSG in frequency control, shown in Fig. 2.

In Fig. 3, $H_{s}$ is the system equivalent inertia time constant, expressed by

where $K$ is the proportion of synchronous generator capacity to the total capacity, and $K_{1}$ is the proportion of wind turbine capacity with VSGs.

The penetration of renewable energy can be obtained and expressed as

where $K_{2}$ represents the proportion of wind turbine capacity without VSGs. After the transfer function operation, the VSG-SFR model can be transformed into the simplified model, shown in Fig. 4.

In the simplified model, the equivalent turbine characteristic coefficient $a_{s}$ and the equivalent governor regulation coefficient $R_{s}$ can be calculated by the formulas expressed as follow

The simplified VSG-SFR model in Fig. 4 can be expressed as an ordinary differential equation, given by

where $t_{g}$ represents the system state of the governor.

The system inertia model under stochastic resources can be formulated as

where $\Delta P$ denotes the stochastic power fluctuation, expressed as

where $P_{\mathrm{G}}$ denotes the power of synchronous generators, $P_{\mathrm{w}}$ is the wind power which is a stochastic resource, and $P_{\mathrm{L}}$ represents the electricity load.

According to (12), (27) and (28), a unified Itô process of the system inertia model can be established and expressed as

Considering the simplified governor model with wind turbine support and linearizing the drift function, a stochastic model of VSG-SFR can be obtained and shown in Fig. 5. According to (26) and (28), the stochastic model of VSG-SFR can be expressed as (30).

These abovementioned variables in (30) are denoted as a vector $\textit{{X}}_{t}=\left[t_{g},\Delta f,P_{\mathrm{w}}\right]^{\mathrm{T}}$, where the superscript “T” represents the transpose operation. The SDEs (29) can be rewritten in a general form, given as

where $\textit{{x}}_{0}$ is the initial value of $\textit{{X}}_{t}$. The state matrix can be expressed by

The diffusion function vector can be expressed by

The complex stochastic characteristic of the wind power $P_{\mathrm{w}}$ can be well expressed in the form of a nonparametric distribution. For nonparametric probability distributions, the diffusion function $\tau\left(\textit{{X}}_{t}\right)$ of SDE (19) is an unknown nonlinear function which cannot be solved by neither direct solution method [16] nor series expansion method [17]. Based on the generalized Itô process model (11), the original complex nonlinear SDE (19) can be converted into a linear combination of $n_{\mathrm{w}}$ linear SDEs corresponding to $n_{\mathrm{w}}$ Gaussian distributions [26]. The i-th SDE is expressed as (34).

Write (34) in a composite form as follow

where $\textit{{B}}_{i}$ and $\textit{{c}}_{i}$ are constant vectors of the i-th SDE, and $\textit{{X}}_{t}^{i}$ is the i-th component of $\textit{{X}}_{t}$.

According to the linearized stochastic theory in [26], rewrite the composite formula (35) as

According to the Itô formula, differentiate $e^{-\textit{{A}}t}d\textit{{X}}_{t}^{i}$, it can be obtained as follow

Substituting (35) into (36), it can be obtained as

By integrating both sides of (38), the solution of (34) can be deduced as

where $e^{\textit{{A}}t}$ is an exponential function of the nth-order square matrix $\textit{{A}}t$.

The system states $\textit{{X}}_{t}^{i}$ can be proved to follow the Gaussian distribution, because $\int_{0}^{t}e^{\textit{{A}}(t-s)}\textit{{B}}_{i}dW_{s}$ follows the Gaussian distribution [27].

The expectations and variances of (39) can be calculated by

where $\mathrm{E}\left\{\textit{{X}}_{t}^{i}\right\}$ is the expectation vector of $\textit{{X}}_{t}^{i}$, $\operatorname{var}\left\{\textit{{X}}_{t}^{i}\right\}$ is the variance matrix of $\textit{{X}}_{t}^{i}$, $\lambda_{k}$ and $\lambda_{j}$ are the k-th and j-th eigenvalue of A, P is a square matrix whose columns are the independent eigenvectors of A, $\Lambda$ is a square matrix whose k-th or j-th diagonal entries are the $\lambda_{k}$ or $\lambda_{j}$, and “$\circ$” is the Hadamard product. The above derivation can be referred to [27].

At time $t$, the i-th sub-Gaussian component PDF of the system dynamic frequency $\Delta f_{t}$ can be described as

where $\mu_{\Delta f,i}$ is the expectation of the i-th sub-Gaussian component of $\Delta f_{t}$, which is actually the second entry of $\mathrm{E}\left\{\textit{{X}}_{t}^{i}\right\}$, and $\sigma_{\Delta f,i}^{2}$ is the variance of the i-th sub-Gaussian component of $\Delta f_{t}$, which is just the second diagonal entry of the $\operatorname{var}\left\{\textit{{X}}_{t}^{i}\right\}$.

The VSG-SFR model proposed in this paper is actually a linear and time-invariant system. The weight of the i-th sub-Gaussian component of the system dynamic frequency $\Delta f_{t}$ remains the same as the weight of the i-th sub-Gaussian component of wind power $P_{\mathrm{w}}$ according to the linear invariance of GMM [28] and the stochastic dynamics theory [29]. The PDF of system dynamic frequency component obtained from each SDE (44) is weighted and integrated to obtain the general probability distribution of the system dynamic frequency $\Delta f_{t}$, expressed as

Accordingly, the CDF of the system dynamic frequency $\Delta f_{t}$ is given as follow

The procedure of the proposed NSA-GIP method for power system dynamic frequency is shown in Fig. 6. in general, there are five steps given as follows

1. Step 1)

   Based on nonparametric probabilistic forecasting [3], the uncertainty of future wind power $P_{\rm{w}}$ can be quantified by quantiles without needs of any specific parametric distribution.

2. Step 2)

   According to the probabilistic forecasting results of Step 1, GMM of wind power $P_{\rm{w}}$ can be established, and the EM algorithm initialized with $k$-means clustering algorithm is used to obtain the parameter set $\theta_{\rm{w}}$.

3. Step 3)

   According to the GMM of wind power $P_{\rm{w}}$, the generalized Itô process model is established as a linear combination of several Gaussian Itô process models.

4. Step 4)

   A simplified VSG-SFR model is established and described via a nonlinear SDE, and the complex nonlinear SDE is decomposed into a linear combination of several linear SDEs based on the generalized Itô process of Step 3.

5. Step 5)

   Each linear SDE is solved to obtain each sub-Gaussian component of the system dynamic frequency.

6. Step 6)

   The sub-Gaussian components of system dynamic frequency obtained from linear SDEs in Step 5 are weighted and integrated to obtain the general probability distribution of the system dynamic frequency.

At first, the accuracy of the proposed VSG-SFR model needs to be verified by comparing with the uniform system frequency response (USFR) model proposed in [25] based on Matlab/Simulink. Meanwhile, the improvement of the proposed model in suppressing frequency changes compared to SFR model also needs to be verified. The VSG-SFR model is simulated as the first case, which only has a thermal machine and a wind turbine. The parameters of the model are derived from the parameter aggregation of the USFR model and shown in Table I.

The VSG-SFR model, USFR model and traditional SFR model without active support of wind power are respectively subject to a constant power disturbance, while changing the proportion of wind power to obtain the system frequency response curves under different renewable energy penetration levels, including 20%, 30%, 40% and 50%, which are shown in Fig. 6. It can be seen from Fig. 7 that no matter how the penetration changes, the VSG-SFR model is always almost identical to the USFR Model, with maximum error ranging from 1.4% to 3.2%, which is sufficient to demonstrate the accuracy of the proposed model. Under the same renewable penetration level, the VSG-SFR model can increase the frequency nadir and quasi-steady frequency and reduce the rate of change of the frequency through the active support of wind power. By comparing the differences of Fig. 7 (a)-(d), it can be found that with the increasing penetration of renewable energy, the VSG-SFR model has better regulation effect on the system dynamic frequency compared with the traditional SFR model. It indicates the validity of the VSG-SFR model for power system with high penetration of renewable energy.

To verify the validity of the proposed method in the VSG-SFR model, the standard variances of the system dynamic frequency deviation over time can be calculated based on Monte Carlo simulation (MCS) and the proposed method. The number of MCS simulations is set to 20000. The standard variance curve of the system dynamic frequency deviation is depicted in Fig. 8. It can be seen from Fig. 8 that the results of the proposed method match well with those of MCS. Moreover, the standard deviation of frequency deviation is very small near the initial time, and its probability distribution is concentrated near the initial value. However, after a certain time, the standard deviation of frequency deviation will tend to be stable, and its probability distribution will converge to a stable distribution.

In order to verify the effectiveness and accuracy of the proposed method, it is necessary to compare whether the PDF and CDF calculated by the proposed method are consistent with those obtained by MCS. Here, the PDF and CDF of the system dynamic frequency deviation at 5s are selected for comparison with those obtained by MCS, which is shown in Fig. 8.

The standard variances of MCS and NSA-GIP are respectively 0.0039 and 0.0040, whose error rate is only 2.56%. Moreover, it can be seen that the PDF and CDF curves of the proposed method match very well with those of MCS.

The validity of the proposed method under different VSG-SFR model parameters needs to be further verified. Since the uncertainty of dynamic frequency with the nonparametric probability distribution, the proportion deviation (PD) [30] is introduced herein to comprehensively measure the accuracy of the proposed method. The proportion deviation of the quantile $q_{x}^{\alpha}$ is defined as

where $\alpha$ represents the nominal proportion level, $N$ is the total number of samples, and $\eta_{i}$ is the indicator function of the i-th sample, described as follow:

where $x_{i}$ is the i-th sample of the system dynamic frequency. Obviously, the closer the quantile deviation is to 0, the more accurate the estimated probability distribution is.

Fig. 9 shows the PD curves of the system dynamic frequency probability distribution under different VSG-SFR model parameters, including inertia time constant $H$, turbine characteristic coefficient $a$, governor regulation coefficient $R$, renewable energy penetration $1$-$K$, virtual inertia time constant $H_{\rm{w}}$, and wind turbine virtual regulation coefficient $\delta_{\rm{w}}$.

It can be seen that the error between the dynamic frequency probability distribution obtained by the proposed method and the reference distribution is within 3% regardless of the model parameters, so the effectiveness of the NSA-GIP method is almost unaffected by parameter changes.

In order to analyze the influence of renewable energy penetration on the proposed method, different cases wind power penetration ranging from 20% to 60% are further studied with an increment of 10%. The comprehensive estimation performance indexes of dynamic frequency uncertainty are shown in Table V. It can be seen that no matter how the renewable energy penetration changes, the maximum PD value can always be maintained below 4% and the W-distance is within 0.0006, which indicates the high accuracy of the proposed method will not be influenced by in the increasing penetration of renewable energy generation.

The influence of the Gaussian component number on the proposed method is shown in Table VI. It can be seen that the more Gaussian component number, the higher the solution accuracy, the longer the time required. However, when the number of Gaussian components is greater than a certain value, the accuracy does not differ much, so it is not necessary to select a very large number of Gaussian components. Furthermore, too many sub-Gaussian components may increase the risk of over-fitting and thus reduce the generalization performance. Therefore, the choice of the component number should consider many factors such as accuracy, efficiency and model generalization performance. In the study, the number of Gaussian components is set to 10.

The time required to complete the IEEE 39-Bus system with different methods is shown in Table VII. All simulations are conducted on a computer with an Intel Core i7-7700 CPU and 16 GB memory. MCS can ensure extremely high accuracy, but its calculation time is too long with more than 500s. DSM has the highest computational efficiency with only 0.51s. Compared with MCS, the SEM calculation time also has a very significant improvement, which only takes 1.28s. However, as aforementioned, the accuracy of DSM and SEM is limited due to the assumption of specific probability distribution model. In contrast, the proposed NSA-GIP method only needs significantly short calculation time 1.24s, while having excellent accuracy, which demonstrates high potential for real-time analysis applications in modern power systems with high penetration of renewable energy.

The above results solidly show that the proposed NSA-GIP method ensures excellent accuracy and calculation efficiency. It can provide an effective support tool to ensure the frequency security of power systems with high penetration of renewables. The probability distribution of stochastic variables is described nonparametric based on the generalized Itô process, which avoids the errors introduced by traditional methods when estimating the probability distribution through parametric assumptions or simple statistical moments, and greatly improves the calculation accuracy under the premise of high calculation efficiency.