CN106300345A - Based on the low-frequency oscillation parameter identification method improving Prony algorithm - Google Patents

Abstract

The invention discloses a kind of based on the low-frequency oscillation parameter identification method improving Prony algorithm, neutral net is used to replace original matrix inversion operation, using complex matrix to be asked as weights, gradient descent method is utilized to adjust weights, solve the low-frequency oscillation parameter of Prony algorithm, it is to avoid complex matrix is asked repeatedly generalized inverse；Use loop iteration simultaneously; calculate simple; degree of accuracy is higher; by the correction of weights; verify cost function whether to accuracy rating within so that calculate more succinct, efficiently; operation for synchronous phasor measuring device brings optimization, and provides reliable initial data and data supporting for functions such as the monitoring of total system electrical network wide area, automation of transformation substations observing and controlling, stability contorting, selfadaptive computation.

Description

Based on the low-frequency oscillation parameter identification method improving Prony algorithm

Technical field

The present invention relates to a kind of low-frequency oscillation parameter identification method based on improvement Prony algorithm, belong to power technology neck Territory.

Background technology

Along with the expanding day of electrical network scale, large sized unit constantly putting into operation in electrical network, generally making of high-speed excitation With, low-frequency oscillation happens occasionally in large-scale interconnected power system.Power system is in oscillatory process, and transmission line of electricity power is back and forth Transmission, have impact on the properly functioning of power system, also directly reduces the transmission capacity of system, make power generation and transmittability Can not get maximum utilization, time more serious, system will lose synchronization.Therefore low-frequency oscillation is to threaten China's interconnected network safety and stability One of major issue run.

Prony algorithm owing to can directly estimate signal frequency, decay, amplitude and initial phase, therefore according to sampled value It is widely used in the identification of low-frequency oscillation of electric power system pattern.But low-frequency oscillation character identification at present is generally real at WAMS main website end Existing, the WAMS main website end PMU data to receiving uses the methods such as Prony algorithm such as to carry out low-frequency oscillation parameter identification.Along with The increase day by day of electrical network scale so that the burden of WAMS main website end increases the weight of day by day, therefore at synchronous phasor measuring device PMU (Phasor Measurment Unit) realizes low-frequency oscillation on-line parameter identification and has great importance.

At present, the application achievement that Prony algorithm obtains in power system, mainly it is reflected in low-frequency oscillation of electric power system inspection Successful Application in terms of survey, in the application prospect that the performance of the aspect such as frequency analysis, Parameter Estimation of Synchronous Machines is certain.Along with grinding That studies carefully deepens continuously and the development of technology, and Prony algorithm is paid close attention to the most widely by obtaining in power system.But at Prony During algorithm application, solving of parameter needs to use inverting of complex matrix, the dimension of the complex matrix that may relate in practice The highest, this brings trouble to calculating, even causes the situation without solving.Meanwhile, when the value of matrix determinant is the least, it is difficult to look for To its inverse matrix.At this time thinking that system is ill, ill-conditioned linear systems is difficult to accurately solve, so to avoid easily as far as possible There is inversion operation during morbid state.

Summary of the invention

In order to solve above-mentioned technical problem, the invention provides a kind of based on the low-frequency oscillation parameter improving Prony algorithm Discrimination method.

In order to achieve the above object, the technical solution adopted in the present invention is:

Low-frequency oscillation parameter identification method based on improvement Prony algorithm, comprises the following steps:

Step 1, collection sampled data y (1), y (2) ..., y (N)；

Wherein, N is sampling number, and y (n) is the n-th sampled value, n ∈ [1, N]；

Step 2, writes out the second moment sample matrix Re of Prony algorithm, determines second moment sample matrix according to sampled data Effective order p of Re and coefficient a₁,a₂,…,a_pLeast-squares estimation；

Step 3, solves proper polynomial 1+a₁z^-1+a₂z^-2+…+a_pz^-pCharacteristic root and sampled-data estimation value；

Wherein, z is Prony discrete function parameter.

Step 4, using the matrix Z that is made up of characteristic root as the sample matrix of neutral net, by by sampled-data estimation value structure The matrix becomeAs the estimation output matrix of neutral net, build complex matrix equationCalculate in complex matrix equation Parameter matrix, i.e. the weight matrix W of neutral net；

Step 5, utilizes parameter matrix to calculate the amplitude of low-frequency oscillation, phase place, frequency and damping ratio, it is achieved low-frequency oscillation mould Formula identification.

Second moment sample matrix Re is,

Wherein, pe is initial exponent number, and value is [N/2],i′∈ [1, pe], j ∈ [1, pe], i ' and j are integer, and y (n '-j) is the n-th '-j sampled value, and y (n '-i ') is that the n-th '-i ' is individual to be adopted Sample value.

The computing formula of sampled-data estimation value is,

$\hat{y}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\hat{y}(n-i)$

Wherein,Being the n-th sampled-data estimation value, i ∈ [1, p], i are integer.

Calculate the parameter matrix in complex matrix equation, i.e. the process of the weight matrix W of neutral net is,

S1, sets up Matrix division:

$\left\{\begin{array}{c}{Z}^{\′}{W}^{\′}-Z^{\′\′}{W}^{\′\′}={\hat{Y}}^{\′}\\ Z^{\′}{W}^{\′\′}+Z^{\′\′}{W}^{\′}=0\end{array}\right.$

Wherein, Z ', Z " being respectively real part and imaginary part, W ', the W of the row vector of Z " are respectively real part and the imaginary part of W,Point Not Wei estimate output real part and imaginary part, i.e.Real part and imaginary part；

S2, selects the initial value W of neutral net weight matrix₀, and resolve into real part and imaginary part；

S3, according to the Matrix division of S1, calculates real part Y ' and imaginary part Y of output "；

S4, by Y ', Y " compare with the real part and imaginary part estimating output respectively, obtain the error of real part and imaginary part arrange to Amount；

S5, squared to the error column vector of real part and imaginary part and, it is thus achieved that cost function；

Cost function is,

$\mathrm{\ξ}\left(\mathrm{\ω}\right)=\frac{1}{2}{E}^{\′}+\frac{1}{2}{E}^{\′\′}$

Wherein,

S6, it is judged that cost function whether within precision threshold, if, then current weight matrix is parameter matrix； If it was not then go to S7；

S7, seeks local derviation by cost function to each element in weight matrix, obtains local derviation matrix；

Local derviation matrix is,

$\left[\begin{array}{cc}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}& \frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}\end{array}\right]$

Wherein, W_i' for the real part of i-th weights, W in weight matrix_i" for the imaginary part of i-th weights in weight matrix；

S8, is adjusted weight matrix according to local derviation matrix, obtains new weight matrix；

In new weight matrix, the real part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′}=W_i^{\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}$

In new weight matrix, the imaginary part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′\′}=W_i^{\′\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}$

Wherein, η is learning rate；

S8, resolves into real part and imaginary part by new weight matrix, goes to step S3.

The computing formula of the amplitude of low-frequency oscillation, phase place, frequency and damping ratio is,

A_i=| W_i|

${\mathrm{\θ}}_i=\arctan \[\frac{\mathrm{Im}\left(W_i\right)}{\mathrm{Re}\left(W_i\right)}\]$

${\mathrm{\σ}}_i=\frac{ln\left|z_i\right|}{\mathrm{\Δ}t}$

$f_i=\frac{\arctan \[\frac{\mathrm{Im}\left(z_i\right)}{\mathrm{Re}\left(z_i\right)}\]}{2\mathrm{\π}\mathrm{\Δ}t}$

Wherein, A_iFor amplitude, θ_iFor phase place, σ_iFor damping ratio, f_iFor frequency, z_iFor ith feature root, when Δ t is for sampling Between be spaced.

The beneficial effect that the present invention is reached: the present invention use neutral net to replace original matrix inversion operation, Complex matrix to be asked, as weights, utilizes gradient descent method to adjust weights, solves the low-frequency oscillation parameter of Prony algorithm, it is to avoid Complex matrix is asked repeatedly generalized inverse；The present invention is mainly loop iteration simultaneously, therefore calculates simple, and degree of accuracy is higher, passes through The correction of weights, verify cost function whether to accuracy rating within so that calculate more succinct, efficiently, for synchronizing The operation of phasor measuring set brings optimization, and be the monitoring of total system electrical network wide area, automation of transformation substations observing and controlling, stability contorting, The functions such as selfadaptive computation provide reliable initial data and data supporting.

Accompanying drawing explanation

Fig. 1 is the flow chart of the present invention.

Fig. 2 is neural network structure figure.

Fig. 3 is the Matlab analogous diagram of original signal and Prony fitted signal.

Detailed description of the invention

The invention will be further described below in conjunction with the accompanying drawings.Following example are only used for clearly illustrating the present invention Technical scheme, and can not limit the scope of the invention with this.

As it is shown in figure 1, low-frequency oscillation parameter identification method based on improvement Prony algorithm, comprise the following steps:

Step 1, collection sampled data y (1), y (2) ..., y (N)；Wherein, N is sampling number, and y (n) is the n-th sampling Value, n ∈ [1, N].

Step 2, writes out the second moment sample matrix Re of Prony algorithm, determines second moment sample matrix according to sampled data Effective order p of Re and coefficient a₁,a₂,…,a_pLeast-squares estimation.

Second moment sample matrix Re is,

Wherein, pe is initial exponent number, and value is [N/2],i′∈ [1, pe], j ∈ [1, pe], i ' and j are integer, and y (n '-j) is the n-th '-j sampled value, and y (n '-i ') is that the n-th '-i ' is individual to be adopted Sample value.

Step 3, solves proper polynomial 1+a₁z^-1+a₂z^-2+…+a_pz^-pCharacteristic root and sampled-data estimation value；Wherein, z is Prony discrete function parameter.

The computing formula of sampled-data estimation value is,

$\hat{y}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\hat{y}(n-i)$

Wherein,Being the n-th sampled-data estimation value, i ∈ [1, p], i are integer.

Step 4, as in figure 2 it is shown, using the matrix Z that is made up of characteristic root as the sample matrix of neutral net, by by adopt The matrix that sample estimated value is constitutedAs the estimation output matrix of neutral net, build complex matrix equationCalculate multiple The weight matrix W of the parameter matrix in matrix equation, i.e. neutral net.

Detailed process is:

S1, sets up Matrix division:

$\left\{\begin{array}{c}{Z}^{\′}{W}^{\′}-Z^{\′\′}{W}^{\′\′}={\hat{Y}}^{\′}\\ Z^{\′}{W}^{\′\′}+Z^{\′\′}{W}^{\′}=0\end{array}\right.$

Wherein, Z ', Z " being respectively real part and imaginary part, W ', the W of the row vector of Z " are respectively real part and the imaginary part of W,Point Not Wei estimate output real part and imaginary part, i.e.Real part and imaginary part；

S2, selects the initial value W of neutral net weight matrix₀, and resolve into real part and imaginary part；

S3, according to the Matrix division of S1, calculates real part Y ' and imaginary part Y of output "；

S4, by Y ', Y " compare with the real part and imaginary part estimating output respectively, obtain the error of real part and imaginary part arrange to Amount；

S5, squared to the error column vector of real part and imaginary part and, it is thus achieved that cost function；

Cost function is,

$\mathrm{\ξ}\left(\mathrm{\ω}\right)=\frac{1}{2}{E}^{\′}+\frac{1}{2}{E}^{\′\′}$

Wherein,

S6, it is judged that cost function whether within precision threshold, if, then current weight matrix is parameter matrix； If it was not then go to S7；

S7, seeks local derviation by cost function to each element in weight matrix, obtains local derviation matrix；

Local derviation matrix is,

$\left[\begin{array}{cc}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}& \frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}\end{array}\right]$

Wherein, W_i' for the real part of i-th weights, W in weight matrix_i" for the imaginary part of i-th weights in weight matrix；

S8, is adjusted weight matrix according to local derviation matrix, obtains new weight matrix；

In new weight matrix, the real part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′}=W_i^{\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}$

In new weight matrix, the imaginary part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′\′}=W_i^{\′\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}$

Wherein, η is learning rate；

S8, resolves into real part and imaginary part by new weight matrix, goes to step S3.

Step 5, utilizes parameter matrix to calculate the amplitude of low-frequency oscillation, phase place, frequency and damping ratio, it is achieved low-frequency oscillation mould Formula identification.

Specific formula for calculation is,

A_i=| W_i|

${\mathrm{\θ}}_i=\arctan \[\frac{\mathrm{Im}\left(W_i\right)}{\mathrm{Re}\left(W_i\right)}\]$

${\mathrm{\σ}}_i=\frac{ln\left|z_i\right|}{\mathrm{\Δ}t}$

$f_i=\frac{\arctan \[\frac{\mathrm{Im}\left(z_i\right)}{\mathrm{Re}\left(z_i\right)}\]}{2\mathrm{\π}\mathrm{\Δ}t}$

Wherein, A_iFor amplitude, θ_iFor phase place, σ_iFor damping ratio, f_iFor frequency, z_iFor ith feature root, when Δ t is for sampling Between be spaced.

In order to further illustrate said method, carry out comparative testing below.

Emulation signal y (t) is made up of three kinds of oscillationg components, and expression formula is as follows:

$\begin{array}{c}y\left(t\right)=0.8e^{-0.1t}\·\cos \left(2\mathrm{\π}\·0.5t+\mathrm{\π}/3\right)\\ +0.4e^{-t}\·\cos \left(2\mathrm{\π}\·1.0t+\mathrm{\π}/4\right)\\ +1.0e^{-0.5t}\·\sin \left(2\mathrm{\π}\·1.5t\right)\end{array}$

Sampled point is 64 points, and employing frequency is 10HZ, and time series 6.4 seconds, learning rate value 0.05, weight matrix is initial Value is null matrix.

In MATLAB, employing the inventive method is to complex matrix equation solution, and compares result, concrete such as table one With two.

Table one is accurately to solve

1	0.0000-0.5000i
		2	0.0000+0.5000i
3	0.2000+0.3464i
		4	0.2000-0.3464i
5	0.1414+0.1414i
		6	0.1414-0.1414i

Table two is the solving result of the present invention

1	0.0000-0.5000i
		2	0.0000+0.5000i
3	0.2004+0.3462i
		4	0.2004-0.3462i
5	0.1414+0.1414i
		6	0.1414-0.1414i

The resultant error of the solution obtained by the present invention is the least as can be seen from the table, and algorithm is simple, it is easy to accomplish, phase Ratio solves generalized inverse, more practical.

As it is shown on figure 3, be the Matlab analogous diagram of original signal and Prony fitted signal, Prony algorithm is used to combine this Bright still can effectively extract oscillating component.

In sum, said method uses neutral net to replace original matrix inversion operation, and complex matrix to be asked is made For weights, utilize gradient descent method to adjust weights, solve the low-frequency oscillation parameter of Prony algorithm, it is to avoid be anti-to complex matrix Ask generalized inverse again；Using loop iteration simultaneously, calculate simple, degree of accuracy is higher, by the correction of weights, verifies cost function Whether to accuracy rating within so that calculating more succinct, efficiently, the operation for synchronous phasor measuring device brings excellent Change, and provide for functions such as the monitoring of total system electrical network wide area, automation of transformation substations observing and controlling, stability contorting, selfadaptive computation Reliable initial data and data supporting.

The above is only the preferred embodiment of the present invention, it is noted that for the ordinary skill people of the art For Yuan, on the premise of without departing from the technology of the present invention principle, it is also possible to make some improvement and deformation, these improve and deformation Also should be regarded as protection scope of the present invention.

Claims (5)

1. based on the low-frequency oscillation parameter identification method improving Prony algorithm, it is characterised in that: comprise the following steps:

Step 1, collection sampled data y (1), y (2) ..., y (N)；

Wherein, N is sampling number, and y (n) is the n-th sampled value, n ∈ [1, N]；

Step 2, writes out the second moment sample matrix Re of Prony algorithm, determines second moment sample matrix Re's according to sampled data Effectively order p and coefficient a₁,a₂,…,a_pLeast-squares estimation；

Step 3, solves proper polynomial 1+a₁z^-1+a₂z^-2+…+a_pz^-pCharacteristic root and sampled-data estimation value；

Wherein, z is Prony discrete function parameter；

Step 4, using the matrix Z that is made up of characteristic root as the sample matrix of neutral net, by by sampled-data estimation value constitute MatrixAs the estimation output matrix of neutral net, build complex matrix equationCalculate the ginseng in complex matrix equation The weight matrix W of matrix number, i.e. neutral net；

Step 5, utilizes parameter matrix to calculate the amplitude of low-frequency oscillation, phase place, frequency and damping ratio, it is achieved low frequency oscillation mode is distinguished Know.

It is the most according to claim 1 based on the low-frequency oscillation parameter identification method improving Prony algorithm, it is characterised in that: Second moment sample matrix Re is,

Wherein, pe is initial exponent number, and value is [N/2],I ' ∈ [1, pe], j ∈ [1, pe], i ' and j are integer, and y (n '-j) is the n-th '-j sampled value, and y (n '-i ') is the individual sampled value of the n-th '-i '.

It is the most according to claim 1 based on the low-frequency oscillation parameter identification method improving Prony algorithm, it is characterised in that: The computing formula of sampled-data estimation value is,

$\hat{y}\left(n\right)=-\underset{i=1}{\overset{p}{\Σ}}{a}_i\hat{y}(n-i)$

Wherein,Being the n-th sampled-data estimation value, i ∈ [1, p], i are integer.

It is the most according to claim 1 based on the low-frequency oscillation parameter identification method improving Prony algorithm, it is characterised in that: Calculate the parameter matrix in complex matrix equation, i.e. the process of the weight matrix W of neutral net is,

S1, sets up Matrix division:

$\left\{\begin{array}{c}{Z}^{\′}{W}^{\′}-Z^{\′\′}{W}^{\′\′}={\hat{Y}}^{\′}\\ Z^{\′}{W}^{\′\′}+Z^{\′\′}{W}^{\′}=0\end{array}\right.$

Wherein, Z ', Z " being respectively real part and imaginary part, W ', the W of the row vector of Z " are respectively real part and the imaginary part of W,0 is respectively Estimate real part and the imaginary part of output, i.e.Real part and imaginary part；

S2, selects the initial value W of neutral net weight matrix₀, and resolve into real part and imaginary part；

S3, according to the Matrix division of S1, calculates real part Y ' and imaginary part Y of output "；

S4, by Y ', Y " to be compared with the real part and imaginary part estimating output respectively, obtains the error column vector of real part and imaginary part；

S5, squared to the error column vector of real part and imaginary part and, it is thus achieved that cost function；

Cost function is,

$\mathrm{\ξ}\left(\mathrm{\ω}\right)=\frac{1}{2}{E}^{\′}+\frac{1}{2}{E}^{\′\′}$

Wherein,

S6, it is judged that cost function whether within precision threshold, if, then current weight matrix is parameter matrix；If Do not exist, then go to S7；

S7, seeks local derviation by cost function to each element in weight matrix, obtains local derviation matrix；

Local derviation matrix is,

$\left[\begin{array}{cc}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}& \frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}\end{array}\right]$

Wherein, W_i' for the real part of i-th weights, W in weight matrix_i" for the imaginary part of i-th weights in weight matrix；

S8, is adjusted weight matrix according to local derviation matrix, obtains new weight matrix；

In new weight matrix, the real part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′}=W_i^{\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′}}$

In new weight matrix, the imaginary part of i-th weights is,

${\stackrel{\‾}{W}}_i^{\′\′}=W_i^{\′\′}-\mathrm{\η}\frac{\∂\mathrm{\ξ}\left(\mathrm{\ω}\right)}{\∂W_i^{\′\′}}$

Wherein, η is learning rate；

S8, resolves into real part and imaginary part by new weight matrix, goes to step S3.

It is the most according to claim 1 based on the low-frequency oscillation parameter identification method improving Prony algorithm, it is characterised in that: The computing formula of the amplitude of low-frequency oscillation, phase place, frequency and damping ratio is,

A_i=| W_i|

${\mathrm{\θ}}_i=\arctan \[\frac{\mathrm{Im}\left(W_i\right)}{\mathrm{Re}\left(W_i\right)}\]$

${\mathrm{\σ}}_i=\frac{ln\left|z_i\right|}{\mathrm{\Δ}t}$

$f_i=\frac{\arctan \[\frac{\mathrm{Im}\left(z_i\right)}{\mathrm{Re}\left(z_i\right)}\]}{2\mathrm{\π}\mathrm{\Δ}t}$

Wherein, A_iFor amplitude, θ_iFor phase place, σ_iFor damping ratio, f_iFor frequency, z_iFor ith feature root, Δ t is between the sampling time Every.