CN106096170A - Wind turbines multivariate failure prediction method based on data-driven - Google Patents

Abstract

The present invention relates to Wind turbines failure prediction method, a kind of Wind turbines multivariate failure prediction method based on data-driven.The present invention solves the problem that existing failure prediction method predictablity rate based on data-driven is low.Wind turbines multivariate failure prediction method based on data-driven, the method is to use following steps to realize: step 1: monitored Wind turbines parts are carried out state data acquisition；Step 2: use 5 moving average methods that characteristic quantity is carried out noise reduction process；Step 3: calculate the degree of association R of characteristic quantity and predicting residual useful life；Step 4: set up multivariate least square method supporting vector machine forecast model；Step 5: regularization parameter γ and nuclear parameter σ to multivariate least square method supporting vector machine forecast model²It is optimized；Step 6: the effectiveness of checking multivariate least square method supporting vector machine forecast model；Step 7: the remaining useful life of prediction Wind turbines parts.The present invention is applicable to Wind turbines failure predication.

Description

Wind turbines multivariate failure prediction method based on data-driven

Technical field

The present invention relates to Wind turbines failure prediction method, a kind of Wind turbines multivariate based on data-driven Failure prediction method.

Background technology

Failure predication refers to predict following fault trend or remaining useful life according to historical data and current data.Wind Group of motors failure predication refers in the case of Wind turbines not yet breaks down, and uses the method for fault reasoning to Wind turbines The following fault that will occur is made prediction.Relative to Fault monitoring and diagnosis, the advantage of failure predication maximum is: if fault The Timing Advance long enough of prediction, and failure cause is enough accurate with location, and field operator just can take accordingly Measure prevents the generation of fault, improves the stability of Wind turbines.But running of wind generating set state has serious mostly Non-linear, the feature such as strong coupling, whole system between the uncertainty of time variation, structure and parameter and multiple-input and multiple-output Huge and complicated, bad environments, Frequent Troubles.These all increase the difficulty of Wind turbines failure predication, so fault is pre- Survey method should possess: can various complex characteristics (such as non-linear, instability, uncertainty etc.) in processing system；Prediction Accuracy rate is high, and rate of false alarm is low；Prediction pre-set time should enough staff's handling failures；Fault location is wanted accurately, it is possible to instruct Staff fixes a breakdown.

Existing Wind turbines failure prediction method includes the following two kinds: Qualitative Forecast Methods and quantitative forecasting technique.Fixed Property Forecasting Methodology, namely Knowledge based engineering method, mainly have rough set, data mining (KDD), support vector machine (SVM), Petri network and specialist system etc..The advantage of Qualitative Forecast Methods is knowledge or the experience that can make full use of expert, but this kind Forecasting Methodology is applicable to qualitative reasoning and has significant limitation in terms of quantitative Analysis.Quantitative forecasting technique includes following two Kind: Forecasting Methodology based on model and Forecasting Methodology based on data-driven.Wherein, Forecasting Methodology based on model refers to utilize Equipment to-be is predicted by the mathematical model having built up, it is achieved the failure prediction of equipment to-be.Based on mould The Forecasting Methodology of type has degree of accuracy height, rate of false alarm is low, be conducive to the advantages such as fault location, but in Wind turbines failure predication It is very restricted and (due to the impact of the factors such as the serious non-linear and various disturbance of Wind turbines self, noise, is difficult to Set up accurate Wind turbines mathematical model).And Forecasting Methodology of based on data-driven is from wind energy turbine set actual operating data, Running of wind generating set status data can be utilized to carry out the hidden feature of excavating equipment.This kind of Forecasting Methodology is widely used, practicality By force, need not founding mathematical models, it has also become the future library of failure predication.

At present, Forecasting Methodology based on data-driven mainly has classical time series method models such as () AR, MA, ARMA, ash Color prediction, regression analysis, artificial intelligence and chaology prediction etc..These methods are single argument and predict (i.e. according to single The inherent Changing Pattern of variable obtains future value, thus carries out Wind turbines failure predication), its predictablity rate is completely dependent on In single variable.But in actual applications, single variable suffers from the impact of multiple variable, and therefore said method is generally deposited In the problem that predictablity rate is low.It is necessary to invent a kind of brand-new Wind turbines failure prediction method for this, existing to solve The problems referred to above that failure prediction method based on data-driven exists.

Summary of the invention

The present invention is to solve the problem that existing failure prediction method predictablity rate based on data-driven is low, it is provided that A kind of Wind turbines multivariate failure prediction method based on data-driven.

The present invention adopts the following technical scheme that realization:

Wind turbines multivariate failure prediction method based on data-driven, the method is to use following steps to realize:

Step 1: monitored Wind turbines parts are carried out state data acquisition, and extracts following special from status data The amount of levying: temporal signatures amount, frequency domain character amount, complexity characteristics amount；

Step 2: use 5 moving average methods that characteristic quantity carries out noise reduction process, thus eliminate between characteristic quantity is random Property impact；Noise reduction process formula is expressed as follows:

$\left\{\begin{array}{c}{x}^{\′}\left(1\right)=\frac{1}{5}\left(3x\right(1)+2x(2)+x(3)-x(4\left)\right)\\ x^{\′}\left(2\right)=\frac{1}{10}\left(4x\right(1)+3x(2)+2x(3)+x(4\left)\right)\\ .\\ .\\ .\\ x^{\′}\left(k\right)=\frac{1}{5}\left(x\right(k-2)+x(k-1)+x(k+1)+x(k+2\left)\right)\\ .\\ .\\ .\\ x^{\′}(n-1)=\frac{1}{10}\left(x\right(n-3)+2x(n-2)+3x(n-1)+4x(n\left)\right)\\ x^{\′}\left(n\right)=\frac{1}{5}(-x(n-3)+x(n-2)+2x(n-1)+3x(n\left)\right)\end{array}\right.---\left(1\right);$

In formula (1): x (k) represents the characteristic quantity data sequence before noise reduction；X'(k) the characteristic quantity data after noise reduction are represented Sequence；K=1,2 ..., n；

In the steady operation period of Wind turbines parts, the status data calculating the root-mean-square value in characteristic quantity corresponding is average Value, and status data meansigma methods is set to standard value, then by root-mean-square value divided by standard value, thus obtain relative root-mean-square Value；

With relative root-mean-square value for label index, set degeneration initiation threshold d of Wind turbines parts_maxAnd failure threshold y_max；

Step 3: calculate the degree of association R of characteristic quantity and predicting residual useful life；Computing formula is expressed as follows:

$R=\frac{\underset{k=1}{\overset{N}{\Σ}}\[x\left(k\right)-\stackrel{\‾}{x}\]\[y\left(k\right)-\stackrel{\‾}{y}\]}{\sqrt{\underset{k=1}{\overset{N}{\Σ}}{\[x\left(k\right)-\stackrel{\‾}{x}\]}^2{\[y\left(k\right)-\stackrel{\‾}{y}\]}^2}}---\left(2\right);$

In formula (2): y (k) represents the residual life data sequence corresponding with x (k)；Represent data sequence x (k) Average；Represent the average of data sequence y (k)；

The characteristic quantity more than 0.6 of the degree of association R with predicting residual useful life is selected to support vector as multivariate least square The input variable of machine forecast model, removes characteristic quantity incoherent with predicting residual useful life and noise characteristic amount simultaneously；Multivariate The input variable of least square method supporting vector machine forecast model includes characteristic quantity following six: root-mean-square value, absolute average, peak Value, kurtosis index, margin index, wavelet packet-envelope sample entropy；

Step 4: set up multivariate least square method supporting vector machine forecast model；Specifically set up process as follows:

Assume univariate time series be X=[x (1), x (2) ..., x (N)], then according to Takens embedding theory, monotropic Amount seasonal effect in time series predictive value is relevant to its front m value, is expressed as follows:

X (k+1)=f (x (k), x (k-1) ..., x (k-m+1)) (3)；

In formula (3): m represents Embedded dimensions, typically take 5～10；K=m, m+1 ..., N；

Choosing front n for training sample, rear N-n is forecast sample, then the input and output sequence of training sample and prediction The input and output sequence of sample is expressed as follows respectively:

$\begin{array}{cc}{X}_{train}=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& \mathrm{...}& x\left(m\right)\\ x\left(2\right)& x\left(3\right)& \mathrm{...}& x(m+1)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x(n-m)& x(n-m+1)& \mathrm{...}& x(n-1)\end{array}\right]& Y_{train}=\left[\begin{array}{c}x(m+1)\\ x(m+2)\\ .\\ .\\ .\\ x\left(n\right)\end{array}\right]\end{array}---\left(4\right);$

$\begin{array}{cc}{X}_{test}=\left[\begin{array}{cccc}x(n-m+1)& x(n-m+2)& \mathrm{...}& x\left(n\right)\\ z(n-m+2)& x(n-m+3)& \mathrm{...}& x(n+1)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x(N-m+1)& x(N-m+2)& \mathrm{...}& x\left(N\right)\end{array}\right]& Y_{test}=\left[\begin{array}{c}x(n+1)\\ x(n+2)\\ .\\ .\\ .\\ x(N+1)\end{array}\right]\end{array}---\left(5\right);$

In formula (4)-(5): X_train、Y_trainRepresent the input and output sequence of training sample；X_test、Y_testRepresent pre-test sample This input and output sequence；

Defined variable L (k) is by M dimensional feature amount x₁(k),x₂(k),...,x_MK () affects, then variables L (k) is expressed as follows:

L (k)=f (x₁(k),x₂(k),...,x_M(k)) (6)；

Formula (6) is substituted into formula (4)-(5), then the input and output sequence of training sample and the input and output of forecast sample Sequence is expressed as follows respectively:

$\begin{array}{cc}{X}_{train}^{\′}=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_2\left(1\right)& \mathrm{...}& x_M\left(1\right)\\ x_1\left(2\right)& x_2\left(2\right)& \mathrm{...}& x_M\left(2\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_1\left(n\right)& x_2\left(n\right)& \mathrm{...}& x_M\left(n\right)\end{array}\right]& Y_{train}^{\′}=\left[\begin{array}{c}L\left(1\right)\\ L\left(2\right)\\ .\\ .\\ .\\ L\left(n\right)\end{array}\right]\end{array}---\left(7\right);$

$\begin{array}{cc}{X}_{test}^{\′}=\left[\begin{array}{cccc}{x}_1(n+1)& x_2(n+1)& \mathrm{...}& x_M(n+1)\\ x_1(n+2)& x_2(n+2)& \mathrm{...}& x_M(n+2)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_1\left(N\right)& x_2\left(N\right)& \mathrm{...}& x_M\left(N\right)\end{array}\right]& Y_{test}^{\′}=\left[\begin{array}{c}L(n+1)\\ L(n+2)\\ .\\ .\\ .\\ L\left(N\right)\end{array}\right]\end{array}---\left(8\right);$

In formula (7)-(8): X'_train、Y'_trainRepresent the input and output sequence of multivariate training sample；X'_test、Y'_test Represent the input and output sequence of multivariable prediction sample；

Step 5: use the regularization to multivariate least square method supporting vector machine forecast model of the heuristic particle cluster algorithm Parameter γ and nuclear parameter σ²It is optimized；Concrete optimization process is as follows:

Step 5.1: following parameter is initialized: population scale, Studying factors, maximum iteration time, each particle Initial position and speed, individual optimal value, global optimum；

Step 5.2: stochastic generation particle position, and particle position coordinate is set to multivariate least square method supporting vector machine The regularization parameter γ of forecast model and nuclear parameter σ², then use training sample pre-to multivariate least square method supporting vector machine Survey model and be trained prediction；Calculate mean square error MSE, and using mean square error MSE as the fitness value of each particle；

Step 5.3: the fitness value of each particle is compared with individual optimal value and global optimum, thereby determines that Global optimum's particle；

Step 5.4: judge whether global optimum's particle meets end condition；If being unsatisfactory for end condition, then according to formula (9)-(10) speed of more new particle and position, and return step 5.2；If meeting end condition, then forward step 5.5 to；

v_id(t+1)=wv_id(t)+c₁r₁(p_id(t)-x_id(t))+c₂r₂(p_gd(t)-x_id(t)) (9)；

x_id(t+1)=x_id(t)+v_id(t+1) (10)；

In formula (9)-(10): w is inertia weight coefficient；c₁、c₂It is two Studying factors, takes c here₁=1.5, c₂= 1.7；r₁、r₂For the random number between [0, l]；x_id(t)、v_idT () represents the Position And Velocity of each particle respectively；p_id(t) table Show the optimal location that up to the present each particle is occurred；p_gdT () represents the optimum that up to the present all particles are occurred Position；

Step 5.5: find the global optimum's particle position meeting end condition, and global optimum's particle position coordinate is made Optimum regularization parameter γ and optimum nuclear parameter σ for multivariate least square method supporting vector machine forecast model²；

Step 6: structure multivariable nonlinearity function, thus verifies multivariate least square method supporting vector machine forecast model Effectiveness；Multivariable nonlinearity function representation is as follows:

$\left\{\begin{array}{c}y=110-1.8{x_1}^3-2x_2-1.5x_1{x}_2-0.8x_1{x}_2{x}_3+N(0,{0.1}^2)\\ x_1~U\[0,3\]+N(0,{0.1}^2)\\ x_2~U\[0,5.8\]+N(0,{0.3}^2)\\ x_3~U\[0,1.5\]\end{array}\right.---\left(11\right);$

In formula (11): x₁、x₂、x₃Represent the input of multivariate least square method supporting vector machine forecast model, by uniformly dividing Cloth function superposes with white noise；Y represents the output of multivariate least square method supporting vector machine forecast model；

Step 7: according to multivariate least square method supporting vector machine forecast model, it was predicted that the residue of Wind turbines parts is effective Life-span；Concrete prediction process is as follows:

Degeneration initiation threshold d according to Wind turbines parts_maxWith failure threshold y_max, not up to lose at relative root-mean-square value Effect threshold value y_maxBefore, have:

||y_k| | ＜ y_max(12)；

Then the remaining useful life of Wind turbines parts is expressed as follows:

$RUL\left(k\right)=min\{p:||{\hat{y}}_k+p||\≥y_{max}\}---\left(13\right);$

In formula (12)-(13): RUL represents the remaining useful life of Wind turbines parts；P represents prediction step.

Compared with existing failure prediction method based on data-driven, wind turbine based on data-driven of the present invention Group multivariate failure prediction method is by setting up multivariate least square method supporting vector machine forecast model, and by by multiple features Measure the input variable as forecast model, it is achieved that the remaining useful life of Wind turbines parts is directly predicted, therefore Its predictablity rate is no longer dependent on single variable, but is together decided on by multiple variablees, thus significantly improves prediction accurately Rate.

The present invention efficiently solves the problem that existing failure prediction method predictablity rate based on data-driven is low, is suitable for In Wind turbines failure predication.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of the present invention.

Fig. 2 is the schematic flow sheet of step 5 in the present invention.

Fig. 3 is data before 8 vibration temporal signatures index noise reductions.

Fig. 4 is data after 8 vibration temporal signatures index noise reductions.

Fig. 5 is the fitness curve of particle cluster algorithm.

Fig. 6 is predicting the outcome of remaining useful life RUL.

Detailed description of the invention

Wind turbines multivariate failure prediction method based on data-driven, the method is to use following steps to realize:

Step 1: monitored Wind turbines parts are carried out state data acquisition, and extracts following special from status data The amount of levying: temporal signatures amount, frequency domain character amount, complexity characteristics amount；

Step 2: use 5 moving average methods that characteristic quantity carries out noise reduction process, thus eliminate between characteristic quantity is random Property impact；Noise reduction process formula is expressed as follows:

$\left\{\begin{array}{c}{x}^{\′}\left(1\right)=\frac{1}{5}\left(3x\right(1)+2x(2)+x(3)-x(4\left)\right)\\ x^{\′}\left(2\right)=\frac{1}{10}\left(4x\right(1)+3x(2)+2x(3)+x(4\left)\right)\\ .\\ .\\ .\\ x^{\′}\left(k\right)=\frac{1}{5}\left(x\right(k-2)+x(k-1)+x(k+1)+x(k+2\left)\right)\\ .\\ .\\ .\\ x^{\′}(n-1)=\frac{1}{10}\left(x\right(n-3)+2x(n-2)+3x(n-1)+4x(n\left)\right)\\ x^{\′}\left(n\right)=\frac{1}{5}(-x(n-3)+x(n-2)+2x(n-1)+3x(n\left)\right)\end{array}\right.---\left(1\right);$

In formula (1): x (k) represents the characteristic quantity data sequence before noise reduction；X'(k) the characteristic quantity data after noise reduction are represented Sequence；K=1,2 ..., n；

In the steady operation period of Wind turbines parts, the status data calculating the root-mean-square value in characteristic quantity corresponding is average Value, and status data meansigma methods is set to standard value, then by root-mean-square value divided by standard value, thus obtain relative root-mean-square Value；

With relative root-mean-square value for label index, set degeneration initiation threshold d of Wind turbines parts_maxAnd failure threshold y_max；

Step 3: calculate the degree of association R of characteristic quantity and predicting residual useful life；Computing formula is expressed as follows:

$R=\frac{\underset{k=1}{\overset{N}{\Σ}}\[x\left(k\right)-\stackrel{\‾}{x}\]\[y\left(k\right)-\stackrel{\‾}{y}\]}{\sqrt{\underset{k=1}{\overset{N}{\Σ}}{\[x\left(k\right)-\stackrel{\‾}{x}\]}^2{\[y\left(k\right)-\stackrel{\‾}{y}\]}^2}}---\left(2\right);$

In formula (2): y (k) represents the residual life data sequence corresponding with x (k)；Represent data sequence x (k) Average；Represent the average of data sequence y (k)；

The characteristic quantity more than 0.6 of the degree of association R with predicting residual useful life is selected to support vector as multivariate least square The input variable of machine forecast model, removes characteristic quantity incoherent with predicting residual useful life and noise characteristic amount simultaneously；Multivariate The input variable of least square method supporting vector machine forecast model includes characteristic quantity following six: root-mean-square value, absolute average, peak Value, kurtosis index, margin index, wavelet packet-envelope sample entropy；

Step 4: set up multivariate least square method supporting vector machine forecast model；Specifically set up process as follows:

Assume univariate time series be X=[x (1), x (2) ..., x (N)], then according to Takens embedding theory, monotropic Amount seasonal effect in time series predictive value is relevant to its front m value, is expressed as follows:

X (k+1)=f (x (k), x (k-1) ..., x (k-m+1)) (3)；

In formula (3): m represents Embedded dimensions, typically take 5～10；K=m, m+1 ..., N；

Choosing front n for training sample, rear N-n is forecast sample, then the input and output sequence of training sample and prediction The input and output sequence of sample is expressed as follows respectively:

$\begin{array}{cc}{X}_{train}=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& \mathrm{...}& x\left(m\right)\\ x\left(2\right)& x\left(3\right)& \mathrm{...}& x(m+1)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x(n-m)& x(n-m+1)& \mathrm{...}& x(n-1)\end{array}\right]& Y_{train}=\left[\begin{array}{c}x(m+1)\\ x(m+2)\\ .\\ .\\ .\\ x\left(n\right)\end{array}\right]\end{array}---\left(4\right);$

$\begin{array}{cc}{X}_{test}=\left[\begin{array}{cccc}x(n-m+1)& x(n-m+2)& \mathrm{...}& x\left(n\right)\\ z(n-m+2)& x(n-m+3)& \mathrm{...}& x(n+1)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x(N-m+1)& x(N-m+2)& \mathrm{...}& x\left(N\right)\end{array}\right]& Y_{test}=\left[\begin{array}{c}x(n+1)\\ x(n+2)\\ .\\ .\\ .\\ x(N+1)\end{array}\right]\end{array}---\left(5\right);$

In formula (4)-(5): X_train、Y_trainRepresent the input and output sequence of training sample；X_test、Y_testRepresent pre-test sample This input and output sequence；

Defined variable L (k) is by M dimensional feature amount x₁(k),x₂(k),...,x_MK () affects, then variables L (k) is expressed as follows:

L (k)=f (x₁(k),x₂(k),...,x_M(k)) (6)；

Formula (6) is substituted into formula (4)-(5), then the input and output sequence of training sample and the input and output of forecast sample Sequence is expressed as follows respectively:

$\begin{array}{cc}{X}_{train}^{\′}=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_2\left(1\right)& \mathrm{...}& x_M\left(1\right)\\ x_1\left(2\right)& x_2\left(2\right)& \mathrm{...}& x_M\left(2\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_1\left(n\right)& x_2\left(n\right)& \mathrm{...}& x_M\left(n\right)\end{array}\right]& Y_{train}^{\′}=\left[\begin{array}{c}L\left(1\right)\\ L\left(2\right)\\ .\\ .\\ .\\ L\left(n\right)\end{array}\right]\end{array}---\left(7\right);$

$\begin{array}{cc}{X}_{test}^{\′}=\left[\begin{array}{cccc}{x}_1(n+1)& x_2(n+1)& \mathrm{...}& x_M(n+1)\\ x_1(n+2)& x_2(n+2)& \mathrm{...}& x_M(n+2)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x_1\left(N\right)& x_2\left(N\right)& \mathrm{...}& x_M\left(N\right)\end{array}\right]& Y_{test}^{\′}=\left[\begin{array}{c}L(n+1)\\ L(n+2)\\ .\\ .\\ .\\ L\left(N\right)\end{array}\right]\end{array}---\left(8\right);$

In formula (7)-(8): X'_train、Y'_trainRepresent the input and output sequence of multivariate training sample；X'_test、Y'_test Represent the input and output sequence of multivariable prediction sample；

Step 5: use the regularization to multivariate least square method supporting vector machine forecast model of the heuristic particle cluster algorithm Parameter γ and nuclear parameter σ²It is optimized；Concrete optimization process is as follows:

Step 5.1: following parameter is initialized: population scale, Studying factors, maximum iteration time, each particle Initial position and speed, individual optimal value, global optimum；

Step 5.2: stochastic generation particle position, and particle position coordinate is set to multivariate least square method supporting vector machine The regularization parameter γ of forecast model and nuclear parameter σ², then use training sample pre-to multivariate least square method supporting vector machine Survey model and be trained prediction；Calculate mean square error MSE, and using mean square error MSE as the fitness value of each particle；

Step 5.3: the fitness value of each particle is compared with individual optimal value and global optimum, thereby determines that Global optimum's particle；

Step 5.4: judge whether global optimum's particle meets end condition；If being unsatisfactory for end condition, then according to formula (9)-(10) speed of more new particle and position, and return step 5.2；If meeting end condition, then forward step 5.5 to；

v_id(t+1)=wv_id(t)+c₁r₁(p_id(t)-x_id(t))+c₂r₂(p_gd(t)-x_id(t)) (9)；

x_id(t+1)=x_id(t)+v_id(t+1) (10)；

In formula (9)-(10): w is inertia weight coefficient；c₁、c₂It is two Studying factors, takes c here₁=1.5, c₂= 1.7；r₁、r₂For the random number between [0, l]；x_id(t)、v_idT () represents the Position And Velocity of each particle respectively；p_id(t) table Show the optimal location that up to the present each particle is occurred；p_gdT () represents the optimum that up to the present all particles are occurred Position；

Step 5.5: find the global optimum's particle position meeting end condition, and global optimum's particle position coordinate is made Optimum regularization parameter γ and optimum nuclear parameter σ for multivariate least square method supporting vector machine forecast model²；

Step 6: structure multivariable nonlinearity function, thus verifies multivariate least square method supporting vector machine forecast model Effectiveness；Multivariable nonlinearity function representation is as follows:

$\left\{\begin{array}{c}y=110-1.8{x_1}^3-2x_2-1.5x_1{x}_2-0.8x_1{x}_2{x}_3+N(0,{0.1}^2)\\ x_1~U\[0,3\]+N(0,{0.1}^2)\\ x_2~U\[0,5.8\]+N(0,{0.3}^2)\\ x_3~U\[0,1.5\]\end{array}\right.---\left(11\right);$

In formula (11): x₁、x₂、x₃Represent the input of multivariate least square method supporting vector machine forecast model, by uniformly dividing Cloth function superposes with white noise；Y represents the output of multivariate least square method supporting vector machine forecast model；

Step 7: according to multivariate least square method supporting vector machine forecast model, it was predicted that the residue of Wind turbines parts is effective Life-span；Concrete prediction process is as follows:

Degeneration initiation threshold d according to Wind turbines parts_maxWith failure threshold y_max, not up to lose at relative root-mean-square value Effect threshold value y_maxBefore, have:

||y_k| | ＜ y_max(12)；

Then the remaining useful life of Wind turbines parts is expressed as follows:

$RUL\left(k\right)=min\{p:||{\hat{y}}_k+p||\≥y_{max}\}---\left(13\right);$

In formula (12)-(13): RUL represents the remaining useful life of Wind turbines parts；P represents prediction step.

In order to verify the pre-of multivariate least square method supporting vector machine forecast model (PSO-MLSSVM model) in the present invention Surveying effect, uniformly generate 220 groups of sequences, remove above 50 groups of stable sequences, middle 130 groups of sequences are used for training, and last 40 Group data are used for testing；

It is respectively adopted PSO-MLSSVM and PSO-LSSVM model to be predicted, first PSO algorithm is carried out parameter initial Change: maximum iteration time is set to 300, and population scale is 30, regularization parameter γ ∈ [100,2000] and nuclear parameter σ²∈ [1, 100]；Model, through training and test, carries out parameter optimization, obtains parameters optimization γ=929.8 and σ²=14.7；

Good and bad in order to specifically evaluate the performance of two kinds of models, introduce mean absolute error (MAE), average relative error (MAPE), the evaluation index such as root-mean-square error (RMSE), the expression formula of each index is as follows:

$\left\{\begin{array}{c}MAE=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}|y_k^{\′}-y_k|\\ MRE=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}\left|\frac{y_k^{\′}-y_k}{y_k}\right|\\ RMSE=\sqrt{\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}{(y_k^{\′}-y_k)}^2}\end{array}\right.---\left(14\right);$

In formula (14): y_kRepresent actual value；y'_kRepresent predictive value；

Obtain two kinds of forecast model evaluation index contrast tables 1, it can be seen that the PSO-MLSSVM model that the present invention proposes is each Item index is superior to PSO-LSSVM, and macro-forecast is respond well；And PSO-LSSVM model is due to by data randomness and noise Impact, it was predicted that effect is poor, demonstrates the multivariate model effectiveness of proposition.

Table 1:

For further illustrating the feasibility of multivariate least square method supporting vector machine forecast model in the present invention, by multivariate Least square method supporting vector machine forecast model is applied in the predicting residual useful life of gearbox of wind turbine bearing: with Shanxi Province The gearbox of wind turbine speed end bearing of wind energy turbine set is object of study, obtains rotating speed in 2014 in on-line monitoring system and exists 1740-1800RPM, power time domain, frequency domain and the vibration number of wavelet packet-envelope sample entropy in the range of 1500-1550KW According to, totally 1190 groups.The degraded data of part time domain vibration index is as shown in Figure 3.As seen in Figure 3: each characteristic quantity with Machine undulatory property is relatively big, and fault threshold is difficult to accurately set, and needs to carry out feature noise reduction process and relative normalized, according to this The fault threshold method to set up of bright proposition, the initial value arranging relative root-mean-square value in conjunction with actual data change trend is 1.3 (figures In the 967th point) and stale value be 3.0 (the 1036th points in figure), such as Fig. 4.Intercept the 967-1036 axle putting totally 70 groups Hold degraded data, predict for bearing remaining useful life.Therefore, it can calculate correspondence according to the threshold value bound arranged Residual life, then utilizes multivariate least square method supporting vector machine forecast model to carry out failure predication.The degeneration number that will intercept According to front 55 groups as training data, rear 15 groups as test data, model parameter after particle cluster algorithm optimization is γ =1822.2 and σ²=20.6, Fig. 5 are the fitness curve of particle cluster algorithm.Predict the outcome such as Fig. 6, lists file names with evaluation index Table 2, relative error is only 7.93% as can be seen from Table 2, and forecast result of model is good, demonstrates multivariate predictive model Effectiveness.

Table 2:

Claims (1)

1. a Wind turbines multivariate failure prediction method based on data-driven, it is characterised in that: the method is to use such as Lower step realizes:

Step 1: monitored Wind turbines parts are carried out state data acquisition, and extracts following feature from status data Amount: temporal signatures amount, frequency domain character amount, complexity characteristics amount；

Step 2: use 5 moving average methods that characteristic quantity carries out noise reduction process, thus eliminate the randomness shadow between characteristic quantity Ring；Noise reduction process formula is expressed as follows:

$\left\{\begin{array}{c}{x}^{\′}\left(1\right)=\frac{1}{5}\left(3x\right(1)+2x(2)+x(3)-x(4\left)\right)\\ x^{\′}\left(2\right)=\frac{1}{10}\left(4x\right(1)+3x(2)+2x(3)+x(4\left)\right)\\ .\\ .\\ .\\ x^{\′}\left(k\right)=\frac{1}{5}\left(x\right(k-2)+x(k-1)+x(k+1)+x(k+2\left)\right)\\ .\\ .\\ .\\ x^{\′}(n-1)=\frac{1}{10}\left(x\right(n-3)+2x(n-2)+3x(n-1)+4x(n\left)\right)\\ x^{\′}\left(n\right)=\frac{1}{5}(-x(n-3)+x(n-2)+2x(n-1)+3x(n\left)\right)\end{array}\right.---\left(1\right);$

In formula (1): x (k) represents the characteristic quantity data sequence before noise reduction；X'(k) the characteristic quantity data sequence after noise reduction is represented； K=1,2 ..., n；

In the steady operation period of Wind turbines parts, calculate the status data meansigma methods that the root-mean-square value in characteristic quantity is corresponding, and Status data meansigma methods is set to standard value, then by root-mean-square value divided by standard value, thus obtains relative root-mean-square value；

With relative root-mean-square value for label index, set degeneration initiation threshold d of Wind turbines parts_maxWith failure threshold y_max；

Step 3: calculate the degree of association R of characteristic quantity and predicting residual useful life；Computing formula is expressed as follows:

$R=\frac{\underset{k=1}{\overset{N}{\Σ}}\[x\left(k\right)-\stackrel{\‾}{x}\]\[y\left(k\right)-\stackrel{\‾}{y}\]}{\sqrt{\underset{k=1}{\overset{N}{\Σ}}{\[x\left(k\right)-\stackrel{\‾}{x}\]}^2{\[y\left(k\right)-\stackrel{\‾}{y}\]}^2}}---\left(2\right);$

In formula (2): y (k) represents the residual life data sequence corresponding with x (k)；Represent the average of data sequence x (k)；Represent the average of data sequence y (k)；

Select pre-as multivariate least square method supporting vector machine with the degree of association R of the predicting residual useful life characteristic quantity more than 0.6 Survey the input variable of model, remove characteristic quantity incoherent with predicting residual useful life and noise characteristic amount simultaneously；Multivariate is minimum Two take advantage of the input variable of SVM prediction model to include characteristic quantity following six: root-mean-square value, absolute average, peak value, Kurtosis index, margin index, wavelet packet-envelope sample entropy；

Step 4: set up multivariate least square method supporting vector machine forecast model；Specifically set up process as follows:

Assume univariate time series be X=[x (1), x (2) ..., x (N)], then according to Takens embedding theory, during single argument Between the predictive value of sequence relevant to its front m value, be expressed as follows:

X (k+1)=f (x (k), x (k-1) ..., x (k-m+1)) (3)；

In formula (3): m represents Embedded dimensions, typically take 5～10；K=m, m+1 ..., N；

Choosing front n for training sample, rear N-n is forecast sample, then the input and output sequence of training sample and forecast sample Input and output sequence be expressed as follows respectively:

$\begin{array}{cc}{X}_{train}=\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& ...& x\left(m\right)\\ x\left(2\right)& x\left(3\right)& ...& x(m+1)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ x(n-m)& x(n-m+1)& ...& x(n-1)\end{array}\right]& Y_{train}=\left[\begin{array}{c}x(m+1)\\ x(m+2)\\ .\\ .\\ .\\ x\left(n\right)\end{array}\right]\end{array}---\left(4\right);$

$\begin{array}{cc}{X}_{test}=\left[\begin{array}{cccc}x(n-m+1)& x(n-m+2)& ...& x\left(n\right)\\ z(n-m+2)& x(n-m+3)& ...& x(n+1)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ x(N-m+1)& x(N-m+2)& ...& x\left(N\right)\end{array}\right]& Y_{test}=\left[\begin{array}{c}x(n+1)\\ x(n+2)\\ .\\ .\\ .\\ x(N+1)\end{array}\right]\end{array}---\left(5\right);$

In formula (4)-(5): X_train、Y_trainRepresent the input and output sequence of training sample；X_test、Y_testRepresent forecast sample Input and output sequence；

Defined variable L (k) is by M dimensional feature amount x₁(k),x₂(k),...,x_MK () affects, then variables L (k) is expressed as follows:

L (k)=f (x₁(k),x₂(k),...,x_M(k)) (6)；

Formula (6) is substituted into formula (4)-(5), then the input and output sequence of training sample and the input and output sequence of forecast sample It is expressed as follows respectively:

$\begin{array}{cc}{X}_{train}^{\′}=\left[\begin{array}{cccc}{x}_1\left(1\right)& x_2\left(1\right)& ...& x_M\left(1\right)\\ x_1\left(2\right)& x_2\left(2\right)& ...& x_M\left(2\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ x_1\left(n\right)& x_2\left(n\right)& ...& x_M\left(n\right)\end{array}\right]& Y_{train}^{\′}=\left[\begin{array}{c}L\left(1\right)\\ L\left(2\right)\\ .\\ .\\ .\\ L\left(n\right)\end{array}\right]\end{array}---\left(7\right);$

$\begin{array}{cc}{X}_{test}^{\′}=\left[\begin{array}{cccc}{x}_1(n+1)& x_2(n+1)& ...& x_M(n+1)\\ x_1(n+2)& x_2(n+2)& ...& x_M(n+2)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ x_1\left(N\right)& x_2\left(N\right)& ...& x_M\left(N\right)\end{array}\right]& Y_{test}^{\′}=\left[\begin{array}{c}L(n+1)\\ L(n+2)\\ .\\ .\\ .\\ L\left(N\right)\end{array}\right]\end{array}---\left(8\right);$

In formula (7)-(8): X'_train、Y'_trainRepresent the input and output sequence of multivariate training sample；X'_test、Y'_testRepresent The input and output sequence of multivariable prediction sample；

Step 5: use the heuristic particle cluster algorithm regularization parameter to multivariate least square method supporting vector machine forecast model γ and nuclear parameter σ²It is optimized；Concrete optimization process is as follows:

Step 5.1: following parameter is initialized: at the beginning of population scale, Studying factors, maximum iteration time, each particle Beginning position and speed, individual optimal value, global optimum；

Step 5.2: stochastic generation particle position, and particle position coordinate is set to the prediction of multivariate least square method supporting vector machine The regularization parameter γ of model and nuclear parameter σ², then use training sample that multivariate least square method supporting vector machine is predicted mould Type is trained prediction；Calculate mean square error MSE, and using mean square error MSE as the fitness value of each particle；

Step 5.3: the fitness value of each particle is compared with individual optimal value and global optimum, thereby determines that the overall situation Optimal particle；

Step 5.4: judge whether global optimum's particle meets end condition；If being unsatisfactory for end condition, then according to formula (9)- (10) speed of more new particle and position, and return step 5.2；If meeting end condition, then forward step 5.5 to；

v_id(t+1)=wv_id(t)+c₁r₁(p_id(t)-x_id(t))+c₂r₂(p_gd(t)-x_id(t)) (9)；

x_id(t+1)=x_id(t)+v_id(t+1) (10)；

In formula (9)-(10): w is inertia weight coefficient；c₁、c₂It is two Studying factors, takes c here₁=1.5, c₂=1.7； r₁、r₂For the random number between [0, l]；x_id(t)、v_idT () represents the Position And Velocity of each particle respectively；p_idT () represents every The optimal location that up to the present individual particle is occurred；p_gdT () represents the optimal location that up to the present all particles are occurred；

Step 5.5: find the global optimum's particle position meeting end condition, and using global optimum's particle position coordinate as many The optimum regularization parameter γ of variable least square method supporting vector machine forecast model and optimum nuclear parameter σ²；

Step 6: structure multivariable nonlinearity function, thus checking multivariate least square method supporting vector machine forecast model is effective Property；Multivariable nonlinearity function representation is as follows:

$\left\{\begin{array}{c}y=110-1.8{x_1}^3-2x_2-1.5x_1{x}_2-0.8x_1{x}_2{x}_3+N(0,{0.1}^2)\\ x_1~U\[0,3\]+N(0,{0.1}^2)\\ x_2~U\[0,5.8\]+N(0,{0.3}^2)\\ x_3~U\[0,1.5\]\end{array}\right.---\left(11\right);$

In formula (11): x₁、x₂、x₃Represent the input of multivariate least square method supporting vector machine forecast model, by being uniformly distributed letter Number superposes with white noise；Y represents the output of multivariate least square method supporting vector machine forecast model；

Step 7: according to multivariate least square method supporting vector machine forecast model, it was predicted that the residue service life of Wind turbines parts Life；Concrete prediction process is as follows:

Degeneration initiation threshold d according to Wind turbines parts_maxWith failure threshold y_max, not up to lost efficacy threshold at relative root-mean-square value Value y_maxBefore, have:

||y_k| | ＜ y_max(12)；

Then the remaining useful life of Wind turbines parts is expressed as follows:

$RUL\left(k\right)=min\{p:||{\hat{y}}_k+p||\≥y_{max}\}---\left(13\right);$

In formula (12)-(13): RUL represents the remaining useful life of Wind turbines parts；P represents prediction step.