CN109212966B - Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method - Google Patents

Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method Download PDF

Info

Publication number
CN109212966B
CN109212966B CN201810922586.5A CN201810922586A CN109212966B CN 109212966 B CN109212966 B CN 109212966B CN 201810922586 A CN201810922586 A CN 201810922586A CN 109212966 B CN109212966 B CN 109212966B
Authority
CN
China
Prior art keywords
working condition
equation
state
mechanical equipment
estimating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810922586.5A
Other languages
Chinese (zh)
Other versions
CN109212966A (en
Inventor
雷亚国
韩天宇
李乃鹏
王彪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201810922586.5A priority Critical patent/CN109212966B/en
Publication of CN109212966A publication Critical patent/CN109212966A/en
Application granted granted Critical
Publication of CN109212966B publication Critical patent/CN109212966B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

A multi-working condition dynamic benchmarking mechanical equipment residual life prediction method comprises the steps of firstly establishing a mechanical equipment degradation state space model comprising a state transition equation and an observation equation, secondly estimating unknown parameters and signal transformation parameters of the model, estimating the parameters of the state transition equation by using a maximum likelihood estimation method based on training sample failure time data, transforming monitoring signals under different operating conditions to transformation parameters of reference working condition monitoring signals by linear interpolation estimation, estimating the parameters of the observation equation by using the transformed signals, secondly dynamically benchmarking the monitoring signals of test samples under different working conditions, estimating the state values of the test samples by using a particle filtering algorithm, and finally calculating an analytic solution of a probability density function of the residual service life of the test samples; the method can dynamically standardize the monitoring signals under multiple working conditions in real time in the residual life prediction process, and is beneficial to improving the residual life prediction precision of mechanical equipment.

Description

Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method
Technical Field
The invention belongs to the technical field of mechanical equipment health management and residual life prediction, and particularly relates to a multi-working-condition dynamic benchmarking method for predicting the residual life of mechanical equipment.
Background
With the technological progress, mechanical equipment is continuously enlarged, complicated and precise, the operation conditions tend to be complicated due to diversification of the use functions of the equipment, faults are easily caused frequently, the operation safety of the equipment is influenced, huge loss is brought to economic benefits, and even the life safety of people is seriously threatened. Therefore, it is important to perform health management and remaining life prediction for mechanical equipment, and to implement preventive maintenance.
In practical production, the multi-working condition mainly has two influences on the prediction of the residual service life of the mechanical equipment, namely, the degradation rate of the mechanical equipment is changed, and amplitude sudden change of a monitoring signal is caused. The prior art only considers the influence of different working conditions on the degradation rate of mechanical equipment, but neglects the influence on the signal amplitude. Sudden changes in signal amplitude easily result in misjudgment of the degradation state process of mechanical equipment, and further result in reduced accuracy of residual life prediction. Therefore, considering the influence of multiple working conditions on the degradation rate and the signal amplitude, the effective identification of the real degradation information of the mechanical equipment from the external interference is of great importance for improving the accuracy of residual life prediction.
Disclosure of Invention
In order to overcome the defects in the prior art, the invention provides a multi-working-condition dynamic benchmarking residual life prediction method for mechanical equipment, which comprises the steps of dynamically converting monitoring signals under different working conditions into a reference working condition, obtaining state space model parameters and signal conversion parameters by using training sample data through methods such as maximum likelihood estimation and the like, evaluating the state value of a test sample by using a particle filter algorithm, and finally predicting the probability distribution of the residual life of the test sample by considering the future operating condition of the equipment, thereby improving the residual life prediction precision of the mechanical equipment.
In order to achieve the purpose, the invention adopts the technical scheme that:
a multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method comprises the following steps:
1) constructing a state space model of mechanical equipment degradation:
considering the influence of the variable working condition on the degradation rate and the signal amplitude, the following state space model is established,
Figure BDA0001764605570000021
Figure BDA0001764605570000022
wherein, the formula (1) is a state transition equation, the formula (2) is an observation equation, and xkIs tkOf time of dayThe state value, namely the recession ratio of the mechanical equipment, is 0 in a healthy state and is 1 in a complete failure state; y iskIs tkThe observation value at each moment is a monitoring index capable of reflecting the state of the mechanical equipment; p is a radical ofkRepresents tkThe operation condition of the equipment is kept at the moment;
Figure BDA0001764605570000023
is a condition p introduced to describe the effect of different conditions on the degradation rate of mechanical equipmentkCoefficient of working condition, η is degradation rate, its value is constant when the working condition is not changed, delta tk＝tk-tk-1A time interval for monitoring signal acquisition; omegakSubject to a mean of 0 and a variance of
Figure BDA0001764605570000024
Normally distributed state transition noise of (1);
Figure BDA0001764605570000025
is a condition p introduced for describing the influence of different conditions on the amplitude of the monitoring signalkThe coefficients are respectively rewritten as a by selecting a working condition as a reference working conditionBAnd bB(ii) a c is the order of the observation equation, and represents the nonlinear characteristic of the degradation trend; upsilon iskSubject to a mean of 0 and a variance of σ2The normally distributed measurement noise of (a);
2) estimating state space model parameters
Figure BDA0001764605570000026
η、
Figure BDA0001764605570000027
σ2、aB、bBC and signal transformation parameter a'p、b′p：
Assuming that the mechanical equipment has P different working conditions and N training samples, preprocessing the obtained data according to the following form: working condition coefficient composition vector under different working conditions
Figure BDA0001764605570000028
Wherein the content of the first and second substances,
Figure BDA0001764605570000029
is a working condition pkThe working condition coefficient of the lower part; time to failure of training samples
Figure BDA00017646055700000210
Operating condition P ═ P (P)1,p2,...,pN) And monitoring signal Y ═ Y1,y2,...,yN) Wherein
Figure BDA00017646055700000211
And
Figure BDA00017646055700000212
respectively representing the working condition and the monitoring signal of the nth sample, wherein the value range of N is 1-N, KnSubscripts to their time to failure;
2.1) estimating State transition equation parameters
Figure BDA00017646055700000213
η、
Figure BDA00017646055700000214
Estimating parameters of state transition equations using maximum likelihood estimation methods
Figure BDA0001764605570000031
η and
Figure BDA0001764605570000032
the estimation process is as follows:
2.1.1) derived from the maximum likelihood estimation method η and
Figure BDA0001764605570000033
the maximum likelihood estimate of (a) is a function of R,
Figure BDA0001764605570000034
wherein, it is made
Figure BDA0001764605570000035
Is the running time of the nth training sample under the p working condition;
2.1.2) substituting equation (3) into equation (4) to obtain a log-likelihood function only for R, and using multidimensional optimization to maximize the function, thereby obtaining an estimate of R
Figure BDA0001764605570000036
Figure BDA0001764605570000037
Wherein, let Ψn＝(ψ1,n,ψ2,n,...,ψP,n)，
Figure BDA0001764605570000038
2.1.3) will
Figure BDA0001764605570000039
Alternative R, substitution of formula (3), gives η and
Figure BDA00017646055700000310
is estimated as a result of
Figure BDA00017646055700000311
And
Figure BDA00017646055700000312
2.2) estimating Signal transformation parameters a'pAnd b'p：
Selecting the working condition with the longest running time as a reference working condition, and converting the monitoring signals under different working conditions into the reference working condition through a conversion algorithm to obtain the conversion relation between the reference working condition and model parameters of other working conditions; signal conversionParameter a'pAnd b'pThe estimation process of (2) is as follows:
2.2.1) establishing the relation between the reference working condition and the monitoring signal under the working condition p according to the formula (6),
Figure BDA00017646055700000313
wherein the content of the first and second substances,
Figure BDA00017646055700000314
for the nth sample at time tkObserving equation values a 'in reference working condition and working condition p respectively'p＝aB/apAnd b'p＝aB(bB-bp) For the transformation parameters to be estimated, aBAnd bBIs a parameter of the observation equation under the reference operating condition, apAnd bpIs the parameter of the observation equation under the working condition p;
2.2.2) finding out all the moments of each sample under the working condition p, obtaining the linear interpolation of the monitoring signal under the working condition in the reference working condition, and recording as an interpolation signal
Figure BDA0001764605570000041
2.2.3) calculating an interpolation signal according to equation (7)
Figure BDA0001764605570000042
And the monitoring signal after conversion
Figure BDA0001764605570000043
The sum of the squared errors of (a) and (b),
Figure BDA0001764605570000044
wherein omegap,nRepresenting the time index set of the nth sample under the working condition p;
2.2.4) substitution of formula (8) into formula (7), and determination of a 'by one-dimensional optimization estimation'pIs expressed as
Figure BDA0001764605570000045
And b 'is obtained by substituting the compound into formula (8)'pIs estimated as a result of
Figure BDA0001764605570000046
Figure BDA0001764605570000047
Wherein, | Ωp,nIs omegap,nLength of (d);
2.2.5) repeating the steps 2.2.1) to 2.2.4), sequentially establishing the relationship between the reference working condition and the monitoring signals of other working conditions, and solving the estimated values of the transformation parameters under P-1 working conditions except the reference working condition;
2.3) estimating the parameters a of the observation equationB，bBC and σ2：
Estimating the parameters of the observation equation under reference conditions, i.e. aB，bBC and σ2The estimation process is as follows:
2.3.1) smoothing the transformed monitoring signal by a local regression algorithm, and recording the smoothed monitoring signal as
Figure BDA0001764605570000048
Figure BDA0001764605570000049
2.3.2) mixing of x 1,n0 and
Figure BDA00017646055700000410
substituted into equation (9), observe the equation parameter aB，bBAnd σ2Are respectively the estimation results of
Figure BDA00017646055700000411
2.3.3) order
Figure BDA0001764605570000051
Wherein the content of the first and second substances,
Figure BDA0001764605570000052
for the nth sample tkEstimation of the operating state values of the mechanical device at the moment,
Figure BDA0001764605570000053
for the failure time after the nth sample transformation, the observation equation order c is estimated according to equation (11),
Figure BDA0001764605570000054
3) estimating the state value of the prediction sample by using online data:
estimating the state value of the prediction sample by using a particle filtering algorithm, wherein the specific method comprises the following steps:
3.1) initialization: at a starting time t0Generating Ns state particles
Figure BDA0001764605570000055
The weight of the particle is
Figure BDA0001764605570000056
3.2) predicting: obtaining a one-step predicted value for each state particle according to the state transfer function by the formula (12),
Figure BDA0001764605570000057
3.3) updating: obtaining tkNew monitor signal y of a momentkThen, if the mechanical equipment is not operated under the reference working condition, the monitoring signal is switched to the reference working condition according to the formula (13),
Figure BDA0001764605570000058
the particle weight is updated and normalized according to equation (14),
Figure BDA0001764605570000059
3.4) resampling: the state particles are resampled Ns times, each particle resampling following
Figure BDA00017646055700000510
To generate a new particle sequence
Figure BDA00017646055700000511
Calculating median of particles
Figure BDA00017646055700000512
As a result of the estimation of the device state value;
4) and (3) predicting the residual life: assuming that the mechanical device performs a predetermined task and that future operating conditions are available, based on this assumption, a probability density function for the remaining life is determined according to equation (15),
Figure BDA0001764605570000061
wherein λ 1 is a failure threshold,
Figure BDA0001764605570000062
l is the future time, tkIs the current time.
The invention has the beneficial effects that:
the invention respectively considers the influence of variable working conditions on the degradation rate and the signal amplitude, and the factors of the degradation rate and the signal amplitude are respectively introduced into a state transfer equation and an observation equation of a state space model, so that the proposed signal transformation algorithm can realize the dynamic benchmarking process of the monitoring signal under different working conditions, thereby effectively reducing the interference of the working condition change on the residual life prediction precision. The method can effectively represent the degradation condition of the mechanical equipment in the operation process in industrial practice, and improve the residual life prediction precision of the mechanical equipment.
Drawings
FIG. 1 is a flow chart of the method of the present invention.
FIG. 2 shows monitoring signals and operating speeds of four bearings according to the embodiment, and FIG. (a) shows monitoring signals and operating speeds of the bearing 1; graph (b) is the bearing 2 monitoring signal and the running speed; graph (c) shows the monitoring signal and the running speed of the bearing 3; and (d) is a bearing 4 monitoring signal and an operation rotating speed.
FIG. 3 is a comparison of the predicted residual life results of four bearings according to two prediction methods of the embodiment, and FIG. (a) is a comparison of the predicted residual life results of the bearing 1; the graph (b) is a comparison of the predicted results of the residual life of the bearing 2; the graph (c) is a comparison of the predicted results of the residual life of the bearing 3; fig. (d) is a comparison of the predicted residual life of the bearing 4.
FIG. 4 is a comparison of the predicted error of four bearings under two prediction methods of the example, wherein (a) is a comparison of the mean values; graph (b) is variance comparison.
Detailed Description
The invention is further elucidated with reference to the figures and embodiments.
Referring to fig. 1, a method for predicting the residual life of mechanical equipment with multi-working-condition dynamic benchmarking includes the following steps:
1) constructing a state space model of mechanical equipment degradation:
considering the influence of the variable working condition on the degradation rate and the signal amplitude, the following state space model is established,
Figure BDA0001764605570000071
Figure BDA0001764605570000072
wherein, the formula (1) is a state transition equation, the formula (2) is an observation equation, and xkIs tkThe state value at the moment, namely the recession ratio of the mechanical equipment, is 0 in a healthy state and 1 in a complete failure state; y iskIs tkThe observation value at each moment is a monitoring index capable of reflecting the state of the mechanical equipment;pkrepresents tkThe operation condition of the equipment is kept at the moment;
Figure BDA0001764605570000073
is a condition p introduced to describe the effect of different conditions on the degradation rate of mechanical equipmentkCoefficient of working condition, η is degradation rate, its value is constant when the working condition is not changed, delta tk＝tk-tk-1A time interval for monitoring signal acquisition; omegakSubject to a mean of 0 and a variance of
Figure BDA0001764605570000074
Normally distributed state transition noise of (1);
Figure BDA0001764605570000075
is a condition p introduced for describing the influence of different conditions on the amplitude of the monitoring signalkThe coefficients are respectively rewritten as a by selecting a working condition as a reference working conditionBAnd bB(ii) a c is the order of the observation equation, and represents the nonlinear characteristic of the degradation trend; upsilon iskSubject to a mean of 0 and a variance of σ2The normally distributed measurement noise of (a);
2) estimating state space model parameters
Figure BDA0001764605570000076
η、
Figure BDA0001764605570000077
σ2、aB、bBC and signal transformation parameter a'p、b′p：
Assuming that the mechanical equipment has P different working conditions and N training samples, preprocessing the obtained data according to the following form: working condition coefficient composition vector under different working conditions
Figure BDA0001764605570000078
Wherein the content of the first and second substances,
Figure BDA0001764605570000079
is a working condition pkThe working condition coefficient of the lower part; time to failure of training samples
Figure BDA00017646055700000710
Operating condition P ═ P (P)1,p2,...,pN) And monitoring signal Y ═ Y1,y2,...,yN) Wherein
Figure BDA00017646055700000711
And
Figure BDA00017646055700000712
respectively representing the working condition and the monitoring signal of the nth sample, wherein the value range of N is 1-N, KnSubscripts to their time to failure;
2.1) estimating State transition equation parameters
Figure BDA00017646055700000713
η、
Figure BDA00017646055700000714
Estimating parameters of state transition equations using maximum likelihood estimation methods
Figure BDA00017646055700000715
η and
Figure BDA00017646055700000716
the estimation process is as follows:
2.1.1) derived from the maximum likelihood estimation method η and
Figure BDA0001764605570000081
the maximum likelihood estimate of (a) is a function of R,
Figure BDA0001764605570000082
wherein, it is made
Figure BDA0001764605570000083
Is the running time of the nth training sample under the p working condition;
2.1.2) substituting equation (3) into equation (4) to obtain a log-likelihood function only for R, and using multidimensional optimization to maximize the function, thereby obtaining an estimate of R
Figure BDA0001764605570000084
Figure BDA0001764605570000085
Wherein, let Ψn＝(ψ1,n,ψ2,n,...,ψP,n)，
Figure BDA0001764605570000086
2.1.3) will
Figure BDA0001764605570000087
Alternative R, substitution of formula (3), gives η and
Figure BDA0001764605570000088
is estimated as a result of
Figure BDA0001764605570000089
And
Figure BDA00017646055700000810
2.2) estimating Signal transformation parameters a'pAnd b'p：
Selecting the working condition with the longest running time as a reference working condition, and converting the monitoring signals under different working conditions into the reference working condition through a conversion algorithm to obtain the conversion relation between the reference working condition and model parameters of other working conditions; signal transformation parameter a'pAnd b'pThe estimation process of (2) is as follows:
2.2.1) establishing the relation between the reference working condition and the monitoring signal under the working condition p according to the formula (6),
Figure BDA00017646055700000811
wherein the content of the first and second substances,
Figure BDA00017646055700000812
for the nth sample at time tkObserving equation values a 'in reference working condition and working condition p respectively'p＝aB/apAnd b'p＝aB(bB-bp) For the transformation parameters to be estimated, aBAnd bBIs a parameter of the observation equation under the reference operating condition, apAnd bpIs the parameter of the observation equation under the working condition p;
2.2.2) finding out all the moments of each sample under the working condition p, obtaining the linear interpolation of the monitoring signal under the working condition in the reference working condition, and recording as an interpolation signal
Figure BDA0001764605570000091
2.2.3) calculating an interpolation signal according to equation (7)
Figure BDA0001764605570000092
And the monitoring signal after conversion
Figure BDA0001764605570000093
The sum of the squared errors of (a) and (b),
Figure BDA0001764605570000094
wherein omegap,nRepresenting the time index set of the nth sample under the working condition p;
2.2.4) substitution of formula (8) into formula (7), and determination of a 'by one-dimensional optimization estimation'pIs expressed as
Figure BDA0001764605570000095
And b 'is obtained by substituting the compound into formula (8)'pIs estimated as a result of
Figure BDA0001764605570000096
Figure BDA0001764605570000097
Wherein, | Ωp,nIs omegap,nLength of (d);
2.2.5) repeating the steps 2.2.1) to 2.2.4), sequentially establishing the relationship between the reference working condition and the monitoring signals of other working conditions, and solving the estimated values of the transformation parameters under P-1 working conditions except the reference working condition;
2.3) estimating the parameters a of the observation equationB，bBC and σ2：
Estimating the parameters of the observation equation under reference conditions, i.e. aB，bBC and σ2The estimation process is as follows:
2.3.1) smoothing the transformed monitoring signal by a local regression algorithm, and recording the smoothed monitoring signal as
Figure BDA0001764605570000098
Figure BDA0001764605570000099
2.3.2) mixing of x1,n0 and
Figure BDA00017646055700000910
substituted into equation (9), observe the equation parameter aB，bBAnd σ2Are respectively the estimation results of
Figure BDA00017646055700000911
2.3.3) order
Figure BDA00017646055700000912
Wherein the content of the first and second substances,
Figure BDA0001764605570000101
for the nth sample tkEstimation of the operating state values of the mechanical device at the moment,
Figure BDA0001764605570000102
for the failure time after the nth sample transformation, the observation equation order c is estimated according to equation (11),
Figure BDA0001764605570000103
3) estimating the state value of the prediction sample by using online data:
estimating the state value of the prediction sample by using a particle filtering algorithm, wherein the specific method comprises the following steps:
3.1) initialization: at a starting time t0Generating Ns state particles
Figure BDA0001764605570000104
The weight of the particle is
Figure BDA0001764605570000105
3.2) predicting: obtaining a one-step predicted value for each state particle according to the state transfer function by the formula (12),
Figure BDA0001764605570000106
3.3) updating: obtaining tkNew monitor signal y of a momentkThen, if the mechanical equipment is not operated under the reference working condition, the monitoring signal is switched to the reference working condition according to the formula (13),
Figure BDA0001764605570000107
the particle weight is updated and normalized according to equation (14),
Figure BDA0001764605570000108
3.4) resampling: the state particles are resampled Ns times, each particle resampling following
Figure BDA0001764605570000109
To generate a new particle sequence
Figure BDA00017646055700001010
Calculating median of particles
Figure BDA00017646055700001011
As a result of the estimation of the device state value;
4) and (3) predicting the residual life: assuming that the mechanical device performs a predetermined task and that future operating conditions are available, based on this assumption, a probability density function for the remaining life is determined according to equation (15),
Figure BDA0001764605570000111
wherein λ 1 is a failure threshold,
Figure BDA0001764605570000112
l is the future time, tkIs the current time.
The rolling bearing is used as a key part in mechanical equipment, is widely applied and is easy to break down, so that the residual life prediction is particularly important to be carried out on the rolling bearing, and preventive maintenance can be carried out in the industry according to the residual life prediction result. To further demonstrate the effectiveness of the method of the present invention, a residual life prediction was developed using vibration acceleration signals obtained from an accelerated degradation experiment of a rolling bearing in conjunction with the method of the present invention and compared to a degradation modeling method (denoted as M1) proposed by LINKAN BIAN et al, of Missippi State university, USA. When the root mean square value of the vibration acceleration of the rolling bearing exceeded 2.2g (g is the gravitational acceleration), the bearing was considered to be failed, and the experiment was terminated. As shown in fig. 2, bearing 1 and bearing 2 operate at 2200rpm and 2600rpm, respectively, and bearing 3 and bearing 4 both operate at varying operating conditions. It can be seen from fig. 2 that the change of the operating condition results in a sudden change of the amplitude of the monitoring signal, the signal randomly fluctuates in the initial stage without a degradation trend, and the monitoring signal gradually shows a degradation trend as the operation time increases. The remaining life prediction is started after the initial prediction point. M1 mixed considers the influence of multiple working conditions on the degradation rate and the signal amplitude, M2 is the method provided by the invention, and the influence of the multiple working conditions on the degradation rate and the signal amplitude is respectively modeled.
In practice, one bearing is selected as the test sample, the other three bearings are used as training samples to estimate the model parameters, and the four bearings are used as the test samples in turn. Fig. 3 shows the predicted residual life of four bearings as test samples. In order to quantitatively evaluate the performance of both methods, the mean and variance of the relative error absolute value ARE of each method were calculated according to equation (18).
Figure BDA0001764605570000113
Wherein, TPreIs the predicted time to failure, TActIs the true time to failure. The calculation results are shown in fig. 4. As can be seen from fig. 3 and 4, the residual life prediction result of the method provided by the present invention is more accurate and stable.
The method for predicting the residual life of the mechanical equipment based on the dynamic benchmarking of the multiple working conditions can be suitable for predicting the residual life of various mechanical equipment. In practical application, an implementer can reasonably determine parameters such as a set corresponding to working conditions, a signal smoothing method and the like according to the operating working conditions of various mechanical equipment. The method provided by the invention is beneficial to improving the accuracy of the residual life prediction of the mechanical equipment. It should be noted that modifications and variations to the method described herein are possible without departing from the inventive concept.

Claims (1)

1. A multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method is characterized by comprising the following steps:
1) constructing a state space model of mechanical equipment degradation:
considering the influence of the variable working condition on the degradation rate and the signal amplitude, the following state space model is established,
Figure FDA0002351296730000011
Figure FDA0002351296730000012
wherein, the formula (1) is a state transition equation, the formula (2) is an observation equation, and xkIs tkThe state value at the moment, namely the recession ratio of the mechanical equipment, is 0 in a healthy state and 1 in a complete failure state; y iskIs tkThe observation value at each moment is a monitoring index capable of reflecting the state of the mechanical equipment; p is a radical ofkRepresents tkThe operation condition of the equipment is kept at the moment;
Figure FDA0002351296730000013
is a condition p introduced to describe the effect of different conditions on the degradation rate of mechanical equipmentkCoefficient of working condition, η is degradation rate, its value is constant when the working condition is not changed, delta tk＝tk-tk-1A time interval for monitoring signal acquisition; omegakSubject to a mean of 0 and a variance of
Figure FDA0002351296730000014
Normally distributed state transition noise of (1);
Figure FDA0002351296730000015
is a condition p introduced for describing the influence of different conditions on the amplitude of the monitoring signalkThe coefficients are respectively rewritten as a by selecting a working condition as a reference working conditionBAnd bB(ii) a c is the order of the observation equation, and represents the nonlinear characteristic of the degradation trend; upsilon iskSubject to a mean of 0 and a variance of σ2The normally distributed measurement noise of (a);
2) estimating a shapeState space model parameters
Figure FDA0002351296730000016
η、
Figure FDA0002351296730000017
σ2、aB、bBC and signal transformation parameter a'p、b′p：
Assuming that the mechanical equipment has P different working conditions and N training samples, preprocessing the obtained data according to the following form: working condition coefficient composition vector under different working conditions
Figure FDA0002351296730000018
Time to failure of training samples
Figure FDA0002351296730000019
Operating condition P ═ P (P)1,p2,...,pN) And monitoring signal Y ═ Y1,y2,...,yN) Wherein
Figure FDA00023512967300000110
And
Figure FDA00023512967300000111
respectively representing the working condition and the monitoring signal of the nth sample, wherein the value range of N is 1-N, KnSubscripts to their time to failure;
2.1) estimating State transition equation parameters
Figure FDA0002351296730000021
η、
Figure FDA0002351296730000022
Estimating parameters of state transition equations using maximum likelihood estimation methods
Figure FDA0002351296730000023
η and
Figure FDA0002351296730000024
the estimation process is as follows:
2.1.1) derived from the maximum likelihood estimation method η and
Figure FDA0002351296730000025
the maximum likelihood estimate of (a) is a function of R,
Figure FDA0002351296730000026
wherein, it is made
Figure FDA0002351296730000027
Figure FDA0002351296730000028
Is the running time of the nth training sample under the p working condition;
2.1.2) substituting equation (3) into equation (4) to obtain a log-likelihood function only for R, and using multidimensional optimization to maximize the function, thereby obtaining an estimate of R
Figure FDA0002351296730000029
Figure FDA00023512967300000210
Wherein, let Ψn＝(ψ1,n,ψ2,n,...,ψP,n)，
Figure FDA00023512967300000211
2.1.3) will
Figure FDA00023512967300000212
The substitution of the R group is carried out,substitution of formula (3) to obtain η and
Figure FDA00023512967300000213
is estimated as a result of
Figure FDA00023512967300000214
And
Figure FDA00023512967300000215
2.2) estimating Signal transformation parameters a'pAnd b'p：
Selecting the working condition with the longest running time as a reference working condition, and converting the monitoring signals under different working conditions into the reference working condition through a conversion algorithm to obtain the conversion relation between the reference working condition and model parameters of other working conditions; signal transformation parameter a'pAnd b'pThe estimation process of (2) is as follows:
2.2.1) establishing the relation between the reference working condition and the monitoring signal under the working condition p according to the formula (6),
Figure FDA00023512967300000216
wherein the content of the first and second substances,
Figure FDA0002351296730000031
for the nth sample at time tkObserving equation values a 'in reference working condition and working condition p respectively'p＝aB/apAnd b'p＝aB(bB-bp) For the transformation parameters to be estimated, apAnd bpIs the parameter of the observation equation under the working condition p;
2.2.2) finding out all the moments of each sample under the working condition p, obtaining the linear interpolation of the monitoring signal under the working condition in the reference working condition, and recording as an interpolation signal
Figure FDA0002351296730000032
2.2.3) calculating the interpolation information according to equation (7)Number (C)
Figure FDA0002351296730000033
And the monitoring signal after conversion
Figure FDA0002351296730000034
The sum of the squared errors of (a) and (b),
Figure FDA0002351296730000035
wherein omegap,nRepresenting the time index set of the nth sample under the working condition p;
2.2.4) substitution of formula (8) into formula (7), and determination of a 'by one-dimensional optimization estimation'pIs expressed as
Figure FDA0002351296730000036
And b 'is obtained by substituting the compound into formula (8)'pIs estimated as a result of
Figure FDA0002351296730000037
Figure FDA0002351296730000038
Wherein, | Ωp,nIs omegap,nLength of (d);
2.2.5) repeating the steps 2.2.1) to 2.2.4), sequentially establishing the relationship between the reference working condition and the monitoring signals of other working conditions, and solving the estimated values of the transformation parameters under P-1 working conditions except the reference working condition;
2.3) estimating the parameters a of the observation equationB，bBC and σ2：
Estimating the parameters of the observation equation under reference conditions, i.e. aB，bBC and σ2The estimation process is as follows:
2.3.1) smoothing the transformed monitoring signal by a local regression algorithm, and recording the smoothed monitoring signal as
Figure FDA0002351296730000039
Figure FDA00023512967300000310
2.3.2) mixing of x1,n0 and
Figure FDA00023512967300000311
substituted into equation (9), observe the equation parameter aB，bBAnd σ2Are respectively the estimation results of
Figure FDA0002351296730000041
2.3.3) order
Figure FDA0002351296730000042
Wherein the content of the first and second substances,
Figure FDA0002351296730000043
for the nth sample tkEstimation of the operating state values of the mechanical device at the moment,
Figure FDA0002351296730000044
for the failure time after the nth sample transformation, the observation equation order c is estimated according to equation (11),
Figure FDA0002351296730000045
3) estimating the state value of the prediction sample by using online data:
estimating the state value of the prediction sample by using a particle filtering algorithm, wherein the specific method comprises the following steps:
3.1) initialization: at a starting time t0Generating Ns state particles
Figure FDA0002351296730000046
The weight of the particle is
Figure FDA0002351296730000047
3.2) predicting: obtaining a one-step predicted value for each state particle according to the state transfer function by the formula (12),
Figure FDA0002351296730000048
3.3) updating: obtaining tkNew monitor signal y of a momentkThen, if the mechanical equipment is not operated under the reference working condition, the monitoring signal is switched to the reference working condition according to the formula (13),
Figure FDA0002351296730000049
the particle weight is updated and normalized according to equation (14),
Figure FDA00023512967300000410
3.4) resampling: the state particles are resampled Ns times, each particle resampling following
Figure FDA00023512967300000411
To generate a new particle sequence
Figure FDA00023512967300000412
Calculating median of particles
Figure FDA00023512967300000413
As a result of the estimation of the device state value;
4) and (3) predicting the residual life: assuming that the mechanical device performs a predetermined task and that future operating conditions are available, based on this assumption, a probability density function for the remaining life is determined according to equation (15),
Figure FDA0002351296730000051
wherein λ 1 is a failure threshold,
Figure FDA0002351296730000052
l is the future time, tkIs the current time.
CN201810922586.5A 2018-08-14 2018-08-14 Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method Active CN109212966B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810922586.5A CN109212966B (en) 2018-08-14 2018-08-14 Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810922586.5A CN109212966B (en) 2018-08-14 2018-08-14 Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method

Publications (2)

Publication Number Publication Date
CN109212966A CN109212966A (en) 2019-01-15
CN109212966B true CN109212966B (en) 2020-04-10

Family

ID=64988549

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810922586.5A Active CN109212966B (en) 2018-08-14 2018-08-14 Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method

Country Status (1)

Country Link
CN (1) CN109212966B (en)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110222855B (en) * 2019-06-12 2021-10-12 中国神华能源股份有限公司 Method and device for processing train wheel degradation data and storage medium
CN112444725B (en) * 2019-09-05 2024-01-26 中车株洲电力机车研究所有限公司 Through hole welding spot contrast accelerated life test method
CN110956339A (en) * 2019-12-17 2020-04-03 上海威派格智慧水务股份有限公司 Flow prediction method
CN111143990B (en) * 2019-12-25 2021-11-16 西安交通大学 Sensitive measuring point selection and fusion machine tool milling cutter residual life prediction method
CN111751724A (en) * 2020-06-24 2020-10-09 湖北文理学院 Motor application working condition information monitoring method and device and readable storage medium
CN112518425B (en) * 2020-12-10 2022-10-04 南京航空航天大学 Intelligent machining cutter wear prediction method based on multi-source sample migration reinforcement learning
CN112855467B (en) * 2021-03-22 2022-01-11 西安交通大学 Wind driven generator reference working condition conversion method
CN117056692B (en) * 2023-10-09 2024-01-23 山东芯赛思电子科技有限公司 Aging prediction method for SiC-based motor driving device
CN117291445B (en) * 2023-11-27 2024-02-13 国网安徽省电力有限公司电力科学研究院 Multi-target prediction method based on state transition under comprehensive energy system

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102542155A (en) * 2011-12-05 2012-07-04 北京航空航天大学 Particle filter residual life forecasting method based on accelerated degradation data
CN102721920A (en) * 2012-06-29 2012-10-10 沈阳工业大学 Prediction device and method for remaining life of operating mechanism of circuit breaker
CN104598734A (en) * 2015-01-22 2015-05-06 西安交通大学 Life prediction model of rolling bearing integrated expectation maximization and particle filter
CN105956236A (en) * 2016-04-22 2016-09-21 西安交通大学 Dual-updating four-factor random degeneration model gear life prediction method
CN106934125A (en) * 2017-02-28 2017-07-07 西安交通大学 A kind of exponential model plant equipment method for predicting residual useful life of trapezoidal noise profile

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106876224B (en) * 2017-01-20 2019-06-11 西安交通大学 The forecasting system and method for impact air pressure suffered by breaker of plastic casing arc extinguishing room housing

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102542155A (en) * 2011-12-05 2012-07-04 北京航空航天大学 Particle filter residual life forecasting method based on accelerated degradation data
CN102721920A (en) * 2012-06-29 2012-10-10 沈阳工业大学 Prediction device and method for remaining life of operating mechanism of circuit breaker
CN104598734A (en) * 2015-01-22 2015-05-06 西安交通大学 Life prediction model of rolling bearing integrated expectation maximization and particle filter
CN105956236A (en) * 2016-04-22 2016-09-21 西安交通大学 Dual-updating four-factor random degeneration model gear life prediction method
CN106934125A (en) * 2017-02-28 2017-07-07 西安交通大学 A kind of exponential model plant equipment method for predicting residual useful life of trapezoidal noise profile

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A New Method Based on Stochastic Process Models for Machine Remaining Useful Life Prediction;Yaguo Lei,Naipeng Li,Jing Lin;《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》;20161231;第65卷(第12期);第2671-2684页 *
Machinery health prognostics: A systematic review from data acquisition to RUL prediction;Yaguo Lei，等;《Mechanical Systems and Signal Processing》;20180531;第800-826页 *
一种新的混合智能预测模型及其在故障诊断中的应用;胡桥，等;《西安交通大学学报》;20050930;第39卷(第9期);第928-932页 *

Also Published As

Publication number Publication date
CN109212966A (en) 2019-01-15

Similar Documents

Publication Publication Date Title
CN109212966B (en) Multi-working-condition dynamic benchmarking mechanical equipment residual life prediction method
CN107145645B (en) Method for predicting residual life of non-stationary degradation process with uncertain impact
Singleton et al. Extended Kalman filtering for remaining-useful-life estimation of bearings
Wu et al. Induction machine fault detection using SOM-based RBF neural networks
EP3719601A1 (en) Abnormality determination system, motor control device, and abnormality determination device
CN109829136B (en) Method and system for predicting residual life of degradation equipment with random jump
CN113947017A (en) Method for predicting residual service life of rolling bearing
CN110631849B (en) Online identification method and system for wear state of friction system
EP3649568B1 (en) Method for automatic detection of physical modes in a modal analysis model
CN108492026B (en) Soft measurement method based on integrated orthogonal component optimization regression analysis
CN112525749A (en) Tribology state online identification method based on friction signal recursion characteristic
CN116007940A (en) Rolling bearing residual life prediction method considering matching failure threshold
CN114755010A (en) Rotary machine vibration fault diagnosis method and system
CN114757087A (en) Tool wear prediction method based on dynamic principal component analysis and LSTM
CN113159088B (en) Fault monitoring and diagnosis method based on multi-feature fusion and width learning
CN113469408A (en) Running state trend prediction method and system for phase modulator
CN110657807B (en) Indoor positioning displacement measurement method for detecting discontinuity based on wavelet transformation
CN115291103A (en) Motor fault diagnosis method based on GR-SWPT wavelet packet algorithm embedded with HD-RCF
CN113822565B (en) Method for graded and refined analysis of time-frequency characteristics of fan monitoring data
Li et al. Multi-scale fractal dimension based on morphological covering for gear fault diagnosis
CN114154736A (en) Rolling bearing service life prediction method based on similarity matching optimization theory
CN114970311A (en) Method for establishing remote module life prediction model and life prediction method
CN111160464B (en) Industrial high-order dynamic process soft measurement method based on multi-hidden-layer weighted dynamic model
Gašperin et al. Condition prognosis of mechanical drives based on nonlinear dynamical models
CN113878613B (en) Industrial robot harmonic reducer early fault detection method based on WLCTD and OMA-VMD

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant