
Research Article


A Method of Least Squares Error with Iterative Resampling for Frequency Tracking in Smart Grid 

Hui Li,
ChunJie Cao
and
Liang Yuan



ABSTRACT

Phase angle and amplitude of a phasor can be provided by method of Least Squares Error (LSE). The LSE adopted iteratively is able to track the frequency and amplitude of power system in steady states and in kinds of nonsteady ones. Iterative LSE with resampling is analyzed in conditions of offnominal input, nominal input with harmonics and decaying Direct Current (DC) offset and additive Gaussian white noise. In the circumstances of frequency and phase step changes, performance of resampling LSE is compared with traditional LSE. Resampling LSE has better performance than the traditional one in frequency tracking ability and can provide less mean square error.





Received: June 04, 2015;
Accepted: July 23, 2015;
Published: September 14, 2015


INTRODUCTION
Frequency is one of the most important situation variables and operation parameter in power system. It is a good indicator of integrity of a power system facing separation and islanding. Phasor Measurement Unit (PMU) provides an indication of lost generation, as an example, a frequency drop of 0.1 Hz is typical in the Western Electricity Coordinating Council for 800 MW generation loss (NAERC., 2010). In the past several decades researchers have paid much attention to frequency measurement and analysis in power engineering. Types of frequency estimation methods have been reported, such as zero crossing (Begovic et al., 1993; Nguyen and Srinivasan, 1984), demodulation technique (Begovic et al., 1993), Newton algorithm (Terzija et al., 1994), Kalman filter (Wood et al., 1985; Routray et al., 2002; Siavashi et al., 2009), prony approach (Lobos and Rezmer, 1997), artificial neural network (Vianello et al., 2010), etc. in timedomain. In the transformdomain, the Discrete Fourier Transform (DFT)/fast fourier transform (Yang and Liu, 2000; Zeng et al., 2011) is widely adopted as well as discrete wavelet transform (Chaari et al., 1996; Huang et al., 1999; Lin et al., 2002; Ren and Kezunovic, 2011).
Frequency can be calculated by Least Squares Error (LSE) with a pseudoinverse of a matrix that is determined by a frequency ω_{0} and samples at time t_{m} (m = 1, 2,…, M, samples at t_{1}~t_{M} constituting a data window) of the input signal (Sachdev and Baribeau, 1979). Using LSE iteratively is able to track system frequency samplebysample (Sidhu and Sachdev, 1998). Its estimation accuracy is determined by the number of iterations per sample interval and by computation ability of hardware in realtime application.
In section II, we firstly presented basic idea of method of LSE, then we put forward a guide line on designing coefficients matrix for signal containing harmonics. In section III, we narrated the process of frequency tracking by resampling LSE. Simulation results of above algorithms in different scenarios were compared and some specialties of iterative LSE were given in section IV. Finally we concluded this study. MATERIALS AND METHODS Traditional LSE method: Suppose the nominal voltage or current signal of a power system is: where, H is the order of highest harmonic and A_{i}, ω_{i} = i×ω_{i}, φ_{i} are the amplitude, angular frequency and initial phase angle of each harmonic. According to rule of Nyquist, least sampling frequency f_{s}≥2×Hf_{1 }should be chosen to eliminate aliasing and thanks to low pass filter, components with frequency higher than Hf_{1} are filtered, so we have:
where, V = [s(t_{1}), s(t_{2}),…, s(t_{2H})]^{T}, in which s(t_{1}), s(t_{2}) and s(t_{2H}) are the samples taken at time t_{1}, t_{2 }and t_{2H} (at least 2H samples are taken); X = [X_{1}, X_{2},…, X_{2H1}, X_{2H}]^{T}, in which X_{1} = A_{1} cosφ_{1}, X_{2} = A_{1} sinφ_{1},… and X_{2H1} = A_{H} cosφ_{H}, X_{2H} = A_{H} sinφ_{H} and the matrix:
in which a_{11} = sin(ω_{1}t_{1}), a_{12} = cos(ω_{1}t_{1}), a_{13} = sin(ω_{2}t_{1}), a_{14} = cos(ω_{2}t_{1}),…, a_{1_2H1} = sin(ω_{H}t_{1}), a_{1_2H} = cos(ω_{H}t_{1}),… and a_{2H_1} = sin(ω_{1}t_{2H}), a_{2H_2} = cos(ω_{1}t_{2H}), a_{2H_3} = sin(ω_{2}t_{2H}), a_{2H_4} = cos(ω_{2}t_{2H}),…, a_{2H_2H} = sin(ω_{H}t_{2H}), a_{2H_2H} = cos(ω_{H}t_{2H}). When, a fault or a disturbance occurs, the current signal consists of exponentially decaying DC offsets in electrical power system (ElRefaie and Megahed, 2010; Benmouyal, 1995; Kezunovic et al., 1992; Sidhu et al., 2003). The decaying rates depend on the time constants determined by the inductive reactance to resistance ratio (X/R ratio) of the system (Diniz de Oliveira et al., 2012; Balamourougan and Sidhu, 2006; Lin and Liu, 2002; Belega and Petri, 2011). The larger the X/R ratio, the slower the DC component decays (Kang et al., 2009; Gu and Yu, 2000; Lee et al., 2009). Signal with the nominal component and a decaying DC offset is represented as:
s(t) = A_{1} cos(ω_{1}t+φ_{1})+A_{dc}e‾^{t/τ} (0≤t≤T_{DC}) 
(3) 
where A_{dc} and τ are the amplitude and time constant of the decaying DC offset component. The T_{DC} is the effecting period of decaying DC offset, since the decaying DC offset exists only in several cycles and then disappears. It is a nonperiodic signal and its frequency spectrum encompasses all the frequencies which cannot be removed by antialiasing low pass filter. Matrix A has to be rewritten according to input signal as:
where, a_{m1} = 1, a_{m2} = t_{m}, a_{m3} = t^{2}_{m}, a_{m4} = sin(ω_{1}t_{m}), at time t_{m } (m≥5) for iterative LSE. In order to obtain X, we have:
X = [A^{T} A]‾^{1} A^{T}V = A^{+}V 
(4) 
We get phasor with amplitude Y^{1} and phase angle θ^{1}:
Y(t_{1}) = Y^{1}e^{jθ1} 
(5) 
Where:
θ^{1} = atan[Y(t_{1})_{image}/Y(t_{1})_{real}] 
(6) 
where, atan is the inverse tangent function. The Y(t_{1})_{image} and Y(t_{1})_{real} are the imaginary part and real part of phasor Y(t_{1}) and:
In order to get Y(t_{2}), we have to sample the input signal at time t_{2} and t_{3}, whose time interval is T_{s}. The data window to get Y(t_{2}) is shifted by one T_{s} from that of Y(t_{1}). So the angular frequency is calculated by f_{1} = (θ^{2}θ^{1})/(2πT_{s}). Generally we have frequency at t_{m}:
f_{1}^{m} = (θ^{m+1}θ^{m})/(2πT_{s}) 
(7) 
Traditional LSE algorithm uses fixed sampling frequency of 720 Hz as specified in Sachdev and Baribeau (1979) and Sidhu and Sachdev (1998). However, frequency of the input may change under different conditions, such as offnominal signal in kinds of transient states or under sudden step changes. In these circumstances one data window may contain noninteger number of samples that leads to spectrum leakage and error of estimation. If we change the rate of sampling according to previous estimated frequencies, there may be integer number of sample in one data window and then the error of frequency would be minimized.
 Fig. 1:  Illustration of shifting windows of a cosine waveform with timetag at starting time of t = 0 sec 
Also the nominal signal may be polluted by harmonics, interharmonics, decaying DC offsets, additive Gaussian with noise, etc. Windowbased data acquiring process is shown in Fig. 1.
Suppose frequency of input signal is an offnominal one. By Eq. 7, we can estimate the frequency iteratively by resetting the matrix A with different estimated frequency f'. In reality, sampling rate is set properly to fulfill the requirement of Nyquist’s rule and the criterion of sampling in Terzija et al. (1994). Secondly, an overdetermined linear set is built to give a more precise estimation with more samples. Steps of iterative LSE method are:
Step 1: 
Sample the signal with frequency f_{s} and form a data window 
Step 2: 
Obtain matrix A according to sampling frequency at time t_{m}, m =1, 2,…, M, where M is the length of data window. M is bigger than length of vector X to construct an overdetermined linear set 
Step 3: 
Obtain X according to Eq. 4 
Step 4: 
Obtain phase angle of phasor Y(t_{1}) and Y(t_{2}) by Eq. 6 and then estimated frequency f' by Eq. 7. Number of iteration is 1 
Step 5: 
Resample the input signal and reset the matrix A according to estimated frequency f' at time t_{m}, m = 1, 2,…, M. Repeat steps 1~4 and get f". Number of iteration is 2 
Step 6: 
If f'f"≤δ, iteration process converge, where δ is a preset minimal number; and also if number iteration approaches the maximal number N_{ite_max}, iteration process stops 
Frequency estimation by iterative resamplingLSE method: In order to track the frequency samplebysample, we have to use the process of the iterative LSE per sampling time interval. The initial frequency to get matrix A in the following interval is also the estimated one of previous interval. Number of iterations per sampling interval, limited by computation ability of a digital processor, influences frequency tracking ability and convergence speed of the estimations process. Method of LSE adopts iteration process to get the offnominal frequency more and more precisely. More iteration is accomplished during one time interval, more accurate we get but more time is consumed. Flow chart of resampling LSE is in Fig. 2. RESULTS Suppose the input signal contains a 3rd order harmonic and a decaying DC offset lasting in several cycles at the beginning of simulation as shown in Eq. 8. Parameters used in simulation are listed in Table 1:
X(t) = A_{1} cos(ω_{1}t+φ_{1})+A_{3} cos (ω_{3}t+φ_{3})+A_{dc}e‾^{t/τ}(0≤t) 
(8) 
We can design the matrix A with the help of Taylor expansion as:
 Fig. 2:  Flow chart of resampling LSE 
Table 1:  Basic parameters utilized in simulation 
 ^{+}The cycle is calculated according to the nominal frequency 
where, a_{m1} = 1, a_{m2} = t_{m}, a_{m3} = t^{2}_{m}, a_{m4} = sin(ω_{1}t_{m}), a_{m5} = cos(ω_{1}t_{m}), a_{m6} = sin(ω_{3}t_{m}), a_{m7} = cos(ω_{3}t_{m}) at time t_{m} (m≥7); ω_{3} = 2πf_{3} and X = [X_{1}, X_{2},…, X_{7}]^{T}, where X_{1} = A_{dc}, X_{2} = A_{dc}/τ, X_{3} = A_{dc}/(2τ^{2}), X_{4} = A_{1} cosφ_{1}, X_{5} = A_{1} sinφ_{1}, X_{6} = A_{3} cosφ_{3}, X_{7} = A_{3} sinφ_{3}.
At first, we choose the best sampling rate and data window for LSE algorithms, since the performance of these algorithms are highly connected with them. Former studies show that lower sampling rate and longer data window help to obtain a better result (Sachdev and Baribeau, 1979; Sidhu and Sachdev, 1998). The f_{input} denotes the frequency of input signal.
In Fig. 3, we find that the number of iterations during each sampling interval is set to be 2 or 3, which is good enough to acquire a high accuracy for resampling LSE and comparing with the traditional one, the resampling LSE method tends to converge more robustly. In Fig. 4, we know that tracking accuracy of resampling LSE in higher than that of the traditional one. The dynamic range of resampling LSE is smaller than that of the traditional one. Figure 5 shows the Mean Square Error (MSE) of two LSE methods with different variables such as time constant and amplitude of decaying DC offset component. When, the amplitude of decaying DC offset is 0.2 p.u. and the time constant is 0.1 T, the MSE are 1.2106 and 0.5222 for traditional LSE and resampling LSE, respectively. Figure 6 tell us that resampling LSE does not work well when signal to noise ratio (SNR) is less than 25~30 dB. In the practice of power system, SNR is always higher than 30 dB. So Gaussian white noise could not do too much harm to resampling LSE in practice.

Fig. 3(ab): 
Performance of traditional and resampling LSE methods with nominal input frequency (Initial sampling rate f_{s} = 720 Hz, data window L = 14, f_{input} = 60 Hz, first 30 estimated frequencies are plotted), (a) Number of iteration N_{ite} = 2 in each sampling interval and (b) Number of iteration N_{ite} = 3 in each sampling interval 

Fig. 4(ab): 
Performance of traditional and resampling LSE methods with offnominal input frequency (Initial sampling rate f_{s} = 720 Hz, data window L = 14, Number of iteration N_{ite} = 2, first 30 estimated frequencies are plotted), (a) f_{input} = 55 Hz and (b) f_{input} = 65 Hz 
 Fig. 5:  Performance of traditional and resampling LSE methods with different decaying DC offset (Initial sampling rate f_{s} = 720 Hz, data window L = 14, Number of iteration N_{ite} = 2, f_{input} = 60 Hz, first 30 estimated frequencies are used for MSE calculation) 

Fig. 6(ab): 
Performance of traditional and resampling LSE methods with noise (Initial sampling rate f_{s} = 720 Hz, data window L = 14, f_{input} = 60 Hz, first 30 estimated frequencies are used for MSE calculation and averaged in 1000 times), (a) Number of iteration N_{ite} = 2 and (b) Number of iteration N_{ite} = 3 
DISCUSSION To study the dynamic tracking ability of LSE, step changes of input signals under various conditions are proposed by IEEE Std. The C37.118 (IEEE Standard, 2006, 2011). We compare the result of resampling LSE with the traditional one in the following parts, numerical comparison are listed in the following tables and figures. Frequency step change:
• 
Scenario 1: Nominal input signal. The input signal contains only the nominal sine signal. The first 1 Hz frequency change occurs at 2.6×T, the second 1 Hz frequency change occurs at 6.9×T 
• 
Scenario 2: Nominal input signal with DC offset and 3rd order harmonics, without noise 
• 
Scenario 3: Nominal input signal with DC offset, 3rd order harmonics and additive Gaussian white noise, SNR = 30 dB 
Table 2:  Performance of traditional and resampling LSE methods with number of iterations (window length L = 19) 
 MSE: Mean square error 
For Table 2 and 3 and Fig. 7ac, we find resampling LSE algorithms has 6~11 dB gains comparing with traditional LSE. Because of iterations in sampling interval, there are still jitters in frequency tracking as shown in above figures. The overshoots appear around the time that step changes happen, however the peak value and MSE of resampling LSE is less than that of traditional LSE. We choose L = 19 as the length of data window which is longer than that of (Sidhu and Sachdev, 1998), since longer data window can minimize the influence of Gaussian white noise and also can improve the accuracy of frequency tracking ability but more computation burden and storage requirement.
Table 3:  Performance of traditional and resampling LSE methods with number of iterations (window length L = 19; 1000 time average) 
 MSE: Mean square error 

Fig. 7(af): 
Comparison of resampling LSE and traditional LSE (Number of iteration N_{ite} = 2), (a) Scenario 1, (b) Scenario 2, (c) Scenario 3, (d) Scenario 4, (e) Scenario 5 and (f) Scenario 6 
Phase step change:
• 
Scenario 4: Nominal input signal. The first 5° phase change happens at 2.6 T, the second 5° phase change happens at 6.9 T 
• 
Scenario 5: Nominal input signal with DC offset and 3rd order harmonics, without noise 
• 
Scenario 6: Nominal input signal with DC offset, 3rd order harmonics and Gaussian white noise. SNR = 30 dB 
Table 4:  Performance of traditional and resampling LSE methods with the number of iterations, frequency is off nominal one (window length L = 19) 
 MSE: Mean square error 
Table 5:  Performance of traditional and resampling LSE methods with the number of iterations, frequency is off nominal one (window length L = 19; 1000 time average) 
 MSE: Mean square error 
Form Table 4 and 5 and Fig. 7df, we can draw the flowing conclusions. There are huge overshot near the point of frequency step change and jitter phenomenon when the input frequency is the offnominal one. There are certain specialties of iterative LSE: Design of matrix A should conform to input signal, in another word, designer has to estimate how many harmonics are contained in the signal after low pass filter. Otherwise, the result of LSE would be error, which is difficult to be revised by iteration. Data window of LSE is defined by the number of unknown parameters in every row of matrix A. A longer data window contains more data which are more easily distorted by dynamic change of power system and it also causes longer delay for reporting PMU to phasor data concentrator. Performance of LSE and convergence speed are affected by designing of matrix A. If several decaying DC offset components are contained simultaneously in the signal, it is a huge problem for designing matrix A. Method of LSE is sensitive to signal distortions and signal step change caused by faults of power system, since it is a timedomain algorithm without abilities of denoising and antialiasing.
CONCLUSIONS
We proposed a method of resampling LSE, whose performance was analyzed and compared with the traditional one, when the incoming signal containing harmonics and decaying DC offset. The key point of resampling LSE is that it resamples the input signal according to frequency calculated in previous sampling interval and the iterations processes are accomplished in each interval, which is limited by the computation ability of digital signal processor chips. Iterative LSE can track the nominal frequency and the offnominal one in presence of additive Gaussian white noise with harmonics. The resampling one is less sensitive to DC offset than traditional one. Resampling LSE shows better performance of frequency tracking comparison with the old one with less MSE and better performance in dynamic states of a power system.
ACKNOWLEDGMENT The authors would like to gratefully acknowledge the financial support from the National Nature Science Foundation of China under grant No. 61071128, Nature Science Foundation of Zhejiang Province under grant No. LY15F010003 and Scientific Research Foundation of Hainan University under grant No. kydq1536.

REFERENCES 
1: NAERC., 2010. Realtime application of synchrophasors for improving reliability. North American Electric Reliability Corporation, (NAERC).
2: Begovic, M.M., P.M. Djuric, S. Dunlap and A.G. Phadke, 1993. Frequency tracking in power networks in the presence of harmonics. IEEE Trans. Power Delivery, 8: 480486. CrossRef  Direct Link 
3: Nguyen, C.T. and K. Srinivasan, 1984. A new technique for rapid tracking of frequency deviations based on level crossings. IEEE Trans. Power Apparatus Syst., 103: 22302236. CrossRef 
4: Terzija, V.V., N.B. Djuric and B.D. Kovacevic, 1994. Voltage phasor and local system frequency estimation using Newton type algorithm. IEEE Trans. Power Delivery, 9: 13681374. CrossRef  Direct Link 
5: Wood, H.C., N.G. Johnson and M.S. Sachdev, 1985. Kalman filtering applied to power system measurements relaying. IEEE Trans. Power Apparatus Syst., 104: 35653573. CrossRef  Direct Link 
6: Routray, A., A.K. Pradhan and K.P. Rao, 2002. A novel Kalman filter for frequency estimation of distorted signals in power systems. IEEE Trans. Instrum. Measur., 51: 469479. CrossRef  Direct Link 
7: Siavashi, E.M., S. Afsharnia, M.T. Bina, M.K. Zadeh and M.R. Baradar, 2009. Frequency estimation of distorted signals in power systems using particle extended Kalman filter. Proceedings of the 2nd International Conference on Power Electronics and Intelligent Transportation System, Vol. 1, December 1920, 2009, Shenzhen, pp: 174178.
8: Lobos, T. and J. Rezmer, 1997. Realtime determination of power system frequency. IEEE Trans. Instrum. Measur., 46: 877881. CrossRef 
9: Vianello, R., M.O. Prates, C.A. Duque, A.S. Cequeira, P.M. da Silveira and P.F. Ribeiro, 2010. New phasor estimator in the presence of harmonics, dc offset and interharmonics. Proceedings of the 14th International Conference on Harmonics and Quality of Power, September 2629, 2010, Bergamo, pp: 15.
10: Yang, J.Z. and C.W. Liu, 2000. A new family of measurement technique for tracking voltage phasor, local system frequency, harmonics and DC offset. IEEE Power Eng. Soc. Summer Meeting, 3: 13271332. CrossRef 
11: Zeng, B., Z. Teng, Y. Cai, S. Guo and B. Qing, 2011. Harmonic phasor analysis based on improved FFT algorithm. IEEE Trans. Smart Grid, 2: 5159. CrossRef 
12: Chaari, O., M. Meunier and F. Brouaye, 1996. Wavelet: A new tool for the resonant grounded power distribution systems relaying. IEEE Trans. Power Delivery, 11: 13011308. CrossRef  Direct Link 
13: Huang, S.J., C.T. Hsieh and C.L. Huang, 1999. Application of morlet wavelets to supervise power system disturbances. IEEE Trans. Power Delivery, 14: 235243. CrossRef 
14: Lin, X., H. Zhang, P. Liu and O.P. Malik, 2002. Wavelet based scheme for detection of torsional oscillation. IEEE Trans. Power Syst., 17: 10961101. CrossRef 
15: Ren, J. and M. Kezunovic, 2011. Realtime power system frequency and phasors estimation using recursive wavelet transform. IEEE Trans. Power Delivery, 26: 13921402. CrossRef 
16: Sachdev, M.S. and M.A. Baribeau, 1979. A new algorithm for digital impedance relays. IEEE Trans. Power Apparatus Syst., 98: 22322240. CrossRef 
17: Sidhu, T.S. and M.S. Sachdev, 1998. An iterative technique for fast and accurate measurement of power system frequency. IEEE Trans. Power Delivery, 13: 109115. CrossRef 
18: ElRefaie H.B. and A.I. Megahed, 2010. A novel technique to eliminate the effect of decaying DC component on DFT based phasor estimation. Proceedings of the IEEE Power Energy Society General Meeting, July 2529, 2010, Minneapolis, MN, pp: 18.
19: Benmouyal, G., 1995. Removal of DCoffset in current waveforms using digital mimic filtering. IEEE Trans. Power Delivery, 10: 621630. CrossRef  Direct Link 
20: Kezunovic, M., P. Spasojevic and B. Perunicic, 1992. New digital signal processing algorithms for frequency deviation measurement. IEEE Trans. Power Delivery, 7: 15631573. CrossRef 
21: Sidhu, T.S., X. Zhang, F. Albasri and M.S. Sachdev, 2003. DiscreteFouriertransformbased technique for removal of decaying DC offset from phasor estimates. IEE Proc. Generat. Trans. Distribut., 150: 745752. CrossRef 
22: Diniz de Oliveira, A., L.R.M. Silva, C.H. Martins, R.R. Aleixo, C.A. Duque and A.S. Cerqueira, 2012. An improved DFT based method for phasor estimation in fault scenarios. Proceedings of the IEEE Power and Energy Society General Meeting, July 2226, 2012, San Diego, CA., pp: 15.
23: Balamourougan, V. and T.S. Sidhu, 2006. A new filtering technique to eliminate decaying dc and harmonics for power system phasor estimation. Proceedings of the IEEE Power India Conference, April 1012, 2006, New Delhi, pp: 16.
24: Lin, Y.H. and C.W. Liu, 2002. A new DFTbased phasor computation algorithm for transmission line digital protection. Proceedings of the Asia Pacific IEEE/PES Transmission and Distribution Conference and Exhibition, October 610, 2002, Yokohama, pp: 17331737.
25: Belega, D. and D. Petri, 2011. Accuracy of a DFT phasor estimator at offnominal frequency in either steady state of transient conditions. Proceedings of the IEEE International Conference on Smart Measurements for Future Grids, November 1416, 2011, Bologna, Italy, pp: 4550.
26: Kang, S.H., D. G. Lee, S.R. Nam, P.A. Crossley and Y.C. Kang, 2009. Fourier transformbased modified phasor estimation method immune to the effect of the DC offsets. IEEE Trans. Power Delivery, 24: 11041111. CrossRef 
27: Gu, J.C. and S.L. Yu, 2000. Removal of DC offset in current and voltage signals using a novel Fourier filter algorithm. IEEE Trans. Power Delivery, 15: 7379. CrossRef  Direct Link 
28: Lee, D.G., Y.J. Oh, S.H. Kang and B.M. Han, 2009. Distance relaying algorithm using a DFTbased modified phasor estimation method. Proceedings of the IEEE Bucharest PowerTech Conference, June 28July 2, 2009, Bucharest, pp: 16.
29: IEEE Standard, 2006. C37.1182005IEEE standard for synchrophasors for power systems. http://ieeexplore.ieee.org/servlet/opac?punumber=10719.
30: IEEE Standard, 2011. C37.118.12011IEEE standard for synchrophasor measurements for power systems. http://ieeexplore.ieee.org/servlet/opac?punumber=6111217.



