Online Evaluation for the POI-Level Inertial Support to the Grid via Ambient Measurements

Abstract

As renewable energy sources like wind and solar power increasingly replace traditional energy sources and are integrated into the power grid, the issue of insufficient system inertia is becoming more apparent. This paper presents an online adaptive time window inertia constant identification method based on ambient measurements to identify the equivalent inertia constant of the time-varying inertia at Point of Interface (POI) level. The proposed method takes advantage of the online inertia estimation and the data-driven equivalent inertia constant identification techniques to simultaneously achieve online tracking and accuracy. With this regard, this paper first describes the inertia providers in modern system. Then, based on the frequency and power data measured by the Phasor Measurement Unit (PMU), this paper provides an improved data-driven equivalent inertia constant identification method. Subsequently, the paper proposes an ambient data smoothing method to cope with the numerical errors and provides, as a byproduct, an adaptive time window inertia constant identification. The adaptive time window is designed to enhance the accuracy of the method. Finally, the feasibility and accuracy of the proposed method of tracking synthetic inertia are validated by the simulation tests based on a grid in northwest China with high renewable energy penetration and a Virtual Power Plant (VPP). The experimental results show that the accuracy of this method is within $5\%$.

1. Introduction

1.1. Motivation

Due to the historical overuse of fossil fuels and the current challenges of global climate change, many countries worldwide are undergoing an energy transition [1]. The shift from fossil fuels to clean, renewable energy sources such as wind and solar power has become an inevitable trend [2]. However, the replacement of Synchronous Generators (SGs) by the Inverter-Based Resources (IBRs), e.g., Wind Generators (WGs) and Photovoltaicss (PVs), leads to the decrease of rotational inertia in the power system, and thus, this affects the frequency stability [3]. The blackout in London on 9 August 2019 has drawn considerable attention to the security concerns associated with renewable energy. The primary cause of this incident was the insufficient system inertia due to the high proportion of renewable energy sources. Consequently, the Rate of Change of Frequency (RoCoF) experienced a significant increase post-disturbance, leading to the separation of the system, and the disconnection of both generators and loads from the grid. This event underscores the critical importance of maintaining adequate system inertia in the presence of high renewable energy penetration to ensure grid stability and prevent such cascading failures.

The existing inertia monitoring system in the majority of the power systems, however, only concerns the rotational inertia of the online SGs. Compared with the rotational inertia, the virtual inertia is featured as adjustable and, therefore, can be regulated according to the needs of the system and/or the owner. The rotational and virtual inertia are abbreviated as the synthetic inertia of the power system [4]. This paper aims to develop an accurate online adaptive time window inertia constant identification method for the time-varying synthetic inertia via the measurements currently available to Transmission System Operators (TSOs). The method proposed in this paper can assist in grid dispatch decisions, lay the foundation for an inertia ancillary services market, and provide a reference for grid planning.

1.2. Literature Review

The widespread implementation of PMUs and Wide Area Measurement System (WAMS) has facilitated real-time inertia estimation for TSOs [5]. Several studies have proposed methods for estimating inertia based on different approaches. References [5,6] proposed a method that estimates the inertia of different regional power grids by analyzing changes in line power during low-frequency electromechanical transient oscillations. However, this approach is more suitable for long-distance interconnected power grids and lacks the ability to accurately monitor inertia changes at the POI level, thus failing to capture the spatial distribution characteristics of inertia. Reference [7] estimated inertia by analyzing the frequency evolution process following significant disturbances in the power system. However, this event-driven method faces challenges in monitoring the inertia levels during normal operating conditions. On the other hand, reference [8] utilized the energy function to identify the equivalent inertia constant of individual inertia providers. Although this data-driven approach shows satisfied accuracy for normal operating conditions, it introduces errors in estimating the equivalent inertia constant of time-varying virtual inertia. Reference [9] proposed a system inertia constant tracking method based on recursive algorithm. However, this method cannot take into account both accuracy and tracking speed. When the tracking speed is high, the resulting inertia constant error tends to increase.

To address the adaptive changes of virtual inertia based on the system’s operating state [10], reference [11] introduced an inertia estimation method that tracks the real-time adaptive changes of virtual inertia using the Frequency Deviation (FD) principle of the power system. This principle, proposed by Prof. F. Milano in 2017 [12], forms the basis of the method. Furthermore, references [11,13] extended the work of [11] by proposing adaptive time window inertia constant identification methods for subsystems connected to the power grid. Although these methods can fast approximately converge to the exact inertia constant, the accuracy of the converged inertia constant is much less than the equivalent inertia identification methods proposed in [5,6,7,14] due to the numerical features. In this context, this paper combines the online inertia estimation method with the equivalent inertia constant identification method to accurately track the POI-level synthetic inertia.

In the field of renewable energy grid integration, the inertia associated with renewable sources is emulated through power electronic devices that simulate the inertia of synchronous machines. This inertia can be dynamically adjusted using advanced control methods. The inertia constant identification method proposed in this paper effectively adapts to variations in renewable energy inertia, enabling accurate and timely measurement of inertia at the POI.

1.3. Contributions

Based on the techniques provided in references [11,14], the specific contributions of the paper are the following:

•

An improved equivalent inertia constant identification method by applying the Elastic Net Regression (ENR) technique;

•

An adaptive-time window inertia constant identification method that can identify the variants of the inertia and accurately identify the equivalent inertia constant at the POI.

The feasibility and the accuracy of the proposed method via the ambient measurements of the power system are validated through the tests on a benchmark system.

1.4. Organization

The remainder of the paper is organized as follows. Section 2 briefly describes the inertia providers in modern system. Section 3 describes the online inertia estimation and the improved data-driven equivalent inertia constant identification methods utilized in this paper. The proposed online adaptive time window inertia constant identification method based on the ambient measurements is described in Section 4. The case study is discussed in Section 5 based on the modified WSCC 9-bus system. Finally, conclusions are drawn, and future work is outlined in Section 6.

2. POI-Level Inertia Providers

2.1. POI-Level Equivalent Inertia Constant

For power system, the concept of the inertia is originally defined as the energy that each rotating motor in the system can use to resist system frequency variations and maintain the rated speed of the motor. Furthermore, generally we use the inertia constant H to represent the inertia. The inertia constant is originally defined by the swing equation that describes the dynamics of the rotor speed of a SG arising from the imbalance between the electromagnetic and mechanical torques [15]. The swing equation in per unit form is as follows:

$$\begin{array}{c}\hfill M_{\mathrm{SG}}{\dot{f}}_{\mathrm{SG}}=p_{\mathrm{m},\mathrm{SG}}-p_{\mathrm{e},\mathrm{SG}}-D_{\mathrm{SG}}(f_{\mathrm{SG}}-f_{\mathrm{grid}}),\end{array}$$

(1)

where $f_{\mathrm{SG}}$ is the frequency of the SG and ${\dot{f}}_{\mathrm{SG}}$ denotes its time derivative, which can be abbreviated as RoCoF; $p_{\mathrm{m},\mathrm{SG}}$ is the mechanical power provided by the turbine; $p_{\mathrm{e},\mathrm{SG}}$ is the electrical power that the SG injects to the grid; $D_{\mathrm{SG}}$ is the damping coefficient; $f_{\mathrm{grid}}$ is the frequency of the grid; and $M_{\mathrm{SG}}$ is the mechanical starting time of the SG that links to the inertia constant $H_{\mathrm{SG}}$ of the machine as follows:

$$\begin{array}{c}\hfill M_{\mathrm{SG}}=2H_{\mathrm{SG}}.\end{array}$$

(2)

In order to deal with inertia deficiency, there are many inertia provider in modern power system. They can be divided into three categories: SG, IBR, and VPP.

2.2. Synchronous Generators

In traditional power systems, SGs are the primary providers of inertia, contributing rotational inertia to maintain frequency stability. The rotor of the SGs possesses rotational inertia, storing kinetic energy during normal operation. When a disturbance occurs in the system, the inertia of the rotor enables the SGs to rapidly and passively convert the stored kinetic energy into electrical power through its power-angle characteristics. This kinetic energy is either released into or absorbed by the system, helping to balance the power mismatch and mitigate frequency fluctuations.

During the inertia response of the SG, inertia support power is always equal to the deviation between mechanical power and electromagnetic power. Furthermore, the main role of inertia in this process can be expressed by the swing equation as Equation (1).

2.3. Inverter-Based Resources

IBR refer to a category of Renewable Energy Resource (RES) that are connected to the grid through the power electronics inverter and can convert DC (Direct Current) to AC (Alternating Current) for various purposes. This category primarily includes RES such as photovoltaic systems, wind turbines, and the energy storage system [16]. Renewable energy is connected to the grid through power electronic devices and controlled by MPPT [17], which results in its decoupling from the grid frequency and is unable to provide frequency support and an inertia response for the grid [18]. In order to improve the frequency response capability of the IBR, the converter control strategy can be improved to provide an inertia response. At present, the control method with the best inertia support effect is Virtual Synchronous Generator (VSG) technology.

VSG, also known as synchronous converter, refers to the introduction of the SG rotor swing equation and electromagnetic transient equation in the converter control, so that it has the external characteristics of SG operation in the grid, and can simulate the SG voltage source characteristic. The inertia response of the VSG must have its own energy source. The photovoltaic VSG needs an additional energy storage device, and the wind turbine VSG can rely on the rotational kinetic energy of the wind rotor. The inertia response of renewable energy sources controlled by VSG is similar to the SG. When a contingency occurs in the system, the output power of the unit mutates, resulting in an imbalance between input power and output power. The inertia energy storage unit will absorb or release energy in time to response to the deviation between input and output power, and the system is given instant inertia to maintain frequency security. In this process, the main role of inertia can also be expressed by the swing equation, shown in Equation (1).

2.4. Virtual Power Plant

The VPP concept refers to the aggregation of several devices, including Distributed Energy Resource (DER), Energy Storage System (ESS), and flexible loads, coordinated to operate as a single generating unit. VPPs consist of devices connected to the grid through power electronic converters in most cases [13]. If the resources aggregated by the virtual power plant are equipped properly control method, the VPP also has external characteristics similar to the SG. ESS in virtual power plants and renewable energy using the Grid-forming method can respond to the system disturbances in time and provide inertia support for the power system. Furthermore, within this inertia response time scale, the way the VPP provides inertia can be regarded as consistent with the synchronous machine, and its mechanism of action follows the swing equation.

The virtual power plant model constructed in this paper includes two wind generation, a photovoltaic, an ESS and several loads. Wind generation, photovoltaic, and ESS are equipped with frequency control to make them have inertia response capability. Then, the method of inertia estimation of virtual power plant is studied.

3. Technology Background: Existing Inertia Estimation Techniques

3.1. Online Inertia Estimation

For the sake of derivation of Equation (1), it is convenient to decompose the mechanical power $p_{\mathrm{m},\mathrm{SG}}$ into the following components:

$$\begin{array}{c}\hfill p_{\mathrm{m},\mathrm{SG}}=p_{\mathrm{PFC}}+p_{\mathrm{SFC}}+p_{\mathrm{UC}},\end{array}$$

(3)

where $p_{\mathrm{PFC}}$ and $p_{\mathrm{SFC}}$ are the active power enforced by the Primary Frequency Control (PFC), e.g., Turbine Governor (TG) of the SG, and Secondary Frequency Control (SFC), e.g., Automatic Generation Control (AGC); $p_{\mathrm{UC}}$ is the power set point allocated by solving the Unit Commitment (UC) problem.

Merging (1) and (3) and differentiating with respect to time, one has:

$$\begin{array}{cc}\hfill M_{\mathrm{SG}}{\ddot{f}}_{\mathrm{SG}}=& {\dot{p}}_{\mathrm{PFC}}+{\dot{p}}_{\mathrm{SFC}}+{\dot{p}}_{\mathrm{UC}}-{\dot{p}}_{\mathrm{e},\mathrm{SG}}\hfill \\ & -D_{\mathrm{SG}}({\dot{f}}_{\mathrm{SG}}-{\dot{f}}_{\mathrm{grid}}).\hfill \end{array}$$

(4)

Considering the time scale of the inertial response of the SG, we can assume that ${\dot{p}}_{\mathrm{PFC}}\approx 0$, ${\dot{p}}_{\mathrm{SFC}}\approx 0$ and ${\dot{p}}_{\mathrm{UC}}\approx 0$. Since $D_{\mathrm{SG}}\ll M_{\mathrm{SG}}$ always holds for the SGs, assuming, in addition, that $D_{\mathrm{SG}}\approx 0$. In this vein, one deduces a simple inertia estimation formula:

$$\begin{array}{c}\hfill H_{\mathrm{SG}}={\displaystyle \frac{{\dot{p}}_{\mathrm{SG}}}{2{\ddot{f}}_{\mathrm{SG}}}}.\end{array}$$

(5)

Equation (5) also works well for IBRs, under the assumption that they are controlled to provide a dynamic response in the same time scale as SGs’ inertial response [19]. In this vein, the subscript in (6) can be removed to show the generality:

$$\begin{array}{c}\hfill H={\displaystyle \frac{\dot{p}}{2\ddot{f}}}.\end{array}$$

(6)

In this case, $\dot{p}$ is the Rate of Change of Power (RoCoP) of the inertia-providing device connected to the grid and $\ddot{f}$ is the second derivative of its frequency. Note that $\dot{p}$ can be estimated based on PMU measurements, and $\ddot{f}$ can be estimated through the techniques provided in [12]. Equation (6) suffers a numerical issue, i.e., the denominator $\ddot{f}$ will be infinitesimal in quasi-steady-state conditions, which results in a severe error in the estimation of H. To cope with this issue, reference [11] proposed an online inertia estimation method by introducing a differential technique to avoid the numerical issue resulting from the division:

$$\begin{array}{c}\hfill T_H{\dot{H}}^{*}=\gamma \left(\ddot{f}\right)(\dot{p}+2H^{*}\ddot{f}),\end{array}$$

(7)

where the superscript * denotes the estimated result and $T_H$ is the time constant for the inertia estimation. The selection of $T_H$ is a trade-off between tracking speed and the numerical stability of the estimator. Furthermore, $\gamma (·)$ in (7) is a function defined as follows:

$$\begin{array}{c}\hfill \gamma \left(x\right)=\left\{\begin{array}{cc}-1,\hfill & x\ge {ϵ}_x,\hfill \\ 0,\hfill & -{ϵ}_x<x<{ϵ}_x,\hfill \\ 1,\hfill & x\le -{ϵ}_x,\hfill \end{array}\right.\end{array}$$

(8)

where ${ϵ}_x$ is a small positive threshold closing to zero aiming at mitigating the impact of measurement noise. This inertia estimation was originally proposed for tracking the inertia of an individual generation during a frequency event. It is expected to converge to the accurate inertia constant at post-contingency status, i.e., $\ddot{f}<{ϵ}_{\ddot{f}}$.

Compared to Formula (6), the online inertia estimation Formula (7) improves the numerical stability while largely decreasing the convergence speed. Therefore, the accuracy of the online inertia estimation is affected by the frequency regulations of PFC and damping. To cope with this issue, Reference [11] extends (7) into the following form:

$$\begin{array}{cc}\hfill T_H{\dot{H}}^{*}& =\gamma \left(\ddot{f}\right)(\dot{p}+2H^{*}\ddot{f}+K_D\dot{f}),\hfill \\ \hfill T_D{\dot{K}}_D^{*}& =\gamma (\Delta f)(\int \dot{p}dt+2H^{*}\dot{f}+K_D^{*}\Delta f),\hfill \end{array}$$

(9)

where $\Delta f=f-f_{\mathrm{ref}}$ and $f_{\mathrm{grid}}\approx f_{\mathrm{ref}}$ holds for normal operating conditions and $K_D^{*}$ is an estimated variable aiming at evaluating the mixing effect of damping and PFC during the inertia estimation period to improve the accuracy of the inertia estimation. As a consequence, (9) can estimate the inertia of a generation device in both normal operating conditions and following a frequency event. Note that since the input of the inertia estimation is implemented in Formula (9), namely the frequency and output power, the result in $H^{*}$ can be a time-varying trajectory converging to the equivalent inertia constant as shown in reference [11].

3.2. Equivalent Inertia Constant Identification Method

The online inertia estimation method introduced in Section 3.1 can track the time-varying inertia. The real-time tracking ability, however, leads to a relatively inferior ability to identify the equivalent inertia constant for a certain inertia level. With this regard, reference [14] provided a data-driven regression-based method to further improve the inertia estimation accuracy while sacrificing the ability of real-time tracking. For the measurement data obtained at time $t_i$, one has:

$$\begin{array}{c}\hfill H\approx {\overline{H}}^{*}={\displaystyle \frac{|\dot{p}\left(t_i\right)+{\xi }_p|}{2|\ddot{f}\left(t_i\right)+{\xi }_f|}}={\displaystyle \frac{|{\dot{p}}^{*}\left(t_i\right)|}{2|{\ddot{f}}^{*}\left(t_i\right)|}},\end{array}$$

(10)

where the absolute values of ${\dot{p}}^{*}$ and ${\ddot{f}}^{*}$ are considered for simplicity without losing information and $\xi$ presents the measurement errors and noises generated from PMUs.

It is relevant to rewrite (10) to solve the numerical issue; thus, one has:

$$\begin{array}{c}\hfill 2{\overline{H}}^{*}|{\ddot{f}}^{*}\left(t_i\right)|-|{\dot{p}}^{*}\left(t_i\right)|=0.\end{array}$$

(11)

One can set up the following two vectors with a set of measurements at time $t_1$, $t_2$, …., $t_{N-1}$, $t_N$:

$$\begin{array}{cc}\hfill \mathit{A}& =\left[\right|{\ddot{f}}^{*}\left(t_1\right)|,|{\ddot{f}}^{*}\left(t_2\right)|,...,|{\ddot{f}}^{*}\left(t_{N-1}\right)|,|{\ddot{f}}^{*}\left(t_N\right){\left|\right]}^{\mathrm{T}},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathit{b}& =\left[\right|{\dot{p}}^{*}\left(t_1\right)|,|{\dot{p}}^{*}\left(t_2\right)|,...,|{\dot{p}}^{*}\left(t_{N-1}\right)|,|{\dot{p}}^{*}\left(t_N\right){\left|\right]}^{\mathrm{T}},\hfill \end{array}$$

(13)

where $\mathrm{T}$ denotes the matrix transposition and $\mathit{A}$ and $\mathit{b}$ present the matrices that consist of the measured data ${\ddot{f}}^{*}$ and ${\dot{p}}^{*}$, respectively. In this vein, the estimation of ${\overline{H}}^{*}$ becomes a Least Square Method (LSM) problem [14]:

$$\begin{array}{c}\hfill \mathrm{Minimize}{(\mathit{A}{\overline{\mathit{H}}}^{*}-\mathit{b})}^2.\end{array}$$

(14)

which admits the well-known solution:

$$\begin{array}{c}\hfill {\overline{H}}^{*}={\left({\mathit{A}}^{\mathrm{T}}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{b}.\end{array}$$

(15)

The LSM is expected to return an “average” value of the equivalent inertia constant within a given time period. This technique, therefore, is prone to over-fitting when the data contain significant noise or has a short sampling period [14].

To address the over-fitting and under-fitting issue in LSM, this paper utilizes the ENR regularization algorithm, which combines the advantages of ridge regularization and lasso regularization, to transform (14) with regularization terms:

$$\begin{array}{c}\hfill L_{\mathrm{T}}\left({\overline{H}}^{*}\right)=\Vert \mathit{A}{\overline{H}}^{*}{-\mathit{b}\Vert }^2+\lambda \rho \Vert {\overline{H}}^{*}{\Vert }_1+{\displaystyle \frac{\lambda (1-\rho )}{2}}{\Vert {\overline{H}}^{*}\Vert }_2^2,\end{array}$$

(16)

where $\lambda$ and $\rho$ are the coefficients of the regular terms, respectively. If $\lambda =0$, the Equation (16) is identical with the least square regression (15). When $\rho =0$, Equation (16) becomes a ridge regression, while it is a Lasso regression when $\rho =1$. When over-fitting occurs, the data that lead to over-fitting can be filtered by adjusting these coefficients.

To minimize the error of the estimated equivalent inertia constant, the following solution is adopted to obtain the optimal result:

$$\begin{array}{c}\hfill {\overline{H}}^{*}=\underset{{\overline{H}}^{*}}{argmin}{L}_{\mathrm{T}}\left({\overline{H}}^{*}\right).\end{array}$$

(17)

where argmin means to minimize the argument. A brief discussion about the key factors to the accuracy and reliability of improved equivalent inertia constant identification method is provided as follows:

• Threshold for input variables: The evolution from (14) to (16) mitigates the over-fitting and under-fitting in LSM, as the ENR regularization algorithm utilized in Section 3.2 can effectively address the over-fitting and under-fitting issue of LSM via regulating $\lambda$ and $\rho$. However, the extreme data pair with significant errors may still affect the accuracy of $H^{*}$. To cope with this issue, the proposed method is fed by the estimated values ${\ddot{f}}^{*}$ and ${\dot{p}}^{*}$ through the PI filters utilized in [11], and the pairs $(|{\ddot{f}}^{*}|,|{\dot{p}}^{*}|)$ with $|{\ddot{f}}^{*}|<{10}^{-5}$ are removed directly.
• Selection of time window: the length of the time window $\tau =t_N-t_1$ presents the amount of data fed to the LSM problem, and thus, this affects the accuracy of (17). The time window located in normal operation implies that the measurement noise may affect the accuracy of the inertia estimation as it is relatively large compared with the RoCoP. On the other hand, the time window located in the evaluation period following a contingency may be impacted by PFC for its effect on boosting the active power outputs at POI responding to the frequency dynamics.
• Change of inertia: The LSM method is expected to return an “average” value of the equivalent inertia constant within a given time window, which cannot identify the occurrence of the change of inertia constant of the device, e.g., the regulation of the control strategies of a VSG and the change of connection status of the synchronous generation units.

4. Proposed POI-Level Online Inertia Evaluation Method

The online inertia estimation method reviewed in Section 3.1 can effectively tracking the real-time inertia but may introduce errors and numerical issues for the converged equivalent inertia constant. While the method proposed in Section 3.2 can accurately deduce the equivalent inertia constant, it cannot identify the time-varying inertia, and thus, it may lead to a fatal error if the equivalent inertia constant is changed within the identification time window, which can occur for the POI-level synthetic inertia. This section proposes an online adaptive time window inertia constant identification method to take the advantages of the above techniques to effectively estimate the POI-level synthetic inertia via the ambient measurements with the ambient data smoothing and inertia variation identification techniques.

4.1. Ambient Data Smoothing

To cope with the numerical errors existed in the trajectory of the estimated inertia obtained by the technique introduced in Section 3, we use the Weighted Least Squares (WLS) method to smooth the trajectory.

The WLS method refines the data smoothing process and adapts it to the characteristics of the data by introducing the concept of weights:

$$\begin{array}{c}\hfill \underset{\mathit{\theta }}{min}\sum _i^n{w}_i{(H^{*}\left(t_i\right)-f(t_i;\mathit{\theta }))}^2,\end{array}$$

(18)

where $w_i$ is the weight of each data; $H^{*}\left(t_i\right)$ is the estimated inertia constant at time $t_i$; $f(t_i;\mathit{\theta })$ is the fitted function; and $\mathit{\theta }$ is the parameter vector of fitted function.

In this paper, the fitting function $f(t_i;\mathit{\theta })$ at a given time period is set to:

$$\begin{array}{c}\hfill {\mathit{H}}_{\mathrm{WLS}}^{*}=k_{\mathrm{WLS}}\mathit{t}+a,\end{array}$$

(19)

where $k_{\mathrm{WLS}}$ and a are the values of $\mathit{\theta }$, which can be obtained by the WLS method:

$$\begin{array}{cc}\hfill k_{\mathrm{WLS}}& ={\left({\mathit{t}}^{\mathrm{T}}\mathit{W}\mathit{t}\right)}^{-1}{\mathit{t}}^{\mathrm{T}}\mathit{W}{\mathit{H}}_{\mathrm{WLS}}^{*},\hfill \\ \hfill a& ={\mathit{H}}_{\mathrm{WLS}}^{*}-k_{\mathrm{WLS}}\mathit{t},\hfill \end{array}$$

(20)

where $\mathit{W}$ is the weight vector. Then, the parameters can be updated based on the residual function $\mathit{R}$ and the Jacobian matrix $\mathit{J}$ until the residual function is relevantly small. One has:

$$\begin{array}{cc}\hfill {\theta }_{k+1}& ={\theta }_k-{({\mathit{J}}^{\mathrm{T}}\mathit{J}+\lambda \mathit{I})}^{-1}{\mathit{J}}^{\mathrm{T}}\mathit{R},\hfill \\ \hfill \mathit{R}& =\mid {\mathit{H}}^{*}-{\mathit{H}}_{\mathrm{WLS}}^{*}\mid ,\hfill \end{array}$$

(21)

where $\lambda$ controls the learning rate during the iterative process.

4.2. Inertia Variation Identification

The online inertia estimation method (9) can track the real-time inertia, and thus, it can detect the variation of the inertia level due to the tuning of virtual inertia or the connection status of the SGs. However, the stochastic perturbations of RESs may introduce numerical fluctuations of the $H^{*}$, which may affect the identification of the inertia variation. In this vein, this subsection proposes a sliding time window method to identify the variation of the inertia constant to avoid the perturbations, as follows:

$$\begin{array}{c}\hfill R_{\mathrm{SC}}={\displaystyle \frac{\Delta S_{H^{*}}}{C_{H^{*}}}},\end{array}$$

(22)

where $\Delta S_{H^{*}}$ is the difference between adjacent time windows; $S_{H^{*}}$ is area of the inertia window within a sliding period $T_{\mathrm{SC}}$; $C_{H^{*}}$ is the perimeter of the inertia window; $R_{\mathrm{SC}}$ is the area-to-perimeter ratio. $S_{H^{*}}$ and $C_{H^{*}}$ are calculated as follows:

$$\begin{array}{cc}\hfill S_{H^{*}}& ={\int }_t^{t_m}{H}^{*}\left(t_n\right)\mathrm{d}t,\hfill \\ \hfill C_{H^{*}}& =\sum _{i=1}^m\sqrt{{(\Delta t)}^2+{(H^{*}\left(t_i\right)-H^{*}\left(t_{i-1}\right))}^2}\hfill \\ \hfill & +H^{*}\left(t\right)+H^{*}\left(t_m\right)+T_{\mathrm{SC}},\hfill \\ \hfill t_m& =t+T_{\mathrm{SC}}=t+m\times \Delta t,\hfill \end{array}$$

(23)

where $\Delta t$ is the sampling rate of the input variable, namely the sampling rate of the ambient measurement in this paper. If $R_{\mathrm{SC}}\le {ϵ}_{\mathrm{SC}}$, one can determine that the inertia level is changed. The selection of ${ϵ}_{\mathrm{SC}}$ should consider the historical datum of the RESs. Through the tests based on the IEC standard stochastic models, we select $T_{\mathrm{SC}}=1$ s and ${ϵ}_{\mathrm{SC}}=0.003$ to obtain a relatively satisfied accuracy.

4.3. Framework of the Online Adaptive Time Window Inertia Constant Identification Method

This subsection proposes an online adaptive time window inertia constant identification method by combining the methods in Section 3.1 and Section 3.2. The diagram of the online adaptive time window inertia constant identification is shown in Figure 1.

The detailed working principles of the system are as follows.

1.

Obtain the active power output and frequency of the synthetic inertia provider at the POI through PMU and estimates $\dot{p}$ and $\ddot{f}$ through the techniques provided in [12].

2.

Utilize the online adaptive time window inertia constant identification method provided in [11] to track the real-time inertia constant of the devices.

3.

Smooth the inertia constant data based on WLS algorithm proposed in Section 4.1.

4.

Judge if the inertia varied according to the technique proposed in Section 4.2.

5.

Once the inertia variation is detected, select a suitable time window following the variation and filter out abnormal values and noise in the frequency and active power output data within this time window.

6.

Finally, the improved inertia constant identification method proposed in Section 3.2 is applied to obtain an accurate equivalent inertia constant at this period.

5. Case Study

In this section, simulation model of a grid in northwest China is served to validate the proposed techniques. The simplified topology of the system is shown in Figure 2. The grid in northwest China consists of 18 buses at 220 kV and 4 buses at 750 kV. The grid can be divided into six zones. Furthermore, it includes nine thermal power plants, six hydropower plants, and four renewable energy plants. The capacity of each generator is shown in

Table 1. The capacity of the generator.

Generator	Capacity [MW]	Generator	Capacity [MW]	Generator	Capacity [MW]
A1	1000	A8	508	B6	540
A2	650	A9	560	C1	632
A3	250	B1	1000	C2	520
A4	830	B2	540	C3	450
A5	508	B3	560	C4	650
A6	650	B4	650
A7	540	B5	650

. Furthermore, all generators are running. The total load in the grid is 10.88 KMW. The RESs penetration is set as 20.7% to simulate the low-inertia system, which is contributed by the renewable energy power stations. In Equation (9), $T_{\mathrm{H}}$ is set to 0.004 and $T_{\mathrm{D}}$ is set to 0.001.

For the system shown in Figure 2, the SG of A4 with $H_{\mathrm{SG}}=6.02\mathrm{MWs}/\mathrm{MVA}$, $D_{\mathrm{SG}}=1.0\mathrm{MWs}/\mathrm{MVA}$, equipped with a TG for PFC, or a VSG with $H_{\mathrm{VSG}}=40.0\mathrm{MWs}/\mathrm{MVA}$, $D_{\mathrm{VSG}}=10\mathrm{MWs}/\mathrm{MVA}$ is connected to D8 for different scenarios. The renewable energy power stations are modeled as VSG. The detailed model of the VSG is provided in [20]. Besides, a VPP is connected to the D14. The VPP includes two wind turbines, a photovoltaic, an ESS, and several loads. The parameters of the VPP are detailed in [13]. We tune the frequency response parameters of the machines in VPP to simulate the inertia variations at POI in reality.

The bus frequency and the active power output of the inertia provider connected to D16 and bus D14 are measured with PMUs with a sampling rate of 100 Hz. The measurement noise of the PMUs is modeled as an Ornstein-Uhlenbeck stochastic process [21]. The uncertainty of wind turbine and photovoltaic output is taken into account when modeling. The dynamic wind speed of the WPP is modeled as a Weibull distribution with exponentially decaying auto-correlation [22]. All simulations are obtained using the Python-based software tool DOME [23]. The ENR algorithm for inertia constant identification is solved by using Python 3.8.10 and Numpy 1.20.0.

5.1. Identification of Equivalent Inertia Constant

Considering that the inertia constant of a SG is fixed and usually known by TSOs, this subsection adopts the estimated result of the equivalent inertia constant for the SG connected to D8 in the grid of northwest China as a reference to select an appropriate $\lambda$ and $\rho$ of (16). Setting appropriate parameters $\lambda$ and $\rho$ can avoid the phenomena of overfitting and underfitting in the algorithm. The influence of different $\lambda$ and $\rho$ values on the estimation results of the equivalent inertia constant for the SG within a 10 s measurement time window is shown in Figure 3.

Firstly, by setting $\rho =0$, the relative error for inertia constant identification is shown in Figure 3a for different $\lambda$ and the optimal selection is $\lambda =1.4\times {10}^{-6}$ with the relative error equal to $0.04\%$. With $\lambda =1.4\times {10}^{-6}$, the relative errors for different $\rho$ are shown in Figure 3b and the optimal selection is $\rho =0.04$ with the relative error equal to $0.01\%$. Also note that when $\lambda =0$, the ENR method (see (16)) becomes LSM method (see (15)), and at this time, the relative error for the inertia constant identification is $0.93\%$, which is much higher than the best performance of the enr. The selection of parameters a and b will affect the accuracy of the ENR algorithm, as indicated by the experimental results. Furthermore, this result further proves the advantages of the ENR algorithm over the LSM.

5.2. Real-Time Inertia Tracking

In this subsection, we validate the functions of the online adaptive time window inertia constant identification method for real-time inertia tracking and identifying the change of inertia of the VSG in the system shown in Figure 2.

The real-time inertia tracking results of the VSG obtained from the online adaptive time window inertia constant identification method are shown in Figure 4. The orange dot line represents the actual equivalent inertia constant for the VSG at D16, which is changed twice at 50 s (from 40 MWs/MVA to 30 MWs/MVA) and 130 s (from 30 MWs/MVA to 50 MWs/MVA). In order to accurately identify the moments of inertia transitions, the WLS algorithm is applied to deal with the estimated inertia (shown as a blue dash line in Figure 4), which effectively smoothed the numerical fluctuations that existed in the original estimated results, shown as the tomato red line.

The technique proposed in Section 4.2 is used to further process the smoothed data by selecting different time intervals in order to detect the change of the inertia. The resulting trajectory is shown in Figure 5. The blue line is the resulting trajectory. The dashed line is the time window when $R_{\mathrm{SC}}$ exceeds the threshold. And the dotted red line indicates the inertia change time. With $T_{\mathrm{SC}}=1$ s and ${ϵ}_{\mathrm{SC}}=0.003$, the system detects two inertia variation periods, namely 44 s to 60 s and 123 s to 140 s. These periods include the actual moments that the equivalent inertia constant of the VSG at D16 changed. This result indicates that the online adaptive time window inertia constant identification method can effectively distinguish the periods for inertia variants and avoids locating the time window for equivalent inertia constant identification at these periods.

5.3. Accuracy of Equivalent Inertia Constant Identification

This subsection aims to investigate the effects of different time windows on the accuracy of equivalent inertia constant identification proposed in Section 4.2. Figure 6 shows the identification results for the VSG at D16 ($H=40$ MWs/MVA) and the SG ($H=6.02$ MWs/MVA) with different time windows. In Figure 6a, all the time windows are located in the after contingency condition, namely the $t_0$ is selected as $0.1$ s later than a sudden load increase. In Figure 6b, all the time windows are located in the normal operation with the stochastic wind as the continuous small disturbances.

According to Figure 6, the proposed method can accurately identify the equivalent inertia constant of the VSG with relative error less than $4\%$ under different conditions with a time window longer than 6 s, and the involvement of the after-contingency evolution can obviously improve the accuracy and shorten the requirement of the time window, as shown in Figure 6a. These results indicate that the proposed method can cope with the impact of PFC during the contingency evolution, while the measurement noise can affect the accuracy of the method in normal operating conditions.

5.4. Inertia Estimation of IBR

This subsection evaluates the overall performance of the proposed online adaptive time window inertia constant identification system. To this aim, we tune the inertia constant of the VSG in the system to simulate the inertia variations in modern power systems. Consider the same scenario with Section 5.3, and set the length of the time window for equivalent inertia identification as 10 s. The proposed online adaptive time window inertia constant identification system automatically picks three time windows to identify the equivalent inertia constant at the POI as shown in Figure 7. In addition to the initial setting of the inertia constants in Section 5, this section modifies the inertia constant values and adds a set of experiments. The corresponding equivalent inertia constant identification results ${\overline{H}}_{\mathrm{VSG}}^{*}$ are shown in

Table 2. The inertia of the IBR.

	Scenario 1			Scenario 2
Time Window	i	ii	iii	i	ii	iii
$H_{\mathrm{VSG}}$[MWs/MVA]	40	30	50	25	32	20
${\overline{H}}_{\mathrm{VSG}}^{*}$[MWs/MVA]	39.18	29.12	49.27	24.37	31.22	19.48
$\overline{\sigma }$[%]	2.05	2.93	1.46	2.52	2.43	2.6

.

As shown in Figure 7, the errors of the inertia constants measured in the three time windows in scenario 1 are 2.05%, 2.93%, and 1.46%, respectively. The errors of the inertia constants measured in the three time windows in scenario 2 are 2.52%, 2.43%, and 2.6%. The results indicate that the proposed online synthetic inertia estimation technique can accurately identify the equivalent inertia constant of the VSG with variant inertia under the normal operating conditions of the power system with high RES penetration.

5.5. Inertia Estimation of VPP

For the system shown in Figure 2, the real-time inertia tracking results of the VPP obtained from the proposed online adaptive time window inertia constant identification system are shown in Figure 8. The blue dashed line represents the actual equivalent inertia constant for the VPP, which is changed at 32 s. The orange line represents the smoothed inertia constant through the WLS algorithm. Then, we use the technique which proposed in Section 4 to detect the change of the inertia. The resulting trajectory is shown in Figure 9. The blue line is the resulting trajectory of VPP. The dashed line is the time window when $R_{\mathrm{SC}}$ exceeds the threshold. And the dotted red line indicates the inertia change time. With $T_{\mathrm{SC}}=1$ s and ${ϵ}_{\mathrm{SC}}=0.05$, the system detects the inertia variation time is 31.3 s.

Consider the scenario with Section 5.4 and set the length of the time window as 10 s. We pick two time windows to identify the equivalent inertia constant at the POI, as shown in Figure 8. In this section, the system’s RES penetration is adjusted by modifying the power of renewable energy units within the grid and the VPP. Two scenarios (22.9% res penetration and 26.5% res penetration) are set up to study the impact of renewable energy penetration on the proposed method. The corresponding equivalent inertia constant identification results ${\overline{H}}_{\mathrm{VPP}}^{*}$ are shown in

Table 3. The inertia of the VPP.

Time Window	i		ii
RES penetration [%]	22.9		26.5
$H_{\mathrm{VPP}}$ [MWs/MVA]	11	3	11	3
${\overline{H}}_{\mathrm{VPP}}^{*}$ [MWs/MVA]	10.60	2.88	10.55	2.90
$\overline{\sigma }$[%]	3.6	4.0	4.1	3.3

.

As shown in Figure 9, the inertia of the virtual power plant measured using the method proposed in this paper is 10.6 MWs/MVA before the change and 2.88 MWs/MVA after the change. The deviations from the actual values are 3.6% and 4.0% when the res penetration reaches 22.9%. The deviations from the actual values are 4.1% and 3.3% when the res penetration reaches 26.5%. The results show that the res penetration have minimal impact. In summary, the proposed online synthetic inertia estimation technique can accurately identify the equivalent inertia constant of the VPP with variant inertia. Besides, we make a comparison to compare the advantages of this work with the latest research in

Table 4. A comparison between the proposed method and the latest research.

Criteria	This Paper Work	Latest Research
Methodology	The proposed method	[8,9]
Accuracy	High accuracy	High accuracy
Real-time	can quickly track the time-varying virtual inertia	[8] no; [9] can quickly track the time-varying virtual inertia but the accuracy will decrease
Computational Efficiency	fast	depend on the data size
Research Subject	SGs, IBRs and VPPs	SGs, and IBRs

.

6. Conclusions

This paper proposes an online adaptive time window inertia constant identification method based on ambient measurements to accurately identify the equivalent inertia constant for the time-varying synthetic inertia at the POI level. This method can be applied to different inertia providers, such as SGs, IBRs, and VPPs. We have verified through SG and VSG experiments that ENR technology can significantly improve the accuracy of the data-driven equivalent inertia constant identification method. Next, based on the method of determine whether the inertia level is changed, the online adaptive time window inertia constant identification method which combine both the real-time inertia tracking technique and the data-driven equivalent inertia constant identification with the inclusion of ambient data smoothing technique is proposed. Simulation results support the proposed online adaptive time window inertia constant identification method by validating its accuracy in tracking variant synthetic inertia. The results indicate that the 10 s window is suitable for inertia constant identification, and the method proposed in this paper has a high success rate, with an error margin of less than 5%.

As a next step, we will aim to collaborate with State Grid Corporation of China or renewable energy companies to test its robustness and accuracy of the proposed method using real data, with the goal of implementing it in actual grids. Additionally, future work will also focus on implementing advanced filtering techniques such as Kalman filtering, to further increase the accuracy of the proposed techniques.