A Robust Algorithm for Real-Time Phasor and Frequency Estimation under Diverse System Conditions

Abstract

This paper presents a comprehensive approach for performing phasor and frequency estimation from voltage and/or current signals of the modern power system. Undesirable components, such as decaying DC, if present in the input signal, are first attenuated using a complex-gain filter. The initial estimates of phasor and frequency are obtained next using the discrete Fourier transform and an improved estimation of signal parameters via rotational invariance technique, respectively. Finally, the accuracy of phasor and frequency estimates are increased based on the identified system condition. Simulations performed to evaluate the proposed approach confirm that it can do fast and accurate estimation of phasor and frequency under diverse operating conditions, making it ideal for wide-area monitoring, protection, and control applications in power systems.

1. Introduction and Background

Optimal energy management in smart grids requires an enhanced visualization of its system parameters, e.g., voltage/current phasors and frequency. This heightened visualization is usually provided by devices capable of producing time-synchronized measurements, such as phasor measurement units (PMUs) in the transmission system and micro-PMUs in the distribution system. To ensure the relevance of the measurements obtained from these devices for monitoring, protection, and control applications, it is necessary that the estimation algorithms used in them are accurate, robust against stray components, computationally efficient, and have low response time [1,2]. Hence, digital signal processing techniques such as discrete Fourier transform (DFT) [3,4,5,6,7,8,9,10,11], least squares (LS) [12,13,14,15], maximum likelihood [16], space vector transform [17], artificial neural networks [18], Hilbert transform [19], Stockwell transform [20], matrix pencil method [21], Kalman filters [22,23], subspace-based methods [24,25], and filter-based methods [26,27] have been proposed recently to estimate phasor and/or frequency under different operating conditions. However, many of the techniques mentioned above suffer from long response time during switching transients [9,13,20], high computational complexity [7,21,24], susceptibility to grid disturbances [12,18,22] and noise [19], lengthy observation window [7,10,11,25,26,27], and performance degradation due to: (i) off-nominal frequency condition [4,5,6,9,12,14,18], and (ii) presence of undesirable components in the input signals, such as decaying DC (DDC) [3,8,15,16,17,23,25]. Thus, there is a genuine need to create robust, real-time phasor and frequency algorithms for synchrophasors as it will ensure a more effective management of the smart grid.

The proposed research satisfies this need and goes beyond the state-of-the-art in the following ways. First, stray components (such as DDC), if present in the input signal, are removed/attenuated by a filter that is designed based on the concept of complex frequencies. The initial values of phasor and frequency are then estimated in parallel from the filtered signal using DFT and an improved version of estimation of signal parameters via rotational invariance technique (ESPRIT) [28], called iESPRIT that significantly reduces the computational burden. Finally, the initial estimates of phasor and frequency are updated based on the identified system condition (normal/abnormal). Performance evaluation of the proposed approach using different test signals indicates that it makes the following salient contributions:

• It efficiently suppresses undesirable components present in the input signal, and can therefore be applied to both voltage and current signals.
• It provides fast, accurate, and robust phasor and frequency estimates under diverse operating conditions (faults, frequency ramps, power swings, harmonics).
• It detects abnormal system conditions quickly (in less than one-cycle).
• It satisfies the requirements of relevant IEEE and IEC Standards (this paper focuses on P-class PMUs which have more stringent speed requirements) [29,30].

2. Proposed Filter for DDC Rejection

Considering the most general setting, let the input signal be described by

$$x\left(t\right)=A_0+\sum _{k=1}^K{A}_{k}cos({\omega }_{k}t+{\theta }_k)+D_0{e}^{-t/\tau }$$

(1)

where $A_0$ denotes a DC offset; K denotes the total number of sinusoidal components; $A_k$ and ${\theta }_k$ denote the magnitude and phase angle of the k-th sinusoid. The angular frequency of the k-th sinusoid is denoted by ${\omega }_k$ and is equal to $2\pi f_k$, where $f_k$ is the frequency of the k-th sinusoid; the nominal frequency and angular frequency are denoted by $f_{nom}$ and ${\omega }_{nom}$, respectively. $D_0$ and $\tau$ are the magnitude and time constant of the DDC component and are called DDC parameters. The above signal can be written in discretized form as

$$x\left[n\right]=A_0+\sum _{k=1}^K{A}_{k}cos({\omega }_{k}n\Delta t+{\theta }_k)+D_0{E}^n$$

(2)

in which $\Delta t$ is the sampling interval and is equal to $1/f_s$ (where $f_s$ is the sampling frequency), n denotes the sample number, and $E\stackrel{\Delta }{=}{e}^{-\Delta t/\tau }$, called the damping factor, depends on the DDC time constant [4]. It is difficult to predetermine the DDC parameters as they depend on fault characteristics, such as, fault inception time/angle, fault location, and fault resistance [4,5,6,31]. However, their presence (usually in the current signals) can cause up to 15% error in the estimated phasor [31]. This can result in a maximum overreach of 115% if the first zone of a distance relay is set at 100% of the line impedance. Thus, DDC suppression is important for safe and reliable operation of the smart grid.

The goals of the proposed filter design are: (i) removing/attenuating the DDC for different values of its time constant, (ii) preventing variation in parameters of the fundamental component during the removal process, and (iii) attaining fast response time (or equivalently low-order characteristic). The frequencies present in (1) are shown in Figure 1, where the frequencies of the DC offset, DDC, and sinusoidal components are located at the origin, negative side of the real axis, and imaginary axis of the s-plane, respectively. The DDC time constant ($\tau$) varies due to fault characteristics (such as fault impedance and fault location) between 0.5 and 5 cycles [31,32]; hence, a corresponding lower bound (${\tau }_l$) and upper bound (${\tau }_u$) are specified. Since the DDC frequency is equal to $-1/\tau$ in the s-plane, the zeros of the proposed filter should lie between $[-1/{\tau }_l,-1/{\tau }_u]$, which is called DDC frequency band (see Figure 1), to filter out the DDC. However, the filter’s zeros will cause phase shifting of the fundamental component. To compensate for this phase-shift, filter’s poles can be added (as was done in [6]). However, doing so will result in: (i) neutralization of the effect of DDC attenuation by the filter’s zeros, (ii) additional error in phasor estimation due to the vicinity of the added poles to the fundamental frequency, and (iii) oscillatory behavior in the estimated magnitude and phase angle. To solve this problem, in this paper, we propose the design of a complex gain filter. The transfer function of the proposed filter is given by

$$G\left(s\right)=K_F\prod _{m=1}^R\left(s-z_m\right)$$

(3)

where the total number of zeros (or equivalently the filter order) is denoted by R which plays a key role in the response time of the proposed filter. The filter gain is expressed by $K_F$ which is complex (i.e., $K_F\stackrel{\Delta }{=}\left|K_F\right|e^{j{\gamma }_F}$). $z_m$s are inserted in the DDC frequency band with the aim of reaching minimum multiplication of their distance to any frequency in this band. This can be formulated as the following optimization problem:

$$\begin{array}{cc}\hfill & \underset{z_m}{min}\prod _{m=1}^R|x-z_m{|}^2\forall x\in \left[-1/{\tau }_l,-1/{\tau }_u\right]\hfill \\ \hfill & \mathrm{subject}\mathrm{to}:-1/{\tau }_l<z_m<-1/{\tau }_u\forall m\in \{1,\cdots ,R\}\hfill \end{array}$$

(4)

Note that this optimization problem ensures that the proposed filter suppresses DDC components regardless of their time-constant values, without any re-tuning.

Next, to compensate for the phase shift due to the added zeros, the phase-angle of the filter gain ${\gamma }_F$ is adjusted as follows (see Figure 1):

$${\gamma }_F+\sum _{m=1}^R{\delta }_m=0$$

(5)

Finally, the magnitude of the filter gain $|K_F|$ is determined with the aim of reaching unit total filter gain at nominal frequency $f_{nom}$. This is given by

$$\left|G\right({\omega }_{nom}=2\pi f_{nom}\left)\right|=1$$

(6)

Since the proposed filter is designed for ${\omega }_{nom}$, its performance may deteriorate at off-nominal frequencies. Under such circumstances, the filter’s parameters can be readjusted as described in Section 6.2. As the number of filter’s zeros increases, the filter can suppress DDC components more effectively, at the cost of increased response time. This trade-off between the proposed-filter’s capacity for DDC attenuation and its response time is investigated in Section 7.1.

3. Phasor Estimation at Off-Nominal Frequencies

DFT is one of the most widely used techniques for phasor estimation due to its simplicity in implementation, low computational burden, harmonic filtering, and noise alleviation [1]. By applying DFT on the filtered signal, denoted by subscript f, we get

$${\acute{\hat{X}}}_1=\frac{2}{N}\sum _{n=r}^{r+N-1}{x}_f\left[n\right]e^{-j\frac{2\pi }{N}n}=\frac{2}{N}\sum _{n=r}^{r+N-1}{x}_f\left[n\right]e^{-jn{\omega }_{nom}\Delta t}$$

(7)

where r is the number of the first sample in the data string, N denotes the number of samples per cycle, and ${\acute{\widehat{X}}}_1$ refers to estimated fundamental phasor at nominal frequency. Note that the DC offset and harmonic components, present in (2), are filtered out by the DFT, while the DDC has already been attenuated by the proposed filter (see Section 2). However, when the frequency deviates from its nominal value $f_{nom}$, the phasor estimated using DFT can be erroneous. This requires additional correction, as illustrated below.

Consider an off-nominal frequency condition in which the angular frequency is different from its nominal value ${\omega }_{nom}$; i.e., let the actual angular frequency be ${\omega }_1$. The DFT of the filtered input signal in this situation can be written as

$${\widehat{X}}_1=\frac{2}{N}\sum _{n=r}^{r+N-1}{x}_f\left[n\right]e^{-jn{\omega }_1\Delta t}$$

(8)

where ${\widehat{X}}_1$ is the estimated phasor by DFT at actual frequency ${\omega }_1$. After performing trigonometric simplifications on (7) and (8), we get

$${\acute{\widehat{X}}}_1=P{\widehat{X}}_1{e}^{jr({\omega }_1-{\omega }_{nom})\Delta t}+Q{\widehat{X}}_1^{*}{e}^{-jr({\omega }_1+{\omega }_{nom})\Delta t}$$

(9)

$$P\stackrel{\Delta }{=}\frac{sin\left(\frac{N({\omega }_1-{\omega }_{nom})\Delta t}{2}\right)}{Nsin\left(\frac{({\omega }_1-{\omega }_{nom})\Delta t}{2}\right)}{e}^{j(N-1)\frac{({\omega }_1-{\omega }_{nom})\Delta t}{2}}$$

(10)

$$Q\stackrel{\Delta }{=}\frac{sin\left(\frac{N({\omega }_1+{\omega }_{nom})\Delta t}{2}\right)}{Nsin\left(\frac{({\omega }_1+{\omega }_{nom})\Delta t}{2}\right)}{e}^{-j(N-1)\frac{({\omega }_1+{\omega }_{nom})\Delta t}{2}}$$

(11)

where P and Q are correction coefficients that are independent of ‘r’ and dependent on N and frequency deviation $\Delta f$. Variations of P and Q for different values of N and $\Delta f$ are shown in

Table A1. Magnitudes and phase angles of correction coefficients for different frequency deviations and sampling rate of 24.

$\mathbf{\Delta }\mathit{f}$ (Hz)	|P|	$\mathbf{\angle }\mathit{P}$ (${}^{\mathbf{\circ }}$)	|Q|	$\mathbf{\angle }\mathit{Q}$ (${}^{\mathbf{\circ }}$)
−5	0.9886	−14.37	0.0434	29.37
−4.5	0.9908	−12.94	0.0390	27.94
−4	0.9927	−11.5	0.0346	26.5
−3.5	0.9944	−10.06	0.0302	25.06
−3	0.9959	−8.62	0.0258	23.62
−2.5	0.9972	−7.19	0.0215	22.19
−2	0.9982	−5.75	0.0171	20.75
−1.5	0.9990	−4.31	0.0128	19.31
−1	0.9995	−2.87	0.0085	17.87
−0.5	0.9999	−1.44	0.0042	16.44
0	1.0000	0	0	15
0.5	0.9999	1.44	0.0042	13.56
1	0.9995	2.87	0.0084	12.12
1.5	0.9990	4.31	0.0125	10.69
2	0.9982	5.75	0.0166	9.25
2.5	0.9972	7.19	0.0206	7.81
3	0.9959	8.62	0.0246	6.37
3.5	0.9944	10.06	0.0285	4.94
4	0.9927	11.5	0.0324	3.5
4.5	0.9908	12.94	0.0363	2.06
5	0.9886	14.37	0.0400	0.62

and

Table A2. Magnitudes and phase angles of correction coefficients for different sampling rates and frequency of 62 (Hz).

N	|P|	$\mathbf{\angle }\mathit{P}$ (${}^{\mathbf{\circ }}$)	|Q|	$\mathbf{\angle }\mathit{Q}$ (${}^{\mathbf{\circ }}$)
12	0.9982	5.5	0.0172	24.5
24	0.9982	5.75	0.0166	9.25
36	0.9982	5.83	0.0164	4.17
48	0.9982	5.876	0.0164	1.62
60	0.9982	5.9	0.0164	0.1
72	0.9982	5.92	0.0164	−0.92
84	0.9982	5.93	0.0164	−1.64
96	0.9982	5.94	0.0164	−2.19
108	0.9982	5.947	0.0164	−2.61
120	0.9982	5.95	0.0164	−2.95

in Appendix A.1. In a typical power system, the variations in the frequency rarely exceed $0.5$ Hz in either direction from its nominal value [33]. As shown in Table A1 and Appendix A.1, Q is negligibly small in comparison to P in the range of $-0.5$ Hz and $+0.5$ Hz. The impact of Q on the accuracy of the estimated phasor by the proposed algorithm is investigated in more detail in Appendix A.2. From the conducted investigation, it is realized that the impact of Q on the accuracy of the estimated phasor by the proposed algorithm is small (in comparison to the total vector error (TVE) limit of 1% [29,30]). Consequently, (9) is simplified and rewritten as follows:

$${\acute{\widehat{X}}}_1\approxeq \underset{\acute{P}}{\underbrace{P{e}^{jr({\omega }_1-{\omega }_{nom})\Delta t}}}{\widehat{X}}_1=\acute{P}{\widehat{X}}_1$$

(12)

Since $\acute{P}$ depends on the actual frequency, it is necessary to estimate the frequency precisely and rapidly for correcting the phasor estimate. The frequency estimation algorithm described in Section 4 can not only improve accuracy of the estimated phasor at off-nominal frequency conditions (see Section 6.2), but also estimate the frequency quickly and accurately for other power system applications.

4. Proposed Frequency Estimation Algorithm

ESPRIT is a well-known subspace-based technique for frequency estimation due to its strong immunity against noise [28]. However, due to its high computational complexity [25], ESPRIT is not suitable for power system applications that require fast estimation of frequency. Hence, an improved version of the ESPRIT, called iESPRIT, is proposed in this paper. The input signal for frequency estimation is the output of the filter that was proposed in Section 2, and can be expressed as

$$x_f\left[n\right]=A_{0f}+A_{1}cos({\omega }_{1}n\Delta t+{\theta }_1)+\sum _{k=2}^K{A}_{kf}cos({\omega }_{k}n\Delta t+{\theta }_{kf})+w_f\left[n\right]$$

(13)

where $w_f\left[n\right]$ denotes an additive noise component with the power of ${{\sigma }_n}^2$. Note that the filter was designed in such a way that it kept the parameters of the fundamental component unchanged. Equation (13) can be equivalently written in the form of $Y\stackrel{\Delta }{=}2K+1$ complex exponentials distorted with noise, called the harmonic model, as follows:

$$\begin{array}{cc}\hfill x_f\left[n\right]& =A_{0f}+\frac{1}{2}{A}_1\left(e^{j\left({\omega }_{1}n\Delta t+{\theta }_1\right)}\right)+e^{-j\left({\omega }_{1}n\Delta t+{\theta }_1\right)}\hfill \\ & +\sum _{k=2}^K\frac{1}{2}{A}_{kf}\left(e^{j\left({\omega }_{k}n\Delta t+{\theta }_{kf}\right)}+e^{-j\left({\omega }_{k}n\Delta t+{\theta }_{kf}\right)}\right)+w_f\left[n\right]=\sum _{k=-K}^K{\alpha }_k{e}^{j{\omega }_{k}n\Delta t}+w_f\left[n\right]\hfill \end{array}$$

(14)

in which ${\alpha }_0=A_{0f}$ and ${\omega }_0=0$; ${\alpha }_1={\alpha }_{-1}^{*}\stackrel{\Delta }{=}\left|{\alpha }_1\right|e^{j{\phi }_1}=\frac{1}{2}{A}_1{e}^{j{\theta }_1}$ and ${\omega }_1=-{\omega }_{-1}$; ${\alpha }_k={\alpha }_{-k}^{*}\stackrel{\Delta }{=}\left|{\alpha }_k\right|e^{j{\phi }_k}=\frac{1}{2}{A}_{kf}{e}^{j{\theta }_{kf}}$ and ${\omega }_k=-{\omega }_{-k}$ for $k=2,\cdots ,K$. For simplicity, we denote the Y${\alpha }_k$s and ${\omega }_k$s for $k=-K,\cdots ,0,\cdots ,K$ with ${\alpha }_1$, ${\alpha }_2$, ⋯, ${\alpha }_Y$ and ${\omega }_1$, ${\omega }_2$, ⋯, ${\omega }_Y$, respectively, and let ${\overline{\omega }}_k\stackrel{\Delta }{=}({\omega }_k/f_s)$, where ${\overline{\omega }}_k$ is thenormalized angular frequency of the k-th exponential component. Then, by using M consecutive samples of $x_f\left[n\right]$, an M-length data vector ${\mathbf{x}}_f\left(n\right)$ is defined as follows:

$$\begin{array}{cc}\hfill & {\mathbf{x}}_f\left(n\right)\stackrel{\Delta }{=}\left[\begin{array}{c}{x}_f\left[n\right]\\ x_f[n+1]\\ ⋮\\ x_f[n+M-1]\end{array}\right]=\underset{{\mathbf{H}}_{M\times Y}}{\underbrace{\left[\begin{array}{ccc}1& \cdots & 1\\ e^{j{\overline{\omega }}_1}& \cdots & e^{j{\overline{\omega }}_Y}\\ ⋮& \ddots & ⋮\\ e^{j(M-1){\overline{\omega }}_1}& \cdots & e^{j(M-1){\overline{\omega }}_Y}\end{array}\right]}}{\Phi }^n\underset{\alpha }{\underbrace{\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_Y\end{array}\right]}}\hfill \\ \hfill & +\underset{{\mathbf{w}}_f\left(n\right)}{\underbrace{\left[\begin{array}{c}{w}_f\left[n\right]\\ w_f[n+1]\\ ⋮\\ w_f[n+M-1]\end{array}\right]}}=\mathbf{H}{\Phi }^n\alpha +{\mathbf{w}}_f\left(n\right)\hfill \\ \hfill & {\Phi }_{Y\times Y}=\mathrm{diag}\left[e^{j{\overline{\omega }}_1},e^{j{\overline{\omega }}_2},\cdots ,e^{j{\overline{\omega }}_Y}\right]\hfill \end{array}$$

(15)

where the Y columns of matrix $\mathbf{H}$ are length-M signal frequency vectors described by

$$\mathbf{h}\left({\overline{\omega }}_k\right)\stackrel{\Delta }{=}{\left[1e^{j{\overline{\omega }}_k}\cdots e^{j(M-1){\overline{\omega }}_k}\right]}^T$$

(16)

for $k=1,\cdots ,Y$. The matrix $\mathbf{H}$ is simplified using (16) as

$${\mathbf{H}}_{M\times Y}=\left[\mathbf{h}\left({\overline{\omega }}_1\right)\mathbf{h}\left({\overline{\omega }}_2\right)\cdots \mathbf{h}\left({\overline{\omega }}_Y\right)\right]$$

(17)

In addition, matrix $\Phi$ is the diagonal matrix of phase shifts and can be expressed as $\Phi =\mathrm{diag}\left[{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_Y\right]$, where ${\varphi }_k=e^{j{\overline{\omega }}_k}$ for $k=1,2,\cdots ,Y$. Since the frequencies of the complex exponentials ${\overline{\omega }}_k$s completely describe the matrix $\Phi$, the estimated frequencies can be obtained by finding $\Phi$.

Let us consider two overlapping sub-matrices within the length-M window vector. In this regard, an $(M-1)\times Y$ matrix ${\mathbf{H}}_{M-1}$, which is the length-($M-1$) sub-matrix of $\mathbf{H}$, is defined as

$${\mathbf{H}}_{M-1}=\left[{\mathbf{h}}_{M-1}\left({\overline{\omega }}_1\right){\mathbf{h}}_{M-1}\left({\overline{\omega }}_2\right)\cdots {\mathbf{h}}_{M-1}\left({\overline{\omega }}_Y\right)\right]$$

(18)

in which ${\mathbf{h}}_{M-1}\left({\overline{\omega }}_k\right)\stackrel{\Delta }{=}{\left[1e^{j{\overline{\omega }}_k}\cdots e^{j(M-2){\overline{\omega }}_k}\right]}^T$ for $k=1,2,\cdots ,Y$. Then, two new $(M-1)\times Y$ matrices which are called ${\mathbf{H}}_1$ and ${\mathbf{H}}_2$ are defined as follows:

$${\mathbf{H}}_1\stackrel{\Delta }{=}{\mathbf{H}}_{M-1}{\Phi }^n\mathrm{and}{\mathbf{H}}_2\stackrel{\Delta }{=}{\mathbf{H}}_{M-1}{\Phi }^{n+1}$$

(19)

Using (15) and (19), we have:

$$\mathbf{H}{\Phi }^n=\left[\begin{array}{ccc}& {\mathbf{H}}_1& \\ *& *& *\end{array}\right]=\left[\begin{array}{ccc}*& *& *\\ & {\mathbf{H}}_2& \end{array}\right]$$

(20)

As seen in (19), the matrices ${\mathbf{H}}_1$ and ${\mathbf{H}}_2$ are related as

$${\mathbf{H}}_2={\mathbf{H}}_1\Phi$$

(21)

Now, a $G\times M$ matrix, ${\mathbf{X}}_f$, called the data matrix, can be obtained from an L-length of input data

$${\mathbf{X}}_f\stackrel{\Delta }{=}\left[\begin{array}{cccc}{x}_f\left[1\right]& x_f\left[2\right]& \cdots & x_f\left[M\right]\\ x_f\left[2\right]& x_f\left[3\right]& \cdots & x_f[M+1]\\ ⋮& ⋮& \ddots & ⋮\\ x_f\left[G\right]& x_f[G+1]& \cdots & x_f[G+M-1]\end{array}\right]$$

(22)

in which $L\stackrel{\Delta }{=}G+M-1$ and G is the total number of data vectors used to construct ${\mathbf{X}}_f$. In conventional ESPRIT, $\Phi$ is calculated from the singular value decomposition (SVD) of ${\mathbf{X}}_f$, as shown below [28]

$${\mathbf{X}}_f=\mathbf{L}\Sigma {\mathbf{U}}^H$$

(23)

where ${\mathbf{L}}_{G\times G}$ and ${\mathbf{U}}_{M\times M}$ are the left and right singular vectors (both of them are unitary), and ${\Sigma }_{G\times M}$ is a diagonal matrix containing singular values ordered in descending magnitude. However, SVD of ${\mathbf{X}}_f$ will have a complexity of $O(G{M}^2+M{G}^2+M^3)$[34], which makes ESPRIT too slow for frequency estimation in P-class PMUs [1,29,30]. To overcome this problem, the proposed iESPRIT estimates an autocorrelation matrix from ${\mathbf{X}}_f$ as shown below

$${\widehat{\mathbf{R}}}_x=\frac{1}{G}{\mathbf{X}}_f^H{\mathbf{X}}_f$$

(24)

where ${\widehat{\mathbf{R}}}_x$ can be expressed in terms of its eigenvalues and eigenvectors using eigen-decomposition (ED) as

$${\widehat{\mathbf{R}}}_x=\mathbf{Q}\Lambda {\mathbf{Q}}^H$$

(25)

where $\Lambda$ is an $M\times M$ diagonal matrix containing the eigenvalues of the ${\widehat{\mathbf{R}}}_x$ in descending order. The columns of $\mathbf{Q}$ are also the corresponding eigenvectors. It is proved in Appendix B that the eigenvectors of ${\widehat{\mathbf{R}}}_x$ (columns of $\mathbf{Q}$) are equal to the columns of $\mathbf{U}$ which is used in ESPRIT (see (23)). Note that the complexity of ED of ${\widehat{\mathbf{R}}}_x$ is $O\left(M^3\right)$[34], which is much lower than the complexity of SVD [35]. The Y largest diagonal elements of $\Lambda$ (i.e., ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_Y$) are equal to $M{|{\alpha }_k|}^2+{{\sigma }_n}^2$ (for $k=1,2,\cdots ,Y$), while the remaining $M-Y+1$ diagonal elements of $\Lambda$ (i.e., ${\lambda }_{Y+1},\cdots ,{\lambda }_M$) have similar values as ${{\sigma }_n}^2$[28]. As a result, the total number of complex exponentials Y can be obtained using the criterion of the relative slope variation between the ${\widehat{\mathbf{R}}}_x$’s eigenvalues as shown in Figure 2. Note that the fundamental frequency corresponds to the largest singular value in $\Lambda$.

In order to find $\Phi$, the matrix $\mathbf{Q}$ is split into two sub-matrices for separation of the signal and noise subspaces

$${\mathbf{Q}}_{M\times M}=\left[{\mathbf{Q}}_S\mid {\mathbf{Q}}_n\right]$$

(26)

where ${\mathbf{Q}}_S$ is an $M\times Y$ sub-matrix that corresponds to the Y largest magnitudes of the eigenvalues (see (25)). Since the elements of the matrix $\mathbf{H}$ defined in (17) completely describe the frequencies of exponential components existing in the input signal, $\mathbf{H}$ similar to ${\mathbf{Q}}_S$ must also lie in the same subspace, called the signal subspace. As a result, an invertible transformation $\mathbf{W}$ exists that maps ${\mathbf{Q}}_S$ into $\mathbf{H}$ as

$$\mathbf{H}={\mathbf{Q}}_S\mathbf{W}$$

(27)

Note that it is not necessary to calculate the transformation $\mathbf{W}$, and it is only formulated as a mapping between the matrices $\mathbf{H}$ and ${\mathbf{Q}}_S$ within the signal subspace. The matrix ${\mathbf{Q}}_S$ can also be split into two smaller $M-1$-dimensional subspaces as was done with the matrix $\mathbf{H}$

$${\mathbf{Q}}_S=\left[\begin{array}{ccc}& {\mathbf{Q}}_1& \\ *& *& *\end{array}\right]=\left[\begin{array}{ccc}*& *& *\\ & {\mathbf{Q}}_2& \end{array}\right]$$

(28)

where both ${\mathbf{Q}}_1$ and ${\mathbf{Q}}_2$ are $(M-1)\times Y$ sub-matrices. Since ${\mathbf{H}}_1$ and ${\mathbf{H}}_2$ correspond to the same subspaces as ${\mathbf{Q}}_1$ and ${\mathbf{Q}}_2$, (27) must also hold for these subspaces, as shown below:

$${\mathbf{H}}_1={\mathbf{Q}}_1\mathbf{W}\mathrm{and}{\mathbf{H}}_2={\mathbf{Q}}_2\mathbf{W}$$

(29)

As seen in (21), the rotational matrix $\Phi$ exists that relates (rotates) ${\mathbf{H}}_1$ to ${\mathbf{H}}_2$. In the same way, a rotational matrix (denoted by $\Psi$ in this paper) must also exist that relates (rotates) ${\mathbf{Q}}_1$ to ${\mathbf{Q}}_2$ as

$${\mathbf{Q}}_2={\mathbf{Q}}_1\Psi$$

(30)

Note that the matrices ${\mathbf{Q}}_1$ and ${\mathbf{Q}}_2$ are known from the ED of ${\widehat{\mathbf{R}}}_x$. Therefore, from (30), the rotational matrix $\Psi$ can be estimated using LS as

$$\widehat{\Psi }={\left({\mathbf{Q}}_1^H{\mathbf{Q}}_1\right)}^{-1}{\mathbf{Q}}_1^H{\mathbf{Q}}_2$$

(31)

Finally, using (21), (29) and (30), we get

$$\widehat{\Phi }={\mathbf{W}}^{-1}\widehat{\Psi }\mathbf{W}$$

(32)

Now, as the diagonal entries of $\widehat{\Phi }$ (${\widehat{\varphi }}_k=e^{j{\widehat{\overline{\omega }}}_k}$ for $k=1,2,\cdots ,Y$) are simply the eigenvalues of $\widehat{\Psi }$, the estimated frequencies can be obtained as

$${\widehat{\overline{\omega }}}_k=\angle {\widehat{\varphi }}_k\mathrm{for}k=1,2,\cdots ,Y$$

(33)

Recall ${\overline{\omega }}_k=2\pi f_k/f_s$; therefore, the estimated frequencies are

$${\widehat{f}}_k=f_s(\angle {\widehat{\varphi }}_k)/\left(2\pi \right)\mathrm{for}k=1,2,\cdots ,Y$$

(34)

5. Fast Detection of System Operating Conditions

The input-signal parameters (magnitude, phase angle, frequency) may abruptly change/oscillate due to power system events, such as faults. Termed abnormal conditions in this paper, these sudden changes can quickly be detected using ${\widehat{\mathbf{R}}}_x$. Let the estimated autocorrelation matrix corresponding to the previous data string be denoted by ${\widehat{\mathbf{R}}}_{x_{old}}$. The new autocorrelation matrix ${\widehat{\mathbf{R}}}_{x_{new}}$ is estimated using (24) when the most recent data string is obtained. We can now create an $M\times M$ matrix ${\Delta \widehat{\mathbf{R}}}_x$ as shown below:

$${\Delta \widehat{\mathbf{R}}}_x\stackrel{\Delta }{=}{\widehat{\mathbf{R}}}_{x_{new}}-{\widehat{\mathbf{R}}}_{x_{old}}$$

(35)

The Frobenius norm of ${\Delta \widehat{\mathbf{R}}}_x$, $\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F$, is given by

$$\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F=\sqrt{\sum _{i=1}^M\sum _{j=1}^M|{\Delta \widehat{\mathbf{R}}}_x{(i,j)|}^2}$$

(36)

A threshold $\zeta$ can now be defined for $\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F$ to determine the system’s current operating condition as follows:

(1)

If $\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F\le \zeta$, it implies that the system operates under normal conditions and the signal parameters have varied slightly. In such conditions, the estimated phasors from previous and current data strings can be used for accurate estimation of frequency (see Section 6.1).

(2)

If $\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F>\zeta$, it implies that the system is faced with an abnormal condition resulting in the signal parameters differing significantly from the previous data string. In such conditions, the proposed iESPRIT must be used for precise estimation of frequency, and subsequent correction of the estimated phasor (see Section 6.2).

The impact of $\zeta$ on the accuracy and computational burden of the proposed approach is investigated in detail in Section 8.6.

6. Accuracy Enhancement of Estimated Phasor and Frequency under Detected Conditions

The main goal of this section is to improve the accuracy of the estimated phasor (Section 3) and frequency (Section 4) based on the system operating condition (normal/abnormal) identified in Section 5.

6.1. Implementation of the Phasor Estimation Algorithm for Frequency Estimation under Normal Conditions

Let the estimated fundamental phasor obtained from the previous N-length data string with sample numbers $\{r,r+1,\cdots ,r+N-1\}$ be denoted by ${\widehat{X}}_{1old}$, which is calculated as follows at ${\omega }_{1old}=2\pi f_{1old}$.

$${\widehat{X}}_{1old}=\frac{2}{N}\sum _{n=r}^{r+N-1}x\left[n\right]e^{-jn{\omega }_{1old}\Delta t}$$

(37)

Similarly, ${\widehat{X}}_{1new}$ is obtained from the current N-length data string with sample numbers $\{r+1,r+2,\cdots ,r+N\}$ at ${\omega }_{1new}=2\pi f_{1new}$ as

$${\widehat{X}}_{1new}=\frac{2}{N}\sum _{n=r+1}^{r+N}x\left[n\right]e^{-jn{\omega }_{1new}\Delta t}$$

(38)

It is proved in Appendix C that, under normal conditions, the parameter $\acute{P}$ (defined in (12)) can be obtained as follows.

$$\acute{P}\approxeq {\widehat{X}}_{1old}/{\widehat{X}}_{1new}$$

(39)

Now, as $\acute{P}$ only depends on N, r, and $\Delta {\omega }_1$, for a pre-determined value of N, $\acute{P}$ can be calculated offline for different values of r and $\Delta {\omega }_1$ in the form of a look-up table. Then, in real-time, the $\Delta {\widehat{\omega }}_1$ can be obtained with high accuracy from its corresponding $\acute{P}$ using the look-up table. Finally, the estimated frequency can be updated during normal conditions as

$${\widehat{f}}_{1new}={\widehat{f}}_{1old}+\Delta {\widehat{\omega }}_1/\left(2\pi \right)$$

(40)

The steps of the proposed frequency estimation algorithm are summarized in Algorithm 1.

Algorithm 1: Proposed Frequency Estimation Algorithm

Input: signal samples, ${\widehat{X}}_{1old}$, ${\widehat{X}}_{1new}$  Output: estimated frequency   Initialization: M and L 1: take L-length of input data 2: create ${\mathbf{X}}_f$, estimate ${\widehat{\mathbf{R}}}_{x_{new}}$, and calculate ${\Delta \widehat{\mathbf{R}}}_x$ using (22), (24) and (35), respectively 3: calculate $\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F$ using (36) 4: save ${\widehat{\mathbf{R}}}_{x_{new}}$ as ${\widehat{\mathbf{R}}}_{x_{old}}$ for the next estimation 5: if ($\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F\le \zeta$) then 6:  calculate $\acute{P}$ using (39) 7:  find $\Delta {\widehat{\omega }}_1$ using the lookup table 8:  estimate frequency using (40) 9:  go to line 18 10: end if 11: if ($\left|\right|{\Delta \widehat{\mathbf{R}}}_x{\left|\right|}_F>\zeta$) then 12:   calculate eigenvectors of ${\widehat{\mathbf{R}}}_x$ 13:   find ${\mathbf{Q}}_S$ by (26), split it into ${\mathbf{Q}}_1$ and ${\mathbf{Q}}_2$ by (28) 14:   estimate $\widehat{\Psi }$ using (31) 15:   calculate eigenvalues of $\widehat{\Psi }$ 16:   estimate frequency using (34) 17: end if 18: return ${\widehat{f}}_1$

6.2. Implementation of the Proposed Frequency Estimation Algorithm for Phasor Estimation under Abnormal Conditions

The first valid estimated frequency realized using Algorithm 1 is obtained after receiving the first L-length of input data; let this frequency be denoted by ${\widehat{f}}_{1,L}$. It will be shown in Section 7.4 that the proposed iESPRIT is able to estimate the frequency precisely even when L is lower than one-cycle-length of input data (i.e., $L<N$), which is a significant achievement especially for protection applications. Therefore, it is possible to correct the estimated phasors by using $\acute{P}$ which is dependent on the accurate frequency ${\widehat{\omega }}_1$ (see (9), (10) and (12)). To do so, we calculate the median ($Mdn$) of a set $\mathfrak{F}$ that contains the previous and current estimated frequencies (i.e., $\mathfrak{F}=\{{\widehat{f}}_{1,L},{\widehat{f}}_{1,L+1},\cdots ,{\widehat{f}}_{1,N}\}$) as shown in (41).

$${\widehat{\omega }}_1=2\pi .Mdn\left(\mathfrak{F}\right)$$

(41)

Employing the median brings additional robustness against noise. The value of the actual angular frequency ${\widehat{\omega }}_1$, obtained from (41), can be used to correct the estimated phasors during abnormal conditions. Moreover, ${\widehat{\omega }}_1$ can also be used to improve the performance of the proposed filter design (see Section 2) by replacing ${\omega }_{nom}$ with ${\widehat{\omega }}_1$ in (5) and (6). In this regard, the filter’s parameters can be calculated offline, for the pre-chosen filter order (see Section 7.1), for different values of ${\omega }_1$ using (5) and (6), and stored in the form of a look-up table. In real time, the filter’s parameters can be adjusted quickly using ${\widehat{\omega }}_1$ and the look-up table.

Finally, note that the proposed frequency estimator itself provides strong immunity against noise and can be used for various power system applications that require accurate and fast frequency measurements, such as under-frequency relaying, power system stabilization/restoration, protection against loss of synchronism, and rapid synthetic inertia control.

Flowchart of the proposed approach for both phasor and frequency estimation is shown in Figure 3 in which normal and abnormal system operating conditions have been distinguished with green dashed and red solid arrows, respectively.

7. Parameter Selection for Implementation

7.1. Choosing Optimal Value for the Filter Order R

To determine the filter order R of the complex gain filter (see Section 2), 100 synthetic test signals containing both fundamental and DDC components are generated. In the test signals, the $D_0$ and $\tau$ (DDC parameters) are randomly selected from a range of $\left[0,1p.u.\right]$ and $\left[0.5,5\right]\left(\mathrm{cycles}\right)$, respectively [5]. Seven different values $\{3,4,6,8,10,15,20\}$ are considered for R. The input signals are passed through the proposed filter and DFT is then performed for estimation of the fundamental phasor. For each filter order, the maximum values of TVE [29,30] and settling time $t_S$ [4], obtained over the 100 test signals are given in

Table 1. Maximum values of TVE and $t_S$ over 100 input test signals for different values of R.

R	3	4	6	8	10	15	20
Max. TVE (%)	0.76	0.67	0.28	0.22	0.14	0.04	0.01
Max.$t_S$ (ms)	20.7	21.7	26.9	30.8	32.9	36.2	38.2

.

As can be seen from Table 1, by increasing R, the DDC is suppressed more effectively resulting in reduction of the maximum TVE at the cost of increasing the maximum $t_S$; i.e., there is a trade-off in choosing the optimal value of the filter order. Table 1 indicates that an acceptable value for the order of the proposed filter would be six since this provides a reasonable balance between DDC rejection (significantly lower maximum TVE than 1%) and response time.

7.2. Choosing Optimal Value for N

The impact of N (or equivalently $f_s$) on the accuracy of the estimated phasor by the proposed algorithm is investigated in this section. In this regard, a synthetic test signal containing a fundamental sinusoidal component that is distorted by DDC, harmonics, and noise is used as an input to the proposed algorithm. Six different values $\{16,32,128,512,1024,16,384\}$ are considered for N. For each value, the fundamental phasor is estimated by the proposed algorithm, and the maximum TVE is calculated, which is given in

Table 2. Impact of N on accuracy of the estimated phasor by the proposed algorithm.

N	16	32	128	512	1024	16,384
$\mathbf{Max}.\mathbf{TVE}$ (%)	0.47	0.43	0.40	0.39	0.39	0.38

.

As can be seen in Table 2, the phasor estimation accuracy of the proposed algorithm converges to its limit value and is negligibly affected by increasing the number of samples per cycle, N (or equivalently the sampling frequency, $f_s$). On the other hand, increasing the sampling frequency causes (i) difficulties in hardware implementation of the phasor and frequency estimation algorithms in numerical relays, and (ii) increased computational burden, which is a major concern for the protection application. Based on Table 2, an appropriate value of N was found to be 32.

7.3. Choosing Optimal Value for M

The value of M plays a key role in computational complexity and robustness of the proposed iESPRIT against noise. As such, seven different values $\{3,4,5,6,8,10,15\}$ are considered for M. A synthetic test signal containing a sinusoidal component (with known frequency) which is distorted by noise is used as an input to the proposed iESPRIT. Note that the value of L (length of data string in (22)) is kept constant at N. The maximum absolute frequency estimation error (denoted by Max. |FE|) [30] and the execution time taken for estimation are presented in

Table 3. Maximum value of |FE| and execution time for different values of M.

M	3	4	5	6	8	10	15
Max. |FE| (mHz)	0.86	0.41	0.34	0.18	0.09	0.03	0.003
Exe. Time (ms)	0.12	0.23	0.47	0.81	1.74	3.05	9.7

for each of the seven values of M. Note that all of the simulations done in this paper have been carried out on a system with Core 2 Duo CPU (2.67 GHz) processor and 4 GB RAM.

It can be concluded from Table 3 that more robustness against noise and higher precision are achieved by choosing greater values for M, at the cost of increasing the execution time due to increasing complexity of calculations. Therefore, there is a trade-off between the precision of estimation and computational burden. From Table 3, an appropriate value of M is found to be four.

7.4. Choosing Optimal Value for L

The value for L (length of data string) plays a key role in convergence speed of the proposed iESPRIT. In this regard, six different values $\{N/4,N/2,3N/4,N,3N/2,2N\}$ are considered for L. Likewise, as we did in Section 7.3, a synthetic signal is used while M is kept constant at four.

The Max |FE| and settling time $t_S$ are given in

Table 4. Maximum value of |FE| and $t_S$ for different values of L.

L	$N/4$	$N/2$	$3N/4$	N	$3N/2$	2N
Max. |FE| (mHz)	1.03	0.33	0.25	0.04	0.007	0.0008
$t_S$ (ms)	6.8	11.4	14.3	17.2	29.1	33.9

for the six values of L. Considering optimal accuracy and convergence speed of the proposed iESPRIT, an appropriate value of L is found to be $N/2$.

Note that the empirically obtained values of the parameters identified in Section 7 are kept constant for all the simulations performed in Section 8.

8. Simulation Results

In the following sub-sections, several static and dynamic test signals that comply with the IEEE and IEC Standards [29,30], as well as power system signals, are used to investigate the performance of the proposed approach. The above-mentioned standards only consider changes in the input signal parameters (i.e., magnitude, phase angle, and frequency). However, in this paper, the performance of the proposed algorithm is also evaluated in the presence of stray components, such as DDC and noise. Four recently published phasor and/or frequency estimation methods, applicable to P-class PMUs, have been selected for comparison purposes: (i) A phasor estimation method presented in [5], which is based on DFT and down sampling (called “DS-DFT”), (ii) a phasor and frequency estimation method presented in [13], which is based on LS and amplitude modulation modeling of the input signal (called “AM”), (iii) a phasor and frequency estimation method presented in [13], which is based on LS and phase modulation modeling of the input signal (called “PM”), and finally (iv) a subspace-based frequency estimation method inspired by the multiple signal classification (MUSIC) technique presented in [24] (called “${\mathrm{A}}^2\mathrm{MUSIC}$”). The nominal frequency $f_{nom}$, sampling frequency $f_s$, and $\zeta$ are set to 60 Hz, 1920 Hz, and ${10}^{-4}$, respectively.

8.1. Performance Evaluation in the Presence of DC Offset, DDC, Harmonics, and Noise in Both Nominal and Off-Nominal Frequency Conditions

The following test signal is used to investigate the performance of the proposed approach for phasor and frequency estimation in the presence of DC offset, DDC, harmonics, and noise, in both nominal and off-nominal frequency conditions

$$x_1\left(t\right)=A_0+D_0{e}^{-t/\tau }+\sum _{k=1}^K\frac{A_1}{k}cos(2\pi k{f}_{1}t+k{\theta }_1)+w\left(t\right)$$

(42)

where the total number of sinusoids K is set equal to 50 [29,30]. Five-hundred test signals are generated by selecting different values for $A_0$, $f_1$ (varied from 58 to 62 Hz [29]), $D_0$ (varied from 0 to $1p.u.$), $\tau$ (varied from 0.5 to 5 cycles), and signal-to-noise ratio (SNR) (varied from 10 to 80 dB) [24]. The fundamental phasor and frequency are estimated by the proposed approach and other understudied algorithms, and the mean, root mean square (RMS), and maximum values of maximum of TVE and |FE| are calculated. From the results shown in

Table 5. Mean, RMS, and maximum values of the indices over the 500 test signals in the presence of undesirable components.

Algorithm	Max. TVE (%)			Max. |FE| (mHz)
	Mean	RMS	Max.	Mean	RMS	Max.
DS-DFT [5]	1.05	2.09	4.86	N/A	N/A	N/A
AM [13]	1.23	2.54	5.78	8.23	11.4	32.4
PM [13]	1.06	2.38	4.97	6.51	9.21	34.6
${\mathrm{A}}^2\mathrm{MUSIC}$ [24]	N/A	N/A	N/A	2.28	3.74	6.63
Proposed Approach	0.58	0.67	0.79	0.56	0.69	1.01

, it is clear that the proposed approach gives higher accuracy for both phasor (Max. TVE) and frequency (Max. |FE|) estimation in comparison to the other methods.

8.2. Performance Evaluation for Dynamic Frequency Ramp Test

In this test, for the same test signals generated in Section 8.1, $f_1$ is ramped to $f_1+5$ (Hz) from 0 to 5 (s) at a rate of 1 Hz/s [3,29,30]. The mean, RMS, and maximum values of the two indices, namely, Max. TVE and Max. |FE|, for the 500 dynamic test signals are given in

Table 6. Mean, RMS, and maximum values of the indices over the 500 dynamic test signals during the frequency ramp.

Algorithm	Max. TVE (%)			Max. |FE| (mHz)
	Mean	RMS	Max.	Mean	RMS	Max.
DS-DFT [5]	1.32	2.59	5.66	N/A	N/A	N/A
AM [13]	1.35	2.67	5.93	9.14	12.6	36.21
PM [13]	1.13	2.52	5.07	6.63	9.42	35.06
${\mathrm{A}}^2\mathrm{MUSIC}$ [24]	N/A	N/A	N/A	2.69	4.17	7.12
Proposed Approach	0.67	0.78	0.89	0.66	0.80	1.07

.

As seen in Table 6, the proposed approach gives better performance for both phasor and frequency estimation under frequency ramps in comparison to the other methods. Apart from the high estimation accuracy, the proposed approach also has faster response time under dynamic conditions which will be explored in more details in Section 8.4.

8.3. Modulation Test

Transmission line faults/outages and bulk load switching often cause rotor-angle oscillations which may lead to power swings and subsequently mal-operation of the protective relays [32]. Generally, power swings are mathematically modeled as a modulated sinusoidal signal [13,29,30,32]

$$x_2\left(t\right)=X_m(1+k_{x}cos\left({\omega }_{m}t\right))\times cos({\omega }_{nom}t+k_{a}cos({\omega }_{m}t-\pi ))$$

(43)

where $k_x$ and $k_a$ denote amplitude and phase modulation coefficients, respectively, both of which are set to 0.1 [29,30]. In this test, the frequency of the input signal remains constant at 60 Hz, but the magnitude and phase-angle oscillate with an angular frequency of ${\omega }_m=2\pi f_m$ where $f_m$ changes from 0.1 Hz to 2 Hz in steps of 0.1 Hz as recommended by the Standards [29,30] for P-class PMUs. The maximum values of TVE and |FE| are given in

Table 7. Maximum TVE and |FE| for the modulation test.

Algorithm	Max. TVE (%)	Max. |FE| (mHz)
DS-DFT [5]	2.37	N/A
AM [13]	0.0096	2.3
PM [13]	0.0008	1.9
${\mathrm{A}}^2\mathrm{MUSIC}$ [24]	N/A	0.95
Proposed Approach	0.43	0.64

for different algorithms.

As can be observed from Table 7, the proposed approach gives better frequency estimation accuracy compared to the other algorithms. Since the algorithms presented in [13] are based on amplitude modulation (AM) and phase modulation (PM) modeling of the input signals, they have the highest phasor estimation accuracy (their Max. TVE is lowest). However, the proposed phasor estimation approach also satisfies the modulation test requirements for P-class PMUs specified in the Standards [29,30].

8.4. Step Change Test

Abrupt changes in the impedance and network configuration during faults may cause current/voltage transients which are mathematically expressed as [13,29,30,32]

$$x_3\left(t\right)=X_m(1+k_{x}u\left(t\right))\times cos({\omega }_{nom}t+k_{a}u\left(t\right))$$

(44)

where $k_x$ and $k_a$ are the magnitudes of the unit step function $u\left(t\right)$ and set to 0.1 and $\pi /18$ rad, respectively, as recommended in the Standards [29,30]. Two well-known indices, namely, response time and maximum over/undershoot, which are defined in [30], are used to evaluate the performance of the different algorithms from the viewpoints of speed and accuracy. The results are given in

Table 8. Performance indices for the step change test.

Algorithm	$k_x$	$k_a$	Phasor		Frequency
			Rsp. (ms)	Max O/U.S. (%) *	Rsp. (ms)
DS-DFT [5]	0.1	0	14.6	4.1	N/A
	0	$\pi /18$	15.5	4.6	N/A
AM [13]	0.1	0	15.3	5	41.3
	0	$\pi /18$	28.6	4	40.7
PM [13]	0.1	0	14.2	5	57.7
	0	$\pi /18$	30.9	4	57.1
${\mathrm{A}}^2\mathrm{MUSIC}$ [24]	0.1	0	N/A	N/A	12.4
	0	$\pi /18$	N/A	N/A	12.7
Proposed Approach	0.1	0	13.9	2.7	11.6
	0	$\pi /18$	14.4	3.3	12.1

* Max O/U.S. values are computed relative to the magnitude of the step [30].

in which the response time and maximum over/undershoot are denoted by Rsp. and Max O/U.S., respectively. Based on the results obtained in Table 8, it is clear that the proposed approach has superior performance in terms of speed and accuracy, and its performance indices can comply with the requirements for P-class PMUs mentioned in the Standards [29,30].

8.5. Performance Evaluation Using Power System Signals

Power system signals of the IEEE 9-bus test system [13] are used to evaluate the performance of the proposed approach for both phasor and frequency estimation when severe oscillations in magnitude, phase angle, and frequency of the input signals take place. Two different scenarios are considered [13]:

Scenario 1: Outage of line 7-8 at $t=0$ s. Voltage of bus 6 is used as the input signal.

Scenario 2:${\mathrm{G}}_2$ (connected to bus 2) trips at $t=50$ ms. Current of line 6-4 is used as the input signal.

The maximum values of TVE and |FE| are given in

Table 9. Maximum TVE (%) and |FE| (mHz) for each of the scenarios *.

Sc.	DS-DFT [5]	AM [13]	PM [13]	${\mathrm{A}}^2\mathrm{MUSIC}$ [24]	Proposed Approach
1	1.11 N/A	2.16 10.3	1.04 6.8	N/A 3.2	0.72 0.81
2	1.26 N/A	2.23 11.2	1.12 7.3	N/A 3.8	0.78 0.85

* For each cell, the first row is Max. TVE and the second row is Max. |FE|.

for the two scenarios. As can be observed from Table 9, the proposed approach provides higher accuracy compared to the other methods. It is also worth mentioning that the abnormal events in both the scenarios were detected by the proposed approach in less than one cycle.

8.6. Computational Burden and Sensitivity Analysis

The threshold $\zeta$ has a significant impact on the accuracy and the computational complexity of the proposed approach (see Section 5). Hence, a sensitivity analysis is performed to investigate the effect of $\zeta$ on the accuracy of the estimated phasors and frequencies and the computational burden of the proposed approach. In this regard, Max. TVE, Max. |FE|, and average execution time (taken for each estimation) are given in

Table 10. Performance indices for different values of $\zeta$.

Test	$\mathit{\zeta }$	Max. TVE (%)	Max. |FE| (mHz)	Execution Time ($\mathsf{\mu }$s)
Steady-State (Section 8.1)	${10}^{-2}$	0.97	5.02	46
	${10}^{-4}$	0.58	0.56	283
	${10}^{-6}$	0.51	0.47	405
Frequency Ramp (Section 8.2)	${10}^{-2}$	1.09	5.17	57
	${10}^{-4}$	0.67	0.66	367
	${10}^{-6}$	0.59	0.58	472
Modulation Test (Section 8.3)	${10}^{-2}$	1.03	5.11	52
	${10}^{-4}$	0.43	0.61	319
	${10}^{-6}$	0.36	0.56	423
Magnitude Step Change (Section 8.4)	${10}^{-2}$	1.12	5.09	54
	${10}^{-4}$	0.69	0.59	326
	${10}^{-6}$	0.60	0.50	446
Phase Angle Step Change (Section 8.4)	${10}^{-2}$	1.14	5.15	56
	${10}^{-4}$	0.71	0.63	339
	${10}^{-6}$	0.63	0.57	468
Scenario 1 (Section 8.5)	${10}^{-2}$	1.16	5.21	59
	${10}^{-4}$	0.72	0.81	370
	${10}^{-6}$	0.64	0.74	484
Scenario 2 (Section 8.5)	${10}^{-2}$	1.25	5.26	63
	${10}^{-4}$	0.78	0.85	391
	${10}^{-6}$	0.71	0.77	510

for different values of $\zeta$ for each of the test signals used in this section.

As can be observed from Table 10, higher accuracy in estimation of phasor and frequency is achieved by choosing smaller values for $\zeta$, at the cost of increasing the execution time (or equivalently computational burden/complexity) i.e., there is a trade-off in choosing the optimal value of the threshold $\zeta$. Table 10 indicates that an acceptable value for $\zeta$ would be ${10}^{-4}$ since this provides a reasonable balance between the estimation accuracy and computational burden. The mean execution time of the proposed approach (over all of the test signals in this section) for each synchrophasor and frequency estimation by the adoption of ${10}^{-4}$ for $\zeta$ was around 0.34214 (ms), which makes it suitable for real-time applications.

8.7. Noise Propagation through the Proposed Algorithm

In this section, the following test signal is used to investigate the propagation of noise through the proposed algorithm:

$$x_4\left[n\right]=A_{1}sin(\frac{2\pi n{f}_1}{f_s}+{\theta }_1)+w\left[n\right]\forall n=1,2,\cdots ,N$$

(45)

Cramér Rao bounds (CRBs) that give lower bounds on the variances of estimated signal parameters (i.e., magnitude, phase angle, and frequency) are obtained for the above-mentioned test signal (45) as follows: [36].

$$\begin{array}{cc}\hfill CRB\left(A_1\right)& =\frac{2{{\sigma }_n}^2}{N}\hfill \\ \hfill CRB\left({\theta }_1\right)& =\frac{4{{\sigma }_n}^2}{{A_1}^2}\frac{2N+1}{N(N-1)}\hfill \\ \hfill CRB\left(f_1\right)& =\frac{24{{\sigma }_n}^2}{{A_1}^{2}N(N^2-1)}{\left(\frac{f_s}{2\pi }\right)}^2\hfill \end{array}$$

(46)

In order to investigate the propagation of noise through the proposed algorithm, 100 test signals are generated by choosing different values for $f_1$ between 58 and 62 (Hz) [29,30]. The variance of the noise in the sinusoidal signal ${{\sigma }_n}^2$ is set to be $0.{01}^2$ [36]. For each of the generated test signals, the signal parameters are estimated by the proposed algorithm. The mean square error (MSE) [37] as well as the CRBs are given in

Table 11. MSEs and CRBs of the estimated parameters by the proposed algorithm.

Parameter	CRB	MSE
$A_1$	$6.25\times {10}^{-6}$	$1.35\times {10}^{-5}$
${\theta }_1$	$2.62096\times {10}^{-5}$	$1.46\times {10}^{-4}$
$f_1$	$6.8458\times {10}^{-3}$	$9.23\times {10}^{-3}$

. As can be seen from Table 11, the values of MSE for all of the estimated parameters are very close to their corresponding CRBs.

9. Conclusions

This paper proposes a fast, accurate, robust, and computationally-efficient approach for estimating phasor and frequency that can be applied to both voltage and current signals without making simplifying assumptions regarding the components present in the signal. The first step of the proposed methodology involves the removal of the DDC that may be present in the input signal using optimal filtering. Care is taken to ensure that the parameters of the fundamental component are left unaltered. In the next step, initial phasor and frequency estimates are obtained in parallel using DFT and the proposed iESPRIT, respectively. Finally, to further improve accuracy, the initial estimates of phasor and frequency are updated, based on the identified system condition.

The simulation results confirm that the proposed approach has high accuracy, fast response time, low computational complexity, and is robust against grid disturbances and off-nominal frequency conditions. Furthermore, the proposed approach fulfills the requirements recommended by the IEEE and IEC standards for P-class PMUs in both steady-state and dynamic conditions. This makes the proposed approach suitable for a variety of power system monitoring, protection, and control applications that need fast and accurate estimates of voltage and current phasors and system frequency.

To provide more robustness against stray components such as harmonics and inter-harmonics, particularly for measurement applications, the combination of more recent versions of DFT such as interpolated DFT (IpDFT) with iESPRIT will be explored in the future.