CN104852392A - Calculation method of sub-synchronous oscillation mode attenuation coefficients based on Prony algorithm - Google Patents

Abstract

The invention relates to a calculation method of the sub-synchronous oscillation mode attenuation coefficients based on Prony algorithm, and mainly aims at solving the problem in extracting key modal parameters as modal frequency and attenuation coefficients from sub-synchronous signals. The calculation method is based on the Prony algorithm, takes shafting rotating speed signal of high sensitivity as input signals, employs least square fit, greatly eliminates influence of noise signals, requires no filtering, reserves all data information, is capable of extracting multiple modal information, and thus, extracts attenuation coefficients of different modes accurately and rapidly.

Description

A kind of sub-synchronous oscillation Modal Decay coefficient calculation method based on Prony algorithm

Technical field

The present invention relates to the Subsynchronous Oscillation of Turbine-generator Set Modal Decay coefficient calculation method be applied in subsynchronous oscillation of electrical power system analysis and simulation.Be specifically related to based on Prony algorithm, calculate sub-synchronous oscillation Modal Decay coefficient, eliminate the noise effect in measuring process, avoid the impact of conventional filter, remain total data information.

Background technology

For carrying out remote, large capacity transmission, the application of high voltage direct current transmission and circuit series capacitor compensation is more and more general, and the sub-synchronous oscillation problem caused thus also becomes increasingly conspicuous.Subsynchronous oscillation of electrical power system is a kind of power oscillation lower than power frequency, and it can cause the fatigue loss of generating set macro-axis, can cause the fracture of turbonator shafting when serious, seriously threatens the safe operation of electric power system.

When sub-synchronous oscillation produces, the mechanical part in electric power system and intercoupling of electric part occurrence dynamics act on and exchange energy.Torsional oscillation mode damping is when vibrating under axle ties up to torsion frequency, and the quantization means of its torsional oscillation rate of decay, comprises mechanical damping system and electrical system damping.When system presents underdamping or negative resistance character, axle system rotor persistent oscillation even oscillation and divergence can be caused, cause shaft train instability, even axle system fracture time serious.Therefore axle system modal damping is assessed very necessary.

Axle system modal damping is one of three elements of sub-synchronous oscillation, and the rate of decay after sign sub-synchronous oscillation is excited is the important criterion judging that whether current system is stable.The calculating of modal damping, is adopted the method for theory calculate comparatively difficult, is generally obtained by field test.Modal damping is relevant to a lot of system parameters, and theory calculate needs a large amount of parameters, and carries out detailed modeling to system, and workload is large, easily makes mistakes.Field test needs the sub-synchronous oscillation of activating system under each subsynchronous frequency, and quantities is large, and adds unit operation risk.

Summary of the invention

After the present invention is mainly used in unit generation sub-synchronous oscillation, the identification problem of mode tach signal rate of decay.Original tach signal after occurring using sub-synchronous oscillation, as input data, adopts Prony algorithm, extracts each Modal Decay coefficient information quickly and accurately.

Usual turbonator shafting has the torsional oscillation mode of multiple frequency, and the generator shaft system tach signal that transducer records is the superposition of multiple mode signals, is shown below:

$x\left(t\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k{e}^{{\mathrm{\σ}}_{k}t}\cos ({\mathrm{\ω}}_{k}t+{\mathrm{\θ}}_k)$

In formula, x is generator shaft system tach signal, and m is mode number, A _kfor the amplitude of a kth mode, σ _kfor the attenuation coefficient of a kth mode, ω _kfor the angular frequency of a kth mode, θ _kfor the phase place of a kth mode, t is sampling time sequence.

The present invention is based on Prony algorithm, from generator shaft system tach signal x, extract the amplitude A of m mode, frequencies omega, attenuation coefficient σ and phase theta information.Concrete steps comprise:

(1) analytic signal is determined: for the transmitting system adopting high voltage direct current transmission or string to mend transmission of electricity, after the disturbances such as grid side generation single phase ground fault, take out the data containing N number of sampling instant in generator shaft system tach signal x, x, be designated as x (1) respectively, x (2) ..., x (g), x (N), the x signal sampling time interval is dt, then the time series t=gdt that g data in x are corresponding, namely $x\left(g\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k{e}^{{\mathrm{\σ}}_k\mathrm{gdt}}\cos ({\mathrm{\ω}}_k\mathrm{gdt}+{\mathrm{\θ}}_k);$

(2) the mode quantity analyzed is determined: set the mode number that will analyze as m, each model frequency is designated as f (1), f (2) ... f (m).

(3) use analytic signal structural matrix, solve least square solution.Structural matrix X, X are the matrix of (N-3m) × (3m+1), by x data allocations in matrix X, and 1st ~ 3m+1 the data of the first behavior x, 2nd ~ 3m+2 the data of the second behavior x, the like, until X matrix is booked.The 1st row getting matrix respectively give vectorial Xb, 2nd ~ 3m+1 row imparting matrix Xa, adopt least square method to solve its solution is designated as a _i;

$X=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& ...& x(3m+1)\\ x\left(2\right)& x\left(3\right)& ...& x(3m+2)\\ x\left(3\right)& x\left(4\right)& ...& x(3m+3)\\ ...& ...& ...& ...\\ x(N-3m)& x(N-3m+1)& ...& x\left(N\right)\end{array}\right]$

$\mathrm{Xa}=\left[\begin{array}{cccc}x\left(2\right)& x\left(3\right)& \·\·\·& x(3m+1)\\ x\left(3\right)& x\left(4\right)& \·\·\·& x(3m+2)\\ x\left(4\right)& x\left(5\right)& \·\·\·& x(3m+3)\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ x(N-3m+1)& x(N-3m+2)& \·\·\·& x\left(N\right)\end{array}\right]$

$\mathrm{Xb}=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ \·\·\·\\ x(N-3m)\end{array}\right];$

(4) structure is with 1, a _i(1), a _i(2) ... a _i(3m) be the equation of coefficient

u ^3m+a _i(1)u ^3m-1+a _i(2)u ^3m-2+…………+a _i(3m-1)u+a _i(3m)＝0

Solve this equation, obtain 3m root u (i), wherein i=1 ~ 3m;

(5) structural matrix A and B,

$A=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ \mathrm{ai}\left(1\right)& \mathrm{ai}\left(2\right)& \·\·\·& \mathrm{ai}\left(3m\right)\\ \mathrm{ai}{\left(1\right)}^2& \mathrm{ai}{\left(2\right)}^2& \·\·\·& \mathrm{ai}{\left(3m\right)}^2\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ \mathrm{ai}{\left(1\right)}^{N-1}& \mathrm{ai}{\left(2\right)}^{N-1}& \·\·\·& \mathrm{ai}{\left(3m\right)}^{N-1}\end{array}\right]B=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ \·\·\·\·\·\·\\ x\left(N\right)\end{array}\right]$

Employing least square method solves its solution is designated as b;

(6) calculating of model frequency f:

$f\left(i\right)=\frac{\arctan \left(\frac{\mathrm{Im}\left(u\left(i\right)\right)}{\mathrm{Re}\left(u\left(i\right)\right)}\right)}{2\mathrm{\πdt}},i=1~3m$

The calculating of Modal Decay factor sigma:

${\mathrm{\σ}}_i=\frac{\ln \left|u\left(i\right)\right|}{\mathrm{dt}},i=1~3m;$

(7) for m mode, the actual value having obtained 3m mode, get distance f (1), f (2), the Frequency point that f (m) is nearest, apart from when identical, get the point that attenuation rate absolute value is less, thus find the attenuation coefficient of corresponding model frequency in 3m point.

The present invention has following useful technique effect:

The present invention analyzes tach signal during unit operation, therefrom obtains Modal Decay coefficient.Adopt the data under unit actual operating mode to input as analysis, the Modal Decay coefficient results obtained thus is more accurate.

The present invention adopts Prony algorithm, using the axle system tach signal of equal interval sampling as analysis data, adopt least square fitting, improve noise resisting ability, extract the frequency of multiple mode, attenuation coefficient, amplitude and phase information without the need to filtering, calculate accurately and rapidly.

Prony method is a kind of time domain approach of Modal Parameter Identification, is to carry out matching equal interval sampling data with the linear combination of exponential function, therefrom can analyze the frequency of signal, decay factor, amplitude and phase place.Prony algorithm is the transform analysis in time domain, there is not the leakage problem in frequency-domain analysis.

Accompanying drawing explanation

Fig. 1 illustrates the sub-synchronous oscillation Modal Decay coefficient calculations flow process based on Prony algorithm;

Fig. 2 illustrates the axle system tach signal to certain unit, the fitted signal obtained according to the present invention and the contrast of original input signal;

Curve line is truly recorded at the scene that Fig. 3 illustrates certain unit, is followed successively by mode one sampled value from top to bottom, mode two sampled value, mode three sampled value, tach signal sampled value;

Fig. 4 illustrates the axle system tach signal to field measurement, the fitted signal obtained according to the present invention and the contrast of original input signal.

Embodiment

According to the present invention, as shown in Figure 1, the step of calculating generator axle system tach signal Modal Decay coefficient is as follows.

(1) analytic signal is determined: for the transmitting system adopting high voltage direct current transmission or string to mend transmission of electricity, after the disturbance of grid side generation single phase ground fault, take out generator shaft system tach signal x, the data containing N number of sampling instant in x, be designated as x (1) respectively, x (2) ..., x (g), x (N), the x signal sampling time interval is dt, then the time series t=gdt that g data in x are corresponding, namely

(2) the mode quantity analyzed is determined: set the mode number that will analyze as m, each model frequency is designated as f (1), f (2) ... f (m);

(3) analytic signal structural matrix is used, solve least square solution: structural matrix X, X is the matrix of (N-3m) × (3m+1), by x data allocations in matrix X, and 1st ~ 3m+1 the data of the first behavior x, 2nd ~ 3m+2 the data of the second behavior x, the like, until be booked by X matrix, the 1st row getting matrix X give vectorial Xb, 2nd ~ 3m+1 the row getting matrix X form matrix Xa, adopt least square method to solve its solution is designated as a _i;

$X=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& ...& x(3m+1)\\ x\left(2\right)& x\left(3\right)& ...& x(3m+2)\\ x\left(3\right)& x\left(4\right)& ...& x(3m+3)\\ ...& ...& ...& ...\\ x(N-3m)& x(N-3m+1)& ...& x\left(N\right)\end{array}\right]$

$\mathrm{Xa}=\left[\begin{array}{cccc}x\left(2\right)& x\left(3\right)& \·\·\·& x(3m+1)\\ x\left(3\right)& x\left(4\right)& \·\·\·& x(3m+2)\\ x\left(4\right)& x\left(5\right)& \·\·\·& x(3m+3)\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ x(N-3m+1)& x(N-3m+2)& \·\·\·& x\left(N\right)\end{array}\right]$

$\mathrm{Xb}=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ \·\·\·\\ x(N-3m)\end{array}\right];$

(4) structure is with 1, ai (1), ai (2) ... the equation that ai (3m) is coefficient

U ^3m+ a _i(1) u ^3m-1+ a _i(2) u ^3m-2+ ... + a _i(3m-1) u+a _i(3m)=0 solves this equation, obtains 3m root u (i), wherein i=1 ~ 3m;

(5) structural matrix A and B,

$A=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ \mathrm{ai}\left(1\right)& \mathrm{ai}\left(2\right)& \·\·\·& \mathrm{ai}\left(3m\right)\\ \mathrm{ai}{\left(1\right)}^2& \mathrm{ai}{\left(2\right)}^2& \·\·\·& \mathrm{ai}{\left(3m\right)}^2\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ \mathrm{ai}{\left(1\right)}^{N-1}& \mathrm{ai}{\left(2\right)}^{N-1}& \·\·\·& \mathrm{ai}{\left(3m\right)}^{N-1}\end{array}\right]B=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ \·\·\·\·\·\·\\ x\left(N\right)\end{array}\right]$

Employing least square method solves its solution is designated as b;

(6) calculating of model frequency f:

$f\left(i\right)=\frac{\arctan \left(\frac{\mathrm{Im}\left(u\left(i\right)\right)}{\mathrm{Re}\left(u\left(i\right)\right)}\right)}{2\mathrm{\πdt}},i=1~3m$

The calculating of Modal Decay factor sigma:

${\mathrm{\σ}}_i=\frac{\ln \left|u\left(i\right)\right|}{\mathrm{dt}},i=1~3m;$

(7) for m mode, the actual value having obtained 3m mode, get distance f (1), f (2), the Frequency point that f (m) is nearest, apart from when identical, get the point that attenuation rate absolute value is less, thus find the attenuation coefficient of corresponding model frequency in 3m point.

For certain turbo generator set, three model frequencies are respectively 13.38Hz, 23.77Hz, 26.51Hz.Generator shaft system tach signal is the superposition of three model frequency signals, and is mixed with high-frequency noise, is shown below: x (t)=5 × cos (2 π × 13.38t) × e ^-0.8t+ 2 × cos (2 π × 23.77t) × e ^-0.5t+ 10 × cos (2 π × 26.51t) × e ^-0.2t+ 0.001 × cos (2 π × 300t)

Sampling interval is 0.001s, and sampling duration is 60 seconds.

Frequency number is 3, then p=3 × 3=9, sampling interval dt=0.001, brings x, p, dt tri-parameters into computational methods of the present invention, obtains following result:

Certain generator unit shaft system tach signal of table 1 adopts the result of calculation of Prony algorithm

Sequence number	Attenuation coefficient	Model frequency
			First group	-74.18936	0
Second group	0	-200
			3rd group	0	200
4th group	-0.20003	26.50999
			5th group	-0.20003	-26.50999
6th group	-0.49996	23.77001
			7th group	-0.49996	-23.77001
8th group	-0.80001	13.37999
			9th group	-0.80001	-13.37999

According to selection principle, the attenuation coefficient absolute value of first group is bigger than normal, can directly cast out.The model frequency of the 4th group, the 6th group and the 8th group is close to actual frequency, and the attenuation coefficient of its correspondence is Trueattenuation coefficient, and through comparing, this value is consistent with original input value, calculates correct.The fitted signal of the solution utilizing the inventive method to obtain compares as Fig. 2 with original input signal, and two curves are basically identical.From above result, the inventive method can obtain frequency and the attenuation coefficient of sub-synchronous oscillation signal quickly and accurately, and effectively can remove noise jamming.

Adopt compute mode attenuation coefficient of the present invention with the performance comparison adopting each mode curve calculation Modal Decay coefficient in prior art:

Certain unit three model frequencies are respectively 13.38Hz, 23.77Hz, 26.51Hz.There is trip accident in certain secondary net side, unit TSR startup separator record ripple, and record curve line as shown in Figure 3.First three curve is mode curve, and the last item curve is original speed curves.Respectively the Fitting Calculation is carried out to first three curve, obtains the attenuation coefficient of each mode as following table:

Table 2 mode curve calculation result

Model frequency	Mode 1	Mode 2	Mode 3
				Modal Decay coefficient	0.2241	0.2827	0.2834

With original tach signal for input data, use the inventive method, calculate each Modal Decay coefficient.Computational process is as follows:

(1) judge according to each mode curve in Fig. 3, the time period of each mode exponentially attenuation trend is 1.3s ~ 10s, and the data in original tach signal in this time period are as analysis data x.

(2) sample frequency is 1000Hz, then sampling interval dt is 0.001s.

(3) in theory, the mode number m of unit is 3, then p=3m=9, brings x, dt, p into computational methods of the present invention, the calculated value substantial deviation actual value obtained, and the curve obtained by calculating value compares with primitive curve and has a long way to go.Analyze reason, site environment is complicated, and noise signal is more, the frequency number contained in tach signal is far above 3, the size of adjustment p, makes to include three theoretical model frequencies in the model frequency calculated, simultaneously matched curve and primitive curve comparatively close.Through adjustment, when p is 150, result of calculation meets the requirements.Adopt the contrast of the result of result of calculation of the present invention and employing mode curve calculation as following table, the contrast of matched curve and original input signal as shown in Figure 4.

Table 3 comparison of computational results

By upper table and Fig. 4, can find out and use the inventive method can obtain each model frequency and Modal Decay coefficient exactly with speed faster.

Claims (1)

1. based on a sub-synchronous oscillation Modal Decay coefficient calculation method for Prony algorithm, turbonator shafting has the torsional oscillation mode of multiple frequency, and generator shaft system tach signal x is the superposition of multiple mode signals, and its expression formula is wherein, m is mode number, A _kfor the amplitude of a kth mode, σ _kfor the attenuation coefficient of a kth mode, ω _kfor the angular frequency of a kth mode, θ _kfor the phase place of a kth mode, t is sampling time sequence, it is characterized in that: described computational methods comprise the following steps:

(1) analytic signal is determined: for the transmitting system adopting high voltage direct current transmission or string to mend transmission of electricity, after the disturbance of grid side generation single phase ground fault, take out generator shaft system tach signal x, the data containing N number of sampling instant in x, be designated as x (1) respectively, x (2) ..., x (g), x (N), the x signal sampling time interval is dt, then the time series t=gdt that g data in x are corresponding, namely $x\left(g\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{A}_k{e}^{{\mathrm{\σ}}_k\mathrm{gdt}}\cos ({\mathrm{\ω}}_k\mathrm{gdt}+{\mathrm{\θ}}_k);$

(2) the mode quantity analyzed is determined: set the mode number that will analyze as m, each model frequency is designated as f (1), f (2) ... f (m).

(3) analytic signal structural matrix is used, solve least square solution: structural matrix X, X is the matrix of (N-3m) × (3m+1), by x data allocations in matrix X, and 1st ~ 3m+1 the data of the first behavior x, 2nd ~ 3m+2 the data of the second behavior x, the like, until be booked by X matrix, the 1st row getting matrix X give vectorial Xb, 2nd ~ 3m+1 the row getting matrix X form matrix Xa, adopt least square method to solve its solution is designated as a _i;

$X=\left[\begin{array}{cccccc}& x\left(1\right)& & x\left(2\right)& ...& x(3m+1)\\ & x\left(2\right)& & x\left(3\right)& ...& x(3m+2)\\ & x\left(3\right)& & x\left(4\right)& ...& x(3m+3)\\ ...& & ...& & ...& ...\\ & x(N-3m)& & x(N-3m+1)& ...& x\left(N\right)\end{array}\right]$

$X=\left[\begin{array}{cccccc}& x\left(2\right)& & x\left(2\right)& ...& x(3m+1)\\ & x\left(3\right)& & x\left(3\right)& ...& x(3m+2)\\ & x\left(4\right)& & x\left(4\right)& ...& x(3m+3)\\ ...& & ...& & ...& ...\\ & x(N-3m+1)& & x(N-3m+2)& ...& x\left(N\right)\end{array}\right]$

$\mathrm{Xb}=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ ...\\ x(N-3m)\end{array}\right];$

(4) structure is with 1, a _i(1), a _i(2) ... a _i(3m) be the equation of coefficient

u ^3m+a _i(1)u ^3m-1+a _i(2)u ^3m-2+…………+a _i(3m-1)u+a _i(3m)＝0

Solve this equation, obtain 3m root u (i), wherein i=1 ~ 3m;

(5) structural matrix A and B,

$A=\left[\begin{array}{cccc}1& 1& ...& 1\\ \mathrm{ai}\left(1\right)& \mathrm{ai}\left(2\right)& ...& \mathrm{ai}\left(3m\right)\\ \mathrm{ai}{\left(1\right)}^2& \mathrm{ai}{\left(2\right)}^2& ...& \mathrm{ai}{\left(3m\right)}^2\\ ...& ...& ...& ...\\ \mathrm{ai}{\left(1\right)}^{N-1}& \mathrm{ai}{\left(2\right)}^{N-1}& ...& \mathrm{ai}{\left(3m\right)}^{N-1}\end{array}\right]$ $B=\left[\begin{array}{c}x\left(1\right)\\ x\left(2\right)\\ x\left(3\right)\\ ......\\ x\left(N\right)\end{array}\right]$

Employing least square method solves its solution is designated as b;

(6) calculating of model frequency f:

$f\left(i\right)=\frac{\arctan \left(\frac{\mathrm{Im}\left(u\left(i\right)\right)}{\mathrm{Re}\left(u\left(i\right)\right)}\right)}{2\mathrm{\πdt}},i=1~3m$

The calculating of Modal Decay factor sigma:

${\mathrm{\σ}}_i=\frac{\ln \left|u\left(i\right)\right|}{\mathrm{dt}},i=1~3m;$

(7) for m mode, the actual value having obtained 3m mode, get distance f (1), f (2), the Frequency point that f (m) is nearest, apart from when identical, get the point that attenuation rate absolute value is less, thus find the attenuation coefficient of corresponding model frequency in 3m point.