Data-driven resilience analysis of power grids

Nomenclature

?

0 Introduction

As the penetration of Renewable Energy Sources (RES)in systems is on an increase,power grid resilience has become a major concern in power engineering.Resilience is the ability of a system to ride through with a disturbance or return to a pre-disturbance state before a system collapses.A critical element of a power system’s resilience is system inertia.Power system inertia represents the resistance that the power grid opposes to changes in its frequency by means of the kinetic energy stored in its rotating masses[1].Although a large number of traditional synchronous generators have been replaced by RES to protect the environment,the intermittent and fluctuating power output of RES has exposed power systems to possible stability issues.Most RES systems do not have rotating equipment with mechanical inertia.RES power electronic systems significantly reduce the system inertia,making it sensitive to disturbances or large transients (load steps).

Reliability issues caused by RES integration may lead to blackouts,thereby negatively impacting electricity users.For example,in 2019,a massive power outage hit Britain,affecting more than one million people.Investigations have shown that this was caused by an increase in the number of wind turbines,greatly reducing the power system inertia and thus decreasing the ability of the system to withstand disturbances and large load transients.Therefore,to evaluate the stability of the grids and implement proactive control before incidents (e.g.,blackout),estimating the inertia of the power system is of great significance.

Extensive studies have been conducted on the inertia of power systems.Conventional inertia is related to the energy restored in the rotating mass,while the inertia of the power system counteracts the rate of change of frequency (ROCOF)after disturbance.The inertia of the entire grid also considers the frequency and voltage dependence of the load and the power-electronic interfaces for generation,load,and storage[2].Existing methods associated with inertia estimation can be classified into three groups:methods based on the swing equation,defamation of swing equation,and trained specific models,which are all summarized below.

According to the swing equation,[3]presents a robust estimation of the generator inertia during disturbance; [4]provides an online algorithm for simultaneous estimation;[5]calculates the inertia offline and selects appropriate time ranges of ROCOF; [2]proposes an algorithm that allows for estimating inertia from ambient frequency and active power signals; [6]utilizes the mean value of generator active power change to estimate a single generator; in [7],a voltage dependency method is proposed to modify the frequency response after disturbance.For the deformed equation,in [8],a polynomial approximation was implemented to fit the oscillatory of frequency change; in [9],the frequency response from a single location was used to estimate the inertia of the Western Electricity Coordination Council (WECC).In [10],phasor measurement units(PMUs) are employed in the power system in Great Britain to provide accurate data for estimation.In [11],a statistical model is formulated as a Gaussian mixture model with temporal dependence encoded as Markov chains to estimate inertia.These methods determine an application time range that reduces the error of the first sample after disturbance.Online methods present a real-time estimation for a power system.However,the major limitation of these existing methods is that the existing methods cannot directly showcase the transient process.Without the trajectory of transient dynamics,determining whether transient states will surpass the resistance of the power system is difficult.

To overcome the limitations of the existing methods,a disturbance-based method is presented herein to estimate inertia,and a Gaussian process regression method is then implemented to predict the trajectory of the power system after disturbance.The main contributions of this study can be summarized as follows:

• An innovative data processing method is proposed for inertia estimation by combining a moving average filter and a median filter.

• R-method is applied to solve the disturbance-based inertia estimation problem.

• Gaussian progress regression is implemented to predict the trajectory after disturbance,focusing on the dynamics of the transient process after disturbance instead of the steady state.

The remainder of this paper is organized as follows.Section 1 introduces the disturbance-based method for estimating power grid inertia.Section 2 presents the Gaussian progress regression method for predicting the power system trajectory after disturbance.A test system is provided in Section 3 to validate the feasibility of our proposed method in estimating system inertia and predicting system dynamic trajectories.Finally,the concluding statements are provided in Section 4.

1 Disturbance-based Method for Estimating Inertia

Based on the swing equation,a disturbance-enabled estimation method is presented to predict the power grid inertia.When a total inertia constant and an average frequency for the whole system are considered (by transforming the system to the Center Of Inertia (COI)reference frame) and the mechanical and electrical powers of the generators are summed up,the swing equation can be employed to describe the electromechanical dynamics of a multi-machine power system [12]:

where H denotes the inertia of the power system in seconds,df(t)/dt is the ROCOF,ΔPm(t) is the mechanical power change,ΔPe(t) is the electrical power change,and ΔP(t) is the total power change of the entire system.

The mechanical power change ΔPm(t) results from the governors’ response to the disturbance.To prevent the imbalance of generation and consumption in a power system,the governors respond to the frequency change,thereby impacting the mechanical and total power changes of the system.

The electrical power change ΔPe(t) is equal to the disturbance and load power changes,which can be defined as

where ΔPdist is the disturbance and ΔPL(t) is the change in load demand.The sign of ΔPdist is assumed to be positive when there is a power deficiency (e.g.,generation loss)caused by disturbance and negative when there is a power surplus (e.g.,disconnection of loads).

In this paper,the total power change is considered as follows:

Here,we examined the behavior of the governor to approach h1(f(t)),a function to describe the relation between the mechanical power change and frequency deviation.In the steady state,the stationary mechanical power change and the frequency deviation have the following relation:

where ΔPm is the mechanical power change,R is the system gain,and Δf is the frequency deviation of the steady state.However,(4) merely describes the relation between the mechanical power change and the frequency change in the steady state.Therefore,in this study,R(t) is used to replace R to accommodate the dynamic changes in frequency.Then,h1(f(t)) can be defined as R(t)Δf(t),and in this case,(1)can be modified as

where Hest is the estimated inertia needed to be obtained.

In (5),df(t)/dt,Δf(t),and ΔPdist are considered to be known,while Hest and R(t) are unknown and required to be calculated.

Before calculating df(t)/dt,Δf(t),and ΔPdist,the data obtained from the measurements must be processed to reduce the noise.In this paper,we present an innovative method that combines a median filter and a moving average filter.

The median filter is a nonlinear filter [13]that has proved to be effective in satisfying the dual requirements of removing impulse noise while preserving rapid signal changes [14,15].L is assumed to be the length of a series of continuous points in the curve.In range L,the median value is selected as the sample value.

The moving average filter establishes a data buffer to store N sampled data sequentially,and the first collected data is discarded after each new data is collected;subsequently,the arithmetic mean or weighted average value (herein,the arithmetic mean is adopted) of N data,including the new data,is obtained.For each sampling,the data block moves once to obtain the sum of the new set of data as well as the average [16].After filtering,the data is much smoother,especially when interference or impulse occurs [17].

In this paper,the median filter is implemented first,and the results are processed by moving the average filter.After filtering the data,df(t)/dt,Δf(t),and ΔPdist can be utilized in the calculation.

To estimate the value of R(t) at certain time tsr,N discrete points around tsr in the ROCOF curve are examined and employed in (5).Therefore,N + 1 unknowns and N + 1 linear equations are formed:

The values of R(t) for the N + 1 selected points can be computed by solving these N + 1 linear equations.The inertia of the power system can be calculated as

The frequency and ΔPdist are obtained from the system,while the application time tsr and range are selected according to the rules below:

• tsr is the first local extreme point of the ROCOF curve after disturbance time tdist.

• The range can be obtained around the first local extreme point (either maximum or minimum) of the ROCOF curve.The application range is N points symmetrically distributed around tsr (i.e.,the range extends from tdist to 2tsr - tdist).

2 Gaussian Process Regression for Predicting Dynamics

After determining the power grid inertia,knowing its dynamics is still necessary.In this section,Gaussian process regression (GPR) is employed to predict the trajectory after a disturbance.GPR models are nonparametric kernel-based probabilistic models that use Gaussian process priors for data regression analysis [18-22].

The GPR model includes noise and a Gaussian process a priori.Its solution is based on Bayesian inference.Without limiting the form of the kernel function,GPR is theoretically a general approximation for any continuous function in the compact space.In addition,GPR can provide the predicted results a posteriori,which has an analytical form when the likelihood is normally distributed.Therefore,GPR is a generalized and analytic probabilistic data-driven method.

We herein present a GPR model from the viewpoint of weight space and function space.The former is based on Bayesian linear regression,and the GPR prediction form is obtained through the mapping of feature space; the latter is based on the Gaussian process,and its equivalent result is obtained from the edge distribution property of the normal distribution.

2.1 The view of weight-space

GPR can be deduced via Bayesian linear regression with a normal distribution [23]-[25].Imagine that we have n data groups:X={X1,X2,…,Xn},y={y1,y2,…,yn}.Assume linear regression is used to describe the relation between outputs and inputs:

where ω is the weight coefficient,∈ is the noise term following a normal distribution, and its Gaussian prior is

We use a Gaussian prior over the weightsThus,the likelihood of Bayesian linear regression is

Subsequently,based on Bayes’ theorem,the posterior is

where

Now,we intend to predict the outputs at the new introduced point X*.The idea is to average out the error and focus on the expected value provided by the function f,through which we can predict f*

Subsequently,the posterior can be integrated to obtain the predict output:

As the regression is chosen to be linear,we assume the function as a linear one.However,in practice,the linear relation is minor.To predict for the majority of nonlinear functions,we can transform the inputs with a nonlinear function Φ(X),X can be mapped into a high-dimensional space:

As Φ(X) is fixed,it has no relation with the weight.Hence,we can directly substitute (15)(16) into the result of Bayesian linear regression and obtain the following result:

where

Then,we utilize the kernel method to transform (8) and obtain the following equation.

where

where k is the kernel function,K is the Gram matrix that equals

2.2 The view of function-space

Imagine we have an observation Dt ={ X t,yt},where the inputs are Xt and the outputs are yt.In GPR,instead of obtaining an accurate function between Xt and yt,we focus on exploring their distributions.

Here,we assumed that f denotes an unknown function that maps inputs x to outputs y:f :x→y,and the relation between x and y can be described as

where ∈ is the additional noise added to the f function and it follows the Gaussian distribution

It is assumed that the term f(x) is distributed as a Gaussian process:f(x)～GP(m(x),k (x,x’)),which is defined by the mean m(x) and covariance function k (x,x’).The mean function reflects the expected value at input x:

and empirically,m(x) is often set to zero.

The covariance function k (x,x’) models the dependence between the function values of the input points x and x':

The function k is called the kernel of the Gaussian process.In supervised learning,a similar value of inputs Xt often corresponds with a close response value yt.The covariance function indicates the two latent variables between two inputs x and x',which determines how x is affected by x'.The covariance function k (x,x') can be defined by different kernel functions,e.g.,squared exponential kernel,exponential kernel,and rational quadratic kernel [26]-[29].In this study,we chose one of the most commonly used kernels,i.e.,the squared exponential kernel,as the covariance function.The squared exponential kernel is defined as

where σl is the characteristic length scale,and σ f is the signal standard.

Theoretically,GPR can predict vectors with infinite size; however,in practice,we only make predictions for definite size.Thus,a multivariate normal distribution with a covariance matrix is used to generate the kernel.Let X* be a n ×n matrix with each row a new input xi*,i=1,2,…,n. X*is defined as follows:

We assume m (x )=0 as it simplifies the computation.Subsequently,we can calculate inputs X* from the multivariate normal distribution f* ～N(0,K(X*,X*)),where f*

As we have the observations we intended to predict new inputs X* by drawing f* from the posterior distributionPrevious observations yt and function f* follow a joint normal distribution:

where K (X t,Xt) is the covariance matrix between all observed points,K (X *,X*) is the covariance matrix between newly observed points,K (X t,X*) is the covariance matrix between observed points and newly introduced points,K (X t,X*) is the covariance matrix between new inputs and previous points,I is the identity matrix,and σ∈2I is the additional noise level of observation (i.e.,the variance of ∈).

Using the standard results,the conditional distributionis then a multivariate normal distribution with meanand a covariance matrix

This posterior is also a GP with mean function

and kernel

To predict f*,we can use the mean function in (24) and kernel in (25).

The kernel often contains hyper-parameters,i.e.,length scale,signal variance,and noise variance,which are unknown and need to be inferred from previous observations [30],[31].Commonly,we calculate the hyperparameters by maximizing the marginal (log) likelihood.Given the observation and hyper-parametersthe marginal likelihood iswhere is the covariance matrix of the outputs yt.

In (26),measures the data fit,i.e.,how well the current kernel parameterization explains the dependent variable,is a complexity penalization term,andlog 2π is a normalization constant [32]-[34].

The marginal likelihood is normally maximized through a gradient-ascent-based optimization tool,i.e.,

where α =Ky-1.

2.3 Flowchart of the algorithm

We herein present the flowchart of the algorithm to provide a clear procedure of the GPR in Fig.1.

Fig.1 Flowchart of the Gaussian Process Regression

3 Numerical Examples

Five typical power grids are designed herein to test and validate the proposed method for estimating power grid inertia and predicting its dynamics.These systems were simulated using Simulink software,and the data obtained were processed in MATLAB.The details of the systems are given in Appendix A.

Fig.2 Test systems for inertia estimation and GPR method

3.1 Results of inertia estimation

In this test,two cases are carried out to verify the feasibility of the disturbance-based inertia estimation method.

(1) Case 1 analysis.

In Case 1,the system shown in Fig.2(a) is used.Ten scenarios were created for evaluating the method to estimate inertia.In those scenarios,the generator’s inertia is 4.9,14.9,24.9,34.9,44.9,54.9,64.9,74.9,84.9,and 94.9 kg·m2.The system’s responses are simulated under the change of power loads.Here,the active power of Load 2 is increased by 10 kW at 6.0 s.The system dynamic responses are shown in Figs.3,4,and 5.Fig.3 shows the active power output from the generator,Fig.4 shows the system frequency changes,and Fig.5 shows the RMS value of voltage changes at bus 12.These data points are used as input data for inertia estimation.The estimations of the overall system inertia are summarized in Table1.

Fig.3 Generator’s active power outputs in Case 1

Fig.4 System frequency changes in Case 1

Fig.5 Bus 12 voltage changes in Case 1

From Figs.3,4,5,and Table1,the following observations can be made:

• As the generator’s inertia increases,the overall system’s response rate decreases,which is consistent with the engineering experience.

• The disturbance-based method can successfully calculate the system inertia based on its dynamic trajectories under disturbances.

• The estimated overall system inertia increases as the generator’s inertia increases.

Table1 Estimations of the overall system inertia in Case 1

Scenarios 1 2 3 4 5 Estimated overall system inertia/kg·m2 1.46 28.29 40.99 58.56 71.37 Scenarios 6 7 8 9 10 Estimated system inertia/kg·m2 89.00 110.82 114.02 148.83 159.37

(2) Case 2 analysis.

In Case 2,the system shown in Fig.2 (b) and (c) are used to demonstrate the impact of topology change on the system inertia.The same load change is introduced here as in Case 1.The system dynamic responses are shown in Figs.6,7,and 8.The estimations of the overall system inertia are summarized in Table2.

Fig.6 Generator’s active power outputs in Case 2

Fig.7 System frequency changes in Case 2

Fig.8 Bus 12 voltage changes in Case 2

From Figs.6,7,8,and Table2,the following observations can be made:

• The system topology will impact the overall system’s inertia.

• As the power system’s total impedance increases(Topology 2),its inertia will also increase.

• When designing a system,the overall system’s inertia should be taken into consideration carefully to improve the system’s dynamic performance.

Table2 Estimations of the overall system inertia in Case 2

Topology 1 Topology 2 Topology 3 Generator’s inertia/kg·m2 54.9 54.9 54.9 Estimated system inertia/kg·m2 89.00 552.95 9.29 Generator’s inertia/kg·m2 94.9 94.9 94.9 Estimated system inertia/kg·m2 159.37 827.56 22.02

(3) Case 3 analysis.

In Case 3,the system shown in Fig.2(d) and (e) is used to demonstrate the impact of the number of generators on the system inertia.The same load change is introduced here as in Case 1.The system dynamic responses are shown in Figs.9,10,and 11.The estimation of the overall system inertia is summarized in Table3.

Fig.9 Generator’s active power outputs in Case 3

Fig.10 System frequency changes in Case 3

Fig.11 Bus 12 voltage changes in Case 3

Table3 Estimations of the overall system inertia in Case 3

Topology 4 Topology 5 Generator’s inertia/kg·m2 54.9 54.9 Number of Generator 1 2 Estimated system inertia/kg·m2 106.61 875.49

Figs.9,10,11,and Table3 show the following observations:

• Number of generators in system will impact the overall system’s inertia.

• As the power system’s total number of generators increases,its inertia will also increase.

3.2 Results of dynamics prediction by GPR

In this test,the GPR method is used to predict the dynamics of the system in Fig.2 (a) subject to disturbances.To apply the GPR method,eight variables were sampled—time,the changes of four power loads,the power outputs of the generator,wind speed,temperature,and irradiance level.

Three cases were tested to verify the feasibility of the GPR method.In Case 4,all eight variables are considered and used to train the kernel-based probabilistic model.Fig.12 shows the comparison between the system’s actual dynamics and prediction results in Case 1,with the loss function equal to 52.54,where the loss function is the l2 norm.In Case 5,all power loads are assumed to have the same change pattern.Fig.13 shows the comparison results,with a loss function of 226.38.In Case 6,random values were used to construct the sample data of the eight variables.Fig.14 shows the comparison results,with the loss function equal to 10,816,590,indicating a prediction failure.

From the above predication and comparison results,we make the following observations:

• The dynamics of power loads are important for predicting system trajectories,which is consistent with engineering experience.

• GPR can predict system dynamics,which verifies the feasibility of GPR in system prediction.

Fig.12 Generator’s active power output prediction via GPR in Case 4

Fig.13 Generator’s active power output prediction via GPR in Case 5

Fig.14 Generator’s active power output prediction via GPR in Case 6

• When very limited or incorrect information or wrong one is learned from the training data,e.g.,random data in Case 5,the GPR method will fail,as shown in Fig.14.Thus,learning and using the most relevant data to predict system dynamics is necessary.

4 Conclusion

A disturbance-based method is presented herein to estimate the power grid inertia,thereby determining the resilient operations of systems.The GPR method is then used to predict the trajectory of the system after disturbance.Simulation results on typical power grids are used to test and validate the proposed data-driven method for estimating the inertia of the system as well a s predicting the dynamics operations of power grids subject to disturbances.

Appendix A

The line impedances in topology 1 of Fig.2 are given in TableA1.The impedances from bus 13 to bus 24,from bus 25 to bus 36,and from bus 37 to bus 48 are the same to those from bus 1 to bus 12.

TableA1 Line impedances between nodes in topology 1 of Fig.1

From To Resistance/(Ω/km)Reactance/(H/km)Length/m 1 2 0.284 0.2202e-3 35 2 3 0.284 0.2202e-3 35 3 4 0.284 0.2202e-3 35 4 5 0.284 0.2202e-3 35 5 6 0.284 0.2202e-3 35 6 7 0.284 0.2202e-3 35 7 8 0.284 0.2202e-3 35 8 9 0.284 0.2202e-3 35 9 10 1.380 0.2175e-3 30 10 11 1.380 0.2175e-3 30 11 12 1.380 0.2175e-3 30

The parameters of the generators of each topology in Fig.2 are the same,and its values are presented in TableA2.

TableA2 Parameters of generator

Parameter Value Nominal Power (Pn/VA) 111.9e3 Voltage 440*sqrt(3)Frequency 60 Stator (Rs/ohm,Ll,Lmd,Lmq/H,Lc/H) 0.26,1.14e-3,13.7e-3,11.0e-3 Field (Rf'/ohm),Llfd'/H) 0.13,2.1e-3 Dampers(Rkd',Llkd',Rkq1',Llkq1') 0.0224,1.4e-3,0.02,1e-3 Friction factor 0 Pole pairs 2 Initial conditions (dw/%,th/deg,ia,ib,ic/A,pha,phb,phc/de,Vf/V)0,-111.48,53.98,53.98,53.98,173.30,-293.30,-53.30,17.89 Active power generation (P/W) 50,000

The load characteristics (Load1,Load2,Load3,Load4,and Load5) of each topology in Fig.2 are the same,and its values are presented in TableA3.

TableA3 Characteristics of load

Parameter Value Load flow Constant Z Nominal frequency/(fn/Hz) 60 Nominal phase-to-phase voltage/(Vn/Vrms) 762 Active power/(P/W) 100,000 Configuration Y(ground)

Declaration of Competing Interest

We declare that we have no conflict of interest.