Performance Analysis of Matrix Pencil Method Applied to High-Resolution Measurement of Supraharmonics

Abstract

Supraharmonics emitted by power electronic equipment will cause electromagnetic interference. However, there is no uniform standard for the measurement of supraharmonics. Unlike non-parametric methods, such as Fourier transform, the matrix pencil method has attracted the attention of researchers in many fields because of its super-resolution characteristics. In this paper, an attempt is made to apply the matrix pencil method to the high-resolution measurement of supraharmonics. Noting that parametric methods offer a promising alternative strategy for high time-resolution analysis of supraharmonics, the matrix pencil method is introduced for high-resolution measurements and estimating the supraharmonic components in grid voltage and current signals. By deriving and analyzing the Cramér–Rao bound of the method, two excellent properties of the matrix pencil method are found: robustness to time-varying signals in supraharmonics measurements and accuracy in frequency localization of components with large amplitudes. Meanwhile, the effectiveness of the matrix pencil method when applied to the high-resolution estimation and measurement of supraharmonics in grid signals is demonstrated by measuring, estimating, and analyzing grid voltage and current signals containing different typical supraharmonics characteristics. This parametric approach opens up a new avenue that allows more signal information to be obtained while maintaining high accuracy, thus enabling the analysis of supraharmonics with higher time resolution.

1. Introduction

As more and more renewable energy generation, electric vehicle charging, and non-linear power loads are connected to the grid, the power electronics of electrical equipment continue to increase. The functional circuit of power electronic electrical equipment will be switched at high frequency, resulting in harmonics interference [1]. As a result, it will distort the voltage and current waveform of the grid. As part of the Electromagnetic Compatibility (EMC) directive, harmonics interference with frequencies below 2 kHz has been subject to regulatory standards and related compliance tests for many years. However, with power electronics technology developing rapidly, the harmonics interference within 2–150 kHz (namely supraharmonics) has increased significantly [2].

For a lot of electrical equipment, the supraharmonics interference will cause the fault of electrical equipment, even reducing the life of electrical equipment. Meanwhile, metering equipment, communication equipment, and other infrastructures will also be affected by supraharmonics interference. For instance, when exposed to supraharmonics interference, the metering results of an energy meter will seriously deviate [3]. Supraharmonics also occupy the same bandwidth as power line communication protocols. They can lead to overheating of reactive power compensation capacitor banks and transformers and the failure of protection devices in the power supply system [4].

Although related issues have attracted attention, the current EMC coordination in the frequency range of 2–150 kHz is still not complete, mainly due to the lack of standardized on-site measurement methods for assessing grid disturbance levels.

In the power system, supraharmonics’ characteristics differ from low-frequency harmonics. Therefore, different methods are needed for their measurement and analysis. Compared with harmonics with frequencies below 2 kHz, supraharmonics are usually non-stationary. The amplitude of supraharmonics will change over time on the millisecond scale, which is generally synchronized with the frequency changes of the power system [5]. Given the complexity of supraharmonics, the measurement and analysis of supraharmonics are usually carried out in the frequency and time domains.

For the measurement of supraharmonics, IEC 61000-4-30 Annex C [6] lists three candidate methods [6]. In the first method, the application scope of the non-overlapping way provided in IEC 61000-4-7 Annex B [7] is extended from 9 kHz to 150 kHz. Specifically, IEC 61000-4-7 defines in Annex B a method for measuring aberrations in the frequency range 2–9 kHz based on a windowed discrete Fourier transform (DFT). The non-overlapping rectangular window is used during the analysis time of 200 ms. The generated 5 Hz frequency resolution interval is grouped into 740 final intervals after clustering, and the frequency resolution of each last interval is 200 Hz [7]. In principle, the supraharmonics can be characterized using traditional tools. Still, supraharmonics is essentially different from harmonics with a frequency below 2 kHz, because the frequency of supraharmonics is not necessarily an integer multiple of the power frequency. Therefore, if the traditional measurement methods are simply used to analyze and calculate the supraharmonics, it will even cause serious spectrum leakage. To address this problem, Ref. [8] proposes a desynchronized processing technique (DPT) as a scheme to evaluate high-frequency distortion (HFD) in power systems. The proposal of using a DPT to assess HFD under experimental conditions is described and demonstrated in ref. [9], where the technique simplifies the measurement hardware.

The second method is the 32 equally spaced measurement method, which is a measurement method with a time interval. Specifically, 32 groups of equal-time data blocks are extracted from the original measured signal with 200 ms, and the data blocks have the same interval of 0.5 ms. Then, the DFT is used to analyze the data blocks, respectively. And 32 sets of characteristic spectral data with a frequency resolution of 2 kHz will be obtained. Because this method only uses 8% of the data in the measured signal, it has small and fast calculation advantages. However, its frequency resolution is low, only 2 kHz, and it cannot accurately measure supraharmonics in the grid. The results cannot be directly compared with the measurement results of other 200 Hz frequency resolution methods. Moreover, because this is a measurement method with a time gap, it may also cause some frequency components in the measured signal to be missed.

The third method is based on the CISPR 16-1-2 method. This method performs DFT using overlapping 20 ms measurement intervals with a Lanczos shape, generating a signal spectrum with a frequency resolution of 50 Hz, which is then aggregated or post-processed using a CISPR 16 detector [10,11]. However, CISPR 16-1-2 is a radio broadcast standard that specifies measurement characteristics for laboratory equipment testing and is not a method for evaluating grid distortions. This method has a 90% overlap of measurement intervals and a large amount of measured data, which is unsuitable for real-time measurements of supraharmonics in the grid.

The above analysis shows that although reference measurements for characterizing supraharmonics emissions from electrical equipment in the laboratory are defined in IEC 61000-4-30, there still needs to be mature measurement methods for supraharmonics identification and estimation. Therefore, the SupraEMI sub-project of the European Metrology Project for Innovation and Research (EMPIR) has been working on strictly characterizing supraharmonics of the power grid under laboratory and field conditions. Through more comprehensive and in-depth research, it is hoped to develop a strict and repeatable standardized measurement software algorithm which can be used to determine the actual distortion limit and uncertainty threshold. In addition to the three candidate methods listed in IEC 61000-4-30 Annex C, the researchers also compared the other three representative methods based on non-parametric models, namely the subsampling approach, wavelet packet decomposition (WPD), and compressive sensing [12]. These methods have different characteristics, involving frequency decomposition, frequency resolution, time resolution, etc.

The subsampling approach method uses a set of analog filters to decompose the measured grid signal into ten sub-signals with a bandwidth of 15 kHz so that a power-quality measuring instrument with a lower sampling rate can be used to analyze the supraharmonics [13]. WPD recursively filters and downsamples the input signals until a 200 Hz frequency resolution is achieved in the 2–150 kHz spectral range [14]. The compressive sensing method can increase the frequency resolution of the measured signal from 2 kHz to 200 Hz while maintaining a time resolution of 0.5 ms [15,16]. However, this method needs to estimate the sparsity of supraharmonics in the measured signal.

In recent years, some scholars have tried to introduce parametric methods to measure supraharmonics. Compared with traditional nonparametric methods, parametric methods are relatively slow due to their higher computational intensity and the need for more computational processing power and memory, which makes them less suitable for large-scale, generalized applications. However, the landscape of computational power is rapidly changing, and the rapid development and advancement of AI technology, in particular, has led to the development of the Graphic Processing Unit (GPU) industry [17]. GPUs were initially designed to handle the computations required to render images and graphics in video games and are essentially designed for parallel computing. They contain hundreds or thousands of cores and can perform many computations simultaneously. This makes them uniquely capable of handling the intensive computational load of parametric approaches [18]. GPU-based technology to implement parametric methods for measuring and evaluating supraharmonics in power grids can provide an effective solution for solving high time-resolution supraharmonic analyses, i.e., so that it can take advantage of the high time-resolution and accuracy achievable with parametric methods, and it can be supported by sufficient computational resource requirements for computationally intensive computations. Therefore, it is now possible to effectively manage the computational effort of parametric methods, making them a viable option for realizing high time-resolution measurements and analyses of supraharmonics in power grids.

In ref. [19], the original waveform was first divided into two frequency bands using the discrete wavelet transform. Then, the Estimation of Signal Parameters using the Rotational Invariance Techniques (ESPRIT) method with sliding window correction and the Nuttal Window’s Sliding Window DFT we used to analyze the low and high frequency bands, utilizing the positive features of each method to minimize the drawbacks and integrate their behaviors. A sliding-window wavelet-modified ESPRIT-based method (SWWMEM) was proposed in ref. [20], which is particularly suitable for the spectrum analysis of waveforms with a vast spectrum. It can be successfully applied to detect spectral components in the 0–150 kHz range introduced in distributed power plants (e.g., wind and photovoltaic systems) and end-user equipment connected to the grid via static converters (e.g., fluorescent lamps).

The matrix pencil method, i.e., a typical parametric method, was proposed by Hua in the early 1990s and introduced into digital signal processing [21,22]. The matrix pencil method belongs to the subspace decomposition algorithm. It can fully use the measured signal information contained in the signal singular value matrix and directly use the orthogonal characteristics of the signal subspace to construct the spectral peak. In terms of estimating the number of frequency components, the matrix pencil method shows better statistical characteristics than the classical non-parametric method, and it does not need to predict the sparsity of frequency components in advance. Some researchers have discussed the above excellent characteristics of the matrix pencil method and have applied it to evaluate harmonics and interharmonics with frequencies below 2 kHz, which achieved good results [23]. The use of the matrix pencil method to measure supraharmonics was proposed in ref. [24], and, based on the analysis that colored noise does exist in the power grid, a supraharmonic high time-resolution measurement method based on matrix pencil with high-order mixed mean cumulant was proposed to improve the estimation accuracy of supraharmonics under colored noise conditions.

However, up to now, in the process of exploring and applying the matrix pencil method, there have been few studies and discussions on the essence of its excellent performance in estimating the different frequency components. Meanwhile, there are few reports on its application to the high-resolution measurement and estimation of supraharmonics. Given this, this paper attempts to introduce the matrix pencil method for the high-resolution measurement and estimation of supraharmonics in the grid voltage and current signals, utilizing its advantages of not needing to predict the sparsity of the different frequency components contained in the measured signals in advance, no iterative computation, and a relatively small computational overhead. By deriving and analyzing the Cramér–Rao bound (CRB) of the method [25], its optimal estimation under different frequency resolutions and time resolutions is explored. Moreover, based on the validation of the measured grid signals containing different typical supraharmonics features, we demonstrate the effectiveness of applying the matrix pencil method to the high-resolution estimation and measurement of supraharmonics in grid signals.

2. Matrix Pencil Method

2.1. Mathematical Model of Supraharmonics

The supraharmonics in the grid can be expressed as a discrete mathematical model, which is given by [15,20]

$$s(n)={\displaystyle \sum _{k=1}^K{\alpha }_k{e}^{-{\beta }_{k}n\Delta t}\cos (n{\omega }_k\Delta t+{\theta }_k)},$$

(1)

where s(n) denotes the value of the discrete signal at the n-th sampling point, and n is an integer indicating the discrete time index or the serial number of the sampling point. k represents the number of supraharmonics components. α_k, w_k, and θ_k represent the amplitude, frequency, and initial phase of the k-th supraharmonics components. β_k is the attenuation coefficient of the k-th supraharmonics components. Δt represents the time interval of equal interval sampling for the voltage and current signals of the measured grid signal. Expressing the above mathematical model in exponential form, we have

$$s(n)={\displaystyle \sum _{k=1}^{2K}{R}_k{e}^{p_{k}n\Delta t}},$$

(2)

where R_k = α_ke^±jθ/2 and p_k = −β_k ± jw_k.

Considering the existence of noise u(n), the supraharmonics components in the grid should be more accurately characterized as

$$x(n)=s(n)+u(n).$$

(3)

2.2. Matrix Pencil Method Applied to Supraharmonics Measurements

The basic idea of the matrix pencil method is to construct a particular matrix (Hankel matrix) and use the special relationship between the matrices to estimate the frequency of the supraharmonics by solving the generalized eigenvalues. To explicitly characterize the physical nature of the supraharmonics time-domain transient-signal-containing noise, that is, how much the supraharmonics components are included in it, the discrete sampling data of the measured time-domain transient signal should be first represented in matrix form and decomposed; then, it can be constructed into the structure of the multiplication of multiple matrices. This is the formation of the so-called matrix pencil. Then, the frequency, amplitude, and initial phase of all the supraharmonics components are extracted from the discrete sampling data of the measured signal by performing singular-value decomposition and using the least squares algorithm to perform the parameter estimation, which minimizes the error and suppresses the noise. Finally, the measured supraharmonics signal’s time-domain waveform is reconstructed as the superimposed combination of the individual supraharmonics components.

For the sampling sequence x(n) of the supraharmonics signal, two Hankel matrices X₁ and X₂ of (N − L) × L are constructed,

$${\mathit{X}}_1=\left[\begin{array}{cccc}x\left(0\right)& x\left(1\right)& \cdots & x\left(L-1\right)\\ x\left(1\right)& x\left(2\right)& \cdots & x\left(L\right)\\ ⋮& ⋮& \ddots & ⋮\\ x\left(N-L-1\right)& x\left(N-L\right)& \cdots & x\left(N-2\right)\end{array}\right],$$

(4)

$${\mathit{X}}_2=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& \cdots & x\left(L\right)\\ x\left(2\right)& x\left(3\right)& \cdots & x\left(L+1\right)\\ ⋮& ⋮& \ddots & ⋮\\ x\left(N-L\right)& x\left(N-L+1\right)& \cdots & x\left(N-1\right)\end{array}\right],$$

(5)

where N represents the number of sampling points. L is a matrix pencil parameter, and its proper selection can reduce the influence of noise. The current research results show that N/3 to 2N/3 is better for L [22]. The Hankel matrices X₁ and X₂ can be represented as

$${\mathit{X}}_1={\mathit{Z}}_L{\mathit{P}\mathit{Z}}_R,$$

(6)

$${\mathit{X}}_2={\mathit{Z}}_L{\mathit{PZZ}}_R,$$

(7)

in which

$${\mathit{Z}}_L=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ e^{p_1\Delta t}& e^{p_2\Delta t}& \cdots & e^{p_{2K}\Delta t}\\ ⋮& ⋮& \ddots & ⋮\\ e^{p_1(N-L-1)\Delta t}& e^{p_2(N-L-1)\Delta t}& \cdots & e^{p_{2K}(N-L-1)\Delta t}\end{array}\right],$$

(8)

$${\mathit{Z}}_R=\left[\begin{array}{cccc}1& e^{p_1\Delta t}& \cdots & e^{p_1(L-1)\Delta t}\\ 1& e^{p_2\Delta t}& \cdots & e^{p_2(L-1)\Delta t}\\ ⋮& ⋮& \ddots & ⋮\\ 1& e^{p_{2K}\Delta t}& \cdots & e^{p_{2K}(L-1)\Delta t}\end{array}\right],$$

(9)

$$\mathit{Z}=diag\left[e^{p_1\Delta t},e^{p_2\Delta t},\cdots ,e^{p_{2K}\Delta t}\right],$$

(10)

$$\mathit{P}=diag\left[R_1,R_2,\cdots ,R_{2K}\right].$$

(11)

There is

$${\mathit{X}}_1-{\lambda \mathit{X}}_0={\mathit{Z}}_L\mathit{P}({\mathit{Z}}_0-\lambda \mathit{I}){\mathit{Z}}_R.$$

(12)

Therefore, the generalized eigenvalues of the matrix pencil X₁ − λX₀ contain the number, frequency, and attenuation coefficient of different frequency components in the supraharmonics signal. Thus, the calculation of various frequency components in the supraharmonics signal is transformed into the generalized eigenvalue of the solution (12), that is,

$${\mathit{X}}_1-{\lambda \mathit{X}}_0={\mathit{X}}_0{}^{+}{\mathit{X}}_1,$$

(13)

where X₀⁺ is the pseudoinverse of X₀. This way, the number, frequency, and attenuation coefficient of different frequency components in the supraharmonics signal can be obtained.

After obtaining the number, frequency, and attenuation coefficient of different frequency components in the supraharmonics signal, the amplitude and initial phase of the specific frequency components in the supraharmonics signal can be obtained by solving the following least squares problem, that is,

$$\left[\begin{array}{c}x(0)\\ x(1)\\ ⋮\\ x(N-1)\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ e^{p_1\Delta t}& e^{p_2\Delta t}& \cdots & e^{p_{2K}\Delta t}\\ ⋮& ⋮& \ddots & ⋮\\ e^{p_1(N-1)\Delta t}& e^{p_2(N-1)\Delta t}& \cdots & e^{p_{2K}(N-1)\Delta t}\end{array}\right]\left[\begin{array}{c}{R}_1\\ R_2\\ ⋮\\ R_{2K}\end{array}\right].$$

(14)

3. Estimation Performance of Matrix Pencil Method

The signal estimation theory can be divided into non-parametric and parametric methods. The latter are different from non-parametric methods, such as power spectrum estimation based on DFT, which does not assume that the data obey a specific probability model. The parametric method assumes that the data follow a probability model with a known structure, but some model parameters are unknown. Parameterized estimation is closely related to the identification of system models. Its primary basis is optimization theory: the estimated parameters should be optimal under specific criteria, and how to obtain the optimal parameter estimation.

Evaluating the estimates of any parameter using a set of sampled samples of the measured signal involves a stochastic process, so it is necessary to consider the estimates as a set of random variables. It is not reasonable to address an estimate in isolation, i.e., it is necessary to know its statistical distribution if the estimate’s accuracy is to be analyzed. A good estimate must be as close to the parameter’s real value as possible. This closeness to the real value can be measured using variances in statistics.

The problem of estimating signal parameters was first proposed in ref. [26]. Later, Reference [27] applied statistical theory to estimate the direction of the arrival of plane waves from signals to linear phased arrays. Fisher and Cramer applied estimation theory to estimate signal parameters [25]. In estimation theory and statistical theory, the Cramér–Rao bound represents the lower bound of the variance of the unbiased estimator of the parameters of the deterministic signal (the signal model is determined, but the model parameters are unknown). The variance of any such estimator is at least as high as the reciprocal of the Fisher information. Equivalently, the CRB denotes an upper bound on the precision of the unbiased estimate (the inverse of the variance). At most, any such estimator’s accuracy is the Fisher information [28].

Unbiased estimators that achieve the CRB are called (completely) valid. This estimation method achieves the lowest possible mean square error and is a minimum variance unbiased (MVU) estimator. It should be noted that there is no unbiased estimate reaching the CRB in some cases.

3.1. Fisher Information

In order to explicitly define the CRB of the matrix pencil method, it is first necessary to deeply understand the concept of the Joint Probability Density Function (JPDF). For a random signal containing N samples, such as x₁, ……, x_N, i.e., denoted as vector x = (x₁, ……, x_N), its JPDF can be defined according to its properties.

Consider x as a complex Gaussian random vector; then, its JPDF can be described in detail according to the definition of ref. [29] as

$$f\left(\mathit{x}\right)=\frac{1}{{\pi }^N\det \left(R_x\right)}{e}^{-{\left(x-E_x\right)}^H{R}_x{}^{-1}\left(x-E_x\right)},$$

(15)

where R_x denotes the autocorrelation matrix of x, and det(·) represents the determinant of the matrix. E_x represents expectation. H denotes the complex conjugate transpose.

A complex Gaussian random variable can be described as a variable with a mean and covariance matrix. Its JPDF reflects the statistical relationship between all the variables, which is a central property of the multidimensional Gaussian distribution. The exact mathematical expression and its detailed derivation process are in the relevant sections in ref. [29].

Once the definition of the JPDF is understood, its relationship with the CRB of the matrix pencil method can be further explored. The CRB provides a lower bound for the estimator’s variance that considers the observed data’s statistical properties, thus providing a theoretically optimal performance benchmark for signal processing and parameter estimation. The specific calculation of the CRB in the matrix pencil method and its impact on the results is an essential part of the study in this paper.

Combining Equations (2) and (3), the JPDF of the supraharmonics can be written as

$$f_{\mathit{x}|\mathit{\eta }}\left(\mathit{x}|\mathit{\eta }\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^N}{e}^{-1/\left(2{\sigma }^2\right){\displaystyle {\sum }_{n=0}^{N-1}{\left|x_n-s_n\right|}^2}},$$

(16)

where |η denotes that the JPDF depends on a vector parameter η. σ is the standard deviation of Gaussian noise. When the matrix pencil method is applied to the supraharmonics measurement, for the k-th supraharmonics components, its amplitude, frequency, initial phase, and attenuation coefficient are unknown parameters, so

$${\mathit{\eta }}^T=\left[{\eta }_1^T,\cdots ,{\eta }_K^T\right],$$

(17)

$${\eta }_k={\left[{\alpha }_k,{\omega }_k,{\theta }_k,{\beta }_k\right]}^T.$$

(18)

Under the condition that the real parameter vector η is given, the quality function V of the sample vector x is defined as the partial derivative of the JPDF’s logarithm Inf(x|η) concerning the real parameter vector η, i.e.,

$$V\left(\mathit{x}\right)=\frac{\partial }{\partial \mathit{\eta }}\mathrm{In}f\left(\mathit{x}|\mathit{\eta }\right)=\frac{\frac{\partial }{\partial \mathit{\eta }}f\left(\mathit{x}|\mathit{\eta }\right)}{f\left(\mathit{x}|\mathit{\eta }\right)}.$$

(19)

Since the mean value of the quality function is zero, its variance is equal to the second-order moment of the quality function, that is, var[V(x)] = E{V²(x)}.

The variance of the quality function is called Fisher information, denoted by J(η) and defined as

$$J\left(\mathit{\eta }\right)=E\left\{{\left[\frac{\partial }{\partial \mathit{\eta }}\mathrm{In}f\left(\mathit{x}|\mathit{\eta }\right)\right]}^2\right\}=-E\left\{\left[\frac{{\partial }^2}{\partial \mathit{\eta }\partial \mathit{\eta }}\mathrm{In}f\left(\mathit{x}|\mathit{\eta }\right)\right]\right\}.$$

(20)

3.2. Cramér–Rao Bound

Let x = (x₁, ……, x_N) be the sample vector drawn from the observed data. When the parameter estimate $\widehat{\mathit{\eta }}$ is an unbiased estimate of the actual parameter η, and both its first-order derivative and its second-order derivative are present, the lower bound on the mean square error (MSE) of the estimate $\widehat{\mathit{\eta }}$ can be characterized by the CRB. This bound represents the minimum variance value that can be achieved by any unbiased estimate given the observations.

The CRB is related to the Fisher Information, which captures the effect of the relationship between the data and the parameter on the accuracy of the parameter estimates. Specifically, the CRB is equal to the inverse of the Fisher information, which means that the larger each element of the Fisher Information is, the smaller the corresponding CRB is and the higher the accuracy of the estimate. The mathematical expression is

$$\mathrm{var}\left(\widehat{\mathit{\eta }}\right)=E\left\{{\left(\widehat{\mathit{\eta }}-\mathit{\eta }\right)}^2\right\}\ge \frac{1}{J\left(\mathit{\eta }\right)}.$$

(21)

The Fisher information J(η) in Equation (21) is defined by Equation (20).

To compute the Fisher information, knowing the probability distribution of the observed data is usually necessary. For continuous data, the elements of the Fisher information can be calculated from the first-order and second-order derivatives of the parameters. This calculation reflects the close relationship between the Fisher information and data distribution and provides a theoretical performance metric for parameter estimation.

A sufficient necessary condition for the equal sign of Equation (21) to hold is

$$\frac{\partial }{\partial \mathit{\eta }}\mathrm{In}f\left(\mathit{x}|\mathit{\eta }\right)=P\left(\mathit{\eta }\right)\left(\widehat{\mathit{\eta }}-\mathit{\eta }\right).$$

(22)

where P(η) is a positive function of η and is independent of sample x₁,……, x_N.

For the unbiased estimation, according to the definition of the CRB, if $\widehat{\mathit{\eta }}$ is an unbiased estimation of η, then the variance of each element $\widehat{\alpha }$_j(j = 1, ……, 4K) of η is not less than the corresponding diagonal term in the inverse of the Fisher information matrix, i.e.,

$$\mathrm{var}\left({\widehat{\alpha }}_j\right)\ge {\left[J^{-1}\right]}_{jj}.$$

(23)

where $\widehat{\alpha }$_j is the estimated value of parameter α_j(j = 1, ……, 4K). [J⁻¹]_jj denotes the j-th diagonal element of the inverse of J.

It is proved in Appendix A that the diagonal block of J⁻¹ can be decomposed into

$$J^{k,k}={\sigma }^2{A}_k^{-1}{H}^{k,k}{A}_k^{-1}.$$

(24)

The meaning of the specific symbol is shown in Appendix A.

Further, the following properties can be obtained from the above proof:

(1) For the k-th supraharmonics component, the CRB for its amplitude α_k, frequency w_k, initial phase θ_k, and attenuation coefficient β_k are independent of the amplitudes of the other supraharmonics components.

It is well-known that the amplitude of some components of the supraharmonics may change on a millisecond time scale, which increases the difficulty of supraharmonics measurements. However, this property also implies that, although there may be some supraharmonics components with dynamically changing amplitudes in the acquired supraharmonics signal, the matrix pencil method is still theoretically capable of accurately estimating and measuring other steady-state supraharmonics components in that supraharmonics signal.

This is because the core idea of the matrix pencil method is based on parameter estimation, which can be adapted to the dynamic changes of the signal, especially for the dynamic changes of the amplitude in the signal. Therefore, even when facing supraharmonics signals with dynamically changing characteristics, the matrix pencil method can accurately measure steady-state supraharmonics components through appropriate parameter adjustment and optimization. This characteristic gives the matrix pencil method greater potential and advantages in applying supraharmonics measurements.

(2) For the k-th supraharmonics component, the CRB for the frequency w_k, the initial phase θ_k, and the attenuation coefficient β_k are inversely proportional to the square of the amplitude α_k².

In other words, the accuracy of the measurements of frequency w_k, initial phase θ_k, and attenuation coefficient β_k increases as the amplitude of the supraharmonics component increases. This property implies that the matrix pencil method can estimate and measure more accurately the frequencies of the supraharmonics components of the grid signals with larger amplitudes that require special attention. This is because the parameter estimation of the matrix pencil method depends on the signal’s amplitude, and signals with larger amplitudes will have more minor estimation errors in their parameters, leading to more accurate measurements.

The above two properties show that the matrix pencil method is excellent for measuring supraharmonics in grid signals. This method can handle dynamic changes in signal amplitude and provide accurate frequency estimates, especially for supraharmonics components with large amplitudes. These properties make the matrix pencil method show significant superiority and application prospects for supraharmonics measurements.

4. Test and Results

4.1. Numerical Simulations

To better-understand the matrix pencil method and its estimation limitations, this paper designs several different synthetic signals based on the typical characteristics of supraharmonics emissions and constructs a simulation model of supraharmonics signals. The first synthetic signal used is described as follows:

K = 2—two supraharmonics components are designed;

α₁ = α₂ = 1—the amplitudes of the two supraharmonics components are equal;

θ₁ = θ₂ = 0—the initial phase of the two supraharmonics components is set to 0;

f₁ = 21,100 Hz.

For this signal containing two supraharmonics components, when the frequency f₁ of the first supraharmonics component is fixed, the frequency f₂ of the second supraharmonics component changes, and the sampling time T changes, the CRB of the matrix pencil method for f₁ is calculated according to Equation (23). The results are shown in Figure 1. The SNR of the added Gaussian white noise is 40 dB.

According to the uncertainty principle of the signal, the frequency resolution of the classical spectral estimation is inversely proportional to the total sampling time of the time-domain signal. That is, for a time-domain signal of length N, if the sampling interval is T_s and the sampling frequency is f_s, then the frequency resolution is f_s/N when using DFT for spectral analysis, i.e.,

$$\Delta f=f_s/N=1/(N{T}_s)=1/T.$$

(25)

where Δf is the frequency resolution; T represents the signal duration used for calculation. It can be seen that there is a certain mutual restriction relationship between the frequency resolution and the signal duration used for analysis. According to IEC 61000-4-30, the frequency resolution of the supraharmonics measurement algorithm is 200 Hz. Based on physical constraints, the limited performance of the time resolution of the supraharmonics signals analyzed using non-parametric methods such as DFT is 5 ms. However, in practice, the change rate of supraharmonics may reach the millisecond scale. Therefore, if the frequency resolution of the supraharmonics measurement algorithm is selected according to the recommendation of IEC 61000-4-30, it may not accurately reflect the dynamic change characteristics of the supraharmonics

In Figure 1, the black line represents the theoretical boundary of DFT. It can be seen that as Δf × T becomes smaller, the CRB of the matrix pencil method becomes larger. Suppose the matrix pencil method can meet the measurement requirements of supraharmonics below the theoretical boundary of the non-parametric method. In that case, it means that the frequency and time resolution performance is improved.

To further explore the influence of the difference in phases on the CRB, under the premise of synthetic signal 1, fix the sampling time to T = 5 ms; the difference in phases of the two supraharmonics components is changed, and the obtained results are shown in Figure 2.

It can be seen from Figure 2 that the CRB of f₁ is affected by the difference in phases between two supraharmonics components when the sampling time is constant. Still, the absolute value of this difference is not significant.

To further explore the theoretical limit of the matrix pencil method, based on synthetic signal 1, the frequency interval f₁ − f₂ = 200 Hz (the frequency resolution Δf recommended by IEC 61000-4-30 is 200 Hz) between the two supraharmonics components was fixed. The CRB of the amplitude α1 of the first supraharmonics component was calculated, and the actual error was obtained according to the matrix pencil method. The results are shown in Figure 3.

It can be seen from Figure 3 that the measurement results of the matrix pencil method are close to the CRB, the error is small, and the trend is the same, indicating that the matrix pencil method has better estimation performance. Further, the difference in the phases Δθ = 0 is fixed, and the sampling time T changes. The results are shown in Figure 4.

From Figure 4, it can be found that even in the case of a short sampling time (T = 2 ms), the calculation results of the matrix pencil method still have high accuracy. That is, high time-resolution analysis was achieved.

To be closer to the actual field situation, the performance of the matrix pencil method under the condition that the measured signal contains multiple supraharmonics components was tested. This paper designs a synthetic signal model containing four supraharmonics components. The frequency and amplitude of each supraharmonics component are provided in

Table 1. Values of parameters in the steady-state supraharmonics simulation model.

Number	Frequency/kHz	Amplitude/a.u.
1	11.9	0.89
2	34.1	0.71
3	43.3	0.58
4	101.5	0.49

. The results obtained through calculation and analysis are shown in Figure 5.

It can be seen from Figure 5 that when the synthetic signal contains multiple supraharmonics components, the matrix pencil method is still close to the CRB and has a low error. This means that even if there are numerous supraharmonics components in the measured signal, the matrix pencil method still has good estimation performance. In addition, it can also be observed that for four supraharmonics components with different amplitudes, the CRB of the amplitude error of each supraharmonics component is the same. This is because the percentage definition was used when defining the amplitude error. According to the derivation of 3.2, it can be seen that for a specific supraharmonic component, the CRB of its amplitude α_k is independent of the amplitude of other supraharmonics components, and the results of Figure 4 are consistent with the theoretical derivation.

To further test the matrix pencil method’s performance in frequency resolution, based on the synthetic signal model in Table 1, the 200 ms duration was selected for analysis, and the frequency of each component was changed by 200 Hz every 50 ms. The calculated results are shown in Figure 6.

Figure 6 shows that the frequency change at every 50 ms can be clearly seen from the results reconstructed using the matrix pencil method, which reflects the better frequency resolution performance of the matrix pencil method.

4.2. Physical Experiment Verification

To further analyze the performance of the matrix pencil method for estimating and measuring the supraharmonics components in the grid, two measured signals were selected for physical experiment verification.

Firstly, the output signal of the PV inverter collected in the field by a laboratory was chosen to be analyzed. The switching frequency of the PV inverter was about 10 kHz, and it generated supraharmonics components, which occurred at the switching frequency and its integer multiples. Its working principle was that the DC power passes through the PV inverter using pulse width modulation (PWM) technology, i.e., it was PWM converted and output, then converted into AC power by an LC filter bank, and finally connected to the low-voltage grid after a transformer was transformed. In the experiments, the output signals of a fully loaded PV inverter were digitally recorded using a YOKOGAVA DL850E recorder [30] and a CAE3N current clamp with a sampling frequency of 500 kHz. To pursue higher time resolution, a rectangular window can be applied to the sampled signal to divide the sampled signal into several shorter time signals. For each signal with a short time, the matrix pencil method was used for analysis, and then the results are jointly presented. In this paper, the rectangular window length was set to 0.5 ms, and a sampling signal of 200 ms was intercepted for analysis. The analysis results are shown in Figure 7a,b.

It can be seen from Figure 7a that the supraharmonics frequency components reconstructed using the matrix pencil method occur at an integer multiple of the switching frequency, which is consistent with the actual situation. Moreover, the matrix pencil method can better-reflect the dynamic time-varying characteristics of the supraharmonics components in the measured signal within milliseconds. That is, the output signal of the photovoltaic inverter fluctuates with time due to the switching frequency.

It can also be found from Figure 7b that the frequency fluctuation range is small for the two supraharmonics components with large amplitudes at 10 kHz and 50 kHz, and the measurement results are relatively stable. This is also consistent with property (2) in Section 3.2: for a specific supraharmonics component, the CRB of frequency w_k is inversely proportional to the amplitude square α_k². For the supraharmonics components with a larger amplitude that need more attention, the matrix pencil method can measure their frequency more accurately. This accurate measurement of frequency also provides a basis for the control of supraharmonics emissions.

To further verify the performance of the matrix pencil method, the measured signal output using a wireless electric-vehicle charging station inverter from a power grid company was selected for physical experiment verification. The three-phase power supply voltage was converted into a DC voltage with a rectifier and then converted into an AC signal with a frequency of about 80 kHz using an inverter and transmitted through a wireless transmitter. The wireless receiver received the AC signal and converted it into a DC signal to charge the electric vehicle. A HIOKI MR8875 recorder and CT9692 current clamp collected the voltage signal at the output end of the inverter (i.e., the input end of the wireless transmitter), and the sampling rate was selected as 200 kHz.

IEC 61000-4-7, in Annex B, defines a method for measuring signal distortion in the 2–9 kHz frequency range based on windowed DFT. Specifically, it uses non-overlapping rectangular windows over a 200 ms measurement interval to obtain frequency analysis results for a 5 Hz frequency interval, which are then clustered. As specified in IEC 61000-4-30 Appendix C, the applicable range of the gapless clustering method in IEC 61000-4-7 Appendix B is extended from 9 kHz to 150 kHz and is divided into 740 groups, with a final frequency resolution of 200 Hz for each group.

$$R_f=\sqrt{{\displaystyle \sum _{x=f-\beta /2+1}^{f+\beta /2}G{\left(x\right)}^2}},$$

(26)

where R_f is the post-clustering amplitude (200 Hz frequency resolution), β is the final frequency resolution (β = 200 Hz), and G(x) is the pre-clustering amplitude (5 Hz frequency resolution).

The above signals were analyzed according to the method defined here and the matrix pencil method, and the results are shown in Figure 8a,b.

As can be seen from Figure 8a, the time-varying dynamic characteristics of the supraharmonics signals are challenging to be reflected by the method recommended by the IEC 61000-4-30 standard because of its low time resolution of only 200 ms. The amplitude of the supraharmonics component reconstructed with the method recommended by the IEC 61000-4-30 standard corresponds to the average value of a 200 ms duration, not the dynamically changing amplitude.

The time-frequency characteristic waveform of the recovery result given in Figure 8b shows that the wireless electric-vehicle charging station inverter outputs a signal containing four supraharmonics components, of which the frequency of the highest amplitude component is about 80 kHz. At the same time, the matrix pencil method can also well-reflect the periodic fluctuation of the output signal of the wireless electric-vehicle charging station inverter over time.

Since the actual physical signal has no reference value in the frequency domain, to evaluate the reconstruction effect of the signal, according to Parseval’s theorem, the total energy of the measured signal in the time domain is equal to the total energy in the frequency domain. That is, there should be:

$${\displaystyle \sum _{n=0}^{N-1}{\left|x(n)\right|}^2}={\displaystyle \sum _{k=0}^{N-1}{\left|X(k)\right|}^2},$$

(27)

where x(n) is the sampling sequence. X(k) represents the reconstructed signal corresponding to the time window.

The energy of the sampled signal was calculated in the time domain every 1 ms, and the energy of the reconstructed signal was calculated in the frequency domain. The comparison results are shown in Figure 9.

It can be seen from Figure 9 that the energy of the signal reconstructed based on the matrix pencil method is close to that of the original signal. After calculation, the energy reconstruction error of the matrix pencil method was only 0.51% in the whole sampling period, which indicates that the method has high recovery accuracy. If the IEC 61000-4-7 method is used, the time resolution needs to reach 5 ms under the premise of a frequency resolution of 200 Hz. In contrast, the matrix pencil method can achieve higher time resolution.

5. Conclusions

The wide bandwidth and millisecond scale dynamics of the supraharmonics impose more stringent requirements on the measurement methods. Current nonparametric methods usually have constraints between frequency resolution and time resolution when dealing with these signals. In contrast, parametric methods typically perform better in frequency resolution and time resolution and can better meet the demands of processing supraharmonics signals. However, most of the current research in parametric methods focuses on the application level and tries to solve practical problems. Studies and discussions on the theoretical performance limits of parametric methods and their performance boundaries in specific cases, such as supraharmonics measurements under different conditions, are still very scarce.

This research gap suggests a need to deepen the theoretical study of parametric methods further to reveal their inherent performance and limitations in dealing with complex signals such as supraharmonics. In addition, more empirical studies are also needed to validate the results of the theoretical studies and provide more instructive suggestions for practical applications. Specifically, the performance of parametric methods under different signal conditions needs to be investigated to find the most suitable methods for handling supraharmonics signals. These studies will contribute to a better understanding and utilization of parametric methods for solving practical problems.

In this paper, the possibility of applying the matrix pencil method to supraharmonics high-resolution measurements is explored in depth. By deriving and analyzing its Cramér–Rao bounds, two excellent properties of the matrix pencil algorithm, namely, robustness to time-varying signals in supraharmonics measurements and accuracy in frequency localization of constituents with large amplitudes, are found to demonstrate the excellent properties of this method in handling supraharmonics high-resolution measurement tasks. The correctness of these theoretical analyses is verified through a series of well-designed numerical simulations and actual physical tests, which are strongly supported by these empirical results. It is also found that the matrix pencil method exhibits extremely high signal-recovery accuracy while maintaining high resolution. This finding further highlights the superiority of the matrix pencil method in handling such complex signals, making it a powerful tool in supraharmonics measurements.

Overall, the research in this paper further expands the understanding of applying the matrix pencil method in supraharmonics measurements and provides a valuable theoretical basis and empirical evidence for further optimization of the method.