CN114781196A - Harmonic detection method based on sparse acquisition model - Google Patents

Harmonic detection method based on sparse acquisition model Download PDF

Info

Publication number
CN114781196A
CN114781196A CN202210698874.3A CN202210698874A CN114781196A CN 114781196 A CN114781196 A CN 114781196A CN 202210698874 A CN202210698874 A CN 202210698874A CN 114781196 A CN114781196 A CN 114781196A
Authority
CN
China
Prior art keywords
harmonic
phasor
frequency
matrix
model
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210698874.3A
Other languages
Chinese (zh)
Inventor
迟傲冰
郭昱甫
郭诗瑶
潘禹飞
曾成碧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan University
Original Assignee
Sichuan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sichuan University filed Critical Sichuan University
Priority to CN202210698874.3A priority Critical patent/CN114781196A/en
Publication of CN114781196A publication Critical patent/CN114781196A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Complex Calculations (AREA)

Abstract

The invention discloses a harmonic detection method based on a sparse acquisition model, which comprises the following steps of: establishing a broadband dynamic signal model, and calculating an attenuation direct current component; estimating harmonic phasor in the broadband dynamic signal: establishing a harmonic phasor estimation model, reconstructing the harmonic phasor estimation model, and solving the harmonic estimation model; simulation analysis: and establishing a test scene, and sequentially establishing a basic performance test scene, a frequency slope test scene, a step change test scene and an interference test scene. The method has higher precision and speed for estimating the broadband dynamic phasor, keeps certain superiority in a comparison algorithm, and simultaneously verifies the dynamic performance and the anti-interference capability of the method.

Description

Harmonic detection method based on sparse acquisition model
Technical Field
The invention relates to the field of electronic information, in particular to a harmonic detection method based on a sparse acquisition model.
Background
With the rapid development of smart grids and new energy power generation, a large number of nonlinear devices are applied to power systems, resulting in the increase of inter-harmonic and DDC components in the grids. On one hand, the harmonic waves and the inter-harmonic waves have frequency spectrum interference on the concerned frequency, and the accurate metering of the electric energy is influenced. On the other hand, aiming at the attenuated direct current component in the electric signal when the power system fails, the traditional algorithm cannot filter the attenuated direct current component, so that the error of the calculation result is large. Therefore, an effective phasor estimation algorithm needs to be found to accurately measure the broadband dynamic phasor.
At present, phasor measurement methods are mainly classified into two types: DFT algorithms (discrete fourier transform) and non-DFT algorithms. The DFT algorithm is to sample the frequency spectrum of a finite time domain discrete signal at equal intervals and discretize a frequency domain function. When the traditional DFT algorithm is adopted to measure phasor, the problems of information redundancy, mutual harmonic interference, information leakage and the like can occur when the actual frequency deviates, and a larger measurement error is caused. Some studies propose to introduce a zero-crossing detection method, a wavelet transform method, an instantaneous value method and the like to optimize the DFT algorithm, but the methods still have large errors at non-nominal frequencies. The related document describes the use of improved second harmonic filtering techniques for monophasic phasor measurements at non-nominal frequencies. The method can integrate uniform sampling and fixed window length phasor measurement without changing the structure of the traditional phasor processing method, and has good performance. However, the phasor estimation algorithm based on DFT under the static model has the problem of poor dynamic performance.
For the bad point that the cracked DTF model is only limited to static signal analysis, some research proposals develop dynamic phasor estimation to different mathematical frameworks to obtain a new generation phasor measurement method based on the model. The method applies a specific mathematical model to define the amplitude and phase changes of the signal in a stable sampling time interval, and gets rid of the limitation of the period in Fourier analysis.
Currently, dynamic harmonic analysis methods, such as taylor-fourier transform (TFT), can achieve phasor estimation in the case of power oscillations, but are susceptible to interference from inter-and higher harmonics. The invention widens the frequency range of the TFT and overcomes the limitation of DFT periodicity. However, the method has many and complex system parameters, and has a problem of large error after multiple fitting, and the measurement error increases with the increase of the order. In order to greatly improve the estimation accuracy, the related literature proposes a discrete fourier transform interpolation method with high resolution. The invention has the capability of recovering the stage information and achieves the noise suppression by blocking the iteration times. However, the DFT sequence in this method can only be resampled at the nominal rate, which undoubtedly adds to the uncertainty of the effect.
Based on the methods, matrix transformation coefficients are obtained through high-order matrix inversion operation. Despite the many effective error reduction algorithms, the effects of higher order inter-harmonics and noise superposition cannot be addressed. Thus, consider incorporating a machine learning algorithm to reduce phasor solution computational complexity. The related literature effectively tracks harmonics for modeling estimation by identifying the most relevant components in the signal, thereby limiting the effects of inter-harmonics. However, this method is only suitable for low-band harmonic vectors, and the phasor estimation parameters cannot be directly obtained because the high-order harmonic signals have time-varying property and the solution of the high-order pseudo-inverse matrix has high complexity. And further considering the field of dynamic harmonic analysis, related documents perform iterative estimation on harmonic parameters by using an STWLS algorithm, and compensate the influence of the harmonic parameters on the estimation. The method can improve the frequency measurement precision, allocate more space for dynamic information, and still be difficult to resist the interference of constraining the inter-harmonic. The related literature proposes O-spline dynamic harmonic analysis, and the method reduces the calculation complexity of the DTTFT in a closed form and provides an optimal data compression algorithm for oscillation. However, under non-ideal conditions, inter-harmonic interference becomes more and more severe as the order increases due to the presence of large noise in the spatial step reconstruction process.
In order to solve the problem of solving the synchronous phasor estimation, a published related document proposes complex wavelet transform based on a Fast Fourier Transform (FFT) algorithm, which is used for segmenting an original waveform and analyzing signal reconstruction. However, in continuous multi-scale analysis, this may result in severe phase distortion and loss of timing information for further signal processing. There is also a published relevant document to obtain frequency, amplitude and phase estimates in harmonics by improving the Taylor weighted least squares algorithm. The method is not suitable for high-order harmonic measurement because the method is only suitable for measuring third-order continuous waves and cannot be used for high-precision harmonic analysis.
Disclosure of Invention
Based on the technical problem, the invention provides a harmonic detection method based on a sparse acquisition model, which comprises the steps of firstly carrying out parameterized modeling on an attenuation direct current component by utilizing a sinc interpolation function, obtaining an accurate DDC component through a least square method, then adopting an improved Taylor-Fourier model for broadband harmonic component estimation, combining a machine learning algorithm, accurately recovering a specific signal by using fewer data points, and further converting a harmonic phasor reconstruction problem into a sparse acquisition model solving problem. Under the regularization frame, the invention expresses the phasor reconstruction problem into a form of a minimization function, and the underdetermined inverse problem is effectively solved by utilizing alternate minimization operation under an SBI system, so that the harmonic phasor estimation value is efficiently and accurately obtained. To prove the rationality of the model, the error range of the estimated values is measured using a cross-entropy objective function. Simulation results show that the method has higher precision and speed for estimating the broadband dynamic phasor, certain superiority is kept in a comparison algorithm, and the dynamic performance and the anti-interference capability of the method are verified.
The technical scheme of the invention is as follows: a harmonic detection method based on a sparse acquisition model comprises the following steps:
a. estimating the attenuation direct current component in the broadband dynamic signal: establishing a broadband dynamic signal model, and calculating an attenuated direct current component;
b. estimating harmonic phasor in the broadband dynamic signal: establishing a harmonic phasor estimation model, reconstructing the harmonic phasor estimation model, and solving the harmonic estimation model;
c. a simulation analysis step: and establishing a test scene, and sequentially establishing a basic performance test scene, a frequency slope test scene, a step change test scene and an interference test scene.
In step a, the step of establishing a broadband dynamic signal model is as follows:
for power system waveforms containing a DC component and a fundamental component
Figure 996151DEST_PATH_IMAGE001
The expression is shown as (1):
Figure 241188DEST_PATH_IMAGE002
(1)
wherein, the first and the second end of the pipe are connected with each other,
Figure 489766DEST_PATH_IMAGE003
is a signal containing fundamental, harmonic and inter-harmonic components,
Figure 150555DEST_PATH_IMAGE004
is a time-varying amplitude of the signal,
Figure 43556DEST_PATH_IMAGE005
is a time-varying phase angle, the reference frequency is
Figure 131597DEST_PATH_IMAGE006
The term "refers to the highest order harmonic component,
Figure 601893DEST_PATH_IMAGE007
is to attenuate the signal of the direct current component,
Figure 66372DEST_PATH_IMAGE008
is the amplitude of the attenuated dc component,
Figure 63147DEST_PATH_IMAGE009
is a time constant;
signals containing fundamental, harmonic and inter-harmonic components
Figure 322090DEST_PATH_IMAGE003
The dynamic phasor model at time t is:
Figure 279682DEST_PATH_IMAGE010
(2)
substituting equation (2) into equation (1) to obtain broadband dynamic signal
Figure 282273DEST_PATH_IMAGE011
The dynamic expression of (2):
Figure 805658DEST_PATH_IMAGE012
(3)
wherein the content of the first and second substances,
Figure 311201DEST_PATH_IMAGE013
is a conjugate operator.
In step a, the step of "calculating the attenuated dc component" is as follows:
attenuating DC component
Figure 818406DEST_PATH_IMAGE014
Represented by the sum of the amplitude and phase time-varying cosine components:
Figure 296792DEST_PATH_IMAGE015
Figure 205842DEST_PATH_IMAGE016
(4)
wherein, the first and the second end of the pipe are connected with each other,
Figure 868905DEST_PATH_IMAGE017
，
Figure 863406DEST_PATH_IMAGE018
，
Figure 145482DEST_PATH_IMAGE019
respectively representing the frequency, amplitude and phase of the nth cosine component, T is the length of an observation interval, and Y is set as3；
Based on the frequency domain sampling theorem, the time domain signal parameterization modeling is carried out on the attenuation direct current component by utilizing a sinc interpolation function, and the dynamic phasor of the DDC component is recorded as
Figure 909039DEST_PATH_IMAGE020
Specifically, it is represented as:
Figure 759314DEST_PATH_IMAGE021
(5)
wherein, the first and the second end of the pipe are connected with each other,
Figure 975532DEST_PATH_IMAGE022
is the sampling frequency, represents the maximum model order, will be set to 2,
Figure 858037DEST_PATH_IMAGE023
is the nth cosine component in the formula (2) at the time
Figure 413784DEST_PATH_IMAGE024
The sample phasor of (c);
for attenuating DC component
Figure 621911DEST_PATH_IMAGE014
Sampling is carried out, and the discrete form of the fitted DDC component in the formula (1) is represented by (6);
Figure 387742DEST_PATH_IMAGE025
(6)
wherein the content of the first and second substances,
Figure 73938DEST_PATH_IMAGE026
representing phasors
Figure 749770DEST_PATH_IMAGE027
Is a column vector of the cosine signal,
Figure 128799DEST_PATH_IMAGE013
as conjugate operators, matrices
Figure 54030DEST_PATH_IMAGE028
Is an element of
Figure 356966DEST_PATH_IMAGE029
At a sample point in the sampling window(s),
Figure 684042DEST_PATH_IMAGE030
representation matrix
Figure 171655DEST_PATH_IMAGE031
，
Figure 849761DEST_PATH_IMAGE029
The definition of (2) is shown as (7):
Figure 205656DEST_PATH_IMAGE032
(7)
determining phasor estimation values by means of least squares
Figure 652818DEST_PATH_IMAGE033
The calculation process is shown in (7) in the matrix
Figure 45753DEST_PATH_IMAGE033
In (1),
Figure 211155DEST_PATH_IMAGE034
for the matrix elements:
Figure 777266DEST_PATH_IMAGE035
(8)
wherein H is a conjugate transposed symbol;
by constructing each cosine component using equation (8)
Figure 157563DEST_PATH_IMAGE029
Followed by
Figure 783716DEST_PATH_IMAGE029
Reconstruction matrix
Figure 374098DEST_PATH_IMAGE036
The amplitude and the phase of each cosine signal are obtained through (5), and then DDC components are estimated
Figure 9478DEST_PATH_IMAGE027
。
In step b, the step of establishing a harmonic phasor estimation model is as follows:
the dynamic phasors are defined using a taylor expansion model of order:
Figure 227970DEST_PATH_IMAGE037
Figure 290604DEST_PATH_IMAGE038
(9)
wherein, the first and the second end of the pipe are connected with each other,
Figure 102702DEST_PATH_IMAGE039
is t =0
Figure 276195DEST_PATH_IMAGE040
K is the taylor expansion order, and T is the observation interval length;
the sinusoidal signals containing harmonics and inter-harmonics take the form of discrete sequences
Figure 552455DEST_PATH_IMAGE041
The sequence length is N, and
Figure 596110DEST_PATH_IMAGE042
and expanding the dynamic phasor method to an actual sequence model through multi-frequency phasor analysis, wherein the Taylor Fourier coefficient expression of each phasor component of the discrete signal is as follows:
Figure 957821DEST_PATH_IMAGE043
(10)
where H is the sampling interval, observation interval length
Figure 872688DEST_PATH_IMAGE044
，
Figure 3455DEST_PATH_IMAGE045
Is the average harmonic and inter-harmonic phasors,
Figure 470208DEST_PATH_IMAGE046
is at phasor frequency
Figure 319215DEST_PATH_IMAGE047
The derivative of the order k of (c) is,
Figure 772193DEST_PATH_IMAGE047
is the fundamental frequency
Figure 23046DEST_PATH_IMAGE048
Integer multiples of (d) represent harmonic frequencies, and non-integer multiples represent inter-harmonic frequencies;
in a size of
Figure 332805DEST_PATH_IMAGE049
The sampling rate of (2) normalizes the frequency of the harmonic signal by a sampling length of
Figure 482158DEST_PATH_IMAGE050
Obtaining the frequency of the harmonic component
Figure 801143DEST_PATH_IMAGE051
The normalized frequency of the Taylor Fourier basis vector is
Figure 844186DEST_PATH_IMAGE052
，
Figure 590425DEST_PATH_IMAGE053
The frequency resolution corresponding to the N coefficient sets is
Figure 210762DEST_PATH_IMAGE054
And the taylor opening order of each i is K, then (10) is written in matrix form:
Figure 333439DEST_PATH_IMAGE055
(11)
wherein the content of the first and second substances,
Figure 230988DEST_PATH_IMAGE056
is a sample column phasor of size
Figure 148128DEST_PATH_IMAGE057
Is a Taylor Fourier-based phasor, e represents noise or measurement error, and x is a length of
Figure 272073DEST_PATH_IMAGE058
The column phasors of (1), describe the current harmonics
Figure 932862DEST_PATH_IMAGE059
When the utility model is used, the water is discharged,
Figure 684917DEST_PATH_IMAGE060
a set of (a);
and (3) phasor solving, namely converting the phasor solving into solving of a pseudo-inverse matrix by using a least square method, and calculating the pseudo-inverse matrix to obtain a corresponding solution when the Euclidean norm is minimum, wherein the solution is the optimal solution under the following constraint conditions:
Figure 38538DEST_PATH_IMAGE061
(12)
wherein, the first and the second end of the pipe are connected with each other,
Figure 633467DEST_PATH_IMAGE062
expressing Euclidean norm, and further obtaining a dynamic phasor coefficient estimation matrix by introducing a Lagrange operator and performing derivation
Figure 97947DEST_PATH_IMAGE063
：
Figure 704508DEST_PATH_IMAGE064
(13)
In step b, the step of reconstructing the harmonic phasor estimation model is as follows:
setting a correction factor
Figure 229031DEST_PATH_IMAGE065
Each group of frequency resolution reaches
Figure 61989DEST_PATH_IMAGE066
The length of the sampling sequence is
Figure 64580DEST_PATH_IMAGE067
The sampling length H is 5 fundamental wave periods, and the harmonic component normalized frequency is
Figure 791227DEST_PATH_IMAGE068
，
Figure 486651DEST_PATH_IMAGE069
And obtaining an improved Taylor Fourier coefficient expression:
Figure 790593DEST_PATH_IMAGE070
(14)
wherein the content of the first and second substances,
Figure 596875DEST_PATH_IMAGE071
，
Figure 443608DEST_PATH_IMAGE072
is expressed as a size of
Figure 44354DEST_PATH_IMAGE073
The matrix of the curvelet transform of (c),
Figure 114554DEST_PATH_IMAGE074
is that
Figure 396631DEST_PATH_IMAGE072
The k-th derivative of (a);
according to the expression (9), an approximate second-order dynamic phasor estimation value and a harmonic frequency are obtained
Figure 160187DEST_PATH_IMAGE075
And frequency rate of change
Figure 259730DEST_PATH_IMAGE076
:
Figure 679210DEST_PATH_IMAGE077
(15)
Figure 905923DEST_PATH_IMAGE078
(16)
Figure 727249DEST_PATH_IMAGE079
(17)
When correcting the factor
Figure 935376DEST_PATH_IMAGE080
When the sensing matrix A and the sparse matrix are in use
Figure 966786DEST_PATH_IMAGE081
Highly uncorrelated, (14) is converted into matrix form, namely:
Figure 590666DEST_PATH_IMAGE082
(18)
wherein, the first and the second end of the pipe are connected with each other,
Figure 328814DEST_PATH_IMAGE083
is a measurement matrix, x represents a length of
Figure 520893DEST_PATH_IMAGE084
For each data block ofIn the present sequence, the sequence of the method,
Figure 914965DEST_PATH_IMAGE081
is the superposition of a small number of phasors extracted after the original sample is segmented,
Figure 732748DEST_PATH_IMAGE085
is dimension of
Figure 263086DEST_PATH_IMAGE086
Each column of the matrix is a base phasor of a discrete fourier transform, e represents a noise signal phasor;
after the sparse sampling measurement model is constructed, matrix blocking is carried out on x by using an Euclidean search algorithm, m data points with optimal matching are found in a search range, and all the data points are combined into a matrix
Figure 813017DEST_PATH_IMAGE087
I is a coordinate of the coefficient matrix, curvelet transformation is carried out on the matrix, an insignificant curvelet coefficient is removed by adopting a curvelet threshold criterion, an effective denoising and compression algorithm is required for the extracted coefficient, curvelet transformation is carried out, an insignificant curvelet coefficient is removed by adopting a curvelet threshold criterion, and an effective denoising and compression algorithm is required for the extracted coefficient;
sparse matrix
Figure 304172DEST_PATH_IMAGE081
Only contains a non-zero significant term set, the phasor of the transformation coefficient column obtained after the significant term set is arranged according to the dictionary order is
Figure 332171DEST_PATH_IMAGE088
Of the transform coefficients
Figure 717016DEST_PATH_IMAGE088
Centered around the zero region and distributed in the form of a fine peak, the signal repeatability is characterized using the distribution characteristics of the transform coefficients, noted as:
Figure 765743DEST_PATH_IMAGE089
(19)
wherein the content of the first and second substances,
Figure 931145DEST_PATH_IMAGE090
represents the L1 norm;
each set of harmonic signal estimates may be inverted according to the transform coefficients:
Figure 497256DEST_PATH_IMAGE091
(20)
wherein, the first and the second end of the pipe are connected with each other,
Figure 736607DEST_PATH_IMAGE092
is that
Figure 628340DEST_PATH_IMAGE093
The inverse operator of (2).
In step b, the step of solving the harmonic estimation model is as follows:
solving a regularized linear least square cost function, and reconstructing a Taylor coefficient matrix x, namely:
Figure 102876DEST_PATH_IMAGE094
(21)
Figure 675940DEST_PATH_IMAGE095
(22)
wherein, the first and the second end of the pipe are connected with each other,
Figure 97694DEST_PATH_IMAGE096
the norm of L2 is shown,
Figure 222645DEST_PATH_IMAGE097
a regularization term is represented as a function of,
Figure 97060DEST_PATH_IMAGE098
is a coefficient of the regularization that,
Figure 473815DEST_PATH_IMAGE099
、
Figure 750076DEST_PATH_IMAGE100
representing local and non-local sparse terms respectively,
Figure 796660DEST_PATH_IMAGE101
is a regularization parameter that balances the two sparse terms;
(22) substituting (21) to obtain:
Figure 158371DEST_PATH_IMAGE102
(23)
wherein the content of the first and second substances,
Figure 807658DEST_PATH_IMAGE103
，
introducing two auxiliary phasors p and q, finally the following scheme is obtained:
Figure 204005DEST_PATH_IMAGE104
(24)
Figure 670758DEST_PATH_IMAGE105
(25)
Figure 519765DEST_PATH_IMAGE106
(26)
Figure 972744DEST_PATH_IMAGE107
(27)
Figure 223596DEST_PATH_IMAGE108
(28)
wherein, the first and the second end of the pipe are connected with each other,
Figure 533355DEST_PATH_IMAGE109
and
Figure 682708DEST_PATH_IMAGE110
the parameters are fixed value parameters, the function is to improve the stability of the algorithm value, and p and q are SBI algorithm auxiliary iteration phasors;
(24) solving the minimization of a strict convex quadratic function, setting the gradient of an objective function as 0, and obtaining a corresponding closed solution:
Figure 1693DEST_PATH_IMAGE111
(29)
wherein the content of the first and second substances,
Figure 44736DEST_PATH_IMAGE112
is of the size
Figure 790975DEST_PATH_IMAGE113
The identity matrix of (a);
the minimization problem of the function is solved by adopting the steepest descent method with the optimal step length, which is expressed as follows:
Figure 411312DEST_PATH_IMAGE114
，
Figure 533989DEST_PATH_IMAGE115
(30)
wherein, the first and the second end of the pipe are connected with each other,
Figure 431538DEST_PATH_IMAGE116
is the gradient of the objective function and,
Figure 348678DEST_PATH_IMAGE117
is the optimum step size for the particular application,
Figure 472623DEST_PATH_IMAGE118
the estimated value is
Figure 133412DEST_PATH_IMAGE119
Using the least mean square error theorem of (25)The problem becomes:
Figure 885467DEST_PATH_IMAGE120
(31)
similarly, the question of (26) is written as:
Figure 239088DEST_PATH_IMAGE121
(32)
here, the
Figure 834017DEST_PATH_IMAGE122
(omission of
Figure 767338DEST_PATH_IMAGE123
）。
Obtaining a set of closed solutions, which are
Figure 701796DEST_PATH_IMAGE124
Observed values seen as noise
Figure 164002DEST_PATH_IMAGE125
Of a certain type, use
Figure 918331DEST_PATH_IMAGE126
Representing the error phasor, the error of each element is respectively
Figure 996621DEST_PATH_IMAGE127
（
Figure 785585DEST_PATH_IMAGE128
) Then assume e is independently distributed over
Figure 153113DEST_PATH_IMAGE129
Mean 0 and variance
Figure 660317DEST_PATH_IMAGE130
From law of large numbers, for arbitrary
Figure 528916DEST_PATH_IMAGE131
All have:
Figure 172387DEST_PATH_IMAGE132
(33)
so as to obtain:
Figure 976395DEST_PATH_IMAGE133
(34)
the transformed error phasor is
Figure 705317DEST_PATH_IMAGE134
，
Figure 128339DEST_PATH_IMAGE135
Figure 891896DEST_PATH_IMAGE136
Representing the error phasor for each element, m being the number of data points for the best match, and the transformation does not change the variance of each group according to the orthogonal nature of the matrix, from which it follows that for each group
Figure 601226DEST_PATH_IMAGE137
Are all independently distributed in
Figure 817443DEST_PATH_IMAGE138
(zero mean and variance of
Figure 965528DEST_PATH_IMAGE139
) Most of them are theoretic and arbitrary
Figure 645908DEST_PATH_IMAGE140
All have:
Figure 588456DEST_PATH_IMAGE141
(35)
the method comprises the following steps:
Figure 229653DEST_PATH_IMAGE142
(36)
to obtain:
Figure 181429DEST_PATH_IMAGE143
(37)
combining the above and the (26) problem to obtain:
Figure 467048DEST_PATH_IMAGE144
(38)
unknown quantity
Figure 846076DEST_PATH_IMAGE145
Is separable in the above formula, and each component is simplified as follows:
Figure 771307DEST_PATH_IMAGE146
(39)
wherein the content of the first and second substances,
Figure 198877DEST_PATH_IMAGE147
refers to the modulus of the corresponding phasor;
under the vector contraction form, the estimated value of x is expressed as:
Figure 791533DEST_PATH_IMAGE148
(40)
after the estimation value of the Taylor coefficient matrix x is obtained, the method is carried out
Figure 403780DEST_PATH_IMAGE149
Harmonic sampling equation to harmonic signals
Figure 19569DEST_PATH_IMAGE150
In the framework of logistic regression, cross entropy was introduced [27 ]]The objective function measures the difference in probability distribution between the estimated and theoretical values x:
Figure 313147DEST_PATH_IMAGE151
(41)
assuming that the error is a binary distribution,
Figure 494729DEST_PATH_IMAGE152
the predicted probability distribution is considered to be closely related to the actual probability distribution, proving that the assumption is consistent with the expected model.
In step c, the step of establishing a basic performance test scene is as follows:
the basic performance of the proposed algorithm is tested, and a broadband dynamic signal model containing DDC components shown in formula (43) is constructed:
Figure 28610DEST_PATH_IMAGE153
(42)
wherein the content of the first and second substances,
Figure 194012DEST_PATH_IMAGE154
for the fundamental frequency, here set at 50Hz,
Figure 963385DEST_PATH_IMAGE155
、
Figure 265053DEST_PATH_IMAGE156
respectively represent the fundamental wave and each harmonic phase angle, and are in
Figure 219103DEST_PATH_IMAGE157
The value of the harmonic times h of the low frequency band is 2-13, the harmonic times h of the high frequency band are 77, 79, 80 and 83 respectively, the sampling frequency is set to be 10kHz, the amplitude of the DDC component is 0.3, 0.4, 0.5 and … 1, and the time constant of the DDC component is arbitrary
Figure 606222DEST_PATH_IMAGE158
Starting from 0.01s, changing from 0.01s to 0.1s in steps;
applying BMP algorithm, the initial value of the iterative phasor is
Figure 179286DEST_PATH_IMAGE159
Regularization coefficient of
Figure 335461DEST_PATH_IMAGE160
Fixed value parameter
Figure 398095DEST_PATH_IMAGE161
，
Figure 82629DEST_PATH_IMAGE162
，
Figure 256122DEST_PATH_IMAGE163
When the zeta value is in the range of [0.05,0.3 ]]Setting the number of the best matched data points in the search window as m;
FFT, Prony, TWLS, SIFE were chosen as comparison algorithms.
In step c, the step of establishing a frequency ramp test scenario is as follows:
the BMP was analyzed for performance under a frequency ramp, providing the following signals:
Figure 470065DEST_PATH_IMAGE164
(43)
wherein the content of the first and second substances,
Figure 703601DEST_PATH_IMAGE165
is the fundamental frequency and takes the value of 50Hz,
Figure 127629DEST_PATH_IMAGE166
is the slope of the fundamental frequency, the value is 1Hz/s,
Figure 104812DEST_PATH_IMAGE167
and
Figure 438841DEST_PATH_IMAGE168
a fundamental phase and a harmonic phase, respectively, the phases being set to
Figure 843278DEST_PATH_IMAGE169
Random numbers uniformly distributed within the range.
In step c, the step of establishing the step change test scene comprises the following steps:
at the start of the test, the amplitude of each component was set to 115% of the initial amplitude, and the phase was changed to
Figure 161127DEST_PATH_IMAGE170
The broadband dynamic signal provided is as follows:
Figure 755050DEST_PATH_IMAGE171
(44)
and setting the sampling rate to be 5kHz and the length of the sampling period to be 5 periods, and evaluating the performance of each algorithm under the condition of step change by using the speed of response time in the standard.
In step c, the step of establishing an interference test scene is as follows:
gaussian white noise with the signal-to-noise ratio of 55dB is introduced into the signal, and specific broadband dynamic signals are as follows:
Figure 5903DEST_PATH_IMAGE172
(45)
wherein the content of the first and second substances,
Figure 518924DEST_PATH_IMAGE173
is the inter-harmonic frequency, and the values thereof are 9652.5Hz, 9751.5Hz, 9850.5Hz, 9949.5Hz, 10048.5Hz and 10147.5Hz in sequence,
Figure 589648DEST_PATH_IMAGE174
and
Figure 970951DEST_PATH_IMAGE175
is the phase of fundamental wave and harmonic wave, and takes the value of
Figure 76310DEST_PATH_IMAGE176
The random number in the range is set to have a sampling frequency of 10kHz and a sampling window length of 5 power frequency periods.
The invention has the beneficial effects that:
a. the method is divided into two parts, wherein the first part obtains accurate DDC components based on a least square method, the second part establishes a Taylor-Fourier model of broadband dynamic phasors, the regularization optimization problem of a sparse acquisition model is solved by utilizing a harmonic vector estimation method, and finally, an estimated value of harmonic phasor measurement is obtained by utilizing a segmented SBI iteration frame, so that the reconstruction of an original signal is realized, the phasor measurement and estimation precision is obviously improved through a simulation test and a performance test, and a reliable theoretical basis can be provided for PMU measurement;
b. the invention firstly uses a sinc interpolation function to carry out parametric modeling on an attenuation direct current component, obtains an accurate DDC component by a least square method, adopts an improved Taylor-Fourier model aiming at broadband harmonic component estimation, combines a machine learning algorithm, accurately recovers a specific signal by fewer data points, further converts a harmonic phasor reconstruction problem into a sparse acquisition model solving problem, under the condition of a regularization frame, the invention expresses the phase reconstruction problem into a form of a minimization function, and effectively solves the underdetermined inverse problem by using alternate minimization operation under an SBI system, thereby efficiently and accurately obtaining a harmonic phasor estimation value, in order to prove the rationality of the model, an error range of the estimation value is measured by using a cross entropy objective function, and a simulation result shows that the estimation of the broadband dynamic phasor has higher precision and speed, certain superiority is kept in the comparison algorithm, and the dynamic performance and the anti-interference capability of the algorithm are verified;
c. the invention provides a new broadband dynamic phasor measurement algorithm, which solves the problem of effective processing of broadband dynamic phasors containing DDC components; based on the idea of regularizing the sparse capture matrix, obtaining a phasor estimation result by solving a sparse regularization problem; simulation and experiment results show that key information of the broadband dynamic phasor is identified, and the calculation complexity is obviously reduced; this represents the ability to obtain more accurate results in a short time, thereby effectively detecting the transient characteristics of a broadband signal; meanwhile, under static and dynamic conditions such as noise interference, inter-harmonic interference, frequency slope and the like, the measurement result can meet the test requirement of the M-level PMU.
Drawings
FIG. 1 is a diagram of the reconstruction effect and the operating time range as a function of parameters;
FIGS. 2[ a ], 2[ b ], 2[ c ] are graphs comparing maximum errors of TVE, FE, RFE under the condition of frequency ramp, FIG. 2[ a ] is a graph comparing maximum errors of TVE under the condition of frequency ramp, FIG. 2[ b ] is a graph comparing maximum errors of FE under the condition of frequency ramp, and FIG. 2[ c ] is a graph comparing maximum errors of RFE under the condition of frequency ramp;
FIG. 3 is a graph of the ratio of run times for the methods at step down;
fig. 4 a, 4 b, and 4 c are maximum error maps of TVE, FE, and RFE under inter-harmonic interference, fig. 4 a is a maximum error map of TVE and RFE under inter-harmonic interference, fig. 4 b is a maximum error map of FE under inter-harmonic interference, and fig. 4 c is a maximum error map of RFE under inter-harmonic interference.
Detailed Description
Embodiments of the present invention will be described in detail below with reference to the accompanying drawings.
Example (b):
1. estimation of attenuated DC components in broadband dynamic signals
1.1 broadband dynamic signal model
For power system waveforms containing a DC component and a fundamental component
Figure 760232DEST_PATH_IMAGE177
The expression is shown as (1):
Figure 318253DEST_PATH_IMAGE178
(1)
wherein
Figure 175350DEST_PATH_IMAGE179
Is a signal containing fundamental, harmonic, and inter-harmonic components.
Figure 213844DEST_PATH_IMAGE180
Is a time-varying amplitude of the signal,
Figure 130985DEST_PATH_IMAGE181
is a time-varying phase angle, the reference frequency is
Figure 113984DEST_PATH_IMAGE182
，
Figure 774773DEST_PATH_IMAGE183
Refers to the highest order harmonic component.
Figure 917041DEST_PATH_IMAGE184
Is to attenuate the signal of the direct current component,
Figure 5083DEST_PATH_IMAGE185
is to attenuate the magnitude of the dc component,
Figure 475379DEST_PATH_IMAGE186
is a time constant.
Signals containing fundamental, harmonic and inter-harmonic components
Figure 939858DEST_PATH_IMAGE187
The dynamic phasor model at time t is:
Figure 687365DEST_PATH_IMAGE188
(2)
substituting equation (2) into equation (1) to obtain broadband dynamic signal
Figure 946308DEST_PATH_IMAGE189
The dynamic expression of (2):
Figure 903900DEST_PATH_IMAGE190
(3)
wherein the content of the first and second substances,
Figure 172070DEST_PATH_IMAGE191
is a conjugate operator.
1.2, attenuated DC component estimation
Attenuating DC component
Figure 23352DEST_PATH_IMAGE192
Can be represented by the sum of the amplitude and phase time-varying cosine components:
Figure 453196DEST_PATH_IMAGE193
Figure 898084DEST_PATH_IMAGE194
(4)
wherein, the first and the second end of the pipe are connected with each other,
Figure 704366DEST_PATH_IMAGE195
，
Figure 82257DEST_PATH_IMAGE196
，
Figure 24281DEST_PATH_IMAGE197
respectively representing the frequency, the amplitude and the phase of the nth cosine component. And T is the observation interval length. To meet the accuracy requirement, the present embodiment sets Y to 3.
Based on the frequency domain sampling theorem, the time domain signal parameterization modeling is carried out on the attenuation direct current component by utilizing a sinc interpolation function, wherein the dynamic phasor of the DDC component is recorded as
Figure 487624DEST_PATH_IMAGE198
Specifically, it is represented as:
Figure 35280DEST_PATH_IMAGE199
(5)
wherein the content of the first and second substances,
Figure 533257DEST_PATH_IMAGE200
is the frequency of the sampling of the samples,
Figure 632800DEST_PATH_IMAGE201
represents the maximum model order, and the present embodiment will
Figure 849018DEST_PATH_IMAGE201
Set to 2.
Figure 934785DEST_PATH_IMAGE202
Is the nth cosine component of the formula (2) at the time
Figure 552849DEST_PATH_IMAGE203
The sample phasor above.
For attenuating DC component
Figure 574025DEST_PATH_IMAGE192
Sampling is performed, and the discrete form of the fitted DDC component in equation (2) can be represented by (5).
Figure 11960DEST_PATH_IMAGE204
(6)
Wherein
Figure 963735DEST_PATH_IMAGE205
Representing phasors
Figure 373988DEST_PATH_IMAGE206
Is the column vector of the cosine signal.
Figure 753017DEST_PATH_IMAGE207
Is a conjugate operator. Matrix of
Figure 6144DEST_PATH_IMAGE208
Is an element of
Figure 496031DEST_PATH_IMAGE209
A sample point in a sampling window.
Figure 26369DEST_PATH_IMAGE210
Representation matrix
Figure 310720DEST_PATH_IMAGE211
，
Figure 801875DEST_PATH_IMAGE212
The definition of (2) is shown as (7):
Figure 95454DEST_PATH_IMAGE213
(7)
determining phasor estimation values by means of least squares
Figure 480298DEST_PATH_IMAGE214
. The calculation process is shown as (9). In a matrix
Figure 935551DEST_PATH_IMAGE215
In (1),
Figure 835373DEST_PATH_IMAGE216
is the matrix element.
Figure 729380DEST_PATH_IMAGE217
(8)
Where H is the conjugate transposed symbol.
Each cosine component being constructed using equation (8)
Figure 296628DEST_PATH_IMAGE218
. Followed by
Figure 126043DEST_PATH_IMAGE218
Reconstruction matrix
Figure 513162DEST_PATH_IMAGE210
. The amplitude and the phase of each cosine signal are obtained through (5), and then DDC components are estimated
Figure 961592DEST_PATH_IMAGE206
。
2. Harmonic phasor estimation in broadband dynamic signals
2.1 establishment of harmonic phasor estimation model
It is contemplated that the taylor fourier method can describe the amplitude and phase transformation over time within the observation interval. Thus, the dynamic phasors are defined using a taylor expansion model of order:
Figure 383346DEST_PATH_IMAGE219
Figure 383664DEST_PATH_IMAGE220
(9)
wherein the content of the first and second substances,
Figure 258079DEST_PATH_IMAGE221
is t =0
Figure 493888DEST_PATH_IMAGE222
K is the taylor expansion order, and T is the observation interval length.
In practice, the sinusoidal signals containing harmonics and inter-harmonics are generally in the form of discrete sequences
Figure 504569DEST_PATH_IMAGE223
The sequence length is N, and
Figure 738104DEST_PATH_IMAGE224
. Through multi-frequency phasor analysis, a dynamic phasor method is expanded to an actual sequence model, Taylor Fourier coefficient expressions of each phasor component of a discrete signal are as follows:
Figure 37499DEST_PATH_IMAGE225
(10)
where H is the sampling interval, observation interval length
Figure 14682DEST_PATH_IMAGE226
。
Figure 967287DEST_PATH_IMAGE227
Is the average harmonic and inter-harmonic phasors,
Figure 371724DEST_PATH_IMAGE228
is at phasor frequency
Figure 158414DEST_PATH_IMAGE229
The k-th order derivative of (c).
Figure 939289DEST_PATH_IMAGE229
Is the fundamental frequency
Figure 658983DEST_PATH_IMAGE230
Integer multiples of (d) represent harmonic frequencies, and non-integer multiples represent inter-harmonic frequencies.
In the size of
Figure 296638DEST_PATH_IMAGE231
The sampling rate of (2) normalizes the frequency of the harmonic signal by a sampling length of
Figure 632941DEST_PATH_IMAGE232
Obtaining the frequency of the harmonic component
Figure 889610DEST_PATH_IMAGE233
The normalized frequency of the Taylor Fourier basis vector is
Figure 994969DEST_PATH_IMAGE234
，
Figure 554258DEST_PATH_IMAGE235
. The frequency resolution corresponding to the N coefficient sets is
Figure 377857DEST_PATH_IMAGE236
. Considering that the taylor opening order of each i is K, (11) can be written in the form of a matrix:
Figure 172638DEST_PATH_IMAGE237
(11)
wherein the content of the first and second substances,
Figure 132504DEST_PATH_IMAGE238
is the sample column phasor. A size of
Figure 377540DEST_PATH_IMAGE239
Is a taylor fourier-based phasor. e represents noise or measurement error. x is a length of
Figure 422857DEST_PATH_IMAGE240
The column phasors of (1), describe the current harmonics
Figure 818066DEST_PATH_IMAGE241
When the utility model is used, the water is discharged,
Figure 835700DEST_PATH_IMAGE242
a collection of (a).
Usually, the phasor solution can be converted into a solution of a pseudo-inverse matrix using a least squares method. The solution corresponding to the minimum Euclidean norm is obtained through the calculation of the pseudo-inverse matrix, and the solution is the optimal solution meeting the following constraint conditions:
Figure 923742DEST_PATH_IMAGE243
(12)
wherein the content of the first and second substances,
Figure 269404DEST_PATH_IMAGE244
representing the euclidean norm. Further, by introducing Lagrange operator and performing derivation, a dynamic phasor coefficient estimation matrix can be obtained
Figure 733883DEST_PATH_IMAGE245
：
Figure 606024DEST_PATH_IMAGE246
(13)
However, as the maximum harmonic order increases, the matrix dimension increases rapidly, and the amount of computation required to solve the multi-element linear equation corresponding to equation (14) increases significantly. For high frequency harmonics, the phasor method based on the least squares method is the least computationally complex. Therefore, the phasor method has a problem of low efficiency in solving the multidimensional matrix. In addition, since the observation time and the frequency resolution are contradictory, the frequency analysis is not accurate enough, and the time resolution is low.
2.2 reconstruction of harmonic phasor estimation model
To obtain more accurate frequency domain results, the present embodiment introduces a correction factor in the harmonic phasor analysis. Each group of frequency resolution reaches
Figure 864967DEST_PATH_IMAGE247
The length of the sampling sequence is
Figure 947193DEST_PATH_IMAGE248
The sampling length H is 5 fundamental periods. The harmonic component normalized frequency can be expressed as
Figure 949784DEST_PATH_IMAGE249
，
Figure 676432DEST_PATH_IMAGE250
。
Improved taylor fourier coefficient expressions are obtained under more precise correction conditions:
Figure 371855DEST_PATH_IMAGE251
(14)
here, the
Figure 613481DEST_PATH_IMAGE252
，
Figure 967233DEST_PATH_IMAGE253
Is expressed as a size of
Figure 876283DEST_PATH_IMAGE254
The matrix of the curvelet transform of (c),
Figure 414712DEST_PATH_IMAGE255
is that
Figure 409212DEST_PATH_IMAGE253
The k-th derivative of (c).
According to the expression (9), an approximate second-order dynamic phasor estimated value and a harmonic frequency can be obtained
Figure 815923DEST_PATH_IMAGE256
And frequency rate of change
Figure 579480DEST_PATH_IMAGE257
:
Figure 554389DEST_PATH_IMAGE258
(15)
Figure 770607DEST_PATH_IMAGE259
(16)
Figure 653112DEST_PATH_IMAGE260
(17)
Since the number of harmonic unknowns is much larger than the observed values, it is assumed that x is compressible in the curvelet transform domain. When correcting the factor
Figure 346874DEST_PATH_IMAGE261
When the sensing matrix A and the sparse matrix are in use
Figure 492684DEST_PATH_IMAGE262
Highly uncorrelated, (14) is converted into matrix form, namely:
Figure 196198DEST_PATH_IMAGE263
(18)
wherein
Figure 944711DEST_PATH_IMAGE264
Is a measurement matrix. x represents a length of
Figure 417281DEST_PATH_IMAGE265
The corresponding sample sequence of each data block of (a),
Figure 999572DEST_PATH_IMAGE262
is the superposition of a small number of phasors extracted after the original sample is segmented.
Figure 924803DEST_PATH_IMAGE266
Is dimension of
Figure 414690DEST_PATH_IMAGE267
Each column of the matrix is a fundamental phasor of a discrete fourier transform. e represents the noise signal phasor.
After the sparse sampling measurement model is constructed by the method, the matrix partitioning is performed on x by using the Euclidean search algorithm based on the sparsity of harmonic frequency domain distribution. Finding m best-matching data points within the search range and combining all data points into a matrix
Figure 554815DEST_PATH_IMAGE268
And i is the coordinate of the coefficient matrix. And performing curvelet transformation on the matrix, removing insignificant curvelet coefficients by adopting a curvelet threshold criterion, wherein the extracted coefficients need effective denoising and compression algorithms. And (4) performing curvelet transformation, and removing insignificant curvelet coefficients by adopting a curvelet threshold criterion, wherein the extracted coefficients need an effective denoising and compression algorithm.
To ensure locally smooth features for harmonic detection, sparse matrices
Figure 104746DEST_PATH_IMAGE269
Only a small set of non-zero valid entries is included. The phasor of the transformation coefficient column obtained after the effective item set is arranged according to the dictionary order is. Due to the transformation of coefficients
Figure 720535DEST_PATH_IMAGE270
The distribution of the transform coefficients is concentrated near the zero region and distributed in the form of a fine peak, so that the repeatability of the signal is characterized by using the distribution characteristics of the transform coefficients, and is recorded as follows:
Figure 14113DEST_PATH_IMAGE271
(19)
wherein
Figure 523591DEST_PATH_IMAGE272
Representing the L1 norm.
Each set of harmonic signal estimates may be inverted according to the transform coefficients:
Figure 978844DEST_PATH_IMAGE273
(20)
wherein
Figure 81929DEST_PATH_IMAGE274
Is that
Figure 913619DEST_PATH_IMAGE275
The inverse operator of (2).
2.3 solution of harmonic estimation model
In order to obtain an estimate with reasonable accuracy and noise robustness, after a harmonic phasor recovery model is established based on the above method, the following taylor coefficient matrix x is reconstructed by solving an optimization problem of a regularized linear least squares cost function, that is:
Figure 293915DEST_PATH_IMAGE276
(21)
Figure 920069DEST_PATH_IMAGE277
(22)
wherein, the first and the second end of the pipe are connected with each other,
Figure 510450DEST_PATH_IMAGE278
representing the L2 norm.
Figure 145831DEST_PATH_IMAGE279
A regularization term is represented as a function of,
Figure 36426DEST_PATH_IMAGE280
is a system of regularizationAnd (4) counting.
Figure 426957DEST_PATH_IMAGE281
、
Figure 239055DEST_PATH_IMAGE282
Representing local and non-local sparse terms, respectively.
Figure 412547DEST_PATH_IMAGE283
Is a regularization parameter that balances the two sparse terms.
(22) Substituting (21) to obtain:
Figure 688808DEST_PATH_IMAGE284
(23)
wherein
Figure 735392DEST_PATH_IMAGE285
。
To solve the above minimization problem, the present embodiment employs an alternative SBI algorithm. Introducing two auxiliary phasors p and q, the following scheme is finally obtained:
Figure 97103DEST_PATH_IMAGE286
(24)
Figure 11970DEST_PATH_IMAGE287
(25)
Figure 142737DEST_PATH_IMAGE288
(26)
Figure 609490DEST_PATH_IMAGE289
(27)
Figure 458498DEST_PATH_IMAGE290
(28)
wherein the content of the first and second substances,
Figure 911476DEST_PATH_IMAGE291
and
Figure 162328DEST_PATH_IMAGE292
the parameter is a fixed value parameter and has the function of improving the stability of the algorithm value. p, q are the SBI algorithm assisted iterative phasors.
Equation (24) is a problem of solving the minimization of a strictly convex quadratic function. Setting the gradient of the objective function to 0, resulting in a corresponding closed solution:
Figure 547786DEST_PATH_IMAGE293
(29)
wherein
Figure 884089DEST_PATH_IMAGE294
Is of size
Figure 937496DEST_PATH_IMAGE295
The identity matrix of (2).
The minimization problem of the function solved by the steepest descent method with the optimal step length can be expressed as follows:
Figure 980538DEST_PATH_IMAGE296
，
Figure 726777DEST_PATH_IMAGE297
(30)
wherein
Figure 347115DEST_PATH_IMAGE298
Is the gradient of the objective function, is the optimal step size.
Figure 469791DEST_PATH_IMAGE299
Taking into account the estimated values
Figure 367340DEST_PATH_IMAGE300
The present embodiment applies the minimum mean square error theorem to change the problem of equation (25) into:
Figure 284481DEST_PATH_IMAGE301
(31)
similarly, the question of (27) can be written as:
Figure 408426DEST_PATH_IMAGE302
(32)
here, the
Figure 69214DEST_PATH_IMAGE303
(omission of
Figure 821269DEST_PATH_IMAGE304
）。
Due to the fact that
Figure 909311DEST_PATH_IMAGE305
The definition of (A) is complicated, and the above formula is difficult to be solved intuitively. This embodiment derives a set of closed solutions by proposing reasonable assumptions. Will be provided with
Figure 441924DEST_PATH_IMAGE306
Of a type of observation viewed as noisy, using
Figure 968720DEST_PATH_IMAGE307
Representing the error phasor, the error of each element is respectively
Figure 637599DEST_PATH_IMAGE308
（
Figure 99804DEST_PATH_IMAGE309
) Then assume e is independently distributed over
Figure 119713DEST_PATH_IMAGE310
Mean 0 and variance
Figure 935353DEST_PATH_IMAGE311
From law of large numbers, for arbitrary
Figure 724317DEST_PATH_IMAGE312
All have:
Figure 357424DEST_PATH_IMAGE313
(33)
so that the following results are obtained:
Figure 333470DEST_PATH_IMAGE314
(34)
the transformed error phasor is
Figure 139752DEST_PATH_IMAGE315
，
Figure 111119DEST_PATH_IMAGE316
Figure 711865DEST_PATH_IMAGE317
Representing the error phasor for each element, and m is the number of data points for the best match. The transformation does not change the variance of each group according to the orthogonal nature of the matrix. It follows therefrom that of each group
Figure 644049DEST_PATH_IMAGE318
Are all independently distributed in
Figure 988443DEST_PATH_IMAGE319
(zero mean and variance of
Figure 565049DEST_PATH_IMAGE320
) Most of them are theoretic and arbitrary
Figure 336696DEST_PATH_IMAGE321
All have:
Figure 756176DEST_PATH_IMAGE322
(35)
the method comprises the following steps:
Figure 638681DEST_PATH_IMAGE323
(36)
it can be easily obtained that:
Figure 256744DEST_PATH_IMAGE324
(37)
combining the problem with the formula (26) to obtain:
Figure 527188DEST_PATH_IMAGE325
(38)
due to unknown quantity
Figure 230702DEST_PATH_IMAGE326
The formula is separable, and each component is simplified into:
Figure 854582DEST_PATH_IMAGE327
(39)
wherein
Figure 592730DEST_PATH_IMAGE328
Refers to the modulus of the corresponding phasor.
Under the vector contraction form, the estimated value of x is expressed as:
Figure 516300DEST_PATH_IMAGE329
(40)
after the estimation value of the Taylor coefficient matrix x is obtained, the method is carried out
Figure 707110DEST_PATH_IMAGE330
A harmonic sampling equation for realizing the harmonic signal
Figure 134680DEST_PATH_IMAGE331
The purpose of accurately detecting the harmonic phasor change is achieved. In order to ensure the accuracy of the model, under the framework of logistic regression, a cross entropy objective function is introduced to measure the probability distribution difference between an estimated value and a theoretical value x:
Figure 727335DEST_PATH_IMAGE332
(41)
assuming that the error is a binary distribution,
Figure 339582DEST_PATH_IMAGE333
the predicted probability distribution can be considered to be very closely related to the actual probability distribution, which demonstrates that the assumption is consistent with the expected model. The algorithm flow is shown in the following table:
Figure 17688DEST_PATH_IMAGE334
3. simulation analysis
This section discusses the result of the broadband dynamic phasor measurement algorithm (BMP) proposed in this embodiment under different test scenarios. The following 4 algorithms are used as comparison algorithms: fast fourier method (FFT), Prony algorithm, Taylor Weighted Least Squares (TWLS), dynamic synchrophasor measurement algorithm (SIFE) based on sinc interpolation function. The test scenario includes a basic performance test, an interference test, a frequency ramp test, and a step change test.
To compare the estimation effects of the different methods, a total phasor error (TVE) was introduced to describe the relative deviation between the theoretical phasor and the estimated phasor. The TVE is closely related to the amplitude error and the phase angle error, but cannot reflect the change of one aspect alone. Therefore, the present embodiment also introduces two indicators: the Frequency Error (FE) and the absolute frequency rate error (RFE) are used to evaluate the phasor estimation effect. In this embodiment, according to the IEEEC37.118.1 standard, the M-type PMU measurement requirement is selected as a standard limit value under different working conditions.
The 5 algorithms all use the same rectangular observation window, the bandwidth of the fundamental frequency is 1Hz, and the length of the sampling window is set as 5 power frequency periods. In order to improve the frequency domain resolution, the present embodiment considers the complexity of the algorithm and the accuracy of the algorithm together, and sets the value to 20.
4.1, basic Performance testing
To test the basic performance of the proposed algorithm, a broadband dynamic signal model containing DDC components is constructed as shown in equation (42):
Figure 983370DEST_PATH_IMAGE335
(42)
wherein, the first and the second end of the pipe are connected with each other,
Figure 430532DEST_PATH_IMAGE336
for the fundamental frequency, here set at 50Hz,
Figure 885784DEST_PATH_IMAGE337
、
Figure 864235DEST_PATH_IMAGE338
respectively represent the fundamental wave and each harmonic phase angle, and are in
Figure 695925DEST_PATH_IMAGE339
Any value within the range. The low-frequency band harmonic times h are 2-13, the high-frequency band harmonic times h are 77, 79, 80 and 83 respectively, and the sampling frequency is set to be 10 kHz. DDC component magnitude
Figure 935277DEST_PATH_IMAGE340
Is 0.3, 0.4, 0.5, … 1. Time constant of DDC component
Figure 827009DEST_PATH_IMAGE341
Starting from 0.01s, it changes from 0.01s to 0.1s in steps.
Iterative phasor initial values when applying the BMP algorithm
Figure 542024DEST_PATH_IMAGE342
Regularization coefficient
Figure 911826DEST_PATH_IMAGE343
. Fixed value parameter
Figure 68001DEST_PATH_IMAGE344
，
Figure 68318DEST_PATH_IMAGE345
，
Figure 942733DEST_PATH_IMAGE346
When used, their zeta value is generally in a range [0.05,0.3 ]]An internal change. The number of data points in the search window that match best is set to m. In the process of matrix blocking, if the number of matched blocks is too large, data points with large noise influence and low matching degree inevitably exist in the block array; otherwise, the influence of the contingency on the construction matrix cannot be avoided. In this embodiment, the maximum iteration number is J, and the higher the iteration number is, the higher the calculation accuracy is, but the calculation cost is significantly increased. In this embodiment, the reconstruction effect and the algorithm running time are analyzed under the condition that the parameters m and J are changed, and the result is shown in fig. 1.
FIG. 1 is a graph of reconstruction effect and run-time range as a function of parameters. As can be seen from fig. 1, as the parameters m and J increase, the total phasor error shows a decreasing trend, but the algorithm runtime increases. When the number of data points m =240 and the number of iterations J is less than 100, changing the number of iterations has less influence on the calculation cost (the change in the running time is small) and the reconstruction effect tends to be stable at this time.
Figure 194854DEST_PATH_IMAGE347
In time, the increasing number of iterations leads to a great increase in computational complexity (large variation in running time), and the TVE also exhibits an unstable downward trend. Therefore, appropriately reducing m can effectively improve the reconstruction effect. Considering both reconstruction performance and runtime, the present embodiment applies the algorithm with the parameter m =240 and the iteration number J = 100.
FFT, Prony, TWLS and SIFE are selected as comparison algorithms, the total phasor error estimation and frequency change rate error estimation results of the comparison algorithms are shown in the following table, and the estimation precision of each algorithm is compared with that of the following table:
Figure 471114DEST_PATH_IMAGE349
Figure 642333DEST_PATH_IMAGE351
Figure 4044DEST_PATH_IMAGE353
as can be seen from the above table, the maximum values of the TVE, FE and RFE indices of the BMP algorithm are 2.79%, 0.096Hz and 2.437Hz/s, respectively. When harmonic components are included, the IEEE standard limits for TVE and FE are 3% and 0.1Hz, respectively, and the RFE limit is 2.7 Hz/s. Compared with other comparison algorithms, the estimation index of the invention completely meets the requirements of the IEEE standard. The result shows that the method still has good detection effect under the condition of the DDC component-containing broadband harmonic wave, and the phasor estimation precision is highest. The BMP uses a sparse distribution of harmonic frequency domain distributions to identify the most relevant components in the signal, which significantly improves the accuracy of the measurement results.
The estimation error of the SIFE method in the vicinity of the low harmonic does not meet the measurement standard of IEEE. This is because the SIFE method is a low pass filter based on baseband signal filtering, and it is difficult to obtain a zero error result. But because the harmonic synchronous phasor can be well estimated due to the wide passband and the wide stopband, the performance of the invention is superior to FFT and TWLS. The FFT measurement results are most affected by spectral leakage, which greatly reduces the accuracy of harmonic parameter identification. The maximum total phasor error exceeds 8%, and the accuracy of the FE and RFE measurement results is not ideal. Under dynamic conditions, the fourier transform model cannot track phasor changes in the observation window, resulting in incorrect phasor evaluation. The TWLS method uses a second order taylor order to fit the signal components. However, the taylor signal model has a large error and a limited accuracy. Increasing the taylor model order can reduce the model error, but the higher the order, the worse the passband performance of the filter. The Prony algorithm uses a parametric model to calculate the signal parameters. However, the estimation order thereof will limit the number of estimated frequency components, and the frequency estimation error will gradually increase.
4.2 frequency ramp test
A power imbalance between the load and the generator may cause the frequency of the broadband signal to decrease with increasing load and increase with increasing input power. To analyze the performance of BMP (the algorithm of this example) in a frequency ramp, the signals provided are as follows:
Figure 777965DEST_PATH_IMAGE354
(42)
wherein the content of the first and second substances,
Figure 174311DEST_PATH_IMAGE355
is the fundamental frequency and takes the value of 50 Hz.
Figure 516431DEST_PATH_IMAGE356
The value of the present embodiment is 1 Hz/s.
Figure 365438DEST_PATH_IMAGE357
And
Figure 693782DEST_PATH_IMAGE358
a fundamental phase and a harmonic phase, respectively, the phases being set to
Figure 944635DEST_PATH_IMAGE359
Random numbers uniformly distributed within the range. The harmonic phasor estimation, frequency estimation and frequency change rate estimation results of the present embodiment and the comparative algorithm are shown in fig. 1. Assume a sampling frequency of 10kHz and a sampling window length of 5 power frequency cycles. Analysis of the signal of (42) using FFT, Prony, TWLS, SIFE as a comparison algorithm, FIG. 2[ a ]]FIG. 2[ b ]]FIG. 2[ c ]]FIG. 2[ a ] is the estimation results generated by each method under the frequency ramp condition]FIG. 2[ b ]]FIG. 2[ c ]]FIG. 2[ a ] is a graph comparing the maximum error of TVE, FE, RFE under the condition of frequency ramp]FIG. 2[ b ] is a graph of maximum error versus TVE under frequency ramp conditions]FIG. 2[ c ] is a graph comparing the maximum error of FE under the frequency ramp condition]For RFE under frequency ramp conditionsMaximum error vs. plot.
As can be seen from fig. 2 a, 2 b, and 2 c, the maximum error of TVE in this embodiment is 1.06%, the maximum error of frequency is 0.0095Hz, and the maximum error of frequency conversion rate is 0.103Hz/s, which all satisfy the IEEE standard requirement [14 ]. The result shows that the method can still keep higher precision when the fundamental frequency changes greatly and the frequency changes linearly, and the estimation effect is superior to other four algorithms. The BMP algorithm adopts a Taylor order to approximate a dynamic signal model, can accurately estimate the frequency and the phase angle, and is minimally influenced by linear frequency change. The method achieves the highest phasor estimation accuracy. In other algorithms, the FFT method cannot track the frequency change in real time under dynamic conditions, and thus the frequency change rate is large. Since TWLS is an estimation algorithm based on a dynamic model, Prony algorithm can accurately extract the low-frequency oscillation characteristic value of the dominant mode, the error calculation results of Prony and TWLS algorithms are smaller than FFT, the error characteristics of parameters such as phase angle and frequency of the algorithms are slightly influenced by frequency change, but the measurement results still cannot meet the requirements of IEEE. The Prony algorithm has difficulty obtaining accurate mode parameters if the time-varying characteristics of the signal are not well understood. The STWLS also cannot accurately estimate the higher harmonic phasor due to the influence of the fundamental component.
4.3 step Change test
In order to simulate fault conditions of sudden changes in the amplitude and phase of the voltage/current signal, it is necessary to simulate the proposed algorithm under these conditions to evaluate the response time and delay. At the start of the test, the amplitude of each component was set to 115% of the initial amplitude, and the phase was changed to
Figure 457656DEST_PATH_IMAGE360
. The broadband dynamic signal provided is as follows:
Figure 528380DEST_PATH_IMAGE361
(43)
the test in this section assumes a sampling rate of 5kHz and a sampling period of 5 cycles. In the standard, the performance of each algorithm under a step change condition is evaluated by the speed of the response time. It is defined as the time interval between the first and last instants greater than a given threshold. The threshold values for the maximum TVE, FE and RFE values were 1.5%, 0.13Hz, and 0.78Hz/s according to the IEEE Standard under the test conditions in this section. The specific simulation results are shown in fig. 3.
As can be seen from fig. 3, the BMP method requires less time to reach the TVE, FE and RFE values of the IEEE standard under the amplitude and phase step change conditions than other algorithms. It can be considered that under the condition of step change, the estimation precision of the invention is the highest and the response time completely meets the standard requirement of P-type PMU. For the evaluation of the entire fourier transform of the step-wise smooth function, the FFT algorithm requires very complex multiplication operations. TWLS performance can be improved by recalculating its coefficients in each report frame, however this greatly increases the computational burden, requiring the computation of the pseudo-inverse. In order to achieve the precision requirement of the SIFE algorithm, a complex matrix of 70 rows and 1400 columns needs to be stored in a memory, and the sub-actual multiplication and the sub-actual addition need to be processed in real time, so that the memory and the processing capacity are very limited. In order to obtain more accurate analysis results by the Prony algorithm, the model order of the method needs to be increased, so that the calculated amount is increased.
4.4, anti-interference test
In general, some inter-harmonics and noise are contained in the power system signal, which seriously affect the estimation of the harmonic phasor. In this section of the test, white gaussian noise with a signal to noise ratio of 55dB was introduced into the signal. Specific broadband dynamic signals are as follows:
Figure 847366DEST_PATH_IMAGE362
(43)
wherein, the first and the second end of the pipe are connected with each other,
Figure 15042DEST_PATH_IMAGE363
is the inter-harmonic frequency, and the values thereof are 9652.5Hz, 9751.5Hz, 9850.5Hz, 9949.5Hz, 10048.5Hz and 10147.5Hz in sequence.
Figure DEST_PATH_IMAGE364
The sum is the phase of the fundamental wave and the harmonic wave, and takes the value of
Figure DEST_PATH_IMAGE365
Random numbers within a range. This section assumes a sampling frequency of 10kHz and a sampling window length of 5 power frequency cycles. The specific simulation results are shown in FIG. 4[ a ]]FIG. 4[ b ]]FIG. 4[ c ]]Shown in FIG. 4[ a ]]FIG. 4[ b ]]FIG. 4[ c ]]FIG. 4 is a graph of maximum error of TVE, FE and RFE under inter-harmonic interference [ a ]]FIG. 4[ b ] is a graph of the maximum error of TVE under inter-harmonic interference]FIG. 4[ c ] is a graph of the maximum error of FE under inter-harmonic interference]Is a maximum error map of RFE under inter-harmonic interference.
As can be seen from FIGS. 4[ a ], 4[ b ] and 4[ c ], the TVEmax of the algorithm of the present embodiment is 2.43%, FEmax is 0.071Hz, and RFEmax is 0.184 Hz/s. The calculation result of each index of the algorithm of the embodiment is superior to that of the other comparison algorithms. Under the condition of the existence of inter-harmonic interference, the IEEE standard limit values of TVE, FE and RFE are respectively 3.5%, 0.2Hz and 3 Hz/s. In the same way, the parameter results of the SIFE algorithms TVEmax, FEmax and RFEmax are respectively 6.64%, 0.115Hz and 7.39Hz/s, the Prony algorithms are 19.73%, 0.136Hz and 9.207Hz/s, the TWLS algorithms are 26.35%, 0.171Hz and 19.634Hz/s, and the FFT algorithms are 49.76%, 17.294Hz and 24.270Hz/s, and the calculation results of all indexes of the invention are superior to other comparison algorithms. Each index of the BMP algorithm meets the requirements of the IEEE standard. In the presence of noise interference in this segment, severe interference and spectral leakage can occur between adjacent harmonics of the FFT method, affecting resolution and accuracy. Under the condition of noise interference in the present segment, the TWLS and SIFE methods can increase the transition band amplitude of their harmonic filters, causing interference between adjacent harmonics, resulting in larger estimation error. In addition, the TWLS algorithm is severely affected by inter-harmonic interference, and it is difficult to accurately estimate each high-frequency component. For high-order components, the model parameters of the Prony algorithm are continuously modified along with the increase of the harmonic order. Therefore, it has a good effect in inter-harmonic phasor and frequency estimation, and can suppress the influence of the inter-harmonic component spectrum leakage to some extent. However, this method requires the order of the dynamic time-varying signal to be estimated in advance.
The BMP algorithm adopts an improved Taylor-Fourier model and a machine learning algorithm, accurately recovers a specific signal by fewer data points, reconstructs harmonic phasors and solves the specific signal by a sparse acquisition model, and effectively improves the reconstruction performance and the anti-noise capability of the algorithm. In addition, the spectral resolution of the reconstruction of the invention is 1Hz, which is beneficial to the accurate detection of the inter-harmonic component.
4. Conclusion
The harmonic detection method based on the sparse acquisition model provided by the embodiment is used for solving the problem of effective processing of the broadband dynamic phasor containing the DDC component. The phasor estimation method is based on the idea of regularizing the sparse capture matrix, and obtains phasor estimation results by solving a sparse regularization problem. Simulation and experiment results show that the method identifies key information of the broadband dynamic phasor, and obviously reduces the computational complexity. This represents the ability to obtain more accurate results in a short time, thereby effectively detecting the transient characteristics of a broadband signal. Meanwhile, under static and dynamic conditions such as noise interference, inter-harmonic interference, frequency slope and the like, the measurement result can meet the test requirements of the M-level PMU.
The above-mentioned embodiments only express the specific embodiments of the present invention, and the description thereof is more specific and detailed, but not construed as limiting the scope of the present invention. It should be noted that, for a person skilled in the art, several variations and modifications can be made without departing from the inventive concept, which falls within the scope of the present invention.

Claims (10)

1. A harmonic detection method based on a sparse acquisition model is characterized by comprising the following steps:
a. estimating the attenuation direct current component in the broadband dynamic signal: establishing a broadband dynamic signal model, and calculating an attenuated direct current component;
b. estimating harmonic phasor in the broadband dynamic signal: establishing a harmonic phasor estimation model, reconstructing the harmonic phasor estimation model, and solving the harmonic estimation model;
c. a simulation analysis step: and establishing a test scene, and sequentially establishing a basic performance test scene, a frequency slope test scene, a step change test scene and an interference test scene.
2. The harmonic detection method based on the sparse acquisition model as claimed in claim 1, wherein in the step a, the step of establishing the broadband dynamic signal model is as follows:
for power system waveforms containing a DC component and a fundamental component
Figure 445225DEST_PATH_IMAGE001
The expression is shown as (1):
Figure 505585DEST_PATH_IMAGE002
(1)
wherein the content of the first and second substances,
Figure 871975DEST_PATH_IMAGE003
is a signal containing fundamental, harmonic and inter-harmonic components,
Figure 992378DEST_PATH_IMAGE004
is a time-varying amplitude of the signal,
Figure 240957DEST_PATH_IMAGE005
is a time-varying phase angle, the reference frequency is
Figure 105007DEST_PATH_IMAGE006
，
Figure 624107DEST_PATH_IMAGE007
Refers to the highest order harmonic component of the signal,
Figure 915411DEST_PATH_IMAGE008
is to attenuate the signal of the direct current component,
Figure 651286DEST_PATH_IMAGE009
is the amplitude of the attenuated dc component,
Figure 53448DEST_PATH_IMAGE010
is a time constant;
signals containing fundamental, harmonic and inter-harmonic components
Figure 191168DEST_PATH_IMAGE011
The dynamic phasor model at time t is:
Figure 387795DEST_PATH_IMAGE012
(2)
substituting equation (2) into equation (1) to obtain broadband dynamic signal
Figure 109501DEST_PATH_IMAGE013
The dynamic expression of (2):
Figure 315354DEST_PATH_IMAGE014
(3)
wherein, the first and the second end of the pipe are connected with each other,
Figure 42002DEST_PATH_IMAGE015
is a conjugate operator.
3. The harmonic detection method based on the sparse acquisition model as claimed in claim 2, wherein in the step a, the step of calculating the attenuated direct current component is as follows:
attenuating DC component
Figure 940688DEST_PATH_IMAGE016
Represented by the sum of the amplitude and phase time-varying cosine components:
Figure 651155DEST_PATH_IMAGE017
Figure 129540DEST_PATH_IMAGE018
(4)
wherein, the first and the second end of the pipe are connected with each other,
Figure 241853DEST_PATH_IMAGE019
，
Figure 541466DEST_PATH_IMAGE020
，
Figure 739230DEST_PATH_IMAGE021
respectively representing the frequency, the amplitude and the phase of the nth cosine component, wherein T is the length of an observation interval, and Y is set to be 3;
based on frequency domain sampling theorem, a sine interpolation function is utilized to carry out time domain signal parameterization modeling on the attenuation direct current component, and the dynamic phasor of the DDC component is recorded as
Figure 21306DEST_PATH_IMAGE022
Specifically, it is represented as:
Figure 988125DEST_PATH_IMAGE023
(5)
wherein, the first and the second end of the pipe are connected with each other,
Figure 228614DEST_PATH_IMAGE024
is the frequency of the sampling of the samples,
Figure 382515DEST_PATH_IMAGE025
indicating the maximum model order, will be set to 2,
Figure 733861DEST_PATH_IMAGE026
is the nth cosine component in the formula (2) at the time
Figure 53722DEST_PATH_IMAGE027
A sample phasor above;
for the attenuation of DC component
Figure 465112DEST_PATH_IMAGE028
Sampling is carried out, and the discrete form of the fitted DDC component in the formula (1) is represented by (6);
Figure 371888DEST_PATH_IMAGE029
(6)
wherein, the first and the second end of the pipe are connected with each other,
Figure 995767DEST_PATH_IMAGE030
representing phasors
Figure 937179DEST_PATH_IMAGE031
Is a column vector of the cosine signal,
Figure 519470DEST_PATH_IMAGE015
as conjugate operators, matrices
Figure 913542DEST_PATH_IMAGE032
Is an element of
Figure 606691DEST_PATH_IMAGE033
At a sample point in the sampling window(s),
Figure 372915DEST_PATH_IMAGE034
representation matrix
Figure 126108DEST_PATH_IMAGE035
，
Figure 7476DEST_PATH_IMAGE036
The definition of (2) is shown as (7):
Figure 504317DEST_PATH_IMAGE037
(7)
determining phasor estimation values by means of least squares
Figure 889162DEST_PATH_IMAGE038
The calculation process is shown in (7) in the matrix
Figure 282097DEST_PATH_IMAGE039
In (1),
Figure 149296DEST_PATH_IMAGE040
for the matrix elements:
Figure 184248DEST_PATH_IMAGE041
(8)
where H is a conjugate transposed symbol;
each cosine component being constructed using equation (8)
Figure 954758DEST_PATH_IMAGE042
Followed by using
Figure 518595DEST_PATH_IMAGE042
Reconstruction matrix
Figure 374555DEST_PATH_IMAGE043
The amplitude and the phase of each cosine signal are obtained through the step (5), and then DDC components are estimated
Figure 213198DEST_PATH_IMAGE044
。
4. The harmonic detection method based on the sparse acquisition model as claimed in claim 1, wherein in the step b, the step of establishing the harmonic phasor estimation model is as follows:
the dynamic phasors are defined using a taylor expansion model of order:
Figure 572636DEST_PATH_IMAGE045
Figure 838532DEST_PATH_IMAGE046
(9)
wherein the content of the first and second substances,
Figure 146236DEST_PATH_IMAGE047
is t =0
Figure 522990DEST_PATH_IMAGE048
K is the taylor expansion order, and T is the observation interval length;
the sinusoidal signals containing harmonics and inter-harmonics take the form of discrete sequences
Figure 2513DEST_PATH_IMAGE049
The sequence length is N, and
Figure 173731DEST_PATH_IMAGE050
and expanding the dynamic phasor method to an actual sequence model through multi-frequency phasor analysis, wherein the Taylor Fourier coefficient expression of each phasor component of the discrete signal is as follows:
Figure 738705DEST_PATH_IMAGE051
(10)
where H is the sampling interval, observation interval length
Figure 653571DEST_PATH_IMAGE052
，
Figure 627081DEST_PATH_IMAGE053
Is the average harmonic and inter-harmonic phasors,
Figure 969201DEST_PATH_IMAGE054
is a phasor frequency of
Figure 21471DEST_PATH_IMAGE055
The derivative of the order k of (c) is,
Figure 975914DEST_PATH_IMAGE055
is the fundamental frequency
Figure 164449DEST_PATH_IMAGE056
Integer multiples of (d) represent harmonic frequencies, and non-integer multiples represent inter-harmonic frequencies;
in the size of
Figure 943050DEST_PATH_IMAGE057
The sampling rate of (2) normalizes the frequency of the harmonic signal by a sampling length of
Figure 217036DEST_PATH_IMAGE058
To obtain the frequency of harmonic component
Figure 739284DEST_PATH_IMAGE059
The normalized frequency of the Taylor Fourier basis vector is
Figure 782327DEST_PATH_IMAGE060
，
Figure 731828DEST_PATH_IMAGE061
The frequency resolution corresponding to the N coefficient sets is
Figure 726067DEST_PATH_IMAGE062
And the taylor opening order of each i is K, then (10) is written in matrix form:
Figure 786427DEST_PATH_IMAGE063
(11)
wherein, the first and the second end of the pipe are connected with each other,
Figure 683975DEST_PATH_IMAGE064
is a sample column phasor of size
Figure 804378DEST_PATH_IMAGE065
Is a Taylor Fourier-based phasor, e represents noise or measurement error, and x is a length of
Figure 318536DEST_PATH_IMAGE066
The column phasors of (1), describe the current harmonics
Figure 917008DEST_PATH_IMAGE067
When the temperature of the water is higher than the set temperature,
Figure 188106DEST_PATH_IMAGE068
a set of (a);
and (3) phasor solving, namely converting the phasor solving into solving of a pseudo-inverse matrix by using a least square method, and calculating the pseudo-inverse matrix to obtain a corresponding solution when the Euclidean norm is minimum, wherein the solution is the optimal solution under the following constraint conditions:
Figure 744989DEST_PATH_IMAGE069
(12)
wherein the content of the first and second substances,
Figure 480864DEST_PATH_IMAGE070
representing Euclidean norm, and further obtaining a dynamic phasor coefficient estimation matrix after introducing a Lagrange operator and derivation
Figure 148606DEST_PATH_IMAGE071
：
Figure 755168DEST_PATH_IMAGE072
(13)。
5. The harmonic detection method based on the sparse acquisition model as claimed in claim 4, wherein in the step b, the step of reconstructing the harmonic phasor estimation model is as follows:
setting a correction factor
Figure 482952DEST_PATH_IMAGE073
Each group of frequency resolution reaches
Figure 939079DEST_PATH_IMAGE074
The length of the sampling sequence is
Figure 144933DEST_PATH_IMAGE075
The sampling length H is 5 fundamental wave periods, and the harmonic component normalized frequency is
Figure 137159DEST_PATH_IMAGE076
，
Figure 35845DEST_PATH_IMAGE077
And obtaining an improved Taylor Fourier coefficient expression:
Figure 480733DEST_PATH_IMAGE078
(14)
wherein, the first and the second end of the pipe are connected with each other,
Figure 224698DEST_PATH_IMAGE079
，
Figure 337011DEST_PATH_IMAGE080
is expressed as a size of
Figure 376904DEST_PATH_IMAGE081
The matrix of the curvelet transform of (c),
Figure 309088DEST_PATH_IMAGE082
is that
Figure 856744DEST_PATH_IMAGE083
The k-th derivative of (a);
according to the expression (9), an approximate second-order dynamic phasor estimated value and harmonic frequency are obtained
Figure 823563DEST_PATH_IMAGE084
And rate of change of frequency
Figure 532893DEST_PATH_IMAGE085
:
Figure 217952DEST_PATH_IMAGE086
(15)
Figure 569299DEST_PATH_IMAGE087
(16)
Figure 623581DEST_PATH_IMAGE088
(17)
When the correction factor is
Figure 34970DEST_PATH_IMAGE089
When the sensing matrix A and the sparse matrix are in use
Figure 676167DEST_PATH_IMAGE090
Highly uncorrelated, (14) is converted into matrix form, namely:
Figure 565626DEST_PATH_IMAGE091
(18)
wherein, the first and the second end of the pipe are connected with each other,
Figure 507037DEST_PATH_IMAGE092
is a measurement matrix, x represents a length of
Figure 292590DEST_PATH_IMAGE093
The corresponding sample sequence of each data block of (a),
Figure 182268DEST_PATH_IMAGE090
is the superposition of a small amount of phasors extracted after the original sample is segmented,
Figure 609838DEST_PATH_IMAGE094
is dimension of
Figure 405756DEST_PATH_IMAGE095
Each column of the matrix is a base phasor of a discrete fourier transform, e represents a noise signal phasor;
after the sparse sampling measurement model is constructed by the method, matrix blocking is carried out on x by utilizing an Euclidean search algorithm, m data points which are optimally matched are found in a search range, and all the data points are combined into one matrix
Figure 158949DEST_PATH_IMAGE096
I is a coordinate of the coefficient matrix, curvelet transformation is carried out on the matrix, an insignificant curvelet coefficient is removed by adopting a curvelet threshold criterion, an effective denoising and compression algorithm is required for the extracted coefficient, curvelet transformation is carried out, an insignificant curvelet coefficient is removed by adopting a curvelet threshold criterion, and an effective denoising and compression algorithm is required for the extracted coefficient;
sparse matrix
Figure 40317DEST_PATH_IMAGE090
Only contains a non-zero significant term set, the transformation coefficient column phasor obtained after the significant term set is arranged according to the dictionary order is
Figure 5999DEST_PATH_IMAGE097
Transformation coefficient
Figure 656423DEST_PATH_IMAGE098
Is focused onThe distribution is near the zero zone and in the form of a fine peak, and the repeatability of the signal is represented by using the distribution characteristics of the transformation coefficient, and is recorded as follows:
Figure 813473DEST_PATH_IMAGE099
(19)
wherein, the first and the second end of the pipe are connected with each other,
Figure 182137DEST_PATH_IMAGE100
represents the L1 norm;
each set of harmonic signal estimates may be inverted according to transform coefficients:
Figure 951510DEST_PATH_IMAGE101
(20)
wherein the content of the first and second substances,
Figure 456441DEST_PATH_IMAGE102
is that
Figure 551436DEST_PATH_IMAGE103
The inverse operator of (2).
6. The harmonic detection method based on the sparse acquisition model as claimed in claim 5, wherein in the step b, the step of solving the harmonic estimation model is as follows:
solving a regularized linear least square cost function, and reconstructing a Taylor coefficient matrix x, namely:
Figure 141817DEST_PATH_IMAGE104
(21)
Figure 714881DEST_PATH_IMAGE105
(22)
wherein the content of the first and second substances,
Figure 310203DEST_PATH_IMAGE106
the norm of L2 is shown,
Figure 310521DEST_PATH_IMAGE107
a regularization term is represented as a function of,
Figure 388198DEST_PATH_IMAGE108
is a function of the regularization coefficients and,
Figure 30532DEST_PATH_IMAGE109
、
Figure 244476DEST_PATH_IMAGE110
representing local and non-local sparse terms respectively,
Figure 681273DEST_PATH_IMAGE111
is a regularization parameter that balances the two sparse terms;
(22) substituting (21) to obtain:
Figure 744782DEST_PATH_IMAGE112
(23)
wherein the content of the first and second substances,
Figure 659648DEST_PATH_IMAGE113
，
introducing two auxiliary phasors p and q, finally the following scheme is obtained:
Figure 993678DEST_PATH_IMAGE114
(24)
Figure 601376DEST_PATH_IMAGE115
(25)
Figure 388067DEST_PATH_IMAGE116
(26)
Figure 841045DEST_PATH_IMAGE117
(27)
Figure 790765DEST_PATH_IMAGE118
(28)
wherein the content of the first and second substances,
Figure 569366DEST_PATH_IMAGE119
and
Figure 108931DEST_PATH_IMAGE120
the parameters are fixed value parameters, the function is to improve the stability of the algorithm value, and p and q are SBI algorithm auxiliary iteration phasors;
(24) solving the minimization of a strict convex quadratic function, setting the gradient of an objective function as 0, and obtaining a corresponding closed solution:
Figure 365600DEST_PATH_IMAGE121
(29)
wherein, the first and the second end of the pipe are connected with each other,
Figure 674222DEST_PATH_IMAGE122
is of the size
Figure 623723DEST_PATH_IMAGE123
The identity matrix of (a);
the minimization problem of the function is solved by adopting the steepest descent method with the optimal step length, which is expressed as follows:
Figure 385006DEST_PATH_IMAGE124
，
Figure 943901DEST_PATH_IMAGE125
(30)
wherein the content of the first and second substances,
Figure 107029DEST_PATH_IMAGE126
is the gradient of the objective function and,
Figure 961853DEST_PATH_IMAGE127
is the optimum step size for the particular application,
Figure 476011DEST_PATH_IMAGE128
the estimated value is
Figure 340061DEST_PATH_IMAGE129
Applying the least mean square error theorem changes the problem of (25) into:
Figure 92117DEST_PATH_IMAGE130
(31)
similarly, the question of (26) is written as:
Figure 649000DEST_PATH_IMAGE131
(32)
here, the
Figure 886340DEST_PATH_IMAGE132
(omit)
Figure 22923DEST_PATH_IMAGE133
）;
Obtaining a set of closed solutions, will
Figure 895064DEST_PATH_IMAGE134
Observed values seen as noise
Figure 357269DEST_PATH_IMAGE135
Of a type using
Figure 314861DEST_PATH_IMAGE136
Representing the error phasor, the error of each element is respectively
Figure 786294DEST_PATH_IMAGE137
（
Figure 778521DEST_PATH_IMAGE138
) Then assume e is independently distributed over
Figure 910162DEST_PATH_IMAGE139
Mean 0 and variance
Figure 620629DEST_PATH_IMAGE140
From law of large numbers, for arbitrary
Figure 364594DEST_PATH_IMAGE141
All have:
Figure 211328DEST_PATH_IMAGE142
(33)
so that the following results are obtained:
Figure 280915DEST_PATH_IMAGE143
(34)
the transformed error phasor is
Figure 213099DEST_PATH_IMAGE144
，
Figure 26334DEST_PATH_IMAGE145
Figure 223179DEST_PATH_IMAGE146
The error phasor representing the error phasor for each element,m is the number of data points of the best match, and the variance of each group is not changed by transformation according to the orthogonal property of the matrix, thereby obtaining the data of each group
Figure 198089DEST_PATH_IMAGE147
Are all independently distributed in
Figure 883148DEST_PATH_IMAGE148
(zero mean and variance of
Figure 234495DEST_PATH_IMAGE149
) Most of them are theoretic, for any
Figure 55820DEST_PATH_IMAGE150
All have:
Figure 201631DEST_PATH_IMAGE151
(35)
the method comprises the following steps:
Figure 842828DEST_PATH_IMAGE152
(36)
to obtain:
Figure 496401DEST_PATH_IMAGE153
(37)
combining the two problems with the (26) problem to obtain:
Figure 172233DEST_PATH_IMAGE154
(38)
unknown quantity
Figure 754524DEST_PATH_IMAGE155
Is separable in the above formula, and each component is simplified as follows:
Figure 883017DEST_PATH_IMAGE156
(39)
wherein the content of the first and second substances,
Figure 576166DEST_PATH_IMAGE157
refers to the modulus of the corresponding phasor;
under the vector contraction form, the estimated value of x is expressed as:
Figure 372084DEST_PATH_IMAGE158
(40)
after the estimation value of the Taylor coefficient matrix x is obtained, the method comprises the following steps
Figure 125276DEST_PATH_IMAGE159
Harmonic sampling equation to harmonic signals
Figure 741065DEST_PATH_IMAGE160
In the framework of logistic regression, cross entropy was introduced [27 ]]The objective function measures the difference in probability distribution between the estimated and theoretical values x:
Figure 473791DEST_PATH_IMAGE161
(41)
assuming that the error is a binary distribution,
Figure 124216DEST_PATH_IMAGE162
the predicted probability distribution is considered to be closely related to the actual probability distribution, proving that the assumption is consistent with the expected model.
7. The harmonic detection method based on the sparse acquisition model as claimed in claim 1, wherein in the step c, the step of establishing a basic performance test scenario is as follows:
the basic performance of the proposed algorithm is tested, and a broadband dynamic signal model containing DDC components shown in formula (43) is constructed:
Figure 48309DEST_PATH_IMAGE163
(42)
wherein the content of the first and second substances,
Figure 151394DEST_PATH_IMAGE164
which is the fundamental frequency, here set at 50Hz,
Figure 920767DEST_PATH_IMAGE165
、
Figure 425698DEST_PATH_IMAGE166
respectively represent the fundamental wave and each harmonic phase angle, and are in
Figure 753649DEST_PATH_IMAGE167
The value of the harmonic frequency h of the low frequency band is 2-13, the harmonic frequency h of the high frequency band is 77, 79, 80 and 83 respectively, the sampling frequency is set to 10kHz, and the amplitude of DDC component is set to be any value in the range
Figure 344030DEST_PATH_IMAGE168
Is taken to be 0.3, 0.4, 0.5, … 1, the time constant of DDC component
Figure 917094DEST_PATH_IMAGE169
Starting from 0.01s, changing from 0.01s to 0.1s in steps;
using the BMP algorithm, the initial value of the iterative phasor is
Figure 542110DEST_PATH_IMAGE170
Regularization coefficient of
Figure 808007DEST_PATH_IMAGE171
Fixed value parameter
Figure 885684DEST_PATH_IMAGE172
，
Figure 262439DEST_PATH_IMAGE173
，
Figure 995425DEST_PATH_IMAGE174
When the zeta value is in the range of [0.05,0.3 ]]The number of data points in the search window that are the best matches is set to m.
8. The harmonic detection method based on the sparse acquisition model as claimed in claim 7, wherein in the step c, the step of establishing a frequency ramp test scenario is as follows:
the BMP was analyzed for performance under a frequency ramp, providing the following signals:
Figure 432223DEST_PATH_IMAGE175
(43)
wherein the content of the first and second substances,
Figure 997196DEST_PATH_IMAGE176
is the fundamental frequency and takes the value of 50Hz,
Figure 177642DEST_PATH_IMAGE177
the slope of fundamental frequency is 1Hz/s,
Figure 777251DEST_PATH_IMAGE178
and
Figure 119370DEST_PATH_IMAGE179
a fundamental phase and a harmonic phase, respectively, the phases being set to
Figure 906061DEST_PATH_IMAGE180
Random numbers uniformly distributed within the range.
9. The harmonic detection method based on the sparse acquisition model as claimed in claim 8, wherein in step c, the step of establishing the step change test scenario is as follows:
at the start of the test, the amplitude of each component was set to 115% of the initial amplitude, and the phase was changed to
Figure 123153DEST_PATH_IMAGE181
The broadband dynamic signal provided is as follows:
Figure 577268DEST_PATH_IMAGE182
(44)
the sampling rate is set to be 5kHz, the length of the sampling period is 5 periods, and in the standard, the performance of each algorithm under the condition of step change is evaluated by the speed of response time.
10. The harmonic detection method based on the sparse acquisition model as claimed in claim 9, wherein in the step c, the step of establishing an interference test scenario is as follows:
gaussian white noise with the signal-to-noise ratio of 55dB is introduced into the signal, and specific broadband dynamic signals are as follows:
Figure 90289DEST_PATH_IMAGE183
(45)
wherein, the first and the second end of the pipe are connected with each other,
Figure 364276DEST_PATH_IMAGE184
is the inter-harmonic frequency, and the values thereof are 9652.5Hz, 9751.5Hz, 9850.5Hz, 9949.5Hz, 10048.5Hz and 10147.5Hz in sequence,
Figure 886524DEST_PATH_IMAGE185
and
Figure 929566DEST_PATH_IMAGE186
is fundamental wave and harmonic waveAnd takes on a value of
Figure 879068DEST_PATH_IMAGE187
The random number in the range is set to have a sampling frequency of 10kHz and a sampling window length of 5 power frequency periods.
CN202210698874.3A 2022-06-20 2022-06-20 Harmonic detection method based on sparse acquisition model Pending CN114781196A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210698874.3A CN114781196A (en) 2022-06-20 2022-06-20 Harmonic detection method based on sparse acquisition model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210698874.3A CN114781196A (en) 2022-06-20 2022-06-20 Harmonic detection method based on sparse acquisition model

Publications (1)

Publication Number Publication Date
CN114781196A true CN114781196A (en) 2022-07-22

Family

ID=82421455

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210698874.3A Pending CN114781196A (en) 2022-06-20 2022-06-20 Harmonic detection method based on sparse acquisition model

Country Status (1)

Country Link
CN (1) CN114781196A (en)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115032453A (en) * 2022-08-08 2022-09-09 四川大学 Multi-frequency dynamic phasor measurement method
CN115343532A (en) * 2022-08-09 2022-11-15 国网福建省电力有限公司 Method for compressing and reconstructing disturbance-containing power quality signal based on compressed sensing
CN115856767A (en) * 2023-02-08 2023-03-28 安徽大学 Reconfigurable intelligent super-surface-assisted wave arrival direction estimation method
CN116150594A (en) * 2023-04-18 2023-05-23 长鹰恒容电磁科技(成都)有限公司 Method for identifying switch element characteristics in spectrum test data
CN117200242A (en) * 2023-11-08 2023-12-08 西安米特电子科技有限公司 Monitoring data processing method and system for intelligent voltage regulating cabinet
CN117639876A (en) * 2024-01-25 2024-03-01 南京理工大学 Linear frequency modulation wave anti-interference DOA estimation method based on space-time modulation super surface

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115032453A (en) * 2022-08-08 2022-09-09 四川大学 Multi-frequency dynamic phasor measurement method
CN115343532A (en) * 2022-08-09 2022-11-15 国网福建省电力有限公司 Method for compressing and reconstructing disturbance-containing power quality signal based on compressed sensing
CN115856767A (en) * 2023-02-08 2023-03-28 安徽大学 Reconfigurable intelligent super-surface-assisted wave arrival direction estimation method
CN116150594A (en) * 2023-04-18 2023-05-23 长鹰恒容电磁科技(成都)有限公司 Method for identifying switch element characteristics in spectrum test data
CN116150594B (en) * 2023-04-18 2023-07-07 长鹰恒容电磁科技(成都)有限公司 Method for identifying switch element characteristics in spectrum test data
CN117200242A (en) * 2023-11-08 2023-12-08 西安米特电子科技有限公司 Monitoring data processing method and system for intelligent voltage regulating cabinet
CN117200242B (en) * 2023-11-08 2024-02-02 西安米特电子科技有限公司 Monitoring data processing method and system for intelligent voltage regulating cabinet
CN117639876A (en) * 2024-01-25 2024-03-01 南京理工大学 Linear frequency modulation wave anti-interference DOA estimation method based on space-time modulation super surface
CN117639876B (en) * 2024-01-25 2024-04-23 南京理工大学 Linear frequency modulation wave anti-interference DOA estimation method based on space-time modulation super surface

Similar Documents

Publication Publication Date Title
CN114781196A (en) Harmonic detection method based on sparse acquisition model
CN108037361B (en) High-precision harmonic parameter estimation method based on sliding window DFT
Boashash et al. Use of the cross Wigner-Ville distribution for estimation of instantaneous frequency
Messina et al. Nonstationary approaches to trend identification and denoising of measured power system oscillations
Yu A multisynchrosqueezing-based high-resolution time-frequency analysis tool for the analysis of non-stationary signals
CN107102255B (en) Single ADC acquisition channel dynamic characteristic test method
CN111289796B (en) Method for detecting subsynchronous oscillation of high-proportion renewable energy power system
Pan et al. A noise reduction method of symplectic singular mode decomposition based on Lagrange multiplier
CN110311392B (en) Electric power system oscillation mode and mode identification method based on SPDMD
CN109659957B (en) APIT-MEMD-based power system low-frequency oscillation mode identification method
Karpilow et al. Characterization of non-stationary signals in electric grids: A functional dictionary approach
Philip et al. An improved Stochastic Subspace Identification based estimation of low frequency modes in power system using synchrophasors
JP2014153354A (en) Method for estimating frequencies and phases in three phase power system
Mei et al. Wavelet packet transform and improved complete ensemble empirical mode decomposition with adaptive noise based power quality disturbance detection
CN112816779B (en) Harmonic real signal parameter estimation method for analytic signal generation
CN112883318A (en) Multi-frequency attenuation signal parameter estimation algorithm of subtraction strategy
Bertocco et al. Frequency tracking for efficient phasor measurement based on a CSTFM model
CN115032453A (en) Multi-frequency dynamic phasor measurement method
CN109254202B (en) Synchronous phasor measurement device applied to power distribution network
CN114184838A (en) Power system harmonic detection method, system and medium based on SN mutual convolution window
Hu et al. Aerodynamic loads correction method based on wavelet packet transform in high-frequency force balance tests
Xingang et al. Supraharmonics measurement algorithm based on CS-SAMP
CN111398735A (en) Transformer substation grounding grid fault detection method based on information entropy
Frigo et al. Compressive sensing Taylor-Fourier and windowing approach for synchronized phasor, frequency and ROCOF measurements
CN116522269B (en) Fault diagnosis method based on Lp norm non-stationary signal sparse reconstruction

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20220722

RJ01 Rejection of invention patent application after publication