Calibrating synchronous-generator-interfaced DG models in microgrids using multiple event data

Calibrating synchronous-generator-interfaced DG models in microgrids using multiple event data

Electrical Power and Energy Systems 120 (2020) 105989 Contents lists available at ScienceDirect Electrical Power and Energy Systems journal homepage...

4MB Sizes 0 Downloads 55 Views

Electrical Power and Energy Systems 120 (2020) 105989

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Calibrating synchronous-generator-interfaced DG models in microgrids using multiple event data☆

T

Zhiwen Wanga, Yijun Chena, Chen-Ching Liub,c, Yin Xud, , Kefei Moe, Kevin P. Schneiderf, Francis K. Tuffnerf, Dan T. Tong ⁎

a

Department of Power Systems, China Electric Power Research Institute (CEPRI), Beijing 100192, China Electrical and Computer Engineering at Virginia Tech (VT), Blacksburg, VA 24061, USA c Electrical Engineering and Computer Science, Washington State University (WSU), Pullman, WA 99164, USA d School of Electrical Engineering, Beijing Jiaotong University (BJTU), Beijing 100044, China e School of Electrical Engineering and Computer Science, Washington State University (WSU), Pullman, WA 99164, USA f Pacific Northwest National Laboratory (PNNL) Located at the Battelle Seattle Research Center in Seattle, WA, USA g U.S. Department of Energy (DOE) Office of Electricity Delivery and Energy Reliability (OE), Washington, DC 20585, USA b

ARTICLE INFO

ABSTRACT

Keywords: Microgrid Distributed generator Parameter estimation Model validation Least square optimization Unscented Kalman filter

Microgrids, consisting of distributed generators (DGs), loads, and energy storage systems, play an important role in future smart grids. In order to evaluate microgrid operations and their dynamic response, it is necessary to develop accurate dynamic models of DGs in microgrids. This paper presents a systematic method to calibrate the model of a synchronous-generator-interfaced DG. Underlying parameters are categorized into two groups, namely, steady-state parameters and time constants, which are estimated in two successive stages using multiple event data. This methodology ensures that the validated model with the estimated parameters are applicable for a set of events. Moreover, assuming that power angle measurements of a synchronous-generator-interfaced DG are unavailable, a modified unscented Kalman filter (UKF) is proposed to facilitate parameter estimation. The effectiveness of the proposed method is validated using field test data of a 2600 kVA diesel generator which is part of an operational microgrid.

1. Introduction FUTURE smart grids will integrate the control and operations of multiple critical infrastructures to increase reliability and resiliency [1]. Because microgrids can be used to restore critical load when utility service is interrupted after major events, they are key components of a resilient power system [2]. In microgrids, a synchronous-generator-interfaced DG, such as a diesel or gas generator, which has black-start capability, is an important resilience source [3,4]. Since the generation capacity of DGs in microgrid is usually small, microgrids do not have sufficient capability to withstand large transients, which can cause instability. Hence, dynamic performance of synchronous-generator-interfaced DGs in load restoration must be carefully evaluated. This requires developing accurate models of DGs. Inaccurate models may lead to unconvincing evaluation results and improper guidance for planning

and operation. Thus, it is important to calibrate the DG models. As the deployment of phasor measurement units (PMUs) in power systems increases, there is an increasing interest in estimating synchronous generator parameters using synchrophasor data [5,6]. The methods in [7–9] adopt a hybrid dynamic simulation (HDS) concept to perform model validation. Simulations with models to be calibrated are performed iteratively. In each step, errors between simulated outputs and PMU data are computed and underlying parameters are updated until the errors reach a local minimum. However, iterative dynamic simulations in these methods are time-consuming. In addition to HDS methods, another popular method is the use of Kalman filters. In [10], the extended Kalman filter (EKF) is applied to calibrate generator models. In [11], an expectation-maximization algorithm is used to enhance noise immunity of the EKF. Due to the use of a linear approximation of generator models by EKF, the level of accuracy and

☆ This work was supported by Office of Electricity Delivery and Energy Reliability (OE), U.S. Department of Energy (DOE) and Pacific Northwest National Laboratory (PNNL). This work was also supported by National Natural Science Foundation of China (grant ID: 5177071141). ⁎ Corresponding author. E-mail addresses: [email protected] (Y. Chen), [email protected] (C.-C. Liu), [email protected] (Y. Xu), [email protected] (K. Mo), [email protected] (K.P. Schneider), [email protected] (F.K. Tuffner), [email protected] (D.T. Ton).

https://doi.org/10.1016/j.ijepes.2020.105989 Received 12 September 2019; Received in revised form 22 January 2020; Accepted 5 March 2020 Available online 20 March 2020 0142-0615/ © 2020 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

robustness is reduced. To improve accuracy and robustness, nonlinear parameter estimators, such as Ensemble Kalman filter (EnKF) [12,13] and the Unscented Kalman filter (UKF) [14–16], are proposed. EnKF and UKF retain full nonlinearity of error statistics through propagation and, therefore, they achieve better estimation results relative to the linearized EKF. Although significant progress of parameter estimation methods has been made [5–16], some issues remain. First, the majority of methods use data from a single event to estimate parameters. As a result, the estimated parameters may only be accurate for a single event. Additionally, when multiple parameters are estimated simultaneously using data from a single event, the Kalman-filter-based methods may not converge [17]. Second, the majority of the reviewed methods [11–16] are implemented at the transmission level, where there are sufficient measurements, including power angle. In contrast, at the distribution level, measurements are typically incomplete. In a distribution system, typical measurements include active/reactive power and voltage at the substation. However, some generator values, such as the power angle, are not measured. For small generators in a microgrid, it is possible to conduct experiments and record performance of the generators under different operational events, such as load pick-up and voltage regulation. Therefore, data from multiple events can be obtained. To this end, this paper develops a systematic method to address those issues in calibrating a synchronous-generator-interfaced DG model in a microgrid. Contributions of this paper are three-fold:

reactance and transient reactance of d and q axes, respectively. Ef is excitation voltage. δ is power angle. Ut is terminal voltage magnitude. Uref is terminal voltage reference. For a diesel generator, a simple voltage exciter can be represented by a first-order differential equation:

TA Ef = (Uref

Ut ) KA

(3)

Ef

where KA is gain, TA is time constant, and Uref is voltage reference. For a synchronous-generator-interfaced DG, there are two measurement equations, also called output equations:

Pe = Ut sin

Qe = Ut cos

E

q

Ut cos X

Eq

d

Ut cos X

d

+ Ut cos

Ut sin X

Ut sin

Ut sin X

Ed (4)

q

E q

d

(5)

Remark 1. Rotor dynamic equations are not involved and, therefore, the power angle δ is not a state variable. It is an input. If governor models are considered, they can be treated in a similar way as an exciter. In such a case, the proposed method is still applicable. The model structure of a synchronous-generator-interfaced DG is determined by Eqs. (1)–(5). To calibrate a synchronous-generator-interfaced machine model, it is necessary to identify underlying para' ' meters, i.e., Tdo , Tqo , TA , Xd , Xq , Xd' , Xq' and KA . The goal of parameter estimation is to find optimal parameter values, with which the model structure can produce outputs that minimize the error when they are compared to measured data. In other words, with the optimal parameter values, the error between reproduced outputs of the model and real measurements is minimized. There are three challenges for estimating these parameters:

(1) Parameter estimation of a synchronous-generator-interfaced DG is formulated as a two-stage problem. Parameters are categorized into two groups: steady-state parameters and time constants. The two types of parameters are estimated in different stages, successively. (2) Multiple event data is used simultaneously to estimate underlying parameters. By doing so, the estimated parameters are applicable for a set of events. A modified UKF is proposed to deal with the case in which measurement of power angle of a synchronous-generatorinterfaced DG is unavailable. (3) In order to verify the effectiveness of the proposed method, a field test using a resistive load bank is designed and performed with the Washington State University (WSU) Microgrid, in which important measurements, including steady-state and dynamic data, are captured. Based on the collected data, models of the operational diesel generators are calibrated and validated.

(1) Underlying parameters can be regarded as extended states. Usually, one may use different kinds of Kalman filters to estimate those extended states. However, when the number of generator parameters increases, Kalman filter may not converge [17]. An example is provided in Fig. 1, where UKF is used to estimate three para' meters, Xd , Xq and Tdo , simultaneously. As shown in Fig. 1, it fails to converge. (2) Typical Kalman filter use measurement data from a single event to estimate parameters. Different event data may result in different estimated parameters, leading to inconsistent results. (3) In microgrids, measurements of some values, i.e., power angles in this paper, of synchronous machines are not available. As a result, many ready-made parameter estimation methods [11–16] do not apply.

The remainder of this paper is organized as follows. Section 2 states the problem of parameter estimation for a synchronous-generator-interfaced machine in a microgrid. The framework of the proposed method is provided in Section 3. Estimation of steady-state parameters and time constants is detailed in Sections 4 and 5, respectively. In Section 6, the effectiveness of the proposed method is validated using the field test data. Conclusions are given in Section 7.

3. Framework In this section, underlying parameters are categorized into two groups, namely, steady-state parameters and time constants. Then, a

2. Problem statement Model validation consists of two phases: determining the model structure and identifying specific parameters. For calibration of a synchronous-generator-interfaced DG, such as a synchronous diesel generator, the model structure can be described by differential equations [18]. In this paper, a synchronous-generator-interfaced DG model is represented by differential equations:

Td0 Eq = Ef

Eq

(X d

X d)

Tq0 Ed =

Ed + (Xq

X q)

Eq

Ut cos X

d

Ut sin X

q

(1)

E

d

(2)

where are d-axis transient voltage, q-axis transient voltage. ' ' They are state variables. Tdo , Tqo are time constants. Xd , Xq , Xd' , Xq' are

Eq' ,

Ed'

Fig. 1. UKF fails to converge when it is used to estimate three parameters. 2

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

two-stage framework for parameter estimation of a synchronous-generator-interfaced DG is established. Details are provided as follows. A. Parameter classification Eqs. (1)–(5) can be formulated in compact forms: (6)

Ti x i = fi (x , p, u, u *)

(7)

ym = hm (x , p, u, u *) where x represents state variables, i.e.,

Eq' ,

and Ef . Ti denotes time ' ' constant, i.e., Tqo , Tdo and TA . p denotes parameters as follows: Xd , Xq , Xd' , Xq' and KA . u represents known inputs, i.e., Ut and Uref . u* represents an unknown input, i.e., δ. ym denotes measurements, i.e., Pe and Qe . Eq. (6) represents Eqs. (1)–(3). Eq. (7) represents Eqs. (4) and (5). When a synchronous-generator-interfaced DG is in a steady state, i.e., x i = 0 , the differential Eq. (6) degenerates to an algebraic equation where Ti is not involved:

Ed'

(8)

0 = fi (x , p, u, u *)

Fig. 3. Procedure of parameter estimation.

That is, in a steady state, only p remains, while Ti disappears. In this sense, p is classified as “steady-state parameter.” Note that Ti only appears in the left-hand term of differential Eq. (6). Therefore, Ti is classified as “time constant.”

residuals of a least-squared optimization. Details of this process are provided in Section 4. Stage 2 In this stage, differential Eq. (6) and output Eq. (7) are needed. Since steady-state Eq. (8) does not apply for dynamic data, it is not needed in stage 2. Suppose that there are multiple dynamic events, such as several load pick-up events. First, in each event, dynamics of a synchronous-generator-interfaced DG will be represented by Eqs. (6,7). Then, all sets of equations are merged to form an aggregated dynamic system. In the aggregated system, time constants Ti are regarded as extended state variables, which are estimated using a modified UKF. Details of the second stage are provided in Section 5.

B. Measurement data classification Measurement data are categorized into two groups: steady-state data and dynamic data. Steady-state data are measurements of active power, reactive power and terminal voltage magnitude when a synchronous-generatorinterfaced DG is in a steady state, i.e., x i = 0 . In practice, when there is no more event (such as load pick-up or a short-circuit fault) for several minutes, measurements can be regarded as steady-state data. On the 0 , i.e., when a synchronous-generator-interfaced other hand, when x i DG is in a dynamic process, measurements are called dynamic data. Steady-state data and dynamic data are illustrated in Fig. 2. C. Two stages of parameter estimation Steady-state data are used to estimate steady-state parameters p in Stage 1. Thereafter, p will be sent into Stage 2, where dynamic data are used to determine time constant Ti . The procedure is shown in Fig. 3. The parameter estimation consists of two stages:

Remark 2. Steady-state parameters p are estimated in stage 1 and, therefore, only a small number of time constants Ti will be estimated using Kalman filter in stage 2. Therefore, the first challenge in Section 2 is overcome. 4. Estimation of steady-state parameters Suppose that there are S steady states. For the sth steady state, there are I equations of fis (x s , p , us, u s ) = 0 , and M equations of hms (x s , p , us, u s ) yms = 0 . Then, there are S × (I + M ) equations in total. Denote the number of x s , p , u s by I, B, and C, respectively. Then, the number of unknown variables (x s, p, u s ) is S × (I + C ) + B . When the number of equations is greater than the number of unknown variables, one can obtain a solution of unknown variables in the least-squared sense. Based on this, if S × (I + M ) is greater than S × (I + C ) + B , least-squared estimation of unknown variables (x s, p, u s ) can be obtained. Specifically, for a synchronous-generator-interfaced DG model (1)(5), I = 3, M = 2, B = 5, C = 1. Then, when S 5, the number of equations will be greater than that of unknown variables, one can estimate p in the least-squared sense. A full least-squared optimization model that determines steady-state parameters p for a synchronous-generator-interfaced DG is given as follows:

Stage 1 In this stage, the output Eq. (7) and steady-state Eq. (8) are needed. Since Eq. (6) degenerates to Eq. (8), it is not needed in stage 1. If the number of steady-state events is large enough, the number of Eqs. (7) and (8) can be higher than the number of parameters. Estimated parameters p are those that minimize the

Pe

dynamic data

steady-state data

t

I

Qe

dynamic data

steady-state data

S

min

~s x s, p, u

s=1

i=1 M

~ s) + fis (x s , p, us, u yms

~ s) hms (x s , p, us , u

(9)

m=1

where

t

f1s = Efs

Fig. 2. Dynamic data and steady-state data in an event. 3

Eq s

(Xd

X d)

E

s

q

Uts cos X

d

s

(10)

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

f2s =

Ed s + (Xq

s f3s = (Uref

y1s

Uts ) KA

X q)

Uts sin

s

X

E

d

s

(11)

q

Efs

(13)

y2s = Qes

(14)

=

h1s = Uts sin

sE q

s

h 2s = Uts cos

s

sE q

(15)

Uts cos s X d

(16)

(22) (23)

d

w1 + wd wD

f (x D , p , u D , u * D )

y1

h (x d, ud , u *d)

v1

yd

= h (x d, ud , u *d)

+ vd

yD

h (xD, uD, u

D)

vD

(24)

yk = h (x k , uk , u *k ) + vk

(25)

A common algorithm for a nonlinear state estimation problem is the UKF, which consists of three steps: sigma point configuration, prediction, and correction. A full description of the UKF is provided in [17]. According to the UKF algorithm, inputs uk are needed in the prediction and the correction steps. When there are unknown inputs u , prior and posterior estimations of states cannot be produced, because the f and h function cannot be computed. In this study, the power angle is unavailable, therefore, the UKF method does not apply. Although there are several methods [21–23] that extend traditional Kalman filter to deal with unknown inputs, they do not apply in this paper. The methods in [21,22] assume that f function is linear (or approximately linear). However, differential Eqs. (1)–(3) of a synchronous-generator-interfaced DG have strong nonlinearity (due to sinusoidal functions). The method in [23] requires that unknown input u does not exist in h function, and f function should have an affine form, i.e., u can be separated from the rest of f function. Since the unknown input u in this paper is power angle δ, functions of f and h (1)–(5) do not meet such requirements. To handle the unknown inputs, a modified UKF is proposed as follows. Note that power angle δ can be expressed as Eq. (26):

where w denotes process noises, and v denotes measurement noises. Let the D dynamic systems represented by (17) and (18) be aggregated into one dynamic system:

f (x d , p, ud, u *d)

x k = f (xk - 1, uk , u *k ) + wk

B. Unscented Kalman Filter with unknown inputs

d

f (x 1, p, u1, u *1)

D

Remark 3. Several methods use data of a single event as inputs for an (I + I)-order Kalman filter to estimate states/parameter, where I is the order of a synchronous machine model, as well as the number of time constants. In theory, the estimated states/parameters are optimal solutions for observed data in the single event. Thereafter, the estimated states/parameters can be validated empirically using different event data. As a comparison, the proposed method uses D event data as inputs for an extended (DI + I)-order Kalman filter. Time constants T are commonly shared in all dynamic events. The estimated parameters T are optimal solutions for the D different event data, which means T values are applicable for a wider range of events. In this regard, the second challenge in Section 2 is overcome.

(18)

y d = h (x d , p, ud, u *d) + v d

d

At this point, if the dynamic states xk of the discrete system can be estimated from measurements, parameters T are identified. Then, the parameter identification is equivalent to a state estimation problem. Note that f and h are nonlinear. Hence, this is a nonlinear state estimation problem.

(17)

Tx d = f (x d , p, ud, u *d) + w d

(21)

where x¯k stands for x ,…, x ,…, x , and T in the k step. For brevity, symbols xk and xk in the subsequent section will be used indiscriminately. Thereafter, Eq. (22) is denoted as:

Once steady-state parameters p are obtained, the next step is to estimate time constants T = [T1…Ti…TI]. Suppose that there are D dynamic events. For the dth dynamic event, let x d , ud , u d denote state variables, known inputs, and unknown inputs, respectively. Thereafter, the model of a synchronous-generator-interfaced DG in the dth dynamic event is represented by Eqs. (17) and (18).

T xD

wd wD 0

f ( x D , p , uD , u * D ) 0

1

A. Extended state-space equation

=

+

yk = h (¯x k , uk , u *k ) + vk

5. Time constant estimation

T xd

f (x d, p , ud, u *d)

x¯ k = f¯ (¯x k - 1, uk , u *k ) + w ¯k

This model is aimed at finding optimal decision variables of xs, p, u s such that the objective function (the square of residuals) is minimized. This is a typical least-squared optimization problem, which can be solved by commercial solvers, e.g., linsqr in MATLAB [19]. It should be noted that due to non-linearity of operating constraints, the true values of generator parameters can be obtained if initial values are sufficiently close to the true values. Generally, there is no theoretical guarantee for identification of all true parameters. This phenomenon is known as “plateau phenomenon” [20]. Although several empirical measures can help ease this problem, further investigations are needed for the future study.

T x1

w1

Given a sampling frequency, the continuous state-space Eq. (21) and output Eq. (20) are converted to discrete forms:

Uts cos s

s s E s d s Ut sin X q

Uts sin

=

TxD T

X d

s s E s d s Ut sin X q

+ Uts cos

T xd

(12)

Pes

f (x1, p , u1, u *1)

T x1

(19)

(20)

u*k

Note that for all D events, time constants are the same, and stay unchanged during the dynamic process. The time constants T can be regarded as extended state variables. Then, an aggregated dynamic system with T can be represented as follows:

k

=

k 1

+(

e )/Fsample

(26)

where Fsample is sampling frequency, e.g., Fsample = 1 kHz. ω is speed of rotor, ωe is electrical frequency of terminal voltage phasor. If the rotor speed deviation (ω − ωe)/ Fsample is small, it can be regarded as pseudo-process noise w’k. Then, Eq. (26) can be rewritten as 4

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

where Rk is covariance of measurement noise vk. Finally, the Kalman Gain Kk is obtained, and the posterior states and the covariance matrix Pk are computed:

follows: k

=

+ wk

k 1

u*k = u *k

1

(27)

+ wk

K k = Ck Sk 1 (46) ~ Pk = Pk K k Sk KkT (47) xk x~k = ~ + Kk [yk uk uk

Eq. (27) shows that that one can use uk 1 as a prior estimation of uk in the prediction step. By doing so, the f function can be computed, and a prior estimation of states xk can be produced. In the correction step, a posterior estimation of uk will be updated. Therefore, the UKF can work with unknown inputs uk . The computation procedure of Unscented Kalman Filter with unknown inputs (UKF_UI) is detailed as follows. First, sigma point of both xk and uk are configured. Then, in the prediction step, the prior estimations ofuk and xk are computed. In the correction step, the posterior estimations of xk and uk are updated based on measurements.

Remark 3. The major differences between the traditional UKF and the proposed UKF_UI are the prior estimation uk and posterior estimation uk . Since the proposed UKF can estimate the unknown inputs, it is able to handle the case where power angle measurement is unavailable. Consequently, the third challenge in Section 2 is overcome. Since the derivation of speed is modeled as pseudo-process noise, the effect of this assumption is reflected in diagonal elements of correlation matrix Q’k-1. Usually, as diagonal elements increase, e.g. the deviation of speed becomes larger and the variance of estimated time constants would increase. It should be noted that if the deviation of speed is large, it cannot be included in the process noise. One can build generator models with different underlying parameters in the simulation environment and check if the assumption is valid during dynamic events.

Step 1 Sigma Point Configuration Suppose that xk-1 is of I dimension, uk 1 is of L dimension. The “sigma points” and a covariance matrix Pk-1 of n dimension are computed as follows: xk 1 x~k0 1 = (28) ~ 0 uk 1 u k 1

x~ki ~ u ki

1

xk uk

=

1

1 1

xk x~kn +1i ~ n + i = uk u k 1

+ ( (n + ) Pk

1

1 )i ,

( (n + ) Pk

1

1 )i ,

i = 1,

, n (29)

i = 1,

, n (30)

= 3 2 n (31) where n = I + L, λ and α are pre-defined parameters of UKF. Step 2 UKF Prediction In the second step, each sigma point obtained in step 1 is used to perform a prediction of state variables xk and unknown inputs uk . There are 2n + 1 predictions, as shown below: ~ i ,u ) x~ki f (~ xki 1, u k 1 k 1 , i = 0, 1, , 2n (32) ~ i = u ~ i u k k 1

6. Experimental results A. Test system and data The Washington State University microgrid includes a 2600 kVA diesel generator. The microgrid is aimed at providing backup power for the campus in case of a failure of utility service. The rated voltage of the diesel generator is 4.16 kV. The voltage reference of the diesel generator can be set at 0.95p.u., 1.00p.u., and 1.05p.u. In order to collect steady-state and dynamic data of multiple events, a rented resistive load bank with four levels (400 kW, 800 kW, 1200 kW, 1600 kW at rated voltage) is installed. The diesel generator can pick up different levels of resistive loads with different voltage reference settings. Measurements of three-phase waveforms of terminal voltage Uabc and output current Iabc are recorded with a sampling frequency Fsample = 1 kHz. The sampling rate plays an important role in Kalman filter. In this paper, sub-transients, e.g. E'd and E'q , of a synchronous machine are considered. Typical time scales of dynamics of E'd and E'q fall between 0.5 s and 10 s [18]. According to Nyquist sampling theorem, the sampling rate should be two times higher than the highest frequency component. Therefore, recorders of 1 kHz sampling rate are used in the field test, which are sufficient to capture dynamics of sub-transients of synchronous machines [24]. Based on waveforms of Uabc and Iabc, the magnitudes of terminal voltage Ut, instantaneous active power Pe, and reactive Qe are computed using the method in Appendix A.

Next, prior estimations of xk and uk are computed as a weighted summation of the 2n + 1 predictions: 2n x~ki x~k Wci ~ i (33) ~ = u uk k i=0 (34)

Wc0 = Wci

n+ 1 , 2(n + )

=

i = 1,

, 2n (35)

A prior covariance matrix is obtained in (36): 2n

~ Pk =

i=0

+ Wc0

=

Wci

=

x~ki ~ uki

Wci

Qk

~ xki ~ u ki

x~k ~ u k

x~k ~ u k

T

(36)

1

Qk

1 2)

+ (1 +

n+ 1 , 2(n + )

i = 1,

(37)

, 2n (38)

where Qk-1 is the covariance of process noises wk, Q’k-1 is covariance of pseudoprocess noises w’k. β is a pre-defined parameter of UKF. Step 3 UKF Correction Step In the third step, the 2n + 1 sigma points should be updated as follows: x k0 x~k = ~ (39) u u 0 k

k

~ xk = ~ u

x ki u ki

k

x kn + i u kn + i

=

~ + ( (n + ) Pk )i , i = 1,

x~k ~ u k

~ ( (n + ) Pk )i , i = 1,

B. Test results of steady-state parameters

, n (40)

There are 4 load levels of the resistive load bank and 3 levels of terminal voltage references. By setting load levels and regulating voltage references, there are 12 steady states in total. Measurements of Uref, Ut, Pe, and Qe in the 12 steady-state events for 1-second are shown in Fig. 4. Steady-state data of the 12 events are divided into a training set (including 7 events) and a test set (including 5 events). In the training set, the number of steady-state equations fis (x s , p , us, u s ) = 0 and hms (x s , p , us, u s ) yms = 0 is 35. The number of unknown variables ( x s , p , u s ) is 33. Mean values of Uref, Ut, Pe, and Qe are fed into the least-squared optimization problem (9), for estimation of Xd , Xq , Xd' , Xq' and KA . Thereafter, the estimated parameters and manufacturer’s parameters are compared in the test set. Test results are shown in

, n (41)

With the new sigma points, one can compute 2n + 1 estimated outputs:

yki = h (x ki , uki, uk) ,

i = 0, 1,

, 2n (42)

Let µk denote a weighted sum of the 2n + 1 estimated outputs

µk =

2n i=0

i k

Wci y ki (43)

At this point, one can compute the measurements covariance matrix Sk and the cross-covariance matrix Ck:

Sk =

Ck =

2n

i=0

2n i=0

Wci [(y ki

Wci

µk )T] + Rk (44)

µk )(y ki

x~ki ~ i u k

x~k ~ u k

(yki

µk ](48)

where yk are measurements in the k step.

µk )T (45)

5

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al. 1.10

1.10

1.05

1.05

0.06

0.95

1.00

0.4

0.04

Qe (p.u.)

1.00

0.5

Pe (p.u.)

Ut (p.u.)

Uref (p.u.)

0.6

0.3 0.2

0.95

0.02

0.1 0.90

0

0.20

0.40 0.60 Time(s)

0.80

1.00

0.90 0

0.20

0.40 0.60 Time(s)

0.80

0

1.00

0

0.20

0.40 0.60 Time(s)

0.80

0 0

1.00

0.20

0.40 0.60 Time(s)

0.80

1.00

Fig. 4. Recorded data of Uref, Ut, Pe, Qe of the diesel generator in 12 steady-state events.

5.5 5.0

Va lue (s)

4.5

4.0 3.5 0

Fig. 5. Estimated values and manufacturer’s values of steady-state parameters.

0

5.00

-0.005

C. Test results of time constants

-0.01 -0.015

There are 14 dynamic events with different load levels and terminal voltage references, details of which are provided in Appendix B. Measurements of Uref, Ut, Pe, and Qe for 5 s are shown in Fig. 7. The first 7 events are used as the training set for estimation of time constants, and the last 7 events are used as the test set for validation. ' Based on the training set, the estimated value of Tdo along with time is shown in Fig. 8, which starts from 4.38 s and converges to 4.65 s. The ' curves of Tqo and TA are similar to Fig. 8. They are omitted due to

Manufacture parameters Estimated parameters 0

5

10 15 No. of equations

20

25

Fig. 6. Residual comparison of steady-state equations in the test set.

1.10

1.05

1.10

0.7

1.05

0.6

1.00

0.95

1.00

Pe (p.u.)

Ut (p.u.)

0.10

0.95 0.90

1.00

2.00 3.00 Time(s)

4.00

5.00

0.85 0

0.5

Qe (p.u.)

Resi du als

4.00

Figs. 5 and 6. According to Fig. 5, the estimated values deviate from the manufacturer’s values of steady-state parameters. From Fig. 6, it can be seen that, in the test set, residuals of the steady-state equations fis and hms with the estimated values are much smaller than those of the manufacturer’s values. This indicates that the estimated parameters can better fit the measurement data to achieve a higher accuracy.

0.005

Uref (p.u.)

2.00 3.00 Time(s)

' Fig. 8. Estimated Tdo in each time step.

0.01

0.90 0

1.00

0.4

0.05

0.3

1.00

2.00 3.00 Time(s)

4.00

5.00

0.2

0

1.00

2.00 3.00 Time(s)

4.00

5.00

Fig. 7. Recorded data of Uref, Ut, Pe, Qe of the diesel generator in 14 dynamic events. 6

0

0

1.00

2.00 3.00 Time(s)

4.00

5.00

Electrical Power and Energy Systems 120 (2020) 105989

0.7

0.7

0.6

0.6

0.5

0.5

Active Power (p.u.)

Active Power (p.u.)

Z. Wang, et al.

0.4 0.3 0.2

Estimated outputs

0.1

0.3

0

1.00

2.00

Time (s)

3.00

4.00

-0.1 0

5.00

0.7

0.6

0.6 Estimated outputs Measurements

0.3 0.2

0

Time (s)

4.00

4.00

5.00

Measurements

0.2

0

3.00

3.00

0.3

0.1

2.00

Time (s)

Reproduced outputs

0.4

0.1

1.00

2.00

0.5

React ive Power (p.u.)

0.4

1.00

Fig. 11. Measured Pe and reproduced Pe in test set.

0.7

0.5

React ive Power (p.u.)

Measurements

0

Fig. 9. Measured Pe and estimated Pe in training set.

-0.1 0

Reproduced outputs

0.2 0.1

Measurements

-0.1 0

0.4

-0.1 0

5.00

1.00

2.00

Time (s)

3.00

4.00

5.00

Fig. 12. Measured Qe and reproduced Qe in test set.

Fig. 10. Measured Qe and estimated Qe in training set. Table 1 Maximum errors of power when UKF_UI converges in training set.

Table 2 Maximum errors of power when UKF_UI converges in test set.

Event No.

Active power error (p.u.)

Reactive power error (p.u.)

Event No.

Active power error (p.u.)

Reactive power error (p.u.)

1 2 3 4 5 6 7

0.0177 0.0177 0.0189 0.0187 0.0182 0.0191 0.0079

0.0286 0.0220 0.0366 0.0318 0.0241 0.0378 0.0154

8 9 10 11 12 13 14

0.0176 0.0180 0.0207 0.0231 0.0182 0.0178 0.0195

0.0135 0.0188 0.0287 0.0351 0.0158 0.0132 0.0205

limitation of space. Their estimated values are 1.05 s and 0.07 s, respectively. During the estimation procedure, the estimated outputs Pe and Qe, produced by UKF_UI using Eq. (43), are compared with measured values. The results are shown in Figs. 9 and 10. It takes about 20 time steps (about 20 ms) for UKF_UI to converge after it starts. Maximum errors of active/reactive power when Kalman filter converges are available in Table 1. It can be seen that the estimated outputs match well with measured values in both Pe and Qe, indicating that estimated results accurately reflect the actual dynamics of the synchronous-generator-interfaced DG. With the estimated steady-state parameters and time constants, one

can reproduce outputs Pe and Qe for the test set. The reproduced outputs and real measurements are compared in Figs. 11 and 12. Maximum errors of active/reactive power are available in Table 2. The errors are small, indicating that the parameters estimated from the training set are applicable to the test set, which validates the proposed method. Evaluation of the estimated parameters’ uncertainty is as important as estimation itself. The evaluation results are provided in Appendix C. 7. Conclusion A two-stage parameter estimation framework is proposed to 7

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

calibrate a synchronous-generator-interfaced DG model. The effectiveness of the proposed method is validated by a field test with the WSU microgrid. With the estimated parameters, simulation with a synchronous-generator-interfaced DG model reproduces the measured curves. Although the proposed method is validated in a field test using a resistive load bank, it can be extended for generation parameter estimation in the real world. For steady-state parameter estimation, when the number of steady-state constraints is greater than that of unknown variables, one can use the proposed method for estimation. For dynamic-state parameter estimation, when active power, reactive power and terminal voltage are captured using PMUs or recorders in similar events such as picking-up loads, the modified UKF can be used to estimate time constants. For future study, as frequency stability is an important issue in microgrids, the proposed method will be used to estimate inertial time constants and calibrate generator governor models.

CRediT authorship contribution statement Zhiwen Wang: Methodology, Writing - original draft. Yijun Chen: Software, Validation. Chen-Ching Liu: Conceptualization, Supervision. Yin Xu: Methodology, Writing - review & editing. Kefei Mo: Software, Validation. Kevin P. Schneider: Resources. Francis K. Tuffner: Resources. Dan T. Ton: Resources. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Computation of Instantaneous Active and Reactive Power from Three-phase Measurements First, voltage, current and their angle in static α-β axes are computed from three-phase measurements: 2

U = 3 (Ua 1 3

U =

0.5Ub

(Ub

2

I = 3 (Ia I =

1 3

0.5Uc )

Uc )

0.5Ib

(Ib

0.5Ic )

Ic ) U

= arctan U

(A-1)

Second, voltage and current in d-q axes are computed:

Ud = U cos Uq =

Id = I cos Iq =

+ U sin

U sin + U cos + I sin (A-2)

I sin + I cos

Finally, instantaneous active power and reactive power are computed based on d-q components of voltage and current:

Pe = 1.5(Ud Id + Uq Iq ) Qe = 1.5(Uq Id Ut =

2

U +U

Ud Iq ) 2

(A-3)

Appendix B Descriptions of 14 dynamic events (see Table B1).

Table B1 Descriptions of 14 dynamic events. No.

Uref (p.u.)

Event

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1.00 0.95 1.05 1.00 0.95 1.05 1.00 0.95 1.05 1.00 1.05 1.00 0.95 1.05

The The The The The The The The The The The The The The

8

generator generator generator generator generator generator generator generator generator generator generator generator generator generator

picks picks picks picks picks picks picks picks picks picks picks picks picks picks

loads loads loads loads loads loads loads loads loads loads loads loads loads loads

from from from from from from from from from from from from from from

0 to780 kW 0 to700 kW 0 to 840 kW 400 to 1140 kW 363 to 1034 kW 441 to 1244 kW 751 to 1120 kW 683 to 1022 kW 829 to 1243 kW 750 to 1576 kW 810 to 1736 kW 1100 to 1475 kW 1002 to 1339 kW 1206 to 1614 kW

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

Appendix C Evaluation of estimated parameters’ uncertainty Since steady-state parameter estimation is formulated as a non-linear problem in Eq. (9), it is impossible to derive analytical formula for distributions of steady-state parameters from the measurements. Therefore, Monte-Carlo simulation is used to estimate distributions of steady-state parameters and time. Note that there are multiple independent measurements of Uref, Ut, Pe, and Qe in the training set (8 steady-state events) for 1-second with a sampling rate 1 kHz, i.e. there are 1000 sets of samples per event. One can use the mth set of samples as inputs to solve the least-squared optimization problem, and obtain estimation of Xd , Xq , Xd' , Xq' and KA . In this paper, this procedure is called the mth point estimation of steady-state parameters. All M (M = 1000) point estimations constitute distributions of steady-state parameters. When the mth point estimation of steady-state parameters is available, one can use UKF_UI to obtain 2n + 1 sigma points of each time constant, which correspond to the mth point estimation. All M point estimations indicate that there are M(2n + 1) sigma points of for each time constant in total. Those sigma points are used to form posterior distributions. Distributions of estimated parameters are shown in Fig. C1.

Fig. C1. Distributions of estimated parameters.

Trans. Smart Grid. early access article. [4] Hossain MA, Pota HR, Hossain MJ, et al. Evolution of microgrids with converterinterfaced generations: challenges and opportunities. Int J Electr Power Energy Syst 2019;109:160–86. [5] Lu S, Hogg BW. Dynamic nonlinear modelling of power plant by physical principles and neural networks. Int J Electr Power Energy Syst 2000;22(1):67–78. [6] Arjona MA, Hernandez C, Cisneros-Gonzalez M, Escarela-Perez R. Estimation of synchronous generator parameters using the standstill step-voltage test and a hybrid Genetic Algorithm. Int J Electr Power Energy Syst 2012;35(1):105–11. [7] Zaker B, Gharehpetian GB, Karrari M, Moaddabi N. Simultaneous parameter

References [1] Samadi P, Mohsenian-Rad H, Schober R, et al. Advanced demand side management for the future smart grid using mechanism design. IEEE Trans Smart Grid Sept. 2012;3(3):1170–80. [2] Shahid MU, Mansoor Khan M, Hashmi K, et al. Renewable energy source (RES) based islanded DC microgrid with enhanced resilient control. Int J Electr Power Energy Syst 2019;113:461–71. [3] Xu Y. et al. DGs for service restoration to critical loads in a secondary network. IEEE

9

Electrical Power and Energy Systems 120 (2020) 105989

Z. Wang, et al.

[8] [9] [10] [11] [12] [13] [14] [15] [16]

identification of synchronous generator and excitation system using online measurements. IEEE Trans on Smart Grid May 2016;7(3):1230–8. Hajnoroozi AA, Aminifar F, Ayoubzadeh H. Generating unit model validation and calibration through synchrophasor measurements. IEEE Trans Smart Grid Jan. 2015;6(1):441–9. Tsai C, et al. Practical considerations to calibrate generator model parameters using phasor measurements. IEEE Trans on Smart Grid Sept. 2017;8(5):2228–38. Huang Z, Du P, Kosterev D, Yang S. Generator dynamic model validation and parameter calibration using phasor measurements at the point of connection. IEEE Trans Power Systems May 2013;28(3):1939–49. Meng D, Zhou N, Lu S, Lin G. An expectation-maximization method for calibrating synchronous machine models. Proc. IEEE PES general meeting, Vancouver, BC. 2013. Huang R, et al. Calibrating parameters of power system stability models using advanced ensemble Kalman filter. IEEE Trans Power Systems May 2018;33(3):2895–905. Fan R, Huang R, Diao R. Gaussian mixture model-based Ensemble Kalman filter for machine parameter calibration. IEEE Trans on Energy Conversion Sept. 2018;33(3):1597–9. Fan L. Least squares estimation and Kalman filter based dynamic state and parameter estimation. Proc. IEEE PES General Meeting, Denver, CO. 2015. Valverde G, et al. Nonlinear estimation of synchronous machine parameters using operating data. IEEE Trans Energy Conversion 2011;26(3):831–9. Aghamolki HG, Miao Z, Fan L, Jiang W, Manjure D. Identification of synchronous

[17] [18] [19] [20] [21] [22] [23] [24]

10

generator model with frequency control using unscented Kalman filter. Electric Power Systems Res 2015;126:45–55. Qi J, Sun K, Wang J, Liu H. Dynamic state estimation for multi-machine power system by unscented Kalman filter with enhanced numerical stability. IEEE Trans Smart Grid March 2018;9(2):1184–96. Kundur P, Balu NJ, Lauby MG. Power system stability and control. New York: McGraw-hill; 1994. Mathworks documentation on linear least-squares problems. [Online]. Available: https://www.mathworks.com/help/optim/ug/lsqlin.html. Choi B, Chiang H. Multiple solutions and plateau phenomenon in measurementbased load model development: issues and suggestions. IEEE Trans Power Syst 2009;24(2):824–31. Pan S, Xiao D, Xing S, et al. A general extended Kalman filter for simultaneous estimation of system and unknown inputs. Eng Struct 2016;109:85–98. Ghahremani E, Kamwa I. Dynamic state estimation in power system by applying the extended Kalman filter with unknown inputs to phasor measurements. IEEE Trans Power Syst Nov., 2011.;26(4):2556–66. Anagnostou G, Pal BC. Derivative-free Kalman filtering based approaches to dynamic state estimation for power systems with unknown inputs. IEEE Trans Power Syst Jan. 2018;33(1):116–30. Huang Z, Scheneider K, Nieplocha J. Feasibility studies of applying Kalman Filter Techniques to power system dynamic state estimation. Proc. The 8th international power engineering conference. 2007.