CN116611609B

CN116611609B - Equipment stability state prediction and evaluation method based on monitoring parameters

Info

Publication number: CN116611609B
Application number: CN202310415583.3A
Authority: CN
Inventors: 倪何; 肖鹏飞; 卓越
Original assignee: China Shipbuilding Group Corp 703 Research Institute; Naval University of Engineering PLA
Current assignee: China Shipbuilding Group Corp 703 Research Institute; Naval University of Engineering PLA
Priority date: 2023-04-18
Filing date: 2023-04-18
Publication date: 2024-01-09
Anticipated expiration: 2043-04-18
Also published as: CN116611609A

Abstract

The invention discloses a method for predicting and evaluating equipment stability state based on monitoring parameters, which is specially used for a power generation turbine of a steam power ship. Decomposing the initial time sequence by adopting an MREMD decomposition method, and providing a non-stationary time sequence single-parameter prediction model based on MREMD decomposition and wavelet threshold noise reduction, so that the prediction precision of equipment is improved; the parameter time sequence of each evaluation index is converted into a fluctuation value sequence, so that positive and negative ideal solutions of each evaluation index of the equipment can be determined, the correlation weight of the weighted evaluation parameter sequence of each evaluation index is obtained by combining an improved TOPSIS algorithm and the existing B-type gray correlation analysis method, the final running stability state score of each prediction trend item is obtained, each evaluation index is adjusted according to the stability state score, and the overall stability of the equipment in a future period can be ensured by timely adjusting the correlation parameters.

Description

Equipment stability state prediction and evaluation method based on monitoring parameters

Technical Field

The invention belongs to the technical field of system stability evaluation, and particularly relates to a device stability state prediction and evaluation method based on monitoring parameters.

Background

The coupling between the parameters related to the running state of the current system is strong, so that equipment is difficult to eliminate, and the cost is high. At present, no mature method for predicting future parameter changes based on the current running state of the system exists, so that state evaluation and maintenance of the system are mainly based on state parameters recorded in the moment of system instability or instability, and in fact, before the instability phenomenon occurs, relevant parameters of the system often have precursors such as state parameter drift, running characteristic change and the like, if the state of equipment cannot be found and adjusted in time, the state of the equipment may be further deteriorated, and finally various faults occur.

For example, the unstable operation of power turbines is particularly evident in extreme variable conditions, cross steam supply and special operating environments of steam powered and nuclear powered vessels. Although the monitoring system ensures the stability of the operation state of the power generation turbine by strictly monitoring the operation state parameters of all the equipment in the system in real time, the problems that the parameter state transition time is longer, the preset control parameters and logic cannot adapt to the technical state change of the equipment and the like still exist due to factors such as thermal inertia and the like. For a power generation turbine in a steam power generation system, the parameters influencing the running condition of the power generation turbine are numerous, the coupling relation is complex, the current monitoring system is difficult to ensure that the power generation turbine can maintain the stable state under all special conditions and disturbance, the power generation turbine lacks corresponding early warning and pre-adjustment capability when running is about to be unstable, and the power generation turbine can only rely on a large number of monitoring personnel to make local adjustment for instability alarms generated when the parameter running in the system is unstable, so that the labor and equipment running cost is huge.

In the aspect of multi-parameter correlation analysis, the entropy weight-ideal solution describes the proximity degree of two evaluation objects through distance measurement, and the method can only show the position relation and can not show the situation change situation among the data sequences of the objects to be evaluated; the gray correlation algorithm describes the correlation by the geometric similarity degree of curves of two study objects, and the method can only reflect the situation change condition of a data curve and cannot reflect the position relation; by adopting the parameter correlation calculation model based on the entropy weight method, the ideal solution and the gray correlation analysis, the similarity of the change trend and the position relation between the parameters can be considered on the basis of determining the importance degree of the parameters.

The prior TOPSIS method constructs a positive ideal solution and a negative ideal solution according to main characteristics of a study object, calculates distances between each evaluation index of the equipment and the positive ideal solution and the negative ideal solution, and then takes the distance degree close to the positive ideal solution and the negative ideal solution as a basis for evaluating the importance degree of the influence factors. The approach degree of two evaluation objects is described through distance measurement, the method can only show the position relation, the situation change situation among the data sequences of the objects to be evaluated can not be shown, and for the data containing non-stationary parameters, the negative ideal solution can not be obtained by using the existing TOPSIS method, and the method can not be applied to the stability evaluation of equipment containing the non-stationary parameters.

It is therefore desirable to provide an improved TOPSIS algorithm by which to combine the two methods to stabilize the operation of a device for a future period of time with device parameters that contain non-stationary parameters, in combination with existing type B gray correlation analysis.

Disclosure of Invention

The invention provides a device stability state prediction and evaluation method based on monitoring parameters, which is used for realizing the prediction and evaluation of the device running stability state containing non-stable parameters, and timely adjusting the running operation of related devices through the stability state score, so that the running stability of the devices in a future period of time is improved.

The invention discloses a device stability state prediction and evaluation method based on monitoring parameters, which comprises the following steps:

s1, analyzing an instability factor affecting the running state of equipment, selecting a plurality of evaluation indexes related to the running stability of the equipment, collecting data sets with the length of N of each evaluation index in a period of time, and forming each data set into an initial time sequence of the running of the evaluation index;

s2, decomposing the initial time sequence by adopting an MREMD decomposition method to obtain n IMF components and a final residual component ；

S3, setting a screening threshold value to screen the IMF components, carrying out wavelet threshold value noise reduction treatment on the IMF components meeting the screening threshold value, recombining to obtain a new IMF component matrix, and processing the new IMF component matrix to obtain a final prediction trend item;

s4, converting the new IMF component matrix to obtain a parameter time sequence of each evaluation index, obtaining a weighted evaluation parameter sequence of each evaluation index by using an entropy weight method, and obtaining positive and negative ideal solutions of each moment point by the weighted evaluation parameter sequence;

s5, respectively calculating the distance relation value between the fluctuation value of each evaluation index and the positive and negative ideal solutions;

s6, respectively calculating gray correlation degrees of fluctuation values and positive and negative ideal solutions of all evaluation indexes by adopting a B-type gray correlation analysis method;

s7, carrying out dimensionless treatment on the distance relation value and the gray correlation value respectively, and fusing the gray correlation value and the distance relation value after the dimensionless treatment into a correlation weight of each evaluation index and the integral running stability state of the equipment;

s8, establishing an operation stability state scoring formula of each evaluation index on the basis of the predicted trend item, obtaining the scoring value of each evaluation index at each moment point in a predicted time period, and establishing an operation scoring formula of the predicted trend item of the whole equipment according to the scoring value and the correlation weight of each evaluation index to obtain the operation scoring of the predicted trend item of the whole equipment;

S9, setting a total reference index, and comparing the operation evaluation score of the overall predicted trend item of the equipment with the total reference index to judge whether adjustment of each evaluation index is needed, so that the operation evaluation score of the overall predicted trend item of the equipment is higher than the total reference index finally.

Further, the step S2 includes the steps of:

s21, utilizing an AR prediction model to perform initial time sequence on each evaluation indexPredicting and extending the left and right endpoints to obtain a new sequence so as to enable +.>Is between two extreme points of the new sequence after extension, wherein the AR prediction model expression is as follows:

（1）

in the above-mentioned method, the step of,the predicted point is the t moment; />The actual values of the parameters at the times t-1, t-2, & gtt-p, & lt/L, & gt>For p+1 real numbers, ">(t=p+1, p+2, & gtis, N) white noise sequence with zero mean;

s22, calculating the average value between adjacent extreme points in the new sequence, wherein the average value point obtained in the new sequence is the initial time sequenceMean point sequence>Wherein k is->The number of the mean value points is k+1, and the +.f is obtained through 3 times spline interpolation>Signal mean sequence>Make->Minus- >Obtain 1 order signal component->The expression is shown as follows:

（2）

s23, calculating the obtained signal componentAs an original signal, the above steps S21 and S22 are repeated for iterative calculation until the signal component after r iterations +.>And the termination condition is satisfied, and the judgment formula of the calculation termination condition is as follows:

（3）

in the above-mentioned method, the step of,for normalizing standard deviation, ++>And->For the signal component h after the r-1 th iteration respectively _1,r-1 The mean point sequence of (t) and the signal component after the r-th iteration +.>The standard deviation of the mean point sequence of (2), P is the conditional probability; definition according to sigma principle->For the initial time sequence->A sense value of the absolute value of the middle extreme point; />And->The mean point of the time series after the initial time series and the r iteration and +.>Ratio sequence of (2); said->、/>And->The formulas of (a) are respectively shown in the following formulas;

（4）

wherein:the z extreme point of the new sequence after extension; />For the initial time sequence->Is a k mean points of (2); />For time series component after the r-1 th iteration +.>Is a k mean points of (2);

s24, willAs an eigenmode function IMF of order 1 ₁ Output of->Minus->Get component->And willThe above steps S21 to S24 are repeated as an original signal until the residual signal becomes a monotone function or no new IMF component is separated, and the residual signal extraction process is as follows:

（5）

Wherein n is the number of resolvable IMF componentsA number;as the final residual component;

s25, after the decomposition, the initial time sequenceCan be expressed as the n IMF components and the final residual component->The sum is as follows:

（6）。

further, the step S3 includes the steps of:

s31, according to the IMF components and the initial time sequenceCorrelation coefficient of->And root mean square errorObtaining the screening threshold value;

s32, screening the n IMF components by using the screening threshold, and dividing the IMF components into IMF components which do not meet the screening threshold and IMF components which meet the screening threshold;

s33, carrying out wavelet threshold denoising treatment on the components meeting the screening threshold, and recombining the residual components, the IMF components which do not meet the screening threshold and the IMF components subjected to denoising treatment into a new IMF component matrix;

s34, performing feature extraction and dimension reduction on the new IMF component matrix by using a singular value decomposition method to obtain a singular value component matrix;

s35, calculating and obtaining the arrangement entropy value of each singular value component matrix by using an arrangement entropy algorithm;

s36, dividing the new IMF component matrix into two types according to the arrangement entropy value by utilizing a K-means clustering algorithm, taking a type of component with a lower entropy value, reconstructing the type of component into an original trend item, carrying out differential stabilization treatment on the original trend item, obtaining an ARIMA model by fitting, and predicting the original trend item by the ARIMA model to obtain a final predicted trend item.

Further, the step S31 includes the steps of:

s311, calculating the IMF components of each order and an initial time sequenceCorrelation coefficient of->And root mean square errorThe expression is shown as follows:

（7）

in the middle ofAnd->Initial time series>And the average value of IMF components of each order, N is the sequence length of each parameter, and N is the number of IMF components;

s312, calculating a screening threshold value, wherein the screening threshold value comprises a correlation coefficient threshold valueAnd mean square error threshold->The expressions are shown below:

（8）

s313, selecting IMF components meeting the following conditions, and carrying out wavelet threshold denoising to obtain denoised IMF components, wherein the judgment formula of the screening conditions is as follows:

（9）。

further, the step S32 includes the steps of:

s321, estimating the noise of each layer in the IMF components meeting the screening conditions through wavelet detail coefficients of each layer, and simultaneously utilizing the coefficientsThe threshold value of the detail coefficient is lowered layer by layer, and the selection criterion formula of the threshold value is as follows:

（10）

in the method, in the process of the invention,a noise threshold that is a layer j wavelet detail coefficient; />The detail coefficient of the wavelet of the j-th layer; media is an intermediate function, namely after each coefficient is arranged in descending order, when the coefficient number is odd, the value of the intermediate number is taken, or when the coefficient number is even, the average value of the intermediate two numbers is taken; n is the number of data points collected by each evaluation index;

S322, performing threshold processing by using a composite threshold function to obtain an IMF component subjected to wavelet noise reduction, wherein the expression of the composite threshold function is as follows:

（11）

in the method, in the process of the invention,the j-th layer detail coefficient after the noise reduction treatment; sign is a sign function, a e [0,1] is an adjustment coefficient, a=0 is a hard threshold function, and a=1 is a soft threshold function.

Further, the step S4 includes the steps of:

s41, sequentially carrying out accumulation reconstruction on an IMF component subjected to noise reduction of each parameter, a component not subjected to noise reduction treatment and a residual component to form a noise-reduced parameter time sequence, converting the noise-reduced parameter time sequence into a fluctuation value sequence according to a setting value of each evaluation index, obtaining an entropy weight value of the fluctuation value sequence by using an entropy weight method, and obtaining a weighted evaluation parameter sequence of each evaluation index by using the entropy weight value and the fluctuation value sequence;

s42, respectively selecting the point with the largest positive fluctuation value of the weighted evaluation parameter sequence of each evaluation index at the same moment as a positive ideal solution, and taking the point with the largest negative fluctuation value as a negative ideal solution, wherein the calculated formulas of the positive ideal solution and the negative ideal solution are as follows:

（12）

In the method, in the process of the invention,for positive understanding, add>Is a negative ideal solution; />，/>N represents the number of evaluation indexes.

Further, the expression of the distance relation between the fluctuation value of each evaluation index sequence and the positive and negative ideal solutions is as follows:

（13）

wherein i=1, 2, …, N; j=1, 2, …, n, n represents the number of evaluation indexes,for the distance relation value between the fluctuation value of the j-th evaluation index and the positive ideal solution,/for the j-th evaluation index>And the relationship value of the distance between the fluctuation value of the j-th evaluation index and the negative ideal solution.

Further, the expression of the operation stability state scoring formula of each evaluation index is as follows:

（14）

wherein j=1, 2, …, n;for the weighted evaluation parameter of the jth evaluation index in the ith point in time, +.>Andrespectively representing the positive maximum fluctuation value and the negative maximum fluctuation value of each evaluation index in the running range of the equipment.

Further, the step S9 includes the steps of:

calculating the running evaluation score of the overall predicted trend item of the equipment, and setting a total reference index;

comparing the operation evaluation score of the overall predicted trend item of the equipment with the total reference index;

when the operation evaluation score of the overall prediction trend item of the equipment is larger than or equal to the total reference index, the equipment is indicated to operate well in a future prediction time period;

When the operation evaluation score of the overall prediction trend item of the equipment is smaller than the total reference index, the operation instability of the equipment in a future prediction time period is indicated, and the evaluation index with the lower score is timely adjusted, so that the operation evaluation score of the overall prediction trend item of the equipment is greater than or equal to the total reference index.

Further, the equipment stability state prediction and evaluation method based on the monitoring parameters is specially used for the power generation steam turbine of the steam power ship.

Compared with the prior art, the invention has the beneficial effects that:

1. the invention provides a non-stationary time sequence single parameter prediction model based on MREMD decomposition and wavelet threshold noise reduction, which improves the prediction precision of the model; in the definition of noise IMF components of each IMF component after MREMD decomposition, a subjective IMF component definition method for removing only a first-order IMF component to achieve the purpose of noise reduction is abandoned, and two indexes of a correlation coefficient and a root mean square error are introduced to realize quantitative definition of the IMF component, so that the calculated error is reduced;

2. the invention combines the improved TOPSIS method with the existing B-type gray correlation degree, and converts the noise-reduced parameter time sequence into a fluctuation value sequence according to the setting value of each evaluation index, so that the positive and negative ideal solutions of each evaluation index of the equipment can be determined, and the problem that the negative ideal solutions cannot be obtained when the existing TOPSIS method is applied to the stability evaluation of the equipment containing non-stable parameters is solved;

3. In the method, an improved distance relation value formula is provided when the positive and negative ideal solutions in the improved TOPSIS method are calculated, and the correlation between the fluctuation value of each evaluation index and the positive and negative ideal solutions is calculated through the formula; in addition, as the distance between the two is smaller, the evaluation index is unstable, the weight occupied by the evaluation index in the stability of the whole equipment is larger, namely the weight of the evaluation index and the corresponding distance relation of the evaluation index are in inverse proportion relation, the negative sign is added in the formula of the distance relation value, so that the real situation can be better reflected, and the obtained result is more accurate;

4. according to the invention, through establishing the equipment operation evaluation score and evaluating the whole predicted trend item of the equipment, the state stability evaluation score curve is generated through software, so that the operation condition of the power generation turbine in a period of time in the future can be intuitively observed, the operation state of the equipment can be observed in time by operation operators, and the labor and the equipment operation cost are saved.

Drawings

FIG. 1 is a flow chart of a single evaluation index timing prediction for a power generation turbine of a steam powered vessel of the present invention;

FIG. 2 shows 7 IMF components and 1 residual component R obtained by decomposing the rotational speed of the power turbine of the present embodiment _n ；

FIG. 3 is an initial time series x of 7 IMF components of the turbine speed for power generation according to the present embodiment ₀ A correlation coefficient threshold screening map of (t);

FIG. 4 is an initial time series x of 7 IMF components of the turbine speed of the power generation turbine of the present embodiment ₀ A mean square error threshold screening graph of (t);

FIG. 5 is a schematic diagram of singular value components of a new IMF component matrix of the power turbine speed of the present embodiment;

FIG. 6 is a graph showing the variation of mutual information values of components of the rotational speed of the power generation turbine with delay time according to the present embodiment;

FIG. 7 is a graph showing the variation of pseudo-neighbor ratios of the various components of the speed of the power turbine according to the present embodiment with the embedding dimension;

FIG. 8 is a schematic diagram of the rotational speed of the power turbine and the original trend thereof according to the present embodiment;

FIG. 9 is a schematic diagram of the rotational speed of the power turbine, the raw trend term and the predicted trend term according to the present embodiment;

FIG. 10 is a diagram showing a positive ideal understanding of the weighted evaluation parameter sequence of the evaluation index according to the present embodiment;

FIG. 11 is a diagram showing a negative idealized solution of a weighted evaluation parameter sequence of the evaluation index according to the present embodiment;

Fig. 12 is a schematic diagram showing evaluation scores of the stability state of the power turbine of the present embodiment over a future period of time.

Detailed Description

The invention is described in further detail below with reference to figures 1 to 12 and the specific examples.

The equipment stability state prediction and evaluation method based on the monitoring parameters is specially used for the power generation steam turbine of the steam power ship.

As shown in fig. 1, a device stability state prediction and evaluation method based on monitoring parameters in this embodiment includes the following steps:

s1, analyzing instability factors affecting the running state of equipment, selecting a plurality of evaluation indexes related to the running stability of the equipment, collecting data sets with the length N of each evaluation index in a monitoring system in a period of time, and forming an initial time sequence x of the running of the evaluation indexes by each data set ₀ (t)。

The power generation turbine unit of the embodiment comprises a turbine body, a steam stop valve, a speed reducer, a speed regulation system, a lubricating oil system, a steam seal air extraction system and a condensate system.

The fluctuation phenomenon of the monitoring parameters such as condenser vacuum and water level, saturated steam pressure, temperature, lubricating oil pressure, lubricating oil temperature and the like is frequently encountered in the use process of the marine power generation turbine, and equipment damage or valve speed closing can be caused when the fluctuation phenomenon is serious.

Typical instabilities in the operation of a power turbine include unstable condenser vacuum, unstable condenser water level, unstable lubricating oil pressure, unstable lubricating oil temperature, unstable steam parameters of the turbine entering and exiting steam, unstable seal pressure and unstable turbine rotational speed.

The power generation steam turbine is provided with a monitoring system, and in the monitoring system, monitoring parameters directly related to the running state of the power generation steam turbine mainly comprise condenser vacuum, front-nozzle steam pressure, front-nozzle steam temperature, rear-nozzle steam pressure, rear-nozzle steam temperature, rear-turbine bypass valve steam pressure, rear-turbine bypass valve steam temperature, turbine gear unit lubricating oil pressure, turbine gear unit lubricating oil temperature, front-rear pressure difference of a turbine lubricating oil filter, power generation steam turbine rotating speed, condenser water level, condenser vacuum degree, 1# to 3# nozzle rear-steam pressure (the power generation steam turbine unit in the embodiment comprises three nozzles), speed valve front-closing steam temperature, nozzle front-steam pressure, speed valve rear-closing steam temperature, turbine unit adjusting stage steam pressure, turbine unit adjusting stage steam temperature, gland sealing pressure and the like. The parameters are the evaluation indexes.

The above parameters are analyzed and a plurality of evaluation indexes related to the running stability of the equipment are selected from the above parameters:

(1) Because the certain type of power generation steam turbine unit adopts saturated steam, the invention has corresponding relation with the steam pressure of corresponding measuring points, and in quantitative analysis, the influence on the running stability state of the power generation steam turbine can be completely converted into the steam pressure, thus the invention does not include the selection range of parameters in analysis.

(2) The parameters such as the front steam pressure of the nozzle of the turbine unit, the rear steam pressure of the bypass valve of the turbine unit and the like have a fixed correlation, if the parameters are selected too much, the redundancy of the parameters is easy to be caused, the calculated amount and the calculated time are increased, and the front steam pressure of the nozzle is an important regulation parameter in the operation process of the power generation turbine, so that only the front steam pressure of the nozzle is reserved in the parameters.

(3) The lubricating oil temperature of the turbine gear unit and the front-rear pressure difference of the turbine lubricating oil filter are fixed relative parameters, and in the operation process of the steam power generation system, the fluctuation range of the lubricating oil temperature of the turbine gear unit is larger than that of the front-rear pressure difference of the turbine lubricating oil filter, so that the lubricating oil temperature of the turbine gear unit is only reserved, and the lubricating oil temperature is called as the lubricating oil temperature in the subsequent calculation.

(4) In the system operation process, the consistency of the steam pressure fluctuation range behind the nozzles 1# to 3# is weaker than that of the operation change trend, so that the three evaluation indexes in the evaluation of the operation state evaluation indexes of the power generation steam turbine are all required to be reserved.

The final evaluation indexes which can be used as the evaluation indexes of the running stability state of the power generation turbine in the embodiment are 10 evaluation indexes, namely the turbine rotating speed, the condenser vacuum, the steam pressure after the No. 1 nozzle to the No. 3 nozzle, the lubricating oil temperature, the steam pressure after the adjusting stage, the steam pressure before the quick valve closing, the condenser water level and the steam seal pressure.

In this embodiment, the data sets with the length N of the 10 evaluation indexes collected by the monitoring system within a period of time are recorded respectively, and each data set is formed into an initial time sequence of the evaluation index operationIn this embodiment, the number of data points collected by each evaluation index is the same, and the number of data points collected by each evaluation index is 940, and s is used as a unit.

S2, decomposing the initial time sequence by adopting an MREMD decomposition method to obtain 7 IMF components and a final residual component。

In this embodiment, step S2 includes the following steps:

（1）

s22, calculating the average value between adjacent extreme points in the new sequence, wherein the average value point obtained in the new sequence is the initial time sequenceMean point sequence>Wherein k is->The number of the mean value points is k+1, and the +.f is obtained through 3 times spline interpolation>Signal mean sequence>Make->Minus->Obtain 1 order signal component->The expression is shown as follows:

（2）

s23, calculating the obtained signal componentAs an original signal, the above steps S21 and S22 are repeated for iterative calculation until the signal component after r iterations +. >And the termination condition is satisfied, and the judgment formula of the calculation termination condition is as follows:

（3）

in the above-mentioned method, the step of,for normalizing standard deviation, ++>And->For the signal component h after the r-1 th iteration respectively _1,r-1 The mean point sequence of (t) and the signal component after the r-th iteration +.>The standard deviation of the mean point sequence of (2), P is the conditional probability; definition according to sigma principle->For the initial time sequence->A sense value of the absolute value of the middle extreme point; />And->Respectively, an initial time sequence and a first time sequenceTime series mean point after r iterations and +.>Ratio sequence of (2); said->、/>And->The formulas of (a) are respectively shown in the following formulas;

（4）

（5）

wherein n is the number of resolvable IMF components, and n is 7 in the embodiment as shown in FIG. 2;as the final residual component;

S25, after the decomposition, the initial time sequenceCan be expressed as the 7 IMF components and the final residual component->The sum is as follows:

（6）。

in the embodiment, the MREMD decomposition methods of formulas (1) - (6) are adopted to decompose the time sequence of the rotation speed of the power generation turbine within 0-940S respectively to obtain 7 IMF components and one residual componentThe results are shown in FIG. 2.

S3, setting a screening threshold to screen the IMF component, and setting an I meeting the screening thresholdPerforming wavelet threshold denoising treatment on the MF component to obtain an IMF component subjected to wavelet denoising, and mixing the IMF component with the IMF component which does not meet the screening threshold and the residual componentAnd (3) recombining to obtain a new IMF component matrix, and processing the new IMF component matrix to obtain a final predicted trend item.

In this embodiment, the step S3 includes the following steps:

s32, screening 7 IMF components by using the screening threshold, and dividing the IMF components into IMF components which do not meet the screening threshold and IMF components which meet the screening threshold;

In the embodiment, the above formulas (7) - (8) are adopted to calculate IMF components and initial signals of each order after the rotation speed of the power generation turbine is decomposedCorrelation coefficient of->Root mean square error->Correlation coefficient screening threshold ∈ ->And root mean square error screening threshold->Wherein the calculated IMF components of each order are combined with the initial signal +.>Correlation coefficient of->Correlation coefficient screening threshold ∈ ->As shown in fig. 3. IMF components of each order and initial signals after the rotation speed of the power generation turbine is decomposedRoot mean square error>Root mean square error screening threshold- >As shown in fig. 4.

From which it can be derived that the formula (9) satisfies the requirementComponent sum->The remaining 5 IMF components do not meet the requirements.

For a pair ofAnd->The components are subjected to wavelet threshold noise reduction processing, the original signals are decomposed into a plurality of layers of approximate coefficients and detail coefficients by utilizing wavelet transformation, and noise information after wavelet transformation is mainly concentrated in the detail coefficients with smaller absolute values, so that the purpose of removing noise can be achieved by setting the detail coefficients with absolute values smaller than a specified threshold to be 0 and reconstructing the residual wavelet coefficients (namely the approximate coefficients obtained by decomposition and the reserved detail coefficients) back to the original signals by utilizing a wavelet inverse transformation method.

The threshold selection and noise reduction of the present embodiment includes the following steps:

estimating noise of each layer in IMF components meeting screening conditions by wavelet detail coefficients of each layer while utilizing coefficientsReducing the threshold value of the detail coefficient layer by layer, thereby preserving as much of the real signal contained in the high frequency component as possible; and carrying out threshold processing by using the composite threshold function to obtain the IMF component after wavelet noise reduction.

And (3) setting a noise threshold value of the j-th layer wavelet detail coefficient and a composite threshold function in the formula (11) by using the formula (10) to perform threshold processing. The signal to noise ratio and root mean square error of the denoised and non-denoised components are shown in table 1.

In the present embodiment, the noise reduction processing will be performedAnd->Component replacement +.>And->And after the components are accumulated with the remaining 5 IMF components which are not denoised, obtaining a historical operation time sequence of the rotating speed of the power generation turbine after denoise treatment. In this embodiment, the signal-to-noise ratio of the time series of the power turbine rotational speed after noise reduction relative to the time series before noise reduction is 81.4374, and the root mean square error is 0.0025. The residual component, the IMF component which does not meet the screening threshold and the IMF component which is subjected to noise reduction treatment are processed>Component sum->The components are reconstructed into a new IMF component matrix.

The same calculation method is used for calculating other 9 evaluation indexes to respectively obtain n IMF components subjected to wavelet threshold noise reduction treatmentAnd residual component->The new IMF component matrix is formed, and will not be described in detail herein.

The invention provides a non-stationary time sequence single parameter prediction model based on MREMD decomposition and wavelet threshold noise reduction, which improves the prediction precision of the model; in the definition of noise IMF components of each IMF component after MREMD decomposition, a subjective IMF component definition method for removing only the first-order IMF component to achieve the purpose of noise reduction is abandoned, and two indexes of a correlation coefficient and a root mean square error are introduced to achieve quantitative definition of the IMF component, so that the calculated error is reduced.

In this embodiment, feature extraction and dimension reduction are performed on the new IMF component matrix by using a singular value decomposition method to obtain a singular value component matrix, including the following steps:

s341, n pieces of processed wavelet threshold noise reductionComponent and residual component->Removing the mean value and then composing the matrix in the form of column vector>And calculates its covariance matrix +.>The expressions of the two matrices are respectively as follows:

（15）

in the method, in the process of the invention,a diagonal matrix consisting of singular values of a covariance matrix B;for n signals->In this embodiment n is 7; wherein->The characteristic values corresponding to covariance matrix B are arranged in descending order; the matrix U is a matrix formed by eigenvectors corresponding to eigenvalues of the covariance matrix B;

s342, combining the matrix U with the matrixMiddle zero bitVector corresponding to the sign value is removed and reconstructed into K singular value components +.>The expression is as follows:

（16）

in the method, in the process of the invention,is the kth singular value component.

In the embodiment, a new IMF component matrix is formed by the rotating speed operation time sequence of the power generation turbine, singular value decomposition is carried out on the IMF component matrix by utilizing formulas (15) - (16), and 8 singular value components are obtained, namelyAs shown in fig. 5.

In this embodiment, the permutation entropy algorithm is used to calculate and obtain permutation entropy values of each singular value component matrix, including the following steps:

S351, reconstructing the phase by using a delay phase space reconstruction methodReconstructing the 1 st singular value component of the plurality of singular value components>Reconstruction into matrix->The expression is:

（17）

in the method, in the process of the invention,the number of sequences after reconstruction, τ is the delay time and m is the embedding dimension.

In the embodiment, during phase space reconstruction, the mutual information method is used to select the optimal delay time, and the mutual information value of each singular value component of the rotation speed of the power generation turbine is within the delay time range of 1-100 secondsThe variation is shown in fig. 6.

Calculating the optimal embedding dimension by adopting a pseudo-neighbor method, wherein the embedding dimension m is within the dimension variation range of 1-10, and the pseudo-neighbor rate of each singular value componentThe variation with embedding dimension m is shown in fig. 7.

S352, arranging each row of elements of the reconstruction matrix in an ascending order to obtain index values of each row of elements, and then forming index vectors by the index values of each row of elements to form an index matrix, wherein the arrangement modes of each row of elements in the index matrix are commonThe possibility that the number of occurrences of the b-th arrangement in the index matrix is counted as +.>Calculate the probability of its occurrence +.>The expression is:

（18）

wherein, the value range of b is (1, K), K is the number of singular value components.

S353, calculating the arrangement entropy of each singular value component, and arranging the entropyThe expression of (2) is as follows:

（19）

in this embodiment, the first extreme point of the time-dependent curve of the mutual information value in fig. 5 is taken as the optimal delay time of each singular value component of the rotation speed of the power generation turbineThe method comprises the steps of carrying out a first treatment on the surface of the Taking the embedding dimension when the pseudo-neighbor rate is less than 5% in FIG. 6 as the optimal embedding dimension of each singular value component +.>The method comprises the steps of carrying out a first treatment on the surface of the Then substituting it into formula (17), and calculating the arrangement entropy value PE of each component by formula (18) and formula (19), the specific values of which are shown in table 2 below:

in this embodiment, the step S36 includes the following steps:

s361, dividing the new IMF component matrix into two types according to the arrangement entropy value by using a K-means clustering algorithm, and randomly selecting the midpoint of the arrangement entropy value sequenceAs a category center, and select and +.>Furthest>The points are another category center, the obtained average value of the two categories of points is used as the new category center to be classified, and then a certain point in the sequence of the arranged entropy values is selected again at will>The calculation is performed as a class center until no more changes occur in both class centers. Finally, the lower class and the higher class of average entropy values are obtained. Wherein the remaining object points x in the sequence of calculated permutation entropy values are respectively associated with +. >Andthe distance of (2) can be +.>Instead, the expression is as follows:

（20）

in the above formula, i=1, …, K is the number of singular value components.

S362, screening out singular value components corresponding to the class with lower average entropy value in the two classes of entropy value points, and reconstructing the singular value components into trend items of the original signalsThe expression is as follows:

（21）

in the method, in the process of the invention, taking a matrix with the mean value expanded for IMF components corresponding to entropy points selected through cluster analysis, and replacing the mean value corresponding to unselected components with a zero vector; />For the screened singular value component matrix, the unselected components are replaced by zero vectors; />To->The non-selected components of the matrix are corresponding to the bitsThe eigenvector is replaced with the zero vector matrix.

In the embodiment, K-means clustering is carried out on the arrangement entropy values of each component of the rotation speed of the power generation turbine by utilizing formulas (20) and (21), and finally selection is carried outAnd->And carrying out trend item reconstruction on the three components to obtain an original trend item, wherein an operation trend chart of the original trend item of the rotation speed of the power generation turbine is shown in fig. 8.

S363, for the trend itemPerforming stability test, and performing differential stabilization treatment on the time sequence with non-stable test result for f times to obtain trend item ++after differential stabilization treatment >The expression of which is shown below,

（22）

in the method, in the process of the invention,the trend item after the differential stabilization treatment is adopted; l is a delay operator; d is the delay order; />Representing the f-time difference.

S364, stabilizing the trend itemPerforming AIC criterion grading, wherein the expression of the AIC criterion grading is as follows:

（23）

in the method, in the process of the invention,the maximum likelihood estimation value is the standard deviation of the model sequence residual error, N is the length of the parameter sequence, r is the number of independent parameters of the model, s is the order of the autoregressive coefficient polynomial, and q is the order of the moving average coefficient polynomial.

S365, obtaining an autoregressive coefficient polynomial order S and a moving average coefficient polynomial order q of the model through AIC test, and fitting an ARIMA (S, f, q) model, wherein the expression is as follows:

（24）

in the method, in the process of the invention,is an AR model coefficient; />Coefficients for an ARIMA model; />Is a white noise sequence; c is a constant; />Is->Is a component of the group.

In this example, the trend term shown in FIG. 8 was checked for stationarity, and Phillips-Perron stationarity was used in this example. The test statistic is suitable for the stationarity test of the heteroscedastic occasion and obeys the limit distribution of the corresponding ADF test statistic. The result of the stability test is a non-stable sequence, and the non-stable sequence is converted into a stable sequence after 1 differential calculation; then, performing a red pool information criterion (AIC) on the sequence subjected to differential stabilization to obtain an autoregressive coefficient polynomial order s and a moving average coefficient polynomial order q which are respectively 3 and 1; after 1 contrast score calculation and LB statistic test, an ARIMA (3, 1, 1) model as shown in formula (24) was fitted.

S366, carrying out LB (Ljung-Box) statistic test on the residual sequence of the ARIMA model, wherein the LB statistic test is represented by the following formula:

（25）

wherein N is an integer and is [1, c ]; c is the residual number;Ljung-Box statistics values for residual sequences; />An estimated value of the model residual error; />Is distributed in a chi-square mode.

In this embodiment, the operation trend of the power generation turbine is predicted at 550s, the operation trends from 400s to 550s are taken as training sets, the difference stabilization processing is performed on the original trend item by using the formula (22), the ARIMA (3, 1) model is obtained by fitting the formulas (23) and (24), the operation trend of 60 seconds in the future between 551-610 s is predicted, and the prediction result of the operation trend of the power generation turbine rotating speed is shown in fig. 9. And (3) performing LB (Ljung-Box) statistic test on the residual sequence of the ARIMA model by using a formula of LB statistic test of the formula (25) so as to meet the requirement.

The operation trend of the other 9 evaluation indexes can refer to the above steps, which are not described herein again, and finally the operation trend prediction sequences of the 10 evaluation indexes are obtained.

S4, converting the new IMF component matrix to obtain a parameter time sequence of each evaluation index, obtaining a weighted evaluation parameter sequence of each evaluation index by using an entropy weight method, and obtaining positive and negative ideal solutions of each moment point by the weighted evaluation parameter sequence.

In this embodiment, the step S4 includes the following steps:

s41, sequentially carrying out accumulation reconstruction on the IMF component subjected to noise reduction of each parameter, the component not subjected to noise reduction treatment and the residual component to form a noise-reduced parameter time sequence, converting the noise-reduced parameter time sequence into a fluctuation value sequence according to the setting value of each evaluation index, obtaining the entropy weight of the fluctuation value sequence by using an entropy weight method, and obtaining the weighted evaluation parameter sequence of each evaluation index by using the entropy weight and the fluctuation value sequence. The entropy weights of the evaluation indexes in this embodiment are respectively: the rotation speed of the power generation turbine is 0.0991, and the steam pressure after the nozzles 0.0906,1# to 3# of the vacuum degree of the condenser is 0.0860, 0.0877 and 0.1200 respectively; 0.1106 for the oil temperature, 0.1042 for the vapor pressure after the conditioning stage, 0.0857 for the vapor pressure after the nozzle, 0.1184 for the condenser water level, 0.0977 for the vapor seal pressure.

The invention combines the improved TOPSIS method with the existing B-type gray correlation degree, and converts the noise-reduced parameter time sequence into the fluctuation value sequence according to the setting value of each evaluation index, so that the positive and negative ideal solutions of each evaluation index of the equipment can be determined, and the problem that the negative ideal solutions cannot be obtained when the existing TOPSIS method is applied to the stability evaluation of the equipment containing non-stable parameters is solved.

（12）

in the method, in the process of the invention,for positive understanding, add>Is negative idealSolving; />，/>N represents the number of evaluation indexes, in this embodiment, N is 10, and the values of positive and negative ideal solutions of the weighted evaluation parameter sequences of the evaluation indexes are finally obtained, and the values of N are 940, and the results are shown in fig. 10 and 11, respectively.

S5, respectively calculating the distance relation value between the fluctuation value of each evaluation index and the positive and negative ideal solutions, wherein the expression is as follows:

（13）

wherein i=1, 2, …, N, the number of N in this embodiment is 940, and the initial time sequenceThe number of N is the same; j=1, 2, …, n, n represents the number of evaluation indexes, in this embodiment, n has a value of 10, < >>For the distance relation value between the fluctuation value of the j-th evaluation index and the positive ideal solution,/for the j-th evaluation index>And the relationship value of the distance between the fluctuation value of the j-th evaluation index and the negative ideal solution.

In the method, an improved distance relation value formula is provided when the positive and negative ideal solutions in the improved TOPSIS method are calculated, and the correlation between the fluctuation value of each evaluation index and the positive and negative ideal solutions is calculated through the formula; in addition, as the distance between the two is smaller, the evaluation index is unstable, the weight occupied by the evaluation index in the stability of the whole equipment is larger, namely the weight of the evaluation index and the corresponding distance relation of the evaluation index are in inverse proportion relation, the negative sign is added in the formula of the distance relation value, so that the real situation can be better reflected, and the obtained result is more accurate.

S6, respectively calculating gray correlation degrees of fluctuation values and positive and negative ideal solutions of all evaluation indexes by adopting a B-type gray correlation analysis method, wherein the method comprises the following steps of:

s61, calculating the overall displacement difference between the fluctuation value of the jth evaluation index at the ith moment and the positive ideal solutionThe expression is as follows:

（26）；

n is the number of data points collected by each evaluation index, and in the embodiment, 940 data points are taken.

S62, respectively calculating the volatility of the j-th evaluation index and the first-order slope difference of the positive ideal solution, wherein the first-order slope difference of the j-th evaluation indexThe expression of (2) is as follows:

（27）

in the above-mentioned method, the step of,a weighted evaluation parameter indicating the (i+1) th time point in the jth evaluation index,/for the (j) th time point>A weighted evaluation parameter indicating an i-th point in time in the j-th evaluation index;

first order slope difference of positive ideal solutionThe expression of (2) is as follows:

（28）

in the above-mentioned method, the step of,representing the maximum forward fluctuation value corresponding to the (i+1) th moment point of n different evaluation indexes; />Representing the maximum forward fluctuation value corresponding to the n different evaluation indexes at the ith moment;

s63, then according to the first-order slope differenceComputing a first order overall slope difference from a first order slope difference of a positive ideal solutionThe expression is as follows:

（29）

Where i=2, 3, …, N, j=1, 2, …, n+1.

S64, calculating a second-order slope difference between the weighted evaluation parameter and the positive ideal solution. Wherein the second order slope difference of the evaluation parameter is weightedThe expression is as follows:

（30）；

second order slope difference of positive ideal solutionThe expression of (2) is as follows:

（31）

s65, according to the second-order slope differenceAnd second order slope difference of ideal>Calculating the second order overall slope difference +.>The expression is as follows:

（32）；

s66, calculating the change rate difference between the weighted evaluation parameter matrix sequence and the positive ideal solutionThe expression is as follows: />

（33）；

Wherein i=3, 4, …, N; j=1, 2, …, n+1;is->Average value of N points in the sample evaluation indexes; />Is the average value of N points in the positive ideal solution;

s67, calculating the association degree of the volatility of the j-th evaluation index and the similarity degree of the positive ideal solution slopeThe expression is shown as follows:

（34）；

s68, calculating the association degree of the volatility of the j-th evaluation index and the similarity of the positive ideal solution change rateThe expression is shown as follows:

（35）

s69, at the association degreeAnd->On the basis of (a) calculating gray correlation degree +_of the j-th evaluation index with respect to the positive ideal solution>The expression is shown as follows:

（36）

Wherein j=1, 2, …, n; omega is an adjusting coefficient, and the value range is (0, 1);

according to the method, the gray correlation degree between the fluctuation value of each evaluation index in the weighted evaluation parameter matrix and the negative ideal solution is calculated。

And S7, carrying out dimensionless treatment on the distance relation value and the gray association degree respectively, and fusing the gray association degree subjected to the dimensionless treatment and the distance relation value into a correlation weight of each evaluation index and the integral running stability state of the equipment.

In this embodiment, the distance relation value is subjected to dimensionless treatmentAnd->As shown in table 3 below.

In this embodiment, the gray correlation degree is obtained after dimensionless treatmentAnd->As shown in table 4 below. />

In this embodiment, after the dimensionless number is obtainedAnd->On the basis of the above, data fusion processing is carried out to obtain a correlation weight value of the volatility of the jth evaluation index and the overall operation stability of the equipment>The expression is as follows:

（37）

wherein, the value range of gamma is [0,1], and the magnitude of gamma determines the curve position relation and geometric shape relation pairIn the present example, the value of γ is 0.5, tableThe distance relation value and the weight value of the type B gray association degree are the same.

In this embodiment, the final weight of each evaluation index is shown below, where the correlation weight of the rotational speed of the power generation turbine is 0.3357, the correlation weight of the vacuum degree of the condenser is 0.2893,1# to 3# of the correlation weights of the steam pressure after the nozzle is 0.3174, 0.3186 and 0.3449, the correlation weight of the oil temperature is 0.5499, the correlation weight of the steam pressure after the adjustment stage is 0.3296, the correlation weight of the steam pressure after the nozzle is 0.3173, the correlation weight of the water level of the condenser is 0.3405, and the correlation weight of the seal pressure is 0.3199.

S8, on the basis of the predicted trend item, establishing an operation stability state scoring formula of each evaluation index, obtaining the scoring value of each evaluation index at each moment point in a predicted time period, and establishing an operation scoring formula of the predicted trend item of the power generation turbine according to the scoring value and the correlation weight of each evaluation index, so as to obtain the operation scoring of the predicted trend item of the whole equipment.

In this embodiment, the expressions of the 10 evaluation index running stability state scoring formulas can be expressed as follows:

（14）

wherein j=1, 2, …, n;for the weighted evaluation parameter of the jth evaluation index in the ith time point of the predicted interval, the predicted interval in this embodiment is 551-610 s, < - >And->And respectively representing the positive maximum fluctuation value and the negative maximum fluctuation value of each evaluation index in the operation range of the power generation turbine.

In this example, each evaluationCorresponding to the indexAnd->The specific numerical value of the current operation evaluation index is determined according to actual operation experience, and is obtained through percentage conversion through setting values of all evaluation indexes, wherein the specific process is the same as that of the current operation evaluation index through percentage conversion into fluctuation value data.

As can be seen from the above-mentioned piecewise function (14), whenWhen the error of (2) is within a maximum fluctuation value of 0.5 times, the corresponding stability state score +.>100 points indicates that the evaluation index is very stable at the ith moment of the prediction interval; when->When the error of (2) is within the range of 0.5-1 times of the maximum fluctuation value, the corresponding stability state score +.>60-100 minutes, which indicates that the evaluation index is relatively stable at the ith moment point of the prediction interval; when->When the error of (2) is outside the maximum fluctuation range of 1 times, the corresponding stability state score +.>And the score is less than 60 minutes, which means that the evaluation index is unstable at the ith moment of the prediction interval, at the moment, the software system can give an alarm in time, and related technicians can adjust related parameters in time according to experience, so that the score is restored to more than 60 minutes. Stability status score of the individual evaluation index +. >The value of (2) is not less than 0.

In this embodiment, the stability status score using 10 evaluation indexesAnd the corresponding correlation weight +.>Establishing a real-time operation evaluation score formula of the power generation steam turbine, wherein the expression is as follows:

（38）

wherein Score is an operation evaluation Score of a predicted trend item of the power generation turbine,is the correlation weight of the fluctuation of the j-th evaluation index and the running stability of the power generation turbine, < ->And (3) scoring the running stability state of the jth evaluation index at the ith moment point of the prediction interval, wherein n is 10.

And S9, setting a total reference index, and comparing the operation evaluation score of the predicted trend item of the power generation turbine with the total reference index to judge whether adjustment of each evaluation index is needed, so that the operation evaluation score of the predicted trend item of the power generation turbine is higher than the total reference index finally.

In this embodiment, the step S9 includes the following steps:

s91, calculating operation evaluation scores of predicted trend items of the power generation steam turbine, and setting total reference indexes, wherein the total reference indexes are 80 minutes in the embodiment;

s92, comparing the operation evaluation score of the predicted trend item of the power generation turbine with the total standard index 80;

And S93, when the operation evaluation score of the predicted trend item of the power generation turbine is greater than or equal to the total reference index, the power generation turbine is well operated in a future predicted time period. In the present embodiment, as shown in fig. 12, the score of the operation evaluation score of the power turbine predicting the trend term between 551s to 610s is greater than 80 points, showing that the power turbine is operating well in the future predicted time period.

When the operation evaluation score of the predicted trend item of the power generation turbine is less than 80 minutes, the power generation turbine is judged to be unstable in operation, at the moment, the software system can give an alarm in time, and review of lower evaluation index scores is returned, and related technicians adjust related parameters in time according to experience, so that the overall real-time operation evaluation score of the power generation turbine is not less than 80 minutes, and the normal operation of the power generation turbine is ensured.

In this embodiment, the software system is Matlab, and each calculation and alarm function is implemented based on the software. Of course other similar software may be used. The operation condition of the power generation turbine in the operation process of the thermodynamic system can be intuitively observed according to the stability state evaluation score diagram of the power generation turbine in a period of time in the future, which is helpful for operation operators to observe the operation state of the condenser in time.

And the software can pre-warn the score of any evaluation index at a certain moment in a future prediction interval, and when the score is less than 60 minutes, relevant technicians timely adjust relevant parameters according to experience to restore the score to more than 60 minutes, so that the safe operation of the software is facilitated. Meanwhile, when the operation evaluation score of the predicted trend item of the power generation turbine is lower than 80, the software system can also give an alarm in time, so that the operation evaluation score of the final predicted trend item is not lower than 80, and the normal operation of the power generation turbine is ensured.

Claims

1. The equipment stability state prediction and evaluation method based on the monitoring parameters is specially used for the power generation steam turbine of the steam power ship and is characterized by comprising the following steps of:

s2, decomposing the initial time sequence by adopting an MREMD decomposition method to obtain n IMF components and a final residual component R _n ；

S3, setting a screening threshold value to screen the IMF components, carrying out wavelet threshold value noise reduction processing on the IMF components meeting the screening threshold value, recombining to obtain a new IMF component matrix, and processing the new IMF component matrix to obtain a final prediction trend item, wherein the method comprises the following steps:

s31, according to the IMF components and the initial time sequence x ₀ Correlation coefficient corr of (t) _i And root mean square error rmse _i Obtaining the screening threshold value;

s36, dividing a new IMF component matrix into two types according to the arrangement entropy value by utilizing a K-means clustering algorithm, taking a type of component with a lower entropy value, reconstructing the type of component into an original trend item, carrying out differential stabilization treatment on the original trend item, obtaining an ARIMA model by fitting, and predicting the original trend item by the ARIMA model to obtain a final predicted trend item;

S4, obtaining a parameter time sequence of each evaluation index by converting the new IMF component matrix, obtaining a weighted evaluation parameter sequence of each evaluation index by using an entropy weight method, and obtaining positive and negative ideal solutions of each moment point by the weighted evaluation parameter sequence, wherein the method comprises the following steps of:

wherein F is ⁺ To positively understand the sequence, F ^- Is a negative ideal solution sequence; n represents the number of evaluation indexes;

S5, calculating a distance relation value between the fluctuation value of each evaluation index and the positive and negative ideal solution respectively, wherein the expression of the distance relation between the fluctuation value of each evaluation index sequence and the positive and negative ideal solution is as follows:

wherein i=1, 2, …, N; j=1, 2, …, n, n represents the number of evaluation indexes,for the distance relation value between the fluctuation value of the j-th evaluation index and the positive ideal solution,/for the j-th evaluation index>The relationship value of the distance between the fluctuation value of the j-th evaluation index and the negative ideal solution;

s8, on the basis of the predicted trend item, establishing an operation stability state scoring formula of each evaluation index, obtaining the scoring value of each evaluation index at each moment point in a predicted time period, and establishing an operation evaluation scoring formula of the predicted trend item of the whole equipment according to the scoring value and the correlation weight of each evaluation index, wherein the expression of the operation stability state scoring formula of each evaluation index is as follows:

Note that when Score _j <At 0, directly define Score _j ＝0

Wherein j=1, 2, …, n; f (f) _ij For the weighted evaluation parameter, delta, of the jth evaluation index at the ith point in time _j，1 And delta _j，2 Respectively representing a positive maximum fluctuation value and a negative maximum fluctuation value of each evaluation index in the equipment operation range;

2. The method for predicting and evaluating the stability status of a device based on monitoring parameters according to claim 1, wherein the step S2 comprises the following steps:

s21, an initial time sequence x of each evaluation index by using an AR prediction model ₀ (t) predicting the left and right endpoints to obtain a new sequence so that x ₀ The end point of (t) is between two extreme points of the extended new sequence, wherein the AR prediction model expression is as follows:

X _t ＝φ ₀ +φ ₁ X _t-1 +φ ₁ X _t-1 +…+φ _p X _t-p +μ _t ；

in the above, x _t The predicted point is the t moment; x is x _t-1 ，x _t-2 ，…，x _t-p The actual values of the parameters at the time t-1, t-2, t-p and phi are respectively shown in the specification ₀ 、φ ₁ 、…、φ _p P+1 real numbers, μ _t (t=p+1, p+2, & gtis, N) white noise sequence with zero mean;

s22, calculating the average value between adjacent extreme points in the new sequence, wherein the average value point obtained in the new sequence is the initial time sequence x ₀ (t) sequence of mean pointsWherein k is x ₀ The number of mean value points of (t) and the number of extreme value points are k+1, and x is obtained through 3 times of spline interpolation ₀ Signal mean sequence m of (t) _1，0 (t) let x be ₀ (t) subtracting m _1，0 (t) obtaining the 1 st order signal component h _1，0 (t) the expression of which is as follows:

h _1,0 (t)＝x ₀ (t)-m _1,0 (t)；

s23, calculatingThe resulting signal component h _1，0 (t) repeating the above steps S21 and S22 as the original signal until the signal component h after r iterations _1，r And (t) a termination condition is satisfied, and the judgment formula of the calculation termination condition is as follows:

in the above, sigma ^* To normalize standard deviation, sigma _r-1 Sum sigma _r For the signal component h after the r-1 th iteration respectively _1,r-1 The mean point sequence of (t) and the signal component h after the r-th iteration _1，r Standard deviation of the mean point sequence of (t), P being a conditional probability; definition according to sigma principleFor an initial time sequence x ₀ A sense value of the absolute value of the extreme point in (t); h is a ₀ And h _r The mean point of the time series after the initial time series and the r iteration and +.>Ratio sequence of (2); the h is ₀ 、h _r And->The formulas of (a) are respectively shown in the following formulas;

wherein: x is X _z The z extreme point of the new sequence after extension;for an initial time sequence x ₀ K mean points of (t); />For the time series component h after the r-1 th iteration _1，r-1 K mean points of (t);

s24, h _1，r (t) as an eigenmode function IMF of order 1 ₁ X is the output of (x) ₀ (t) subtracting IMF ₁ Obtaining component R ₁ And R is taken as ₁ The above steps S21 to S24 are repeated as an original signal until the residual signal becomes a monotone function or no new IMF component is separated, and the residual signal extraction process is as follows:

wherein n is the number of resolvable IMF components; r is R _n As the final residual component;

s25, after the decomposition, the initial time sequence x ₀ (t) can be expressed as the n IMF components and the final residual component R _n The sum is as follows:

3. the method for predicting and evaluating the stability status of a device based on monitoring parameters according to claim 2, wherein the step S31 comprises the following steps:

s311, calculating the IMF components of each order and an initial time sequence x ₀ Correlation coefficient corr of (t) _i And root mean square error rmes _i The expression is shown as follows:

in the middle ofAnd->Respectively initial time series x ₀ (t) and the average value of IMF components of each order, wherein N is the sequence length of each parameter, and N is the number of IMF components;

S312, calculating a screening threshold value, wherein the screening threshold value comprises a correlation coefficient threshold value lambda _IMF And a mean square error threshold delta _IMF The expressions are shown below:

λ _IMF ≥λ

δ _IMF ≤6。

4. a method for predicting and evaluating the stability status of a device based on monitored parameters as claimed in claim 3, wherein said step S32 comprises the steps of:

s321, estimating noise of each layer In IMF components meeting screening conditions through wavelet detail coefficients of each layer, and simultaneously reducing thresholds of the detail coefficients layer by utilizing coefficients In (j+1), wherein a selection criterion formula of the thresholds is as follows:

wherein lambda is _j A noise threshold that is a layer j wavelet detail coefficient; d, d _(j) The detail coefficient of the wavelet of the j-th layer; medium is an intermediate value function, and N is the number of data points collected by each evaluation index;

in the method, in the process of the invention,the j-th layer detail coefficient after the noise reduction treatment; sign is a sign function, a.epsilon.0, 1]The adjustment coefficient is a hard threshold function when a=0, and the soft threshold function when a=1.

5. The method for predicting and evaluating the stability status of a device based on a monitored parameter as set forth in claim 4, wherein said step S9 comprises the steps of: