Low-Latency, Three-Phase PMU Algorithms: Review and Performance Comparison

Featured Application

Algorithms for synchrophasor, frequency and rate of change of frequency estimation are essential to implement effective Phasor Measurement Units. In particular, the three-phase characteristics of AC power systems can be favorably exploited for measuring the positive sequence contribution. The present paper introduces the mathematical framework of the Space Vector transformation and compares three estimation algorithms based on this approach, characterized by reporting latencies lower than two nominal cycles, thus suitable for relaying and control applications.

Abstract

Phasor Measurement Units are the most advanced instruments for power network monitoring, since they allow phasors, frequency and rate of change of frequency (ROCOF) to be measured in predetermined time instants with respect to an absolute time reference. The employed estimation algorithm plays a key role in overall performance under off-nominal conditions; the challenge to be faced is combining high steady-state accuracy with fast responsiveness in dynamic conditions, small reporting latency and reduced computational burden. Under regular operation, AC power networks are weakly unbalanced three-phase systems. Based on this consideration, the recent literature has proposed native three-phase estimation algorithms that effectively exploit this property to accurately identify the positive sequence synchrophasor, frequency and ROCOF. In this respect, the present paper describes three among the most promising three-phase algorithms based on the Space Vector transformation. By means of numerical simulations, it compares the achieved performance in terms of response time and estimation error both under steady-state and dynamic conditions. All the considered approaches enable a flexible design that allows balancing accuracy and responsiveness. For this analysis, the reporting latency has been limited to about one and half nominal cycles, i.e., 30 ms at 50 Hz; the P-class algorithm suggested by IEC/IEEE Std 60255-118-1 has also been included as comparison benchmark.

1. Introduction

In recent years, modern power networks are experiencing an ever increasing integration of renewable energy sources and distributed generation, characterized by higher volatility and faster dynamics [1,2]. In this context, the measurement infrastructure plays a crucial role in guaranteeing the operability and continuity of service of the power network [3]. In particular, an optimal tradeoff between high measurement accuracy and low reporting latency is required in many monitoring and control applications, like fault detection [4] or load shedding [5].

In such challenging scenario, Phasor Measurement Units (PMUs) might represent a promising solution as they provide estimates of synchrophasors, frequency and rate of change of frequency (ROCOF) associated to the fundamental component in remote nodes of the grid, which are synchronized to Coordinated Universal Time (UTC) [6]. The recent IEEE/IEC Std 60255-118-1-2018 [7] (denoted as IEC Std for the sake of brevity), defines the operative requirements for PMUs in terms of measurement accuracy, latency and transient performance. The measurement accuracy is quantified by means of three indices, namely Total Vector Error (TVE), Frequency Error (FE), and ROCOF Error (RFE). The reporting latency, instead, represents the time delay after which measurement data is available at the PMU output with respect to the corresponding reporting instant (timestamp). Transient performance is measured by the response time required by the algorithm output to settle down after fast variations in amplitude or phase and by such response shape parameters. In particular, the IEC Std introduces two performance classes: M for high-accuracy monitoring applications, and P for low-latency applications where high responsiveness is recommended, such as relaying and fast control.

The requirements of the IEC Std, as well as the superseded IEEE standards [8,9], have pushed the development of highly sophisticated algorithms for synchrophasor, frequency and ROCOF estimation, capable of rejecting spurious contributions as well as accounting for the dynamics of the fundamental component [10,11,12]. In particular, the traditional definition of phasor has been extended to the novel concept of dynamic synchrophasor, whose amplitude and phase angle may vary as functions of time [13,14].

In practice, though, an effective employment of PMU measurements requires complying with stringent latency limits [15]. Most of the existing approaches rely on segmenting the acquired signal and processing each segment via Discrete Fourier Transform (DFT) [16] or Taylor-Fourier Transform (TFT) [17], in order to extract the parameters associated to the fundamental component. In general, processing a longer segment enables better disturbance rejection and wider measurement bandwidth, but it also results in increased computational cost, higher latency and slower response. Therefore, an optimal tradeoff between accuracy and responsiveness is hard to find and a practical limit of two nominal cycles has been achieved through a combined application of Hilbert transform and TFT, at the expense though of higher computational complexity and sensitivity to spurious interferences [18]. A detailed review of the main synchrophasor estimation approaches is beyond the scope of this paper and can be found in [6,19].

The Fortescue transformation is often used for an effective analysis of three-phase AC networks: it allows expressing a sinusoidal three-phase signal in terms of its symmetrical components [20]. A peculiar characteristic of real-world power systems is that they are weakly unbalanced during regular operation: in this case, the positive sequence component alone contains most of the information associated to the fundamental term of a three-phase quantity. In the recent literature, several PMU algorithms exploits this symmetry property by using the Clarke/Park transformation [21] or Principal Component Analysis and Maximum Likelihood estimators [22]. In this context, the Space Vector (SV) -based approach [23] can be also adopted for reducing computational burden without sacrificing performance. Instead of dealing with each phase separately, the SV with respect to a reference frame rotating at the system nominal frequency is computed. Digital filters whose frequency responses can be customized to reach specific performance goals [24] are employed to extract the contribution of the positive sequence component from the complex-valued SV signal. It is worth noting that the SV approach can be also used in conjunction with other well-known methods for estimating synchrophasor, frequency and ROCOF such as interpolated DFT (IpDFT) [25] or Taylor-Fourier (TF) filtering [26], thus leveraging the three-phase symmetry of power systems for reduced computational cost and improved accuracy.

In this paper, we carry out a rigorous and extensive performance assessment among three SV-based estimation algorithms, specifically designed for high responsiveness applications, with maximum latency set to two nominal cycles. To the best of Authors’ knowledge, this review represents the first attempt of a comparative analysis for highlighting peculiarities and strengths of the different methods. For this purpose, we analyze the algorithms’ performance in terms of response time and measurement accuracy, both under steady-state and dynamic conditions; the P-class reference algorithm provided by the IEC Std is considered as benchmark. In this framework, the paper main contributions consist in: (i) a reasoned analysis of the most suitable setting for an optimal tradeoff between accuracy and latency, (ii) a thorough characterization of the algorithms’ performance not only in the entire IEC Std test set, but also under other conditions including unbalance and harmonic distortion at off-nominal frequency.

The paper is organized as follows. In Section 2, we introduce the considered measurement algorithms. Section 3 presents the results of the performance comparison. In Section 4, we provide some closing remarks and outline the future steps of the research activity.

2. Measurement Algorithms

In this paper, we present a comparison among three algorithms for estimating positive sequence synchrophasor, frequency and ROCOF starting from a three-phase signal. The considered techniques are the SV filtering algorithm (SV-F) [24], the SV Taylor-Fourier algorithm (SV-TF) [26] and the SV IpDFT algorithm (SV-IpDFT) [25]. All of them feature high flexibility: their parameters can be tuned in order to meet specific requirements in terms of accuracy and latency according to the application. In this respect, we focus on fast response, fixing the algorithm latency at about one and half nominal cycles. The P-class algorithm suggested by the IEC Std, which is used as reference to compare the performance obtained with the three aforementioned algorithms, has a different latency, but always below two cycles. In the following, the common ground and useful definitions are firstly introduced, then all the considered algorithms are explained in detail.

2.1. Common Background on Synchrophasor Measurements

Let us consider a three-phase power system characterized by the rated frequency $f_0$, corresponding to the angular frequency ${\omega }_0$. The synchrophasor approach is based on the following quasi-steady-state model of the generic phase p signal $x_p\left(t\right)$ (with $p\in \{a,b,c\}$), which may represent the voltage or current waveform measured in a node as a function of a shared time coordinate t:

$$\begin{array}{c}{x}_p\left(t\right)=x_{1,p}\left(t\right)+d_p\left(t\right)\\ x_{1,p}\left(t\right)=\sqrt{2}{X}_p\left(t\right)cos\left({\theta }_p\left(t\right)\right)=\sqrt{2}{X}_p\left(t\right)cos({\omega }_{0}t+{\phi }_p\left(t\right))\end{array}$$

(1)

Basically, the generic phase p waveform consists of a main term $x_{1,p}\left(t\right)$, which is assumed to be a modulated sinewave whose peak amplitude $\sqrt{2}{X}_p\left(t\right)$ and phase angle ${\phi }_p\left(t\right)$ are slowly-varying with respect to the rated frequency. $d_p\left(t\right)$ represents a second contribution containing the other components that may be present in the signal, which are considered as disturbances (e.g., harmonics and interharmonics); its magnitude is supposed to be considerably lower than that of $x_{1,p}\left(t\right)$, while its spectral content near $f_0$ is assumed to be negligible. Starting from these considerations, the corresponding phase p (dynamic) synchrophasor ${\overline{X}}_p\left(t\right)$ is defined as:

$${\overline{X}}_p\left(t\right)\triangleq X_p\left(t\right)e^{j{\phi }_p\left(t\right)}$$

(2)

In the following, overbar is used to highlight complex-valued quantities; for a lighter notation, a previously defined complex variable written without overbar denotes its magnitude. From an operative point of view, the synchrophasor is a possibly time-varying complex-valued number associated to an electric quantity of an AC power system that, under the previous assumptions, is able to carry its most relevant informative content. Using the Euler’s formula, an expression of $x_{1,p}\left(t\right)$ where the synchrophasor explicitly appears can be easily obtained:

$$x_{1,p}\left(t\right)=\frac{\sqrt{2}}{2}\left[{\overline{X}}_p\left(t\right)e^{j{\omega }_{0}t}+{\overline{X}}_p^{\ast }\left(t\right)e^{-j{\omega }_{0}t}\right]$$

(3)

where ${}^{\ast }$ denotes the complex conjugate operator. From another point of view, $x_{1,p}\left(t\right)$ has been decomposed into the sum of two counter-rotating vectors in the complex plane: one with positive, the other with negative angular speed, which is often known as image. It is worth noting that when assuming sinusoidal steady-state conditions at the rated frequency, the synchrophasor corresponds to the usual phasor. Considering once again the waveform model (1), the frequency $f\left(t\right)$ (corresponding to the angular frequency $\omega \left(t\right)$) and its rate of change ROCOF$\left(t\right)$ are defined as:

$$f\left(t\right)=\frac{1}{2\pi }\frac{d{\theta }_p\left(t\right)}{dt}=\frac{1}{2\pi }\frac{d{\phi }_p\left(t\right)}{dt}+f_0=\Delta f\left(t\right)+f_0$$

(4)

$$\mathrm{ROCOF}\left(t\right)=\frac{df\left(t\right)}{dt}=\frac{1}{2\pi }\frac{d^2{\theta }_p\left(t\right)}{d{t}^2}=\frac{1}{2\pi }\frac{d^2{\phi }_p\left(t\right)}{d{t}^2}$$

(5)

$\Delta f\left(t\right)$ is also introduced: it represents the frequency deviation with respect to rated conditions, which is the rotational speed of the synchrophasor in the complex plane divided by $2\pi$. The target of a PMU is sampling the three-phase waveform of an electrical quantity and, by means of a proper algorithm, extracting the corresponding synchrophasors, frequency and ROCOF at discrete and fixed time instants $t_r=i{T}_{RR}$ (with i integer, see [7]), hence multiples of the reporting interval $T_{RR}$. Since $T_{RR}$ is generally an integer multiple of $T_0=1/f_0$, ${\theta }_p\left(t_r\right)={\phi }_p\left(t_r\right)$ when phase angle wrapping is considered. The algorithm should be able to track dynamic variations as well as rejecting disturbances.

The study of three-phase systems is more efficient when it is performed by using the symmetrical components. Having available the phase a, b and c synchrophasors, the corresponding positive, negative and zero sequence synchrophasors (${\overline{X}}_{+}\left(t\right)$, ${\overline{X}}_{-}\left(t\right)$ and ${\overline{X}}_0\left(t\right)$, respectively) can be obtained as follows by means of the Fortescue transformation (here the unitary formulation is considered):

$$\left[\begin{array}{c}{\overline{X}}_{+}\left(t\right)\\ {\overline{X}}_{-}\left(t\right)\\ {\overline{X}}_0\left(t\right)\end{array}\right]=\left[\begin{array}{c}{X}_{+}\left(t\right)e^{j{\phi }_{+}\left(t\right)}\\ X_{-}\left(t\right)e^{j{\phi }_{-}\left(t\right)}\\ X_0\left(t\right)e^{j{\phi }_0\left(t\right)}\end{array}\right]=\frac{1}{\sqrt{3}}\left[\begin{array}{ccc}1& \overline{\alpha }& {\overline{\alpha }}^2\\ 1& {\overline{\alpha }}^2& \overline{\alpha }\\ 1& 1& 1\end{array}\right]\left[\begin{array}{c}{\overline{X}}_a\left(t\right)\\ {\overline{X}}_b\left(t\right)\\ {\overline{X}}_c\left(t\right)\end{array}\right]$$

(6)

where $\overline{\alpha }\triangleq e^{j2\pi /3}$. Among the symmetrical components, the positive sequence term has by far the highest magnitude, since three-phase systems are weakly unbalanced during regular operation, in particular as far as the transmission grid is considered; for this reason, many applications rely on a positive sequence representation of the power system.

It is worth noting that processing the single-phase signals while using the definition (4) may lead to different frequency measurements for each of the phases; conversely, as from [7], a unique frequency (and ROCOF) value should be provided for each three-phase quantity. On the one hand, this can be obtained through additional processing (e.g., averaging) on the three frequency estimates [27]; on the other hand, computing a shared frequency value reminds that, from a physical point of view, a three-phase quantity is not a mere set of three independent single-phase quantities. In this regard, a combined processing of the three signals would enable a direct estimate of frequency and ROCOF other than the symmetrical components of the synchrophasors, and in particular of the positive sequence term. Moreover, native three-phase PMU algorithms may exploit that, under regular operation, three-phase quantities are weakly unbalanced: this may be beneficial to improve the quality of the estimates while reducing the computational cost.

2.2. IEC/IEEE 60255-118-1 P-Class Reference Algorithm

The IEC Std suggests a measurement model for synchrophasor estimation and two PMU algorithms, one for class P, intended for applications requiring fast response, and another one for class M, designed for measurement applications. Both of them are based on the same concept of quadrature demodulation, which is briefly summarized in the following. Let us substitute (3) in the decomposition (1) of the phase p waveform $x_p\left(t\right)$ and let us multiply by $e^{-j{\omega }_{0}t}$, namely a unitary magnitude vector whose rotational speed in the complex plane corresponds to the rated angular frequency of the power system. The signal ${\overline{x}}_{p,d}\left(t\right)$ is thus obtained, which results:

$${\overline{x}}_{p,d}\left(t\right)=\frac{\sqrt{2}}{2}{\overline{X}}_p\left(t\right)+\frac{\sqrt{2}}{2}{\overline{X}}_p^{\ast }\left(t\right)e^{-j2{\omega }_{0}t}+d_p\left(t\right)e^{-j{\omega }_{0}t}$$

(7)

From the previous expression, ${\overline{x}}_{p,d}\left(t\right)$ is a complex-valued signal consisting of three contributions. Except for a constant scale factor equal to $\sqrt{2}/2$, the first term coincides with the quantity of interest, i.e., the synchrophasor, that—based on the assumptions introduced in Section 2.1—is expected to have a spectral content restricted at very low frequencies. The second term is generated by the counter-rotating image component in (3), and is represented by the same scale factor multiplied by the synchrophasor conjugate and an unitary vector, that is rotating clockwise in the complex plane with angular speed equal to $2{\omega }_0$; therefore, its frequency content is confined around $-2f_0$. The last term is the disturbance multiplied by an unitary vector rotating in the complex plane with angular speed equal to ${\omega }_0$; reminding the hypothesis about $d_p\left(t\right)$, the spectral content of this third contribution is virtually entirely located away from zero. Thanks to the spectral separation of these terms, the synchrophasor can be extracted from ${\overline{x}}_{p,d}\left(t\right)$ through proper low-pass filtering. In this respect, it should be stressed that it is particularly important to somehow attenuate the second contribution (namely that involving the conjugate of the synchrophasor): it represents a very large unwanted component, since it has the same magnitude as the low-frequency term to be extracted.

After having explained the basic concepts, we focus on the details of the P-class algorithm reported in [7] (from here on IEC-P algorithm for the sake of brevity), which will be considered as a benchmark for the comparison among the different low-latency techniques; its block diagram is reported in Figure 1.

Let us suppose that the phase quantities have been sampled with rate $f_s=M_C{f}_0$ (the positive integer $M_C$ is the number of samples per nominal cycle) corresponding to the sampling interval $T_s=1/f_s$; it is also assumed that $T_{RR}=Q{T}_s$, with Q positive integer. For each phase p, a first estimate ${\hat{\overline{X}}}_p^{\prime }$ of the synchrophasor in the generic reporting instant $t_r=i{T}_{RR}$ is obtained by means of the following expression, which corresponds to multiplication by the rotating exponential and finite impulse response (FIR) filtering:

$${\hat{\overline{X}}}_p^{\prime }\left(t_r\right)={\hat{X}}_p^{\prime }\left(t_r\right)e^{j{\hat{\phi }}_p\left(t_r\right)}=\frac{\sqrt{2}}{G}\left({\overline{x}}_{p,d}\ast w_{tri}\right)\left(t_r\right)=\frac{\sqrt{2}}{G}\sum _{n=-M_C+1}^{M_C-1}{x}_p(t_r-n{T}_s)w_{tri}\left[n\right]e^{-j{\omega }_0(t_r-n{T}_s)}$$

(8)

where * indicates the convolution sum; from here on, hat is used to denote estimated values. Filter coefficients $w_{tri}\left[n\right]$ correspond to the $2M_C-1$ sample triangular window:

$$w_{tri}\left[n\right]=1-\frac{\left|n\right|}{M_C},\mathrm{for}n=-M_C+1,\dots ,M_C-1$$

(9)

and the constant G is its DC gain:

$$G=\sum _{n=-M_C+1}^{M_C-1}{w}_{tri}\left[n\right]=M_C$$

(10)

It is worth noting that the low-pass filter defined by (9) is linear-phase, non causal, with zero group delay. Let us assume sinusoidal steady-state operation with frequency equal to $f_0$: ${\overline{X}}_p\left(t\right)$ in (7) is constant, while the term involving ${\overline{X}}_p^{\ast }\left(t\right)$ rotates with angular speed $-2{\omega }_0$. It can be easily proven that the employed filter has frequency response zeros located at multiples of $f_0$. Hence, under these conditions it is able to completely cancel out the term rotating at $-2{\omega }_0$. Since the DC gain of the triangular window FIR filter defined as in (9) is $G=M_C$, the estimate provided by (8) is theoretically exact.

In this algorithm, a three-phase definition of frequency is employed in order to extract a unique value from the three phases. In particular, the frequency deviation $\Delta f$ corresponds to the rotational speed (in hertz) of the positive sequence synchrophasor, whose estimate ${\hat{\overline{X}}}_{+}^{\prime }$ is obtained by applying the unitary Fortescue transformation (6):

$${\hat{\overline{X}}}_{+}^{\prime }\left(t_r\right)={\hat{X}}_{+}^{\prime }\left(t_r\right)e^{j{\hat{\phi }}_{+}\left(t_r\right)}=\frac{{\hat{\overline{X}}}_a^{\prime }\left(t_r\right)+\overline{\alpha }{\hat{\overline{X}}}_b^{\prime }\left(t_r\right)+{\overline{\alpha }}^2{\hat{\overline{X}}}_c^{\prime }\left(t_r\right)}{\sqrt{3}}$$

(11)

From a physical point of view, this definition of frequency corresponds to the electrical speed of the air gap magnetic field in an ideal, symmetric three-phase machine whose winding magnetomotive force is purely positive sequence and expressed by the synchrophasor ${\overline{X}}_{+}$.

Frequency deviation and ROCOF measurements are obtained as first and second order discrete-time derivatives of the estimated phase angle ${\hat{\phi }}_{+}$. For this purpose, let us firstly compute ${\hat{\phi }}_{+}$ in the time instants $t_r-S{T}_s$ and $t_r+S{T}_s$, where $S\le Q$ is a positive integer; after that, frequency and ROCOF estimates result:

$$\hat{f}\left(t_r\right)=\widehat{\Delta f}\left(t_r\right)+f_0=\frac{{\hat{\phi }}_{+}(t_r+S{T}_s)-{\hat{\phi }}_{+}(t_r-S{T}_s)}{4\pi S{T}_s}+f_0$$

(12)

$$\widehat{\mathrm{ROCOF}}\left(t_r\right)=\frac{{\hat{\phi }}_{+}(t_r+S{T}_s)+{\hat{\phi }}_{+}(t_r-S{T}_s)-2{\hat{\phi }}_{+}\left(t_r\right)}{2\pi {\left(S{T}_s\right)}^2}$$

(13)

Therefore, (12) and (13) correspond to the application of two linear-phase FIR filters $D_S^1$ and $D_S^2$ to ${\hat{\phi }}_{+}/\left(2\pi \right)$, respectively. Now, let us suppose that the input signals $x_p\left(t\right)$ are pure sinewaves whose frequency f is different from $f_0$. In this case, the phase p synchrophasor is:

$${\overline{X}}_p\left(t\right)=X_p{e}^{j((\omega -{\omega }_0)t+{\phi }_{ph,p})}={\overline{X}}_{ph,p}{e}^{j(\omega -{\omega }_0)t}$$

(14)

where ${\overline{X}}_{ph,p}$ represents the usual phasor at the frequency f. Hence, the synchrophasor is a complex number having the same magnitude as the phasor, but its phase angle is time-varying, since it rotates in the complex plane with angular speed equal to the difference between $\omega$ and ${\omega }_0$. In particular, the real and imaginary parts of ${\overline{X}}_p\left(t\right)$ are quadrature sinewaves whose frequencies are equal to the frequency deviation $\Delta f$. From (7), the signal to be filtered is made of two contributions:

$${\overline{x}}_{p,d}\left(t\right)=\frac{\sqrt{2}}{2}{\overline{X}}_{ph,p}{e}^{j(\omega -{\omega }_0)t}+\frac{\sqrt{2}}{2}{\overline{X}}_{ph,p}^{\ast }{e}^{-j(\omega +{\omega }_0)t}$$

(15)

Introducing $W_{tri}\left(j\omega \right)$ as the frequency response of the previously defined triangular window filter (real-valued and identical to the amplitude response), the evaluated synchrophasor results:

$${\hat{\overline{X}}}_p^{\prime }\left(t_r\right)=\frac{W_{tri}(j2\pi \Delta f)}{G}{\overline{X}}_p\left(t_r\right)+\frac{W_{tri}\left(j(\omega +{\omega }_0)\right)}{G}{\overline{X}}_{ph,p}^{\ast }{e}^{-j(\omega +{\omega }_0)t_r}$$

(16)

reminding that $2\pi \Delta f=\omega -{\omega }_0$. The first thing to be noticed is that the term proportional to ${\overline{X}}_{ph,p}^{\ast }$ is attenuated but not completely rejected by the filter, since its rotational speed does not correspond to a zero of the filter. Furthermore, even neglecting this effect, the magnitude of the synchrophasor estimate is biased: in fact $W_{tri}(j2\pi \Delta f)\ne G$ when $\Delta f\ne 0$. This effect is significant even for small values of the frequency deviation. However, having measured $\widehat{\Delta f}$, it can be compensated since the filter response is known; a more accurate synchrophasor estimate ${\hat{\overline{X}}}_p$ can thus be obtained. In this respect, ref. [7] adopts the following approximated expression instead of using the exact filter response:

$${\hat{\overline{X}}}_p\left(t_r\right)=\frac{{\hat{\overline{X}}}_p^{\prime }\left(t_r\right)}{A\left(\widehat{\Delta f}\left(t_r\right)\right)}=\frac{{\hat{\overline{X}}}_p^{\prime }\left(t_r\right)}{sin\left(\pi \frac{f_0+1.62\widehat{\Delta f}\left(t_r\right)}{2f_0}\right)}$$

(17)

As mentioned before, the positive sequence synchrophasor estimate has key importance in power systems. Let us obtain its expression under off-nominal frequency conditions; the substitution of (16) into (11) leads to:

$${\hat{\overline{X}}}_{+}\left(t_r\right)=\frac{W_{tri}(j2\pi \Delta f)}{A\left(\widehat{\Delta f}\left(t_r\right)\right)}{\overline{X}}_{+}\left(t_r\right)+\frac{W_{tri}\left(j(\omega +{\omega }_0)\right)}{A\left(\widehat{\Delta f}\left(t_r\right)\right)}{\overline{X}}_{-}^{\ast }\left(t_r\right)e^{-j2{\omega }_0{t}_r}$$

(18)

As in (16), also ${\hat{\overline{X}}}_{+}\left(t_r\right)$ is disturbed by the presence of a component which rotates clockwise in the complex plane with angular speed $\omega +{\omega }_0$. Based on the assumption of weakly unbalanced system, such disturbance has a much lower relative magnitude, since $W_{tri}(j2\pi \Delta f$) is close to G, whereas $X_{-}/X_{+}$ (often called unbalance factor) is well below one. Therefore, the positive sequence synchrophasor estimate is considerably less affected by this phenomenon. Finally, it is worth stressing that ${\hat{\overline{X}}}_{+}$ does not depend on the zero sequence component.

The overall algorithm latency is given by $L=(M_C-1)T_s+S{T}_s$, which keeps into account the filtering process and the phase-angle derivative computation. It is thus always $T_0\le L<2T_0$, thus $S\le M_C$ in order to comply with the latency limit for P class algorithms.

2.3. Space Vector Algorithm

PMU algorithms based on the SV transformation of the phase quantities and filtering of the resulting signal when decomposed in rectangular and polar coordinates have been recently proposed by the authors in [23,24]. The starting point is the usual measurement model (1), here rewritten by using the vector notation in order to consider the three phases simultaneously:

$${\mathbf{x}}_{abc}\left(t\right)={\mathbf{x}}_{1,abc}\left(t\right)+{\mathbf{d}}_{abc}\left(t\right)$$

(19)

where:

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{1,abc}\left(t\right)=\left[\begin{array}{c}{x}_{1,a}\left(t\right)\\ x_{1,b}\left(t\right)\\ x_{1,c}\left(t\right)\end{array}\right]& \hfill & \hfill {\mathbf{d}}_{abc}\left(t\right)=\left[\begin{array}{c}{d}_a\left(t\right)\\ d_b\left(t\right)\\ d_c\left(t\right)\end{array}\right]\end{array}$$

(20)

From the three-phase waveforms, it is possible to compute the SV ${\overline{x}}_{SV}\left(t\right)$ in a rotating reference frame whose instantaneous angular position is $\beta \left(t\right)$; it results:

$${\overline{x}}_{SV}\left(t\right)=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}1& \overline{\alpha }& {\overline{\alpha }}^2\end{array}\right]{\mathbf{x}}_{abc}\left(t\right)e^{-j\beta \left(t\right)}={\overline{x}}_{S{V}_0}\left(t\right)e^{-j\beta \left(t\right)}$$

(21)

where ${\overline{x}}_{S{V}_0}$ is the SV in a stationary reference frame characterized by $\beta =0$. Substituting (19) into (21) while using (3) and reminding (6) leads to:

$$\begin{array}{c}{\overline{x}}_{SV}\left(t\right)={\overline{X}}_{+}\left(t\right)e^{j({\omega }_{0}t-\beta \left(t\right))}+{\overline{X}}_{-}^{\ast }\left(t\right)e^{-j({\omega }_{0}t+\beta \left(t\right))}+{\overline{d}}_{SV}\left(t\right)\\ {\overline{d}}_{SV}\left(t\right)=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}1& \overline{\alpha }& {\overline{\alpha }}^2\end{array}\right]{\mathbf{d}}_{abc}\left(t\right)e^{-j\beta \left(t\right)}={\overline{d}}_{S{V}_0}\left(t\right)e^{-j\beta \left(t\right)}\end{array}$$

(22)

Three terms appear: the first one is related to the positive sequence synchrophasor, the second depends on the negative sequence synchrophasor, the last is produced by disturbances. Of course, the term due to the positive sequence synchrophasor is expected to be considerably higher with respect to the others. Furthermore, it is worth noting that the zero sequence synchrophasor does not affect the SV signal. Now, let us choose $\beta \left(t\right)={\omega }_{0}t$, namely performing the SV transformation in a reference frame whose angular speed corresponds to the rated angular frequency. In this case, the positive sequence synchrophasor produces a very low frequency contribution; instead, the term related to the negative sequence synchrophasor rotates with an angular speed close to $2{\omega }_0$ in clockwise direction. Reminding that the spectral content of $d_p\left(t\right)$ is assumed to be located away from $f_0$, the spectrum of ${\overline{d}}_{SV}\left(t\right)$ is well separated from that of the positive sequence synchrophasor. This suggests that an estimate of ${\overline{X}}_{+}\left(t\right)$ can be extracted from ${\overline{x}}_{SV}\left(t\right)$ by means of proper low-pass filtering, without explicitly computing per phase synchrophasors. Adopting the same three-phase definition of frequency deviation already employed by the IEC-P algorithm, $\Delta f$ and ROCOF represent (except for a scale factor $2\pi$) the first and second order derivatives of the positive sequence synchrophasor phase angle. Therefore, their estimates can be obtained by filtering ${\overline{x}}_{SV}$, computing its phase angle and performing numerical differentiations.

The architecture of this approach corresponds to the block diagram reported in Figure 2; its practical implementation requires first of all to properly sample the phase signals with rate $f_s$. From the acquired data, applying (21) with $\beta \left(t\right)={\omega }_{0}n{T}_s$ the samples of the SV signal ${\overline{x}}_{SV}$ are obtained. A first low-pass filtering stage H with unit DC gain is applied to the real and imaginary parts of ${\overline{x}}_{SV}$ in order to remove most of ${\overline{d}}_{SV}$ as well as the term produced by the negative sequence. Both FIR and infinite impulse response (IIR) filters can be employed, even if best performance and greater flexibility is obtained with FIR design. The samples of ${\overline{x}}_{SV,f}$ are thus obtained and they have to be further processed in order to compensate for the group and phase delay introduced by H; it is extremely simple when a linear-phase FIR filter is adopted.

After that, ${\overline{x}}_{SV,f}$ is decomposed into its magnitude $x_{SV,f}$ and argument ${\phi }_{SV,f}$. The estimate ${\hat{\phi }}_{+}\left(t_r\right)$ of the positive sequence synchrophasor phase angle in the generic reporting instant $t_r$ is obtained by applying to ${\phi }_{SV,f}$ a linear-phase FIR low-pass filter P (characterized by unit DC gain, designed to attenuate the residual impact of ${\overline{d}}_{SV}$ and of the oscillating term produced by ${\overline{X}}_{-}$) and compensating its delay. Its bandwidth should be sufficient to properly follow the phase angle dynamics, which are typically rather slow.

Similarly, frequency deviation and ROCOF measurements in the reporting instants ($\widehat{\Delta f}\left(t_r\right)$ and $\widehat{\mathrm{ROCOF}}\left(t_r\right)$, $t_r=Q{T}_s$) are computed by filtering ${\phi }_{SV,f}$ with linear-phase, partial band first and second-order FIR differentiators (F and R respectively) and compensating the introduced delays. Their frequency responses should match those of ideal first and second order differentiators only near zero (just slow frequency variations have to be accurately tracked), while they should have considerably lower magnitudes at higher frequencies in order to achieve good disturbance rejection and noise immunity.

An estimate ${\hat{X}}_{+}^{\prime }$ of the positive sequence synchrophasor magnitude can be obtained by applying a low-pass filter M to $x_{SV,f}$ while taking into account the introduced delay. Similarly to P, M should have unitary DC gain and should be designed to reject the infiltration of ${\overline{d}}_{SV}$ as well as the term proportional to ${\overline{X}}_{-}$. Its bandwidth should be large enough to properly follow amplitude modulations, which are fairly slow in practical power systems.

However, it should be noticed that ${\hat{X}}_{+}^{\prime }$ may be significantly biased under off-nominal frequency conditions, similarly to what happens with the IEC-P algorithm discussed in the previous subsection. Assuming perfectly sinusoidal, positive sequence input with angular frequency $\omega \ne {\omega }_0$, ${\overline{x}}_{SV,f}$ results:

$${\overline{x}}_{SV,f}\left(n{T}_s\right)=H(j2\pi \Delta f){\overline{X}}_{+}\left(n{T}_s\right)=H(j2\pi \Delta f){\overline{X}}_{ph,+}{e}^{j(\omega -{\omega }_0)n{T}_s}$$

(23)

where ${\overline{X}}_{ph,+}$ is the positive sequence phasor while filter H is assumed to be a linear-phase FIR filter having amplitude response $H\left(j\omega \right)$ For a linear-phase filter H, the relationship between its frequency response $\overline{H}\left(j\omega \right)$ and amplitude response $H\left(j\omega \right)$ is: $\overline{H}\left(j\omega \right)=H\left(j\omega \right)e^{-j\omega {\tau }_H}$, where ${\tau }_H$ is the group delay. According to the procedure described in the previous lines, it is easy to derive the magnitude estimate ${\hat{X}}_{+}=H(j2\pi \Delta f)X_{+}$, reminding that the DC gain of M is unitary by assumption. Similarly to what explained in Section 2.2, the effect of filter H can be removed since its response is known and a measurement of the frequency deviation is available. Hence a better estimate ${\hat{X}}_{+}$ of the positive sequence synchrophasor magnitude is:

$${\hat{X}}_{+}\left(t_r\right)=\frac{{\hat{X}}_{+}^{\prime }\left(t_r\right)}{H\left(j2\pi \widehat{\Delta f}\left(t_r\right)\right)}$$

(24)

One of the main advantages of this SV-based approach is its flexibility: achieved performance strongly depends on the filters H, M, P, F and R, which can be tuned in order to reach predetermined goals in terms of accuracy, latency and responsiveness. It is worth highlighting that closed-form expressions provided in [24] allow predicting the results of the P- and M-class compliance tests prescribed by [7], thus substantially helping the design of the aforementioned filters.

The algorithm latency in this case is $L=((N_H+max\{N_M,N_P,N_F,N_R\})/2-1)T_s$, where $N_H$, $N_M$, $N_P$, $N_F$, $N_R$ are the number of taps of the filters H, M, P, F and R, respectively.

2.4. Space Vector Taylor-Fourier Algorithm

In [13], the TF approach to synchrophasor estimation was introduced. Basically, it is based on a measurement model which is a truncated Taylor expansion of the phase p synchrophasor around the generic reporting time $t_r$ in order to model its time evolution. Model parameters are obtained through least-squares fitting of a sliding window of the collected samples, thus corresponding to FIR filtering. Synchrophasor, frequency and ROCOF estimates can be easily obtained from the model parameters. It is worth noting that a frequency estimate for each phase is obtained. Since according to [7] a unique frequency value for each three-phase quantity must be provided, a possible solution is computing the average between the three different frequency measurements.

In [26], the technique has been extended for being applied to the SV signal, thus directly estimating positive and negative sequence synchrophasors. This enables higher design flexibility with respect to a per phase approach, which results in better performance and lower computational burden. Furthermore, it directly provides frequency and ROCOF values according to the three-phase definition adopted by the previously presented algorithms.

The block diagram of the method is reported in Figure 3; the starting point is the SV ${\overline{x}}_{S{V}_0}\left(t\right)$ on a stationary reference frame, which is obtained by applying the SV transformation (21) with $\beta =0$ to the three-phase signals. Under the usual assumptions, it results:

$${\overline{x}}_{S{V}_0}\left(t\right)={\overline{X}}_{+}\left(t\right)e^{j{\omega }_{0}t}+{\overline{X}}_{-}^{\ast }\left(t\right)e^{-j{\omega }_{0}t}+{\overline{d}}_{S{V}_0}\left(t\right)$$

(25)

As mentioned above, the TF approach relies on a measurement model obtained through truncated Taylor expansions of the synchrophasors around the reporting instant. Therefore, considering the positive and negative sequence synchrophasors in the neighborhood of $t_r$:

$$\begin{array}{c}{\overline{X}}_{+}\left(t\right)\approx \sum _{k=0}^{K_{+}}{\overline{X}}_{+}^{\left(k\right)}\left(t_r\right)\frac{{(t-t_r)}^k}{k!}\\ {\overline{X}}_{-}\left(t\right)\approx \sum _{k=0}^{K_{-}}{\overline{X}}_{-}^{\left(k\right)}\left(t_r\right)\frac{{(t-t_r)}^k}{k!}\end{array}$$

(26)

where ${\overline{X}}_{+}^{\left(k\right)}\left(t_r\right)$ and ${\overline{X}}_{-}^{\left(k\right)}\left(t_r\right)$ are the kth order derivatives of the positive and negative sequence synchrophasor model at the reporting instant $t_r$; considering $k=0$, they correspond to the two synchrophasors. Expansion orders $K_{+}$ and $K_{-}$ are in general different: this additional degree of freedom is enabled thanks to the SV-based approach, while it would have not been possible if the conventional implementation were adopted.

The expressions (26) can be used to write a model for the SV signal in the neighborhood of the reporting instant:

$${\overline{x}}_{S{V}_0}\left(t\right)\approx {\overline{X}}_{+}\left(t\right)e^{j{\omega }_{0}t}+{\overline{X}}_{-}^{\ast }\left(t\right)e^{-j{\omega }_{0}t}$$

(27)

Now, let us suppose that the phase signals have been sampled with rate $f_s$; the samples of the SV signal ${\overline{x}}_{S{V}_0}$ can be easily computed. Furthermore, let us consider a sliding window made of $N_w$ samples (with $N_w$ odd) of ${\overline{x}}_{S{V}_0}$ and centered on $t_r$; these samples can be arranged in a vector ${\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)$. Using (27) it is possible to write the corresponding vector of $N_w$ samples obtained from the model, having assumed that its parameters (namely the synchrophasor derivatives) are constant within the window. Adopting vector notation, this leads to:

$${\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)\approx \left[\overline{\mathbf{\Phi }}{\mathbf{A}}_{+}{\overline{\mathbf{\Phi }}}^{\mathrm{H}}{\mathbf{A}}_{-}\right]·{\overline{\mathbf{p}}}_{\pm }\left(t_r\right)={\overline{\mathbf{B}}}_{\pm }·{\overline{\mathbf{p}}}_{\pm }\left(t_r\right)$$

(28)

where:

$$\begin{array}{cc}\hfill & {\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)=\left[\begin{array}{c}{\overline{x}}_{S{V}_0}(t_r+\frac{N_w-1}{2}{T}_s)\\ ⋮\\ {\overline{x}}_{S{V}_0}\left(t_r\right)\\ ⋮\\ {\overline{x}}_{S{V}_0}(t_r-\frac{N_w-1}{2}{T}_s))\end{array}\right]\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\overline{\mathbf{p}}}_{\pm }\left(t_r\right)={\left[{\overline{X}}_{+}^{\left(0\right)}\left(t_r\right),\cdots ,{\overline{X}}_{+}^{\left(K_{+}\right)}\left(t_r\right),{\overline{X}}_{-}^{\left(0\right)\ast }\left(t_r\right),\cdots ,{\overline{X}}_{-}^{\left(K_{-}\right)\ast }\left(t_r\right)\right]}^{\mathrm{T}}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & {\mathbf{A}}_{+}=\left[\begin{array}{ccccc}1& \frac{N_w-1}{2}{T}_s& \frac{{\left(\frac{N_w-1}{2}{T}_s\right)}^2}{2}& \cdots & \frac{{\left(\frac{N_w-1}{2}{T}_s\right)}^{K_{+}}}{k!}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& 0& 0& & 0\\ ⋮& ⋮& ⋮& & ⋮\\ 1& -\frac{N_w-1}{2}{T}_s& \frac{{\left(-\frac{N_w-1}{2}{T}_s\right)}^2}{2}& \cdots & \frac{{\left(-\frac{N_w-1}{2}{T}_s\right)}^{K_{+}}}{k!}\end{array}\right]\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & \overline{\mathbf{\Phi }}=\left[\begin{array}{ccccc}{e}^{j{\omega }_0\frac{N_w-1}{2}{T}_s}& & & & \\ & \ddots & & & \\ & & 1& & \\ & & & \ddots & \\ & & & & e^{-j{\omega }_0\frac{N_w-1}{2}{T}_s}\end{array}\right]e^{j{\omega }_0{t}_r}\hfill \end{array}$$

(32)

and ${\mathbf{A}}_{-}$ is obtained from (31) by replacing $K_{+}$ with $K_{-}$; the superscripts ${}^{\mathrm{H}}$ and ${}^{\mathrm{T}}$ indicate the Hermitian transpose and the transpose operators, respectively.

${\overline{\mathbf{p}}}_{\pm }\left(t_r\right)$ is the vector of the model parameters in the reporting instant. It can be estimated by minimizing the Euclidean norm of the vector of the differences between the samples of the SV signal and those obtained from the model, hence:

$${\widehat{\overline{\mathbf{p}}}}_{\pm }\left(t_r\right)=\underset{{\overline{\mathbf{p}}}_{\pm }\left(t_r\right)}{arg\; min}∥{\overline{\mathbf{B}}}_{\pm }·{\overline{\mathbf{p}}}_{\pm }\left(t_r\right)-{\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)∥$$

(33)

which corresponds to an ordinary least squares (LS) problem whose solution is:

$$\begin{array}{c}{\widehat{\overline{\mathbf{p}}}}_{\pm }\left(t_r\right)=\overline{\mathbf{H}}·{\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)\\ \overline{\mathbf{H}}={\left({\overline{\mathbf{B}}}_{\pm }^{\mathrm{H}}{\overline{\mathbf{B}}}_{\pm }\right)}^{-1}{\overline{\mathbf{B}}}_{\pm }^{\mathrm{H}}\end{array}$$

(34)

where $\overline{\mathbf{H}}$ is a complex filter bank, namely its hth row ${\overline{\mathbf{H}}}_{h,\ast }$ contains the complex-valued coefficients of the FIR filter that permits obtaining the hth element of ${\widehat{\overline{\mathbf{p}}}}_{\pm }$. In the same fashion, $\overline{\mathbf{H}}\left(j\omega \right)$ denotes the vector-valued frequency response of the filter bank; its hth element is the frequency response of the FIR filter defined by the coefficients ${\overline{\mathbf{H}}}_{h,\ast }$. Therefore, the estimated positive sequence synchrophasor is:

$${\hat{\overline{X}}}_{+}\left(t_r\right)={\overline{\mathbf{H}}}_{1,\ast }·{\overline{\mathbf{x}}}_{S{V}_0}\left(t_r\right)$$

(35)

Estimates can also be obtained by using a weighted LS (WLS) method [28] where the weights are given, for instance, by the squared coefficients of a cosine window. Up to now, the WLS estimator has never been applied in conjunction with the SV-TF approach.

As previously mentioned, the TF expansion applied to the complex-valued SV signal permits selecting different expansion orders for the positive and negative sequence synchrophasor, which is not possible when considering the conventional implementation in the real-valued phase signals: this additional degree of freedom in the design of TF filters enables better performance. In general, increasing the positive sequence synchrophasor expansion order $K_{+}$ is beneficial, since the model is able to better represent its dynamics. On the contrary, when the focus is measuring the positive sequence synchrophasor, choosing $K_{-}=K_{+}$ often leads to an overparametrized model. In this case, the best choice is $K_{-}<K_{+}$, since it increases the robustness with respect to noise and disturbances (which are bettered filtered thanks to the “stiffer” underlying model), without significant drawbacks.

According to the three-phase definition of frequency deviation and ROCOF, they can be obtained, with the formulas reported in [13], from the estimated derivatives of the positive sequence synchrophasor at $t_r$, that is as:

$$\begin{array}{cc}\hfill \widehat{\Delta f}\left(t_r\right)& =\frac{1}{2\pi }\frac{\Im \left[{\hat{\overline{X}}}_{+}^{\left(1\right)}\left(t_r\right)·{\hat{\overline{X}}}_{+}^{\ast }\left(t_r\right)\right]}{{\left|{\hat{\overline{X}}}_{+}\left(t_r\right)\right|}^2}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \widehat{\mathrm{ROCOF}}\left(t_r\right)& =\frac{1}{\pi }\left(\frac{\Im \left[{\hat{\overline{X}}}_{+}^{\left(2\right)}\left(t_r\right){\hat{\overline{X}}}_{+}^{\ast }\left(t_r\right)\right]}{2{\left|{\hat{\overline{X}}}_{+}\left(t_r\right)\right|}^2}-\frac{\Re \left[{\hat{\overline{X}}}_{+}^{\left(1\right)}\left(t_r\right){\hat{\overline{X}}}_{+}^{\ast }\left(t_r\right)\right]\Im \left[{\hat{\overline{X}}}_{+}^{\left(1\right)}\left(t_r\right){\hat{\overline{X}}}_{+}^{\ast }\left(t_r\right)\right]}{{\left|{\hat{\overline{X}}}_{+}\left(t_r\right)\right|}^4}\right)\hfill \end{array}$$

(37)

From the above equations, it is clear that an order $K_{+}\ge 2$ is needed to estimate frequency and ROCOF.

The latency of the algorithm depends only on the length of the window and thus, in this case, is $L=\frac{N_w-1}{2}{T}_s$.

2.5. Space Vector IpDFT Algorithm

IpDFT algorithms based on proper weighting windows [29,30] have been widely used to measure spectral components and frequencies of electrical signals [31]. More recently, several PMU techniques exploiting this approach have been proposed in the literature [16,32,33]. In order to obtain a synchrophasor estimate, they require observing the phase signal $x_p\left(t\right)$ over a time interval centered around the reporting instant $t_r$ and corresponding to an integer number $C\ge 2$ of rated cycles; the sampling frequency $f_s$ is multiple of the rated frequency $f_0$, so that $M_C$ samples per nominal cycle are collected. The underlying signal model is very similar to (1), but in this case the amplitude and frequency (defined as in (4)) are assumed to be constant within the analyzed time interval. Therefore, it results:

$$\begin{array}{cc}\hfill x_p\left(t\right)& =\frac{\sqrt{2}}{2}\left[{\overline{X}}_{ph,p}\left(t_r\right)e^{-j\omega \left(t_r\right)t_r}{e}^{j(\omega \left(t_r\right)-{\omega }_0)t}{e}^{j{\omega }_{0}t}+{\overline{X}}_{ph,p}^{\ast }\left(t_r\right)e^{j\omega \left(t_r\right)t_r}{e}^{-j(\omega \left(t_r\right)-{\omega }_0)t}{e}^{-j{\omega }_{0}t}\right]+d_p\left(t\right)\hfill \\ \hfill & =\frac{\sqrt{2}}{2}\left[{\overline{X}}_p\left(t_r\right)e^{j{\omega }_{0}t}+{\overline{X}}_p^{\ast }\left(t_r\right)e^{-j{\omega }_{0}t}\right]+d_p\left(t\right)\hfill \end{array}$$

(38)

The first step of the IpDFT algorithm is applying a proper tapering window to the collected samples. Assuming that spectral interference is negligible, there should be at least two DFT bins produced by the signal component rotating with angular frequency $\omega$ (namely the one related to the synchrophasor to be estimated) which falls under the main lobe of the window. Since its shape is known, the ratio between the magnitudes of these components allows estimating the frequency deviation $\Delta f$ and, in turns, also the synchrophasor. In particular, the knowledge of $\Delta f$ allows compensating the effect of scalloping loss which may undermine the evaluated magnitude under off-nominal frequency conditions.

As mentioned above, the method properly works if the two considered DFT bins are not affected by spectral interference: this may be produced by the disturbance $d_p\left(t\right)$, but also by the image component (namely the counter-rotating term appearing in (38)): specific care has to be taken, since it has the same magnitude as the term to be evaluated. This effect can be reduced by selecting a proper window (e.g., Rife-Vincent type I, characterized by maximum asymptotic decay of the sidelobes) and increasing the observation interval, even if the latter expedient results in several drawbacks. First of all, algorithm latency and computational burden are increased; secondly, the assumption of having constant signal parameters over a longer time window becomes harder to be met, thus it may seriously jeopardize the achieved dynamic performance.

Having estimated the three synchrophasors of the phase quantities, the positive sequence synchrophasor ${\overline{X}}_{+}$ is evaluated through the Fortescue transformation (6). From the previous considerations, the estimates ${\hat{\overline{X}}}_p$ of the phase synchrophasors may contain disturbances due to the infiltration of the image components. Similarly to what explained in Section 2.2, this effect is significantly attenuated in the positive sequence synchrophasor measurement ${\hat{\overline{X}}}_{+}$ thanks to the weak unbalance level of real-world three-phase quantities. However, explaining how this cancellation occurs is troublesome, since the IpDFT algorithm is inherently nonlinear and closed-form expressions cannot be easily derived. Finally, it is also worth noting that applying the IpDFT algorithm to the phase signals produces three (different) frequency estimates, which have to be processed (averaged) in order to obtain a unique value.

Similarly to the TF approach, the IpDFT algorithm can be favorably applied to the SV signal ${\overline{x}}_{S{V}_0}$ in a stationary reference frame, as in the block diagram reported in Figure 4; this enables a direct estimation of the positive sequence synchrophasor. As in Section 2.4, ${\overline{x}}_{S{V}_0}$ is computed from the phase waveforms by means of the SV transformation (21) with $\beta =0$. Its model, which also in this case is assumed to be valid over the C nominal cycle interval centered around the reporting instant, results by applying to (38) the SV transformation on the same reference frame, hence:

$$\begin{array}{cc}\hfill {\overline{x}}_{S{V}_0}\left(t\right)& ={\overline{X}}_{ph,+}\left(t_r\right)e^{j\omega \left(t_r\right)t_r}{e}^{j(\omega \left(t_r\right)-{\omega }_0)t}{e}^{j{\omega }_{0}t}+{\overline{X}}_{ph,-}^{\ast }\left(t_r\right)e^{-j\omega \left(t_r\right)t_r}{e}^{-j(\omega \left(t_r\right)-{\omega }_0)t}{e}^{-j{\omega }_{0}t}+{\overline{d}}_{S{V}_0}\left(t\right)\hfill \\ \hfill & ={\overline{X}}_{+}\left(t_r\right)e^{j{\omega }_{0}t}+{\overline{X}}_{-}^{\ast }\left(t_r\right)e^{-j{\omega }_{0}t}+{\overline{d}}_{S{V}_0}\left(t\right)\hfill \end{array}$$

(39)

The next step is considering the $N_w=C{M}_C$ samples of the SV signal collected around $t_r$ and, having applied a smoothing window defined by the sequence $w\left[n\right]$ ($n=0,\dots ,N_w-1$), computing the DFT bins ${\overline{X}}_{S{V}_0}(k,t_r)$ whose corresponding frequencies are close to $f_0$, namely whose indexes k are near to C. Let us assume that, as typically happens in practical implementations, $M_C$ and thus also $N_w$ are even numbers.

$${\overline{X}}_{S{V}_0}(k,t_r)=\sum _{n=0}^{N_w-1}w\left[n\right]{\overline{x}}_{S{V}_0}(t_r+(n-N_w/2)T_s)e^{-j2\pi \frac{kn}{N_w}}$$

(40)

The key underlying assumption of IpDFT algorithms is that long-range leakage is negligible for the considered indexes k. When comparing (39) with (38), it is evident that this assumption is more easily met as long as the IpDFT algorithm is applied to the SV signal instead of a phase waveform. In fact, thanks to low unbalance level of power systems, in the first case the disturbance produced by the counter-rotating component is considerably lower in relative value.

Therefore, neglecting spectral interference, it results ($\omega =\omega \left(t_r\right)$ in the following):

$$\begin{array}{cc}\hfill {\overline{X}}_{S{V}_0}(k,t_r)& \approx {\overline{X}}_{ph,+}\left(t_r\right)\overline{W}(k-\lambda )e^{-j\omega N_w{T}_s/2}\hfill \\ \hfill & ={\overline{X}}_{ph,+}\left(t_r\right)\overline{W}(k-\lambda )e^{-j\pi \lambda }\hfill \end{array}$$

(41)

where $\overline{W}\left(\gamma \right)$ is the discrete-time Fourier transform of $w\left[n\right]$ computed in the generalized bin $\gamma$, while $\lambda =f{N}_w{T}_s$. Now, let us suppose to have employed a periodic window having $w\left[0\right]=0$: it is worth noting that this property holds true for many widely employed windows, even including those belonging to the Rife-Vincent class I family. In this case, it is possible to write:

$${\overline{X}}_{S{V}_0}(k,t_r)\approx {(-1)}^k{\overline{X}}_{ph,+}\left(t_r\right)W(k-\lambda )$$

(42)

where $W\left(\gamma \right)$ is the amplitude response of the window. Let us assume that $k_0$ is the closest integer to $\lambda$, namely the index of the highest DFT component. Writing $\lambda =k_0+\delta \left(t_r\right)$ ($\left|\delta \right(t_r\left)\right|<1/2$), $\delta$ can be estimated from the ratio between the magnitudes of the two largest bins with the following equations:

$$\begin{array}{cc}\hfill \xi \left(t_r\right)& =\frac{\left|{\overline{X}}_{S{V}_0}(k_0+l\left(t_r\right),t_r)\right|}{\left|{\overline{X}}_{S{V}_0}(k_0,t_r)\right|}\approx \frac{\left|W(l\left(t_r\right)-\delta \left(t_r\right))\right|}{\left|W(-\delta \left(t_r\right))\right|}=g\left(\delta \left(t_r\right)\right)\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \hat{\delta }\left(t_r\right)& =g^{-1}\left(\xi \left(t_r\right)\right)\hfill \end{array}$$

(44)

where $l\left(t_r\right)=sign(\left|{\overline{X}}_{S{V}_0}(k_0+1,t_r)\right|-\left|{\overline{X}}_{S{V}_0}(k_0-1,t_r)\right|$).

The frequency is then estimated as:

$$\hat{f}\left(t_r\right)=f_0+\frac{\hat{\delta }\left(t_r\right)}{N_w{T}_s}$$

(45)

Using $\hat{\delta }$, the positive sequence synchrophasor in the reporting instant is obtained through scalloping loss compensation and phase shift of the highest DFT component:

$${\hat{\overline{X}}}_{+}\left(t_r\right)={(-1)}^k\frac{{\overline{X}}_{S{V}_0}(k_0,t_r)}{W\left(\hat{\delta }\left(t_r\right)\right)}$$

(46)

Finally, the ROCOF is obtained by using the previously described algorithm to perform frequency measurements in the time instants $t_r-S{T}_s$ and $t_r+S{T}_s$ while computing a discrete-time derivative, thus applying the linear-phase FIR filter $D_S^1$ (introduced in Section 2.2) to the frequency estimates:

$$\widehat{\mathrm{ROCOF}}\left(t_r\right)=\frac{\hat{f}(t_r+S{T}_s)-\hat{f}(t_r-S{T}_s)}{2T_s}.$$

(47)

As in Section 2.2, the algorithm latency is given by half the window length plus the delay required by the discrete-time derivative used to compute ROCOF, that is $L=(N_w/2-1+S)T_s=(C{M}_C/2+S-1)T_s$.

3. Simulation Results

3.1. Implementation of the Algorithms

The target of this section is comparing the capabilities of the previously discussed PMU algorithms in estimating the positive sequence synchrophasor, frequency and ROCOF under different conditions. The tests waveforms follow the standard [7], with reference to P compliance class and maximum reporting rate, but additional cases are also included. Accuracy of the estimates is quantified in terms of total vector error (TVE), absolute frequency error (|FE|) and absolute ROCOF error (|RFE|) defined, with reference to the generic measurement instant $t_r$, as:

$$\begin{array}{c}\mathrm{TVE}\left(t_r\right)=\frac{|{\hat{\overline{X}}}_{+}\left(t_r\right)-{\overline{X}}_{+}\left(t_r\right)|}{X_{+}\left(t_r\right)}\\ |\mathrm{FE}\left(t_r\right)|=|\hat{f}\left(t_r\right)-f\left(t_r\right)|\\ |\mathrm{RFE}\left(t_r\right)|=|\widehat{\mathrm{ROCOF}}\left(t_r\right)-\mathrm{ROCOF}\left(t_r\right)|\end{array}$$

(48)

The considered algorithms have been applied with the same sampling rate $f_s=200·f_0$; a reporting interval equal to the sampling interval (that is $Q=1$, $T_{RR}=T_s$) is employed for a more accurate comparison. In order to have similar latencies, so that the comparison is more significant, all the methods except for IEC-P have been implemented to use three nominal cycles for each synchrophasor estimation. Filter design has key importance in the performance achieved by the SV-F algorithm. The first filter H is a 181 tap lowpass linear-phase equiripple having 2 Hz passband frequency with $2·{10}^{-3}$ passband ripple, 50 Hz stopband frequency and $0.03$ stopband ripple. Filters M and P are identical 421 tap lowpass equiripple linear-phase filters with 2 Hz passband frequency, $0.01$ passband ripple, 50 Hz stopband frequency and $0.03$ stopband ripple. All the coefficients of H, M and P have been scaled in order to have unit DC gain; this ensures theoretically zero error when the input is sinusoidal and positive sequence, also under off-nominal frequency. Filter F is a 421 tap linear-phase, equiripple partial-band differentiator. Its coefficients have been scaled so that the output is one when the input is $x\left[n\right]=n{T}_s$; it results in theoretically zero frequency estimation error when input is sinusoidal, positive sequence but frequency is different from its rated value. Finally, filter R is a 421 tap, linear-phase, equiripple 2nd order partial-band differentiator. Its coefficients have been scaled in order to return unitary output when input is $x\left[n\right]={\left(n{T}_s\right)}^2/2$; this corresponds to zero RFE in the presence of a purely positive sequence input exhibiting a frequency ramp. The SV-TF estimator has been designed in order to exploit the additional flexibility enabled by the native three-phase approach: $K_{+}=3$ guarantees a good modeling of the time-dependency of the positive sequence synchrophasor, while $K_{-}=1$ is generally sufficient to prevent the infiltration of the negative sequence term in the positive sequence synchrophasor estimate. When considering the IEC-P and SV-IpDFT methods, $S=1$ is employed, namely the finite differences in the reporting instant $t_r$ are computed from the estimates in $t_r-T_s$ and $t_r+T_s$. Finally, a periodic Hann window is used for the implementation of the SV-IpDFT algorithm.

3.2. Tests under Steady-State Conditions

These tests require applying three-phase waveforms whose parameters, which are the positive sequence phasor, the frequency and the ROCOF, have constant values (it implies $\mathrm{ROCOF}=0$). The target is evaluating the capabilities of the PMU algorithms in returning accurate estimates (ideally constant and equal to those present in the test waveform) in the presence of different disturbing factors.

3.2.1. Off-Nominal Frequency

The first set of tests is focused on testing the accuracy achieved by the different algorithms when the input is sinusoidal, positive sequence, but the frequency differs from its rated value. In particular, different frequencies have been considered. Figure 5 reports the maximum TVE under steady-state conditions for all the algorithms; frequency spans the range $[48,52]\mathrm{Hz}$ with a step of $0.2\mathrm{Hz}$. Since the negative sequence input is zero, positive sequence synchrophasor estimates are not affected by spectral interference. In particular, TVE values are negligible for the SV-F and SV-IpDFT, since they provide a theoretically exact compensation of the scalloping loss, which is the only error source under these conditions. TVE is not negligible (albeit it is still very low) as far as the SV-TF (<$7.1·{10}^{-3}$) and IEC-P (<$4.5·{10}^{-3}$) methods are concerned. The reason is that the SV-TF technique does not compensate for scalloping loss due to the employed filters, while IEC-P uses an approximate formula, thus leading to its mitigation but not to a complete cancellation. All the algorithms result in negligible RFE values (in the order of ${10}^{-8}\mathrm{Hz}/\mathrm{s}$ or even lower), while the |FE| has a meaningful value only for SV-TF ($0.078\mathrm{mHz}$ at $48\mathrm{Hz}$). In this case, the reason is undermodeling of the positive sequence synchrophasor because of its truncated Taylor expansion and the consequent error infiltration in (36), which is used for obtaining the frequency estimate. It is worth noting that both FE and TVE values can be reduced by obtaining the filter bank from a WLS solution (e.g., weights corresponding to the coefficients of a Hann window), as previously mentioned in Section 2.4.

3.2.2. Harmonic Disturbances

Another important set of tests is focused on assessing the impact of harmonic disturbances on positive sequence synchrophasor, frequency and ROCOF measurements obtained with the different algorithms. SV-based PMU algorithms (SV-F, SV-TF and SV-IpDFT) estimate the positive sequence synchrophasor, frequency and ROCOF by processing the SV signal. Therefore, in order to understand their behavior, it is significant to write the expression of ${\overline{x}}_{SV}$ under steady-state conditions in the presence of harmonic disturbances. It results from (22) where:

$$\begin{array}{ccc}\hfill {\overline{X}}_{+}\left(t\right)={\overline{X}}_{ph,+}{e}^{j(\omega -{\omega }_0)t}& \hfill & \hfill {\overline{X}}_{-}\left(t\right)={\overline{X}}_{ph,-}{e}^{j(\omega -{\omega }_0)t}\end{array}$$

(49)

$${\overline{d}}_{S{V}_0}\left(t\right)=\sum _{h=2}^{\infty }\left[{\overline{X}}_{ph,h,+}{e}^{jh\omega t}+{\overline{X}}_{ph,h,-}^{\ast }{e}^{-jh\omega t}\right]$$

(50)

${\overline{X}}_{ph,h,-}$ is the negative sequence phasor of the fundamental; ${\overline{X}}_{ph,h,+}$ and ${\overline{X}}_{ph,h,-}$ are, respectively, the positive and negative sequence phasors of the hth order harmonic. Therefore, positive sequence terms produce complex exponential contributions to ${\overline{x}}_{SV}$ which rotate with angular speed $h\omega$, hence in positive direction; conversely, negative sequence components are responsible for complex exponential terms which rotate with angular speed $-h\omega$, thus in negative direction. Zero sequence terms do not appear in ${\overline{x}}_{SV}$.

As mentioned before, SV based PMU algorithms apply some kind of bandpass filtering to ${\overline{x}}_{S{V}_0}$ in order to extract only the components whose rotational speeds are close to ${\omega }_0$, thus representing the centre of the passband. It is worth noting that, for a given harmonic order h, negative sequence components are easier to be filtered with respect to positive sequence terms, since they are $2{\omega }_0$ farther with respect to the middle of the passband. Similar consideration applies also to the IEC-P method. In fact, with some manipulations, it is possible to show that the provided positive sequence synchrophasor estimate corresponds to the application of the triangular FIR filter to the SV signal whose instantaneous angular position of the reference frame is $\beta \left(t\right)={\omega }_{0}t$.

In practical applications, the three-phase set is very close to be symmetric: this means that phase b and c waveforms are obtained by time-shifting that of phase a by $-1/\left(3f\right)$ and $1/\left(3f\right)$, respectively. Under this assumption, it can be shown that:

• harmonics whose orders are multiple of three are purely zero sequence contributions;
• harmonics whose orders are given by $h=3k+1$ are purely positive sequence contributions;
• harmonics whose orders are given by $h=3k-1$ are purely negative sequence contributions.

In order to test the performance of the different algorithms in the presence of harmonic disturbances, three-phase test signals characterized by a positive sequence contribution at frequency f with a superimposed symmetric hth order harmonic (h ranging from 2 to 50) having $1\%$ relative amplitude have been applied. For each test signal and algorithm, maximum TVE, |FE| and |RFE| values over a 2-s time interval have been evaluated and compared. It is worth reminding that because of the aforementioned considerations, injected harmonics whose orders are multiple of 3 do not affect the estimates. A particular situation for methods IEC-P and SV-IpDFT occurs when frequency is equal to its rated value. In fact, they apply a frequency shift of $-f_0$ to the SV signal and then either a two-cycle triangular or a three-cycle Hann window having frequency response zeros located at multiples of $f_0$, thus in correspondence of the harmonics, which are fully removed. As a result, positive sequence synchrophasor, frequency and ROCOF estimates obtained with these techniques are intrinsically immune to harmonic disturbances when $f=f_0$. Under the same conditions, the SV-F method provides very good synchrophasor and frequency estimates, with maximum TVE and |FE| below $2·{10}^{-3}\%$ and $0.2\mathrm{mHz}$; |RFE| is significant, reaching a maximum value of $0.39\mathrm{Hz}/\mathrm{s}$, which is still compliant with the P-class limit of the IEC/IEEE standard. The SV-TF method results in the highest errors, with maximum TVE of $3.6·{10}^{-2}\%$ (thus well below the $1\%$ limit prescribed by the IEC/IEEE standard) in the presence of the $2\mathrm{nd}$ order harmonic, and then decaying for the higher-order components thanks to the response of the TF filter. Conversely, |FE|, exceeds $13\mathrm{mHz}$ at the $4\mathrm{th}$ order harmonic, namely significantly higher with respect to the $5\mathrm{mHz}$ limit; |RFE| almost reaches $0.6\mathrm{Hz}/\mathrm{s}$, also in this case higher than the corresponding limit. However, harmonic rejection capability of the SV-TF approach can be dramatically improved by applying weighting coefficients to the samples, such as those of a Hann window.

Harmonic rejection capability has been tested also under off-nominal frequency (not requested by the IEC/IEEE standard), namely considering $f=49\mathrm{Hz}$; results are reported in the following Figure 6, Figure 7 and Figure 8.

The first consideration is that from the previous graphs the envelope of the errors reflects the response of the filters employed by the different algorithms. The somewhat jagged trend occurs since, according to their harmonic order, superimposed disturbances in the SV signal have alternatively positive angular speed, negative angular speed, or zero magnitude. Under these conditions, it is clear that IEC-P and SV-IpDFT are no longer able to fully reject harmonic disturbance. When looking at the TVE values (Figure 6), the SV-IpDFT algorithm still results in excellent accuracy: maximum is about $3·{10}^{-4}\%$ at the $4\mathrm{th}$ order harmonic; the IEC-P method reaches a larger maximum TVE slightly above $2·{10}^{-3}\%$, which is still pretty low. As for the SV-F algorithm, performance is very close to that achieved at nominal frequency. The same applies to the SV-TF algorithm, with the highest error now occurring at the $4\mathrm{th}$ order harmonic. From Figure 7, maximum frequency errors for SV-F and SV-TF algorithms are similar to those with $f=f_0$, while the SV-IpDFT approach guarantees the best results; this means that the adopted Hann window is very effective in preventing spectral interference. Frequency estimates provided by the IEC-P method are affected by significant errors, but they remain below $2.5\mathrm{mHz}$. Finally, when analyzing the results in terms of RFE, values obtained with the SV-F algorithm are also in this case very similar to those measured with nominal frequency; maximum |RFE| of the SV-TF technique slightly increases to almost $0.8\mathrm{Hz}/\mathrm{s}$ at the $4\mathrm{th}$ order harmonic. SV-IpDFT ensures the best performance also in terms of ROCOF estimate, while the |RFE| values corresponding to the IEC-P algorithm are extremely large. The impact of harmonic disturbances affecting the positive sequence synchrophasor estimates, which are no more suppressed by the triangular window, is highly magnified by the second order differentiator that allows computing the ROCOF.

3.2.3. Unbalance

Asymmetries between the three phases may be present in real-world applications: thus, assuming purely sinusoidal steady-state conditions, the three-phase input to the PMU is not purely positive sequence, but it may contain also negative and zero sequence contributions. From the considerations of Section 3.2.2, it can be noticed that, for all the algorithms, the zero sequence contribution is inherently rejected, while the negative sequence term results as a complex exponential disturbance, whose rotational speed is $-(\omega +{\omega }_0)$, affecting the positive sequence synchrophasor, frequency and ROCOF estimates. Therefore, the target of these tests is assessing the capability of the different algorithms to provide accurate estimates also in the presence of three-phase asymmetry; for the sake of completeness, frequency f is swept in the range $[48,52]\mathrm{Hz}$, as in the previous tests. Accuracy is quantified by the maximum TVE, |FE| and |RFE| evaluated on a 2-s time interval for each test condition and PMU algorithm. It is worth noting that the IEC-P and SV-IpDFT methods are intrinsically immune to the infiltration of the negative sequence component when $f=f_0$, thanks to the zeros of the employed triangular or Hann window. Also the SV-TF method results in virtually zero error in this case: in fact, the underlying model produces zeros located at $-f_0$ in the frequency responses in the resulting complex-valued filters which are applied to the SV signal in order to obtain the positive sequence synchrophasor and its derivatives.

First of all, amplitude unbalance has been considered: phase a magnitude is set $10\%$ larger with respect to the others. TVE is still very low for all the algorithms and frequency values (below $1.2·{10}^{-2}\%$) and hence not shown, while results in terms of FE and RFE are much more significant; they are reported in Figure 9 and Figure 10 as functions of f. As mentioned above, the negative sequence acts as a disturbance in the SV signal located at frequency $-(f+f_0)$ whose relative magnitude is $3.2\%$ with respect to the positive sequence component; it has to be rejected by the filtering stages of the different algorithms. As expected, when $f=f_0$ all the algorithms except for SV-F result in zero error and exhibit a monotonic increase of errors as long as f deviates from $f_0$. For significant values of $\Delta f$, SV-F becomes the most accurate algorithm, both for frequency and ROCOF measurement. The behavior of SV-IpDFT is fairly good as far as the frequency estimate (|FE| below $1\mathrm{mHz}$), while |RFE| values are significant, exceeding $0.5\mathrm{Hz}/\mathrm{s}$; the reason is the spectral interference produced by the negative sequence term (it becomes larger as the difference between f and $f_0$ increases) that the algorithm assumes to be negligible, whose impact is magnified by the numerical differentiation. The triangular window adopted by the IEC-P is not very effective in reducing the oscillation produced by the negative sequence component in the positive sequence synchrophasor phase angle estimate under off-nominal frequency conditions. When numerical first and, even more, second order discrete differentiation is applied to obtain frequency and ROCOF estimates, the magnitude of this oscillation increases noticeably: it is not surprising that the IEC-P method achieves the worst |RFE| values. The SV-TF method achieves the highest frequency error because of the response of the filter used for obtaining the positive sequence synchrophasor first-order derivative, which in turns is used for frequency estimation. Better results would have been obtained by designing a filter with a wider stopband around $-f_0$. ROCOF measurement is fairly good (second only to the SV-F method), but it is worth highlighting that significantly lower |FE| and |RFE| can be obtained by adopting, as mentioned before, Hann weighting in filter design; values of $0.06\mathrm{mHz}$ and $0.01\mathrm{Hz}/\mathrm{s}$ at $52\mathrm{Hz}$ are achieved in a worst-case scenario.

Similar results, confirming the above discussion, can be found also in case of phase unbalance, thus reducing the angle of phase a by ${10}^{°}$ with respect to purely sinusoidal positive sequence conditions. Since the unbalance level is higher in this case ($5.6\%$), the errors are expected to be higher too. TVE increases but it is still very small (below $1.8·{10}^{-2}\%$); Figure 11 and Figure 12 show the |FE| and |RFE| values. It is worth highlighting that results show that the magnitudes of all the errors increase proportionally to the unbalance level; this has been extensively studied in in [34] for IEC-P, SV-F, and TF-based algorithms and analytical expression to predict their values are also available. This consideration applies with excellent accuracy also for the SV-IpDFT estimator, even though it is intrinsically nonlinear. Once again, the performance of SV-TF can be enhanced if the filters are obtained with the WLS approach (weights corresponding to the squares of the Hann window coefficients) and |FE| and |RFE| thus become lower than $0.12\mathrm{Hz}$ and $0.03\mathrm{Hz}/\mathrm{s}$, respectively, again reflecting the stronger impact of the disturbance with respect to the previous unbalance test case.

3.2.4. Wideband Noise

Tests with additive uniform white noise permit to verify the immunity of the algorithms to disturbances that are not narrowband (such as unbalance and harmonics discussed in Section 3.2.2 and Section 3.2.3) but correspond instead to wideband noise which may be due to different sources, such as thermal noise of analog signal conditioning stages and quantization noise of digitization stages. Purely sinusoidal, positive sequence input at nominal frequency has been considered and independent zero-mean uniform noise has been superimposed to the samples. This permits to highlight the impact of the noise alone having considered different yet realistic signal-to-noise ratios (SNRs). TVE, FE and RFE have been evaluated over a 20-s test duration, and their Root Mean Square (RMS) values have been computed. Maximum errors are not considered in this case, since they have intrinsically a considerably worse statistical behavior. Figure 13 reports the impact of noise on the RMS value of TVE. It is clear that the impact depends on the equivalent noise bandwidth of the filters or cascades of filters used for synchrophasor estimation; as expected, the error increases linearly with the noise level. In this respect, SV-IpDFT (characterized by the narrowest equivalent bandwidth) achieves the lowest error values. On the opposite side, the SV-TF technique suffers from the largest TVEs: it is somewhat expected, since TF filters have been designed with dynamic performance in mind. In general, increasing the expansion order widens the equivalent bandwidth and thus also the sensitivity to noise. Furthermore, it is interesting to notice that the aforementioned Hann weighting in SV-TF leads, in the presence of wideband noise, to an error which is $25\%$ higher, meaning that the improved characteristics of dynamic synchrophasor tracking have, as a side effect, lower noise immunity.

Figure 14 and Figure 15 report the RMS errors for frequency and ROCOF, respectively. Also in this case, performance depends on the equivalent noise bandwidths (this applies to linear but also nonlinear algorithms trough small-signal linearization) that characterize the different methods for estimating frequency and ROCOF, and as expected the values of the errors are proportional to the noise level. It can be noticed that the differences among the algorithms under test are larger when RMS FE and RFE are considered instead of the RMS TVE.

When looking at the RMS FE, the SV-F algorithm results in the lowest error values, while the other techniques achieve similar performance (SV-TF is marginally the worst). Focusing on the RMS RFE, it is clear that IEC-P and SV-IpDFT suffer from the largest errors: the reason is that ROCOF is obtained with first (SV-IpDFT) or even second-order (IEC-P) discrete-time derivatives, whose outputs are notably highly sensitive to the presence of input noise. The SV-TF with the selected expansion orders provides the best ROCOF estimate in the presence of noise, with SV-F resulting in slightly higher RMS RFE. Also in this case, Hann weighting for SV-TF leads to significantly worse errors, with a growth up to $78\%$ and $176\%$ for |FE| and |RFE|, respectively. This result highlights how the values of the adopted weighting coefficients must be carefully considered: they can be very helpful for synchrophasor estimation algorithm design (as proven in previous sections), but only if high SNR is expected.

3.3. Tests under Dynamic Conditions

One of the breakthroughs of PMUs is that they enable monitoring the time evolution of the power system. For this purpose, the employed algorithms should be able to track the value of synchrophasor, frequency and ROCOF also when the parameters defining the input signal are time-varying. Tests under dynamic conditions are aimed at evaluating the accuracy in these scenarios.

3.3.1. Amplitude and Phase-Angle Modulations

Electrical signals in the presence of power system oscillations can be modeled as in (1). Therefore, the tracking capability of PMUs under these conditions can be assessed by applying three-phase balanced sinewaves with sinusoidal amplitude or phase angle modulation.

$${\mathbf{x}}_{abc}\left(t\right)=\sqrt{2}X(1+k_{x}cos\left(2\pi f_{m}t\right))\left[\begin{array}{c}cos({\omega }_{0}t+k_{a}cos(2\pi f_{m}t-\pi ))\\ cos({\omega }_{0}t-\frac{2}{3}\pi +k_{a}cos(2\pi f_{m}t-\pi ))\\ cos({\omega }_{0}t+\frac{2}{3}\pi +k_{a}cos(2\pi f_{m}t-\pi ))\end{array}\right]$$

(51)

where $k_x$ and $k_a$ are the amplitude and phase-angle modulation indexes, respectively, X is the unmodulated RMS amplitude and $f_m$ is the modulation frequency. Applying the SV transformation to (51), it results:

$${\overline{x}}_{SV}\left(t\right)=\sqrt{3}X(1+k_{x}cos\left(2\pi f_{m}t\right))e^{j(k_{a}cos(2\pi f_{m}t-\pi )+{\omega }_{0}t-\beta \left(t\right))}={\overline{X}}_{+}\left(t\right)e^{j({\omega }_{0}t-\beta \left(t\right))}$$

(52)

Tests are performed with different values of $f_m$, which is swept from 0 (unmodulated) to $2\mathrm{Hz}$ with $0.1\mathrm{Hz}$ step in order to evaluate the measurement bandwidths of the different techniques. For each algorithm and condition, accuracy is expressed by the maximum values of TVE, |FE| and |RFE| over a 10-s test duration.

Let us firstly consider amplitude modulated signals, with $k_x=0.1$ and $k_a=0$; in this case, the algorithms are required to extract a bandpass signal from ${\overline{x}}_{SV}$, whose spectral content is located between $f_0\pm f_m$. This is evident for SV based algorithms, but it is true also for the IEC-P method as long as the positive sequence synchrophasor, frequency and ROCOF measurements are considered. The SV-TF immediately appears as perfectly suited to cope with this test, since TF filters are designed with dynamic performance in mind; in particular, the filter that allows positive sequence synchrophasor measurement exhibits a flat frequency response around $f_0$ (due to multiple null derivatives according to the selected expansion order). Figure 16 reports the TVE values for all the algorithms; as expected, errors monotonically increase with the modulation frequency and, furthermore, it is clearly noticeable that SV-TF provides the most accurate synchrophasor measurement. IEC-P and SV-F achieve very close error values, but it should be highlighted that the results obtained with SV-F strongly depend on the filter design, which is obtained as a tradeoff between steady-state and dynamic performance according to the specific needs. It is not surprising that the highest error values are reached by the SV-IpDFT method, which relies on a steady-state model of the SV signal.

Frequency errors are either small or negligible (<$2·{10}^{-3}\mathrm{mHz}$) for all algorithms but SV-IpDFT, which shows higher |FE| up to $13.25\mathrm{mHz}$ due to its underlying steady-state assumption that is not able to effectively cope with amplitude modulated input signals. As a consequence, |RFE| achieved by SV-IpDFT is also huge (it can easily go beyond $10\mathrm{Hz}/\mathrm{s}$, when $f_m$ gets larger than a few tenths of hertz). Conversely, the other algorithms result in negligible RFE values (order of magnitude ${10}^{-5}\mathrm{Hz}/\mathrm{s}$, or less).

More varied is the performance of the methods under phase angle modulations, as shown in Figure 17, Figure 18 and Figure 19. Again, the dynamic model which the SV-TF is based on results in the most accurate synchrophasor estimations: the TF filters preserve the bandpass characteristics of the SV signal also in this case. As for the other algorithms, TVE values are low, but they show how the synchrophasor tracking properties can be significantly different. In particular, SV-IpDFT shows the highest errors because of the steady-state signal representation, which shows its major weakness under dynamic conditions. The same considerations hold true also for frequency estimation: the SV-TF is also in this case the less sensitive to the sinusoidal frequency oscillation (see Figure 18). The other algorithms show similar errors, with the SV-F approach being the least accurate, even if better dynamic performance would be enabled by different filter design. As far as ROCOF measurements is concerned, SV-IpDFT is again more prone to errors for the aforementioned reasons. The other algorithms are quite similar and it is interesting to notice that SV-TF is no longer the best one, mainly because, when sinusoidal variations are involved, the truncation of the Taylor expansion limits the tracking capabilities. A higher expansion order might help reducing the RFE values.

3.3.2. Frequency Ramp

This test requires applying a three-phase positive sequence sinusoidal waveform whose frequency varies from $48\mathrm{Hz}$ to $52\mathrm{Hz}$ (between the instants $t_{start}=3\mathrm{s}$ and $t_{end}=7\mathrm{s}$) with constant ROCOF equal to $1\mathrm{Hs}/\mathrm{s}$, starting from and reaching steady-state conditions. Performance is evaluated in terms of TVE, |FE| and |RFE|, but measurements whose reporting instants are before $t_{start}+2/f_0$ (thus two nominal cycles after the start of the ramp) or after $t_{end}-2/f_0$ (that is two nominal cycles before the end of the ramp) are excluded from the analysis, as prescribed by [7]. Under this condition, the TVE values resulting from the algorithms IEC-P, SV-F and SV-IpDFT, which are not intrinsically based on a dynamic synchrophasor model, are significantly degraded with respect to off-nominal frequency conditions, as it can be noticed from Figure 20. On the contrary, the SV-TF method, whose filters are derived from a Taylor expansion of the synchrophasor, shows almost the same performance as that reported in Figure 5. In fact, the dynamics corresponds to a slowly-varying fundamental frequency which can be embedded by the TF synchrophasor model. The SV-TF method reaches the highest |FE| values, which are very close to those obtained under off-nominal frequency conditions (see Section 3.2.1); also the corresponding maximum |RFE| is not negligible (always <$0.06\mathrm{Hz}/\mathrm{s}$), since the Taylor expansion truncation and the linearization used for ROCOF computation (see (37)) become more significant when higher-order derivatives are involved. It is interesting to notice that thanks to the normalization of filter R coefficients, SV-F has an ideal behaviour in the presence of constant ROCOF, and thus leading to virtually zero RFE in the presence of a frequency ramp. ROCOF measurement provided by the SV-IpDFT algorithm has a small spike (resulting in $0.025\mathrm{Hz}/\mathrm{s}$ RFE) when frequency reaches $50\mathrm{Hz}$. The reason is that the frequency estimation error changes sign exactly at the nominal frequency, and this affects the ROCOF measurement, which is obtained by numerical differentiation.

3.3.3. Step Tests

In these tests, the applied input signal is sinusoidal and positive sequence with rated frequency, but it exhibits a step variation of either its magnitude or phase angle in a predetermined time instant $t_{step}$. The target is assessing the behavior of the PMU algorithms in the presence of fast transients, such as those triggered by switching and faults. As prescribed by [7], step magnitudes are $10\%$ in magnitude and ${10}^{°}$ in phase; performance are evaluated thus assessing the times required by TVE, FE and RFE to return below the corresponding steady-state limits (response times), the delay time and the magnitude of overshoot/undershoot which may be present in the estimates.

Figure 21 shows the magnitude of the positive sequence synchrophasor estimated by the different algorithms in the presence of a $+10\%$ amplitude step, starting from $1\mathrm{p}.\mathrm{u}.$ magnitude. All the algorithms have negligible delay time thanks to a proper compensation of the group delay: halfway of the step is reached within 1 sample ($100\mathsf{\mu }\mathrm{s}$) from the step instant $t_{\mathrm{step}}=1\mathrm{s}$. The magnitude estimate obtained with the SV-TF method exhibits undershoot in the pre-transition region and overshoot in the post-transition region (as evident from Figure 21) with a maximum of $8.08\%$ of the step size. Such effect, which is due to the characteristics of the TF filters, can be mitigated by using Hann weighting as discussed above. All the other algorithms do not show significant under/overshoot: they are easily within the limit indicated by the standard ($5\%$ of the step magnitude).

The TVE, FE and RFE response times for all the algorithms reported in

Table 1. Response times for $+10\%$ magnitude step test.

Algorithm	Response Time [ms]
	TVE	Frequency	ROCOF
IEC-P	21.7	0.0	0.0
SV-F	21.6	0.0	0.0
SV-TF	22.6	52.3	58.7
SV-IpDFT	28.4	52.5	58.1

. It is worth highlighting that SV-F and IEC-P have zero error both in frequency and ROCOF since in this case phase angle is decoupled from amplitude estimate. SV-TF has response times that strongly depend on the order of expansion and on the possible use of weights, which can easily lead to maximum |FE| and |RFE| lower than the corresponding steady-state limits: in that case, frequency and ROCOF response times becomes equal to zero.

Figure 22 reports the positive sequence synchrophasor phase angle estimated by the considered algorithms when a $-{10}^{°}$ phase step is applied. Analogously to amplitude step test, overshoot and undershoot are noticeable only for SV-TF algorithm, which leads to a maximum value of $8.1\%$ with respect to the step magnitude, thus above the limit of the standard. Response times are reported in

Table 2. Response times for $-{10}^{\circ }$ phase step test.

Algorithm	Response Time [ms]
	TVE	Frequency	ROCOF
IEC-P	26.4	39.8	40.0
SV-F	31.0	56.8	60.0
SV-TF	53.0	60.0	60.0
SV-IpDFT	34.6	51.0	56.2

, where SV-TF shows the longest TVE response time because of the secondary lobes present in the time evolution of TVE (see Figure 23) that go beyond the $1\%$ steady-state limit prescribed by [7]. The Hann weighting allows improving SV-TF dynamic response, leading to a TVE response time of $15.0\mathrm{ms}$ and to a maximum under/overshoot of $4.4\%$. In addition, it also allows reducing the frequency and ROCOF response times, which, as reported in Table 2, are very close to the window length for all the algorithms (namely two and three nominal cycles for IEC-P and the other methods, respectively). In fact, FE and RFE promptly react to the abrupt phase change as soon as the step enters in the sample window, resulting in large deviations from actual frequency and ROCOF (up to about $1\mathrm{Hz}$ and $80\mathrm{Hz}/\mathrm{s}$ for |FE| and |RFE|, respectively).

4. Conclusions

The paper has presented a thorough comparison among the performance achieved by four low-latency algorithms for positive sequence synchrophasor, frequency and ROCOF estimation. Three of them are designed directly on the SV signal model, while the fourth one is the reference P-class method of the latest IEC/IEEE standard, which is used as a benchmark. The paper has shown how, depending on the actual conditions of the three-phase signals and on the specific assumptions of each method, results can differ noticeably. It is extremely important to notice that the algorithms are characterized by parameters that need to be finely tuned according to the desired performance target. It is also clear that techniques designed to track the synchrophasor dynamics better preserve the bandpass characteristics of the positive sequence synchrophasor, but at the expense of a lower immunity to disturbances. The paper has highlighted that, depending on the expected level of narrowband or wideband disturbances, some methods can be prone to errors. Countermeasures, such as different weighting functions have been also presented, while keeping in mind that a tradeoff between bandwidth and immunity, between fast response to abrupt changes and small artifacts in the response itself, is always to be sought. The presented algorithms have, in general, good performance, and peculiar characteristics, but some unexpected outcomes have been highlighted by the analysis (e.g., those about the RFE under phase modulation). The performed comparison allows understanding in detail the potentialities of these methods that are simple and lightweight and thus likely to be implemented in commercial PMUs for three-phase systems.