IpDFT-Tuned Estimation Algorithms for PMUs: Overview and Performance Comparison

Featured Application

Real-time measurement of amplitude, phase, frequency and rate of change of frequency (ROCOF) of current and voltage AC signals in power transmission and distribution systems.

Abstract

The Interpolated Discrete Fourier Transform (IpDFT) is one of the most popular algorithms for Phasor Measurement Units (PMUs), due to its quite low computational complexity and its good accuracy in various operating conditions. However, the basic IpDFT algorithm can be used also as a preliminary estimator of the amplitude, phase, frequency and rate of change of frequency of voltage or current AC waveforms at times synchronized to the Universal Coordinated Time (UTC). Indeed, another cascaded algorithm can be used to refine the waveform parameters estimation. In this context, the main novelty of this work is a fair and extensive performance comparison of three different state-of-the-art IpDFT-tuned estimation algorithms for PMUs. The three algorithms are: (i) the so-called corrected IpDFT (IpDFTc), which is conceived to compensate for the effect of both the image of the fundamental tone and second-order harmonic; (ii) a frequency-tuned version of the Taylor Weighted Least-Squares (TWLS) algorithm, and (iii) the frequency Down-Conversion and low-pass Filtering (DCF) technique described also in the IEEE/IEC Standard 60255-118-1:2018. The simulation results obtained in the P Class and M Class testing conditions specified in the same Standard show that the IpDFTc algorithm is generally preferable under the effect of steady-state disturbances. On the contrary, the tuned TWLS estimator is usually the best solution when dynamic changes of amplitude, phase or frequency occur. In transient conditions (i.e., under the effect of amplitude or phase steps), the IpDFTc and the tuned TWLS algorithms do not clearly outperform one another. The DCF approach generally returns the worst results. However, its actual performances heavily depend on the adopted low-pass filter.

1. Introduction

The Phasor Measurement Units (PMUs) are instruments designed to measure amplitude, phase, frequency, and Rate Of Change Of Frequency (ROCOF) of voltage and current AC signals in crucial locations of transmission and distribution systems [1]. While up to a few years ago PMUs were mainly conceived to detect anomalous operating conditions in transmission systems [2], nowadays they are regarded as essential devices to monitor the actual operating conditions of smart distribution grids in real-time, e.g., for island maneuvering [3], faulted line detection and identification [4], short-term voltage stability monitoring [5], and state estimation [6,7,8,9].

The key features of PMUs compared with other measurement devices for power system applications are:

• the capability to measure all quantities at times synchronized with the Universal Coordinated Time (UTC);
• the capability to stream real-time data to the so-called Phasor Data Concentrators (PDCs) with reporting rates ranging between about 10 frame/s and 100 frame/s;
• high measurement accuracy, as specified in the IEEE/IEC Standard 60255-118-1:2018, in the following simply referred to as the “Standard” for the sake of brevity [10];
• fast responsiveness (in the order of a few power line cycles) in order to trigger adequate protection schemes when severe phasor magnitude or phase changes occur due to faults or grid topology changes [11]. In this respect, the IEEE/IEC Standard classifies PMUs as protection-oriented (P Class) or measurement-oriented (M Class) devices. The former ones can indeed be less accurate than the latter, but of course they must be considerably faster.

Quite importantly, in the case of active distribution systems, the maximum Total Vector Error (TVE) of PMU measurements has to be even smaller than the Standard limits, i.e., in the order of 0.1% [12], and with phase uncertainty contributions not greater than a few mrad [13].

The need to meet the specifications shortly summarized above required remarkable research efforts in the area of signal processing and metrology. In particular, both flexible testbeds for PMU algorithm characterization [14], and new PMU calibrators were developed [15,16].

Among the many estimation algorithms for PMUs proposed over the last few years, the Interpolated Discrete Fourier Transform (IpDFT) is definitely one of the most popular, due to its low computational complexity and its good capability to compensate for the effect of off-nominal frequency deviations [17]. Such deviations might otherwise heavily affect the waveform parameter estimation based on the classic Discrete Fourier Transform (DFT) [18]. Even if the use of the IpDFT in signal processing is well established (e.g., for multi-tone parameter estimation [19]), its use in power systems monitoring applications was proposed only a few years ago. Since then, a number of IpDFT variants and multi-step algorithms based on the IpDFT stage have been developed to overcome the inherent limitation of the basic algorithm. In particular, the classic IpDFT algorithm can be used to obtain a preliminary estimate of some quantities of interest (most notably the fundamental frequency of the analyzed signal) prior to applying a further algorithm to return more accurate estimates of synchrophasor amplitude, phase, frequency, and ROCOF.

For instance, in [20], a three-point IpDFT algorithm is proposed to compensate for the infiltration of the spectral image of the fundamental component on the estimated parameters. However, this algorithm returns only a frequency estimate. This is indeed the most common case [21,22,23], although a dynamic ROCOF estimator based on IpDFT and Kalman filtering has been recently presented in [24]. The same idea also underlies the so-called enhanced IpDFT [25], which mitigates the detrimental effect of the image component on the estimated synchrophasor parameters by applying the IpDFT iteratively. In [26] this approach is further extended to compensate for the effect of harmonics and inter-harmonic components. The so-called corrected IpDFT (IpDFTc) consists of a two-step IpDFT designed to remove the spectral leakage due to both the image of the fundamental tone and the second-order harmonic [27]. In fact, these steady-state components are often the most critical ones for synchrophasor estimation, especially over short observation intervals. Quite importantly, the IpDFTc algorithm is not iterative. Therefore, it requires a lower computational burden than the estimator proposed in [26].

In [28] and [29] the fundamental frequency estimated through a classic IpDFT is used to tune the matrix coefficients of real-valued weighted least-squares estimators obtained from the Taylor’s series expansion of the dynamic synchrophasor model. In the original version of such estimators, the matrix coefficients are computed at the nominal frequency. As a consequence, the estimates are significantly affected by off-nominal frequency deviations. The IpDFT-based frequency tuning instead significantly improves the estimation accuracy of both the so-called Taylor-Fourier filter (which does not include the harmonics in the signal model) [30], and the Taylor-Fourier Transform (that instead includes a number of harmonics in the signal model) [31].

The preliminary fundamental frequency estimation based on the IpDFT can be also used to perform a better direct down-conversion of the voltage or current waveforms, before extracting synchrophasor amplitude and phase through low-pass filtering [32].

To the best of the Authors’ knowledge, the three IpDFT-tuned algorithms described above are among those that exhibit excellent performance in terms of both accuracy and computational complexity. Unfortunately, the results reported in the scientific literature do not allow a fair and easy comparison of their actual performances because different observation interval lengths or testing conditions are considered. Thus, this paper aims at presenting a thorough performance comparison of the IpDFTc estimator [27], the tuned Taylor Weighted Least Squares (TWLS) algorithm [28], and a tuned implementation of the frequency down-conversion and low-pass filtering (DCF) technique [32].

The rest of the paper is structured as follows. In Section 2, the theoretical background of the algorithms considered is briefly presented. In Section 3, the simulation results achieved in the steady-state, dynamic and transient testing conditions described in the IEEE/IEC Standard 60255-118-1:2018 are shown and commented on. Section 4 presents some final comments on the considered IpDFT-tuned algorithm compliance, advantages and disadvantages. Section 5 presents some further experimental results in borderline operating conditions compliant with the Standard EN 50160:2010 [33]. Finally, Section 6 concludes the paper.

2. Theoretical Background

In this Section, firstly a brief overview of the signal model used for synchrophasor, frequency, and ROCOF estimation is described. Then, the basic IpDFT algorithm used as a preliminary stage of all the techniques under test is recalled. Finally, the theoretical background of the three IpDFT-tuned estimators introduced in Section 1 is presented in a common framework.

2.1. Signal Model and Basic IpDFT Algorithm

Since the current or voltage AC waveform parameters change with time, a general AC waveform model is:

$$y\left(t\right)=s\left(t\right)+\eta \left(t\right)=a_1\left(t\right)\cos \left(2\pi f_1\left(t\right)t+{\varphi }_1\left(t\right)\right)+\eta \left(t\right)$$

(1)

where $s\left(t\right)=a_1\left(t\right)\cos \left(2\pi f_{1}t+{\varphi }_1\left(t\right)\right)$ is the fundamental sinewave (with amplitude, phase, and frequency at a given time t represented by a₁(t), $\varphi$₁(t), and f₁(t), respectively) and η(t) includes all possible disturbances, such as harmonics, inter-harmonics, and wideband noise. The sine-wave dynamic phasor is defined as $p\left(t\right)\stackrel{\mathrm{def}}{=}{a}_1\left(t\right)e^{j{\varphi }_1\left(t\right)}$. It is called a synchrophasor when it is referred to timestamped UTC instants. If waveform (1) is sampled by a PMU at a rate f_s and an odd number of samples M = 2N_h + 1 is acquired, the reference time t_r for synchrophasor estimation (with r being integer) can be chosen exactly in the middle of the r-th digitized data record, which is given by:

$$y_r\left(m\right)=a_1\left(rL+m\right)\cos \left(2\pi \frac{f_1\left(rL+m\right)}{f_s}\left(rL+m\right)+{\varphi }_1\left(rL+m\right)\right)+\eta \left(rL+m\right)$$

(2)

for $m=-N_h,\dots ,0,\dots ,N_h$, where L is the distance (in samples) between subsequent data records. Note that, in order to ensure a continuous data acquisition, L can range from 1 (if the observation intervals shift by 1 sample at a time) to M (when consecutive disjoint observation intervals are considered).

Thus, the objective of the estimation algorithms for PMUs is to determine the values of parameters a₁(t_r), ϕ₁(t_r) and f₁(t_r) and ${\frac{d{f}_1\left(t\right)}{dt}|}_{t_r}$ at the chosen UTC reference time within the data record (2). Such data have to be transferred to the PDC at a given Reporting Rate RR. For this reason, in practice it is often reasonable to choose $L=\frac{f_s}{RR}$. In the following, the waveform parameters of interest within the r-th data record will be shortly denoted as a_1,r, ϕ_1,r f_1,r and ROCOF_1,r, respectively.

The ratio between the fundamental AC frequency f_1,r and the sampling rate f_s can be expressed as:

$$\frac{f_{1,r}}{f_s}=\frac{{\nu }_{1,r}}{M}=\frac{J_1+{\delta }_{1,r}}{M},$$

(3)

where ν_1,r = J₁ + δ_1,r is the sinewave frequency expressed in bins, i.e., the number of acquired waveform cycles. Expression (3) results from the sum of an integer part J₁ and a fractional part δ_1,r (with −0.5 ≤ δ_1,r < 0.5). While J₁ typically does not change from record to record, the value of δ_1,r changes as a function of time due to frequency fluctuations, and becomes zero only if coherent sampling is performed. The sampling rate is selected so that the following condition is fulfilled:

$$\frac{f_n}{f_s}=\frac{J_1}{M-1}.$$

(4)

where f_n is the nominal sinewave frequency (i.e., 50 Hz or 60 Hz). Notice that the value of J₁ coincides also with the number of cycles acquired at the nominal frequency.

In practice, a sinewave is usually sampled non-coherently (i.e., δ_1,r ≠ 0). Therefore, the discrete spectrum of the acquired signal is affected by the spectral leakage due to the finite duration of the considered observation interval. The Maximum Sidelobe Decay (MSD) windows are often employed to reduce the impact of leakage on the estimated parameters [34,35]. The expression of a generic MSD window is:

$$w\left(m\right)=\frac{1}{2^{2H-2}}\left[{\mathrm{C}}_{2H-2}^{H-1}\cos \left(\frac{2\pi }{M}m\right)+2{{\displaystyle \sum }}_{h=1}^{H-1}{\mathrm{C}}_{2H-2}^{H-h-1}\cos \left(\frac{2\pi h}{M}m\right)\right], m=-N_h,\dots ,0,\dots ,N_h$$

(5)

where H ≥ 2 is the number of window terms and ${\mathrm{C}}_n^q=\frac{n !}{\left(n-q\right) ! q !}$. In particular, the 2-term MSD (or Hann) window will be used in this paper. Since the DFT of the r-th windowed data record is given by

$$Y_{r,w}\left(k\right)\stackrel{\mathrm{def}}{=}{{\displaystyle \sum }}_{m=-N_h}^{N_h}{y}_r\left(m\right)·w\left(m\right)e^{-j\frac{2\pi }{M}km},$$

(6)

the preliminary estimates of the sine-wave parameters ν₁, a₁, and ϕ₁ returned by the IpDFT algorithm are based on the DFT sample Y_r,w(J₁) and the larger of its two adjacent spectral samples, i.e., Y_r,w(J₁ − 1) or Y_r,w(J₁ + 1). In particular, it can be shown that [17]:

$${\widehat{\nu }}_{10,r}=J_1+{\widehat{\delta }}_{10,r}$$

(7)

with

$${\widehat{\delta }}_{10,r}=\frac{\left(H-1+i\right){\alpha }_{1,r}-H+i}{{\alpha }_{1,r}+1}$$

(8)

in which ${\alpha }_{1,r}=\frac{\left|Y_{r,w}\left(J_1+i\right)\right|}{\left|Y_{r,w}\left(J_1+i-1\right)\right|}$ with i = 0 if |Y_r,w(J₁ − 1)| > |Y_r,w(J₁ + 1)| or i = 1 if |Y_r,w(J₁ − 1)| < |Y_r,w(J₁ + 1)| and

$${\widehat{a}}_{10,r}=\frac{2\left|Y_{r,w}\left(J_1\right)\right|}{W\left({\widehat{\delta }}_{10,r}\right)},{\widehat{\varphi }}_{10,r}=angle\left\{Y_{r,w}\left(J_1\right)\right\}.$$

(9)

Note that in the estimator ${\widehat{a}}_{10,r}$, the DTFT of the window sequence (5) for M >> 1 and |λ| << M, is given with high accuracy by [36]

$$W\left(\lambda \right)\stackrel{\mathrm{def}}{=}{\displaystyle \sum }_{m=-N_h}^{N_h}w\left(m\right)e^{-j2\pi \frac{\lambda }{M}m}=\frac{M\sin \left(\pi \lambda \right)}{2^{2H-2}\pi \lambda }\frac{\left(2H-2\right)!}{{{\displaystyle \prod }}_{h=1}^{H-1}\left(h^2-{\lambda }^2\right)}.$$

(10)

which is a real-valued function due to the window even symmetry.

2.2. IpDFTc Algorithm

The IpDFTc estimator improves the performance of the basic IpDFT algorithm as it greatly reduces the effect on the estimated parameters of the spectral leakage due to both the image of the fundamental component and the second-order harmonic [27]. In the IpDFTc algorithm, the estimator of the number of observed signal cycles at the reference time within the r-th collected record is:

$${\widehat{\nu }}_{1,r}={\widehat{\nu }}_{10,r}-{\Delta }_c{\nu }_{1,r},$$

(11)

where ${\widehat{\nu }}_{10,r}$ is returned by (8) and the term ${\Delta }_c{\nu }_{1,r}$ is given by [27]:

$$\begin{array}{cc}\hfill {\Delta }_c{\nu }_{1,r}=& {\Delta }_c{\delta }_{1,r}=(H+{(-1)}^i{\widehat{\delta }}_{10,r})[\frac{{(-1)}^{i}2{\widehat{\nu }}_{10,r}}{2{\widehat{\nu }}_{10,r}-{\widehat{\delta }}_{10,r}+{(-1)}^{i+1}H}\frac{W(2{\widehat{\nu }}_{10,r}-{\widehat{\delta }}_{10,r})}{W\left({\widehat{\delta }}_{10,r}\right)}\cos (2{\widehat{\varphi }}_{10,r})+\hfill \\ \hfill & \frac{{(-1)}^i({\widehat{\nu }}_{20,r}-{\widehat{\nu }}_{10,r})}{{\widehat{\nu }}_{20,r}-{\widehat{\nu }}_{10,r}+{\widehat{\delta }}_{10,r}+{(-1)}^{i}H}\frac{{\widehat{a}}_{20,r}}{{\widehat{a}}_{10,r}}\frac{W({\widehat{\nu }}_{20,r}-{\widehat{\nu }}_{10,r}+{\widehat{\delta }}_{10,r})}{W\left({\widehat{\delta }}_{10,r}\right)}\cos ({\widehat{\varphi }}_{20,r}-{\widehat{\varphi }}_{10,r})].\hfill \end{array}$$

(12)

In (12) the values of ${\widehat{a}}_{10,r},$ ${\widehat{\varphi }}_{10,r}$ result from (9), namely from the preliminary IpDFT, whereas ${\widehat{a}}_{20,r}$ and ${\widehat{\varphi }}_{20,r}$ are the amplitude and the phase estimators of the possible second-order harmonic at the same reference time t_r. Let ${\tilde{y}}_r(·)$ be the residual signal obtained from (2) after subtracting the fundamental sinewave reconstructed by using the IpDFT estimated parameters, i.e.,

$${\tilde{y}}_r\left(m\right)=y_r\left(m\right)-{\widehat{a}}_{10,r}\cos \left(2\pi \frac{{\widehat{\nu }}_{10,r}}{M}\left(rL+m\right)+{\widehat{\varphi }}_{10,r}\right),m=-N_h,\dots ,0,\dots ,N_h.$$

(13)

If ${\tilde{Y}}_{r,w}(·)$ is the windowed DFT of (13) obtained using the H-term MSD window [i.e., ${\tilde{y}}_{r,w}\left(m\right)={\tilde{y}}_r\left(m\right)·w\left(m\right)$], the amplitude and phase estimators of the second-order harmonic are obtained by applying the IpDFT algorithm to (13), that is:

$${\widehat{a}}_{20,r}=\frac{2\left|{\tilde{Y}}_{r,w}\left(J_2\right)\right|}{W\left({\widehat{\delta }}_{20,r}\right)},{\widehat{\varphi }}_{20,r}=angle\left\{{\tilde{Y}}_{r,w}\left(J_2\right)\right\}.$$

(14)

In (14), J₂ denotes the integer part of the number of second harmonic cycles, while ${\widehat{\delta }}_{20,r}=\frac{\left(H-1+i\right){\alpha }_{2,r}-H+i}{{\alpha }_{2,r}+1}$ is the corresponding off-nominal frequency deviation estimator expressed in bins. Note that in ${\alpha }_{2,r}=\frac{\left|{\tilde{Y}}_{r,w}\left(J_2+i\right)\right|}{\left|{\tilde{Y}}_{r,w}\left(J_2+i-1\right)\right|}$ the index i is equal to 0 if $\left|{\tilde{Y}}_{r,w}\left(J_2-1\right)\right|>\left|{\tilde{Y}}_{r,w}\left(J_2+1\right)\right|$ or it is equal to 1 if $\left|{\tilde{Y}}_{r,w}\left(J_2-1\right)\right|<\left|{\tilde{Y}}_{r,w}\left(J_2+1\right)\right|$. Therefore, the second harmonic frequency, expressed in bins, is given by ${\widehat{\nu }}_{20,r}=J_2+{\widehat{\delta }}_{20,r}$. In order to compensate for the effect of the spectral leakage caused by both the fundamental image component and the second harmonic, the following expression is applied to estimate the sine-wave amplitude at time t_r:

$${\widehat{a}}_{1,r}={\widehat{a}}_{10c,r}-{\Delta }_c{a}_{1,r}$$

(15)

where

• ${\widehat{a}}_{10c,r}=\frac{2\left|Y_{r,w}\left(J_1\right)\right|}{W\left({\widehat{\delta }}_{10,r}-{\Delta }_c{\delta }_{1,r}\right)}$ (with ${\widehat{\delta }}_{10,r}$ and ${\Delta }_c{\delta }_{1,r}$ given by (8) and (12), respectively) is the amplitude estimator of the fundamental sinewave after compensating for the fundamental frequency estimation bias due to the second harmonic and
• ${\Delta }_c{a}_{1,r}={\widehat{a}}_{10,r}\frac{W\left(2{\widehat{\nu }}_{10,r}-{\widehat{\delta }}_{10,r}\right)}{W\left({\widehat{\delta }}_{10,r}\right)}\cos \left(2{\widehat{\varphi }}_{10,r}\right)+{\widehat{a}}_{20,r}\frac{W\left({\widehat{\nu }}_{20,r}-{\widehat{\nu }}_{10,r}+{\widehat{\delta }}_{10,r}\right)}{W\left({\widehat{\delta }}_{10,r}\right)}\cos \left({\widehat{\varphi }}_{20,r}-{\widehat{\varphi }}_{10,r}\right)$ is the amplitude correction term [27].

Similarly, the IpDTFc synchrophasor phase estimator at time t_r is

$${\widehat{\varphi }}_{1,r}={\widehat{\varphi }}_{10,r}-{\Delta }_c{\varphi }_{1,r},$$

(16)

where ${\Delta }_c{\varphi }_{1,r}=-\frac{W\left(2{\widehat{\nu }}_{10,r}-{\widehat{\delta }}_{10,r}\right)}{W\left({\widehat{\delta }}_{10,r}\right)}\sin \left(2{\widehat{\varphi }}_{10,r}\right)+\frac{{\widehat{a}}_{20,r}}{{\widehat{a}}_{10,r}}\frac{W\left({\widehat{\nu }}_{20,r}-{\widehat{\nu }}_{10,r}+{\widehat{\delta }}_{10,r}\right)}{W\left({\widehat{\delta }}_{10,r}\right)}\sin \left({\widehat{\varphi }}_{20,r}-{\widehat{\varphi }}_{10,r}\right)$ and ${\widehat{\varphi }}_{10,r}$ results from (10).

Finally, the fundamental frequencies and the ROCOF are estimated respectively as:

$${\widehat{f}}_{1,r}=f_n+\frac{f_s}{M}{\widehat{\delta }}_{1,r}.$$

(17)

and

$$\widehat{ROCO{F}_{1,r}}=\left|{\widehat{f}}_{1,r}-{\widehat{f}}_{1,r-1}\right|·\frac{f_s}{L},$$

(18)

in which if ${\widehat{f}}_{1,r}$ and ${\widehat{f}}_{1,r-1}$ are the fundamental frequency estimates in two subsequent data records and L/f_s represents the time interval between the corresponding UTC reference times.

Summarizing, the steps of the IpDFTc estimation algorithm are orderly reported in

Table 1. Steps of the IpDFTc estimation algorithm for PMUs.

Step 1: Apply the IpDFT algorithm based on the Hann window to the r-th data record yr(·) and determine the estimators ${\widehat{\nu }}_{10,r}, {\widehat{a}}_{10,r},\mathrm{and}{\widehat{\varphi }}_{10,r}$ by using (7) and (9). Step 2: Compute the residual signal ${\tilde{y}}_r(·)$ by removing the estimated sinewave from the acquired samples yr(·). Step 3: Apply the IpDFT algorithm based on the Hann window again to the residual signal ${\tilde{y}}_r(·)$ and determine the second harmonic parameters ${\widehat{\nu }}_{20,r}, {\widehat{a}}_{20,r}, \mathrm{and}{\widehat{\varphi }}_{20,r}.$ Step 4: Compute the new sine-wave parameter estimates ${\widehat{\nu }}_{1,r}, {\widehat{a}}_{1,r},\mathrm{and}{\widehat{\varphi }}_{1,r}$ by using (11), (15), and (16). Step 5: Compute $f_{1,r}$ and ${\widehat{ROCOF}}_{1,r}$ by using (17) and (18). Step 6: Return the estimated parameters.

.

2.3. Tuned Real-Valued TWLS Algorithm

The fundamental component of (1) can be rewritten as:

$$s\left(t\right)=a_1\left(t\right)\cos \left(2\pi f_{1}t+{\varphi }_1\left(t\right)\right)=c_1\left(t\right)\cos \left(2\pi f_{1}t\right)-s_1\left(t\right)\sin \left(2\pi f_{1}t\right),$$

(19)

where c₁(t) ₌ a₁(t) cos ${\varphi }_1\left(t\right)$ and s₁(t) ₌ a₁(t) sin ${\varphi }_1\left(t\right)$. To implement the real-valued TWLS algorithm [28], the two quadrature components c₁(t) and s₁(t) are approximated by their respective Taylor’s series about the reference time t_r truncated to the K-th order term (K ≥ 2). As a result,

$$c_1\left(t\right)\cong c_1\left(t_r\right)+c_1^{\left(1\right)}\left(t_r\right)\Delta t+\frac{c_1^{\left(2\right)}\left(t_r\right)}{2!}\Delta t^2+\dots +\frac{c_1^{\left(K\right)}\left(t_r\right)}{K!}\Delta t^K$$

(20)

and

$$s_1\left(t\right)\cong s_1\left(t_r\right)+s_1^{\left(1\right)}\left(t_r\right)\Delta t+\frac{s_1^{\left(2\right)}\left(t_r\right)}{2!}\Delta t^2+\dots +\frac{s_1^{\left(K\right)}\left(t_r\right)}{K!}\Delta t^K$$

(21)

where Δt = t − t_r and $c_1^{\left(k\right)}\left(t_r\right)$, $s_1^{\left(k\right)}\left(t_r\right)$, k = 1, 2,…, K, are the k-th order derivatives of c₁(t) and s₁(t), respectively, evaluated at the time instant t_r. The original Taylor-based weighted approach described in [30] assumes that the sinewave frequency coincides with the nominal one, i.e., f₁ = f_n in (19). However, possible off-nominal frequency deviations may significantly affect the accuracy of the estimation results. This is why a preliminary IpDFT-based fundamental frequency estimation is recommended in [29].

Let us consider the r-th data record collected in the observation interval $\left[\frac{rL-N_h}{f_s},\frac{rL+N_h}{f_s}\right]$ centered at the reference time $t_r=\frac{rL}{f_s}$. If (20) and (21) are replaced into (19) and the resulting expression is evaluated and equated to each sample of the data record, a linear system of M = 2N_h + 1 equations with 2K unknowns and coefficients $c_1\left(t_r\right)=c_{1,r}^{\left(0\right)},c_1^{\left(k\right)}\left(t_r\right)=\left(k!f_s^k\right)c_{1,r}^{\left(k\right)},s_1\left(t_r\right)=s_{1,r}^{\left(0\right)},s_1^{\left(k\right)}\left(t_r\right)=\left(k!f_s^k\right)s_{1,r}^{\left(k\right)},$ for k = 1, 2, …, K, can be built. If the unknown variables are gathered into vector $P_r\triangleq {\left[c_{1,r}^{\left(K\right)}\dots c_{1,r}^{\left(1\right)} c_{1,r}^{\left(0\right)}-s_{1,r}^{\left(0\right)}-s_{1,r}^{\left(1\right)}\dots -s_{1,r}^{\left(K\right)}\right]}^T$ and ${\widehat{f}}_{10,r}={\widehat{\nu }}_{10,r}\frac{f_s}{M}$ [with ${\widehat{\nu }}_{10,r}$ being given by (7)] is used to estimate f₁ in (19), the coefficients of the linear system can be arranged into a matrix V_K,r whose elements are [28]

$$\begin{array}{c}{\left(v_{K,1}\right)}_{lq}={\left(l-N_h\right)}^{K+1-q}\cos \left(2\pi {\widehat{\nu }}_{10,r}\left(N_h-l\right)\right),l=1,\dots , N_h,q=1,\dots , K+1\hfill \\ {\left(v_{K,2}\right)}_{lq}=l^{K+1-q}\cos \left(2\pi {\widehat{\nu }}_{10,r}l\right),l=1,\dots , N_h,q=1,\dots , K+1\hfill \\ {\left(v_{K,3}\right)}_{lq}=-{\left(l-N_h\right)}^{q-1}\sin \left(2\pi {\widehat{\nu }}_{10,r}\left(N_h-l\right)\right),l=1,\dots ,N_h,q=1,\dots , K+1\hfill \\ {\left(v_{K,4}\right)}_{lq}=l^{q-1}\sin \left(2\pi {\widehat{\nu }}_{10,r}l\right),l=1,\dots , N_h,q=1,\dots , K+1\hfill \end{array}$$

(22)

Thus, the weighted least squares approach returns the following estimate of vector $P_r$:

$$\widehat{P_r}={\left(V_{K,r}^T{W}^{T}W{V}_{K,r}\right)}^{-1}{V}_{K,r}^T{W}^{T}W{y}_r$$

(23)

where y_r represents the column vector including the samples of the r-th data record (2) and W is a M × M diagonal matrix. The nonzero elements of W are the coefficients of the window function w(·) that is used to smooth the approximation error near the ends of each observation interval. It is worth recalling that if K = 2 the coefficients of $\widehat{P_r}$ can be estimated not only in the time domain, but also in the frequency domain, as summarized in Appendix A. Due to the preliminary fundamental frequency tuning, the elements of $V_{K,r}$ have to be recomputed for every new record of collected data. For this reason, the complexity of the tuned-TWLS algorithm is generally greater than the computational burden of the IpDFTc estimator.

If the pairs of the real and the imaginary parts of the elements of vector $\widehat{P_r}$ returned by (23) are combined into complex numbers (i.e., ${\widehat{p}}_{k,r}={\widehat{c}}_{1,r}^{\left(k\right)}+j{\widehat{s}}_{1,r}^{\left(k\right)}$ for $k=0,1,2$), the synchrophasor magnitude and phase as well as the fundamental frequency and the ROCOF evaluated at time t_r can be obtained respectively as

$${\widehat{a}}_1=\left|{\widehat{p}}_{0,r}\right|,{\widehat{\varphi }}_1=angle\left\{{\widehat{p}}_{0,r}\right\},$$

(24)

$${\widehat{f}}_{1,r}=f_n+\frac{f_s}{2\pi }\frac{Im\left\{{\widehat{p}}_{1,r}{\widehat{p}}_{0,r}^{*}\right\}}{{\left|{\widehat{p}}_{0,r}\right|}^2} \mathrm{and}$$

(25)

$${\widehat{ROCOF}}_{1,r}=\frac{f_s^2}{\pi }\left[\frac{Im\left\{{\widehat{p}}_{2,r}{\widehat{p}}_{0,r}^{*}\right\}}{{\left|{\widehat{p}}_{0,r}\right|}^2}-\frac{Re\left\{{\widehat{p}}_{1,r}{\widehat{p}}_{0,r}^{*}\right\}Im\left\{{\widehat{p}}_{1,r}{\widehat{p}}_{0,r}^{*}\right\}}{{\left|{\widehat{p}}_{0,r}\right|}^4}\right].$$

(26)

Ultimately, the steps of the tuned real-valued TWLS estimation algorithm are orderly reported in

Table 2. Steps of the tuned real-valued TWLS estimation algorithm for PMUs.

Step 1: Apply the IpDFT algorithm based on the Hann window to the r-th data record yr(·) and determine the estimator ${\widehat{\nu }}_{10,r}$ from (7). Step 2: Compute the elements of the vector $\widehat{P_r}$ through (23). Step 3: Compute the estimators ${\widehat{a}}_{1,r}\mathrm{and}{\widehat{\varphi }}_{1,r}$ by (24), and the estimators of ${\widehat{f}}_{1,r}$ and ${\widehat{ROCOF}}_{1,r}$ with (25) and (26), respectively. Step 4: Return the estimated parameters.

.

2.4. Tuned Frequency Down-Conversion and Low-Pass Filtering (DCF)

As briefly outlined in the Introduction, the DCF approach is one of the most classic techniques for synchrophasor estimation. Indeed, it was mentioned in past editions of the Standard, too. The underlying idea is quite simple and it is very similar to the direct frequency down-conversion technique implemented in many radio receivers. If the samples of the r-th data record (2) are mixed with two digital quadrature sinewaves of unit amplitude and frequency f_1,r, the following in-phase (I) and quadrature (Q) components of the input signal are obtained from the prosthaphaeresis formulas:

$$\begin{array}{cc}\hfill y_r^I\left(m\right)=& y_r\left(m\right)\cos \left(2\pi \frac{f_{1,r}}{f_s}m\right)=\hfill \\ \hfill & =\frac{a_1\left(rL+m\right)}{2}\left[\cos \left({\varphi }_1\left(rL+m\right)\right)+\cos \left(4\pi \frac{f_{1,r}}{f_s}m+{\varphi }_1\left(rL+m\right)\right)\right]\hfill \\ \hfill & +{\eta }^I\left(rL+m\right)\hfill \\ \hfill y_r^Q\left(m\right)=& -y_r\left(m\right)\sin \left(2\pi \frac{f_{1,r}}{f_s}m\right)=\hfill \\ \hfill & =\frac{a_1\left(rL+m\right)}{2}\left[\sin \left({\varphi }_1\left(rL+m\right)\right)+\sin \left(4\pi \frac{f_{1,r}}{f_s}m+{\varphi }_1\left(rL+m\right)\right)\right]\hfill \\ \hfill & +{\eta }^Q\left(rL+m\right)\hfill \end{array}$$

(27)

where ${\eta }^I(·)$ and ${\eta }^Q(·)$ are the contributions due to harmonics, inter-harmonics and wideband noise resulting from the frequency down-conversion of term $\eta (·)$ in (2). Expression (27) shows clearly that, as a result of mixing, the power of the fundamental sinewave is evenly split into two components: one at DC, and another around 2f_1,r. Similarly, every harmonic and inter-harmonic component is split into pairs of contributions shifted by ±f₁,_r with respect to their respective original frequency. It is worth emphasizing that the frequency of the digital quadrature oscillators can be hardly set exactly equal to the fundamental frequency f_1,r. Therefore, in practice, either it is set equal to the nominal value f_n (as suggested in the Standard) or it is extracted from the input signal. This can be done, for instance, by estimating f_1,r with a phase-locked loop or with a preliminary IpDFT algorithm, as done in this paper [32].

Now, let $h\left(m\right)$ be the impulse response of a digital low-pass filter with a a null phase response and with magnitude response $\left|H\left(e^{j\omega }\right)\right|$ close to 2 in the passband (assuming a cutoff frequency $\omega <2\pi \frac{f_n}{f_s}$), and to 0 in the stopband, respectively. If this filter is applied in parallel to both $y_r^I\left(m\right)$ and $y_r^Q\left(m\right)$, it follows that

$$\begin{array}{c}{\tilde{y}}_r^I\left(m\right)=y_r^I\left(m\right)⊛h\left(m\right)=a_1\left(rL+m\right)\cos \left({\varphi }_1\left(rL+m\right)\right)+e_r^I\left(m\right)\hfill \\ {\tilde{y}}_r^Q\left(m\right)=y_r^I\left(m\right)⊛h\left(m\right)=a_1\left(rL+m\right)\sin \left({\varphi }_1\left(rL+m\right)\right)+e_r^Q\left(m\right).\hfill \end{array}$$

(28)

In (28) $⊛$ denotes the convolution operator, while $e_r^I\left(m\right)$ and $e_r^Q\left(m\right)$ are the error terms affecting the baseband I/Q components due to the filter non-idealities. These are due to (i) the frequency fluctuations of the quadrature sinewaves generated by the internal digital oscillator used for frequency down-conversion; (ii) the ripples affecting the gain of the low-pass filter in the passband; and (iii) the imperfect attenuation of both the double-frequency terms and all the disturbances lying in the stopband. However, if the filter is suitably designed, such error terms can be kept under control. Thus, from (28), it follows that synchrophasor magnitude and phase at time t_r can be easily estimated as

$${\widehat{a}}_{1,r}=\left|{\tilde{y}}_r^I\left(m\right)+j{\tilde{y}}_r^Q\left(m\right)\right|,{\widehat{\varphi }}_{1,r}=atan2\left({\tilde{y}}_r^Q\left(m\right),{\tilde{y}}_r^I\left(m\right)\right)$$

(29)

where atan2(y,x) computes the principal value of the phase of a complex number x + jy between [−π, π]. Moreover, the sinewave fundamental frequency and its ROCOF can be estimated by applying the first-order and second-order backward Euler differences over subsequent data records, i.e.,

$${\widehat{f}}_{1,r}=\left({\widehat{\varphi }}_{1,r}-{\widehat{\varphi }}_{1,r-1}\right)\frac{f_s^{}}{2\pi L}$$

(30)

and

$${\widehat{ROCOF}}_{1,r}=\frac{f_s^2}{2\pi L^2}\left({\widehat{\varphi }}_{1,r}-2{\widehat{\varphi }}_{1,r-1}+{\widehat{\varphi }}_{1,r-2}\right).$$

(31)

Of course, the performances of the DCF approach strongly depend on the features of the adopted low-pass filter as well as on the accuracy of the frequency tuning of the sinewave used for down-conversion. In fact, if the signal nominal frequency f_n rather than the estimated value of f_1,r is used for down-conversion, the overall estimation accuracy can be heavily affected by possible off-nominal frequency deviations of the acquired signal, unless the filter passband is enlarged accordingly [37]. On the contrary, the frequency tracking of the fundamental component can be exploited to place the zeros of the filter frequency response at harmonic frequencies to mitigate their influence [32]. A custom approach to design optimal differentiator filters for PMUs is presented for instance in [38], but it is applied to the Taylor-Fourier Transform and not to the DCF method.

In the next Section, the results obtained with two alternative low-pass filters are reported. The steps of the tuned DCF estimation technique are instead summarized in

Table 3. Steps of the tuned DCF estimation algorithm for PMUs.

Step 1: Apply the IpDFT algorithm based on the Hann window to the r-th data record yr(·) and determine the estimator ${\widehat{\nu }}_{10,r}$ from (7). Step 2: Mix the samples of the data record yr(·) with two quadrature sinewaves of unit amplitude and frequency ${\widehat{f}}_{1,r}={\widehat{\nu }}_{10,r}·f_s$/M to obtain the I/Q components of the collected signal. Step 3: Apply a suitable low-pass filter to the I/Q components to remove disturbances and the double-frequency terms. Step 4: Compute the estimators ${\widehat{a}}_{1,r}\mathrm{and}{\widehat{\varphi }}_{1,r}$ by (29), and the estimators ${\widehat{f}}_{1,r}$ and ${\widehat{ROCOF}}_{1,r}$ with (30) and (31), respectively. Step 5: Return the estimated parameters.

.

3. Estimation Results and Performance Analysis

In this Section, the performances of the IpDFT-tuned algorithms presented in Section 2 are compared under multiple testing conditions. In particular, the maximum values of the following parameters are computed, i.e.,

• the Total Vector Errors (TVEs), namely the magnitude of the error vectors of the estimated synchrophasor relative to the actual synchrophasor values at the same UTC reference times;
• the Frequency Errors (FE), i.e., the differences between the estimated fundamental frequencies and the actual ones and
• the Rate of change of Frequency Errors (RFEs) namely the differences between the estimated ROCOF values and the actual ones.

Such parameters are estimated assuming L = 120 and a reporting rate RR = 50 readings/s. This is the largest mandatory rate reported in the Standard [10]. To this purpose, the results of the IpDFT-tuned estimation algorithms under study are compared over 1000 data records in the P Class and M Class steady-state and dynamic testing conditions specified in [10]. The initial phase values of all possible disturbances (e.g., harmonics, inter-harmonics or modulating tones) are varied randomly with a uniform probability density function in the range [0, 2π]. Then, the TVE, |FE| and |RFE| response times as well as the corresponding delay times and the maximum overshoots due to amplitude or phase steps in the worst-case P Class and M Class simulated transient conditions reported in [10] are summarized. All the simulated P Class and M Class tests are performed assuming the following common settings, i.e., nominal fundamental frequency f_n = 50 Hz, sampling frequency f_s = 6 kHz and observation interval length ranging from J₁ = 2 or 4 nominal power line cycles for P Class tests to J₁ = 6 or 8 nominal power line cycles for M Class tests.

Such observation interval values were chosen as they provide a good trade-off between estimation accuracy and responsiveness, although with actual performances that may differ considerably from case to case. The preliminary IpDFT stage (common to all algorithms) is based on the 2-term MSD (or Hann) window defined in [35] with two DFT interpolation points [17]. Quite importantly, the TVE, |FE|, and |RFE| values achievable with the IpDFTc algorithm in steady-state tests with off-nominal frequency and harmonics are much lower than those obtained with the classic IpDFT estimator, as confirmed by Figure 1. In the other testing conditions instead, the results of the IpDFTc and the IpDFT almost coincide. Therefore, in the rest of this paper the latter ones are omitted for the sake of brevity. In Figure 1 the curves with different line styles and colors refer to the maximum TVE, |FE| and |RFE| values returned by the classic IpDFT and the IpDFTc estimators as a function of the off-nominal frequency deviation. Such maxima are computed over 1000 6-cycle-long or 8-cycle-long records, assuming that the first 50 harmonics (all of them with amplitude equal to 10% of the fundamental component) affect the collected waveform. The benefits of the compensation of the leakage due to both the second-order harmonic and the fundamental image component are evident (particularly in the 6-cycle case) and justify the adoption of the IpDFTc approach in the following analysis. Further results (not shown for the sake of brevity) confirm that the IpDFTc algorithm can meet the P Class limits reported in the Standard even over 3-cycle observation intervals, whereas the classic IpDFT estimator cannot.

3.1. P Class and M Class Steady-State and Dynamic Testing Conditions

The main P Class and M Class steady-state and dynamic testing conditions specified in [10] are briefly recalled below.

• Fundamental frequency off-nominal static deviations within ±2 Hz (P Class) or ±5 Hz (M Class);
• Harmonics from the 2nd to the 50th (taken one at a time) with amplitude set to 1% (P Class) or 10% (M Class) of the fundamental with the same off-nominal frequency deviations mentioned above;
• Amplitude modulation (AM) or phase modulation (PM) with modulation indexes equal to k_a = 0.1 or k_p = 0.1 rad, respectively, and frequency of the modulating tones equal to f_a = f_p = 2 Hz (P Class) or 5 Hz (M Class);
• Linear increment/decrement of the fundamental frequency at rates of ±1 Hz/s within ±2 Hz (P Class) or ±5 Hz (M Class);
• White zero-mean Gaussian noise with variance such that the Signal-to-Noise Ratio (SNR) is 60 dB (this testing condition is not mentioned in the IEEE/IEC Standard, but it has a high practical relevance due to the strong impact of wideband noise on frequency and above all ROCOF estimation [39]);
• Finally, in the M Class case only, a single inter-harmonic tone of amplitude equal to 10% of the fundamental and frequency changing within the bands [10 Hz, 25 Hz] or [75 Hz, 100 Hz] is generated for testing, assuming that the off-nominal fundamental frequency deviations are comprised within [−2.5 Hz, 2.5 Hz] and that RR = 50 readings/s.

Table 4. Maximum TVE, |FE|, and |RFE| values obtained with the estimation algorithms for PMUs described in Section 2 over 2-cycle (a) or 4-cycle (b) observation intervals. All the testing conditions, except the one with the additive wideband noise, are specified in the IEEE/IEC Standard 60255-118-1:2018 for P Class PMUs [10]. The limits mentioned in the Standard (when defined) are also reported. The values exceeding such limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded.

(a)
Test Type	TVEmax(%)					|FE|max(mHz)					|RFE|max(Hz/s)
	Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
				LSF	STD				LSF	STD				LSF	STD
Freq. dev. only (±2 Hz)	1	0.01	0.00	0.38	0.17	5	8.9	17.7	13.5	7.2	0.4	0.2	0.0	8.1	4.2
Freq. dev.+ 1% harmonics	1	0.01	0.74	0.38	0.17	5	10.6	141	14.1	7.5	0.4	0.3	54.7	8.4	4.8
AM (fa = 2 Hz, ka = 10%)	3	0.04	0.00	0.36	0.05	60	8.0	0.0	0.0	0.0	2.3	0.1	0.0	0.0	0.0
PM (fp = 2 Hz, kp = 0.1 rad)	3	0.04	0.00	0.37	0.03	60	1.8	0.5	1.6	0.7	2.3	0.0	0.0	0.6	0.0
Freq. ramp 1 (±2 Hz @ ±1 Hz/s)	1	0.02	0.00	0.40	0.19	10	5.4	10.8	13.9	7.2	0.4	0.1	0.0	8.0	4.2
AWGN (SNR = 60 dB)	-	0.03	0.04	0.36	0.03	-	14.1	24.3	2.7	2.5	-	0.6	1.4	0.3	0.3
(b)
Test Type	TVEmax(%)					|FE|max(mHz)					|RFE|max(Hz/s)
	Lim.	IpDFTc	TunedTWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
				LSF	STD				LSF	STD				LSF	STD
Freq. dev. only (±2 Hz)	1	0.00	0.03	0.06	0.16	5	0.1	0.1	0.4	6.3	0.4	0.0	0.0	0.3	3.8
Freq. dev.+ 1% harmonics	1	0.00	0.01	0.06	0.16	5	0.1	1.4	0.4	6.5	0.4	0.0	0.3	0.3	4.0
AM (fa = 2 Hz, ka = 10%)	3	0.15	0.00	0.09	0.15	60	14.3	0.0	0.0	0.1	2.3	0.2	0.0	0.0	0.3
PM (fp = 2 Hz, kp = 0.1 rad)	3	0.16	0.00	0.06	0.12	60	2.1	1.8	0.4	2.0	2.3	0.0	0.0	0.0	0.0
Freq. ramp 1 (±2 Hz @ ±1 Hz/s)	1	0.07	0.03	0.06	0.21	10	0.1	0.1	0.5	6.3	0.4	0.0	0.0	0.4	3.8
AWGN (SNR = 60 dB)	-	0.02	0.03	0.09	0.02	-	2.9	3.5	2.9	2.2	-	0.2	0.3	0.3	0.1

¹ Results obtained considering an exclusion interval of 40 ms (i.e., twice the maximum mandatory reporting period at 50 Hz) at the beginning and at the end of the frequency ramp.

a,b shows the maximum values of TVE, |FE| and |RFE| (over 1000 records) defined as reported in the Standard and obtained over 2-cycle (a) or 4-cycle (b) observation intervals in the P Class testing conditions summarized above. The P Class limits mentioned in the Standard (when defined) are also reported. The values exceeding such limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded. In the case study considered, both the IpDFTc and the tuned TWLS algorithm rely on the Hann window. In the DCF case instead, two different sets of results obtained using two different filters are shown. The first set of DCF results is obtained with 2-cycle-long and 4-cycle-long linear-phase Finite Impulse Response (FIR) filters designed by minimizing the least-squares errors with respect to the frequency response of an ideal low-pass filter. These filters will be shortly labeled as Least-Squares Filters (LSF) in the following. The second set of DCF results is obtained instead with the simple filter suggested in Annex D of the Standard [10], i.e., a filter with a triangular impulse response denoted simply as standard (STD) P Class filter in the rest of this manuscript. Quite importantly, also the classic equiripple FIR filters resulting from Parks and McClellan minimax optimization were considered in the analysis. However, they led to worse performances in almost all testing conditions. Therefore, the results obtained with such filters are omitted for the sake of brevity. This is probably due to the fact that, unlike the LSF and STD filters, the constant ripple magnitude in the stopband cannot be further reduced unless the order of the filter is increased. As a consequence, possible crucial disturbances (e.g., large high-order harmonics) with frequencies close to the position of the peaks of the stopband ripples are not attenuated enough.

The results in Table 4a,b can be summarized as follows.

• As far as synchrophasor estimation is concerned, all algorithms return maximum TVE values compliant with the Standard even over observation intervals so short as two nominal power line cycles. However, the IpDFTc and the tuned-TWLS algorithms provide the lowest TVE values. In particular, the IpDFTc algorithm is generally preferable in steady-state testing conditions and in the presence of wideband noise, whereas the tuned TWLS estimator is superior when dynamic disturbances (i.e., modulations and frequency ramps) are considered.
• The estimation of the fundamental frequency is quite problematic. If just two cycles are observed, no estimator is able to meet the P Class |FE| limits reported in the Standard in all the testing conditions. In fact, further simulations showed that at least 3-cycle observation intervals are needed. In this scenario, the DCF approach with the STD filter seems to be globally preferable. If 4-cycle observation intervals are used, all estimators return results compliant with the Standard. The |FE| values are indeed comparable in most cases. However, the tuned TWLS approach appears slightly better than the others (and particularly much better than the IpDFTc under the effect of modulations).
• Finally, as far as ROCOF estimation is concerned, the IpDFTc estimator is the only one that is able to meet the P Class |RFE| limits reported in the Standard by using both 2-cycle and 4-cycle observation intervals. The tuned TWLS approach over 2-cycle intervals returns accurate results in most cases, but it is strongly affected by the presence of harmonics. On the other hand, if 4-cycle observation intervals are used, the IpDFTc, the tuned TWLS, and the DCF approach with the LSF filter exhibit comparable performances.

Table 5. Maximum TVE, |FE| and |RFE| values obtained with the estimation algorithms for PMUs described in Section 2 over 6-cycle (a) and 8-cycle (b) observation intervals. All the testing conditions, except the one with the additive wideband noise, are specified in the IEEE/IEC Standard 60255-118-1:2018 for M Class PMUs [10]. The limits mentioned in the Standard (when defined) are also reported. The values exceeding such limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded.

(a)
Test Type	TVEmax(%)					|FE|max(mHz)					|RFE|max(Hz/s)
	Lim.	IpDFTc	TunedTWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
				LSF	STD				LSF	STD				LSF	STD
Freq. dev. only (±5 Hz)	1	0.00	0.00	0.33	0.13	5	0.0	0.0	17.0	8.4	0.4	0.0	0.0	9.8	4.7
Freq. dev.+ 10% harmonics	1	0.00	0.02	0.34	0.13	5	0.1	0.2	21.1	10.2	0.4	0.0	0.3	12.1	6.5
AM (fa = 5 Hz, ka = 10%)	3	2.36	0.07	0.49	0.89	300	140	5.2	1.5	2.3	14	4.2	0.1	3.4	8.1
PM (fp = 5 Hz, kp = 0.1 rad)	3	2.12	0.07	0.32	0.69	300	70.4	59.0	9.6	36.4	14	2.1	1.7	1.5	1.9
Freq. ramp 1 (±5 Hz @ ±1 Hz/s)	1	0.15	0.00	0.29	0.20	10	0.0	0.1	1.2	9.0	0.4	0.0	0.0	0.1	4.0
AWGN (SNR = 60 dB)	-	0.02	0.03	0.33	0.02	-	1.3	1.6	14.1	1.0	-	0.1	0.1	9.1	0.1
10% inter-harmonics	1.3	0.08	2.32	0.57	0.54	10	25.5	100	57.6	70.9	-	2.5	23.6	12.6	12.3
(b)
Test Type	TVEmax(%)					|FE|max(mHz)					|RFE|max(Hz/s)
	Lim.	IpDFTc	TunedTWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
				LSF	STD				LSF	STD				LSF	STD
Freq. dev. only (±5 Hz)	1	0.00	0.00	0.34	0.11	5	0.0	0.00	0.7	9.9	0.4	0.0	0.0	0.4	5.7
Freq. dev.+ 10% harmonics	1	0.00	0.00	0.34	0.12	5	0.0	0.2	0.7	11.4	0.4	0.0	0.0	0.4	6.5
AM (fa = 5 Hz, ka = 10%)	3	3.92	0.21	0.43	1.46	300	164.7	9.9	0.5	4.2	14	4.9	0.1	2.2	17.8
PM (fp = 5 Hz, kp = 0.1 rad)	3	3.51	0.20	0.35	1.21	300	119.2	100.3	4.1	59.9	14	3.6	2.8	0.2	2.4
Freq. ramp 1 (±5 Hz @ ±1 Hz/s)	1	0.26	0.00	0.37	0.20	10	0.1	0.1	1.2	9.9	0.4	0.0	0.0	0.4	5.8
AWGN (SNR = 60 dB)	-	0.02	0.02	0.35	0.20	-	0.9	1.2	0.8	0.7	-	0.1	0.1	0.5	0.0
10% inter-harmonics	1.3	0.02	0.02	0.38	0.12	10	3.6	6.5	7.7	7.8	-	0.4	0.3	1.3	3.0

¹ Results obtained considering an exclusion interval of 140 ms (i.e., 7 times the maximum mandatory reporting period at 50 Hz) at the beginning and at the end of the ramp.

a,b shows the maximum values of TVE, |FE| and |RFE| (computed over 1000 records) obtained over 6-cycle (a) or 8-cycle (b) observation intervals in the M Class testing conditions listed at the beginning of this Section. In this case, longer observation intervals are considered as they are needed to ensure a higher measurement accuracy. The M Class limits mentioned in the Standard (when defined) are also shown. The values exceeding such limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded. Again, both the IpDFTc and the tuned TWLS algorithm rely on a Hann window. The LSF FIR filters used with the DCF approach is again based on least-squares frequency response error minimization assuming 6-cycle-long and 8-cycle-long impulse responses, respectively. On the contrary, the M Class STD filter suggested in Annex D of the Standard does no longer exhibit a triangular impulse response. In fact, it results from the classic window design method: the non-causal sinc-shaped impulse response of the ideal low-pass filter is indeed modulated by a Hamming window.

The results in Table 5a,b can be summarized as follows.

• As far as synchrophasor estimation is concerned, the maximum TVE values basically confirm the conclusions drawn in the P Class case. The IpDFTc algorithm is generally preferable in steady-state testing conditions and in the presence of wideband noise, whereas the tuned TWLS estimator is superior when dynamic disturbances (i.e., modulations and frequency ramps) are considered. However, in the M Class case, such differences are exacerbated. The tuned TWLS approach is not able to meet the Standard limits in the presence of out-of-band inter-harmonics over 6-cycle intervals, but it meets the Standard requirements over 8-cycle intervals. Quite interestingly, the accuracy of the tuned TWLS algorithm over longer intervals under dynamic conditions slightly degrades due to the larger Taylor’s series approximation errors of the phasor model at the edges of the observation intervals [29]. The IpDFTc estimator is not able to meet the Standard limits over 8-cycle intervals when amplitude or phase oscillations perturb the fundamental component. Only the DCF approach is able to meet the Standard limits with both interval lengths, although in some testing conditions it exhibits a slightly worse accuracy than the other algorithms.
• The |FE| values exhibit a trend similar to the TVE ones. The main difference is in the case of the DCF approach whose behavior is more heterogeneous. In this case, ensuring compliance with the Standard requirements is much more challenging. Indeed, no estimator is able to meet the Standard limits in the case of inter-harmonic interference over 6-cycle intervals. Over 8-cycle intervals instead, the IpDFTc estimator returns generally the best results except in the case of amplitude and phase modulations, which are more effectively mitigated by the DCF approach with the LSF filter.
• Finally, in the case of ROCOF estimation, almost all algorithms except the DCF approach with the STD filter exhibit good results. They are indeed compliant with the limits reported in the Standard (when defined). Globally, the IpDFTc estimator is preferable over 6-cycle intervals, whereas the tuned TWLS algorithm is better over 8-cycle intervals in accordance with the TVE trend.

3.2. Results in Transient Conditions

To compare the behavior of the three considered IpDFT-tuned estimation algorithms in transient conditions, both the response times associated with synchrophasor, frequency and ROCOF estimation, and the measurement delay times are computed in the worst-case conditions specified in the Standard [10]. In this case, the testing conditions are: either amplitude step changes equal to ±10% of the nominal amplitude or phase step changes of ±10°. It is worth recalling that the response times are defined as the temporal intervals between the instants at which, as a result of a step input change, the measured values of synchrophasor, frequency or ROCOF exceed a specified TVE, |FE| or |RFE| threshold (i.e., 1%, 5 mHz or 0.4 Hz/s, respectively, in the case of P Class PMUs, or 1%, 5 mHz or 0.1 Hz/s in the case of M Class PMUs) and the times after which the measured values steadily remain under such thresholds.

During the transients, the maximum over- or undershoots affecting amplitude or phase estimates (depending on the type of applied input step) are also of interest, as they are required to be smaller than 5% or 10% of the input step size for P Class and M Class PMUs, respectively.

The measurement delay time is defined as “the time interval between the instant that a step change is applied to the input of a PMU and the measurement time when the stepped parameter achieves a value that is half-way between the initial and final steady-state values” [10].

In the present study, the P Class and M Class maximum response times and the maximum measurement delay times (obtained by shifting sample-by-sample the instants at which the steps occur) are shown in

Table 6. Maximum response times (expressed in nominal power line cycles) associated with synchrophasor, frequency and ROCOF estimation (a) and maximum delay times (b), under the influence of the worst-case amplitude or phase step changes specified in the IEEE/IEC Standard 60255-118-1:2018 for P Class PMUs using 2-cycle or 4-cycle observation intervals. The limits specified in the Standard are also shown. The values exceeding the limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded.

(a)
Test Type	Obs. Inter. Length [cycles]	Synchrophasor Resp. Time (cycles)					Frequency Resp. Time (cycles)					ROCOF Resp. Time (cycles)
		Lim	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
					LSF	STD				LSF	STD				LSF	STD
±10% amp. step	2	2	1.05	0.56	1.23	1.18	4.5	1.84	1.88	9.15	2.76	6	2.89	1.80	9.32	2.98
	4	2	1.99	0.94	1.04	2.29	4.5	3.56	3.66	2.73	4.45	6	4.68	3.43	7.97	4.91
±10° phase step	2	2	1.25	1.00	1.35	1.43	4.5	1.85	1.88	9.25	2.88	6	2.88	1.84	9.35	2.99
	4	2	2.35	1.10	1.20	2.68	4.5	3.47	3.61	4.33	4.74	6	4.64	3.58	8.18	4.98
(b)
	Test Type		Obs. Inter. [cycles]			|Delay Time| (ms)
						Lim.		IpDFTc	Tuned TWLS		DCF
											LSF		STD
	±10% amp. step		2			5		1.50	1.67		3.00		1.80
			4			5		1.50	1.50		1.80		1.80
	±10° phase step		2			5		1.83	2.00		4.70		2.80
			4			5		1.83	1.83		2.00		2.50

a,b over 2-cycle and 4-cycle observation intervals and in

Table 7. Maximum response times (expressed in nominal power line cycles) associated with synchrophasor, frequency and ROCOF estimation (a) and maximum delay times (b), under the influence of the worst-case amplitude or phase step changes specified in the IEEE/IEC Standard 60255-118-1:2018 for M Class PMUs using 6-cycle or 8-cycle observation intervals. The limits specified in the Standard are also shown. The values exceeding the limits are highlighted in bold, whereas the cells containing the best results obtained in each test are shaded.

(a)
Test Type	Obs. Inter. (cycles)	Synchrophasor Resp. Time (cycles)					Frequency Resp. Time (cycles)					ROCOF Resp. Time (cycles)
		Lim	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF		Lim.	IpDFTc	Tuned TWLS	DCF
					LSF	STD				LSF	STD				LSF	STD
±10% amp. step	6	7	2.96	1.36	1.68	1.98	14	5.17	5.37	5.78	4.14	14	6.43	4.90	7.33	12.27
	8	7	3.93	1.78	1.81	2.60	14	6.71	6.99	4.37	4.96	14	8.09	6.17	8.98	8.91
±10° phase step	6	7	3.49	1.56	1.94	2.32	14	5.13	5.31	7.13	6.36	14	6.36	5.22	7.33	12.28
	8	7	4.63	2.04	2.17	3.06	14	6.66	6.96	7.11	7.23	14	8.05	6.78	8.98	8.91
(b)
	Test Type		Obs. Inter. [cycles]			|Delay Time| (ms)
						Lim.		IpDFTc	Tuned TWLS		DCF
											LSF		STD
	±10% amp. step		6			5		1.50	1.50		1.50		1.50
			8			5		1.50	1.50		2.20		2.00
	±10° phase step		6			5		1.83	1.83		4.80		4.50
			8			5		1.83	1.83		3.00		4.70

a,b over 6-cycle and 8-cycle observation intervals, respectively.

The results in both cases are quite consistent and can be summarized as follows:

• The limits reported in the Standard are met in almost all cases with a few exceptions.
• The maximum synchrophasor and ROCOF response times of the TWLS estimator are always considerably shorter than those of the IpDFTc estimator. However, the IpDFTc frequency estimation response times are generally the shortest ones, except under the effect of amplitude steps when 4-cycle observation intervals are considered. In fact, this is the sole testing condition in which the DCF approach outperforms the other estimators. Furthermore, the over- and undershoots obtained with the IpDFTc estimator are generally smoother than those of the other algorithms. In any case, the over- and undershoots of all considered algorithms are well below the limits reported in the Standard. Therefore, they are not reported for the sake of brevity.
• As far as the measurement delay times are concerned, the tuned TWLS and the IpDFTc estimators exhibit similar performances, quite better than those of the DCF method. Moreover, they are well below the limits reported in the Standard [10].

4. Discussion

Table 8. Overall comparison between the IpDFT-tuned estimation algorithms under study in the P Class (a) and M Class (b) PMU testing conditions reported in the IEEE/IEC Standard 60255-118-1:2018. Letter C stands for “full compliance” with all the limits reported in the testing conditions considered in Section 3 for RR = 50 frame/s, whereas “not compliant” (NC) means that at least one of the Standard limits is exceeded. The cells related to the best performing algorithm are shaded.

(a)
Estimation Algorithm	No. Cycles	TVE	|FE|	|RFE|	Phasor Resp. Time	Freq. Resp. Time	ROCOF Resp. Time	Delay Time
IpDFTc	2	C	NC	C	C	C	C	C
	4	C	C	C	NC	C	C	C
Tuned TWLS	2	C	NC	NC	C	C	C	C
	4	C	C	C	C	C	C	C
DCF with LSF	2	C	NC	NC	C	NC	NC	C
	4	C	C	C	C	C	NC	C
DCF with STD	2	C	NC	NC	C	C	C	C
	4	C	NC	NC	NC	NC	C	C
(b)
Estimation Algorithm	No. Cycles	TVE	|FE|	|RFE|	Phasor Resp. Time	Freq. Resp. Time	ROCOF Resp. Time	Delay Time
IpDFTc	6	C	NC	C	C	C	C	C
	8	NC	C	C	C	C	C	C
Tuned TWLS	6	NC	NC	C	C	C	C	C
	8	C	C	C	C	C	C	C
DCF with LSF	6	C	NC	NC	C	C	C	C
	8	C	C	C	C	C	C	C
DCF with STD	6	C	NC	NC	C	C	C	C
	8	C	NC	NC	C	C	C	C

a,b summarizes the results of the performance analysis in the P Class and M Class testing conditions reported in Section 3. The estimation algorithms with TVE, |FE| and |RFE| values and with response times (in the case of amplitude or phase step changes) below the limits reported in the Standard in all the considered testing conditions are denoted as Compliant (C). If at least one of the TVE, |FE|, |RFE| or response time limits is exceeded, an estimator is labeled as Not Compliant (NC). Again, the algorithm with the best performances in estimating each parameter is highlighted by shading the corresponding cell in Table 8a,b. The criteria adopted to identify the best estimator of each parameter for a given observation interval length are the following:

• an algorithm able to return fully compliant TVE, |FE|, |RFE| or response time values is considered to be better than another one exceeding even a single limit, regardless of the estimation errors associated with the individual tests;
• if two algorithms exhibit a different number of C and NC tests, the one with a larger number of compliant results in different tests is considered the better one, regardless of the estimation errors in the individual tests;
• if two algorithms exhibit the same amount of compliant or non-compliant tests, the algorithm with a higher number of more accurate results shown in Table 4, Table 5, Table 6 and Table 7 is regarded as the better one;
• if two algorithms exhibit the same amount of compliant or non-compliant tests and they do not clearly outperform one another (e.g., because results are very similar or because the number of testing conditions in which they are better is the same), the estimators are considered to have an equivalent performance.

From the analysis of Table 8a it results clearly that in the P-Class testing conditions full compliance is achieved only with the Tuned TWLS estimator over 4-cycle intervals. However, this is also the most demanding algorithm from the computational point of view. The IpDFTc estimator requires a lower computational burden and it is globally the best one over 2-cycle intervals, although it suffers from excessive |FE| values in steady-state testing conditions. Over 4-cycle intervals instead, despite its high accuracy, the main problem of the IpDFTc algorithm is the excessive response time duration affecting synchrophasor estimation after step changes.

In the M-Class testing conditions summarized in Table 8b the IpDFTc estimator is globally the best one over 6-cycle observation intervals, but it again suffers from excessive |FE| values when out-of-band inter-harmonics are considered. On the contrary, the Tuned TWLS algorithm is confirmed to be both fully compliant and the best estimator over 8-cycle observation intervals. Quite interestingly, the DCF technique with the LSF filter over 8-cycle intervals ensures full compliance in all tests considered, despite its accuracy and responsiveness not being excellent in any test.

5. Experimental Results in a Practical Case Study

To validate the simulation-based performance comparison reported above in a case study of practical interest, 1000 2-cycle, 4-cycle, 6-cycle and 8-cycle data records of a waveform generated by a function generator Agilent 33220A were collected through a 12-bit data acquisition system (DAQ) NI6023E. Each generated waveform exhibits nominal normalized amplitude, initial phase chosen at random in the range [0, 2π], fundamental frequency within [47 Hz, 52 Hz] and a Total Harmonic Distortion (THD) of about 9.4%. The THD level purposely slightly exceeds the limit reported in the Standard EN 50160:2010 [33]. In particular, according to that Standard, the range [47 Hz, 52 Hz] shall never be exceeded, although in practice the system frequency is required to lie within [49.5 Hz, 50.5 Hz] for 95% of the time of a week. The THD instead must be not greater than 8% with the first 25 harmonics that have to be smaller than a value ranging from 0.5% to about 6% of the fundamental, depending on the harmonic order. In the current study, just the harmonics from the 2nd to the 7th ones are included. Their amplitude is equal to 2%, 5%, 1%, 6%, 0.5% and 5% of the fundamental, respectively. The values of amplitude, phase, frequency and ROCOF evaluated at the reference times were estimated a posteriori through a multi-harmonic sine-fitting procedure [40].

The bar diagram in Figure 2 shows the maximum TVE, |FE| and |RFE| values computed over all data records and associated with the IpDFTc algorithm, the tuned TWLS estimator and the DCF technique using either the LSF or the STD filter, respectively. Given that, as shown in Section 3 and Section 4, the best results as well as the full compliance of some algorithms are obtained over 4-cycle intervals in the P Class case and over 8-cycle intervals in the M Class case, respectively, the results obtained over 2-cycle and 6-cycle intervals are omitted here for the sake of brevity. Moreover, the TVE, |FE| and |RFE| values obtained over 2 cycle intervals are very large due to the large spectral leakage caused by the 2nd and 3rd harmonic.

The results in Figure 2 confirm that, despite the presence of a variety of uncertainty contributions (most notably the limited spectral purity of the adopted waveform generator, the limited effective resolution of the DAQ and the uncertainty associated with the reconstruction of the reference values of amplitude, phase, frequency and ROCOF) the behavior of the estimation algorithms under test is consistent with the simulation results shown in Table 4b and Table 5b. In particular, it is evident that:

• the IpDFTc estimator exhibits the best accuracy both over 4-cycle and 8-cycle observation intervals, as expected since it is natively conceived to counteract the effect of the 2nd harmonic;
• the tuned TWLS returns slightly worse results, but it is worth recalling that it behaves better in dynamic and transient testing conditions;
• the DCF technique is outperformed with either filter although the results obtained using a 4-cycle-long LSF are not much worse than those achieved with the IpDFTc and the tuned TWLS estimators. Note that the TVE values of the DCF technique with the LSF over 8-cycle intervals are quite worse than those over 4 cycle-intervals. However, this behavior is consistent with the corresponding simulation results shown in Table 4b and Table 5b and this is probably due to the larger bandwidth of the LSF with an 8-cycle-long impulse response compared with the LSF with a 4-cycle-long impulse response. Indeed, while the latter filter is conceived for P Class applications, the former one is optimized not to perturb stronger amplitude and phase oscillations in the M Class case. Thus, the LSF with an 8-cycle-long impulse response makes the DCF technique more sensitive to noise (e.g., injected by the experimental tested).

6. Conclusions

Recent research results have shown that if a preliminary fundamental frequency measurement is performed to compensate for the effect of static off-nominal frequency deviations, the accuracy of some of the most common estimation algorithms for PMUs can be greatly improved. Till now, to the best of Authors’ knowledge, no clear performance comparisons between estimation algorithms exploiting the benefits of a preliminary frequency measurement have been reported in the scientific literature. In this paper, the performances of three alternative techniques for the estimation of AC waveform amplitude, phase, frequency and ROCOF based on an a prior 2-point Interpolated Discrete Fourier Transform (IpDFT) stage are evaluated in the very same P Class or M Class testing conditions specified in the IEEE/IEC Standard 60255-118-1:2018 and using data records of the same size. The algorithms considered in this study are: the corrected IpDFT (IpDFTc), the tuned Taylor Weighted Least Squares (TWLS) estimator, and the frequency Down-Conversion and low-pass Filtering (DCF) technique with two alternative types of low-pass filter. The performance comparison shows that the IpDFTc and the tuned TWLS methods are generally more accurate and exhibit shorter transients than the DCF technique, although the latter one still fulfills the Standard requirements in most testing conditions. Quite importantly, the behavior of the DCF algorithm is strongly affected by the chosen low-pass filter and it could be further improved by adopting custom optimization strategies for FIR filter design. The IpDFTc technique is preferable over short intervals and it overcomes the others when the analyzed waveform is affected by steady-state harmonic and inter-harmonic disturbances. However, it is considerably worse than the tuned TWLS estimator in the presence of amplitude or phase oscillations and it is not able to ensure compliance to the IEEE/IEC Standard in a few testing conditions. On the contrary, the tuned TWLS estimator is much more effective in dynamic conditions, it exhibits shorter synchrophasor and ROCOF response times (although with greater overshoots) and, assuming a reporting rate of 50 readings/s, it returns fully compliant results in the P Class and M Class tests described in the IEEE/IEC Standard over 4-cycle and 8-cycle observation intervals, respectively.