Harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing

Abstract

The present invention discloses a harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing, comprising the following steps: (1) constructing a sensor receiving signal into a parallel factor trilinear model; (2) compressing the three-dimensional matrix of the received signal according to the compressed sensing theory; (3) using a trilinear alternating least squares algorithm to perform cyclic iterative decomposition on the compressed model until convergence to obtain an estimate of the parameter matrix; (4) constructing a sparse recovery problem in frequency estimation, and sequentially using an orthogonal matching pursuit method to obtain frequency estimates. The present invention is applicable to harmonic and interharmonic signals in power systems and is easy to implement in actual use; the algorithm of the present invention has a small amount of calculation and has a low capacity requirement for data storage; the algorithm of the present invention can obtain paired frequency estimates with the help of a parallel factor model without the need for additional pairing.

Description

Technical Field

[0001] The present invention relates to the field of harmonic measurement in power systems, and in particular to a harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing. Background Art

[0002] With the widespread application of nonlinear devices such as power electronics in power systems, a large number of nonlinear loads are connected to the power grid, resulting in increasingly serious harmonic and interharmonic pollution in the power grid. Various faults and accidents caused by harmonics and interharmonics continue to occur, seriously affecting the safety of the power grid and the quality of power. Therefore, it is necessary to control harmonics and interharmonics, and accurate and effective harmonic and interharmonic parameter estimation is the premise and important guarantee for harmonic and interharmonic control.

[0003] The classic power system harmonic and interharmonic estimation method is mainly based on Fourier transform. This type of algorithm has the advantages of fast operation speed and easy hardware implementation, but the algorithm has two disadvantages: first, the frequency resolution is limited, and only integer harmonic parameter estimation can be achieved, and interharmonic estimation cannot be achieved; second, the algorithm requires synchronous sampling, but interharmonics often exist in actual power grids, and it is difficult to achieve synchronous sampling. When sampling asynchronously, the spectra of each harmonic and interharmonic will interfere with each other, causing serious frequency leakage and fence effect, resulting in failure of harmonic parameter estimation. Windowing and interpolation are effective methods to overcome Fourier transform spectrum leakage, but their resolution is limited. Wavelet transform can also overcome the shortcomings of Fourier transform, but when estimating interharmonics, there is a mixing phenomenon, resulting in low resolution and inability to identify similar frequencies.

[0004] In order to solve the inherent defects of the classical Fourier transform, modern spectral estimation theory has been applied to the estimation of harmonic and interharmonic parameters in power systems. Modern spectral estimation theory is developed on the concept of linear algebraic subspace. The basic principle is to use the eigenvalue decomposition theory of the signal autocorrelation matrix to decompose the feature space of the received signal into signal subspace and noise subspace, and use their correlation characteristics to achieve super-resolution parameter estimation. The parallel factor (PARAFAC) method was first proposed in physiology for the decomposition of multidimensional matrices. In recent years, this method has been gradually applied to signal parameter estimation. The parallel factor method usually requires a cyclic iterative decomposition of the received signal model to obtain an estimate of the parameter matrix. This process is generally complex, especially when the number of received signal samples is large. Compressed sensing (CS) is a technology used to analyze high-dimensional data. With the help of compressed sensing ideas, high-dimensional data can be compressed to reduce the sampling data and reduce the data storage space. Therefore, combining compressed sensing theory with the parallel factor model can effectively reduce the computational complexity of the parallel factor method and the requirements for data storage space. Summary of the invention

[0005] 1. Technical problems solved:

[0006] The technical problem to be solved by the present invention is to provide a harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing, which can obtain frequency estimation with performance close to that of the traditional parallel factor method, and at the same time its computational complexity is lower than that of the traditional parallel factor method.

[0007] (II) Technical solution:

[0008] In order to solve the above technical problems, the present invention provides a harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing, comprising the following steps:

[0009] (1) The sensor receiving signal is constructed as a parallel factor trilinear model;

[0010] (2) According to the compressed sensing theory, the three-dimensional matrix of the received signal is compressed;

[0011] (3) Using the trilinear alternating least squares algorithm, the compressed model is iteratively decomposed until convergence to obtain an estimate of the parameter matrix;

[0012] (4) Construct the sparse recovery problem in frequency estimation and use the orthogonal matching pursuit method to obtain frequency estimation.

[0013] Preferably, step (1) is specifically:

[0014] The power system frequency signal (voltage or current) containing noise, power frequency, harmonics and interharmonic components at the receiving end can be expressed as

[0015]

[0016] Where: M is the number of sinusoidal components including fundamental wave, harmonics and interharmonics; α m is the amplitude of the mth sinusoidal component; ω m is the angular frequency of the mth sinusoidal component; is the phase of the mth sinusoidal component; e(t) is the noise signal.

[0017] Using Euler's formula, it can be transformed into

[0018]

[0019] Where: A m is the amplitude of the harmonic signal, which is a complex constant; ω m is the normalized frequency to be estimated; is the normalized initial phase, uniformly distributed in the interval (-π,π); w is an uncorrelated complex Gaussian white noise with zero mean. If the noise component is not considered, it can be simplified to

[0020]

[0021] In the formula

[0022] In order to estimate the frequency of the received signal, a multi-delay output signal model of the receiving end is constructed. The P (P>2) discrete sinusoidal signals at the receiving end can be expressed as

[0023]

[0024] Where: {D P ,p=1,2,…,p-1} is a multi-delay output; {w p (n),p=0,1,…,P-1} is an uncorrelated complex Gaussian white noise with zero mean.

[0025] According to the received signal, the following vector group is constructed:

[0026]

[0027] Where: Q is the vector length; [·] T represents the matrix transpose. Without considering the influence of noise, the signal in the formula can be expressed as

[0028]

[0029] Where:

[0030]

[0031]

[0032]

[0033] A is a frequency estimation matrix with Vandermonde characteristics. Define the delay matrix H and the source matrix S:

[0034]

[0035]

[0036] Construct a set of matrices:

[0037]

[0038] It can be expressed as

[0039] Z p =ADiag p+1(H)S, p = 0, 1, …, P - 1

[0040] where: Diag i (·) is a diagonal matrix constructed using the i-th row of the matrix; Z P represents the p-th slice in the time-delay direction.

[0041] Considering the influence of noise, the observation model becomes

[0042]

[0043] where W P is the received noise, which is a Q×L matrix.

[0044] By rearranging the equation, a trilinear model can be obtained:

[0045]

[0046] p = 1, 2, …, P; l = 1, 2, …, L

[0047] where: a q,m is the (q, m) element of matrix A, h p,m is the (p, m) element of matrix H, s m,l is the (m, l) element of matrix S. Utilizing the symmetry of the trilinear model, the observation models in other directions can be obtained by rearrangement:

[0048]

[0049]

[0050] Preferably, step (2) is specifically:

[0051] To reduce the complexity of the parallel factor decomposition algorithm, the three-dimensional matrix of the received data is compressed into a smaller matrix where Q' < Q, P' < P, L' < L. The compressed matrix is defined as The above compressed matrix needs to satisfy the RIP property and can be obtained through Tucker3 decomposition or partial random matrices such as Gaussian random matrices, Bernoulli random matrices, or local Fourier matrices.

[0052] Then the compressed three-dimensional matrix can be expressed as

[0053] Z' p' = A'Diag p'+1 (H')S', p' = 0, 1, …, P' - 1

[0054] where A' = U H A, H' = V HH, S' = SW.

[0055] Preferably, step (3) is specifically:

[0056] According to the trilinear alternating least squares algorithm (TALS), the least squares update of the parameter matrices A', H', S' is

[0057]

[0058]

[0059]

[0060] in is a compressed signal slice containing noise; using the above updating method, the final estimated values ​​of the parameter matrices A', H', S' are obtained by cyclic iteration.

[0061] Preferably, step (4) is specifically:

[0062] Define the mth column of matrix A' as a'(ω m ), then Where a(ω m ) is the mth column of matrix A, is the scale factor; definition as a sampling grid of all frequencies with R>>M, used to construct a complete dictionary for frequency estimation if Rewrite the above formula as where e p is a P×1 sparse matrix. Using the L0 norm constraint, we can get ω m The estimate is

[0063]

[0064] (III) Beneficial effects

[0065] Compared with the prior art, the beneficial effects of the present invention are:

[0066] (1) The present invention is applicable to harmonic and interharmonic signals in power systems and is easy to implement in practical use;

[0067] (2) The algorithm of the present invention has a small amount of calculation and requires a low capacity for data storage;

[0068] (3) The algorithm of the present invention can obtain frequency estimation of pairing with the help of parallel factor model without the need for additional pairing. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 It is a flow chart of the present invention.

[0070] Figure 2 It is a schematic diagram of compression of the receiving signal model of the present invention.

[0071] Figure 3 Schematic diagram comparing the complexity of the parallel factor algorithm of the present invention and the traditional one.

[0072] Figure 4 This is a distribution diagram of the frequency estimation results of the present invention.

[0073] Figure 5 This is a comparison chart of the frequency estimation performance of the present invention and two commonly used algorithms.

[0074] Figure 6 This is a comparison chart of the frequency estimation performance of the present invention under different sampling numbers.

[0075] Figure 7 This is a comparison chart of the frequency estimation performance of the present invention under different numbers of sensors.

[0076] Figure 8 This is a distribution diagram of dynamic frequency estimation results of the present invention. DETAILED DESCRIPTION

[0077] The present invention will be further explained below in conjunction with the accompanying drawings.

[0078] The present invention proposes a method for estimating the frequency of harmonics and interharmonics in a power system based on compressed sensing parallel factor decomposition. The method first performs multi-stage delay reception on the signal in order to estimate the frequency of the received signal, and constructs the multi-stage delay output signal into a parallel factor trilinear model. In order to reduce the complexity of decomposing the model, the algorithm combines the compressed sensing theory to compress the model, and then decomposes the model using a trilinear alternating least squares algorithm to obtain an estimate of the signal parameter matrix, and obtains the frequency estimates of the fundamental and harmonic waves from the obtained parameter matrix. During the decomposition process, the parameter matrix has the same column ambiguity, so the frequency estimates finally obtained are paired. According to simulation results, the frequency estimation performance of the algorithm proposed in the present invention is close to that of the traditional parallel factor method, and is superior to some commonly used modern spectral estimation algorithms, and its complexity is lower than that of the traditional parallel factor method.

[0079] (1) Signal model

[0080] The power system frequency signal (voltage or current) containing noise, power frequency, harmonics and interharmonic components at the receiving end can be expressed as

[0081]

[0082] Where: M is the number of sinusoidal components including fundamental wave, harmonics and interharmonics; α m is the amplitude of the mth sinusoidal component; ω m is the angular frequency of the mth sinusoidal component; is the phase of the mth sinusoidal component; e(t) is the noise signal.

[0083] Using Euler's formula, it can be transformed into

[0084]

[0085] Where: Am is the amplitude of the harmonic signal, which is a complex constant; ω m is the normalized frequency to be estimated; is the normalized initial phase, uniformly distributed in the interval (-π,π); w is an uncorrelated complex Gaussian white noise with zero mean. If the noise component is not considered, it can be simplified to

[0086]

[0087] In the formula

[0088] In order to estimate the frequency of the received signal, a multi-delay output signal model of the receiving end is constructed. The P (P>2) discrete sinusoidal signals at the receiving end can be expressed as

[0089]

[0090] Where: {DP,p=1,2,…,p-1} is the multi-delay output; {wp(n),p=0,1,…,P-1} is the uncorrelated complex Gaussian white noise with zero mean.

[0091] According to the received signal, the following vector group is constructed:

[0092]

[0093] Where: Q is the vector length; [·]T represents the matrix transpose. Without considering the influence of noise, the signal in the formula can be expressed as

[0094]

[0095] Where:

[0096]

[0097]

[0098]

[0099] A is a frequency estimation matrix with Vandermonde characteristics. Define the delay matrix H and the source matrix S:

[0100]

[0101]

[0102] Construct a set of matrices:

[0103]

[0104] Then Z can be expressed as

[0105] Z p = ADiag p+1 (H)S,p = 0, 1, …, P - 1

[0106] where: Diagi(·) is a diagonal matrix constructed using the i-th row of the matrix; Z P represents the p-th slice in the time-delay direction.

[0107] Considering the influence of noise, the observation model becomes

[0108]

[0109] where WP is the received noise, which is a Q×L matrix.

[0110] By rearranging the above equation, a trilinear model can be obtained:

[0111]

[0112] p = 1, 2, …, P; l = 1, 2, …, L

[0113] where: a q,m is the (q, m) element of matrix A, h p,m is the (p, m) element of matrix H, s m,l is the (m, l) element of matrix S. Utilizing the symmetry of the trilinear model, the observation models in other directions are obtained by rearrangement:

[0114]

[0115]

[0116] (2) Model Compression

[0117] To reduce the complexity of the parallel factor decomposition algorithm, the three-dimensional matrix of the received data is compressed into a smaller matrix where Q’ < Q, P’ < P, L’ < L. The compressed matrix is defined as The above compression matrix needs to satisfy the RIP property and can be obtained through Tucker3 decomposition or a partial random matrix such as a Gaussian random matrix, a Bernoulli random matrix or a local Fourier matrix.

[0118] Then the compressed three-dimensional matrix can be expressed as

[0119] X' p' =A'Diag p'+1 (H')S',p'=0,1,…,P'-1

[0120] where A'=U H A, H'=V H H, S' = SW.

[0121] (3) Model decomposition

[0122] according to

[0123]

[0124] Get the least squares fit

[0125]

[0126] From the above formula, the least squares update of A' is

[0127]

[0128] in and are the estimated values ​​of S' and H' obtained in the previous iteration respectively.

[0129] Depend on

[0130]

[0131] The least squares fit is obtained as

[0132]

[0133] From the above formula, the least squares update of H' is

[0134]

[0135] Depend on

[0136]

[0137] The least squares fit is obtained as

[0138]

[0139] From the above formula, the least squares update of S' is

[0140]

[0141] The decomposition process of the TALS algorithm for the trilinear model has been given above, and the identifiability of the algorithm is given below. q=1,2,…,Q, where A is a Vandermont matrix if

[0142] k sT +min(Q+k H ,2M)≥2M+2

[0143] Then A, H and S are unique with respect to column transformation and scale transformation.

[0144] Consider the first row of the matrix A and the matrix H are all vectors of 1, and the normalized by trilinear decomposition Normalization and normalization as follows:

[0145]

[0146] Where П is the permutation matrix.

[0147] The signal in this article is an independent source Independent delay signal k H =min(P,M), and rewrite the above formula as

[0148] min(M,L)+min(Q+min(P,M),2M)≥2M+2

[0149] When L≥M, P≥M, M≥2, the identifiability is met, and M≤min(L,P). Obviously, the above parameters satisfy L≥M, P≥M, and the number of identifiable harmonics is min(L,P).

[0150] (4) Frequency estimation

[0151] Define the mth column of matrix A' as a'(ω m ), then Where a(ω m ) is the mth column of matrix A, is the scale factor; definition as a sampling grid of all frequencies with R>>M, used to construct a complete dictionary for frequency estimation if Rewrite the above formula as where e p is a P×1 sparse matrix. Using the L0 norm constraint, we can get ω m The estimate is

[0152]

[0153] (5) Algorithm complexity analysis and Cramer-Rao bound (CRB)

[0154] a) Complexity Analysis

[0155] For the proposed frequency estimation algorithm based on compressed sensing parallel factorization, the complexity is O[Q'L'P'M+3M 3 +M 2 (Q'L'+Q'P'+L'P'+L'+Q'+P')], while for the traditional parallel factor algorithm, the complexity is O[PQLM+3M 3 +M 2 (QL+QP+LP+P+Q+L)]. In order to intuitively compare the computational complexity of the two algorithms, the above model is substituted for comparison, considering Q=12, P=4, M=3, Q'=P'=L'=0.5. Figure 4 The computational complexity comparison between the traditional parallel factor algorithm and the algorithm proposed in this invention under different sampling numbers L is given. Figure 4 It can be seen from the figure that the computational complexity of the proposed algorithm is greatly reduced due to the compression of the received signal.

[0156] b) Cramer-Rao bound

[0157] definition

[0158] B=H⊙A

[0159] For the frequency estimation of harmonic signals, the Cramer-Rao bound is

[0160]

[0161] ⊙ denotes the Khatri-Rao product, Denotes the Hadamard product.

[0162] in, b m is the mth column of matrix B, in 1 3×3 is a 3×3 matrix with all elements set to 1.

[0163] (6) Simulation results

[0164] We assume that the harmonic signal received by the sensor is

[0165] x(t)=0.2cos(2π·45t)+cos(2π·50t)+0.2cos(2π·75t)+e(t)

[0166] This signal contains three frequency components: 50Hz power frequency, 45 and 75Hz interharmonics, and e(t) is Gaussian white noise.

[0167] L represents the number of signal samples, P is the number of delay levels, and Q is the number of sensors. In the simulation, we use the root mean square error (RMSE) of signal frequency estimation to evaluate the parameter estimation performance of the proposed algorithm, which is defined as

[0168]

[0169] in, is ω in the i-th Monte Carlo simulation m I is the total number of simulations, and we take I=100 in the following simulation.

[0170] Simulation 1: Figure 4 This is the distribution diagram of the frequency estimation results of the algorithm proposed in the invention when the signal-to-noise ratio is 20dB. The size of the parallel factor model in the algorithm is defined as Q×P×L. From the figure, it can be seen that the algorithm can still effectively identify fixed frequencies at a lower signal-to-noise ratio.

[0171] Simulation 2: Figure 5 The frequency estimation performance comparison chart of the proposed algorithm, the traditional parallel factor algorithm, and the ESPRIT algorithm. It can be seen from the figure that the frequency estimation performance of the proposed algorithm is close to the traditional parallel factor algorithm, and is better than the ESPRIT algorithm.

[0172] Simulation 3: Figure 6 The frequency estimation performance of the proposed algorithm under different sampling numbers L. Set Q = 12, P = 4, sampling number L = 50, L = 100, L = 150, L = 200. It can be seen from the figure that the frequency estimation performance of the proposed algorithm will improve as L increases.

[0173] Simulation 4: Figure 7 The frequency estimation performance of the proposed algorithm under different sensor numbers Q. Set L = 200, P = 4, and the number of sensors Q = 8, Q = 12, Q = 16, and Q = 20. It can be seen from the figure that the frequency estimation performance of the proposed algorithm will improve as Q increases.

[0174] Simulation 5: Figure 8 The performance of the proposed algorithm under the condition of dynamic frequency tracking. The power system simulation signal is

[0175] x(t)=0.2cos(2π·f1t)+cos(2π·f2t)+0.2cos(2π·f3t)+e(t)

[0176]

[0177] The figure shows the tracking performance of dynamic frequency. It can be seen that the algorithm can track dynamic frequency and realize dynamic frequency tracking and estimation.

[0178] The present invention proposes a method for estimating the frequency of harmonics and interharmonics in a power system based on compressed sensing parallel factor decomposition. First, the sensor receiving signal is constructed as a parallel factor trilinear model, and then the three-dimensional matrix of the received signal is compressed according to the compressed sensing theory. Finally, the compressed model is iteratively decomposed using a trilinear alternating least squares algorithm to obtain an estimate of the parameter matrix, and the frequency estimate of the signal is obtained from the obtained parameter matrix.

[0179] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing, characterized in that: The steps include: 1) Construct a parallel factor trilinear model based on the sensor's received signal model; 2) constructing a compression matrix to compress the constructed parallel factor trilinear model; 3) Using the trilinear alternating least squares algorithm, the compressed model is iteratively decomposed until convergence to obtain an estimate of the compressed parameter matrix; 4) Construct the sparse recovery problem in harmonic and interharmonic frequency estimation, and use the orthogonal matching pursuit method to obtain the harmonic and interharmonic frequency estimates; The specific steps of step (1) are: Define the signal model Z at the sensor receiving end p for: Z p =ADiag p+1 (H)S,p=0,1,L,P-1; Where: A is the frequency estimation matrix, H is the delay matrix, S is the source matrix, Diag p+1 (H) is a diagonal matrix constructed using the p+1th row of the delay matrix H, p is the signal model Z P The pth slice in the delay direction, P is the signal model Z p The number of delays; Considering the noise W P The sensor receiving end signal model Z p becomes Using the symmetry of the trilinear model, the signal models in the other two directions are obtained by rearrangement: Among them: Diag l (S T ) is to use the source matrix transpose S T A diagonal matrix constructed from the lth row of T is the transpose of the frequency estimation matrix A, L is the number of samples; Diag q (A) is a diagonal matrix constructed using the qth row of matrix A, H T is the transpose of the delay matrix H, Q is the number of sensors; The specific steps of step (2) are: Define three compression matrices: where Q’ < Q, P’ < P, L’ < L, and all three compression matrices satisfy the RIP property; the signal model at the sensor receiving end is compressed through the Tucker decomposition of the third-order tensor to obtain: Z’ p’ =A’Diag p’+1 (H’)S’,p’=0,1,L,P’-1; Where the parameter matrix A' = U H A, H'=V H H, S' = SW, similarly, we get the compressed signal Two other rearrangement matrices: X’ l’ =H’Diag l’ (S’ T )A’ T ,l’=1,2,L,L’, Y’ q’ =S’ T Diag q’ (A’)H’ T ,q’=1,2,L,Q’; The specific steps of step (4) are: Define the mth column of the parameter matrix A' as a'(ω m ), then: in Q' <Q,U H represents the conjugate transpose of U, a(ω m ) is the mth column of the frequency estimation matrix A, is the scale factor; definition As the sampling grid of all harmonic and interharmonic frequencies and R>>M, M represents the number of frequencies contained in the signal, it is used to construct a complete dictionary Ω for harmonic and interharmonic frequency estimation: in represents the rth column of the complete dictionary Ω, 1 <r<R; if a'(ω m ) is rewritten as: where e r is a P×1 sparse matrix, where P is the signal model Z p The number of delays; The frequency estimates of harmonics and interharmonics obtained by using the L0 norm constraint are: where || || F represents the Frobenius norm.

2. The harmonic and interharmonic frequency estimation method based on parallel factor compressed sensing according to claim 1, characterized in that: The specific steps of step (3) are: According to the trilinear alternating least squares algorithm, the least squares of the parameter matrices A', H', S' are updated as follows: in It is a compressed signal slice containing noise; the final estimated value of the parameter matrix A', H', S' is obtained by loop iteration.

Images