An Optimization of Least-Square Harmonic Phasor Estimators in Presence of Multi-Interference and Harmonic Frequency Variance

Abstract

The wide application of power electronic devices brings an increasing amount of undesired harmonic and interharmonic tones, and accurate harmonic phasor estimation under a complex signal input is an important task for smart grid applications. In this paper, an optimization of least-square dynamic harmonic phasor estimators, considering multi-interference and harmonic frequency variance, is proposed. A comprehensive error index (CEI) composed of the fundamental-leakage-led harmonic amplitude estimation error, harmonic mutual interference, out-of-band interference, and harmonic frequency deviation is employed. The largest CEI part of least-square algorithms using three different signal decomposition models is analyzed for the first time, and variables to reduce this error component are then introduced using singular value decomposition. With the CEI and defined variables, a minimum-error estimation of harmonic phasors under various interference and harmonic frequency change is discussed. Numerical tests are performed, and the test results show that after the proposed optimization is applied to least-square algorithms, the harmonic phasor estimation errors are considerably reduced, especially for low-order harmonics. We also show the possibility of choosing desired optimal phasor filter design by balancing the measurement accuracy and data latency. For example, when the window length is set to three nominal cycles, the proposed optimization can yield both good accuracy and fast measurement speed for estimating harmonic phasors under multi-interference and harmonic frequency variance.

1. Introduction

The development of low-carbon power systems and smart grids greatly promotes the application of power electronic equipment [1,2]. Although the power electronic device improves the efficiency and flexibility of power utilization, it also brings an increasing amount of undesired harmonic and interharmonic tones [1]. For example, considerable harmonic and interharmonic components have been observed in the current of DC arc furnaces [3]. These undesired components distort the signal waveform, deteriorate power quality, cause oscillations, and even lead to power system failures [4,5,6]. Measuring harmonic components is the first and critical step to locating harmonic sources, and then yields the possibility of suppressing harmonic components by using active filters [7]. Additionally, the accurate harmonic phasor estimation is of great importance for some advanced applications, e.g., identification of a distribution topology [8], high impedance fault location [9], and arcing fault detection [10], etc. However, the increasingly complex signal environment hinders the accurate estimation of harmonic phasors.

The signal complexity and measurement algorithm are two core elements ruling the performance of a harmonic phasor estimator. The standards for synchrophasor estimation, e.g., IEC/IEEE 60255-118-1 [11], summarize a variety of complex signals into several typical characteristics, such as amplitude deviation, harmonic interference, out-of-band interference (OBI), amplitude modulation, phase modulation, as well as frequency deviation or linear change. Similar conditions apply for harmonic measurements, and usually two major targets, i.e., the accuracy and response speed, are most concerned in phasor estimation. The accuracy relies on the suppression of various interference and how well the harmonic frequency change is tracked. The speed concern expects as short time window length as possible to get the estimation, which is, however, likely to be detrimental to the estimation accuracy [12].

Classical synchrophasor estimators based on DFT employ window function or interpolation to help improve the performance, and the approaches work well under dynamic conditions. However, they are less effective for the harmonic phasor estimation [13,14], because they can not provide enough compensation in the wider frequency deviation range. To improve the estimation accuracy for harmonic phasors, researchers have developed several frequency-domain algorithms. For example, the IEC 61000-4-7 recommended a grouping harmonic method based on DFT and Parseval’s theorem to calculate the harmonic amplitude [15], but it is unable to calculate the harmonic phase. Reddy et al. employed compressive sensing to estimate dynamic harmonic phasor, yielding a high frequency resolution at the expanse of a heavier computation burden [16].

Time-domain methods have also been developed to improve the dynamic performance of harmonic phasor estimations. The Taylor signal model was adopted using least-square or Kalman to lower the passband ripple, e.g., Taylor–Fourier transform (TFT) [17], discrete-time TFT filter banks designed with O-Splines [18], and Taylor–Kalman filter [19]. Choosing an appropriate window function also reduces the passband ripple, and increases the attenuation at mutual harmonic interference frequencies [20]. The passband ripple and, hence, the estimation error of these dynamic algorithms, however, increase with the harmonic order [21]. To address this issue, research works have explored harmonic filters with passbands that increase with the harmonic order, and, in this case, the passband ripple stops growing with the harmonic order. For example, sinc interpolation function [21], and frequency-domain sampling theorem [22] have been proposed to replace the Taylor series function to reach the above targets. However, these least-square algorithms suffer from large transition band gain [23], and, hence, are sensitive to the possible OBI tones, e.g., a communication signal at 178 Hz around the third harmonic [15], some interharmonics at 330 Hz, 380 Hz, and 530 Hz, respectively, around the sixth, seventh, and tenth harmonics in an electric power system [24]. Increasing the window length can mitigate the negative effect from OBI, but it may not meet the fast response requirement when measuring strong dynamic signal in active distribution networks [25]. In [26], we proposed to enable the TFT algorithm with improved suppression of OBI tones without increasing window length. However, the proposal is less practical because it does not consider the harmonic frequency deviation and dynamic change.

For the suppression of interference, especially the OBI, harmonic phasor filters with notches at the interference frequencies are required. A key step in such an approach is to detect the component number and frequency from signal samples, which in practice can be realized by employing, e.g., the IpDFT [27], Prony [4], matrix pencil [28], or ESPRIT [1] methods. The least-square is finally applied to estimate the harmonic phasors, and usually satisfying performance, even under dynamic and OBI conditions, can be achieved. The drawback of the notch filter algorithms is that they produce a heavy real-time computation burden, and the estimation error can increase considerably when the signal component number and frequency are not well identified.

The above research works tend to focus on harmonic signals with a variant frequency or OBI tone. For example, the least-square algorithms have good dynamic performance but weak suppression of interharmonics while the SVDHPE is the opposite [26]. Therefore in this paper, we aim to keep the good dynamic performance of least-square algorithms, and enable them with good suppression of OBI, fundamental leakage, and harmonic mutual interference under a short time window, e.g., three nominal cycles. First, in Section 2, the comprehensive error index (CEI) is employed to indicate the harmonic amplitude estimation error caused by the multi-interference and harmonic frequency deviation. The largest source contributing to the CEI for harmonic phasor filters using three different least-square algorithms is analyzed. Then in Section 3, new variables are introduced to the least-square equations by singular value decomposition (SVD), and a minimum error optimization is delivered based on the CEI and the defined variables. Next, in Section 4, the proposed minimum error estimation of harmonic phasors under various interference and harmonic frequency change is tested using numerical tests. Finally, the applicability of the proposed optimization is analyzed in Section 5, and the conclusion is drawn in Section 6.

To summarize, the novelty and contributions of this paper mainly include: (1) The interharmonic is found as the largest error source for the harmonic phasor estimation with least-square algorithms, and hence some variables to suppress interharmonic interference and reduce CEI are introduced; (2) An optimization to minimize the CEI by changing the introduced variables is proposed. The optimal design, compared to other approaches, takes both the different interference suppressions and the dynamic accuracy as optimization targets; (3) We conduct a detailed theoretical analysis and design a set of typical scenarios to fully evaluate the proposed algorithm in terms of accuracy, response performance, and computational burden.

2. The Comprehensive Error Index, CEI

In this section, we employ the comprehensive error-index to indicate the performance of harmonic phasor filters under the multi-influence of interference and harmonic frequency deviation. It begins with the signal model without losing generality as

$$\begin{array}{c}\hfill x\left(t\right)=\sum _{h=1}^H(a_h\cos (2\pi h{f}_{1}t+{\varphi }_h)+a_{h,i}\cos (2\pi f_{h,i}t+{\varphi }_{h,i}))=\sum _{h=1}^H{x}_h\left(t\right)+\sum _{h=1}^H{x}_{h,i}\left(t\right),\end{array}$$

(1)

where the signal contains the fundamental $x_1\left(t\right)$, harmonic $x_h\left(t\right)$, and OBI component $x_{h,i}\left(t\right)$. The OBI tone means the signal component resulting in spectral aliasing when the phasor data are transferred to the master station, and in (1), there are many OBI tones around the fundamental and harmonics. The fundamental’s frequency, amplitude, and phase are $f_1$, $a_1$, and ${\varphi }_1$. The harmonic $x_h\left(t\right)$ has an integer multiple frequencies of the fundamental, i.e., $h{f}_1$, with an amplitude $a_h$ and phase ${\varphi }_h$. The OBI component $x_{h,i}\left(t\right)$ has a frequency pertaining to $[(h-1)f_0,h{f}_0-0.5f_{\mathrm{r}}]\cup [h{f}_0+0.5f_{\mathrm{r}},(h+1)f_0]$, where nominal frequency $f_0=50$Hz, and reporting rate [11] $f_{\mathrm{r}}=50$frames/s in this paper. The amplitude of the OBI component is $a_{h,i}$, while its phase is ${\varphi }_{h,i}$. According to [11], the synchrophasor and harmonic phasor are defined as $p_h\left(t\right)=a_h{\mathrm{e}}^{\mathrm{j}{\varphi }_h\left(t\right)}$. Then, the time-domain form of the fundamental and harmonic signals can be rewritten as

$$\begin{array}{c}\hfill x_h\left(t\right)=\frac{1}{2}{p}_h\left(t\right){\mathrm{e}}^{\mathrm{j}2\pi h{f}_{1}t}+\frac{1}{2}{p}_h^{*}\left(t\right){\mathrm{e}}^{-\mathrm{j}2\pi h{f}_{1}t},\end{array}$$

(2)

where ${}^{*}$ denotes the conjugate operation. Considering that $p_h\left(t\right)$ changes slowly with time, then from (2), the Fourier transform of $x_h\left(t\right)(h\in [1,H])$ is obtained as

$$\begin{array}{c}\hfill X_h\left(\omega \right)=\pi (p_h\left(\omega \right)\otimes \delta (\omega -h{\omega }_1)+p_h^{*}\left(\omega \right)\otimes \delta (\omega +h{\omega }_1)),\end{array}$$

(3)

where the angular frequency $\omega =2\pi f$ and ${\omega }_1=2\pi f_1$; $\delta (·)$ represents the unit impulse function; ⊗ denotes the convolution operation. Suppose that the harmonic phasor filter has an ideal magnitude-frequency response as follows

$$\begin{array}{c}\hfill L_h^{\mathrm{i}}\left(\omega \right)=\left\{\begin{array}{c}1,\omega \in 2\pi [h(f_0-v),h(f_0+v)],\hfill \\ 0,\mathrm{others},\hfill \end{array}\right.\end{array}$$

(4)

where v represents the fundamental frequency deviation, i.e., the filter passband, and hence the harmonic frequency deviation is $hv$. Then, the positive frequency part of $X_h\left(\omega \right)$ in (3) can be obtained by multiplying the $h-$order harmonic phasor filter and the Fourier transform of signal $x\left(t\right)$, i.e., $X\left(\omega \right)$, as

$$\begin{array}{c}\hfill X\left(\omega \right)L_h^{\mathrm{i}}\left(\omega \right)=\pi p_h\left(\omega \right)\otimes \delta (\omega -h{\omega }_1).\end{array}$$

(5)

The corresponding time-domain form is

$$\begin{array}{c}\hfill x\left(t\right)\otimes l_h^{\mathrm{i}}\left(t\right)=\frac{1}{2}{p}_h\left(t\right){\mathrm{e}}^{\mathrm{j}2\pi h{f}_{1}t},\end{array}$$

(6)

Let the signal $x\left(t\right)$ be sampled at frequency $f_{\mathrm{s}}$, and the center of every time window be tagged as the initial time $t=0$, then the harmonic amplitude and phase can be calculated by

$$\begin{array}{cc}\hfill a_h\left(t_r\right)=& 2|x\left(n{T}_{\mathrm{s}}\right)\otimes l_h^{\mathrm{i}}\left(n{T}_{\mathrm{s}}\right)|,\hfill \\ \hfill {\varphi }_h\left(t_r\right)=& \angle (x\left(n{T}_{\mathrm{s}}\right)\otimes l_h^{\mathrm{i}}\left(n{T}_{\mathrm{s}}\right)),\hfill \end{array}$$

(7)

where $t_r$ denotes the Coordinated Universal Time (UTC) attached to the center time of the r-th time window; $n\in [-N_h,\cdots N_h]$ represents the data number in a time window; $T_{\mathrm{s}}=1/f_{\mathrm{s}}$ is the sampling interval.

Considering that the implemented harmonic filter often has non-zero gain out of the passband, with interference, the filtering produces a harmonic phasor estimation as

$$\begin{array}{c}\hfill x\left(t\right)\otimes l_h^{\mathrm{r}}\left(t\right)=x\left(t\right)\otimes l_h^{\mathrm{i}}\left(t\right)+x\left(t\right)\otimes (l_h^{\mathrm{r}}\left(t\right)-l_h^{\mathrm{i}}\left(t\right))=\frac{1}{2}{p}_h\left(t\right){\mathrm{e}}^{\mathrm{j}2\pi h{f}_{1}t}+o\left(t\right),\end{array}$$

(8)

where $l_h^{\mathrm{r}}\left(t\right)$ denotes the implemented filter. Compared with (6), the undesired term $o\left(t\right)$ will cause harmonic phasor estimation error, and its Fourier transform is

$$\begin{array}{cc}\hfill O\left(\omega \right)=& \pi (p_h\left(\omega \right)\otimes \delta (\omega -h{\omega }_1)(L_h^{\mathrm{r}}\left(\omega \right)-L_h^{\mathrm{i}}\left(\omega \right))+p_h^{*}\left(\omega \right)\otimes \delta (\omega +h{\omega }_1)L_h^{\mathrm{r}}\left(\omega \right)\hfill \\ \hfill & +\sum _{i=1,i\ne h}^{2H}{p}_i\left(\omega \right)\otimes \delta (\omega -{\omega }_i)L_h^{\mathrm{r}}\left(\omega \right)+\sum _{i=1,i\ne h}^{2H}{p}_i^{*}\left(\omega \right)\otimes \delta (\omega +{\omega }_i)L_h^{\mathrm{r}}\left(\omega \right)),\hfill \end{array}$$

(9)

where $L_h^{\mathrm{r}}\left(\omega \right)$ is the magnitude-frequency response of $l_h^{\mathrm{r}}\left(t\right)$; ${\omega }_i$ represents the signal component in the signal $x\left(t\right)$ after omitting the $h-$order harmonic. For the $h-$order harmonic phasor filter, some indices used in the following text are defined as

$$\begin{array}{cc}\hfill {\delta }_{h,h}^{\mathrm{p}}& =\max |L_h^{\mathrm{r}}\left(h{\omega }_1\right)-L_h^{\mathrm{i}}\left(h{\omega }_1\right)|,\hfill \\ {\delta }_{h,h}^{\mathrm{n}}\hfill & =\max |L_h^{\mathrm{r}}(-h{\omega }_1)|,\hfill \\ {\delta }_{h,i}^{\mathrm{p}}\hfill & =\max |L_h^{\mathrm{r}}\left({\omega }_i\right)|,\hfill \\ {\delta }_{h,i}^{\mathrm{n}}\hfill & =\max |L_h^{\mathrm{r}}(-{\omega }_i)|,\hfill \end{array}$$

(10)

where $\max \left|Y\right|$ denotes the maximum absolute value of Y. The four indices defined in (10) respectively represent the maximum gain error between $L_h^{\mathrm{r}}\left(\omega \right)$ and $L_h^{\mathrm{i}}\left(\omega \right)$ in the passband, the conjugate range of the passband, the positive frequency range of the interference, and the negative frequency range of the conjugate components for the $h-$order harmonic phasor filter. Specifically, for the h-order harmonic, $h{\omega }_1\in 2\pi h[f_0-v,f_0+v]$. For the fundamental and $i-$order ($i\in [1,H],i\ne h$) harmonic, ${\omega }_i\in 2\pi i[f_0-v,f_0+v]$. For OBI tones, i.e., $i\in [H+1,2H]$, ${\omega }_i\in 2\pi [(i-H-1)f_0,(i-H)f_0-0.5f_{\mathrm{r}}]\cup 2\pi [(i-H)f_0+0.5f_{\mathrm{r}},(i-H+1)f_0]$. Note that the OBI frequency range containing the passband of the $h-$order harmonic filter needs to be eliminated, i.e., when $i=h-1+H$, ${\omega }_i\in 2\pi [(i-H-1)f_0,(i-H)f_0-0.5f_{\mathrm{r}}]$, and when $i=h+1+H$, ${\omega }_i\in 2\pi [(i-H)f_0+0.5f_{\mathrm{r}},(i-H+1)f_0]$). Combining (9) and (10), the maximum amplitude of $o\left(t\right)$, i.e., the inverse Fourier transform of $O\left(\omega \right)$, is

$$\begin{array}{c}\hfill {\left|o\left(t\right)\right|}_{\max }=\frac{1}{2}\left(\right|p_h\left(t\right)|{\delta }_{h,h}^{\mathrm{p}}+|p_h^{*}\left(t\right)|{\delta }_{h,h}^{\mathrm{n}}+\sum _{i=1,i\ne h}^{2H}\left(\right|p_i\left(t\right)|{\delta }_{h,i}^{\mathrm{p}}+|p_i^{*}\left(t\right)|{\delta }_{h,i}^{\mathrm{n}}\left)\right).\end{array}$$

(11)

As a result, a comprehensive error-index, standing the $h-$order harmonic amplitude estimation error attached to the implemented filter $l_h^{\mathrm{r}}\left(t\right)$, is defined as [29]

$$\begin{array}{cc}\hfill {\mathrm{CEI}}_h& =\frac{|{\widehat{a}}_h-a_h|}{a_h}=\frac{\left|2\right|x\left(t\right)\otimes l_h^{\mathrm{r}}\left(t\right)|-a_h|}{a_h}\le \frac{{2\left|o\left(t\right)\right|}_{\max }}{a_h}\hfill \\ \hfill & \le {\delta }_{h,h}^{\mathrm{p}}+{\delta }_{h,h}^{\mathrm{n}}+\sum _{i=1,i\ne h}^{2H}(k_i{\delta }_{h,i}^{\mathrm{p}}+k_i{\delta }_{h,i}^{\mathrm{n}}),\hfill \end{array}$$

(12)

where $\widehat{Y}$ denotes the estimation of Y; $k_i=a_i/a_h$ is a parameter related to the amplitudes of harmonic and interharmonic tones. Equation (12) shows that ${\mathrm{CEI}}_h$ results from the joint action of non-ideal filtering, and in the presence of interference and harmonic parameter change. ${\mathrm{CEI}}_h$ can be divided into three parts: $({\delta }_{h,h}^{\mathrm{p}}+{\delta }_{h,h}^{\mathrm{n}})$ is related to the non-ideal harmonic parameters such as the harmonic frequency deviation; ${\sum }_{i=1,i\ne h}^H(k_i{\delta }_{h,i}^{\mathrm{p}}+k_i{\delta }_{h,i}^{\mathrm{n}})$ is with the fundamental leakage and harmonic mutual interference, and ${\sum }_{i=H+1}^{2H}(k_i{\delta }_{h,i}^{\mathrm{p}}+k_i{\delta }_{h,i}^{\mathrm{n}})$ with the OBI tones.

With the CEI defined, it can then be used to analyze the errors of the existing dynamic harmonic phasor estimators. Here we take as an example the recently developed least-square algorithm based on the frequency domain sampling theorem (LS-FST). The LS-FST analysis begins with the decomposition model of the discretized harmonic phasor [22], i.e.,

$$\begin{array}{c}\hfill p_h\left(n{T}_{\mathrm{s}}\right)=\sum _{k=-K}^K{p}_{h,k}{\mathrm{e}}^{\mathrm{j}2\pi k\Delta f_{h}n{T}_{\mathrm{s}}},\end{array}$$

(13)

where the harmonic model order $Q=2K+1$; $p_{h,k}$ denotes the k-th sample of the discrete time Fourier transform of $p_h\left(n{T}_{\mathrm{s}}\right)$, and $\Delta f_h$ is the sample interval. With this decomposition, the $N=2N_h+1$ signal samples in a time window for signal $s\left(n{T}_{\mathrm{s}}\right)={\sum }_{h=1}^H{x}_h\left(n{T}_{\mathrm{s}}\right)=\frac{1}{2}{\sum }_{h=1}^H(p_h\left(n{T}_{\mathrm{s}}\right){\mathrm{e}}^{\mathrm{j}2\pi h{f}_{1}n{T}_{\mathrm{s}}}+p_h^{*}\left(n{T}_{\mathrm{s}}\right){\mathrm{e}}^{-\mathrm{j}2\pi h{f}_{1}n{T}_{\mathrm{s}}})$ can yield the following overdetermined system of equations, i.e.,

$$\begin{array}{c}\hfill \mathit{s}={\left[\begin{array}{c}{\mathit{E}}_1\\ ⋮\\ {\mathit{E}}_H\\ {\mathit{E}}_1^{*}\\ ⋮\\ {\mathit{E}}_H^{*}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cccccc}{\mathit{F}}_1& & & & & \mathbf{0}\\ & \ddots & & & & \\ & & {\mathit{F}}_H& & & \\ & & & {\mathit{F}}_1^{*}& & \\ & & & & \ddots & \\ \mathbf{0}& & & & & {\mathit{F}}_H^{*}\end{array}\right]\left[\begin{array}{c}{\mathit{p}}_1\\ ⋮\\ {\mathit{p}}_H\\ {\mathit{p}}_1^{*}\\ ⋮\\ {\mathit{p}}_H^{*}\end{array}\right]=\mathit{EFp}=\mathit{Gp},\end{array}$$

(14)

where the column vector $\mathit{s}\in {\mathbb{R}}^{N\times 1}$ contains N samples of signal $s\left(t\right)$, i.e., $s(-N_h{T}_{\mathrm{s}}),\cdots ,$$s\left(N_h{T}_{\mathrm{s}}\right)$; the diagonal matrix ${\mathit{E}}_h\in {\mathbb{C}}^{N\times N}$ has the diagonal elements as ${\mathrm{e}}^{\mathrm{j}2\pi h{f}_{1}n{T}_{\mathrm{s}}}$ ($-N_h\le n\le N_h,1\le h\le H$); the matrix ${\mathit{F}}_h\in {\mathbb{C}}^{N\times Q}$ and its $(k+K+1)$-th column elements are ${\mathrm{e}}^{\mathrm{j}2\pi k\Delta f_{h}n{T}_{\mathrm{s}}}(-K\le k\le K)$; the column vector ${\mathit{p}}_h\in {\mathbb{C}}^{Q\times 1}$ consists of the spectrum sampling parameters to be solved, i.e., $p_{h,-K}/2,...,p_{h,K}/2$; ${\mathit{Y}}^{\mathrm{T}}$ denotes the transpose of the matrix $\mathit{Y}$. The least-square algorithm is employed to estimate the unknown vector $\mathit{p}$, i.e.,

$$\begin{array}{c}\hfill \widehat{\mathit{p}}={\left({\left(\mathit{EF}\right)}^{\mathrm{H}}\left(\mathit{EF}\right)\right)}^{-1}{\left(\mathit{EF}\right)}^{\mathrm{H}}\mathit{s}={\mathit{G}}^{+}\mathit{s},\end{array}$$

(15)

where ${\mathit{Y}}^{\mathrm{H}},{\mathit{Y}}^{-1},{\mathit{Y}}^{+}$ represent the conjugate transpose, inverse, and generalized inverse operations for the matrix $\mathit{Y}$; ${\mathit{G}}^{+}\in {\mathbb{C}}^{2HQ\times N}$. Combining (13) and (15), the h-order harmonic phasor estimation can be obtained by

$$\begin{array}{c}\hfill {\widehat{p}}_h\left(0\right)=\sum _{k=-K}^K{\widehat{p}}_{h,k}=\sum _{i=(h-1)Q+1}^{hQ}\widehat{\mathit{p}}(i,1),\end{array}$$

(16)

where $\widehat{\mathit{p}}(i,1)$ represents the i-th element in the column vector $\widehat{\mathit{p}}$. Based on (15) and (16), the filter coefficients for h-order harmonic phasor are given as

$$\begin{array}{c}\hfill {\mathit{d}}_h=\sum _{i=(h-1)Q+1}^{hQ}{\mathit{G}}^{+}(i,:),\end{array}$$

(17)

where ${\mathit{G}}^{+}(i,:)\in {\mathbb{C}}^{1\times N}$ denotes the i-th row elements of matrix ${\mathit{G}}^{+}$. For harmonic filter obtained by (17), three components are separately indicated following the three terms of CEI${}_h$ defined in (12).

$$\begin{array}{cc}\hfill {\mathrm{EP}}_{h,1}& =100\%\times ({\delta }_{h,h}^{\mathrm{p}}+{\delta }_{h,h}^{\mathrm{n}})/{\mathrm{CEI}}_h,\hfill \\ \hfill {\mathrm{EP}}_{h,2}& =100\%\times \left(\sum _{i=1,i\ne h}^H(k_i{\delta }_{h,i}^{\mathrm{p}}+k_i{\delta }_{h,i}^{\mathrm{n}})\right)/{\mathrm{CEI}}_h,\hfill \\ \hfill {\mathrm{EP}}_{h,3}& =100\%\times \left(\sum _{i=H+1}^{2H}(k_i{\delta }_{h,i}^{\mathrm{p}}+k_i{\delta }_{h,i}^{\mathrm{n}})\right)/{\mathrm{CEI}}_h.\hfill \end{array}$$

(18)

Comparing three components with a numerical example helps to understand the relative contributions to the CEI of LS-FST.

Table 1. Each part proportion of the harmonic amplitude estimation error obtained by different least-square algorithms.

Order	LS-FST (%)			LS-Taylor (%)			LS-Sinc (%)
h	EP${}_{\mathit{h},1}$	EP${}_{\mathit{h},2}$	EP${}_{\mathit{h},3}$	EP${}_{\mathit{h},1}$	EP${}_{\mathit{h},2}$	EP${}_{\mathit{h},3}$	EP${}_{\mathit{h},1}$	EP${}_{\mathit{h},2}$	EP${}_{\mathit{h},3}$
2	0.004	1.751	98.245	0.011	4.292	95.697	0.005	2.062	97.933
3	0.009	1.836	98.155	0.026	4.398	95.576	0.010	2.128	97.862
4	0.018	2.137	97.844	0.055	5.053	94.892	0.018	2.474	97.507
5	0.033	2.608	97.359	0.104	6.098	93.799	0.031	3.007	96.961
6	0.056	3.227	96.717	0.182	7.455	92.363	0.049	3.702	96.250
7	0.090	4.014	95.897	0.299	9.152	90.549	0.072	4.574	95.354
8	0.138	4.985	94.877	0.468	11.194	88.338	0.104	5.632	94.264
9	0.207	6.156	93.637	0.702	13.584	85.714	0.143	6.897	92.959
10	0.304	7.479	92.217	1.024	16.240	82.736	0.193	8.308	91.499
11	0.490	8.868	90.642	1.572	18.964	79.464	0.303	9.770	89.972
12	0.916	10.063	89.021	2.749	21.300	75.951	0.566	11.025	88.409
13	2.864	7.042	90.093	7.377	14.362	78.262	2.053	7.785	90.162

shows an example, and the parameters included in this example are: the nominal frequency $f_0=50$Hz; the sampling frequency $f_{\mathrm{s}}=10$kHz; the time window $T_{\mathrm{w}}=3/f_0$; the signal model order $Q=3$; the maximum harmonic order $H=13$ [22]; the frequency interval $\Delta f_h=0.5h$Hz [22]; the fundamental frequency $f_1=f_0$; the amplitudes of the fundamental, harmonic, and OBI are $a_1=1$pu, $a_h=0.1$pu, and $a_{h,i}=0.01$pu [11], respectively. It can be seen from Table 1 that the EP${}_{h,3}$, coming from the OBI components, takes over 90% weight for all harmonic orders less than 13, and is the largest error contribution for LS-FST harmonic phasor estimators.

The above analysis fits other least-square algorithms. Table 1 also shows the proportion of each CEI part obtained by the other two least-square algorithms, respectively, based on Taylor expansion (LS-Taylor) [17], and sinc interpolation function (LS-Sinc) [21]. The analysis is under the same test example as above for the LS-FST. A similar result, i.e., EP${}_{h,3}$ takes the largest error contribution, is found.

3. Proposed Optimal Design of Dynamic Harmonic Phasor Filter

Since the major error component of least-square algorithms is from OBI tone, reducing the response gain among the OBI frequency range should be an effective way to reduce the CEI${}_h$ and improve the measurement accuracy of least-square algorithms. In this paper, the core idea is to use singular value decomposition to lower the OBI band gain, the third part of the CEI${}_h$. Our previous work has shown that the SVD can efficiently increase the attenuation to OBI for fundamental and harmonic phasor estimation based on LS-Taylor [26,30], and here we implement the approach in harmonic phasor estimation based on LS-FST, LS-Sinc, and LS-Taylor as part of the new improvement target. In the following exploration of SVD potential in harmonic phasor estimation, the latest LS-FST is still used as an example. Only the key formulas are given here, and the detailed derivation process can be found in [26,30].

The analysis begins with decomposing ${\mathit{F}}_h$ by SVD, i.e.,

$$\begin{array}{c}\hfill {\mathit{F}}_h={\mathit{U}}_h{\mathsf{\Sigma }}_h{\mathit{V}}_h^{\mathrm{H}},\end{array}$$

(19)

Using ${\mathit{F}}_h$ in (19), the elements from row $(h-1)Q+1$ to row $hQ$ of the matrix ${\mathit{G}}^{+}$ in (15) can be written as

$$\begin{array}{c}\hfill {\mathit{G}}_h^{+}={\mathit{V}}_h{\mathsf{\Sigma }}_h^{-1}(\sum _{j=1}^H{\mathit{M}}_{h,j}{\mathit{U}}_j^{\mathrm{H}}{\mathit{E}}_j^{*}+\sum _{j=H+1}^{2H}{\mathit{M}}_{h,j}{\mathit{U}}_{(j-H)}^{\mathrm{T}}{\mathit{E}}_{(j-H)}),\end{array}$$

(20)

where ${\mathit{G}}_h^{+}\in {\mathbb{C}}^{Q\times N}(1\le h\le H)$; ${\mathit{M}}_{h,j}\in {\mathbb{C}}^{Q\times Q}(1\le h,j\le 2H)$ is formed by the intersection elements of $(h-1)Q+1$ to $hQ$ rows and columns in the matrix $\mathit{M}$, $\mathit{M}={\left({\mathit{U}}^{\mathrm{H}}{\mathit{E}}^{\mathrm{H}}\mathit{EU}\right)}^{-1}$, and $\mathit{M}\in {\mathbb{C}}^{2HQ\times 2HQ}$; $\mathit{U}=\mathrm{diag}({\mathit{U}}_1,\cdots ,{\mathit{U}}_H,{\mathit{U}}_1^{*},\cdots ,{\mathit{U}}_H^{*})$ and $\mathrm{diag}(·)$ represents a matrix formed by putting matrices in the parentheses on the diagonal; $\mathit{E}$ is defined in (14). Then, the filter coefficients for $h-$order harmonic phasor in (17) can be rewritten as ${\mathit{d}}_h={\sum }_{i=1}^Q{\mathit{G}}_h^{+}(i,:)$. Referring to [26,30], the weight pair ${\mathit{V}}_h{\mathsf{\Sigma }}_h^{-1}$ is chosen to change the filtering performance by introducing variables $y_{h,k_{\mathrm{e}}}$, and then the $h-$order ($h\in [2,H]$) dynamic harmonic phasor filter is designed as

$$\begin{array}{c}\hfill {\overline{\mathit{d}}}_h=\sum _{k_{\mathrm{e}}=1}^Q\frac{{\mathit{w}}_h(1,k_{\mathrm{e}})}{y_{h,k_{\mathrm{e}}}·{\mathsf{\Sigma }}_h(k_{\mathrm{e}},k_{\mathrm{e}})}{\mathit{L}}_h(k_{\mathrm{e}},:),\end{array}$$

(21)

where ${\mathit{w}}_h(1,k_{\mathrm{e}})$ is the $k_{\mathrm{e}}$-th element of the row vector ${\mathit{w}}_h={\sum }_{i=1}^Q{\mathit{V}}_h(i,:)$ (${\mathit{w}}_h\in {\mathbb{C}}^{1\times Q}$), and ${\mathit{V}}_h(i,:)\in {\mathbb{C}}^{1\times Q}$ denotes the i-th row elements of matrix ${\mathit{V}}_h$; ${\mathsf{\Sigma }}_h(k_{\mathrm{e}},k_{\mathrm{e}})$ is the $k_{\mathrm{e}}$-th element in the diagonal of ${\mathsf{\Sigma }}_h$; ${\mathit{L}}_h(k_{\mathrm{e}},:)\in {\mathbb{C}}^{1\times N}$ denotes the $k_{\mathrm{e}}$-th row elements of matrix ${\mathit{L}}_h={\sum }_{j=1}^H{\mathit{M}}_{h,j}{\mathit{U}}_j^{\mathrm{H}}{\mathit{E}}_j^{*}+{\sum }_{j=H+1}^{2H}{\mathit{M}}_{h,j}{\mathit{U}}_{(j-H)}^{\mathrm{T}}{\mathit{E}}_{(j-H)}({\mathit{L}}_h\in {\mathbb{C}}^{Q\times N})$.

Both Refs. [26,30] show the appropriate variable values can reduce the transition band gain, and hence the estimation error from OBI, with the increase of passband ripple, and the error caused by frequency deviation. Whether the CEI${}_h$ can be reduced depends on their joint effects, as well as the harmonic mutual interference. Therefore, we formulate the following minimum optimization with the defined CEI${}_h$ and a total of Q variables introduced with SVD to explore the LS-FST performance improvement under multi-interference and harmonic frequency change.

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}{y}_{h,k_{\mathrm{e}}}\end{array}}{\mathbf{Min}}& {\mathrm{CEI}}_h\hfill \\ \hfill \mathbf{s}.\mathbf{t}.& \forall y_{h,k_{\mathrm{e}}}>0,\hfill \\ \hfill & c\ge Q,\hfill \end{array}$$

(22)

where CEI${}_h$ is calculated by (12) and (10) for the $h-$order harmonic phasor filter designed by (21); $c=T_{\mathrm{w}}{f}_0$ denotes the number of nominal cycles in a time window. The first constraint is used to keep computation stability. The second one is obtained by the following two points: (1) Equation (14) should have a least-square or a unique solution, which requires the number of equations greater than or equal to the number of the unknown sampling parameters, i.e., $cM\ge 2HQ$ where $M=f_{\mathrm{s}}/f_0$ is the sample number in a nominal cycle; (2) The Nyquist sampling theorem requires $0.5M{f}_0\ge H{f}_0$, namely $m=M/2H\ge 1$. Merging two conditions yields $c\ge Q/m$. Since this relationship should always hold, with the largest value on the right at $m=1$, $c\ge Q$. The optimization aims to minimize the CEI${}_h$ with appropriate variable values. For LS-Taylor and LS-Sinc algorithms, a similar optimization problem can be obtained by replacing the signal decomposition model in (13) with the corresponding basis function.

Here a two-step numerical search method is employed to find optimal variable values and the minimum CEI${}_h$. The first step is to observe the CEI${}_h$ change as a function of each $y_{h,k_{\mathrm{e}}}$ (fixing other parameters). With this step, only the variables that help to reduce CEI${}_h$ value are set adjustable, and others are fixed at one. Then, the adjustable variables are enumerated to reach an optimal combination for minimizing CEI${}_h$.

So far, the optimization for the three least-square algorithms is achieved under the given time window c and signal model order Q ($Q=2K+1$ for LS-FST and LS-Sinc; $Q=K+1$ for LS-Taylor). The following is an example to visually display the optimization results. For a simple case, the time window length is set to three nominal cycles, and the signal model order is set to 3 for the least-square algorithms, i.e., $c=Q=3$. Then, three variables, $y_{h,1},y_{h,2},$ and $y_{h,3}$, are introduced. Figure 1 shows their relationships with the CEI${}_h$ while the other two variables are fixed to unchanged one. It can be seen that: (1) The individual change of $y_{h,1}$ toward either direction from one will increase CEI${}_h$, (2) CEI${}_h$ is independent of $y_{h,2}$, and (3) $y_{h,3}$ can reduce CEI${}_h$ when it is greater than one. The optimal $y_{h,3}$ value for each harmonic is chosen at the valley of the green curve. For the same least-square algorithm, the optimal values for different harmonics are slightly different, 1.56–1.66, 1.46–1.59, and 1.56–1.64 for the LS-FST, LS-Taylor, and LS-Sinc, respectively.

The CEI${}_h$ change with and without optimization for the three least-square algorithms is plotted in Figure 2. The results agree with the analysis: (1) EP${}_{h,2}$ and EP${}_{h,3}$ caused by fundamental leakage, mutual harmonic interference, and OBI components are well suppressed after optimization; (2) the EP${}_{h,1}$ related to the filter passband ripple increases, but the overall CEI${}_h$, especially for $h\in [2,10]$, is reduced considerably. This observation is further verified by the frequency responses in Figure 3, i.e., the optimized harmonic filters reduce transition band and stop band gain, but enlarge passband ripple.

It is interesting to know whether the proposed optimization in (22) always has the effect of CEI${}_h$ reduction and harmonic phasor estimation accuracy improvement under different least-square order Q and window length c. Here we show numerical examples in the range {3, 4, 5, 6, 7} for both parameters. $Q=$ 1, and 2 are not considered because they yield only one variable $y_{h,1}$ for optimization, and changing $y_{h,1}$ can only increase the CEI${}_h$ value. Since $c\ge Q$, c can not go lower than 3. The upper limit, 7, is considering the length limit set in IEC/IEEE standard [11]. Figure 4 displays the results of three least-square algorithms with and without the proposed optimization in three columns.

As we can see from Figure 4, the red-line values are always less than the black lines at the same configuration of least-square order and window length, which checks the validation of the proposed optimization. The optimized CEI${}_h$ varies with the harmonic order h, and in most cases, increases with h. The reason is that the passband ripple effect becomes higher at high-order harmonics. Another interesting conclusion observed from Figure 4 is that when keeping Q constant, a larger window length yields a higher optimized CEI${}_h$ value. Considering $c\ge Q$, it suggests using $c=Q$ as the optimal when the window length varies in different applications. The last row of subplots in Figure 4 summarizes the CEI${}_h$ at the optimal $c=Q$. There is no doubt that good CEI reductions are achieved at long window lengths, e.g., $c=5,6,7$, and these optimizations can be used in measurement class schemes. For the proposal, however, we observe better improvement for odd c than for even c. Therefore, a special interest result is for $c=3$: it can achieve a much lower CEI${}_h$ value than $c=4$ and uses the smallest response latency for phasor update. In this case, i.e., $c=3$, the optimization can be used in protection-scheme applications.

To conclude, the diagram of the proposed optimal algorithm for dynamic harmonic phasor filters is shown in Figure 5. As the black boxes show, Equations (1)–(12) construct the CEI${}_h$ for any non-ideal harmonic filter under multi-interference and parameter deviation; The green boxes show that Equations (13)–(18) analyze the largest CEI source (i.e., the OBI tones denoted by EP${}_{h,3}$) for the harmonic filters obtained by the three least-square (LS) algorithms; In Equations (19)–(21) displayed in the blue boxes, variables $y_{h,k_{\mathrm{e}}}$ are introduced to reduce CEI${}_h$ with the help of the SVD, where the relationships between CEI${}_h$ and the introduced variables for different least-square algorithms are shown in Figure 1; Then, Equation (22) yields the maximum CEI${}_h$ reduction by changing the introduced variables, and Figure 2 and Figure 3 show an example of the maximum CEI${}_h$ reduction and filtering response improvement; Finally, we evaluate the optimization performance under different configurations of two parameters, i.e., signal model order Q and window length c as in Figure 4, and recommend $c=Q=3$ for realizing both good accuracy and fast measurement speed.

4. Results

In this section, typical scenarios referring to [11] are carried out to test the estimation accuracy and response speed of the optimized harmonic phasor measurement algorithms. To keep the text concise, we take the case with window length $c=3$ as an example. Conventional LS-FST [22], LS-Taylor [17], and LS-Sinc [21] are used as input algorithms. Their performances with and without the proposed optimization, as well as two state-of-the-art works, the matrix pencil for estimating harmonics and interharmonics (HI-MP) [28], and the ESPRIT for measuring harmonic phasors [1], are compared. The test settings are identical for all eight algorithms: the nominal frequency $f_0=50$Hz; the sampling frequency $f_{\mathrm{s}}=10$kHz; the reporting rate $f_{\mathrm{r}}=50$frames/s; the signal model order $Q=3$. The frequency interval is $0.5h$Hz for LS-FST [22], and the basic frequency band is 0.575 $h$ Hz for LS-Sinc [21] where h is the harmonic order. The total vector error (TVE) and response time defined in [11] are used as indicators. The mean TVE of the signal parameter variation (magnitude, frequency variation, etc.) over 100 runs is plotted in the following graphs.

4.1. Harmonic Amplitude Test

The first four steady tests employ the signal in (1) to test the eight algorithms under different harmonic amplitude, noise interference (white noise is added to signal in (1)), interharmonic amplitude, and harmonic frequency deviation.

As shown in (1), the signal contains the fundamental tone, harmonics with an order from 2 to H, and each component is accompanied by an OBI tone. The fundamental has frequency $f_1=f_0$, amplitude $a_1=1$pu, and phase ${\varphi }_1$ randomly chosen within $[-\pi ,\pi ]$. The harmonic frequency is set to an integer multiple of the fundamental, the maximum harmonic order $H=$ 13 [22], the harmonic amplitudes increase from 8% to 12% of the fundamental in a step of 0.5% in each implementation, i.e., $a_h\in [0.08,0.12]$pu, and the harmonic phases ${\varphi }_h\in [-\pi ,\pi ]$. The OBI tones have frequencies $f_{h,i}=h{f}_0-0.5f_{\mathrm{r}}$ (In the OBI tone frequency range, this is the worst case), amplitudes $a_{h,i}=0.01$ pu, and phases ${\varphi }_{h,i}\in [-\pi ,\pi ](h\in [1,H])$. The fundamental amplitude, maximum harmonic order, OBI frequencies, and all component phases have the same settings as here and it will not be repeated in the following tests.

The mean TVEs of the three least-square algorithms with and without the proposed optimization, as well as HI-MP and ESPRIT methods, are shown in Figure 6. The test results show that the proposed method has a better performance than the HI-MP and ESPRIT algorithms. The main reason is that the number and frequencies of signal components are not accurately identified for the comparing algorithms, causing significant harmonic phasor estimation errors.

Figure 6 also shows that the three least-square algorithms have very similar performances under multi-interference before optimization. Note that this similarity of the three least-square algorithms remains the same in the following tests, and the reason is that they own a comparable ability to suppress OBI components. In Figure 6, after optimization, the total errors caused by the fundamental leakage, 11 harmonic mutual interference (except the estimated harmonic from the 12 harmonics when $H=$ 13), and 13 OBI components are below 5.9%, 3.2%, and 4.3% for LS-FST, LS-Taylor, and LS-Sinc, respectively. The mean TVE reductions within a harmonic order of 10 are over 0.03, 0.043, and 0.038 for the three least-square algorithms. However, the TVEs of optimized LS-FST and LS-Sinc increase with the harmonic order. This is because (1) Their deviations of the filter gain from the ideal value of one at the set nominal frequency become larger with the increase of the harmonic order, which does not exist in the optimized LS-Taylor, and (2) The interharmonic amplitude related parameter $k_i\in [0.1/0.12,0.1/0.08]$ in this test is different from the value used in the optimization $k_i=0.1$.

4.2. Noise Interference Test

This test employs the signal in (1) added with white noise to evaluate the noise immunity of the eight algorithms. The signal settings include: the frequencies of the fundamental and harmonics are $f_1=f_0$ and $h{f}_1$, respectively; the amplitudes of the harmonic and OBI are $a_h=0.1$pu and $a_{h,i}=0.01$pu, respectively. The noise intensity is characterized by the signal-to-noise ratio (SNR) defined as

$$\begin{array}{c}\hfill \mathrm{SNR}=10log\frac{a_1^2}{2{\beta }^2},\end{array}$$

(23)

where $a_1=1$pu represents the fundamental amplitude in (1), and $\beta$ denotes the standard deviation of the white noise.

For each test, the SNR increases from 50 dB to 80 dB with a step of 5 dB, and the average TVEs across different SNRs and 100 running times are shown in Figure 7. After optimization, the ceiling values for average TVEs are reduced from 7.72%, 7.72%, and 7.73% to 4.4%, 3.2%, and 3.7% for 2–10 harmonic phasors obtained by LS-FST, LS-Taylor, and LS-Sinc algorithms, respectively. Moreover, the optimization achieves higher accuracy than the HI-MP and ESPRIT algorithms, and shows good robustness to noise interference.

4.3. OBI Amplitude Test

This test is to estimate the performance of the eight algorithms under different interharmonic amplitudes. The test signal in (1) is set as follows: the fundamental frequency $f_1=f_0$; the harmonic frequencies are $h{f}_1$; the OBI amplitudes $a_{h,i}\in [0.001,0.02]$pu; and all the harmonic amplitudes $a_h=0.1$pu.

In each test, $a_{h,i}$, the amplitudes of all OBI increase from 1% to 20% of the harmonic in a step of 1%. The optimized LS-FST, LS-Taylor, and LS-Sinc respectively achieve a TVE reduction over 0.03, 0.045, and 0.038 for 2–10 harmonics. The results in Figure 8 show that the optimization still produces good accuracy improvement, and achieves a much lower estimation error than the HI-MP and ESPRIT algorithms, while the ceiling OBI amplitude $a_{h,i}$ = 0.02 used in the test is much larger than that $a_{h,i}$ = 0.01 in the optimization analysis. This fact proves the proposed optimization is less sensitive to the OBI amplitude increase.

It is observed that the estimation error curves of the three least-square algorithms with the proposed optimization in Figure 6, Figure 7 and Figure 8 appear to be the same curve trend as the optimized subplots in Figure 2. This is because for high-order harmonic, the error EP${}_{h,1}$ caused by passband ripple takes a larger proportion, and EP${}_{h,3}$ related to the transition band becomes smaller. The proposed optimization realizes CEI reduction by decreasing EP${}_{h,3}$, and hence the CEI reduction becomes smaller with the harmonic order increase. Namely, the proposal yields better performance improvement for low harmonics than for high harmonics.

4.4. Harmonic Frequency Deviation Test

Figure 2 implies that the passband ripples of the optimized harmonic filters are bigger than the original conditions. This will lead to larger errors caused by harmonic frequency deviation in the harmonic phasor estimation. Therefore, this part tests the proposed harmonic estimators and the compared algorithms when there are harmonic frequency deviations and multiple instances of interference. The signal in (1) is also used in this test, and the settings are given as the fundamental frequency $f_1$ changes within [49.5, 50.5] Hz [22] with a step of 0.1 Hz in each test; accordingly, the harmonic frequencies $h{f}_1$ increase from 49.5h Hz to 50.5h Hz in a step of 0.1h Hz; the harmonic and OBI amplitudes are $a_h=0.1$pu and $a_{h,i}=0.01$pu, respectively.

As seen in Figure 9, the proposed optimization still performs a goodly decrease in the mean TVEs, and produces much lower error than the HI-MP and ESPRIT algorithms. The values are always below 4.2%, 4.5%, and 3.9% for the optimized LS-FST, LS-Taylor, and LS-Sinc, respectively. Moreover, we compare the maximum TVEs and the theoretical values in Figure 2. Most differences are less than 0.02 for all three least-square algorithms, which proves the test results are in good agreement with the theoretical analysis.

4.5. Harmonic Amplitude Modulation Test

The following three tests are carried out to check the algorithms’ ability to deal with dynamic conditions, i.e., harmonic modulation in amplitude and phase, and harmonic frequency ramp. The signal for the harmonic amplitude modulation test is

$$\begin{array}{c}\hfill x\left(t\right)=\sum _{h=1}^H{a}_h(1+0.1\cos \left(2\pi f_{\mathrm{m}}t\right))\cos (2\pi h{f}_{1}t+{\varphi }_h)+\sum _{h=1}^H{a}_{h,i}\cos (2\pi f_{h,i}t+{\varphi }_{h,i}),\end{array}$$

(24)

where the fundamental and harmonic amplitudes, i.e., $a_1=1$pu and $a_h=0.1$pu are modulated with level $0.1$, and frequency $f_{\mathrm{m}}\in [0.1,2]$ Hz [11] in a step of 0.1 Hz for each test; the OBI amplitudes $a_{h,i}=0.01$pu; the frequencies of the fundamental and harmonics are $f_1=f_0$ and $h{f}_1$, respectively.

Figure 10 shows the optimized LS-FST, LS-Taylor, and LS-Sinc ensure a TVE reduction over 0.028, 0.041, and 0.035 for 2–10 harmonics, respectively, and achieve more accurate estimation than HI-MP and ESPRIT methods. This observation indicates the proposed optimization still has good performance under dynamic signal condition.

4.6. Harmonic Phase Modulation Test

Equation (25) gives the signal used in the harmonic phase modulation test. The fundamental phase ${\varphi }_1$ and harmonic phases ${\varphi }_h$ are all added with a modulation component having level $0.1h$ and frequency $f_{\mathrm{m}}\in [0.1,2]$Hz [11]. Other settings are as follows: the fundamental has frequency $f_1=f_0$ and amplitude $a_1=1$pu; the harmonics own frequencies $h{f}_1$ and amplitudes $a_h=0.1$pu; the OBI tones have amplitudes $a_{h,i}=0.01$pu.

$$\begin{array}{c}\hfill x\left(t\right)=\sum _{h=1}^H{a}_h\cos (2\pi h{f}_{1}t+0.1h\cos (2\pi f_{\mathrm{m}}t-\pi )+{\varphi }_h)+\sum _{h=1}^H{a}_{h,i}\cos (2\pi f_{h,i}t+{\varphi }_{h,i}).\end{array}$$

(25)

In this test, the mean TVEs resulting from the phase modulation and multi-interference are always below 5.8%, 3.2%, and 4.1% for the three least-square algorithms with the proposed optimization. Figure 11 also indicates that the optimized TVEs are respectively reduced by more than 0.033, 0.043, and 0.04 compared with the original LS-FST, LS-Taylor, and LS-Sinc, and are much lower than the TVEs of HI-MP and ESPRIT methods for 2–10 harmonics.

4.7. Harmonic Frequency Ramp Test

Harmonic frequency ramp may result in a bigger estimation error than frequency deviation when covering the same frequency range. Therefore, the test evaluates the eight algorithms under harmonic frequency ramp and various interference. The test signal is

$$\begin{array}{cc}\hfill x\left(t\right)=& \sum _{h=1}^H{a}_h\cos (2\pi h{f}_{1}t+\pi h{R}_{\mathrm{f}}{t}^2)+\sum _{h=1}^H{a}_{h,i}\cos \left(2\pi f_{h,i}t\right),\hfill \end{array}$$

(26)

where $a_1=1$pu, $a_h=0.1$pu, and $a_{h,i}=0.01$pu; the fundamental frequency changes linearly from 49.5 Hz to 50.5 Hz in one second, i.e., $R_{\mathrm{f}}=1$Hz/s, and correspondingly the harmonic frequencies increase from 49.5h Hz to 50.5h Hz at the ramp $h{R}_{\mathrm{f}}=h$Hz/s.

Figure 12 shows that the optimized LS-FST, LS-Taylor, and LS-Sinc maintain the mean TVEs at a low level. The TVE reduction is over 0.043, 0.046, and 0.047 for harmonics with an order from 2 to 10. Moreover, the error values and volatility of the optimized least-square method are smaller than those of the HI-MP and ESPRIT algorithms.

As seen in the seven tests above, the optimized least-square algorithms achieve good accuracy improvement compared with the original form, and more accurate estimation than the HI-MP and ESPRIT methods in most conditions, which proves the proposal is valid both in steady and dynamic conditions. However, the proposal fails to estimate the thirteenth harmonic phasor for LS-FST. There are two reasons for this failure: (1) The error proportion EP${}_{h,3}$ caused by OBI tones decreases with harmonic order, which means the CEI reduction becomes small for high-order harmonic filters; (2) The interharmonic amplitude related parameter $k_i$ in the tests is different from the value used in the optimization, e.g., $k_i\in [0.083,0.125]$ in the first test, while $k_i=0.1$ in the optimization. In this case, the optimization should be implemented with the appropriate maximum harmonic order H and amplitude parameter $k_i$.

4.8. The Sampling Frequency Test

In this test, the sampling frequency is ramped from 5 kHz to 10 kHz in steps of 1 kHz. The signal in (1) added with a 60 dB white noise is employed, and the signal settings include the following: the fundamental frequency has a deviation of 0.2 Hz, i.e., $f_1=50.2$Hz, and hence the harmonic frequency is $50.2h$Hz; the amplitudes of the harmonic and OBI are $a_h=0.1$pu and $a_{h,i}=0.01$pu, respectively.

The mean TVEs of the second and thirteenth harmonic phasors are shown in Figure 13. The results for 3–12 harmonics have similar trends. It is observed that the three least-square algorithms with the proposed optimization maintain better robustness over different sampling frequencies compared to the HI-MP and ESPRIT methods.

4.9. Harmonic Amplitude and Phase Step Tests

The step tests are implemented to obtain the response times of the optimized algorithms. The test signal includes the fundamental and one harmonic, and they perform steps together, i.e.,

$$\begin{array}{c}\hfill x\left(t\right)=(1+k_{\mathrm{a}}b\left(t\right))\cos (2\pi f_{0}t+k_{\mathrm{p}}b\left(t\right))+0.1(1+k_{\mathrm{a}}b\left(t\right))\cos (2\pi h{f}_{0}t+k_{\mathrm{p}}b\left(t\right)),\end{array}$$

(27)

where the step level configurations $k_{\mathrm{a}}=0.1,k_{\mathrm{p}}=0$ and $k_{\mathrm{a}}=0,k_{\mathrm{p}}=-\pi /18$ are taken for the $h-$order harmonic amplitude and phase step tests, respectively; the harmonic order h increases from 2 to 13 [22] for both step tests; the step function $b\left(t\right)=0,t<0$ and $b\left(t\right)=1,t>0$. To calculate the response times, the errors caused by the non-ideal gain at the nominal harmonic frequency (shown in Figure 3) for the optimized LS-FST and LS-Sinc are compensated.

As shown in

Table 2. The response times of the LS-FST, LS-Taylor, and LS-Sinc with and without the proposed optimization under harmonic amplitude step test.

Order	LS-FST		LS-Taylor		LS-Sinc
h	Original	Optimized	Original	Optimized	Original	Optimized
2	58.0	47.8	58.0	48.6	58.0	47.9
3	58.3	41.9	58.2	43.9	58.3	42.5
4	58.0	36.1	57.8	40.6	58.0	37.6
5	58.3	32.7	58.1	38.8	58.2	34.4
6	57.7	27.3	57.7	35.5	57.7	29.5
7	42.4	23.3	42.3	35.2	42.4	26.8
8	26.0	16.6	26.0	32.7	26.0	21.8
9	25.2	14.1	25.2	32.9	25.2	19.1
10	24.4	24.0	24.3	30.6	24.3	15.0
11	24.2	27.1	24.0	31.2	24.1	25.9
12	22.0	26.4	22.0	30.2	22.0	25.5
13	31.7	29.2	31.8	33.3	33.3	29.1

and

Table 3. The response times of the LS-FST, LS-Taylor, and LS-Sinc with and without the proposed optimization under harmonic phase step test.

Order	LS-FST		LS-Taylor		LS-Sinc
h	Original	Optimized	Original	Optimized	Original	Optimized
2	19.7	32.1	19.7	33.1	19.7	32.3
3	19.2	30.5	19.2	32.8	19.2	31.2
4	19.0	28.7	19.0	32.8	19.0	30.0
5	19.1	27.7	19.0	32.5	19.0	28.9
6	19.1	25.9	19.1	32.5	19.1	27.7
7	19.1	24.5	19.1	32.4	19.1	26.5
8	19.1	23.1	19.1	32.4	19.1	25.1
9	19.2	21.8	19.1	32.4	19.1	24.1
10	19.2	20.6	19.2	32.4	19.1	22.9
11	19.3	19.4	19.3	32.7	19.2	22.0
12	19.6	18.4	19.6	33.0	19.5	21.2
13	47.4	18.2	44.7	46.3	46.3	20.8

, the optimized algorithms have a fast response speed, and the most response times are below two nominal cycles, no matter in the harmonic amplitude or phase step test.

5. Discussion

5.1. Accuracy Analysis under Different Window Lengths

This part analyzes the relationship between the time window length and estimation accuracy of the three least-square algorithms with and without the proposed optimization to determine the application environment. The harmonic frequency deviation test signal with a 40 dB [31] white noise is employed in this section. When the time window length increases from three to seven nominal cycles [11], the the estimation errors and the improvement degrees of the proposed optimization are respectively shown in the left and right subplots of Figure 14. As seen in the right subplots, the proposal produces larger CEI reduction under c = 3, 5, and 7 than c = 4 and 6. Wherein, the response latency of c = 3 is the shortest among c = 3, 5, and 7. Therefore, for applications that require fast response and good accuracy, we recommend c = 3.

5.2. Computational Burden Analysis of the Proposed Optimization

Since the proposed optimization is implemented offline, the last Equation in (15) produces the mainly real-time computation for the three least-square algorithms with and without the proposed optimization [32]. This means the six algorithms have comparable time complexity, as shown in

Table 4. The computation analysis of the LS-FST, LS-Taylor, and LS-Sinc with and without the proposed optimization.

Algorithms	Time Complexity	Execution Time (Min–Max ms)
LS-FST (original)	$O\left(2N\right(H-1\left)\right)$	0.425–0.457
LS-FST (optimized)	$O\left(2N\right(H-1\left)\right)$	0.420–0.459
LS-Taylor (original)	$O\left(2N\right(H-1\left)\right)$	0.423–0.464
LS-Taylor (optimized)	$O\left(2N\right(H-1\left)\right)$	0.426–0.459
LS-Sinc (original)	$O\left(2N\right(H-1\left)\right)$	0.428–0.475
LS-Sinc (optimized)	$O\left(2N\right(H-1\left)\right)$	0.426–0.456

. Moreover, Table 4 also shows the minimum and maximum execution times for a single implementation of the six algorithms when the time window increases from three to seven nominal cycles. The test signal has the same settings as the harmonic frequency deviation test, while the fundamental frequency $f_1$ is fixed to $f_0$. Each execution time is the average value of 10,000 running using MATLAB R2022a on a computer with a 3.2 GHz processor and 16 GB RAM. It is observed these algorithms require light online computation, and have the potential for real-time applications.

6. Conclusions

Harmonic phasor estimation has fundamental significance for addressing harmonic-related issues and developing high-level applications. However, the undesired interharmonic tones and the harmonic parameter variance make it a great challenge to estimate harmonic phasors accurately and fast. In this paper, the largest error source for the comprehensive error index (CEI) caused by multi-interference, and harmonic frequency variance is analyzed. By using singular value decomposition, we introduce variables able to reduce the major component of CEI for least-square algorithms. Based on the CEI and variables, an optimization is constructed aiming to minimize the harmonic phasor estimation error CEI under various interference and harmonic frequency changes. Given the recommended window length of three nominal cycles, the mean TVEs of 2–10 harmonics can be reduced by more than half, and the response latencies of most harmonics are less than two nominal cycles for all three least-square algorithms. In the future, a further accuracy improvement for high-order harmonic phasors under harmonic frequency deviation and dynamic conditions will be preferred.