A dynamic window length-based synchronous phasor measurement method and storage medium

Abstract

The application relates to a dynamic window length-based synchronous phasor measurement method and a storage medium, the measurement method comprising the following steps: obtaining PMU sampling data, extracting an original sampling sequence, denoising and filtering the original sampling sequence, obtaining a pretreated signal and the component number of the pretreated signal; determining a fundamental frequency based on the pretreated signal and the component number; determining a window dynamic length according to the fundamental frequency; and realizing synchronous phasor calculation by using a discrete Fourier transform method based on the window dynamic length. Compared with the prior art, the application has the advantages of improving measurement accuracy and the like, and provides accurate data and information for power grid fault diagnosis, voltage stability monitoring, low-frequency oscillation analysis, transient stability analysis and control and the like.

Description

Technical Field

[0001] This invention relates to the field of power systems, specifically to a synchronous phasor measurement method, and more particularly to a synchronous phasor measurement method and storage medium based on dynamic window length. Background Technology

[0002] The integration of large-scale distributed power sources and electric vehicles into the power grid significantly increases its complexity and insecurity. To meet the requirements for real-time monitoring, fault diagnosis, and rapid control of the power grid, it is necessary to measure the grid's phasors. A phasor measurement unit (PMU) is used for measuring the grid's phasors. The implementation of phasor measurement by the PMU is becoming increasingly important for network monitoring of the power grid system. Measurement-based phasor measurement algorithms have high accuracy requirements, which means that the algorithms need to be able to withstand strong interference.

[0003] Currently widely used synchronous phasor measurement methods include zero-crossing detection, Kalman filtering, and the traditional Discrete Fourier Transform (DFT). Zero-crossing detection is simple in principle; after data is generated, a zero-crossing circuit detects the signal's zero-crossing moment, and then compares it with a reference moment to calculate the phase angle. However, zero-crossing detection has poor real-time performance and is easily affected by harmonics and noise; moreover, its measurement accuracy is very poor when the system frequency fluctuates. The core of Kalman filtering lies in establishing a state-space description of the signal and performing iterative calculations on individual sampled data. Its advantages include no need for a data window and good real-time performance. However, its disadvantage is that ordinary Kalman filtering cannot quickly track sudden changes in system parameters after the system reaches a steady state. Traditional DFT, when the system experiences disturbances and enters a dynamic state, will experience asynchronous sampling, introducing errors in amplitude and phase calculations. Furthermore, the error is proportional to the degree of asynchronous phenomenon, resulting in spectral leakage and the picket fence effect, thus reducing measurement accuracy. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a synchronous phasor measurement method and storage medium with high measurement accuracy based on dynamic window length.

[0005] The objective of this invention can be achieved through the following technical solutions:

[0006] A method for synchronous phasor measurement based on dynamic window length includes the following steps:

[0007] Acquire PMU sampling data, extract the original sampling sequence, denoise and filter the original sampling sequence to obtain the preprocessed signal and the number of components of the preprocessed signal;

[0008] The fundamental frequency is determined based on the preprocessed signal and the number of components.

[0009] The dynamic length of the window is determined based on the fundamental frequency, and the synchronous phasor is calculated using the discrete Fourier transform method based on this dynamic length of the window.

[0010] The original sampled sequence is denoised using singular value decomposition to obtain a denoised signal, and the denoised signal is then subjected to local weighted regression smoothing filtering to obtain the preprocessed signal.

[0011] Furthermore, the denoising of the original sampled sequence using singular value decomposition specifically includes:

[0012] Based on the original sampling sequence with a finite length N, construct an L×K Hankel matrix, where L is a pre-selected value and K = N - L + 1;

[0013] Perform singular value decomposition on matrix H to obtain the original singular value sequence. Then, calculate the average of the (2i-1)th and 2ith singular values ​​in sequence, i = 1, 2, ..., to obtain a new singular value sequence.

[0014] Calculate the curvature value of each new singular value in the new singular value sequence, where the largest curvature C is... k The inflection point represents the transition point between useful signal and noise. Its subscript k represents the position of the inflection point. The number of components M is determined based on the position of the inflection point, M = k-1.

[0015] Let r = 2M, extract the first r singular values ​​and their corresponding left and right singular vectors from the original singular value sequence, and calculate the denoising matrix H. de ;

[0016] matrix H de The denoised signal is obtained by performing anti-angle averaging.

[0017] Furthermore, the range of L is N / 3≤L≤N / 2.

[0018] Furthermore, the local weighted regression smoothing filter applied to the denoised signal specifically involves:

[0019] By selecting a local waveform and its corresponding position parameters in the denoised signal through a filtering window, selecting a filtered point position, solving the LWRS problem, obtaining the local waveform within the fitting filtering window, and then obtaining the filtered value at the filtered point position.

[0020] Move the filtering window and repeat the above filtering process until the entire signal sequence has been filtered, thus obtaining the preprocessed signal.

[0021] Furthermore, the fundamental frequency is determined using variational mode decomposition or empirical mode decomposition.

[0022] Furthermore, when using the variational mode decomposition method, the preprocessed signal is decomposed according to the number of components, and the final Fourier transform sequence of the preprocessed signal is obtained through iterative solution. The Fourier transform sequence is then subjected to inverse Fourier transform to obtain the fundamental signal and the corresponding fundamental frequency.

[0023] Furthermore, during the iterative solution process, the iteration stops when the Fourier transform sequence between two iterations satisfies the following formula or when the maximum number of iterations is reached:

[0024]

[0025] In the formula, ε is the given convergence tolerance. Let represent the Fourier transform sequence of the I-th iteration, and let the subscript c represent the c-th component.

[0026] Furthermore, the original sampling sequence is denoised using wavelet transform to obtain a denoised signal, and the denoised signal is then subjected to local weighted regression smoothing filtering to obtain the preprocessed signal.

[0027] The present invention also provides a computer-readable storage medium including one or more programs executable by one or more processors of an electronic device, the one or more programs including instructions for performing the synchronous phasor measurement method based on the described dynamic window length.

[0028] Compared with the prior art, the present invention has the following beneficial effects:

[0029] 1. This invention is based on PMU sampling data. By denoising the data, the impact of high-intensity background noise on the accurate calculation of synchronous phasors is reduced. At the same time, by determining the fundamental signal frequency, the window length of DFT calculation is dynamically changed to reduce the impact of spectral leakage and picket fence effect, thereby improving the calculation accuracy of synchronous phasors under asynchronous sampling.

[0030] 2. This invention uses a combined denoising filtering method of SVD and LWRS, which can improve the accuracy of subsequent measurements and effectively determine the number of components.

[0031] 3. The present invention adopts a frequency determination method based on VMD, which produces clear and reliable signals, providing a foundation for improving the accuracy of synchronous phasor measurement.

[0032] 4. This invention can provide accurate data and information for power grid fault diagnosis, voltage stability monitoring, low-frequency oscillation analysis, transient stability analysis and control. Attached Figure Description

[0033] Figure 1 This is a schematic diagram of the process of the present invention. Detailed Implementation

[0034] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. These embodiments are implemented based on the technical solution of the present invention, providing detailed implementation methods and specific operating procedures. However, the scope of protection of the present invention is not limited to the following embodiments.

[0035] This embodiment provides a synchronous phasor measurement method based on dynamic window length, such as... Figure 1 As shown, it includes the following steps:

[0036] Step S1: Acquire PMU sampling data, extract the original sampling sequence, denoise and filter the original sampling sequence to obtain the preprocessed signal and the number of components of the preprocessed signal.

[0037] In this embodiment, this step uses singular value decomposition to denoise the original sampling sequence to obtain a denoised signal, and then performs local weighted regression smoothing filtering on the denoised signal to obtain the preprocessed signal. Specifically, this step includes:

[0038] S101. The original sampled sequence is denoised using singular value decomposition.

[0039] Let x[n] = [x1 x2 … x N ] is an original sampling sequence of finite length N, which is constructed as an L×K Hankel matrix H.

[0040]

[0041] The parameter L affects the noise reduction effect of the subsequent Singular Value Decomposition (SVD). Generally, it is more appropriate to take N / 3≤L≤N / 2, and the value of K is N-L+1.

[0042] Performing singular value decomposition on matrix H yields:

[0043]

[0044] In the formula, d is the number of non-zero singular values ​​of matrix H; σ i Let u be the i-th singular value of matrix H arranged in descending order; i and v i These correspond to σ i The left and right singular vectors.

[0045] According to the theory of singular value decomposition and the matrix best approximation theorem in the sense of Frobeious norm, the useful signal, i.e., the fundamental and harmonic signals, is mainly reflected by the first r larger singular values, while noise is represented by the remaining smaller singular values. Therefore, the denoising effect of SVD depends on whether the choice of r is appropriate.

[0046] The average of the (2i-1)th and 2ith (i = 1, 2, ...) singular values ​​is used to obtain a new sequence of singular values. The curvature value corresponding to each new singular value is calculated using equation (3).

[0047]

[0048] In the formula, C m For λ m The corresponding curvature; λ m K is the m-th new singular value; m For λ m The smaller of the absolute values ​​of the forward and backward differences.

[0049] Maximum curvature C k This represents the inflection point between the useful signal and noise, and its subscript k represents the location of the inflection point. Therefore, the number of signal components can be determined as follows:

[0050] M = k-1 (4)

[0051] Therefore, the value of r can be determined as:

[0052] r = 2M (5)

[0053] The SVD denoising matrix H is obtained by calculating the first r original singular values ​​and their corresponding left and right singular vectors using equation (2). de .

[0054] matrix H de By performing anti-angle averaging, we can obtain the SVD-denoised signal y[n] = [y1 y2…y N ].

[0055] In other implementations, other denoising methods besides SVD, such as WT, can be used to achieve denoising.

[0056] S102. Perform local weighted regression smoothing filtering on the denoised signal. Select a local waveform and its corresponding position parameters in the denoised signal through a filtering window, select a filtering point position, solve the LWRS problem, obtain the local waveform within the fitting filtering window, and then obtain the filtering value at the filtering point position; move the filtering window and repeat the above filtering process until the entire signal sequence is filtered, that is, obtain the preprocessed signal.

[0057] Under strong noise, the noise has already contaminated the useful signal, and the first r singular values ​​are contributed by both the useful signal and the noise. Therefore, further SVD denoising is needed for the signal y[n]=[y1 y2…y N Locally weighted regression smoothing (LWRS) filtering is performed to eliminate the influence of residual noise.

[0058] Let y p [n] = [y p1 y p2 …y pn [ ] is a local waveform of length n of the SVD denoised signal y[n], and its corresponding position parameter is l. p [n] = [l p1 l p2 …l pn ], l ps The location of the filtered point is set to the midpoint of the local waveform in this embodiment. Let the parameters of the t-th order polynomial be {α}. j ,j=0,1,...,t},w i For (l) pi y pi If the fitting weight parameters for the points are given, then the LWRS problem can be expressed as:

[0059]

[0060] w i Values ​​of the cubic kernel function:

[0061]

[0062] In the formula, d(l ps ) is (l ps y ps The farthest distance from α to other points in the filter window. The weighted least squares solution for α is:

[0063]

[0064] In the formula, L p Let Y be a position parameter matrix of size n×t. p Let W be a local waveform value matrix of size n×1, and let W be an n×n weighted diagonal matrix.

[0065] After obtaining the value of α, the local waveform within the filtering window can be fitted to obtain l. ps The filtered value z of the point ps Move the filtering window and repeat the above filtering process until the entire signal sequence has been filtered, thus obtaining the signal z[n] = [z1 z2 ... z...]. N ].

[0066] Step S2: Determine the fundamental frequency based on the preprocessed signal and the number of components. When using the variational mode decomposition method, the preprocessed signal is decomposed according to the number of components. Through iterative solution, the final Fourier transform sequence of the preprocessed signal is obtained. The inverse Fourier transform of the Fourier transform sequence is then performed to obtain the fundamental signal and its corresponding fundamental frequency.

[0067] Signal z[n]=[z1 z2…z N The window contains fundamental and harmonic components. To determine the dynamic length of the window, the frequency of the fundamental wave must be determined. Variational Mode Decomposition (VMD) is essentially the process of solving a variational problem. Its algorithm mainly includes constructing and solving the variational problem. During decomposition, the number of components needs to be set; in SVD denoising, we determined the number of components as C = M = r / 2.

[0068] First, initialize the parameters. and iteration number I, express The Fourier transform of u1, u2, ..., u2, u3, where the superscript 1 indicates the initial values ​​of the first iteration. c This represents the C components decomposed from the original structure; express The Fourier transform of ω1, ω2, ..., ω2, where the superscript 1 indicates the initial values ​​of the first iteration. c Indicates the center frequency of the corresponding component; Represents {λ 1 The Fourier transform of}, where the superscript 1 indicates the initial value of the first iteration, and λ represents the Lagrange multiplication operator.

[0069] Then update based on (9) and (10) and

[0070]

[0071] In the formula, α is the set penalty factor.

[0072]

[0073] Then update

[0074]

[0075] In the formula, τ is the set noise tolerance.

[0076] When the result between two iterations satisfies equation (12) or I reaches the maximum number of iterations I. max When the time comes, the VMD algorithm stops.

[0077]

[0078] In the formula, ε is the given convergence tolerance.

[0079] For the final result Performing the inverse Fourier transform, the actual part is {u c (t)}, where the one closest to 50Hz is the fundamental signal, and its frequency f0 is what we need.

[0080] In other implementations, signal decomposition methods other than VMD, such as EMD, can be used for signal decomposition.

[0081] Step S3: Determine the dynamic length of the window based on the fundamental frequency, and use the discrete Fourier transform method to calculate the synchronous phasor based on the dynamic length of the window.

[0082] After determining the fundamental frequency f0, the actual number of sampling points per cycle can be calculated:

[0083]

[0084] In the formula, f s Where N is the sampling frequency. i Represents N a The integer part of N; f Represents N a The fractional part.

[0085] According to the principle of DFT, the DFT at time t0 can be written as:

[0086]

[0087] Let P = N i -1+N f In the formula and It can be calculated using the following formula:

[0088]

[0089] In the formula, T s The sampling period.

[0090] Since P is not an integer, z(t0+PT) s ) is not in the sampling sequence, therefore N can be used f ·z(t0+N i T s If we replace ) with , then in equation (15) It can be rewritten as:

[0091]

[0092] According to equations (15) and (16), the real and imaginary parts of the DFT at time t0 are as follows:

[0093]

[0094]

[0095]

[0096] By X R and X I The magnitude and phase angle of the phasor at time t0 can then be calculated:

[0097]

[0098]

[0099] If the above methods are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0100] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A method for synchronous phasor measurement based on dynamic window length, characterized in that, Includes the following steps: Acquire PMU sampling data, extract the original sampling sequence, denoise and filter the original sampling sequence to obtain the preprocessed signal and the number of components of the preprocessed signal; The fundamental frequency is determined based on the preprocessed signal and the number of components. The dynamic length of the window is determined based on the fundamental frequency, and the synchronous phasor is calculated using the discrete Fourier transform method based on this dynamic length of the window. The fundamental frequency is determined by variational mode decomposition, specifically by decomposing the preprocessed signal according to the number of components, obtaining the final Fourier transform sequence of the preprocessed signal through iterative solution, and performing an inverse Fourier transform on the Fourier transform sequence to obtain the fundamental signal and its corresponding fundamental frequency.

2. The synchronous phasor measurement method based on dynamic window length according to claim 1, characterized in that, The original sampled sequence is denoised using singular value decomposition to obtain a denoised signal, and the denoised signal is then subjected to local weighted regression smoothing filtering to obtain the preprocessed signal.

3. The synchronous phasor measurement method based on dynamic window length according to claim 2, characterized in that, The denoising of the original sampled sequence using singular value decomposition specifically includes: Based on having a finite length N From the original sampling sequence, construct a L × K Hankel matrix, L For a pre-selected value, K = N - L +1; For matrix H Perform singular value decomposition to obtain the original singular value sequence, and then sequentially extract the second singular value from it. i -1 and the 2nd i The average value is calculated from the singular values. i =1,2,…, to obtain a new sequence of singular values; Calculate the curvature value of each new singular value in the new singular value sequence, where the largest curvature is... It represents the inflection point between useful signal and noise, and its subscript... k This represents the location of the inflection point, and the number of components is determined based on the location of this inflection point. , ; make Extract the first few idiomatic values ​​from the original singular value sequence. r The denoising matrix is ​​obtained by calculating the sum of the singular values ​​and their corresponding left and right singular vectors. H de ; matrix H de The denoised signal is obtained by performing anti-angle averaging.

4. The synchronous phasor measurement method based on dynamic window length according to claim 3, characterized in that, L The range of values ​​is N / 3≤ L ≤ N / 2.

5. The synchronous phasor measurement method based on dynamic window length according to claim 2, characterized in that, The local weighted regression smoothing filter for the denoised signal is specifically performed as follows: By selecting a local waveform and its corresponding position parameters in the denoised signal through a filtering window, selecting a filtered point position, solving the LWRS problem, obtaining the local waveform within the fitting filtering window, and then obtaining the filtered value at the filtered point position. Move the filtering window and repeat the above filtering process until the entire signal sequence has been filtered, thus obtaining the preprocessed signal.

6. The synchronous phasor measurement method based on dynamic window length according to claim 1, characterized in that, During the iterative solution process, the iteration stops when the Fourier transform sequence between two iterations satisfies the following formula or when the maximum number of iterations is reached: In the formula, For a given convergence tolerance, Indicates the first I The Fourier transform sequence of the next iteration, index c Indicates the first c Each component.

7. The synchronous phasor measurement method based on dynamic window length according to claim 1, characterized in that, The original sampling sequence is denoised using wavelet transform to obtain a denoised signal, and the denoised signal is then subjected to local weighted regression smoothing filtering to obtain the preprocessed signal.

8. A computer-readable storage medium, characterized in that, Includes one or more programs executable by one or more processors of an electronic device, the one or more programs including instructions for performing the synchronous phasor measurement method based on dynamic window length as described in any one of claims 1-7.

Images