An adaptive filters based PMU algorithm for both steady-state and dynamic conditions in distribution networks

Electrical Power and Energy Systems 117 (2020) 105714

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

An adaptive filters based PMU algorithm for both steady-state and dynamic conditions in distribution networks

T

Yinfeng Wanga, Chao Lua, , Innocent Kamwab, Chen Fangc, Ping Lingc ⁎

a

Department of Electrical Engineering, Tsinghua University, Beijing, China Hydro-Québec Research Institute (IREQ), Varennes, Quebec, Canada c Shanghai Electric Power Research Institute, State Grid Corporation of China, Shanghai, China b

ARTICLE INFO

ABSTRACT

Keywords: Phasor estimation ESPRIT Finite impulse response filter Taylor Fourier transform Transient recognition Distribution network

Although synchronized phasor measurement unit (PMU) technology is widely applied in transmission networks, conventional PMU algorithms suffer from new challenges (e.g. large frequency deviation, multiple interferences and dynamic behaviors, etc.) if directly applied to distribution networks. On the basis of predecessors' research, this paper presents a comprehensive PMU algorithm, using adaptive filters to simultaneously satisfy the highprecision and dynamic-response requirements of phasor estimation in distribution networks. In this algorithm, out-of-band (OOB) interferences that are difficult to be filtered out, are specifically considered in both steadystate and dynamic signal modelling processes. Based on the models, a finite impulse response (FIR) filter and a Taylor–Fourier transform (TFT) filter that both contains notches at the OOB interferences are elaborately designed and implemented in two parallel channels. With the support of a super-resolution frequency analysis method for signal components, the filters’ central frequency and notches’ position can be online adaptively adjusted. Finally, the outputs are switched between the two-channel estimations according to a designed transient recognition strategy. Test results in IEEE standard scenarios and some other realistic but severe scenarios demonstrate the effectiveness and reliability of this proposed algorithm.

1. Introduction With the continuously growing penetration of distributed generations (DGs), electric vehicles (EVs) and other interactive loads, certain new features gradually emerge in distribution networks. Such features include intermittency and randomness of power sources, changeable operating modes of power grids, diversity and interactivity of loads, and potential islanding morphology, all of which complicate the operation and control of distribution networks. With the development of wide area measurements (WAMS) technology, phasor measurement units (PMUs) are widely used in bulk power systems, such as state estimation, dynamic monitoring, security assessment, stability control and other fields [1]. Theoretical research and application results have demonstrated its effectiveness and advancement. Compared to the traditional measurement methods, such as supervisory control and data acquisition (SCADA), PMUs have higher measurement accuracy, better dynamic performance, and unique global synchronization capability, which could help to address the new challenges in distribution networks. Typical applications based on the synchrophasor can be primarily divided into three categories, i.e., monitoring, diagnostic and

⁎

control, as discussed in [2–4], showing that the PMUs have good development prospects at the distribution level. It is well-known that the phasor estimation algorithm determines the precision and quality of PMU data under various measurement conditions, which is the cornerstone of all types of PMU-based technologies and applications. However, the previous algorithms designed for transmission network are not suitable for distribution network applications. Actually, measurement scenarios in distribution networks are more complicated. Firstly, distribution network lines are so short that the voltage phasor angle differences (tenths of a degree or less) are typically two orders of magnitude smaller than those in transmission networks, demanding the improvement of the angular resolution [2,5]. Even worse, there are more non-linear loads (e.g., rectifier, arc furnace, and frequency converter), distributed renewable generations (wind power, photovoltaic) and energy storages in distribution networks, which introduce a large number of interharmonics and strong noise [6–8]. Besides, distribution networks are closer to the customers than the transmission networks, resulting in more sophisticated dynamic behaviours [9,10], e.g. waveform oscillation, voltage flicker and frequency change. In summary, the signals in distribution networks have

Corresponding author. E-mail address: [email protected] (C. Lu).

https://doi.org/10.1016/j.ijepes.2019.105714 Received 15 October 2018; Received in revised form 21 October 2019; Accepted 18 November 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

more interferences and dynamic changes, which makes it difficult to estimate the phasor accurately. However, the short-line structure requires a higher phasor measurement precision at the same time. As we know, the IEEE std. C37.118.1-2011 [11] and its amendment [12] (uniformly referred as the IEEE stds. below) has defined some error limitations for the PMU algorithms that are classified into M-class and P-class algorithms, respectively corresponding to the steady-state and dynamic conditions. However, requirements of synchrophasor estimation in distribution network appear excessive to the IEEE stds., which means that higher accuracy is expected under the conditions where several interferences and dynamic behaviors may exist simultaneously. However, current researches about PMU algorithm still respectively focus on two aspects according to the IEEE stds.: measurement accuracy and dynamic performance. The discrete Fourier transform (DFT) is the most widely used PMU algorithm for its simple principle and harmonics suppression ability. However, inherent fence effect and spectrum aliasing will degrade its accuracy under some asynchronous sampling conditions. To improve the accuracy of traditional DFT based algorithms, several efforts have been utilized in recent years [13–15]. In [13], an adaptive samplinginterval method was utilized to reduce leakage error whenever frequency deviation occurs. But it is so complex that it is difficult to be implemented on the hardware. Ref. [14] adaptively composed the errors of the DFT algorithm based on a simple frequency tracking method in frequency deviation scenarios. However, the expressions of phase error and amplitude error are invalid when there are several interharmonics. In [15], a two-term window-based method was proposed to minimize the detrimental effects of the spectrum leakage due to offnominal frequency deviations. However, it cannot deal with low-order harmonic distortion or dynamic variation. In order to enhance dynamic-response ability of phasor estimating process, Ref. [16] employed Kalman filter (KF) to approximate the dynamic phasor, which illustrated instantaneous and synchronous attribute in case of smooth oscillations. However, its convergence speed is a little slow in transient changes. In [17], an extended complex Kalman filter (EKF) was designed to fast track the instantaneous frequency changes of the power system signals. Its stability and accuracy under multiple disturbances and fast variations should be further investigated. Some algorithms based on Taylor–Fourier Transform (TFT) and its improved forms that were presented in [18,19], also aimed to well track the signal dynamic behaviours. Ref. [18] used the Taylor approximation and weighted least squares (WLS) method based on classical windows to design a maximally flat differentiators with linear phase response, which can obtain accurate and effective estimates of dynamic phasor. In [19], the Taylor and WLS method were improved based on the frequency estimate returned by a classical interpolated DFT algorithm. In [20,21], recursive wavelet transform (RWT) was presented as an instantaneous and accurate tool to increase the phasor estimation efficiency. Among them, Ref. [20] proposed a real-time phasor estimation approach based on a newly constructed recursive RWT. Ref. [21] designed a special mother wavelet function, which can reduce computational complexity and response time compared to conventional RWT methods. Nevertheless, the stability of the RWT methods in more sharp disturbances should be further evaluated, just like the EKF methods. In summary, reliabilities of the above algorithms need to be further verified, especially under multiple interference conditions or fast change conditions, which may be more in accordance with the actual situation in distribution network. Besides, these methods are difficult to balance the steady-state accuracy and dynamic response requirements in distribution network. To improve the performance under multiple interference conditions, some researchers took harmonics or interharmonics into signal models to improve the anti-interference abilities of traditional TFT and EKF filters [22–26]. In [22], an extended TFT with complex envelope at each harmonic was employed to design a set of maximally flat filters for dynamic harmonic estimation. The research mainly focused on

amplitude estimation, while the accuracy of phase estimation was not mentioned. Similarly, Ref. [23] designed the TFT filter at each harmonic by fixed polynomials, which can avoid solving the whole pseudoinverse problem, thus reduces the computational burden. In [24], compressive sensing method and TFT method were utilized in a joint way to limit the impact of harmonic and interhamonic interferences. However, an iterative searching process was required to identify the interference components. In [25,26], the interharmonic was additionally detected and filtered out by adding a notch to the EKF filter. However, only one interharmonic could be identified and rejected, lacking to the resolution of the frequency scanning method. Although the above algorithms have shown their excellent abilities to supress interferences, their performance under the conditions with multiple harmonics, and especially with several interharmonics simultaneously, should be further investigated. As stated earlier, most of the current algorithms were generally dedicated to the steady-state or dynamic scenarios, but could not simultaneously comply with the requirements in both steady-state and dynamic conditions. To break this limitation, some researchers creatively proposed some cascaded or parallel filter structures [27–29]. In [27], on the premise of preserving the performance of the original TFT algorithm, an observation window length adaptively reduction approach based on fast changes detection was proposed to improve TVE response speed. However, the detection thresholds were heuristically set and did not consider the potential slow changes of interference level in distribution network, which may cause that the monitor indices exceed the thresholds even in steady state, resulting in invalid detection. Ref. [28] put forward two kinds of cascaded filters using for phasor estimation, respectively based on TickTock and Asymmetric structural variants. Entire and partial length of the filters are separately employed for P class and M class applications. The results of both filters were output, instead of choosing the results that matches the current scenario better. In [29], two different lengths and orders of TFT-WLS filters were integrated into a unified structure and their outputs were chosen according to a very simple transient detector. This algorithm can basically comply with the challenging requirements imposed for both P and M class PMUs. However, it doesn’t specially consider the harmful effects of out-of-band harmonics/interharmonics that are near to the fundamental component. Besides, the filters’ characteristics cannot be adjusted online to eliminate the unexpected interferences, due to the lack of real-time analysis of signal frequency components. In summary, these algorithms have not yet specially considered the detection and modelling of the harmonics and out-of-band (OOB) interferences in the filter design procedures. Besides, most of them are designed to meet the requirements of the IEEE stds. for transmission network. The above reasons lead to the fact that the multiple interferences rejection abilities of these algorithms cannot be fundamentally guaranteed. Based on the researches of PMU in transmission network, some pioneers have designed several kinds of PMU prototypes recently, aiming to monitor the frequency and synchrophasor at distribution network level [4,5,30–33]. In [4], a high-precision PMU based on hardware platform of power quality analyser, named Micro-PMU, was designed and applied in some scenarios of distribution network, while its PMU algorithm was not publicly introduced. However, the researches in [5,30,31] mainly focused on the low-cost hardware implementations. Ref. [32,33] developed a frequency record device (FDR) employed in residential voltage level. The algorithm utilized in the FDR was dedicated to improving the measurement accuracy under strong harmonic and noise distortion condition. This algorithm consists of six modules in series, resulting in long response time. Moreover, the demand of simultaneously dealing with multiple harmonic and OOB interferences still have not drawn enough attention in these PMUs for distribution network. In order to deal with the above problems, this paper proposes a comprehensive two-channel switching PMU algorithm based on the forefathers’ valuable research. Reasonable coordination and targeted 2

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

D channel: To improve the response speed and track the dynamic changes, a shorter filter with only two-cycle length in the D channel is employed to calculate the dynamic phasors, whose filtering properties are also available to be optimized online as the filter in S channel. Frequency analysis: To our knowledge, this is the first time to apply a frequency analysis method rather than prior knowledge to determine the frequencies of multiple unknown interharmonics in phasor estimating process. The results are utilized to tune the filters’ characteristics in the two channels online. Transient recognition: A transient recognition strategy is developed to detect the fast changes, such as step change, sudden modulation and unexpected harmonic/interharmonic appearance. The recognition results of the signal states are used for guiding the output switching activities. Based on the framework presented, the proposed PMU algorithm is attempted to excavate its potential for complying with the specified scenario requirements at the distribution level and achieving a balance between measurement accuracy and dynamic response. Its details are described in the following sections.

improvement of the previously presented methodologies in [24,27–29] have been made to design this comprehensive phasor estimation algorithm, which can be effectively utilized under multiple interferences and variations conditions in distribution network. The innovations of this algorithm are mainly embodied in the following three aspects: (a) Multiple harmonics/ interharmonics that are close to the fundamental wave are specially considered in both steady-state and dynamic signal modelling, which theoretically ensures that the designed filters in both channels can effectively remove the interferences’ influence on the fundamental phasor estimation; (b) The Cosh window that has narrow main lobe and low side lobe is firstly employed to improve traditional FIR and TFT filters for phasor estimation. Besides, some notches at the identified harmonics/interharmonics are added to the filters’ transition band, which is intended to further suppress all of these interferences, not only one interharmonic. (c) The signal frequency components including harmonics and interharmonics are analysed online based on a super-resolution spectrum estimation method, fourth-order ESPRIT, the results of which are utilized to adjust filters’ central frequencies and notches’ position, ensuring that the unexpected interferences can be identified and filtered out adaptively.

3. Signal modeling and filter design The waveform of interest for most studies in transmission networks is generally a clean sinusoidal signal that corresponds to the phasor definition. However, as mentioned before, the signals at distribution level generally contain many harmonics/interharmonics, although most of the harmonics/interharmonics could be filtered out with appropriate filters, and those close to the fundamental frequency in the passband or transition band of the filters are not efficiently suppressed. The interference over a range from below the passband to the second harmonic is specially called the out-of-band (OOB) interference in the IEEE stds. However, note that the interferences considered in the signal models here are not limited to the definitions in the IEEE stds. The interharmonics and some low-order harmonics (e.g., 2nd harmonic) that probably cause larger estimation errors are all included. The choice of these interferences depends on the amplitude-frequency response characteristics of the filters. Further, in contrast to the models in [25,26] that only can filter out one interharmonic component, several interharmonics are taken into account in this paper.

The rest of this paper is organized as follows. An overall framework of the proposed algorithm is shown in Section 2. Section 3 presents the signal modelling and adaptive filters. Two assistant functions, namely, transient recognition and frequency analysis, are introduced in Section 4. By integrating those techniques, the algorithm implementation is given and subsequently verified based on a series of tests in Section 5. Comparative studies and relevant discussions are also included. The study’s conclusions and several possible future research directions are presented in Section 6. 2. Framework of the proposed algorithm The proposed PMU algorithm is mainly composed of two parallel phasor estimating channels, namely, steady-state channel (S channel) and dynamic channel (D channel), which are respectively utilized to process the steady-state and dynamic signals. The framework of this algorithm is shown in Fig. 1. Except for the two channels, there are another two assistant modules in this algorithm: frequency analysis and transient recognition. The functions of the four modules are roughly introduced as follows, while the technical details involved in these modules will be explained in detail in Section 3 and Section 4. S channel: The steady-state phasors are estimated based on a fourcycle filter implemented in this channel. The central frequency of this filter is adjusted in real time, and some notches corresponding to the interferences near the fundamental component are adaptively added to the filter.

3.1. Signal modelling 3.1.1. Steady-state signal model In steady state, large frequency deviation may be mainly caused by the local power imbalance in distribution networks, such as in the electrical island mode, which is always assumed to be constant. Thus, the steady-state signal model is established as

xs (k ) = Ac cos(2 fc kT + +

c)

Aih cos(2 fih kT +

+

Am cos(2 mfc kT +

ih ) +

m)

Ab cos(2 fb kT +

b)

+ xNoise (k ) (1)

where

• x (k) represents the steady-state signal, superimposed with Gaussian s

• • • Fig. 1. Framework of the proposed algorithm. 3

noise (i.e., xNoise(k)), harmonics, interharmonics and the OOB components that involves the interferences of special interest separately. The fixed sampling interval is T. k is the sampling index. fc is the fundamental frequency, which equals the normal frequency f0 plus a deviation △f; Am is the amplitude of the harmonic, and φm is the phase. m is the harmonic order of this signal. When m is c, the harmonic is the fundamental frequency component. Aih, φih and fih correspond to the amplitude, phase and frequency values of the interharmonic respectively, with the subscript ih representing the interharmonic order.

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

• A ,φ

b b and fb indicate the possible OOB waves’ amplitude, phase and frequency, respectively, with b denoting the OOB order.

3.1.2. Dynamic signal model Under dynamic conditions, amplitudes and phases of the measurement signals are no longer constant but probably vary linearly or nonlinearly over time. However, special attention should be paid to the dynamic behaviours of the fundamental and the OOB component because the former is used to estimate the phasor, and the latter has a significant effect on the phasor estimation. To approximate the variations, the two components’ amplitudes and phase angles are linearly expanded at the middle point of the data window, utilizing a secondorder Taylor series. This model is shown as

xd (k ) = Ac (k )cos(2 fc kT + +

c (k ))

Aih cos(2 fih kT +

ih )

+

Ab (k )cos(2 fb kT +

+ m 1

Am cos(2 mfc kT +

b (k ))

m)

Fig. 2. Amplitude-frequency characteristics of these window functions.

+ xNoise (k ) (2)

filters because of its stability and linear-phase properties. Considering that the shorter reporting delay is recommended in the IEEE stds., a shorter filter length is preferred. However, a shorter filter length will result in a wider main lobe, wrapping more OOB components in. Obviously, this approach represents another difficult trade-off. In this paper, a four-cycle FIR filter is chosen, and the time delay is to be compensated. Although the Cosh window has shown a better performance than the other windows, there is still no guarantee that the OOB components can be effectively supressed. Thus, these components are paid special attention in the filter design process. To design a FIR filter, an attenuation threshold (e.g., 100 dB) is first determined to ensure that most of the interferences can be effectively suppressed. Next, the components, whose attenuation levels are below this threshold according to the amplitude-frequency response curve of the Cosh window in Fig. 2, are all treated as the OOB interferences here. As described in the steady-state signal model, the signal samples can be written as a linear combination of fundamental and OOB components at the middle position of the data-window, omitting the components whose attenuation are over the threshold.

where

• A (k) = a •

2 2 c 0 + a1kT + a2(kT) + o[(kT) ] is the expansion of the fundamental amplitude. Similarly, the phase is expended as φc(k) = b0 + b1kT + b2(kT)2 + o[(kT)2], where o[(kT)2] is the Peano’s remainder. Ab(k) = ab0 + ab1kT + ab2(kT)2 + o[(kT)2] is the OOB amplitude, and φb(k) = bb0 + bb1kT + bb2(kT)2 + o[(kT)2] refers to the OOB phase angle.

Note that the frequencies of harmonics, interharmonics and the specially defined OOB components in the steady-state and dynamic models are estimated based on a carefully selected ESPRIT based frequency analysis method that has shown a super-resolution ability. The procedures of method selection and validation will be described in the Section 4. 3.2. Adaptive filters design

xs ( - KT )

Based on the signal models above, filters with linear-phase property have been designed in this part. The window function based filter deign strategy utilized here is simple and practical to cope with the drawbacks of the DFT-based algorithms. However, the window function should be carefully selected because there is an unavoidable trade-off between the main-lobe width and side-lobe attenuation. To find a compromise solution, the window function with relatively narrow main lobe and low side lobe becomes a better choice. An adjustable window function named three-parameter Cosh window recently proposed in [34] has shown a better performance than the others:

cosh ac 1 wcosh (n) =

( )

2n 2 N 1

cosh(ac )

Xs =

xs (KT )

=

c KT

e

- j bi KT

1 ej

ej

1

c KT

ej

bi KT

bi KT

1 e-j

bi KT

ej

c KT

1 e-j

c KT

pc pbi p¯bi p¯c

m

, |n|

xs (0)

e-j

= Cs ps (4)

where is the fundamental angular frequency, and c = 2 fc bi = 2 fbi (i = 1, 2, … Nb) is the angular frequency of each OOB component. The fundamental component can be written in complex Ac cos(2 fc kT + c ) = 0.5A c e j c e j2 fc kT + 0.5Ac exponential form: e j c e j2 fc kT (k = −K, −K + 1, … K). Similarly, the OOB components Abi cos(2 fbi kT + bi) = can also be express as pc = 0.5Abi e j bi e j2 fbi kT + 0.5Abi e j bi e j2 fbi kT .In these forms, 0.5Ac e j n represents the steady-state fundamental phasor, and pbi = 0.5Abi e j bi represents the OOB component. p¯c and p¯bi are their conjugate values. Xs is the vector composed by a series of signal samples. ps is the phasor vector. Cs is the coefficient matrix whose column vector represents the Fourier basis corresponding to the fundamental or OOB frequency. Similar to Ref. [18,19,29],we also adopts a window weighted least square (WLS) method to calculate the phasors in the vector ps, which can be directly derived as

K (3)

where K = (N − 1)/2 is the half of data-window length. The cosh means hyperbolic cosine function. ac and ρm are the adjustable parameters. N is the order of the window. The frequency response characteristic of the well-adjusted Cosh window (ac = 7.0, ρm = 1.2) is compared with several common windows in Fig. 2. Clearly, the Cosh window has lower side lobe than the others when they have the same lengths and the 3-dB bandwidths. For this advantage, this well-adjusted Cosh window is applied in the following filters design. 3.2.1. Steady-state FIR filter It is well-known that the FIR filter is one of the most popular digital 4

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

estimate the fundamental phasor and its derivatives at the reference time

(5)

ps = (CHs WHWCs) 1Cs W· X s = Fs·X s

where W represents the weight matrix, whose diagonal elements are composed of the values of the four-cycle well-adjusted Cosh window. Fs is the product of the pseudoinverse matrix, coefficient matrix and weight matrix. The symbol H represents the Hermite operator. Note that the OOB phasors in the vector ps will be calculated and will be attenuated to some extent because of the window function characteristic in the Fig. 2. However, our efforts primarily focus on the fundamental phasor. Thus, the first row related to the fundamental phasor estimation in the matrix Fs is chosen as the coefficient vector of the FIR filter.

In particular, the TFT filter has a unique advantage that the fundamental frequency and its rate of change of frequency (ROCOF) are involved in the phasor derivatives and can be derived simultaneously as

Ac = a 0 = 2 |pc(0) | c

3.2.2. Dynamic TFT filter When the signal is under fast dynamic change condition, its amplitudes or the phases can no longer be assumed to be constant, even though the window isn’t too long. Therefore, unlike the steady-state filter, the variations should also be considered in the filter design. Similarly, the OOB components are also included. Based on the dynamic signal model, the TFT method is utilized to improve the filter as

xd ( KT ) Xd =

c (- K , c )

c (- K , bi )

c¯ (- K , bi )

c¯ (- K , c )

c (0, c )

c (0, bi )

c¯ (0, bi )

c¯ (0, c )

c (K , c )

c (K , bi )

c¯ (K , bi )

c¯ (K , c )

=

xd (0) xd (KT )

fc = f0 +

ROCOF =

= Cd pd (7)

where Xd is the vector of dynamic signal samples. pd is the vector of the phasors and their derivatives, whose elements are

pc = pc(0) , pc(1) , pbi =

pbi(0) ,

pbi(1) ,

1 (2) p 2 c

T

1 (2) p 2 bi

, p¯c = p¯ (0) c , p¯ (1) c , T

, p¯c =

p¯ (0)

bi ,

p¯ (1)

1 (2) p¯ c 2

(8)

c (k, c ) = [e j c¯ (k, c ) = [e

c kT ,kTe j c kT ,(kT ) 2e j c kT ],

- j c kT ,kTe - j c kT ,(kT ) 2e - j c kT ]

(10)

c (k, bi) = [e j bi kT ,kTe j bi kT ,(kT ) 2e j bi kT ], c¯ (k , bi) = [e - j bi kT ,kTe - j bi kT ,(kT )2e - j bi kT ] k = - K , - K + 1,

, 0,

K - 1, K

(11) The vector pd can also be derived using the WLS method in a similar form

pd = (CHd WHWCd) 1Cd W ·X d = Fd· X d

jb0}

jb0 }

(17)

+ 2a 0 ·b12

Im{pc (2) e jb0} b2 = 2 2 a0

(18)

a1·b1 (19)

Obviously, the numbers of the harmonics and interharmonics contained in the signal models and their frequencies are generally not a priori information. Most of the research efforts [22–26] always assumed that the frequency composition of the signal had already been known before the filter design. However, such assumptions may lead to a contradiction: On the one hand, if the frequency component actually contained in the signal is omitted, the phasor estimation error will increase; On the other hand, including too many frequency components that do not exist will increase the computational burden and even cause numerical stability problems. Therefore, a reliable signal frequency analysis method is urgently needed. Recently, the modern spectrum estimation (MSE) based methods, such as multiple signal classification (MUSIC) and estimating signal parameters via rotational invariance technique (ESPRIT), have shown their advantages of super resolution and good stability in frequency analysis [35,36]. Compared to the traditional Fast Fourier Transform (FFT) based method, the MSE based method (MUSIC) has stronger ability to distinguish the interharmonics that are close to harmonics (as shown in Fig. 3). Compared with MUSIC method, the ESPRIT method can get the signal frequency components by calculating the generalized eigenvalues of the matrix, rather than scanning the whole frequency domain. Therefore, its computational efficiency is greatly improved. In this paper, a fourth-order autocorrelation function of the signal is used according to Ref. [37].

(9)

pc(2) are

pc(1)

2 Im{pc(1) e b1 = f0 + 2 2 a0

(16)

4.1. Frequency analysis

T

where and the zero-order, first-order and second-order derivatives of the fundamental phasor, respectively. pbi(0) , pbi(1) and pbi(2) are the derivatives of the OOB phasors. p¯ (0) c denotes the conjugate value of pc(0) . p¯c (1) and p¯ (2) c are the conjugate values of pc(1) and pc(2) , respectively. Besides, p¯bi(0) , p¯bi(1) and p¯bi(2) are the conjugate values of the different derivatives of the OOB phasor. T is the transpose operator. The matrix Cd consists of the Taylor series basis of different frequency components at each sample

pc(0) ,

jb0}

4. Frequency analysis and transient recognition

T

1 (2) p¯ bi bi , 2

(15)

where Ac and c are the amplitude and phase angle of fundamental phasor, respectively. a0 , a1and a2 are the zero-order, first-order and second-order derivatives of the amplitude.b0 , b1 and b2 are the zeroorder, first-order and second-order derivatives of the phase angle. f0 is the fundamental frequency rating. fc is its estimated value. ROCOF is the rate of change of the fundamental frequency. Re{} represents the operation that the real part of the complex number is selected. Im{} refers to selecting the imaginary part. By considering the OOB interferences in the two filters designing procedures, notches related to these components are adaptively placed on the filters’ passband and transition band to better eliminate their influences, which will be introduced concretely in the next section. To reduce the computational burden, the filter coefficients under most conditions can be calculated offline and conveniently retrieved online.

pbi

p¯ 1

(14)

pc(0)

a2 = Re{pc(2) e

p1

p¯ bi

= b0 =

a1 = 2 Re{pc(1) e

(6)

Hs = Fs (1, :)

(13)

Hd = Fd (1 : 3, :)

(12)

where W is also the diagonal matrix with a two-cycle Cosh window, and Fd is the product of the pseudoinverse matrix, coefficient matrix and weight matrix. The first three rows of Fd can be used as the filter to 5

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

Uniform Phase angle time stamp 4-cycle FIR compensation window 20ms

Raw phasor

20ms

Latest sampling point Reporting latency

2-cycle TFT window Fig. 7. The schematic diagram of phasor compensation in the S channel.

Fig. 3. Spectral analysis results of the MSE based and FFT based methods.

Fig. 8. Verification results of the frequency analysis function.

R 4x ( ) =

xs ( ) xs (t

)]2

)[x s (t

t=

Then the fourth-order autocorrelation values are employed to form the Hankel Matrix Y:

Fig. 4. Frequency analysis errors of the three MSE based methods.

Maximum Frequency Error (Hz)

0.45

2-cycle data 3-cycle data 4-cycle data 6-cycle data

0.4

0.35 0.3

Y=

0.4

0.2 0.13

0.05 0

0.02

0 0 0

80

0.01 0 0

0.02 0.02 0.01

70 60 Noise Level (dB)

R 4x (M 1) R 4x (M ) R 4x ( L + M

1)

(21)

The detail process is introduced in Ref. [37]. The fourth-order autocorrelation function is used instead of the measured signal to suppress the influence of Gaussian noise in the frequency detection accuracy. Three common MSE based methods, Root-MUSIC, traditional ESPRIT (TR-ESPRIT) and the fourth-order ESPRIT (FO-ESPRIT) based on sixcycle simulated data, are simply compared here. The test signal is also simulated as the example in the Fig. 3. The average frequency analysis errors of these three methods under different levels of Gaussian noise are shown in Fig. 4. After comparison and analysis, the FO-ESPRIT method is found to have better accuracy than the other methods. Therefore, it is employed in the frequency analysis module in Fig. 1 to overcome the quantity limitation of interharmonic detection due to the interpolation FFT method in the Refs. [26,27]. As shown in the Fig. 3, the FO-ESPRIT that is a kind of MSE method. Thus, it is less sensitive to the length of analysis window than FFT-based methods. Therefore, a shorter window that is less than six cycles seems to be an alternative. However, considering the influence of noise, the window length cannot be too short. Because the fourth-order autocorrelation quantity is estimated based on the finite sampling data. Too short window will result in degradation of noise suppression capability. The frequency analysis errors based on four different lengths of windows under several noise levels are shown in Fig. 5. Therefore, a 6-cycle FO-ESPRIT method is employed for online frequency analysis in this paper, as a trade-off between the response time and the estimation accuracy under noisy conditions.

0.25

0.1

R 4x (0) R 4x (1) R 4x (L)

0.27

0.15

(20)

0.09 0.07 0.02

50

Fig. 5. The frequency analysis errors using different lengths of data under several noise levels.

4.2. Transient recognition

Fig. 6. Illustration of the algorithm implementation procedure.

The phasors estimated based on the FIR filter are committed to 6

Electrical Power and Energy Systems 117 (2020) 105714

2.5

0.3

2

0.25 Probability (%)

Time error (ms)

Y. Wang, et al.

1.5 1 0.5 0

0.2 0.15 0.1 0.05

0

200

400

600

800

0

1000

0

0.5

1

1.5

2

2.5

3

Number of the Monte Carlo test

Time error (ms)

(a) Time errors of the Monte Carlo test

(b) Probability distribution of the time errors

Fig. 9. The time errors and their probability distributions of the step change identification.

data window, then there will also be an obvious difference between the errors based on the first half and the second half of the data window. Thus, two other indices based on the above analysis are defined to further improve the reliability of the transient recognition:

Table 1 TVE compliance tests under steady-state conditions. Test condition

Frequency range Harmonic distortion Out-of-band interference*

Range

± 5.0 Hz 10%, each harmonic up to 50th 10% of amplitude

TVE limit (%)

Max TVE (%) RWT

EKF

TCS

1.0 1.0

1.4 × 10−1 5.9 × 10−2

1.2 × 10−1 1.7 × 10−9

3.4 × 10−6 1.6 × 10−7

1.3

−1

1.1 × 10

7.2 × 10

−2

Indice3 =

k= K K

Indice4 =

−7

3.2 × 10

having higher precision under normal conditions. However, the result obtained from the TFT filter will be a better choice to realize fast dynamic response when there are some abrupt changes. These changes include the transient processes that cause the signal to enter dynamic state or a new steady state from the current steady state, such as step changes, modulation or linear variations on signal amplitude/phase/ frequency, etc. Therefore, a recognition strategy is required to detect the occurrence time of these changes and guide the selection of appropriate outputs. Several indices based on experience are set to comprehensively identify the moment when the signal state suddenly changes. Generally, the amplitude and phase angle will become changeable rather than constant, when the signal enters a certain dynamic state from the original steady state, such as amplitude/phase /frequency modulations. Therefore, the first-order derivatives of the amplitude and phase from the dynamic phasor estimation procedure are used as parts of the transient recognition criteria

{Indice1 >

| x (k )

0

x (k )|2

k= 0

x (k )|2

| x (k )

k= K

(23)

(K + 1)·|A1 |

th1}

{Indice2 >

th2}

{Indice3 >

th3}

{Indice4 >

th4}

(24)

Transition occurs

where σth1, σth2, γth2 and γth2 are the thresholds of these indices. Generally, the harmonics, interharmonics or noise level probably vary slowly during a short-term. The levels of the indices in Formula (22) and (23) are almost assumed to keep unchanged. However, these parameters may change significantly after a long-time accumulation. This phenomenon is very likely to happen in distribution networks, which means that if the thresholds are configured as constant all the time, the detection criterion in Formula (24) has a risk of failure. Therefore, the thresholds should be adaptively updated in a regular time or when the signal re-enters a new steady state after some dynamic events: thi

=

i·max{avg[Indicei (1:

N )], … , avg[Indicei (9N + 1:10N )] }

i = 1, 2, 3 or 4

(22)

Indice2 = |b1|

|x (k ) x (k )|2 (2K + 1)·|A1 |

where x(k) is the actual value at sampling point k, and x (k ) is the corresponding predicted value. Finally, the four indices are synthetically employed in the detection criterion of the transient events

* For the 50 Hz system, the range of out-of-band frequency is from 10 Hz to 100 Hz.

Indice1 = |a1|,

K

(25)

where some average values of each index at serval steady-state periods are calculated, then the maximal one is multiplied with a margin factor (i.e., αi) to determine its threshold. The thresholds are applied only to detect the initial moment of the dynamic behaviours. When the indices satisfy the criterion in Formula (24), the D channel will be selected and locked to output, until the indices are all less their thresholds.

If the signal frequency and the sampling interval are known, the signal can be reconstructed based on the phasor estimation. The steadystate estimation results are utilized for their better anti-interference ability. The error between the real values and the reconstructed values is lower than a setting threshold in most cases, except that there are some sudden changes. Additionally, if a sudden change occurs in the

Table 2 FE and RFE compliance tests under steady-state conditions (Range settings are the same as Table 1). Test condition

Frequency range Harmonic distortion Out-of-band interference

FE(mHz)/RFE (Hz/s) limit

5.0/0.1 25/None 10/None

Max FE (mHz)/RFE (Hz/s) RWT

EKF

TCS

28/0.2 2.3/0.1 62/0.3

2.8/0.06 4.8 × 10−10/3.5 × 10−11 0.3/0.1

2.5 × 10−4/6.0 × 10−6 1.3 × 10−6/9.4 × 10−9 1.4 × 10−5/3.3 × 10−7

7

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

Table 3 Phasor performances under dynamic conditions. Test condition

Range

Modulation

TVE limit (%)

0.1 × amplitude, 0.1 to 5 Hz 0.1 rad, 0.1 to 5 Hz Rate: ± 1 Hz/s, Range: ± 2 Hz

Frequency ramp

TVE response time (ms)/Max overshoot or undershoot Step changes Amplitude ± 10% Angle ± 10°

Max TVE (%) RWT

EKF

TCS

3.0 3.0 1.0

0.2 2.8 0.7

0.6 1.4 0.3

1.0 × 10−2 1.0 × 10−2 3.9 × 10−5

40/5.0% 40/5.0%

38/0.5% 74/0.2%

18/0.4% 163/11%

10/2.0% 10/3.4%

Table 4 Frequency and ROCOF performances under dynamic conditions (Range settings are the same as Table 3). Test condition

FE(mHz)/RFE (Hz/s) limit

Modulation

0.1 × amplitude, 0.1 to 5 Hz 0.1 rad, 0.1 to 5 Hz Rate: ± 1 Hz/s, Range: ± 2 Hz

Frequency ramp

Frequency/ROCOF response time (ms) Step changes Amplitude ± 10% Angle ± 10°

Max FE (mHz)/RFE (Hz/s) RWT

EKF

TCS

60/2.3 60/2.3 10/0.4

0.7/0.1 20/2.6 1.4/0.9

16/2.3 40/4.4 17/3.2

2.0/1.3 7.3/0.8 6.3 × 10−3/4.5 × 10−3

90/120 90/120

83/90 146/152

88/103 98/124

29/30 18/19

Fig. 11. Property of the FIR filter under the multiple-interference condition. Fig. 10. Impacts of noise on these three PMU methods.

changes, and then achieving a new steady state become more important. Therefore, the higher-order and shorter-window TFT filter is required. Based on the above considerations, the two designed filters are integrated in the framework in Fig. 1, together with the frequency analysis method and transient recognition strategy, forming the comprehensive phasor estimation algorithm. Its whole implementation procedure is illustrated in Fig. 6. First, six-cycle signal samples are down-sampled by two for the frequency analysis based on the FO-ESPRIT method. This module is implemented in a separate thread and update the results every 10 cycles (i.e. 200 ms), which means that the signal frequency components can be calculated independently, avoiding task conflict with the phasor estimation. When there is an abrupt change detected, the frequency analysis method is shut down until the signal return to steady state. This design mechanism is set to keep the frequency analysis from making mistakes in sudden change or dynamic scenarios. Next, the fundamental frequency together with the OOB frequencies are fed back to adaptively adjust the filters in S channel and D channel. Moreover, the amplitudes of interferences, such as out-of-band harmonics and interharmonics, are calculated by WLS method. This process is performed every 10 cycles, except when a transient event occurs. During the transient event period, the results of signal frequency

Table 5 Comparisons of the three methods under the sever condition. Methods

Max TVE (%)

Max FE (mHZ)

Max RFE (Hz/s)

RWT EKF TCS

1.03 4.73 0.05

13.2 55.9 1.45

0.15 39.8 0.13

5. Implementation and Verification 5.1. Algorithm implementation As mentioned above, the steady-state filter is designed based on assumption that the amplitude and phase of synchrophasor are unchanged in the 4-cycle window. It is mainly designed to estimate the synchrophasor in cases that the signal is stationary or slowly varies, which are the most common scenes in distribution network. In these scenes, accurate estimation of phasors is the primary task. However, if the state of the signal suddenly changes, for example sudden dynamics or unexpected interference occurrence, responding quickly to the 8

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

New steady state

interharmonic arrives 0.7

0.8

0.9

1

OOB Flag

1

0

1.1 0 1

0.5

1.2

1.3

0.8

0.9

1.4

1.5

Without OOB With OOB

Lock frequency analysis module 0.7

State flag

0.5 0

2

0 Steady state 1 Dynamic/Transient

Frequency analysis flag

State Flag

1

Identify OOB components 1

1.1

1.2

1.3

1.4

1.5

0

Phase error (°)

TVE %

0.01

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

0.01

Frequency error (Hz)

TVE (%)

0.02

0 -0.01 0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

Fig. 12. The processes of an unexpected interharmonic occurrence.

0 0.6 2

Dynamic arrives 0.8

modulation 1

1.4 1.5 1.6

0

-50

-50

-100

-100

-150

Filter s characteristic before 1.3s

-250

-300

1

Lock frequency analysis module

0.5

First analysis after this dynamic

0.8

1

1.2

1.4

1.6

1.8

0.8

1

1.2

1.4

1.6

1.8

0.8

1

1.2

1.4

1.6

1.8

0.8

1

1.2 Time (s)

1.4

1.6

1.8

0.2 0.1

0.02 0 -0.02 -0.04 0.6 0.4 0.2 0

-0.2 0.6

-150

Filter s characteristic after 1.3s

-200 -250 -300

-350 -200 -150 -100 -50

1.8

0 Stop frequency analysis 1 Start frequency analysis

1.5

0 0.6 0.3

1.2

From the Fig. 7, we can see that the distance between the raw phasor and the uniform time stamp is only 20 ms. As mentioned before, the amplitudes of the phasors estimated by the FIR filter can be

0

-200

0.5

Fig. 14. The algorithm response processes when an unexpected amplitude modulation occurs.

Attenuation (dB)

Attenuation (dB)

components remain unchanged until the signal returns to steady state again. Then, the dynamic and steady-state phasors are estimated by the 2cycle TFT filter and 4-cycle FIR filter respectively, at a report rate of 100 frames/s. The time stamp is set at the middle point of the TFT filter window (reference time for Taylor expansion), which is one cycle before the latest available sampling point. It is worth noting that the raw phasor directly estimated by the FIR filter in S channel, lags behind the dynamic phasor by one cycle, because it corresponds to the middle point of the 4-cycle filter window. In order to synchronize the two estimation phasors, the raw phasor should be compensated to the uniform time stamp. This process is shown in Fig. 7.

1

0 0.6 0.04

Rocof error (Hz/s)

Time (s)

0 Steady state 1 Dynamic/Transient Dynamic leaves Amplitude

1.5

-125

0 50 100 150 200 250 300 Frequency (Hz)

-350 -200 -150 -100 -50

125

0 50 100 150 200 250 300 Frequency (Hz)

Fig. 13. The amplitude-frequency characteristics of S-channel filter before and after 1.3 s. 9

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

Fig. 15. Signal dynamic behavior sequence diagram.

a 1-cycle low-pass filter. Finally, the transition recognition strategy designed in the Section 4 is employed to detect the sudden dynamic and transient events every 4 points (i.e. a time interval of 1.25 ms) online. A state flag is set to control the frequency analysis process and switch the phasor estimation outputs between the two available channels. The thresholds of the indices are adaptively updated in steady-state condition as in Formula (25). 5.2. Tests and verifications 5.2.1. Frequency analysis function verification As aforementioned, frequency analysis is a necessary part of the proposed algorithm framework. Furthermore, it can also play a role in power quality analysis in distribution networks, providing the frequencies and amplitudes of the harmonics and interharmonics every 0.2 s. This ability can meet the requirement in IEC std. 61000-4-7 [38]. To verify it, a test signal with multiple frequency components and a 60dB Gaussian noise in Formula (27) was tested in Fig. 8. The frequency analysis results are updated with a 200 ms period. The testing results clearly demonstrate that the frequency analysis method can distinguish the harmonics and interharmonics and accurately estimate their amplitudes, although they are very close to each other, for example the harmonic that is 250 Hz and the interharmonic that is 287.5 Hz in the time range from 0 to 10 s. Besides, this figure also shows this method’s online frequency components tracking capacity, when the signal components’ frequencies and amplitudes change at 10 s.

Fig. 16. Results of the test signal based on the TCS algorithm. Table 6 Profiling algorithm’s one execution cycle on ThinkPad T430/64- bit Windows 7 Intel I5-3230 + 2.60Ghz platform. Calculating Thread 1 Phasor estimation Transient Detection Output selection Total time

1.0 cos(50 × 2 t + /3) + 0.1 cos(150 × 2 t ) + 0.06 cos(200 × 2 t ) + 0.1 cos(250 × 2 t ) + 0.12cos(287.5 × 2 t) + 0.08 cos(350 × 2 t ), t < 10s x (t ) = 0.9 cos(49 × 2 t + /3) + 0.12 cos(147 × 2 t ) + 0.08 cos(196 × 2 t ) + 0.08 cos(245 × 2 t ) + 0.1cos(281.75 × 2 t) + 0.08 cos(343 × 2 t ), t 10s

Calculating Thread 2 61.84us 16.97us 3.95us 82.76us

Frequency analysis & coefficients update

2.14 ms

2.14 ms

5.2.2. Transient recognition strategy verification According to the definitions of indices defined in Formula (22) and Formula (23), they may be influenced by the positions of the potential abrupt changes in the data window. Therefore, an abrupt change such as amplitude step that randomly occurs in different positions in the data window is simulated to show the effectiveness of the transient recognition strategy based on the designed strategy. Monte Carlo method is used to randomly generate 1000 locations where the amplitude step change may occur. These positions obey uniform distribution. The errors between the identification value and the true value of the step change occurrence time and their distribution are shown in Fig. 9. From this figure, it can be seen that the time errors of the step change identification in about 82.0% test cases are within 1.0 ms. Besides the maximum error does not exceed 2.5 ms, which is two times of the intervals of online transient recognition. These results demonstrate that the designed transient recognition strategy can correctively

approximated to be unchanged within this 20 ms in steady-state scenarios. Therefore, it is only necessary to compensate the phase of the raw phasor to the time stamp position, which is performed by a common operation Stamp

=

Raw

+

=

Raw

+ 2 f· t

(27)

(26)

where Stamp is the phase angle at the time stamp, Raw is the phase angle of the raw phasor estimated by the FIR filter. t is equal to 20 ms. f is the fundamental frequency, which is assumed to be constant within the time interval of compensation. Besides, it is clearly seen that the reporting latency of the FIR filter at steady state is 20 ms. The same is true for the reporting latency of TFT filter. Moreover, the frequency and ROCOF in the D channel can be obtained together with the dynamic phasor as shown for Formula (14) to Formula (19). Differently, the frequency and ROCOF in the S channel are derived based on the steady-state phasor angle and then filtered by 10

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

locate this step change, no matter where it happens.

however, its errors under lower noise levels are larger than the other two algorithms. After comprehensive comparison, it is concluded that the proposed TCS algorithm is better at noise rejection and can make the TVE, FE and RFE sufficiently small over most noise levels.

5.2.3. Tests under the IEEE std. Conditions In this part, the proposed Two-Channel Switching PMU algorithm (TCS) was tested under the specified conditions in the IEEE stds. based on the MATLAB. Moreover, two types of representative algorithms, namely, the EKF based algorithm [26] and the RWT transform based algorithm [12], were also evaluated for comparative studies. The sampling rate is set as 6400 Hz, and the reporting rate is 100 frames/s for the 50 Hz system. In addition, the frequency and ROCOF estimation were also analysed here for their potential applications in the master station, such as fault diagnosis and frequency monitoring. Under steady-state conditions, the testing requirements for M-class PMU in the IEEE stds. are preferred because they are much stricter than the requirements for P-class PMU. For example, a wider frequency range of 50 ± 5.0 Hz rather than 50 ± 2.0 Hz is necessary to adequately evaluate the performance of the proposed algorithm, especially in some unexpected scenarios at distribution level, such as an isolated island of electricity. Therefore, all the M-class PMU testing conditions in the IEEE stds are employed here, and the results of total vector error (TVE) are listed in Table 1. Besides, the frequency error (FE) and ROCOF error (RFE) are listed in Table 2. The results that exceed the limits are written in bold form. According to Table 1, the TVEs of the three algorithms are very small under these steady-state conditions, which are considerably lower than the limits for the M-class algorithm in the IEEE stds. The TCS algorithm has similar TVEs under all three conditions, because the FIR filter’s central frequency and stopband can be adaptively adjusted according the given signal’s property. Table 2 shows that the required specifications of FE and RFE can only be achieved by the EKF and the TCS. Besides, the TCS performances are even better when the signal suffers from frequency deviations or OOB interferences. Therefore, we can conclude that the TCS has a better capability of frequency offset correction and OOB interference suppression. Generally, the dynamic behaviors of the signal in electric power system are mainly classified into three categories: amplitude/phase modulation, frequency ramp and amplitude/phase step. Their mathematical representations are expressly defined in the IEEE stds. These dynamic conditions are simulated to evaluate the three algorithms’ compliance with the P-class requirements under dynamic conditions. The results are listed in Tables 3 and 4, where these undesirable ones are also marked with bold. It is seen from the Table 3 that the TVEs of the TCS algorithm under modulation and frequency ramp conditions are much less than the TVEs of the other two algorithms. The TVE response times of this algorithm under both amplitude/phase step change conditions are about 10 ms. As shown in Table 4, the frequency errors and ROCOF errors of this algorithm can totally meet the critical regulations. The frequency and ROCOF response speeds of the TCS are algorithm also superior to the others. The results under dynamic conditions further demonstrate the proposed algorithm’s high accuracy and fast response.

5.2.5. Multiple interference tests As mentioned above, the signals at the distribution level may contain several harmonics, interharmonics and a notable level of noise. Even worse, a large frequency deviation may also occur at the same time that is more severe but realistic than the test conditions in the IEEE stds. To verify the algorithms under this situation, we constructed a signal in Formula (28). The estimation errors of phasor, frequency and ROCOF are presented in Table 5.

x (t ) = cos(2 ft ) + 0.06 cos(0.5 × 2 ft + /4) + 0.1 cos(3 × 2 ft ) + 0.06 cos(3.25 × 2 ft ) + 0.05 cos(5 × 2 ft + /9) + 0.05 cos(7 × 2 ft ) + xNoise 60dB (t ), f = 45Hz

(28)

It is obvious that the multiple-interference signal poses a challenge for the PMU algorithms, especially for the EKF algorithm, which only has the ability to address one interharmonic. However, the TCS algorithm is proven to be more efficient in processing the complicated signals because it benefits from its adaptive filters, which are adjusted according to the online measurements of the fundamental and OOB frequencies. This property is revealed in Fig. 11, where the filter is added with some notches at the OOB frequencies and their conjugate frequencies, consequently improving its interference rejection ability. 5.2.6. Unexpected interference/dynamic occurrence tests To verify the coordination of the several modules in the proposed algorithm, some tests under conditions that an unexpected interference or an unexpected dynamic sudden occurs, are carried out here. (1) Unexpected interharmonic test For example, a 10% 2.5th interharmonic occurs at 1.0 s. The whole change processes of TVEs, phase errors and some flags are respectively shown in Fig. 12. As soon as the interharmonic arrives, the state flag is set to 1 until 1.1 s. During this period, the frequency module is locked. After this period, the signal frequency components are re-analyzed after 0.2 s (10 cycles) and the filters coefficients are updated. As shown in the figure, the TVEs and phase errors during 1.1 s to 1.3 s are larger than those after 1.3 s. The response time is about 300 ms, after which the effect of interharmonic interference can be eliminated by updating the filters coefficients. The filter’s characteristics in S channel before 1.3 s and after 1.3 s are shown in Fig. 13. It can be seen that a pair of new notches is added on the FIR filter to eliminate the unexpected interharmonic. Although it is not possible to respond immediately to interference component changes, it has been significantly improved compared to traditional phasor algorithms that do not consider the online variation of signal interference components. These results demonstrate that the proposed algorithm can effectively and adaptively deal with unexpected interferences online.

5.2.4. Noise rejection test The signal noise in distribution networks generally contains the signal inherent noise and the measurement noise (including channel noise and quantization noise). The signal noise is usually assumed to obey a Gaussian distribution; therefore, it is distributed uniformly over the entire frequency band, making it difficult to be removed. Regrettably, the noise that may reduce the signal-to-noise ratio (SNR) to approximately 60 dB [32] is not considered at all in the IEEE stds. However, the noise rejection effects of the PMU algorithms are wellassessed in this paper. The test results based on a pure cosine signal with different levels of noises are shown in Fig. 10. We can directly observe that the performances become worse as the noise level increases in general. Among the algorithms considered, the RWT algorithm is relatively insensitive to the noise level changes;

(2) Unexpected dynamic test To illustrate the algorithm response when an unexpected dynamic occurs, the amplitude modulation is chosen as a representative dynamic behavior. The amplitude modulation suddenly appears at 1.0 s and lasts for 0.5 s. The whole processes of the measured errors of the parameters such as synchrophasor, frequency, ROCOF and state flag, are shown in Fig. 14. The state flag is set to 1 from 1.0 s to 1.54 s (i.e. about 2 cycles after the dynamic behavior is over). The frequency module is locked during 11

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al.

this period. The first time frequency analyzing results are obtained at the time 1.74 s, which is about ten cycles after the state flag is reset to zero. Overall, the TVEs of synchrophasor, the frequency errors and the ROCOF errors are generally small in the whole process, except for the transitional stages, which shows validity and reliability of the proposed algorithm under sudden dynamic occurrence condition.

and ROCOF under conditions containing multiple interferences, dynamic behaviors and unexpected interference/dynamic. This algorithm is now at the stage of development and testing on a specialized System on Chip (SoC) based device. In our future project, the PMUs utilizing this algorithm will be produced and further demonstrated in a real distribution network.

5.2.7. Multiple dynamic behaviour tests Finally, a signal with a series of dynamic behaviors is designed to further validate the adaptive adjustment and switching capability of the TSA algorithm. Several harmonics, interharmonics and 65-dB noise are also taken into consideration. The sequence diagram of the dynamic behaviors and their settings are shown in Fig. 15. The results are compared with the actual values in Fig. 16. The results in Fig. 16 include TVE, amplitude and phase errors of the proposed algorithm. In the first figure, it is indicated that the TVEs are generally lower than 0.1%, whenever the signal is under frequency ramp, modulation or step change condition. The second, fourth and fifth figures show that the amplitude, phase, frequency and ROCOF of synchrophasor can well track the actual values, although the testing signal contains various disturbances. In the third figure, the phase errors are less than 0.05° except for the moments of state transition, which satisfies the measurement requirements at the distribution level. Finally, the last figure illustrates the good consistency between the state flag and the event process in the Fig. 15, which demonstrates that the changes can be effectively recognized and the outputs can be correctly switched by the proposed TCS algorithm.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported in part by National Key Research and Development Program of China (2017YFB0902800) and in part by National Natural Science Foundation of China (51677097, U1766214). Appendix A. Supplementary material Supplementary data to this article can be found online at https:// doi.org/10.1016/j.ijepes.2019.105714. References [1] Lu C, Shi B, Wu X, et al. Advancing China’s smart grid: Phasor measurement units in a wide-area management system. IEEE Power Energy Mag 2015;15(3):60–71. [2] Meier AV, Culler D, McEachern A, et al. Micro-synchrophasors for distribution systems. Proc. IEEE PES Innovative Smart Grid Technol. (ISGT), Washington, DC, USA. Feb 2014. p. 1–5. [3] Arnold DB, Roberts C, Ardakanian O, Stewart EM. Synchrophasor data analytics in distribution grids. Proc. IEEE PES Innovative Smart Grid Technol. (ISGT), Washington, DC, USA. 2017. p. 2472–8152. [4] Meier AV, Stewart E, Mceachern A. Precision micro-synchrophasors for distribution systems: a summary of applications. IEEE Trans Smart Grid 2017;8(6):2926–36. [5] Romano P, Paolone M, Chau T, et al. A High-performance, Low-cost PMU Prototype for Distribution Networks based on FPGA. Proc. IEEE PES Powertech, Manchester, UK. 2017. [6] Li Y, Saha T, Krause O, et al. An inductively active filtering method for powerquality improvement of distribution networks with nonlinear loads. IEEE Trans Power Del 2013;28(4):2465–73. [7] Li L, Ma M, Xu B, et al. Power-quality of distribution networks with high penetrated intermittent distributed generation: a survey. Proc. IEEE Power System Technology (POWERCON), Chengdu, China. 2014. p. 2933–9. [8] Wen H, Zhang J, Yao W, et al. FFT-based amplitude estimation of power distribution systems signal distorted by harmonics and noise. IEEE Trans Indus Inform 2017;14(4):1447–55. [9] Dahal S, Mithulananthan N, Saha T. An approach to control a photovoltaic generator to damp low frequency oscillations in an emerging distribution system. Proc. IEEE PES General Meeting, Detroit, MI, USA. 2011. p. 1–8. [10] Naidu SR, De Andrade GV, Da Costa EG. Voltage sag performance of a distribution system and its improvement. IEEE Trans Ind Appl 2012;48(1):218–24. [11] IEEE Standard for Synchrophasor Measurements for Power Systems, IEEE Standard C37.118.1-2011, 2011. [12] IEEE Standard for Synchrophasor Measurements for Power Systems— Amendment 1: Modification of Selected Performance Requirements, IEEE Standard C37.118.1a2014, Apr. 2014, pp. 1–25. [13] Benmouyal G. An adaptive sampling-interval generator for digital relaying. Trans Power Del 1989;4(3):1602–9. [14] Wang MH, Sun YZ. A practical, precise method for frequency tracking and phasor estimation. IEEE Trans Power Del 2004;19(4):1547–52. [15] Macii D, Petri D, Zorat A. Accuracy analysis and enhancement of DFT-based synchrophasor estimators in off-nominal conditions. IEEE Trans Instrum Meas 2012;61(10):2653–64. [16] De la O Serna JA, Rodriguez J. Instantaneous dynamic phasor estimates with Kalman filter. Proc. IEEE PES Gen. Meeting. 2010. p. 1–6. [17] Dash P, Jena R, Panda G. An extended complex Kalman filter for frequency measurement of distorted signals. IEEE Trans Instrum Meas 2000;49(4):746–53. [18] Platas-Garza MA, De la O Serna JA. Dynamic phasor and frequency estimates through maximally flat differentiators. IEEE Trans Instrum Meas 2010;59(7):1803–11. [19] Belega D, Fontanelli D, Petri D. Dynamic phasor and frequency measurements by an improved Taylor weighted least squares algorithm. IEEE Trans Instrum Meas 2015;64(8):2165–78. [20] Ren J, Kezunovic M. Real-time power system frequency and phasors estimation using recursive wavelet transform. IEEE Trans Power Del 2011;26(3):1392–402. [21] Rahmati A, Adhami R, Dimassi M. Real-time electrical variables estimation based

5.2.8. Profiling CPU time This algorithm is realized based on Visual Studio 2008 software. The computing time results of each part of the algorithm are listed in Table 6. Three voltage channels and three current channels (3U3I) are simulated. The frequency analysis task is handled in a single thread to avoid affecting the phasor estimation. The steady-steady and dynamic algorithms using a lookup table approach only take a short amount of time. Although some additional computational burdens are incurred by transient detection and output selection, they also require only a very short amount of time. As shown in the Table 6, the time consumed by the phasor estimation is 82.76us, which is much less than the inherent latency of the filters in the proposed algorithm. As discussed in the Section 5.1, the reporting latency of this algorithm is about 20 ms as a whole, which meets the requirements of P-Class PMU (ie. 20 ms for the reporting rate at 100/s) and much less than the M-Class PMU limitation (ie. 50 ms for the reporting rate at 100/s). 6. Conclusion Based on the application contexts and the special requirements of the synchrophasor measurements in distribution networks, this paper presented a proposed two-channel switching PMU algorithm to balance the high-precision and fast-response demands. The OOB interferences are specifically considered in the steady-state and dynamic signal modeling. Next, the two adaptive filters are well-designed and online adjusted by adding notches to suppress the OOB interferences; this task is a common challenge for traditional algorithms. Two additional functions: a super-resolution frequency analysis method and a transient recognition strategy with thresholds self-adaptively update, are also designed to support the two parallel channels, respectively for the filters’ adjustment and appropriate estimates’ selection. Finally, the effectiveness of the proposed algorithm is adequately tested not only according to the IEEE stds. but also under some more severe but realistic scenarios. The results have shown the proposed algorithm’s ability on online frequency analysis and transient detection. Besides, the results have also demonstrated that this algorithm can accurately estimate and quickly track not only the fundamental phasor but also the frequency 12

Electrical Power and Energy Systems 117 (2020) 105714

Y. Wang, et al. on recursive wavelet transform. Int J Electr Power Energy Syst 2015;68(6):170–9. [22] Platas-Garza MA, De la O Serna JA. Dynamic harmonic analysis through TaylorFourier transform. IEEE Trans Instrum Meas 2011;60(3):804–13. [23] Platas-Garza MA, De la O Serna JA. Polynomial implementation of the TaylorFourier transform for harmonic analysis. IEEE Trans Instrum Meas 2014;63(12):2846–54. [24] Bertocco M, Frigo G, Narduzzi C, Muscas C, Pegoraro PA. Compressive sensing of a Taylor-Fourier multifrequency model for synchrophasor estimation. IEEE Trans Instrum Meas 2015;64(12):3274–83. [25] Kamwa I, Samantaray SR, Joos G. Wide frequency range adaptive phasor and frequency PMU algorithms. IEEE Trans Smart Grid 2014;5(2):569–79. [26] Chakir M, Kamwa I, Huy HL. Extended C37.118.1 PMU algorithms for joint tracking of fundamental and harmonic phasors in stressed power systems and microgrids. IEEE Trans Power Del 2014;29(3):1465–80. [27] Castello P, Lixia M, Muscas C, et al. Adaptive Taylor-Fourier synchrophasor estimation for fast response to changing conditions. Proc. IEEE Instrumentation and Measurement Technology Conference (I2MTC), Graz, Austria. 2012. [28] Roscoe AJ, Abdulhadi IF, Burt GM. P and M class phasor measurement unit algorithms using adaptive cascaded filters. IEEE Trans Power Del 2013;28(3):1447–59. [29] Castello P, Liu J, Muscas C, et al. A fast and accurate PMU algorithm for P + M class measurement of synchrophasor and frequency. IEEE Trans Instrum Meas 2014;63(12):2837–45. [30] Carta A, Locci N, Muscas C. A PMU for the measurement of synchronized harmonic phasors in three-phase distribution networks. IEEE Trans Instrum Meas

2009;58(10):3723–30. [31] Wei Y, Wang K, Chen Y. A single-phase Phasor Measurement Unit for smart distribution systems. Proc. IEEE Power Electronics and Motion Control, Hefei, AH, CHN. 2016. p. 3559–65. [32] Zhan L, Liu Y, Culliss J, et al. Dynamic single-phase synchronized phase and frequency estimation at the distribution level. IEEE Trans Smart Grid 2015;6(4):2013–22. [33] Liu Y, You S, Yao W, et al. A Distribution Level Wide Area Monitoring System for the Electric Power Grid–FNET/GridEye. IEEE Access 2017;5:2329–38. [34] Kaplan H, Nacaroglu A. A new window function for FIR filter design. Proc. IEEE Signal Processing and Communications Applications, Malatya, TUR. 2015. p. 2458–61. [35] Janik P, Waclawek Z. MUSIC algorithm for estimation of parameters of signals in power system’. Proc. IEEE Environment and Electrical Engineering, Rome, Italy. 2015. p. 2236–40. [36] Jain SK, Jain P, Singh SN. A fast harmonic phasor measurement method for smart grid applications. IEEE Trans Smart Grid 2017;8(1):493–502. [37] Jin T, Liu S, Flesch CC. Mode identification of low-frequency oscillations in power systems based on fourth-order mixed mean cumulant and improved TLS-ESPRIT algorithm’. IET Gener Transm Distrib 2017;11(15):3739–48. [38] IEC Standard for Electromagnetic compatibility (EMC)-Part 4-7: Testing and Measurement Techniques-General guide on harmnics and interharmonices measurements and instrumentation, for power supply systems and equipment connected thereto, IEC Standard 61000-4-7-2009 (Edition 2.1), 2009.

13