Decentralized Robust Dynamic State Estimation in Power Systems Using Instrument Transformers

This paper proposes a decentralized method for estimation of dynamic states of a power system. The method remains robust to time-synchronization errors and high noise levels in measurements. Robustness of the method has been achieved by incorporating internal angle in the dynamic model used for estimation and by decoupling the estimation process into two stages with continuous updation of measurement-noise variances. Additionally, the proposed estimation method does not need measurements obtained from phasor measurement units; instead, it just requires analog measurements of voltages and currents directly acquired from instrument transformers. This is achieved through statistical signal processing of analog voltages and currents to obtain their magnitudes and frequencies, followed by application of unscented Kalman filtering for nonlinear estimation. The robustness and feasibility of the method have been demonstrated on a benchmark power system model.


Index Terms-Decentralized, time-synchronization error, internal angle, statistical signal processing, dynamic state estimation (DSE), pseudo-input, unscented Kalman filtering (UKF), discrete-time Fourier transform (DFT), Hanning-window, instrument transformers, phasor measurement unit (PMU).
NOMENCLATURE α difference of rotor angle and stator voltage phase in rad. 0 denotes a zero matrix (or vector) of appropriate size. χ denotes a state sigma point. γ denotes a measurement sigma point. g,ḡ discrete and continuous forms of differential functions, resp. h a column vector of the system algebraic functions. I denotes an identity matrix of appropriate size. K the Kalman gain matrix. P denotes a covariance matrix or a cross-covariance matrix. u , v column vectors of pseudo-inputs and process noise, resp.
V , V m analogue stator voltage and its magnitude, resp., in p.u. V d , V q d-axis and q-axis stator voltages, resp., in p.u. V r , V ref AVR-filter voltage and AVR-reference voltage, resp., in p.u. W DFT of Hann window function. X d , X q d-axis and q-axis synchronous reactances, resp., in p.u. X d , X q d-axis and q-axis transient reactances, resp., in p.u. X d , X q d-axis and q-axis subtransient reactances, resp., in p.u. X l armature leakage reactance in p.u. Y denotes a sinusoidal signal with harmonics and noise.

I. INTRODUCTION
A DISTURBANCE in a power system (such as a fault) can initiate spontaneous oscillations in the power-flows in transmission lines. These oscillations grow in magnitude within few seconds if they are undamped or poorly damped. This can lead to loss in synchronism of generators or voltage collapse, ultimately resulting in wide-scale blackouts. The power blackout of August 10, 1996 in the Western Electricity Co-ordination Council region is a famous example of blackouts caused by such oscillations [1], [2]. In order to monitor and control such oscillations and related dynamics which cause instability, the operating state of the system needs to be estimated in real-time, with update rates which are in time scales of ten milliseconds or less (as the time constants associated with such oscillations are not more than ten milliseconds), and this realtime estimation of operating state is known as dynamic state estimation (DSE) [3]- [12]. DSE is a fast growing and widely researched field, and it lays the foundation for a new generation of control methods which can prevent blackouts by dynamically stabilizing the system.
The dynamic states which are estimated and obtained as outputs from DSE algorithms are angles, speeds, voltages and fluxes of the rotors of all the generators in the power system. The inputs which are given to DSE algorithms are some measurable time-varying quantities such as voltage and current of the stator, and some measurable time-invariant quantities such as resistances, reactances, inertia and other constants for the generator. The constant quantities are measured beforehand, and are used as parameters in DSE algorithms.
A generator's voltage, current and power are sinusoidal quantities, and since each sinusoid has a magnitude and a phase (which are together known as a phasor), these quantities can either be represented as sine waves, or as phasors. The conversion of sine waves to phasors is done by phasor measurement units (PMUs). During this conversion, PMUs provide a common reference angle to the phase of the phasor. This is required because a power system is a rotational system (because of the rotational parts of the generators), and every rotational system needs to have a reference angle which is common for all the angles in the system. This common reference angle is provided by PMUs using a common time source for synchronization which is obtained using time clock of global positioning system (GPS) [13]. One important dynamic state which also requires this common reference angle is the rotor angle of a machine. Thus, in order to estimate the rotor angle any DSE algorithm available in power system literature requires synchronized measurements obtained using PMUs [3]- [12].
One problem with time synchronization is that it has associated noise and synchronization-errors [13]. Synchronizationerrors increase the total vector error (TVE) of PMU measurements. For instance, in [14] it was demonstrated that a time synchronization error of just 10 μs can make the TVE greater than 1% for PMUs, even though this synchronization error is much less than the maximum allowable error of 31.6μs as per IEEE standard [15], [16]. As synchronized measurements are used for DSE, these errors can get propagated to the estimated states and deteriorate the overall accuracy and robustness of estimation. It is also not possible to completely eliminate time synchronization as it is inherently required for estimation of rotor angles. This leads to the main idea of the paper: although time synchronization is needed for estimation of rotor angle, it is not needed for estimation of other dynamic states, such as rotor speed, rotor voltages and fluxes, as these states are not defined with respect to a common reference angle. Thus, if the dynamic model which is used for estimation can be modified in such a way that rotor angle is replaced with another angle which does not require time-synchronization, then this can minimize the effects of synchronization on accuracy and robustness of estimation.
The proposed method provides an algorithm for DSE which realizes the above idea. This is done by modifying the estimation model to estimate a relative angle (which does not require synchronization) instead of rotor angle. One such angle is the difference between the rotor angle and the generator terminal voltage phase, also known as the internal angle of the generator. As the rotor angle and the voltage phase have a common reference angle, this reference angle gets cancelled in the difference of the two quantities. Thus, the internal angle, rotor speed, voltages and fluxes can be estimated using the modified estimation model without requiring any synchronized measurements. These dynamic states can then be utilized for decentralized control of the generator [17]- [21]. It should be noted that if the estimation of rotor angle is specifically required then it can be indirectly estimated as the sum of the estimated internal angle and the measured terminal voltage phase obtained using PMU.
To further elaborate, the novel contributions and advantages of the proposed method are enumerated as follows.
r All the dynamic states are estimated without any timesynchronization by incorporating internal angle in estimation model, which in turn ensures robustness of the method to synchronization errors.
r The error in phasor measurements considered in several existing methods of DSE is much less than 1% TVE (see, for example, [3], [4], [9] and [11]). 1% TVE is the permissible error in PMU measurements as per IEEE standard [15], [16], and, hence, these methods of DSE do not consider realistic errors in measurements. The proposed method considers and remains accurate for varying levels of errors in measurements -from 0.1% to 10%. Also none of the currently available methods take into account GPS synchronization errors.
r As synchronization is not required for estimation of the states, DSE for these states can be performed using the analogue measurements directly acquired from current transformers (CTs) and voltage transformers (VTs). This is particularly beneficial for decentralized control purposes, such as in [17]- [21]. r A dual-stage estimation process has been proposed in which interpolated discrete-time Fourier transform (DFT) [22] and unscented Kalman filtering (UKF) [23] have been combined as two stages of estimation. The DFT stage dynamically provides estimates of means and variances of the inputs required by the UKF stage, and this continuous updation of variances is one of the reasons for noiserobustness of the proposed method. In existing methods of DSE for power systems, only static estimates of measurement variances are provided to the estimator.
r Analytical expressions have been obtained for the means and variances of the parameter estimates of a sinusoidal signal (which are given as input to the UKF stage from the DFT stage). Most of these expressions are currently not available in literature. Rest of the paper is organized as follows. Section II specifies the decoupled equations which are used in the proposed method. Section III describes the process for estimation of magnitude, phase and frequency of the analogue signals of terminal voltage and current, while Section IV explains how these estimates can be further used for DSE using unscented Kalman filtering. Section V presents simulations to demonstrate the developed estimation method. Section VI concludes the paper.

II. POWER SYSTEM DYNAMICS IN A DECOUPLED FORM
A power system consists of a wide variety of elements, including generators, their controllers, transmission lines, transformers, relays and loads. All these elements are electrically coupled to each other, and, therefore, in order to define a power system using dynamic mathematical equations, complete knowledge of the models, states and parameters of all these constituent elements is required. Acquiring this knowledge in real-time is not feasible as power systems span wide geographic regions, which are as large as a country, or even a continent. Therefore, it is a practical necessity to represent the dynamic equations of power system in a decoupled form, so that the real-time estimation of dynamic states can be conducted in a decentralized manner. Such a decoupling of system equations can be achieved if a generator and its controller(s) is considered as a decentralized unit, and the stator terminal voltage magnitude, V m , and its phase, θ V , are treated as 'inputs' in the dynamic equations, instead of considering them as algebraic quantities or measurements. This concept of 'pseudo-inputs' for decoupling the equations was introduced in [9], and has been used in the proposed method as well.
In order to estimate the internal angle (which is the difference between the rotor angle and the voltage phase) instead of estimating the rotor angle (as explained in Section I), the decoupled equations and the pseudo-inputs for a generator also get altered. The altered decoupled equations are given by (1)- (11), derived using the subtransient model of machines with four rotor coils in each machine, known as IEEE Model 2.2 [24]. In these equations, the altered pseudo-inputs are V m and voltage frequency, f V , and i refers to the system's ith machine, 1 ≤ i ≤ M . Slow dynamics of the speed-governor have been ignored in this model (although they can also be added, if required). Also, model of a static automatic voltage regulator (AVR) is included with the model of each machine.
The above equations can be written in the following composite state-space form which will be used for DSE (here pseudoinputs are denoted by u i , and the process noise and the noise in pseudo inputs have been included, and denoted by v i and w i , respectively).

III. INTERPOLATED DFT BASED ESTIMATION
Several methods have been proposed in literature for estimating the parameters of a sinusoidal signal, but most of these methods are computationally expensive and, hence, are not suitable for real-time applications [25]. Recently, an interpolated DFT based estimation method was proposed in [22] and was shown to be both fast and accurate enough for real-time control applications in power systems. This method has been further developed in this section for finding the estimates of frequency, magnitude and phase of the fundamental components of measurements obtained from CTs and PTs.
The fundamental component of a sinusoidal signal can be extracted by multiplying the signal with a suitable window function which eliminates other harmonics and higher frequency components in the signal, followed by finding its DFT. One such function is Hanning window function given by h k = sin 2 ( π k N ), and if this function is multiplied with N samples of an analogue signal Y (t) sampled at a frequency f s , then DFT of the product is given by Z(λ) as follows [22].
The key concept in interpolated-DFT based estimation is to approximate W (λ) with the following expression, provided that N >> 1 and λ << N [22], [26].
whereŶ m ,θ andf denote the estimates of Y m , θ and f , respectively. As (16) has three unknowns (which areŶ m ,θ andf ), three distinct equations are required to estimate these unknowns. This can be done by choosing any three distinct values of λ in (16) (say λ = 1, λ = 2 and λ = 3). The obtained values ofŶ m ,θ andf will have associated estimation errors which will depend on N and on the values of λ which are used for generating the three distinct equations. More precisely, these estimation errors are inversely proportional to N 4 [22], and, hence, N should be as large as practically feasible. In this paper N is taken to be in the order of 10 3 , as this is the highest order for N for which interpolated-DFT can run on a state-of-the-art DSP processor without overloading it [22] (The DSP processor used in [22] is a TMS320C6713 with 225 MHz clock rate and 264 KB onchip-RAM. Overloading refers to overall processor usage of above 95%.). Also, for a given N , the estimation errors are minimized if the choices for λ are taken as λ = 0, λ = 1 and λ = 2, provided thatf N f s < 2.1; otherwise, for 2.1 <f N f s < 3, the errors are minimized if the choices are λ = 1, λ = 2 and λ = 3 [22]. The value off N f s should not be greater than 3 as then the delay in obtaining the estimated values becomes too large (that is, more than two cycles, or more than 0.04 s for a 50 Hz power system), and at the same time it should not be too small as then the accuracy of estimation is diminished [22]. In this paper an intermediate value off N f s ≈ 1.5 has been taken and, hence, the former choices of λ = 0, λ = 1 and λ = 2 are applicable.
Remark 1:f N f s is an unknown quantity as f needs to be estimated. But because of power system operational requirements [27], f should remain within 5% of the base system frequency, f 0 (which is usually 50 Hz or 60 Hz), and, hence, if N and f s are chosen such that f 0 N f s = 1.5, thenf N f s ≈ 1.5. The 3 equations which are obtained by putting λ = 0, λ = 1 and λ = 2 in (16) can be written in matrix form as follows. ⎡ Equation (17) implies that the product of a square-matrix and a column vector is equal to a zero vector, when both the matrix and the vector have non-zero elements. This can only happen if the columns of the matrix are linearly dependent, that is, the determinant of the matrix is zero, given as follows.
Simplification of the above determinant givesf as follows.
θ can be obtained by substituting the above value off back into (16) and eliminatingŶ m . To do this, the equation which is obtained by putting λ = 0 in (16) is divided by the equation obtained by putting λ = 1 in (16), which comes as follows.
Solving for e jθ using (20) gives the following expression.
Y m = 8πZ 0 / N Be jθ + Ce −jθ (22) where B, C, and e jθ are given by (20) and (21). Remark 2: It should be noted thatf ,θ andŶ m are real quantities, but they are obtained as functions of complex quantities (given in the right hand sides (RHSs) of (19), (21) and (22), respectively). Hence, these quantities will have negligible but finite imaginary parts associated with them because of finite computational accuracy of any computational device. Thus, during implementation, the imaginary parts should be ignored and only the real parts of RHSs should be assigned tof ,θ orŶ m . Also, asf andŶ m are strictly positive, absolute values of real parts of respective RHSs should be assigned to them.
It was found in [22] that the variance of the above estimate of f in (19) is approximately twice the minimum possible variance which is theoretically achievable using any unbiased estimator (known as Cramer-Rao bound (CRB) [28]). CRB for frequency estimation of a sinusoidal signal has been derived in [28] and is given by CRB(f ) (in Hz 2 ) as follows.
where σ 2 Y is the variance of noise in Y (in p.u.). CRBs forŶ m and θ have been derived in Appendix A, and are given by CRB(Ŷ m ) (in p.u.) and CRB(θ) (in rad 2 ), respectively, as follows.
Following the statistical analysis given in [22] for finding the variance off , the variances ofŶ m andθ are found to be approximately two and six times the above CRBs in (24), respectively; and hence, the estimated variances off ,θ andŶ m are given byσ 2 f (in p.u.),σ 2 Y m (in p.u.) andσ 2 θ (in rad 2 ), respectively, as follows.
where CRB(f ), CRB(Ŷ m ) and CRB(θ) are given by (23) and (24). Estimates of means and variances obtained above are given as inputs to the UKF stage, as detailed in the next section. Remark 3: The advantage of obtaining the analytical expressions forσ 2 f ,σ 2 Y m andσ 2 θ in (23)- (25) is that these variances can be continuously updated and provided to the dynamic estimator (which is the UKF stage) along withf ,θ andŶ m , thereby improving the accuracy of dynamic state estimation.

IV. UNSCENTED KALMAN FILTERING
UKF is a nonlinear method for obtaining dynamic state estimates of a system. It employs the basic idea that performing DSE is easier if the distribution of state estimates is transformed, than if the system model itself is transformed through linearization. System linearization requires computation of Jacobian matrices and is a mathematically challenging task for a high order power system model, especially if it needs to be done at every iteration.
Since linearization is not required in UKF, and, moreover, it has higher accuracy and similar computational speeds as that of linear methods of DSE [8], UKF has been used for performing DSE in this paper.
UKF is a discrete method and, hence, the system given by (12) needs to be discretized before UKF can be applied to it. Discretizing (12) at a sampling period T 0 , by approximatinġ x i with (x ik − x ik )/T 0 , gives the following equation (where k andk represent the kth and (k − 1)th samples, respectively).
In the above model,V ik m andf ik V (found using the DFT method) are used in the pseudo-input vector u ik as follows.
UKF also requires a measurement model besides the above process model. The estimates of active power, P ik e (defined by (10)), and stator current magnitude, I ik m (defined by (11)), which are obtained using the DFT method are used as measurements for UKF. After incorporating the measurement noise, w ik , the measurement model is given as follows. (28) u ik and y ik are estimated quantities and have finite variances which need to be included in the process and measurement models, respectively. This is done by including w ik and w ik in the models as the following zero-mean noises. (30) where P ik w and P ik w denote the covariance matrices of w ik and w ik , respectively. In order to find the estimates and variances in (27) putting Y (t) = I i (t). As P ik e = V ik m I ik m cos (θ ik V − θ ik I ) (from (10)) and the mean values and variances of V ik m , I ik m , θ ik V and θ ik I are known, the mean value of P ik e (denoted asP ik e ) and its estimated variance (denoted asσ 2 P i k e ) can be represented in terms of these known quantities, and have been obtained as follows (here it should be noted that by definition θ ik V and θ ik I lie in the interval (−π/2, π/2], hence, they should be 'unwrapped' by adding or subtracting suitable multiples of π to them, in order to find cos (θ ik V − θ ik I )).
Thus, the four quantities which are required by the UKF stage from the DFT stage are u i , y i , P ik w and P ik w , given by (27)- (30). These quantities should be updated every T 0 s, as this is the sampling period of the UKF stage. Also, in (26), both x ik and w ik are unknown quantities and can be combined together as a composite state vector X ik with a composite covariance matrix P ik X defined as follows.
Here P ik x is the covariance matrix of x ik , and P ik xw is the crosscovariance matrix of x ik and w ik . With the above definition, the model in (26)-(28) is redefined as follows.
With (33) as model and x i0 as steady state estimate of x ik and with the knowledge of g i , h i , u i , y i , P ik w , P ik w and the process noise covariance matrix, P ik v , the filtering equations of UKF for kth iteration and ith unit are given as follows [9].
else reinitializeŵ ik and P ik w in (32) according to (29), leaving rest of the elements inX ik and P ik X unchanged. STEP 2: Generate sigma points

V. CASE STUDY
A model 16-machine, 68-bus benchmark test system (see Fig. 1) has been used for the case study and MATLAB-Simulink (using ode45 solver) running on Windows 7 has been used for its modeling and simulation. A detailed description of the system (including various parameters) is given in [2] or [29]. Static AVRs are used in all the machines, and their parameters are given in [17].
The robust dynamic state estimator (developed in Section II-Section IV) runs at the location of each generation unit, and provides dynamic state estimates for the unit. The measurements which are required by the estimator are V (t) and I(t), and are generated by adding noise to the simulated analogue values of terminal voltage and current of the unit. As explained in Section III, N , f 0 and f s are taken as 1200, 50 Hz and 40000 Hz, respectively. The sampling period of UKF stage, T 0 , is taken as 0.01 s, as explained in [9] or [30]; and thus, the estimates obtained from the DFT stage are also updated every 0.01 s. Also, P ik v is found as described in [30]. For comparison with the proposed estimator, another UKF based dynamic state estimator which uses PMU measurements (given in [9]) also runs at each unit's location and is termed as DSE-with-PMU. Estimate of the internal angle in case of DSE-with-PMU is obtained by subtracting the measurement of terminal voltage phase from the estimate of rotor angle.
The measurement error for the robust DSE method is the percentage error in the analogue signals of V (t) and I(t), while the measurement error for DSE-with-PMU method is the total vector error (TVE) in the phasor measurements of terminal voltage and current. As the measurement errors for the two estimators are of two different kinds, these methods can not be directly compared for same noise levels. Nevertheless, performance of the two methods for standard measurement errors can be compared, as specified by IEEE [15], [16], [31], and IEC [32]. As mentioned in these standards, the measurement error in CTs/VTs should be less than 3%, while the standard error for PMUs is 1% TVE. Hence, in the base case for comparison, the  measurement error for robust DSE is taken as 3%, while for the DSE-with-PMU method, it is taken as 1% TVE.
The system starts from a steady state in the simulation. Then at t = 1 s, a disturbance is created by a three-phase fault at bus 54 and is cleared after 0.18 s by opening of one of the tie-lines between buses 53-54. The simulated states, along with their estimated values for the base case for one of the units (the 13th unit), have been plotted in Figs. 2 and 4. Corresponding estimation errors, which is the difference of estimated and simulated values, have also been plotted in Figs. 3 and 5.    It can be seen in Figs. 2-5 that for robust DSE the plots of estimated values almost coincide with those of the simulated values and the estimation errors are low, but for the DSE-with-PMU method the difference between the simulated and estimated values is apparent and the estimation errors are much higher. This shows that the proposed method performs accurately with standard measurement errors in CTs/VTs, while DSE-with-PMU fails to do so with standard errors in PMU measurements.
Robustness of the proposed method has been tested against varying noise levels in measurements. Fig. 6 shows the estimation results for ω 13 for two more cases: in the first case the noise levels are one-third the base case, while in the second case the noise levels are thrice the base case. Also, root mean squared errors (RMSEs) for varying error levels have been calculated and tabulated in Tables I and II for the two methods. It can be observed that the performance of the proposed method remains robust to errors up to 3%, and even for 10% error, its performance deteriorates only to a small extent. On the other hand, DSE-with-PMU method does not perform accurately for error levels above 0.3% TVE, that is, it is not accurate for 1% TVE and 3% TVE.
The proposed method has also been tested in presence of non-Gaussian noises. This testing has been done by including three different colored noises in the measurements: pink noise, blue noise and violet noise. The estimation results in presence of colored noises for both the proposed method and the DSE-with-PMU method have been presented in Fig. 7, Tables III and IV. It   TABLE II  ROOT MEAN SQUARE ERRORS FOR DSE-WITH-  can be observed from these figure and tables that the proposed method remains robust to non-Gaussian noises as well, while the DSE-with-PMU method gives inaccurate estimation results. It can also be observed that the proposed method has higher estimation errors for pink noise than for blue or violet noises. Computational feasibility of the proposed method can be inferred from the fact that the entire simulation, including simulation of the power system, with two estimators at each machine, runs in faster-than-real-time on MATLAB-Simulink running on Windows 7 on a personal computer with Intel Core 2 Duo, 2.0 GHz CPU and 2 GB RAM. The expression 'fasterthan-real-time' here means that 1 second of the simulation takes less than 1 second of processing time. Also, the total execution time for all the operations for the proposed method for one time step (that is for one iteration) is 0.44 millisecond. Specifically, execution time for the proposed DFT stage is 0.11 ms, while that for UKF is 0.33 ms (for both the proposed method and the PMU method). Thus, the method can be easily implemented using current technologies as the update rate required by the proposed method is 10 milliseconds.

VI. CONCLUSION
A method for dynamic state estimation in power systems has been presented which works using analogue measurements from instrument transformers in order to make the estimation robust to time-synchronization errors. The method is also robust to a wide range of measurement noise which can be encountered in state-of-the-art instrument transformers and has practical computational requirements for real-time operation. This has been achieved using a two-stage estimation algorithm based on interpolated DFT and unscented Kalman filtering. The authors believe that the method will pave the way for fast adoption of methods of dynamic state estimation.

APPENDIX A DERIVATION OF CRAMER RAO BOUND (CRB) FOR PARAMETER ESTIMATION OF A SINUSOIDAL SIGNAL
Let the N samples of a sinusoidal signal Y , sampled at a sampling frequency of f s , be given as follows.
Y k = Y m sin φ k + k , φ k = 2πkf f s +θ, ∀k = 1, 2, . . . , N (34) Here, k is the noise in Y k , and the variance of k is σ 2 Y . The set of parameters which need to be estimated for (34) is Θ = {Θ 1 , Θ 2 , Θ 3 } = {Y m , θ, f}. A lower bound on the variance of any unbiased estimator of Θ is given by the CRB [33], which is found using the inverse of the information matrix, I(Θ). For (34), the (i, j)th element of the matrix I(Θ) is given as follows.