Subscribe Now Subscribe Today
Research Article

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

Hui Li, Chun-Jie Cao and Liang Yuan
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

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 non-steady ones. Iterative LSE with resampling is analyzed in conditions of off-nominal 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.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Hui Li, Chun-Jie Cao and Liang Yuan, 2015. A Method of Least Squares Error with Iterative Resampling for Frequency Tracking in Smart Grid. Journal of Applied Sciences, 15: 1093-1102.

DOI: 10.3923/jas.2015.1093.1102

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


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 time-domain. In the transform-domain, 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 pseudo-inverse of a matrix that is determined by a frequency ω0 and samples at time tm (m = 1, 2,…, M, samples at t1~tM constituting a data window) of the input signal (Sachdev and Baribeau, 1979). Using LSE iteratively is able to track system frequency sample-by-sample (Sidhu and Sachdev, 1998). Its estimation accuracy is determined by the number of iterations per sample interval and by computation ability of hardware in real-time 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.


Traditional LSE method: Suppose the nominal voltage or current signal of a power system is:


where, H is the order of highest harmonic and Ai, ωi = i×ωi, φi are the amplitude, angular frequency and initial phase angle of each harmonic. According to rule of Nyquist, least sampling frequency fs≥2×Hf1 should be chosen to eliminate aliasing and thanks to low pass filter, components with frequency higher than Hf1 are filtered, so we have:

V = A×X

where, V = [s(t1), s(t2),…, s(t2H)]T, in which s(t1), s(t2) and s(t2H) are the samples taken at time t1, t2 and t2H (at least 2H samples are taken); X = [X1, X2,…, X2H-1, X2H]T, in which X1 = A1 cosφ1, X2 = A1 sinφ1,… and X2H-1 = AH cosφH, X2H = AH sinφH and the matrix:

in which a11 = sin(ω1t1), a12 = cos(ω1t1), a13 = sin(ω2t1), a14 = cos(ω2t1),…, a1_2H-1 = sin(ωHt1), a1_2H = cos(ωHt1),… and a2H_1 = sin(ω1t2H), a2H_2 = cos(ω1t2H), a2H_3 = sin(ω2t2H), a2H_4 = cos(ω2t2H),…, a2H_2H = sin(ωHt2H), a2H_2H = cos(ωHt2H). 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) = A1 cos(ω1t+φ1)+Adce‾t/τ (0≤t≤TDC)

where Adc and τ are the amplitude and time constant of the decaying DC offset component. The TDC is the effecting period of decaying DC offset, since the decaying DC offset exists only in several cycles and then disappears. It is a non-periodic signal and its frequency spectrum encompasses all the frequencies which cannot be removed by anti-aliasing low pass filter. Matrix A has to be rewritten according to input signal as:

where, am1 = 1, am2 = tm, am3 = t2m, am4 = sin(ω1tm), at time tm (m≥5) for iterative LSE. In order to obtain X, we have:

X = [AT A]‾1 ATV = A+V

We get phasor with amplitude Y1 and phase angle θ1:

Y(t1) = Y1ejθ1


θ1 = atan[Y(t1)image/Y(t1)real]

where, atan is the inverse tangent function. The Y(t1)image and Y(t1)real are the imaginary part and real part of phasor Y(t1) and:

In order to get Y(t2), we have to sample the input signal at time t2 and t3, whose time interval is Ts. The data window to get Y(t2) is shifted by one Ts from that of Y(t1). So the angular frequency is calculated by f1 = (θ21)/(2πTs). Generally we have frequency at tm:

f1m = (θm+1m)/(2πTs)

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 off-nominal signal in kinds of transient states or under sudden step changes. In these circumstances one data window may contain non-integer 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 time-tag at starting time of t = 0 sec

Also the nominal signal may be polluted by harmonics, inter-harmonics, decaying DC offsets, additive Gaussian with noise, etc. Window-based data acquiring process is shown in Fig. 1.

Suppose frequency of input signal is an off-nominal 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 over-determined 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 fs and form a data window
Step 2: Obtain matrix A according to sampling frequency at time tm, m =1, 2,…, M, where M is the length of data window. M is bigger than length of vector X to construct an over-determined linear set
Step 3: Obtain X according to Eq. 4
Step 4: Obtain phase angle of phasor Y(t1) and Y(t2) 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 tm, 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 Nite_max, iteration process stops

Frequency estimation by iterative resampling-LSE method: In order to track the frequency sample-by-sample, 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 off-nominal 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.


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) = A1 cos(ω1t+φ1)+A3 cos (ω3t+φ3)+Adce‾t/τ(0≤t)

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, am1 = 1, am2 = tm, am3 = t2m, am4 = sin(ω1tm), am5 = cos(ω1tm), am6 = sin(ω3tm), am7 = cos(ω3tm) at time tm (m≥7); ω3 = 2πf3 and X = [X1, X2,…, X7]T, where X1 = Adc, X2 = -Adc/τ, X3 = Adc/(2τ2), X4 = A1 cosφ1, X5 = A1 sinφ1, X6 = A3 cosφ3, X7 = A3 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 finput 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(a-b):
Performance of traditional and resampling LSE methods with nominal input frequency (Initial sampling rate fs = 720 Hz, data window L = 14, finput = 60 Hz, first 30 estimated frequencies are plotted), (a) Number of iteration Nite = 2 in each sampling interval and (b) Number of iteration Nite = 3 in each sampling interval

Fig. 4(a-b):
Performance of traditional and resampling LSE methods with off-nominal input frequency (Initial sampling rate fs = 720 Hz, data window L = 14, Number of iteration Nite = 2, first 30 estimated frequencies are plotted), (a) finput = 55 Hz and (b) finput = 65 Hz

Fig. 5:
Performance of traditional and resampling LSE methods with different decaying DC offset (Initial sampling rate fs = 720 Hz, data window L = 14, Number of iteration Nite = 2, finput = 60 Hz, first 30 estimated frequencies are used for MSE calculation)

Fig. 6(a-b):
Performance of traditional and resampling LSE methods with noise (Initial sampling rate fs = 720 Hz, data window L = 14, finput = 60 Hz, first 30 estimated frequencies are used for MSE calculation and averaged in 1000 times), (a) Number of iteration Nite = 2 and (b) Number of iteration Nite = 3


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. 7a-c, 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(a-f):
Comparison of resampling LSE and traditional LSE (Number of iteration Nite = 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. 7d-f, 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 off-nominal 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 time-domain algorithm without abilities of de-noising and anti-aliasing.


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 off-nominal 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.


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.

1:  NAERC., 2010. Real-time 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: 480-486.
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: 2230-2236.
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: 1368-1374.
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: 3565-3573.
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: 469-479.
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 19-20, 2009, Shenzhen, pp: 174-178.

8:  Lobos, T. and J. Rezmer, 1997. Real-time determination of power system frequency. IEEE Trans. Instrum. Measur., 46: 877-881.
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 26-29, 2010, Bergamo, pp: 1-5.

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: 1327-1332.
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: 51-59.
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: 1301-1308.
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: 235-243.
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: 1096-1101.
CrossRef  |  

15:  Ren, J. and M. Kezunovic, 2011. Real-time power system frequency and phasors estimation using recursive wavelet transform. IEEE Trans. Power Delivery, 26: 1392-1402.
CrossRef  |  

16:  Sachdev, M.S. and M.A. Baribeau, 1979. A new algorithm for digital impedance relays. IEEE Trans. Power Apparatus Syst., 98: 2232-2240.
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: 109-115.
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 25-29, 2010, Minneapolis, MN, pp: 1-8.

19:  Benmouyal, G., 1995. Removal of DC-offset in current waveforms using digital mimic filtering. IEEE Trans. Power Delivery, 10: 621-630.
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: 1563-1573.
CrossRef  |  

21:  Sidhu, T.S., X. Zhang, F. Albasri and M.S. Sachdev, 2003. Discrete-Fourier-transform-based technique for removal of decaying DC offset from phasor estimates. IEE Proc. Generat. Trans. Distribut., 150: 745-752.
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 22-26, 2012, San Diego, CA., pp: 1-5.

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 10-12, 2006, New Delhi, pp: 1-6.

24:  Lin, Y.H. and C.W. Liu, 2002. A new DFT-based phasor computation algorithm for transmission line digital protection. Proceedings of the Asia Pacific IEEE/PES Transmission and Distribution Conference and Exhibition, October 6-10, 2002, Yokohama, pp: 1733-1737.

25:  Belega, D. and D. Petri, 2011. Accuracy of a DFT phasor estimator at off-nominal frequency in either steady state of transient conditions. Proceedings of the IEEE International Conference on Smart Measurements for Future Grids, November 14-16, 2011, Bologna, Italy, pp: 45-50.

26:  Kang, S.H., D. G. Lee, S.R. Nam, P.A. Crossley and Y.C. Kang, 2009. Fourier transform-based modified phasor estimation method immune to the effect of the DC offsets. IEEE Trans. Power Delivery, 24: 1104-1111.
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: 73-79.
CrossRef  |  Direct Link  |  

28:  Lee, D.G., Y.J. Oh, S.H. Kang and B.M. Han, 2009. Distance relaying algorithm using a DFT-based modified phasor estimation method. Proceedings of the IEEE Bucharest PowerTech Conference, June 28-July 2, 2009, Bucharest, pp: 1-6.

29:  IEEE Standard, 2006. C37.118-2005-IEEE standard for synchrophasors for power systems.

30:  IEEE Standard, 2011. C37.118.1-2011-IEEE standard for synchrophasor measurements for power systems.

©  2021 Science Alert. All Rights Reserved