Introduction

Due to the lack of recommended specific algorithms to estimate phasor in I E E E S t d.C37.118, phasor estimation has attracted lots of attentions recently [1]. Phasor estimation is a significant key of wide area monitoring and protecting in power systems. Fast and precise estimation is also necessary for accurate decision in power system control. Dynamic phasor application is not limited to P M U. For example, there are some utilizations in power system simulator programs [2]. Recent developments, particularly the emerging of power electronics based equipment like F A C T S devices, clarified an absence of suitable definition in the typical power system analysis methods which have considered the sinusoidal signal with constant amplitude and phase. For such components (power electronic based components) a full time domain simulation is needed due to incomplete concept of phasor. The concept of time varying phasor (dynamic phasor) has been proposed in [3] for the first time to overcome this problem. This concept has several advantages compared to time-based simulation. For example, it noticeably decreases the simulation time as advantage, but as a disadvantage, increases the number of variables and equations.

Several literatures discussed new algorithms of dynamic phasor estimation. In [4], a new method based on adaptive complex band pass filter was proposed to estimate phasor. Xianing et al. [5] proposed a method based on an angle-shifted energy operated to extract the instantaneous amplitude. An integrated phasor and frequency estimation using a Fast Recursive Gauss Newton algorithm was proposed in [6]. A method based on modified Fourier transform to eliminate DC offset was suggested in [7]. A phasor estimation algorithm based on the least square curve fitting technique was presented in [8] for the distorted secondary current due to CT saturation. In [9], an innovative approach was proposed to estimate the phasor parameters including frequency, magnitude and angle in real time based on a newly constructed recursive wavelet transform. Reference [10] discussed phasor and frequency estimations under transient system conditions: electromagnetic and electro-mechanic. Maximally flat differentiators [11] and phasorlet [12] are other new methods for dynamic phasor estimation. Mai et al., [13]; Serna and Martin [14]; Serna [15] proposed modified forms of earlier methods.

Historically, Guass invented least square method and used it as estimator technique [16]. He suggested that the most appropriate value for the unknown parameter is the most probable one, which is the sum of the square of the observed and the computed values difference. Although Kalman filter is proposed fifty years ago, it is still one of the most important and common data fusion algorithms today. The great success of the Kalman filter is result of its low computational requirement, recursive property and its optimal estimation capability with Gaussian error [17]. The least square and Kalman filter based methods are discussed in this paper, as two general types of phasor estimation. Six specific methods based on these two types have been selected in this study which three of them are based on least square and others are based on Kalman filter.

Method1) Traditional method: This algorithm is based on zeroth-order Taylor expansion and least square to estimate phasor [18].

Method2) Fourier Taylor method: This method is based on second-order Taylor expansion and least square to approximate dynamic phasor [18].

Method3) Shank method: The idea of this method is based on consecutive delays of unit response (digital filter design theory) and least square method to estimate dynamic phasor [19].

Method4) Kalman Taylor method: The main concern of mentioned three methods is delay. In the next three methods, in contrast with the priors, Kalman filter is used as an alternative observer to address the delay challenge [20].

Method5) Fourier Kalman Taylor method: The main idea of this method is based on introducing augment state space which can overcome harmonic infiltration problem [21].

Method6) Modified Kalman Taylor method: The main contribution of this method is to modify modeling process of state space to decrease error bound [22].

The six concepts of algorithms are discussed as different common starting points in a unified manner. The main purpose of this paper is to review and provide a framework in order to compare past and future algorithms in this area.

Dynamic Phasor estimation

Consider a sinusoidal quantity with time-varying amplitude and phase given by:

$$\begin{array}{@{}rcl@{}} S(t)=a(t)cos(2\pi.f_{1}t+\phi(t)) \end{array} $$
(1)

where a(t) and ϕ(t) are amplitude and phase angle of S(t) respectively. f 1 is the frequency of the signal. p(t) is dynamic phasor (complex envelope) that is defined as:

$$\begin{array}{@{}rcl@{}} \ p(t)=a(t)e^{j\phi(t)} \end{array} $$
(2)

By substituting (2) in (1), S(t) can also be written as:

$$\begin{array}{@{}rcl@{}} \ S(t)=\frac{1}{2}\left(p(t)e^{j2.\pi.f_{1}.t}+p^{*}(t)e^{-j2.\pi.f_{1}.t}\right) \end{array} $$
(3)

* represents conjugating operator. In order to estimate dynamic phasor p(t), Taylor series of p(t) at t=0 is used as:

$${} \left\{\begin{array}{l} p(t)=p_{0}+p_{1}t+p_{2}t^2+\ldots+p_{k}t^{k}\\ \left(p_{0}=p(0), p_{1}={p^{\prime}}(0), p_{2}={p^{\prime\prime}}(0)/2,\ldots, p_{k}=p^{k}(0)/k!\right) \end{array}\right. $$
(4)

where the coefficients of the series (P 0, P 1,P 2,..., P k ,) are the derivatives of the dynamic phasor at the observation interval center. All six mentioned methods are similar until this step and differences come to show then.

Method1) traditional method

S(t) can be written based on zeroth-order Taylor polynomial of p(t) as:

$$\begin{array}{@{}rcl@{}} p(t)=\frac{1}{2}\left(p_{0}e^{j2.\pi.f_{1}.t}+p^{*}_{0}e^{-j2.\pi.f_{1}.t}\right) \end{array} $$
(5)

where p 0 and \({p^{*}_{0}}\) are constant term and its conjugated term of Taylor series of p(t). This truncated model can be used in any interval observation like T. The signal S(t) is sampled N 1 times in one period of fundamental frequency (T 1), so interval observation size will be N=(T/T 1)N 1. By substituting N1=2N h +1, (−N h ←n→+N h ) in (5), (6) will be resulted.

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} S(0) \\ \vdots \\ S(N_{h}) \\ \vdots \\ S(n) \\ \vdots \\ S(N-1) \\ \end{array}\right] = \left[ \begin{array}{cc} e^{-j\omega_{1}N_{h}} & e^{j\omega_{1}N_{h}} \\ \vdots & \vdots \\ 1 & 1 \\ \vdots & \vdots\\ e^{\ j\omega_{1}n} & e^{j\omega_{1}n} \\ \vdots & \vdots \\ e^{j\omega_{1}N_{h}} & e^{-j\omega_{1}N_{h}} \\ \end{array}\right].\left(\frac{1}{2}\right) \left[\begin{array}{c} p_{0} \\ p^{*}_{0} \end{array}\right] \end{array} $$
(6)

From left side, first matrix is named S, second matrix is named B (0) and the third one is \({\hat {P}^{(0)}}\) and ω 1=2π/N 1 corresponds to the fundamental angular frequency. The best estimation \({\hat {P}^{(0)}}\), is obtained by least square method as:

$$\begin{array}{@{}rcl@{}} \hat{P}^{(0)}=\left({B^{(0)}}^{H}B^{(0)}\right)^{-1}{B^{(0)}}^{H}.S \end{array} $$
(7)

where H is the Hermitian transpose operator. The estimated time-varying amplitude (\({\hat {a}(t)}\)) and time-varying phase (\({\hat {\phi }(t)}\)) are:

$$\begin{array}{@{}rcl@{}} \left\{\begin{array}{c} \hat{a}(t)=2|\hat{p}_{0}|\\ \hat{\phi}(t)=\angle\hat{p}_{0} \end{array}\right. \end{array} $$
(8)

Method2) Fourier Taylor method

The phasor has been assumed a constant amplitude and phase in previous method which is inappropriate for power system during oscillation like power swing, so time-varying amplitude and time-varying phase are better models in this condition. Based on explained restriction, second method (Fourier Taylor method) has been proposed. Difference between first and second methods is in the usage of higher terms in Taylor expansion. Using first three terms of polynomial in estimation process, S(t) can be written as:

$$\begin{array}{@{}rcl@{}} S(t)&=&\frac{1}{2}\left\{\left(p_{0}+p_{1}t+p_{2}t^{2}\right)e^{j2.\pi.f_{1}.t}\right.\\ &&\quad\left.+\left (p^{*}_{0}+p^{*}_{1}t+p^{*}_{2}t^{2}\right)e^{-j2.\pi.f_{1}.t}\right\} \end{array} $$
(9)

where P 2, P 1, P 0, \({p^{*}_{0}}\), \({P^{*}_{1}}\) and\({ P^{*}_{2}}\) are coefficients of second-order Taylor series and their conjugated, respectively. N linear equations are created as (10):

$$\begin{array}{@{}rcl@{}}{} \hspace{80pt} \left[\begin{array}{c} S(0) \\ \vdots \\ S(N_{h}) \\ \vdots \\ S(n) \\ \vdots \\ S(N-1) \\ \end{array}\right] = \left[\begin{array}{cccccc} N^{2}_{h}e^{-j\omega_{1}N_{h}} & - N_{h}e^{-j\omega_{1}N_{h}} & e^{-j\omega_{1}N_{h}} & e^{j\omega_{1}N_{h}} & -N_{h}e^{j\omega_{1}N_{h}} & N^{2}_{h}e^{j\omega_{1}N_{h}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 1 & 1 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ n^{2}e^{j\omega_{1}n} & ne^{j\omega_{1}n} & e^{\ j\omega_{1}n} & e^{-j\omega_{1}n} & ne^{-j\omega_{1}n} & n^{2}e^{-j\omega_{1}n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ N^{2}_{h}e^{j\omega_{1}N_{h}} & N_{h}e^{j\omega_{1}N_{h}} & e^{j\omega_{1}N_{h}} & e^{-j\omega_{1}N_{h}} & N_{h}e^{-j\omega_{1}N_{h}} & N^{2}_{h}e^{-j\omega_{1}N_{h}} \\ \end{array}\right] \left[\begin{array}{c} p_{2}/2 \\ p_{1}/2 \\ p_{0}/2 \\ p^{*}_{0}/2 \\ p^{*}_{1} /2 \\ p^{*}_{2}/2 \end{array}\right] \end{array} $$
(10)

From left side, first matrix is named S, second matrix is named B (2) and the third one is named \({\hat {P}^{(2)}}\). The best estimation \({\hat {P}^{(2)}}\) is obtained by least square as:

$$\begin{array}{@{}rcl@{}} \hat{P}^{(2)}=({B^{(2)}}^{H}B^{(2)})^{-1}{B^{(2)}}^{H}.S \end{array} $$
(11)

The relationships between the estimated coefficients, time-varying amplitude, time-varying phase and their derivatives, are given by:

$$ \left\{\begin{array}{l} \hat{a}(t)=2|\hat{p}_{0}|\\ \hat{\phi}(t)=\angle\hat{p}_{0}\\ \hat{a}^{\prime}(t)=Re\left\{\hat{p}_{1}e^{-j\hat{\phi}(t)}\right\}\\ \hat{\phi}^{\prime}(t)=\frac{1}{\hat{a}(t)}Img\left\{\hat{p}_{1}e^{-j\hat{\phi}(t)}\right\}\\ \hat{a}^{\prime\prime}(t)=2Re\left\{\hat{p}_{2}e^{-j\hat{\phi}(t)}\right\}+\hat{a}(t)\left[\hat{\phi}^{\prime}(t)\right]^{2}\\ \hat{\phi}^{\prime\prime}(t)=\frac{1}{\hat{a}(t)}2Img\left\{\hat{p}_{2}e^{-j\hat{\phi}(t)}\right\}+\hat{a}^{\prime}(t)\left[\hat{\phi}^{\prime}(t)\right] \end{array}\right. $$
(12)

It is clear from (12) that first and second derivatives of phasor can be calculated by this method.

Method3) Shank method

Digital filters can be directly designed based on least square in Z domain. Shank method is one of these direct filter designs. In this method, the parameters are computed based on the least square criterion. Measurement data are considered as a unit response of digital filter in this method [23]. Just like previous method (Fourier Taylor), S(t) can be written based on second-order Taylor polynomial as (9). Then in discrete time: (sampling time= τ)

$$\begin{array}{@{}rcl@{}} S(t)&=&\frac{1}{2}\left\{\left(\rho_{0}+\rho_{1}.n+\rho_{2}.n^{2}\right)e^{jn\theta_{0}}\right.\\ &&\left.\quad+ \left(\rho^{*}_{0}+\rho^{*}_{1}.n+\rho^{*}_{2}.n^{2}\right)e^{-jn\theta_{0}}\right\} \end{array} $$
(13)
$$\begin{array}{@{}rcl@{}} \rho_{0}=p_{0}, \rho_{1}=p_{1}\tau, \rho_{2}=p_{2}\tau^{2} \end{array} $$

where ρ 2, ρ 1, ρ 0, \({\rho ^{*}_{0}}\), \({\rho ^{*}_{1}}\) and \({ \rho ^{*}_{2}}\) are coefficients of second-order Taylor series and their conjugates in discrete time. While P 2, P 1, P 0, \({p^{*}_{0}}\), \({P^{*}_{1}}\) and \({ P^{*}_{2}}\) are coefficients of second-order Taylor series and their conjugates in continuous time. By applying z transform to truncated Taylor polynomial (13) we have:

$$\begin{array}{@{}rcl@{}} S(z)&=\!&\frac{1}{2}\left\{\left(\frac{\rho_{0}}{1-e^{j\theta_{0}}z^{-1}}\!\right)+\!\left(\frac{\rho^{*}_{0}}{1-e^{-j\theta_{0}}z^{-1}}\!\right)\right.\\ &&\left(\!\frac{\left(\rho_{1}+\rho_{2}/2\right)e^{j\theta_{0}}z^{-1}}{\left(1-e^{j\theta_{0}}z^{-1}\right)^{2}}\right)+\left(\frac{\left(\rho^{*}_{1}+\rho^{*}_{2}/2\right)e^{-j\theta_{0}}z^{-1}}{\left(1-e^{-j\theta_{0}}z^{-1}\right)^{2}}\right)\\ &&\left.\left(\frac{\rho_{2}e^{j2\theta_{0}}z^{-2}}{\left(1-e^{j\theta_{0}}z^{-1}\right)^{3}}\right)+\left(\frac{\rho^{*}_{2}e^{-j2\theta_{0}}z^{-2}}{\left(1-e^{-j\theta_{0}}z^{-1}\right)^{3}}\right)\right\} \end{array} $$
(14)

where z is transformation operator. θ 0=(2.π/N 1) is the sampling angle of fundamental frequency. This can be reduced to rational form by some mathematical operations as:

$$\begin{array}{@{}rcl@{}} S(z)=\frac{\sum_{k=0}^{k=5}b_{k}z^{-k}}{2\left(1-e^{j\theta_{0}}z^{-1}\right)^{3}{\left(1-e^{-j\theta_{0}}z^{-1}\right)^{3}}} \end{array} $$
(15)

According to (15), there are two triple poles at \({e^{j\theta _{0}}}\phantom {\dot {i}\!}\) and \(\phantom {\dot {i}\!}{e^{-j\theta _{0}}}\)and b k coefficients (k=0,1,..,5) include phasor information (ρ 2, ρ 1, ρ 0, \({\rho ^{*}_{0}}\), \({\rho ^{*}_{1}}\) and\({ \rho ^{*}_{2}}\)). This information could be extracted by Shank method. Since poles are determined so locating zeros is the aim of this part. Separating poles from zeros in (15) produces two transfer functions as shown in (16).

$$ \left\{\begin{array}{l} H_{1}(z)=\frac{0.5}{2\left(1-e^{j\theta_0}z^{-1}\right)^{3}\left(1-e^{-j\theta_0}z^{-1}\right)^{3}}\\ H_{2}(z)={\sum\limits_{k=0}^{k=5}b_{k}}z^{-k} \end{array}\right. $$
(16)

Based on Fig. 1 and by considering v(n) as impulse response of H 1 in time domain, (17) is created as:

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} S(0) \\ S(1) \\ \vdots \\ S(N_{1}-1) \\ \end{array}\right] =V. \left[\begin{array}{c} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \\ b_{4} \\ b_{5} \end{array}\right] \end{array} $$
(17)
Fig. 1
figure 1

Poles and zeros separation mentioned in (16)

Full size image

where the left sidematrix is S, the right side matrix is B and the middle one is:

$$\begin{array}{@{}rcl@{}} V= \left[\begin{array}{cccc} \nu(0) & 0 & \cdots & 0 \\ \nu(1) & \nu(0) & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ \nu(N_{1}-1) & \nu(N_{1}-2) & \cdots & \nu(N_{1}-6) \\ \end{array}\right] \end{array} $$

The best estimation of B using least square method is calculated as:

$$\begin{array}{@{}rcl@{}} \hat{B}=\left(V^{H}V\right)^{-1}V^{H}.S \end{array} $$
(18)

The mentioned three methods were based on least square. An important point about least square observer is its delay. It means that dynamic phasor is tracked with delay which will be shown in Section “Simulation results” later. To overcome this problem, next methods (Kalman filter based methods) have been proposed in literatures.

Method4) Kalman Taylor method

Kalman filter is an outstanding method to compute state variables recursively and instantaneously. Regardless to previous methods, next methods are based on state space model and Kalman filter. State space is a complete model for analyzing dynamic system. In this model, state value at each sample time is calculated by its value at previous sample. State space model among with Kalman filter are used to estimate phasor in these methods. The main advantage of Kalman filter based methods is their instantaneous tracking property. State transition matrix can be obtained from the derivatives of p(t). Suppose t 0=(n−1).τ and t=(n).τ, are two consecutive samples where τ is sampling time. The Kth- order Taylor series and its derivatives are:

$$ \left\{ \begin{array}{l} p(t)=p(t_{0})+{p^{\prime}}(t_{0})\tau+{p^{\prime\prime}}(t_{0})\frac{\tau^{2}}{2!}+\ldots+{p^{(k)}}(t_{0})\frac{\tau^{k}}{k!}\\ p^{\prime}(t)=p^{\prime}(t_{0})+{p^{\prime\prime}}(t_{0})\tau+\ldots+{p^{(k)}}(t_{0})\frac{\tau^{(k-1)}}{(k-1)!}\\ \vdots \\ p^{(k)}(t)=p^{(k)}(t_{0}) \end{array}\right. $$
(19)

where p ′(t), p ′′(t),…,p (k)(t) are derivatives of p(t) in time domain. Based on (19), state equations are written as (20):

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} p(t)\\ p^{\prime}(t)\\ p^{\prime\prime}(t)\\ \vdots\\ p^{(k)}(t)\\ \end{array}\right] = \left[\begin{array}{ccccc} 1 & \tau & \frac{\tau^{2}}{2!} & \cdots & \frac{\tau^{k}}{k!} \\ 0 & 1 & \tau & \cdots & \frac{\tau^{(k-1)}}{(k-1)!} \\ 0 & 0 & 1 & \cdots & \frac{\tau^{(k-2)}}{(k-2)!} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 0 & 0 & 0 & 0 & 1 \\ \end{array}\right] \left[\begin{array}{c} p(t_{0})\\ p^{\prime}(t_{0})\\ p^{\prime\prime}(t_{0})\\ \vdots\\ p^{(k)}(t_{0}) \end{array}\right] \end{array} $$
(20)

From left side, the first matrix is named P M(t) which is state vector at time t, the second one named ϕ(τ) which is state transition matrix and the third one named P M(t 0) which is state vector at time t 0. After state equations, measurement equations are obtained. Based on (3) and (20), S(t) is presented by:

$$ \left\{\begin{array}{ll} S(t)=Re\left\{h^T.PM(t).e^{j2\pi.f_{1}t}\right\}=Re\left\{h^Tr(t)\right\}\\ h^{T}=[1~ 0 ~ 0 \cdots 0]\\ r(t)=PM(t).e^{2\pi.f_{1}t} \end{array}\right. $$
(21)

h T=[1 0⋯ 0] is used to extract first component of state vector (P M(t)). As it can be seen in (21), the term \({e^{j.2\pi.f_{1}.t}}\phantom {\dot {i}\!}\) is produced in measurement equations. So new variable r(t), rotated vector, is introduced. The rotated state equations based on r(n) and its conjugates can be written in discrete time as:

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} r(n)\\ r^{*}(n)\\ \end{array}\right] = \left[\begin{array}{cc} \phi(\tau)\varphi_{1} & 0 \\ 0 & \phi(\tau)\varphi^{*}_{1} \\ \end{array}\right] \left[\begin{array}{c} r(n-1) \\ r^{*}(n-1) \\ \end{array}\right] \end{array} $$
(22)

where from left side, the second matrix is named R ϕ(τ) which is rotated state transition matrix and the first and the third matrix are named X(n) and X(n−1) which are rotated state vectors at sample n and n−1 respectively. φ 1 is the sampling factor (\(\phantom {\dot {i}\!}{\varphi _{1} =e^{j\theta _{1}}}\)) of fundamental frequency where θ 1=2.π/N 1 and N 1 is sample number in one fundamental period. S(t) is reconfigured based on r(t) as:

$$\begin{array}{@{}rcl@{}} S(n)=\frac{1}{2}\left[h^{T}h^{T}\right] \left[\begin{array}{cc} r (n) \\ r^{*}(n) \\ \end{array}\right] \end{array} $$
(23)

In (23), from left side, the first matrix is named S(n) and second one is named M S which is final measurement matrix. Finally (22) and (23) are considered as final state and final measurement equations. Kalman filter is applied to these two equations in two steps in order to estimate phasor. Predicting and updating steps are as:

Prediction step:

$$ \left\{\begin{array}{l} X^{-}(n)=R.\phi(\tau).X(n-1)\\ p^{-}(n)=R.\phi(\tau).p(n-1).R.\phi^{H}(\tau)+\Gamma\Gamma^T.\sigma^{2}_{v} \end{array}\right. $$
(24)

where X(n−1) is rotated state vector at (n−1)th sample and X −(n) is its prediction in n th sample. p −(n) is prior error covariance and \({\sigma ^{2}_{v}}\) is the variance of model error. Noise is assumed to affect only the rotated state vector despite of its derivatives thus considered as (h T.h T). In (24) H is the Hermitian transpose operator.

Update step:

$$ \left\{ \begin{array}{l} K(n)=p^{-}(n).MS^T.\left(MS.p^{-}(n).MS^T+\sigma^{2}_{w}\right)\\ X(n)=X^{-}(n)+K(n).\left(S(n)-MS.X^{-}(n)\right)\\ p(n)=(I-K(n).MS).p^{-}(n) \end{array}\right. $$
(25)

where:

X(n) is rotated state vector at sample n.

K(n), Kalman gain, reveals how much modification is needed for state variables based on measurement.

\({\sigma ^{2}_{w}}\) is measurement noise variance created by sensors.

p(n) is posterior error covariance. And

I is the unit matrix.

These Kalman equations make it possible to calculate X(n); Therefore dynamic phasor can be calculated based on estimated value of X(n).

Method5) Fourier Kalman Taylor method

It is clear that Kalman filter works successfully if input signal is matched with the model which the filter is designed based on. In previous method (Kalman Taylor) Kalman filter is designed based on (1) containing only the fundamental frequency. In the cases that input signal is contaminated by harmonic, it is expected that Kalman filter not work properly. So the complete modeling of input signal is necessary to guarantee the accurate operation of filter. Based on mentioned reason, the complete model of main signal is considered as:

$$\begin{array}{@{}rcl@{}} S(t)&=&\{a_{0}(t)cos(2.\pi.f_{0}.t+\phi_{0}(t))\\ &&+a_{1}(t)cos(2.\pi.f_{1}.t+\phi_{1}(t))\\ +&\cdots&+a_{N-1}(t)cos(2.\pi.{(N-1)}f_{1}.t+\phi_{N-1}(t))\} \end{array} $$
(26)

where:

N is sample number in fundamental period.

a 0(t) and ϕ 0(t) are DC amplitude and phase.

a 1(t) and ϕ 1(t) are fundamental amplitude and phase.

a N−1(t) and ϕ N−1(t) are amplitude and phase of (N−1)th harmonic.

f 0 is zero (DC) frequency and f 1 is the fundamental frequency of the signal.

Based on complete model (26), the transition matrix ϕ(τ) is extended to include all harmonics as (27). So it is expected to have individual dynamic phasor for each harmonic (p 0(t)=a 0(t).e j ϕ 0(t)), (p 1(t)=a 1(t).e j ϕ 1(t)), ⋯, (p N−1(t)=a N−1(t).e j ϕ N−1(t)). In this condition, fundamental phasor is free from harmonics that demonstrates the superiority of Fourier Kalman Taylor method compared to Kalman Taylor method.

$$\begin{array}{@{}rcl@{}} \Omega(\tau) = \left[\begin{array}{cccc} \phi(\tau)e^{j0\theta} & 0 &\cdots & 0 \\ 0 & \phi(\tau)e^{j1\theta} &\cdots & 0 \\ 0 & 0 &\cdots & 0 \\ 0 & 0 & 0 & 0 \\ \vdots & \vdots &\cdots & \vdots \\ 0 & 0 &\cdots & \phi(\tau)e^{j(N-1)\theta} \\ \end{array}\right] \end{array} $$
(27)

Rest of dynamic phasor estimation in this method is the same as previous method.

Method6) Modified Kalman Taylor method

Last Kalman based method is obtained by modifying the modeling process of Kalman Taylor method (4th method) in order to decrease error estimation. This model represents a more accurate dynamic of system. As mentioned in (21), the term \(\phantom {\dot {i}\!}{e^{j2\pi.f_{1}.t}}\) has been created in measurement equations. So new variable r(t), rotated vector, has been introduced based on multiplying\(\phantom {\dot {i}\!}{e^{j2\pi.f_{1}.t}}\) by p(t) which produces the rotated state equations written as (22). However this process does not completely express the dynamic of system due to time-varying term \(\phantom {\dot {i}\!}{e^{j2\pi.f_{1}.t}}\). Dynamic behavior of this term (derivatives of this term) has not been considered in the state equation. Therefore, in this method (6th method) consecutive derivatives of \(\phantom {\dot {i}\!}{(e^{j2\pi.f_{1}.t}.p(t))}\) are utilized to produce more accurate state equation. By means of consecutive derivatives we have:

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} r(t) \\ r^{\prime}(t) \\ r^{\prime\prime}(t) \\ \vdots \\ r^{(k)}(t) \\ \end{array}\right] = G. \left[\begin{array}{c} p(t_{0}) \\ p^{\prime}(t_{0}) \\ p^{\prime\prime}(t_{0}) \\ \vdots \\ p^{(k)}(t_{0}) \end{array}\right] \end{array} $$
(28)
$$\begin{array}{@{}rcl@{}} G=e^{j\omega_{1}t} \left[\begin{array}{ccccc} 1 & 0 & 0 & \cdots & 0 \\ j\omega_{1} & 1 & 0 & \cdots & 0 \\ {(j\omega_{1})}^{2} & j\omega_{1} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ {(j\omega_{1})}^{k} & {(j\omega_{1})}^{k-1} & {(j\omega_{1})}^{k-2} & \cdots & 1\\ \end{array}\right] \end{array} $$

where ω 1 is angular frequency of fundamental component. p ′(t),p(t),⋯,p (k)(t) are derivatives of p(t) in time domain and r ′(t),r ″(t),⋯,r (k)(t) are derivatives of r(t). Supposed that t 0=(n−1).τ and t=(n).τ are two consecutive samples and τ is the sampling interval, so:

$$\begin{array}{@{}rcl@{}} R(t)=e^{{jw}_{1}\tau}.G.\phi(\tau).G^{-1}.R(t_{0}) \end{array} $$
(29)

where ϕ(τ) is the matrix described in (20). By considering Q=G.ϕ(τ).G −1, the form of modified state space in discrete time is:

$$\begin{array}{@{}rcl@{}} \left[\begin{array}{c} R(n) \\ R^{*}(n) \\ \end{array}\right] = \left[\begin{array}{cc} e^{{jw}_{1}\tau}.Q & 0 \\ 0 & e^{-{jw}_{1}\tau}.Q^{*} \\ \end{array}\right] \left[\begin{array}{c} R(n-1) \\ R^{*}(n-1) \\ \end{array}\right] \end{array} $$
(30)

Kalman filter is used in method 6 as methods 4 and 5, so the rest of dynamic phasor process is similar to these methods.

Simulation results

First, the test signal which is common in oscillating conditions has been used to examine the proposed methods. Consider the test case as:

$$\begin{array}{@{}rcl@{}} S(t)=a(t)~cos\left(2.\pi.f_{1}.t+\phi(t)\right) \end{array} $$
(31)

where

$$\left\{\begin{array}{l} a(t)=a_{0}+\left(a_{1}cos\left(2.\pi.f_{a}t\right)\right)\\ \phi(t)=\phi_{0}+\left(\phi_{1}cos\left(2.\pi.f_{\phi}.t\right)\right)\\ a_{0}=\phi_{0}=1,a_{1}=\phi_{1}=0.1\\ f_{a}=f_{\phi}=5,N_{1}=16 \end{array}\right. $$

The signal is sampled at 960H z, so 16 samples are obtained over a window of 16.66 m s, which corresponds to one period of the 60 H z. \({\sigma ^{2}_{v}}\) and \({\sigma ^{2}_{w}}\) values are 1×10−2 and 1×10−4 respectively. The oscillation of main signal, shown in Fig. 2, is perceptible around the fundamental frequency. The dynamic phasor is estimated using second order Taylor model in all methods except method 1 which is based on zeroth-order. Figures 3 and 4 show amplitude and phase estimation of dynamic phasor respectively. In these figures the dashed and solid lines represent ideal (real) components and their estimates. According to the figures, it is clear that the main difference between least square based methods (methods 1, 2 and 3) and Kalman based ones (methods 4, 5 and 6) is the estimation delay due to utilization of data window in least square. As a first result, Kalman filter based methods are able to provide instantaneous estimations which are promising result in wide area protection field and synchrophasor application (P M U). An essential attribute of these applications is their synchrony that is provided by Kalman based methods.

Fig. 2
figure 2

Main signal

Full size image
Fig. 3
figure 3

Amplitude estimation

Full size image
Fig. 4
figure 4

Phase estimation

Full size image

As the second result, dynamic phasor concept (methods 2, 3, 4, 5, 6) compared to traditional one (method 1) is more flexible in oscillating conditions. In method 1, a slight distortion appears at estimated amplitude (Fig. 3) and phase (Fig. 4) while this distortion is disappeared in other methods. This improvement is caused due to relaxing amplitude and phase in dynamic phasor model. Total Vector Error (T V E) criterion detects phasor magnitude and angle estimation error, defined as:

$$\begin{array}{@{}rcl@{}} TVE=|\frac{X_{r}-Xe}{X_{r}}| \end{array} $$
(32)

where X r and X e are real and estimated values. Figure 5 depicts the total vector error of all six methods. In order to represent more clearly, first ten cycles has been shown in this figure. These results indicate that the high estimation errors of least square based methods (methods 1, 2 and 3) are mainly due to their one cycle delay. T V E index is not a useful index to compare traditional and dynamic phasor concept because the delay causes high error values in T V E. It is apparent that the low estimation error of Kalman filter based methods (methods 4, 5 and 6) are due to their instantaneous estimations. The value of T V E is achieved approximately 6×10−2 by the method 4. Even though method 5 has been designed to deal with harmonic conditions, this design increases the estimation error. Method 6 provides least error (approximately 4×10−2) which validates the applied modification in modeling process in this method.

Fig. 5
figure 5

Total Vector Error (TVE)

Full size image

Another feature of dynamic phasor concept compared to traditional one is its ability to calculate the derivatives of the phasor. Based on traditional model, it is obvious that estimating the phasor speed and acceleration are impossible by method 1. However, it is possible to obtain estimations of the first and second derivatives of the phasor with the second order Taylor model (all methods except method 1), which are shown in Fig. 6. In these figures the dashed and solid lines represent ideal (real) derivatives and their estimates.

Fig. 6
figure 6

First derivative of phase

Full size image

According to Fig. 6, it is observed that phasor derivative estimations are not as accurate as the phasor estimation (Fig. 4) which indicates the elimination of higher terms in Taylor expansion. These derivatives have two important roles. First, they reduce error estimation as shown in simulation results; Second, they are able to calculate frequency and detect faults and power swings. It is the superiority of dynamic phasor compared to traditional concept of phasor.

In order to clarify this capability, consider a disturbance which occurs in a power system. It is important for us to be discovered immediately to take accurate actions. Power systems make use of distance relays in transmission lines to detect this condition. A distance relay is a device that measures the apparent impedance as an index of distance from the relay location. The power swing is a consequence of a severe disturbance like line fault, loss of generator unit and switching heavy load and creates large fluctuations (just like dynamic phasor condition) of active and reactive power between two areas of a power system. Power swing affects the distance relay behavior and causes its malfunction. Fast detection of power swing is interested in distance protection of transmission lines. Several methods have been proposed to solve this problem till now [24–29]. However the detection based on first and second derivatives of dynamic phasor can be a novel method and makes this aim accessible.

Lack of comprehensive indices to explain the discrepancy of different methods motivated us to establish a framework for comparing presented methods. Twelve indices that are utilized to form a complete benchmark in the paper, are:

  • T V E to examine error bound

  • Step amplitude-phase benchmark tests to analyze dynamic response of the methods in amplitude and phase step condition

  • Step frequency benchmark tests to analyze dynamic response of the methods in frequency step condition

  • Frequency response to demonstrate the delay of the methods

  • Histogram tool to examine RMS error of amplitude estimation

  • Signal taken from a PMU to check presented methods in practical conditions

  • Harmonic and DC offset infiltration

  • Derivatives of amplitude and phase

  • Transient monitor index

  • Computation time

  • Sampling number

  • Noise infiltration

This benchmark is shown in Fig. 7. As mentioned in I E E E standard, the exact algorithm used by P M U in non-steady state condition is beyond of standard scope. However, some simple tests are proposed to evaluate this condition. Two benchmark tests are described in standard as: Magnitude-phase step and Frequency step.

Fig. 7
figure 7

Outline of benchmark for comparison

Full size image

Magnitude-phase step test

To investigate the dynamic response of presented methods, dynamic benchmark based on amplitude-phase step is considered. The test has the form as:

$$ \left\{\begin{array}{ll} a(n)=1, \phi(n)=0 & {0<n<10N_{1}}\\ a(n)=(1+0.9)/2, \phi(n)=\pi/4 &{n=10N_{1}}\\ a(n)=0.9, \phi(n)=\pi/2 & {10N_{1}\leq n} \end{array}\right. $$
(33)

It is 10 % magnitude step and 90° phase step. According to Figs. 8 and 9, the estimated amplitude and phase track their real values accurately after transient period. Method 5 shows the longest transient period which indicates the presence of extremely close poles to unit circle in the z plane among the other methods. Settling time in methods 1 and 2 is twice as method 3 which is dependent on their observation window. The high overshoot value of Kalman based methods is because of their instantaneous behavior which estimate based on previous sample behalf of a samples window. Another reason of this transient response comes from Taylor model which is more appropriate for smooth signals and not sudden changes in signals. This test results show further investigations are needed to improve these transient responses. A possible solution is to add feedback path in observer space state in order to make the dominate poles away from unit circle in z plane.

Fig. 8
figure 8

Amplitude estimation (magnitude and phase step test)

Full size image
Fig. 9
figure 9

Phase estimation (magnitude and phase step test)

Full size image

Frequency step test

The second test waveform is 5 H z frequency step used to evaluate transient response in frequency step condition.

$$ \left\{\begin{array}{ll} S(n)=cos(2\pi.f_{1}\Delta t) & {0<n<10N_{1}}\\ S(n)=cos(2\pi.(f_{1}+5)\Delta t) & {10N_{1}\leq n} \end{array}\right. $$
(34)

Transient responses of magnitude and phase estimation in subjected to frequency step condition are similar to magnitude-phase step condition, which are shown in Figs. 8 and 9 respectively. The main contribution of dynamic phasor can be easily evaluated by this test as shown in Fig. 10. This figure shows the estimation of phase derivative obtained from all six presented methods. The first derivative of the phase is related to frequency (multiplied by 1/(2π)). According to Fig. 10, + 5 H z frequency step is tracked by all methods except method 1. This negative aspect comes from serious limitation of traditional phasor concept which considers phasor as constant amplitude and phase. Therefore the output of this method (method 1) is zero. However, all other methods have frequency and R O C O F (rate of change of frequency) tracking feature. Power system frequency measurement has been in use since the advent of alternating current generator and systems. A number of techniques for measuring power system frequency have been published in technical literatures [30–34]. The frequency estimation of a power signal has wide applications in the control and protection of power systems. Furthermore increasing penetration of distribution generator makes the electrical network to estimate the power signal parameters such as frequency accurately.

Fig. 10
figure 10

Frequency estimation (frequency step test)

Full size image

Frequency measurement is an important parameter in power system operation because it indicates the dynamic balance between power generation and its consumption. To protect power system and detect islanding situation and its time, the frequency and its rate of change are utilized as indicators. Synchronized phasor measurements offer an opportunity to measure power system frequency. Dynamic phasor estimation, presented in this paper, as an inner part of P M U can be utilized in this field.

Frequency response

The frequency analysis illustrates the relationship between input and output in frequency domain. It is done by stimulating the estimator with complex exponential sequence (\(\phantom {\dot {i}\!}e^{jk \omega _{0}}\), where ω 0 is fundamental angular frequency) as:

$$ \left\{ \begin{array}{ll} H(z)=\frac{\hat {P}(z)}{S(z)} &\\ z=e^{jk\omega_0} & {-\infty<k<\infty} \end{array}\right. $$
(35)

where S(z) and \({\hat p(z)}\) are z transforms of input (main signal) and output (estimated dynamic phasor). The frequency responses of all methods are shown in Figs. 11 and 12 (first row includes the magnitude responses and second row includes phase responses). As shown in Figs. 11 and 12, magnitude responses are unsymmetrical. This is explained by the fact that they belong to complex filters. It is apparent that when the input signal corresponds to a constant amplitude and phase signal, all methods work properly with a gain equals to one p.u in magnitude response and group delay equals to zero in phase response.

Fig. 11
figure 11

Frequency response (method 1, 2, 3)

Full size image
Fig. 12
figure 12

Frequency estimation (Frequency response (method 4, 5, 6)

Full size image

However method1 does not work appropriately when the input frequency deviates from nominal frequency even in very little amount. According to Fig. 11 (method1), magnitude response gain decreases sharply before and after fundamental frequency. This frequency domain analysis simply shows the drawback of method1 in oscillating condition. In contrast to method1, other methods show different behaviors around fundamental frequency. The main feature of an appropriate filter is to provide a gain by the value of two around positive and a zero-gain in negative fundamental frequency. This is mainly because the signal model corresponds to two rotator components: first one rotates at positive and second one rotates at negative fundamental frequency. Therefore, complete elimination of negative one and complete pass of positive one is the main duty of the filters (presented methods).

This feature could be accessible by dynamic phasor concept as shown in Figs. 11 and 12. The flatter interval around fundamental frequency leads to lower distortion in phasor estimation. The more persistent flat gain of method 2 compared to method 1 validates the strength of dynamic phasor concept in oscillating conditions. These figures also help to explain the behavior of the estimations when other frequency components are present in the input signal. Methods 1, 2 and 5 provide zero-gain in non fundamental component which demonstrates their ability to remove harmonic in the output.

However methods 3, 4 and 6 do not have this capability. So these methods have difficulties in harmonic conditions especially when subjected to DC component which is common in fault time. In addition, a useful result is attainable based on phase response of output. According to Fig. 11 (second row), all least square based methods (methods 1, 2 and 3) show non-zero around fundamental frequency point which indicates non-zero group delay when subjected to oscillation.

However Kalman filter based methods (methods 4, 5 and 6) show instantaneous estimation property because of their constant zero phase around nominal frequency. The great advantage of these methods is that they can be truly employed in P M U due to their synchronization in nanosecond scale.

It is worthy to note that good performances of all methods except method1, in off nominal frequency are limited around nominal frequency (small vicinity), but because frequency never deviates from the nominal value more than a few m H z in real conditions, the performances of these methods are acceptable. It is also worth noting that, noise effect in dynamic phasor concept is higher than traditional concept which is shown in noise simulation section. This fact can be explained by increase the flat gain length in the interval of fundamental frequency that makes methods more sensitive to noise.

Error bounds

In order to determine the error bound of phasor estimation in all six methods, test signals with variable envelopes are used. 150 amplitude oscillation signals were generated over 10 cycles under the signal model, represented in:

$$\begin{array}{@{}rcl@{}} a(t)=1+0.5\sum_{i=1}^{i=3}e^{-t/\tau{i}}cos(2\pi.f_{1}t) \end{array} $$
(36)

where the time constants (τ i ) were generated by a uniform random process in the interval of [20, 40] cycles. In a similar way, the three frequencies were randomly generated in the intervals of [1, 3], [3, 5] and [5, 7] H z. The error is calculated by:

$$\begin{array}{@{}rcl@{}} rms-error=\sqrt{\sum{(a(n)-\hat{a}(n))}} \end{array} $$
(37)

where a(n) and \({\hat a(n)}\) are real and estimated amplitudes respectively. In statistics, histogram is a graphical representation of the data distribution. It is an estimate of the probability distribution of a continuous variable. Histogram tool creates disjoint categories (known as bins) and counts the number of observations that fall into each of these bins. Thus, by considering N as the total number of observations, K as the total number of bins, and m k as the number of observations in k t h bin, the histogram meets \({\sum _{k=1}^{k=K}}\).In this study, 20 bins and 150 observations have been considered.

Histogram tool distributes the bins along the x-axis between the minimum and maximum values of the error. Thus, by analyzing histogram plot, minimum and maximum errors of methods are attained and the error with the most distribution (the most probable) is also obtainable. Figure 13 shows the histograms of the errors attained by six methods. According to this figure, it is observed that the lowest error is related to method 5 which ranges between 2×10−3 and 8×10−3. Method 3 is also in the second level. The magnitude errors of Kalman based methods (methods 4, 5 and 6) are reasonable compared to least square based methods (methods 1, 2 and 3). The higher errors of least square based methods (methods 1, 2 and 3) correspond to their delays which range between 0.04 and 0.13. It is better to compare all methods without taking account the delay, so the superiority of dynamic phasor concept is more perceivable.

Fig. 13
figure 13

RMS error of amplitude estimation

Full size image

Harmonic infiltration (power system test case)

In real conditions, power system signals may be polluted by harmonics or DC component which is necessary to be examined. In order to analyze oscillated signal along with harmonic, quite identical signal examined in [18] is used in this section. The data were taken from a P M U, installed in a substation. A three-phase fault is created in one of power system lines at 0.1 s e c. The fault is cleared at 0.2 s e c by opening breakers at the ends of faulted line. Line removing causes a swing condition with frequency 5 H z. Measured signal consists of fundamental component (60 H z) and fifth harmonic component (300 H z,1.5 % of fundamental). This signal was sampled at 32 samples per cycle.

All six methods are applied to main signal shown in Fig. 14. Results of estimated amplitude based on six methods are illustrated in Fig. 15. The ideal amplitude components of the phasor are represented in dashed lines and their estimates are shown in solid ones. The discussion on this section is very close to frequency response section and the results are almost similar to those results shown in Figs. 11 and 12. As mentioned before, in frequency response section, least square based methods (methods 1, 2 and 3) create one cycle delay in output due to non-zero flat gain around the fundamental frequency but Kalman filter based methods estimate immediately. These delayed and immediate estimates are clearly shown in Fig. 15. Moreover, this experimental test shows the effect of harmonic component on presented methods. Based on frequency response section, methods 3, 4 and 6 provide non-zero gains in harmonic component. According to Fig. 15, these methods are not able to remove fifth harmonic from the output. It is easy to see that under harmonic condition, input’s harmonics infiltrate to output in methods 3, 4 and 6. However, other methods are able to remove harmonics.

Fig. 14
figure 14

Main signal ([18])

Full size image
Fig. 15
figure 15

Amplitude estimation

Full size image

By analyzing the derivative of amplitude shown in Fig. 16, it is possible to distinguish fault from power swing as explained earlier. This derivative provides beneficial information about power system state. It is obvious that method 1 does not provide this information since it has been built on traditional phasor concept (constant amplitude and phase). Abrupt change at 0.1 sec in estimated derivative obtained by all methods except method 1, demonstrates occurrence of fault in power system. Fault detection process is done by setting appropriate threshold and calculating index in every sample. First derivative of phase and second derivatives of amplitude and phase are also valuable in diagnosis process similar to first derivative of amplitude.

Fig. 16
figure 16

First derivative of phase

Full size image

Noise infiltration

It is useful to analyze the estimation error when noise level changes. A criterion named transient monitor (TM) is used in this section [19]. TM is an index to monitor sudden change of measured signal and is calculated by:

$$ \left\{\begin{array}{l} \left[t_{n}\right]=\left[S_{n}-\hat{S}_{n}\right]\\ TM= \sum\limits_{n=r-N_1}^{n=r}|t_{n}| \end{array}\right. $$
(38)

where S n is the real sample, measured by P M U and \({\hat S_{n}}\) is recomputed sample of S n obtained by dynamic phasor. Estimation error t n , is difference of these two quantities. T M is calculated by estimation error in every sample and can be used as a quality measure of phasor estimation. Figure 17 shows the T M as a function of noise variance. As the noise variance increases, higher errors are obtained by all methods. So there is an upward trend when noise level increases. According to Fig. 17, there is a stable trend for least square based methods (method 1, 2 and 3) before critical point with variance about 10−2.

Fig. 17
figure 17

Transient monitor versus noise variance

Full size image

After this point, transient monitor increases steeply by increasing error variance. Generally Kalman based methods (methods 4, 5 and 6) show lower values compared to least square based methods according to Fig. 17. It is also observed that method 6 provides minimum T M index before critical point. However after critical point (Variance = 10−2) method 1 shows minimum T M index because the dynamic phasor concept is more sensitive to noise due to its flat gain around fundamental frequency. This flat gain allows the noise to be transferred into the output (estimated dynamic phasor).

It is also interesting to analyze the relation between T M behavior and sampling number per cycle which is shown in Fig. 18. Based on the figure, all methods except method 3 generally depict downward trend by increasing the sampling frequency because the estimation accuracy increases by high sampling rate. The different behavior of method 3 can be explained by increasing its noise sensitivity at high sampling rate.

Fig. 18
figure 18

Transient monitor versus sampling frequency

Full size image

Simulation time

All presented methods are examined in computation time when sampling frequency is increased due to its importance in practical cases. This index is utilized to evaluate the capability of methods to be employed off-line or on-line. So this comparison has been done and the results are shown in Fig. 19. A 50-cycle window (computed sample number = 50∗ sample per cycle) is considered in this section and results are obtained by Processor: I n t e l(R) C o r e(T M)2D u o C P U T9550. According to the figure, by increasing the sample number per cycle, the calculation times of C P U move up in all methods. Method 5 (Taylor-Kalman-Fourier) is appropriate to estimate dynamic phasor under harmonic condition, but this matter increases the computation time according to Fig. 19. It is mainly due to its state space model which includes all of harmonic. As it is apparent, the best method in high sampling rate is method 4 and the worst one is method 5.

Fig. 19
figure 19

Computation time versus sampling frequency

Full size image

Conclusion

Phasor and its derivatives are very useful tools to improve the simulation programs of power system and its control. The constant amplitude and phase (traditional definition) impose serious restriction on monitor and control of power system, thus dynamic phasor attracted lots of attentions nowadays. Six dynamic estimation methods have been classified based on least square and Kalman filter; then they have been explained and compared in different test cases. All methods have both advantages and disadvantages. The defect of least square based methods is that they create a delay in estimation process but Kalman filter based methods provide instantaneous estimation. The advantage of all method except 1 is that the speed (first derivative) and acceleration (second derivative) of phasor could be calculated. Lack of comprehensive indices to explain the discrepancy of different method is motivating to establish a framework in order to compare presented methods. Twelve indices are utilized to form a complete benchmark in the paper including: TVE to examine error bound, amplitude-phase step benchmark test to analyze dynamic response of the methods in amplitude and phase step condition, frequency step benchmark test to analyze dynamic response of the methods in frequency step condition, frequency response to demonstrate the delay of the methods, histogram tool to examine RMS error of amplitude estimation, signal taken from a PMU to check presented methods in practical conditions, harmonic and DC offset infiltration, derivatives of amplitude and phase to detect the fault occurrence time, transient monitor index, computation time, sampling number and noise infiltration. Simulation results show that the the dynamic phasor concept gives advantageous results in slow frequency oscillation.