Improved Sage–Husa Unscented Kalman Filter for Harmonic State Estimation in Distribution Grid

Abstract

In distribution grids with the large-scale integration of renewable energy sources and energy storage systems, power signals are often contaminated with time-varying noise and frequency deviation caused by low-frequency inertia. To achieve an accurate dynamic harmonic state estimate (HSE), a novel method based on an improved Sage–Husa unscented Kalman filter (ISHUKF) is proposed. Considering the frequency deviation, a nonlinear filter model for power signal is proposed, and a UKF is used to address the nonlinear estimation. A Sage–Husa noise estimator is incorporated to enhance the robustness of the UKF-based HSE against the time-varying noise. Additionally, the noise covariance of the Sage–Husa algorithm is modified to ensure the rapid convergence of the estimation. Then, the performance of the proposed method is validated using an IEEE 14-node system. Finally, the method is applied to evaluate the harmonic states of grid-connected inverter faults in real-world scenarios. The simulation and experiment results demonstrate that the proposed method provides an accurate dynamic HSE even in the presence of time-varying noise and frequency deviation.

1. Introduction

Distributed energy sources are being integrated into power systems on an increasing scale. Nonlinear devices, including grid-connected converters, electromagnetic saturated motors, and impact loads, serve as harmonic sources that inject harmonics into the public grid, thereby impairing the normal operation of the power system [1,2,3,4,5]. To mitigate these harmonics, an accurate harmonic state evaluation (HSE) is essential [6,7,8,9].

The current methods for HSE can be broadly categorized into two types: non-parameter HSE methods and parameter-based HSE methods. Non-parameter HSE methods utilize transforms to map harmonics into a specific domain, where harmonic estimation is performed based on the domain coefficients. Fourier transform (FT), wavelet transform (WT), and Hilbert–Huang transform (HHT) are typical non-parameter HSE techniques [10,11,12,13,14,15,16,17,18,19]. FT methods project harmonics into the frequency domain, allowing for easy harmonic estimation through Fourier coefficients. To mitigate the effects of frequency deviation, windowed interpolation FT has been proposed [10,11]. WT methods map harmonics into the wavelet domain, where each harmonic is analyzed in different frequency bands and evaluated separately [12,13,14]. HHT methods, on the other hand, project harmonics into the Hilbert domain and use empirical mode decomposition (EMD) to estimate the harmonics [15,16,17]. Non-parameter HSE methods require a sufficient length of data to achieve accurate harmonic estimation due to the resolution limitations of the estimation process. Longer data lengths provide more accurate results but also lead to longer estimation delays. For example, when using a 0.2 s window for FT-based HSE, the estimation reflects data from 0.2 s prior to the present moment [18]. Consequently, non-parameter HSE methods are not well-suited for dynamic HSE in distribution grids, where real-time or near-real-time estimation is required.

Parameter-based HSE methods involve the construction of models for the harmonics, from which the harmonic state can be evaluated based on the model parameters. Typical parameter-based HSE methods include stochastic model approaches such as autoregressive moving average (AR/ARMA) and Prony, as well as sinusoidal model methods such as the estimation of signal parameters via rotational invariance (ESPRIT), an adaptive linear combiner (ADALINE), and a Kalman filter (KF) [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. The AR/ARMA method utilizes least-squares estimation to identify the parameters of a suitably chosen model for harmonic detection [19,20,21]. AR/ARMA methods can provide accurate HSE even in the presence of transient disturbances. The Prony method extracts harmonic parameters by solving a set of linear equations for the coefficients of the recurrence relation that the signal satisfies. Prony is capable of detecting all harmonic attributes and the damping factor in the presence of oscillations, but it is sensitive to noise interference [22,23]. The ESPRIT method employs a sinusoidal exponential model and decomposes the covariance matrix into the noise and signal subspaces. Once the frequency components are known, the harmonic state can be determined from the eigenvalues of the signal subspace. References [24,25] applies fast Fourier transform to evaluate harmonic order, followed by an exact model-order ESPRIT algorithm for harmonic estimation. Reference [24] introduces a joint sliding-window ESPRIT and FT scheme for harmonic estimation in power systems, where frequency deviation is accounted for by synchronizing the FT window. ESPRIT, combined with a filter bank, is employed in [26]. However, ESPRIT methods are limited to the frequency domain and are primarily suited for stationary HSE applications. The ADALINE method assumes the signal is a linear combination of harmonics and employs an artificial neural network to estimate harmonic attributes [27,28,29,30]. ADALINE is robust against noise and exhibits good convergence when the harmonic model is appropriately constructed. The KF is widely used for parameter estimation in noisy environments [31,32,33,34]. KF methods are resilient to noise and can provide accurate dynamic HSE, provided prior information about the noise is available.

In distribution grids, the harmonic currents injected into the power system exhibit dynamic characteristics due to varying load levels and the power balance of renewable energy sources. Additionally, large-scale grid-connected converters introduce significant time-varying noise to the grid, primarily due to the high-frequency switching associated with power conversion. Moreover, as the penetration of power converters increases, the frequency inertia of the distribution grid diminishes, leading to larger frequency deviations than before. These factors collectively complicate the accurate HSE in practical applications.

To further enhance the accuracy and robustness of dynamic HSE in the presence of time-varying noise and frequency deviations, this paper proposes an HSE approach based on the improved Sage–Husa unscented Kalman filter (ISHUKF). The key merits and contributions of this technique are as follows:

(1)

Frequency deviation is an inevitable source of interference in HSE. This paper designs a nonlinear harmonic model and employs the UKF to achieve accurate HSE despite the presence of frequency deviations, leveraging the capabilities of the nonlinear model.

(2)

Time-varying noise is another factor that impairs the accuracy of HSE. This paper integrates a Sage–Husa noise estimator into the UKF to minimize the noise influence on HSE.

(3)

A novel strategy is introduced to enhance the convergence of the UKF, thereby improving the dynamic performance of HSE.

The structure of the paper is as follows: Section 2 presents the proposed method, Section 3 evaluates the new method through simulations in a distribution grid, Section 4 demonstrates its performance in field harmonic state estimation, and Section 5 concludes the work.

2. The Proposed Method

In this section, the harmonics are non-linearly modeled and the ISHUKF is proposed to process the model for HSE.

2.1. Harmonic Model

To simultaneously estimate the amplitudes and phase angles of each harmonic, the harmonic signal y_k at time k can be expressed as follows:

$$y_k={\displaystyle \sum _{h=1}^N{a}_k^h}\cos ({\omega }_k^{h}k\mathrm{\Delta }t+{\theta }_k^h)$$

(1)

In the equation, h is the harmonic order; N is the total number of harmonics; $\mathrm{\Delta }t$ denotes the sampling interval; ${\omega }_k^h$ represents the frequency of the h-th harmonic, $a_k^h$ is the amplitude of the h-th harmonic; and ${\theta }_k^h$ is the phase angle of the h-th harmonic. The essence of the HSE method lies in estimating the amplitude and phase at the target frequency. However, if the frequency is fixed, the HSE method cannot produce accurate results for harmonics with varying frequencies. Hence, the frequency has to be set as a state variable and be estimated during the HSE:

$$\left\{\begin{array}{l}{x}_k^1=a_k^1\\ x_k^2={\omega }_{k}k\mathrm{\Delta }t+{\theta }_k^1\\ ⋮\\ x_k^{2h-1}=a_k^h\\ x_k^{2h}=h{\omega }_{k}k\mathrm{\Delta }t+{\theta }_k^h\\ ⋮\\ x_k^{2N-1}=a_k^N\\ x_k^{2N}=N{\omega }_{k}k\mathrm{\Delta }t+{\theta }_k^N\\ x_k^{2N+1}={\omega }_k\end{array}\right.$$

(2)

Then, the harmonics can be modeled as

$$\left\{\begin{array}{l}{x}_k=f(x_{k-1})+q_{k-1}\\ y_k=h(x_k)+r_k\end{array}\right.,$$

(3)

where k is sample index, the system Gaussian noise q follows the $N(0,Q)$ distribution, and the measurement Gaussian noise r follows the $N(0,R)$ distribution. Then, the state variables can be obtained by

$$\left[\begin{array}{l}{x}_{k+1}^1\\ x_{k+1}^2\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_{k+1}^{2h-1}\\ x_{k+1}^{2h}\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_{k+1}^{2N-1}\\ x_{k+1}^{2N}\\ x_{k+1}^{2N+1}\end{array}\right]=\left[\begin{array}{l}1\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}1\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }t\\ ⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 1\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}1\dots 0\hspace{0.17em}\hspace{0.17em}h\mathrm{\Delta }t\\ ⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}1\dots 0\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 1\hspace{0.17em}\hspace{0.17em}N\mathrm{\Delta }t\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}1\end{array}\right]\left[\begin{array}{l}{x}_k^1\\ x_k^2\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_k^{2h-1}\\ x_k^{2h}\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_k^{2N-1}\\ x_k^{2N}\\ x_k^{2N+1}\end{array}\right].$$

(4)

Because the sample interval $\mathrm{\Delta }t$ is much smaller than the harmonic number N, the Equation (4) can be rewritten as

$$\left[\begin{array}{l}{x}_{k+1}^1\\ x_{k+1}^2\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_{k+1}^{2N-1}\\ x_{k+1}^{2N}\\ x_{k+1}^{2N+1}\end{array}\right]=\left[\begin{array}{l}1\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}1\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\\ ⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}1\dots 0\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 1\hspace{0.17em}\hspace{0.17em}0\\ 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}0\dots 0\hspace{0.17em}\hspace{0.17em}1\end{array}\right]\left[\begin{array}{l}{x}_k^1\\ x_k^2\\ \hspace{0.17em}\hspace{0.17em}⋮\\ x_k^{2N-1}\\ x_k^{2N}\\ x_k^{2N+1}\end{array}\right].$$

(5)

The state function and measurement function for the proposed method are

$$\left\{\begin{array}{l}f(x_{k+1})=E{x}_k\\ h(x_k)={\displaystyle \sum _{h=1}^N{x}_k^{2h-1}\cos (x_k^{2h})}\end{array}\right.,$$

(6)

where the identity matrix E is the state matrix. In Equation (6), the h(.) represents the nonlinear function of x, which necessitates the use of a nonlinear method for the HSE.

2.2. UKF for HSE

The UKF is a nonlinear version of the KF with the accuracy of a second-order Taylor series approximation [35]. To implement the nonlinear transformation $y=h(x)$, where the state variable x has a mean value of $\overline{\mathit{x}}$ and covariance of $P_x$, a set of $2n+1$ sigma points $\left\{{\chi }_i\right\}$ is constructed as

$${\chi }_i=\left\{\begin{array}{l}\overline{x}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=0\\ \overline{x}+{\left(\sqrt{(n+\lambda )P_x}\right)}_i\hspace{0.17em}\hspace{0.17em}i=1,2,\cdots ,n\\ \overline{x}-{\left(\sqrt{(n+\lambda )P_x}\right)}_i\hspace{0.17em}\hspace{0.17em}i=n+1,\cdots ,2n\end{array}\right..$$

(7)

Then, the corresponding weights for the sigma points are

$$\left\{\begin{array}{l}{W}_{(0)}^m={\displaystyle \frac{\lambda }{n+\lambda }}\\ W_{(0)}^c={\displaystyle \frac{\lambda }{n+\lambda }}+(1+\beta -{\alpha }^2)\\ W_i^m=W_i^c={\displaystyle \frac{1}{2(n+\lambda )}}\hspace{0.17em}\hspace{0.17em}i=1,\cdots ,2n\end{array}\right.,$$

(8)

where the $\lambda ={\alpha }^2(n+\kappa )-n$ is a scaling parameter that controls the distance between the two sigma points. The constant $\alpha$ is a factor that specifies the spread of the sigma points, $\beta$ is a factor to incorporate prior knowledge of the distribution of $\overline{x}$, and $\beta =2$ is represents the nominal distribution. The $\kappa$ is an auxiliary scaling parameter for the spread of the sigma points, which can be set to 0 [35]. Here, the $\alpha$, $\beta$ and $\kappa$ are set as 0.001, 2 and 0, respectively.

The detailed steps of the UKF algorithm are as follows:

(1)

Initialization, k = 0, ${\widehat{x}}_0$ and $P_0$ are

$$\left\{\begin{array}{l}{\widehat{x}}_0=E[x_0]\\ P_0=E\left[\right(x_0-{\widehat{x}}_0\left)\right(x_0-{\widehat{x}}_0{)}^{\mathrm{T}}]\end{array}\right.,$$

(9)

where $E[.]$ is the expect function.

(2)

Time Update

Calculate the Sigma points at time k, and construct the matrix:

$${\chi }_{i,k-1}=\left\{\begin{array}{l}{\widehat{x}}_{k-1}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} ,i=0\\ {\widehat{x}}_{k-1}+{(\sqrt{(n+\lambda )P_{k-1}})}_i\hspace{0.17em}\hspace{0.17em},i=1,2,\dots ,n\\ {\widehat{x}}_{k-1}-{(\sqrt{(n+\lambda )P_{k-1}})}_i\hspace{0.17em}\hspace{0.17em},i=n+1,n+2,\dots ,2n\end{array}\right..$$

(10)

(3)

Prediction

Perform the transformation of the state equation on the set of Sigma points, and then predict the statistical quantities at time k:

$${\widehat{\chi }}_{i,k\left|k-1\right.}=f({\chi }_{i,k-1})+q_{k-1},$$

(11)

$${\widehat{x}}_{k\left|k-1\right.}={\displaystyle \sum _{i=0}^{2n}{W}_i^m}{\widehat{\chi }}_{i,k\left|k-1\right.},$$

(12)

$$P_{k\left|k-1\right.}={\displaystyle \sum _{i=0}^{2n}{W}_i^c}\left({\widehat{\chi }}_{i,k\left|k-1\right.}-{\widehat{x}}_{k\left|k-1\right.}\right)\cdot {\left({\widehat{\chi }}_{i,k\left|k-1\right.}-{\widehat{x}}_{k\left|k-1\right.}\right)}^{\mathrm{T}}+Q_{k-1}.$$

(13)

Perform the nonlinear transformation $h(\cdot )$ on the ${\chi }_{i,k-1}$ in the measurement equation, and then obtain the statistic $y_{i,k\left|k-1\right.}$ as

$$y_{i,k\left|k-1\right.}=h({\chi }_{i,k-1})+r_{k-1}.$$

(14)

Then, use the statistic $y_{i,k\left|k-1\right.}$ to recover the measurement prediction ${\widehat{y}}_{k\left|k-1\right.}$ as

$${\widehat{y}}_{k\left|k-1\right.}={\displaystyle \sum _{i=0}^{2n}{W}_i^m{y}_{i,k\left|k-1\right.}}.$$

(15)

and the auto-covariance $P_{yy}$ and cross-covariance $P_{xy}$ are calculated as follows:

$$P_{yy}={\displaystyle \sum _{i=0}^{2n}{W}_i^c(y_{i,k\left|k-1\right.}-{\widehat{y}}_{k\left|k-1\right.})(y_{i,k\left|k-1\right.}-{\widehat{y}}_{k\left|k-1\right.}){)}^{\mathrm{T}}}+R_{k-1},$$

(16)

$$P_{xy}={\displaystyle \sum _{i=0}^{2n}{W}_i^c({\chi }_{i,k-1}}-{\widehat{x}}_{k\left|k-1\right.}){(y_{i,k\left|k-1\right.}-{\widehat{y}}_{k\left|k-1\right.})}^{\mathrm{T}}.$$

(17)

(4)

Measurement Update

Calculate the filter gain $K$:

$$K_k=P_{xy}{P}_{yy}^{-1}.$$

(18)

Then, the updated state estimate $x_k$ and covariance estimate $P_k$ at time step k are

$${\widehat{x}}_{k\left|k\right.}={\widehat{x}}_{k\left|k-1\right.}+K_k(y_k-{\widehat{y}}_{k\left|k-1\right.}),$$

(19)

$$P_k=P_{k\left|k-1\right.}-K_k{P}_{yy}{K}_k^{\mathrm{T}}.$$

(20)

According to the filter model, the h-order harmonic state at time k is evaluated as

$$\left\{\begin{array}{l}{a}_k^h=x_k^{2h-1}\\ {\theta }_k^h=x_k^{2h}-h{x}_k^{2N+1}k\mathrm{\Delta }t\end{array}\right..$$

(21)

2.3. Sage–Husa Noise Estimation

In traditional UKF methods, the noise covariance Q and R in Equation (15) and Equation (18) are fixed during the estimation process. This limitation prevents the accurate stratification of HSE in distribution grids, especially due to the presence of dynamic noise. The Sage–Husa estimation method is introduced to adjust the covariance values according to the noise level [36], thereby improving the accuracy of HSE when using the UKF.

(1)

The covariance of system noise is evaluated as

$$\begin{array}{ll}{\widehat{Q}}_k& =(1-d_k){\widehat{Q}}_{k-1}+d_k(K_k{e}_k{e}_k^T{K}_k^T+P_k-E_{k-1}{P}_{k-1}{E_{k-1}}^T)\\ & =(1-d_k){\widehat{Q}}_{k-1}+d_k(K_k{e}_k{e}_k^T{K}_k^T+P_k-P_{k-1})\end{array},$$

(22)

where $e_k=y_k-{\widehat{y}}_{k\left|k-1\right.}$, and $d_k=(1-b)/(1-b^{k+1})$ is the weighted coefficient and b ⋲ [0.95, 0.99] is the forgetting factor.

(2)

The covariance of measurement noise is estimated as

$${\widehat{R}}_k=(1-d_k){\widehat{R}}_{k-1}+d_k{e}_k{e}_k^{\mathrm{T}}.$$

(23)

2.4. Improved Sage–Husa Noise Estimation

In Equation (22), in the system with dynamic noise, the estimation of ${\widehat{Q}}_k$ is prone to filter divergence because it tends to fall into a nonpositive semidefinite. Then, adjustment factor $\mu$ is introduced to modify the estimation of ${\widehat{Q}}_k$ without losing critical information from the original Sage–Husa estimation [37]. Then the ${\widehat{Q}}_k$ is given by

$${\widehat{Q}}_k=(1-d_k){\widehat{Q}}_{k-1}+d_k(K_k{e}_k{e}_k^T{K}_k^T+{({\mu }_k)}^p(P_k-P_{k-1})),$$

(24)

where p is initialized as 0, and the adjustment factor $\mu$ is calculated as

$${\mu }_k=\exp (-{\displaystyle \frac{\left|tr(P_k-P_{k-1})\right|}{\left|K_k{e}_k{e}_k^T{K}_k^T\right|}}),$$

(25)

where tr(.) is the trace operation. At every step k, the ${\widehat{Q}}_k$ is monitored. If the $\min (eig({\widehat{Q}}_k))<0$, the ${\widehat{Q}}_k$ is not a positive semidefinite, and then p = p + 1 in Equation (26) until the $\min (eig({\widehat{Q}}_k))>0$, which ensures that the ${\widehat{Q}}_k$ is a positive semidefinite. Here, the min(.) returns the minimum value of (.) and eig(.) returns the eigenvalue vector of (.).

2.5. Harmonic State Estimation by Improved SHUKF (ISHUKF)

The HSE process based on the new method is as follows:

(1)

Sample the power signals.

(2)

Use the UKF to process the power signals and estimate the harmonics state by the state variable vector x according to Equation (21).

(3)

Use Sage–Husa noise to estimate covariances of system noise Q and measurement noise R, respectively, for the UKF iteration.

(4)

Monitor the Q, and use Equation (27) to make sure that the ${\widehat{Q}}_k$ is a positive semidefinite.

(5)

Refresh the Q and R, then move forward to Step (2).

(6)

Stop the whole process if the HSE is over.

The flowchart of the new method is shown in Figure 1.

The computational complexity of the UKF is O((2N + 1)³). Similarly, the one of Equation (26) is also O((2N + 1)³). Therefore, the computational complexity of the new method is O(2(2N + 1)³), which is directly related to the number of state variables of the filter model. In application, the size of the filter model can be calibrated for the balance of interested HSE accuracy and computation consumption.

Pseudocode of the new method

1. Build the nonlinear filter model for the interested harmonics (Equation (6)). 2. Initialize the parameters: r, q, P, $\alpha$, $\beta$, $\kappa$ and b. 3. Obtain power signal $y_k$. 4. Generate sigma points for the UKF estimation (Equations (7) and (8)). 5. Perform UKF to estimate $y_k$ (Equation (7) to Equation (21)). 6. Monitor the noise by improved Sage–Husa estimator (Equation (23) to Equation (25)). 7. Return to 5 to modify the Q and R of the UKF.

3. Simulation Analysis

In this section, the new method is applied for HSE in a simulated distribution grid system. The settings of the new method are given in

Table 1. Parameter settings of the new method.

UKF						Sage–Husa
r	q	P	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\kappa }$	b
0.01	0.01	0	0.001	2	0	0.98

. Dynamic HSE methods such as ADALINE and KF are used for the ISHUKF comparisons. The sampling rate is 4 kHz.

3.1. Harmonic State Estimation in Distribution Grid Model

To verify the performance of the ISHUKF for HSE, this paper constructs the IEEE 14-node system for verification, as shown in Figure 2. This system has generators connected to nodes 1, 2, 3, 6, and 8 for power supply. The specific parameters of the system components such as generators, transformers, buses, and branches can be found in [38]. Apart from the harmonic sources at nodes 10 and 13, the rest are non-harmonic source buses. The harmonic sources are modeled as current sources, and their harmonic spectrum is shown in

Table 2. Harmonic source spectrum table.

Harmonic Order	Amplitude/p.u.	Phase/°
1	1.0000	0.00
5	0.1824	−55.68
7	0.1190	−84.11
11	0.0573	−143.56
13	0.0401	−175.58

. Measurement devices are installed at nodes 2, 4, 5, 6, 7, 8, 11, and 12 to capture the harmonic at each node.

In the simulation, the harmonic impedance matrix of the system is established using the BIBC (bus injection current-branch current)–BCBV (branch current-bus voltage) method, as described in [39]. The ISHUKF is applied for the estimation of harmonic currents injected into the grid with simulated noise.

To evaluate the accuracy of HSE, the root-mean-square error (RMSE) is introduced as the evaluation metric. The RMSE can measure the deviation between the estimated and true values of the harmonic current magnitude and phase angle, and the RMSE for the h-order harmonic evaluation is given by

$$RMS{E}_h=\sqrt{{\displaystyle \frac{1}{L}}{\displaystyle \sum _{l=1}^L{({\widehat{y}}_{l,h}-y_{l,h})}^2}},$$

(26)

where l is the sample index and L is the data length. The ${\widehat{y}}_{l,h}$ and $y_{l,h}$ represent the estimated and true values of the magnitude or the phase angle of the h-th harmonic current respectively. A smaller RMSE indicates a more accurate HSE.

To evaluate the HSE in the presence of noise, the standard deviation (STD) is used as a metric. A method with smaller STD is less affected by noise, indicating a smoother result with reduced noise contamination. The standard deviation is given by

$$ST{D}_h=\sqrt{{\displaystyle \frac{1}{L}}{\displaystyle \sum _{l=1}^L{({\widehat{y}}_{l,h}-{\overline{y}}_{l,h})}^2}},$$

(27)

where ${\overline{y}}_{l,h}$ represents the mean values of the magnitude of the h-th harmonic current.

3.2. Harmonic State Estimation with Dynamic Noise

To verify the estimation accuracy of the ISHUKF with dynamic noise interference, the 0.2 s harmonic current is divided into five segments, each lasting 0.04 s. The signal-to-noise ratios (SNRs) of the segments are 50 dB, 45 dB, 40 dB, 35 dB, and 30 dB, respectively, from the first to the last segment. The noise is Gaussian with a mean of zero, and the system frequency is 50 Hz. The harmonic current and the noise affecting the harmonic are shown, respectively, in Figure 3.

The RMSE and STD comparisons of the three methods are presented in

Table 3. The RMSE comparison of the HSE by the three methods.

Harmonic Order	Harmonic Parameters	Methods
		ADALINE		KF		ISHUKF
		RMSE	STD	RMSE	STD	RMSE	STD
5th	Amplitude (p.u.)	0.0096	0.0084	0.0091	0.0087	0.0073	0.0072
	Phase (degree)	2.717	2.672	2.697	2.681	2.189	2.220
7th	Amplitude (p.u.)	0.0078	0.0076	0.0073	0.0077	0.0053	0.0053
	Phase (degree)	4.465	4.064	4.173	4.380	3.358	3.390
11th	Amplitude (p.u.)	0.0067	0.0066	0.0066	0.0069	0.0038	0.0039
	Phase (degree)	10.767	7.582	9.42	7.749	6.001	6.013
13th	Amplitude (p.u.)	0.0066	0.0067	0.0066	0.0066	0.0034	0.0035
	Phase (degree)	10.381	9.588	10.885	10.365	4.854	4.064

. Both the ADALINE and KF are freed linear estimators, with harmonic frequencies set based on the system frequency. As a result, the HSE errors arise solely from the noise. The ADALINE and KF yield similar RMSE and STD values for harmonic estimation since both methods rely on the least-mean-squares rule. However, because the parameters of ADALINE and KF, such as the learning rate and noise variance, are fixed, the estimation of the two methods cannot be adjusted to the change of noise. In contrast, the ISHUKF uses the Sage–Husa algorithm to estimate the noise and tunes the parameter of the UKF at every sample of the current, which provides more accuracy (smaller RMSE) and is less noise-affected (smaller STD) than the two methods. Moreover, the improved Q calculation helps to prevent filter divergence when noise conditions change, further contributing to the superior performance of the proposed method.

The performance comparison of the three methods for the 5th, 7th, 11th and 13th harmonics are shown in the Figure 4, Figure 5, Figure 6, and Figure 7, respectively. From these figures, it is evident that the estimation curve of the ISHUKF exhibits less noise. This is primarily due to the underlying principle of the UKF. When the noise is significant, the UKF places greater trust in the estimated state variables (${\widehat{x}}_{k\left|k-1\right.}$) in Equation (19), and the estimation ${\widehat{x}}_{k\left|k\right.}$ predominantly relies on the filter’s prediction. Conversely, when the noise is low, the UKF relies more on the measurement ($y_k-{\widehat{y}}_{k\left|k-1\right.}$), and the estimation is primarily influenced by the measurement. The Kalman gain ($K_k$) controls the weighting of the measurement. Because the Kalman gain is adjusted at each time step k, based on the noise covariance Q and R, the ISHUKF is less affected by noise compared to the other two methods.

4. Experiment

To test the performance of the method proposed in this paper on a hardware platform, we have established a setup consisting of the following components, as shown in Figure 8. The platform includes a 600 V DC power supply, a three-phase two-level inverter (simulating the harmonic source), a 40 kW programmable three-phase AC power source (simulating the distribution grid, with adjustable frequency), a programmable load, a rapid control prototyping (RCP) unit, and a computer with CPU intel i7-8750H memory DDR4 16G (for HSE and result display) (ThinkPad, Beijing, China). The computer downloads the program to the inverter controller via the RCP unit. Meanwhile, the signal conditioning circuit collects the three-phase current signals, which are then transmitted to the computer through the RCP unit for the HSE and display. The current sampling rate is 4 kHz, and the filtering capacitor and inductor are 50 uF and 1 mH, respectively. When the inverter switches are open-circuited, all three phase currents become distorted due to the condition where the sum of the currents equals zero, leading to the injection of significant harmonics into the grid. Since the harmonic orders of the currents are low during the fault condition, the DC component, fundamental component, second harmonic, and third harmonic are modeled. The computational complexity of this process is O(2 × 11³).

Table 4. Key parameters of the experimental settings.

System frequency (Hz)	50 ± 0.5
Injected current (A)	1 to 7.5
Connected grid line voltage (V)	380
Sampling rate (Hz)	4000
Interested harmonic order	Second and third

presents the key parameters of the experiments.

4.1. HSE for Single Switch Open-Circuited Harmonics

This experiment evaluates the new method for dynamic HSE under frequency deviation. The load current is initially 5 A and increases to 7.5 A at 0.3 s, then decreases back to 5 A at 0.5 s. The upper arm switch of Phase A is open-circuited at 0.4 s. The system operates at 49.5 Hz, controlled by the programmable AC source. The current signals of the three-phase inverter and the estimated frequency by ISHUKF are shown in Figure 9. Prior to the fault, the harmonics of the current are mainly from switches by pulse width modulation and very small. After the fault occurs, intense harmonics appear in the faulty phase, as half of the current waveform is missing [40]. Also, because the sum of the three phase currents is zero, the currents of the unfaulty phases also exhibit harmonics. The runtime of the ISHUKF is 0.08 s, which is much less than 0.6 s, demonstrating the real-time capability of the proposed method.

The estimation comparisons for the second harmonic and third harmonic of each phase are shown in Figure 10, Figure 11, and Figure 12, respectively, because the two harmonic and third harmonic are the two largest harmonics in the three-phase currents.

In the Figure 10, Figure 11 and Figure 12, both ADALINE and KF exhibit oscillations for all the estimation of all harmonics, as these two methods are designed for a 50 Hz system. When there is frequency deviation, the harmonic frequencies do not align with the model assumptions of these methods. As a result, even in the absence of harmonics (before the fault), the outputs from both methods show oscillatory behavior. In contrast, the ISHUKF incorporates frequency as part of the state variables and estimates the frequency at each step k (See Figure 9). The estimated frequency is then used for the HSE, which ensures that frequency deviations do not affect the performance of the ISHUKF.

The ISHUKF requires approximately 0.01 s (half a cycle) for the dynamic HSE, which is close to that of ADALINE and about 0.003 s faster than that of the KF. Both ADALINE and KF use least-mean-squares estimation. The weight updates in ADALINE are managed by a neural network, making it robust against noise interference. On the other hand, the state variables of the KF are prone to divergence when the filter model is large and with noise. Hence, the KF takes a longer time for the HSE to keep the filter convergence. In contrast, the ISHUKF is less affected by the noise and provides fast HSE, similar to ADALINE.

4.2. HSE for Multiple Switch Open-Circuited Harmonics

This experiment aims to evaluate the performance of the proposed method for steady HSE under conditions of low SNR. In this scenario, the upper arms of Phase A and Phase B are open-circuited, with the current set to 1 A. Since two phases are affected by faults, the direct component of the fault current is higher than that of a single-phase fault. The peak value of the fault current is approximately 2.6 A, which is about half of the value observed in Section 4.1. While the noise remains nearly unchanged, the currents in this case are reduced. Figure 13 shows the three-phase currents.

Table 5. The standard deviation comparison of the HSE by the three methods.

Harmonic Order	ADALINE	KF	ISHUKF
Second	0.0046	0.0056	0.0039
Third	0.0041	0.0046	0.0035

presents STDs of the three methods for HSE, and Figure 14 shows the detailed performance of the Phase A HSE.

In Table 4, both ADALINE and KF show similar estimation results, as both methods are based on least-mean-squares estimation. However, due to the presence of noise, the true harmonic amplitudes are distorted, leading to larger standard deviations (STDs) in the estimations. ADALINE adjusts the weights of the neurons for HSE using gradient-based optimization, while KF uses Kalman gain, which is more susceptible to saturation, to perform a similar function. As a result, the ADALINE method responds more quickly to noise than KF. In contrast, the ISHUKF avoids Kalman gain saturation because the process and measurement noise covariances (Q and R) are updated at each step k. Additionally, by estimating the noise, the ISHUKF is better able to identify the true harmonic values amidst the noise, leading to more accurate HSE with smaller STDs compared to ADALINE and KF.

5. Conclusions

In order to accurately perform the dynamic harmonic state estimation of the power system under the influence of time-varying noise and frequency deviation, this paper proposes a dynamic HSE based on an improved Sage–Husa UKF. Because the frequency is estimated, the filter model of HSE becomes nonlinear, and the UKF is used to handle this nonlinearity. The Sage–Husa noise estimator is integrated to evaluate the noise covariance, ensuring the accuracy of the UKF estimation even in the presence of noise. Furthermore, the improved algorithm maintains the system noise covariance as a positive semidefinite, allowing the UKF to remain accurate despite time-varying noise. The proposed method is validated through simulations of the IEEE 14-node system for HSE with time-varying noise, as well as through experiments on a faulty grid-connected inverter for field HSE under frequency deviation, load changes and noise.

The proposed method can be flexibly adapted for harmonic estimation by adjusting the size of the filter model and can also be applied in various areas, such as harmonic source localization and power quality evaluation. Additionally, many types of power equipment generate harmonics under faulty conditions. The new method also holds potential for fault detection applications.