CN109740255B - Life prediction method based on time-varying Markov process - Google Patents

Life prediction method based on time-varying Markov process Download PDF

Info

Publication number
CN109740255B
CN109740255B CN201910001112.1A CN201910001112A CN109740255B CN 109740255 B CN109740255 B CN 109740255B CN 201910001112 A CN201910001112 A CN 201910001112A CN 109740255 B CN109740255 B CN 109740255B
Authority
CN
China
Prior art keywords
time
state
probability
varying
formula
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
CN201910001112.1A
Other languages
Chinese (zh)
Other versions
CN109740255A (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.)
Beijing Jiaotong University
Original Assignee
Beijing 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 Beijing Jiaotong University filed Critical Beijing Jiaotong University
Priority to CN201910001112.1A priority Critical patent/CN109740255B/en
Publication of CN109740255A publication Critical patent/CN109740255A/en
Application granted granted Critical
Publication of CN109740255B publication Critical patent/CN109740255B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

The invention provides a life prediction method based on a time-varying Markov process, which comprises the steps of establishing a time-varying Markov model, calculating initial forward probability, backward probability and a conditional probability formula of a state sequence by utilizing an improved forward-backward algorithm, reestimating model parameters, obtaining the expectation of the stay time of equipment in the state according to the determined current running state of the equipment, solving the stay time of the equipment in each state, and calculating the residual life of the stay time of the current state.

Description

Life prediction method based on time-varying Markov process
Technical Field
The invention relates to the field of reliability, in particular to a life prediction method for a rolling bearing of rail transit.
Background
The Security Region (SR) is a quantitative model describing the overall safe and stable operation Region of the system from the perspective of the domain, and the relative relationship between the security Region boundary and the system operation point can provide the quantitative operation Safety margin and the optimal control information of the system under different conditions.
In the prior art, the concept of a security domain is expanded to the field of rail transit, a complete security domain estimation method framework is provided, and beneficial attempts are made for rolling bearings. The research divides the bearing state characteristic variable space into two regions, a safety region and a non-safety region, and directly corresponds the state characteristic data to different categories such as 'normal' and 'fault' according to the state information. However, most bearing performance is a gradual degradation process, and between normal and fault, there are multiple safety sub-domains, and there is some probabilistic transfer relationship between the different safety sub-domains. And with the development of sensor technology, the observation of the operation process of a system or key components in the field of industrial engineering is more and more common, so that a large number of time series are generated, and the time series are often multivariable mutual influence and mutual connection. Therefore, the bearing degradation model based on the safety domain estimation theory is expanded to multiple safety sub-domains, and the automatic positioning of the boundary of the safety domain and the non-safety domain and the critical position of each safety sub-domain is realized by adopting the time sequence fuzzy clustering algorithm.
Disclosure of Invention
In order to achieve the purpose, the reliability analysis and service life prediction method based on the time-varying Markov process mainly comprises four steps: the method comprises the steps of degradation characteristic selection, fuzzy security domain division, time-varying Markov process modeling, reliability evaluation and residual life prediction, and a technical route diagram is shown in figure 1. The invention specifically adopts the following technical scheme:
(1) extracting a full-life operation characteristic vector, segmenting based on a multivariate fuzzy segmentation algorithm, and performing sample time warping by adopting a dynamic time warping algorithm to finish sample data preprocessing;
(2) dividing the states of the components by adopting a fuzzy security domain algorithm;
(3) establishing a time-varying Markov model, and calculating an initial forward probability, a backward probability and a conditional probability formula of a state sequence by using an improved forward-backward algorithm; and reestimating the model parameters; obtaining the expectation of the stay time of the equipment in the state i according to the determined current running state i of the equipment, and solving the stay time of the equipment in the states i +1, i +2, …, N; calculating the dwell time of the current state i as
Figure BDA0001933654120000011
The remaining life of;
firstly, determining the residence time distribution of each state according to a sample set, and estimating model parameters by adopting an improved forward-backward algorithm:
1) when the model M is known as the (pi, A, ST) parameter, the system stays in the state i at the time tIs at interval of T*Probability of (2)
Figure BDA0001933654120000021
When t is 1, the initial forward probability is:
Figure BDA0001933654120000022
Figure BDA0001933654120000023
Figure BDA0001933654120000024
when T is 2,3, …, T, the forward probability recurrence formula is:
Figure BDA0001933654120000025
Figure BDA0001933654120000026
Figure BDA0001933654120000027
2) the residence time of state i at time t is
Figure BDA0001933654120000028
The state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
Figure BDA0001933654120000029
When T is T, the initial backward probability formula is:
Figure BDA00019336541200000210
Figure BDA00019336541200000211
Figure BDA00019336541200000212
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
Figure BDA00019336541200000213
Figure BDA00019336541200000214
based on the forward probability formula and the backward probability formula, the state vector S ═ (S) can be derived1,s2,...,sT)，s1,s2,...,sTThe conditional probability formula of epsilon { i,1 is more than or equal to i and less than or equal to N } is
Figure BDA0001933654120000031
3) Parameter estimation
Obtaining a parameter re-estimation formula in a time-varying semi-Markov process model based on a forward-backward probability formula
Figure BDA0001933654120000033
ζt(i, j, st) for a given model M and a state sequence of S ═ S1,s2,...,st)，s1,s2,...,stWhen e is { i,1 is more than or equal to i and less than or equal to N }, the system is transferred after the retention time of the state i is stProbability to state j.
γt(i, st) given the model M and the state sequence S, the probability that the system stays in state i for a time st at time t,
Figure BDA0001933654120000038
setting etat(i, st) is the probability of the system staying at the state i for the time st at the moment t,
Figure BDA0001933654120000039
by
Figure BDA0001933654120000034
In a clear view of the above, it is known that,
Figure BDA0001933654120000035
Figure BDA0001933654120000036
represents the desired number of transitions from state i to state j while the system remains in state i for a time st
Figure BDA0001933654120000037
Represents the desired number of transitions of the system from state i when the system stays in state i for a time st, and therefore aijThe formula for reevaluation of (st) is
Figure BDA0001933654120000041
The probability of residence time in state i in model M obeys the mean value μiVariance is
Figure BDA0001933654120000042
The distribution of the gaussian component of (a) is,
Figure BDA0001933654120000043
Figure BDA0001933654120000044
4) reliability determination and update
The failure rate function of the device is λ (t), the life distribution function is F (t), the density function is F (t), F (t) ═ F' (t), and the reliability function is r (t), then F (t) + r (t) ═ 1, λ (t) ═ F (t)/r (t) exist.
When the total number of samples is M, M (t) represents the number of samples with faults before the time t, Δ M (t) represents the number of samples with faults in the time interval (t, t + Δ t), and if,
Figure BDA0001933654120000045
to pair
Figure BDA0001933654120000046
With R (t)>0, therefore, λ (t) can be considered as an estimate of the (t, t + Δ t) interval fault condition probability.
L denotes the life cycle of the device, rul (t) denotes the remaining life expectancy from t to system failure at time t when the device is operating normally,
Figure BDA0001933654120000047
the known apparatus is at time t*Entering a state i, and when the stay time of the state i is st, the residual life of the state i is equal to the sum of the effective residual stay time of the operation state i and the stay time of other residual states,
wherein the effective remaining dwell time in the operating state i
Figure BDA0001933654120000048
Figure BDA0001933654120000049
For the apparatus at time t*Entering a state i, and when the stay time of the state i is st, operating the effective residual stay time of the state i; ruliFor the desired dwell time of the device in state i, ruli＝μi；
The device being at time t*Entering a state i, and when the stay time of the state i is st, the remaining life is as follows:
Figure BDA0001933654120000051
drawings
FIG. 1 is a technical roadmap for the present invention.
FIG. 2 is a comparison of multiple features of a bearing.
Fig. 3 is a graph of feature data segmentation effects.
Detailed Description
And (I) extracting a full-life operation characteristic vector, segmenting based on a multivariate fuzzy segmentation algorithm, and performing sample time normalization by adopting a DTW algorithm to finish sample data pretreatment.
The non-extensive entropy Tsallis entropy is an expression of entropy containing an index q proposed by ConstantinoTsallis in Brazil physicist in 1988, is used for describing the disorder degree of a system based on a non-extensive statistical mechanism generated by Boltzmann-Gibbs theory, has good description performance on the system with long-range interaction, long-term memory influence or system evolution in space and time of a multi-fractal sample, and has wide application in various fields of physics, chemistry, biology, medicine, economy, geophysical and the like. The Tsallis entropy has the limitation of concavity, stability and entropy generation rate in unit time, however, the entropy of shanon entropy, Renyi entropy and the like which has very rich application in the field of bearing analysis does not have the characteristics, and is reflected in the characteristic extraction of the bearing, and the signal is not well represented. Therefore, the invention extracts the Tsallis entropy of the vibration signal as a data feature, which is defined as:
Figure BDA0001933654120000052
the Tsallis entropy is to find an unambiguous q value (generally not equal to 1) to characterize as little as possible a phenomenon that is neither regular nor completely chaotic or random. Through a plurality of experiments, the variable characteristics of the bearing can be better described when q is 2.
(II) partitioning the state of the component by adopting a fuzzy security domain algorithm
The invention expands a bearing degradation model based on a safety domain estimation theory to multiple safety sub-domains, and adopts a time sequence fuzzy clustering algorithm to realize the automatic positioning of the boundary of the safety domain and the non-safety domain and the critical position of each safety sub-domain.
Let the sample sequence be
Figure BDA0001933654120000053
Where N is the data length of time series X. If xiFor multiple samples, xi＝{xij I j 1, 2.., M }, M being the dimension of the sample.
One segment of X may represent a time-continuous series of sample points za,za+1,...,zbAnd is marked as S (a, b) { a ≦ i ≦ b }. Suppose that the sample sequence X can be divided into c segments of non-overlapping portions, denoted as
Figure BDA0001933654120000061
Wherein a is1＝1,bc＝m,ak＝bk-1+1, i.e. the presence of segment boundaries s1＜s2＜…scSo that Sk(sk-1+1,sk)。
Time series segmentation can be regarded as acquiring clusters of misaligned segment constraints along a time axis, and generally, clustering processing is performed on acquired time series data points under a continuous condition according to similarity of the acquired time series data points.
Defining a sample sequence segmentation objective function as
Figure BDA0001933654120000062
It is generally defined as the distance between the true data points of the time series and the data points of a fitting function (typically a linear or polynomial function) to the sample series.
Figure BDA0001933654120000063
I.e. by minimizing the objective function, the zone boundary a is calculatedk,bkAnd k is 1, 2.., c, and an optimal segmentation position is obtained. Currently, dynamic programming and various heuristic algorithms are commonly used to minimize the objective function.
Figure BDA0001933654120000064
Wherein the content of the first and second substances,
Figure BDA0001933654120000065
is the cluster center of the kth security sub-domain;
Figure BDA0001933654120000066
is xiDistance to cluster center;
βk(ti) Is xiMembership function belonging to the k-th security sub-domain.
In practical situations, fuzzy boundaries exist between security sub-domains, which are not suitable for using fragile membership functions. Therefore, the invention adopts the multivariate Gaussian membership function to solve the problem of fuzzy boundary.
Similar to the improved Gath-Geva clustering algorithm, the method adopts a multivariate mixed Gaussian function as a sequence fitting function of a clustering prototype, and obtains the optimal region division by minimizing the sum of squares of weighted distances between data points and the center of the clustering prototype. Namely:
Figure BDA0001933654120000067
μi,krepresenting sample data points
Figure BDA0001933654120000068
A degree of membership in the kth security sub-domain, k ═ 1, 2., c; m is [1, ∞) is a weighting index that determines the fuzzy clusters that are generated, typically m is 2; d2(zi,ηk) Representing a distance between the sample data point and the cluster prototype; etakIs the clustering prototype function of the kth security sub-domain, here a multivariate mixed gaussian function.
Assuming that the sample sequence obeys the expectation of vkThe covariance matrix is FkMultivariate Gaussian distribution of (1), p (z)iEta) represents the probability density function that the sample data point belongs to the c security sub-domains.
Figure BDA0001933654120000071
Wherein the content of the first and second substances,
Figure BDA0001933654120000072
p(ηk) It is the unconditional clustering probability that,
Figure BDA0001933654120000073
ηkfor the cluster prototype function of the kth safety sub-domain, ηk＝{p(ηk),vk,Fk|k＝1,2,...,c}。
Distance function D of Gath-Geva clustering algorithm according to probability theory2(zi,ηk) And ziDegree of membership p (z) in the kth security sub-domaini|ηk) In inverse proportion, and a time variable t in the sample dataiAnd a characteristic variable xiIndependently of each other, there are:
Figure BDA0001933654120000074
wherein alpha iskIs the initial probability of clustering;
Figure BDA0001933654120000075
the distance between the time variable of the ith sample data point and the center of the clustering time variable;
Figure BDA0001933654120000076
is the distance of the feature variable in the ith sample data point from the center of the clustered feature,
Figure BDA0001933654120000077
for the cluster center of the kth safety sub-domain in the feature space, r is the feature variable distance norm AkIs determined.
Distance norm AkThere are many ways of calculating (A), the Mahalanobis distance is a good choice and can be used to adjust the correlation between variables, where Ak＝FkWherein is FkIs a fuzzy covariance matrix of a multivariate gaussian distribution,
Figure BDA0001933654120000081
to facilitate covariance matrix FkThe strong correlation between variables must be eliminated. Principal Component Analysis (PCA) can perform a series of operations and transformations on high-dimensional data, and eliminates the correlation among the high-dimensional data while retaining the original variable information as much as possible, thereby realizing the dimension reduction.
Let covariance matrix FkWith q non-zero eigenvalues (in descending order) λ k,l1,2, q and their corresponding feature vectors,
is provided with
Figure BDA0001933654120000082
For sample dataCharacteristic variable xiReducing to q dimension by PCA algorithm to obtain
Figure BDA0001933654120000083
Wherein
Figure BDA0001933654120000084
According to the ppca (probabilistic principal component analysis) algorithm,
Figure BDA0001933654120000085
wherein R iskIs any one of the q x q orthogonal transformation matrices,
Figure BDA0001933654120000086
so far, the fuzzy security domain automatic partitioning algorithm is converted into an optimization problem:
an objective function:
Figure BDA0001933654120000087
Figure BDA0001933654120000088
constraint conditions are as follows:
Figure BDA0001933654120000089
Figure BDA00019336541200000810
the problem can be solved by using an alternating optimization AO (optimization) algorithm and a Picard iteration algorithm.
The method comprises the following basic steps:
initialization:
giving the segmentation number of the sample sequence X and the space dimension of the feature vector reserved by the PCA algorithm, and selecting the proper stopping condition epsilon, epsilon>0, initialization
Figure BDA0001933654120000091
And (3) cyclic calculation: for the
Computing clustering prototype ηkThe parameters of (2):
clustering initial probability:
Figure BDA0001933654120000092
clustering center:
Figure BDA0001933654120000093
wherein
Figure BDA00019336541200000913
And (3) updating the weight:
Figure BDA0001933654120000096
and (3) updating the variance:
Figure BDA0001933654120000097
distance norm:
Figure BDA0001933654120000098
calculating the clustering prototype center and standard deviation parameter of the sample sequence time parameter:
Figure BDA00019336541200000914
clustering and merging:
for two adjacent security sub-domains So,SpAnd comparing the similarity of the two to determine whether the combination is needed. Since the PCA algorithm is used in the foregoing, the PCA similarity factor based method proposed by Krzanowski is adopted as one of the merging criteria,
Figure BDA00019336541200000911
wherein, Uo,q,Up,qAre respectively a security sub-domain So,SpThe first q principal components of the feature vector.
Another merging criterion is the Security sub-Domain So,SpDistance between feature vector cluster centers:
Figure BDA00019336541200000912
because the clustering process is carried out in the whole range of the sample sequence, the fuzzy decision algorithm is adopted to measure the clustering similarity of each safety subdomain in the whole range, and the whole similarity matrix of the decision process is H ═ Ho,p|1≤o≤c,1≤p≤c},
Figure BDA0001933654120000101
When h is generatedo,o+1Above a set threshold, the safety sub-field So,So+1Merging is carried out until the algorithm termination condition epsilon is reached.
Establishing a time-varying Markov model, calculating an initial forward probability, a backward probability and a conditional probability formula of a state sequence by using an improved forward-backward algorithm, and reestimating model parameters; obtaining the expectation of the stay time of the equipment in the state i according to the determined current running state i of the equipment, and solving the stay time of the equipment in the states i +1, i +2, …, N; calculating the dwell time of the current state i as
Figure BDA0001933654120000106
The remaining life of the battery.
And the division of the component security domains provides scientific basis for selecting the state quantity of the Markov process. We consider that the transition probability between states is not only related to the state the system is in, but also to the time the system stays in the current state. The traditional HSMM model albeitThe Markov chain of the HMM model is improved, namely the state dwell time obeys general distribution, but the influence of the state dwell time on the state transition probability is not considered, so that a time-varying transition state is introduced into the Markov chain model to obtain a state transition matrix (sojourn time based transition probabilities) related to the dwell time, and a reliability evaluation model for realizing the time-varying Markov process is recorded as M ═ M (pi, A, ST), wherein pi is the initial state probability distribution pi ═ pi { [ pi ] pi { (pi, A, ST) }i},1≤i≤N，πi＝P[s1＝i]I is more than or equal to 1 and less than or equal to N, and N is a state number; a is a time-varying state transition matrix
Figure BDA0001933654120000102
Wherein
Figure BDA0001933654120000103
Indicating that the residence time of the device in state i at time t is
Figure BDA0001933654120000104
The state transition probability of the system transitioning from state i to state j; ST (ST)iThe maximum dwell time at state i. At the same time, the state transition probability satisfies
Figure BDA0001933654120000105
ST is the state dwell time distribution, assuming a gaussian-compliant distribution.
Firstly, determining the residence time distribution of each state according to a sample set, and estimating model parameters by adopting an improved Baum-Welch algorithm.
Improved forward-backward algorithm
(1) When the model M is known as the (pi, A, ST) parameter, the stay time of the system in the state i at the time T is T*Probability of (2)
Figure BDA0001933654120000111
When t is 1, the initial forward probability is:
Figure BDA0001933654120000112
Figure BDA0001933654120000113
Figure BDA0001933654120000114
when T is 2,3, …, T, the forward probability recurrence formula is:
Figure BDA0001933654120000115
Figure BDA0001933654120000116
Figure BDA0001933654120000117
(2) the residence time of state i at time t is
Figure BDA0001933654120000118
The state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
Figure BDA0001933654120000119
When T is T, the initial backward probability formula is:
Figure BDA00019336541200001110
Figure BDA00019336541200001111
Figure BDA00019336541200001112
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
Figure BDA00019336541200001113
Figure BDA00019336541200001114
based on the forward probability formula and the backward probability formula, the state vector S ═ (S) can be derived1,s2,...,sT)，s1,s2,...,sTThe conditional probability formula of epsilon { i,1 is more than or equal to i and less than or equal to N } is
Figure BDA00019336541200001115
(3) Parameter estimation
Based on a forward-backward probability formula, a parameter re-estimation formula in a time-varying semi-Markov process model can be obtained
Figure BDA0001933654120000122
ζt(i, j, st) for a given model M and a state sequence of S ═ S1,s2,...,st)，s1,s2,...,stWhen the e is larger than or equal to { i,1 and smaller than or equal to i and smaller than or equal to N }, the probability that the system is transferred to the state j after the retention time of the state i is st.
Let gammat(i, st) given the model M and the state sequence S, the probability that the system stays in state i for a time st at time t,
Figure BDA0001933654120000123
setting etat(i, st) is the probability of the system staying at the state i for the time st at the moment t,
Figure BDA0001933654120000129
by
Figure BDA0001933654120000124
In a clear view of the above, it is known that,
Figure BDA0001933654120000125
Figure BDA0001933654120000126
represents the desired number of transitions from state i to state j while the system remains in state i for a time st
Figure BDA0001933654120000127
Represents the desired number of transitions of the system from state i when the system stays in state i for a time st, and therefore aijThe formula for reevaluation of (st) is
Figure BDA0001933654120000128
Assuming that the probability of residence time in state i in model M obeys the mean value μiVariance is
Figure BDA0001933654120000131
The distribution of the gaussian component of (a) is,
Figure BDA0001933654120000132
Figure BDA0001933654120000133
(4) reliability determination and update
If the failure rate function of the device is λ (t), the lifetime distribution function is F (t), the density function is F (t), F (t) ═ F' (t), and the reliability function is r (t), F (t) + r (t) ═ 1, λ (t) ═ F (t)/r (t) exist.
When the total number of samples is M, M (t) represents the number of samples with faults before the time t, Δ M (t) represents the number of samples with faults in the time interval (t, t + Δ t), and if,
Figure BDA0001933654120000134
to pair
Figure BDA0001933654120000135
With R (t)>0, therefore, λ (t) can be considered as an estimate of the (t, t + Δ t) interval fault condition probability.
The life cycle of the equipment is represented by L, RUL (t) represents the residual life expectancy from t to system failure when the equipment operates normally at the moment t,
Figure BDA0001933654120000136
the known apparatus is at time t*Entering a state i, and when the stay time of the state i is st, the residual life of the state i is equal to the sum of the effective residual stay time of the operation state i and the stay time of other residual states,
wherein the effective remaining dwell time in the operating state i
Figure BDA0001933654120000137
Figure BDA0001933654120000138
For the apparatus at time t*Entering a state i, and when the stay time of the state i is st, operating the effective residual stay time of the state i; ruliTo be provided withExpected residence time ready at state i, ruli＝μi。
Thus, the device is at time t*Entering a state i, and when the stay time of the state i is st, the remaining life is as follows:
Figure BDA0001933654120000141
example 1
The rolling bearing is a universal mechanical part which is most widely applied in rail transit vehicles, but the failure occurrence rate is high, only 10% -20% of the rolling bearings can reach the design life according to statistics, and the rolling bearing often causes abnormal machine performance due to various reasons such as fatigue, abrasion, strain, electric corrosion, fracture, gluing and the like in the use process and cannot work normally. Therefore, the bearing is taken as an example in this chapter to verify the scientificity of the method.
The invention adopts a French FEMTO-ST institute, and is designed and realized by AS2M department, and is specially used for testing and verifying the bearing total life data collected on a PRONOSTIA experimental platform. PRONOSTIA consists of three main components: a rotating part, a loading part (exerting radial force on the test bearing) and a measuring part.
The rotating part comprises an asynchronous motor with a gearbox and its two shafts, one of which is close to the motor and the second of which is placed on the running side of the incremental encoder. The motor power, 250W, transmits rotational motion through the gearbox, which allows the motor to reach a rated speed of 2830 revolutions, thereby providing a rated torque while maintaining the speed of the layshaft at a speed of less than 2000 revolutions per minute. The compliant and rigid couplings are used to create a shaft support bearing to which the connection for transmitting rotational motion is created by the motor. The bearing support shaft guides the bearing through its inner race. This remains fixed to the shaft, with a shoulder on the right and a threaded locking ring on the left. The shaft made of one piece is held by two pillows and their gearwheel. The two clamps allow the shaft to be longitudinally blocked between the two pillows. A human machine interface allows an operator to set speed, select directional rotation of the motor and set monitoring parameters, such as instantaneous temperature of the motor as a percentage of the maximum use temperature.
And a loading part: the assembly of the components is divided into distinct and identical aluminum plates, with partial isolation from the instrument portion by a thin layer of polymer. The aluminum plate supports a pneumatic jack, a vertical shaft and a lever arm thereof, a force sensor, a test bearing of a clamping ring, a support test bearing shaft, two bearing blocks and a large-scale oversized bearing thereof. The force from the jack is first amplified by the lever arm and then indirectly applied to the outer race of the test ball bearing via its clamping ring. This loading section forms the core of the overall system. In fact, the radial force reduction determines the service life of the bearing to be 4000N by setting its value to the maximum dynamic load of the bearing (see appendix a.1). The load is generated by a force actuator consisting of a pneumatic jack, the supply pressure of which is provided by a digital electro-pneumatic regulator.
The running condition of the bearing is measured by radial force applied to the bearing, the rotating speed of the bearing and torque on the bearing, and the sampling frequency is 100 Hz; the degradation of the bearing is measured by a vibration sensor and a temperature sensor, vibration signals are collected by two accelerometers placed at an angle of 90 degrees, the model DYTRAN 3035B is respectively placed on a horizontal axis and a vertical axis, the sampling frequency is 25.6kHz, the sampling is carried out once every 10 seconds, the sampling time is 0.1s every time, and the sampling time comprises 2560 sample points; the temperature sensor is a platinum RTD PT100 resistance temperature detector, is positioned in a hole close to an outer ring of the bearing, the sampling frequency is 10Hz, and 600 points are sampled per minute.
Samples that did not contain temperature data and were significantly misleading in temperature data were excluded.
(1) Feature extraction
Because the bearing data comprises two sets of vibration acceleration and one set of temperature data, entropy is calculated on the bearing vibration acceleration data, and the temperature data is aligned with the vibration acceleration data. In fig. 2, Tsallis entropy, kurtosis, Shannon entropy and Renyi entropy of q ═ 2 are extracted from the same life cycle signal of the bearing, and it can be known from signal representation that Shannon entropy and Renyi entropy are not monotonous in the life cycle signal and have large fluctuation; the trend of the kurtosis value is not obvious enough, the middle part is larger than the value after the fault occurs, and the performance degradation of the bearing cannot be accurately expressed by the three characteristics; the Tsallis entropy not only shows a monotonous characteristic in the whole life cycle, but also shows a descending trend when the bearing performance starts to degrade (at 1500. sample number in the figure), and is reflected in the descending trend, and the numerical stability is better, so that compared with kurtosis, Shannon entropy and Renyi entropy, the Tsallis entropy not only can reflect the degradation trend of the bearing performance, but also is very stable in performance, and a failure threshold value is determined by using the following steps.
The bearing temperature signals are resampled to the same length of the vibration characteristic signals, the maximum value is extracted, the three groups of signals are respectively normalized and converted into a dimensionless expression.
(2) Data segmentation results
The feature data in the normalized data set is automatically segmented, and the input features are [ tsalis 1, tsalis 2, Temperature ], and the obtained result is shown in fig. 3.
Therefore, the segmentation results of Bearing1_1 are [0,360], [361,1800], [1801,2740] and [2741,2803], namely the degradation process of the Bearing can be divided into four stages, namely a running-in stage, a normal operation stage, a degradation failure stage and a complete failure stage are considered to be consistent with actual use, and the effectiveness and scientific rationality of the multivariate fuzzy segmentation algorithm are fully verified. We treat the complete failure phase as a non-secure domain; the running-in stage, the normal operation stage and the degradation failure stage are regarded as security domains, and the three stages are respectively marked as a security sub domain 1, a security sub domain 2 and a security sub domain 3.
(3) Secure domain state partitioning
Training sample alignment: different bearings have different degradation tracks which are not identical, and the degradation amounts are different, so when the samples are regularly aligned, the samples are not regulated according to the degradation amounts, but according to the variation trend of the degradation amounts, namely the slopes of all points of the degradation tracks.
Based on the segmentation position of Bearing1_1, the segmentation position of Bearing1_2 is ([1,289], [290,392], [393,844], [845,871 ]).
However, in many cases, not all bearings will experience these four degradation stages, and in the samples we tested, there were Bearing2_4 that failed directly without significant degradation stages, and Bearing2_5 that was in degradation stages until failure. Taking Bearing2_4 as an example, the validity of the algorithm is verified.
Break-in phase [1:402], degenerate failure phase [403,740], and full failure phase [741,751] of Bearing2_ 4. It was therefore concluded that: the dynamic time warping algorithm based on the slope still has good matching performance on samples which do not completely accord with the degradation rule, and negative effects on degradation rule summarization due to forced alignment are avoided.
The results of the segmentation for the different bearing samples are shown in table 1.
Table 1 total sample segmentation results
Figure BDA0001933654120000161
(4) Time varying Markov process model
In the time-varying Markov process model, the state transition probability is assumed to obey a mixed Gaussian distribution probability density function, the state stay probability distribution function obeys a single Gaussian distribution probability density function, and the number of states N is 4 according to a multivariate fuzzy segmentation algorithm. The maximum iteration step number in the training process is set to 1000, and the error convergence of the algorithm is 1 × e-5。
According to the bearing security domain evaluation result, an initial matrix of mutual conversion between each security sub domain and each non-security domain and the average residence time of each security domain can be obtained, and the results are shown in tables 2 and 3.
TABLE 2 bearing security domain initial transformation matrix
Figure BDA0001933654120000171
TABLE 3 mean residence time of bearings in each safety domain
Figure BDA0001933654120000172
By utilizing the full-life data of the rolling bearing, a 4-state time-varying half Markov process fault rate prediction model is obtained, and even if mechanical parts are in the same operation state, the transition probability between the operation states and the mean value and the variance of the residence time of each operation state are different due to different residence times in the operation state. Tables 4 and 5 show the mean and variance of the operating state transition probability and state dwell time when the component is in the security sub-domain 1 with a dwell time of 10. Tables 6 and 7 show the mean and variance of the operating state transition probability and state dwell time when the component is in the security sub-domain 2 with a dwell time of 4.
TABLE 4
Figure BDA0001933654120000173
Figure BDA0001933654120000181
TABLE 5
Figure BDA0001933654120000182
TABLE 6
Figure BDA0001933654120000183
TABLE 7
Figure BDA0001933654120000184
For the bearing used in the experiment, the effective life is the life value of the safety subdomains 1,2 and 3, and the bearing completely fails when reaching the non-safety domain state.
Compared with the HMM, the predicted service life of the algorithm provided by the invention has higher accuracy, and accords with the gradual decline trend of the service life along with the time.

Claims (1)

1. A life prediction method based on a time-varying Markov process is characterized by comprising the following steps:
(1) extracting a full-life operation characteristic vector, segmenting based on a multivariate fuzzy segmentation algorithm, and performing sample time warping by adopting a dynamic time warping algorithm to finish sample data preprocessing;
(2) dividing the states of the components by adopting a fuzzy security domain algorithm;
(3) establishing a time-varying Markov model, and calculating a conditional probability formula of initial forward probability, backward probability and a state sequence by using an improved forward-backward algorithm; reestimating the time-varying Markov model parameters; obtaining the expectation of the stay time of the component in the state i according to the determined current operation state i of the component, and solving the stay time of the component in the states i +1, i +2, …, N; calculating the dwell time of the current state i as
Figure FDA00028305664700000112
The remaining life of;
Figure FDA00028305664700000110
representing the dwell time of the system in state i at time t; t is*Refers to the time that the system has been in residence,
Figure FDA00028305664700000111
refers to the dwell time at time t when the system has been in state i; t represents the maximum possible residence time;
firstly, determining the residence time distribution of each state according to a sample set, and estimating time-varying Markov model parameters by adopting an improved forward-backward algorithm:
1) when the time-varying Markov model M is known as (pi, A, ST) parameter, pi is initial state probability distribution, A is time-varying state transition matrix, and the staying time of the system in the state i at the moment T is T*Probability of (2)
Figure FDA0002830566470000011
When t is 1, the initial forward probability is:
Figure FDA0002830566470000012
Figure FDA0002830566470000013
Figure FDA0002830566470000014
when T is 2,3, …, T, the forward probability recurrence formula is:
Figure FDA0002830566470000015
Figure FDA0002830566470000016
Figure FDA0002830566470000017
2) the dwell time of state i at time t is
Figure FDA0002830566470000018
The state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
Figure FDA0002830566470000019
When T is T, the initial backward probability formula is:
Figure FDA0002830566470000021
Figure FDA0002830566470000022
Figure FDA0002830566470000023
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
Figure FDA0002830566470000024
Figure FDA0002830566470000025
based on the forward probability formula and the backward probability formula, the state vector S ═ (S) can be derived1,s2,...,sT)，s1,s2,...,sTThe conditional probability formula of epsilon { i,1 is more than or equal to i and less than or equal to N } is
Figure FDA0002830566470000026
3) Parameter estimation
Obtaining a parameter re-estimation formula in a time-varying Markov model in a time-varying half Markov process based on a forward-backward probability formula
Figure FDA0002830566470000027
ζt(i, j, st) is given a time-varying Markov model M and the state sequence is S ═ S1,s2,...,st)，s1,s2,...,stWhen the e is larger than or equal to { i,1 and is less than or equal to N }, the probability that the system is transferred to the state j after the retention time of the state i is st;
γt(i, st) given the time-varying Markov model M and the state sequence S, the probability that the system stays at state i for time st at time t,
Figure FDA0002830566470000031
setting etat(i, st) probability of the system staying in state i for time st at time t when no model and state sequence are given,
Figure FDA0002830566470000032
by
Figure FDA0002830566470000033
In a clear view of the above, it is known that,
Figure FDA0002830566470000034
Figure FDA0002830566470000035
represents the desired number of transitions from state i to state j while the system remains in state i for a time st
Figure FDA0002830566470000036
The presentation system is inThe desired number of transitions from state i, hence a, for the system at state i dwell time stijThe formula for reevaluation of (st) is
Figure FDA0002830566470000037
The probability of residence time in state i in the time-varying Markov model M obeys a mean value of μiVariance is
Figure FDA0002830566470000038
The distribution of the gaussian component of (a) is,
Figure FDA0002830566470000039
Figure FDA00028305664700000310
4) reliability determination and update
The failure rate function of the part is λ (t), the life distribution function is F (t), the density function is F (t), F (t) ═ F' (t), and the reliability function is r (t), F (t) + r (t) ═ 1, λ (t) ═ F (t)/r (t);
when the total number of samples is M, M (t) represents the number of samples with faults before the time t, Δ M (t) represents the number of samples with faults in the time interval (t, t + Δ t), and if,
Figure FDA0002830566470000041
to pair
Figure FDA0002830566470000042
R (t) > 0, therefore, lambda (t) can be regarded as the estimation of the fault condition probability in the (t, t + delta t) interval;
l denotes the life cycle of the component, RUL (t) denotes the time tThe part is operating normally, from t to the residual life expectancy of the system failure,
Figure FDA0002830566470000043
the known component enters the state i at the time t and, at a dwell time st of the state i, has a residual life equal to the sum of the effective residual dwell time in the operating state i and the dwell time in the other residual states,
wherein the effective remaining dwell time in the operating state i
Figure FDA0002830566470000044
Figure FDA0002830566470000045
The method comprises the steps that a component enters a state i at a time t, and when the stay time of the state i is st, the effective residual stay time of the component in an operation state i is obtained; ruliDesired dwell time of the component in state i, ruli＝μi；
The remaining life of the component at time t into state i and at the dwell time st of state i is:
Figure FDA0002830566470000046
CN201910001112.1A 2019-01-02 2019-01-02 Life prediction method based on time-varying Markov process Active CN109740255B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910001112.1A CN109740255B (en) 2019-01-02 2019-01-02 Life prediction method based on time-varying Markov process

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910001112.1A CN109740255B (en) 2019-01-02 2019-01-02 Life prediction method based on time-varying Markov process

Publications (2)

Publication Number Publication Date
CN109740255A CN109740255A (en) 2019-05-10
CN109740255B true CN109740255B (en) 2021-02-09

Family

ID=66363140

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910001112.1A Active CN109740255B (en) 2019-01-02 2019-01-02 Life prediction method based on time-varying Markov process

Country Status (1)

Country Link
CN (1) CN109740255B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112016866B (en) * 2019-05-31 2023-09-26 北京京东振世信息技术有限公司 Order data processing method, device, electronic equipment and readable medium
TWI706149B (en) * 2019-12-04 2020-10-01 財團法人資訊工業策進會 Apparatus and method for generating a motor diagnosis model
CN111260504B (en) * 2020-02-11 2020-11-17 南京瀚元科技有限公司 Intelligent power grid monitoring method and system and intelligent power grid controller

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011156587A2 (en) * 2010-06-09 2011-12-15 Daiichi Sankyo, Inc. Methods and systems for anticoagulation risk-benefit evaluations
CN108763644A (en) * 2018-04-24 2018-11-06 招商新智科技有限公司 A kind of method and apparatus of highway electromechanical equipment durability analysis

Also Published As

Publication number Publication date
CN109740255A (en) 2019-05-10

Similar Documents

Publication Publication Date Title
CN109740255B (en) Life prediction method based on time-varying Markov process
CN109800487B (en) Rail transit rolling bearing service life prediction method based on fuzzy security domain
CN106484949B (en) Momenttum wheel fail-safe analysis and method for predicting residual useful life based on degraded data
Shen et al. A new intermediate-domain SVM-based transfer model for rolling bearing RUL prediction
Said et al. Machine learning technique for data-driven fault detection of nonlinear processes
US20210064988A1 (en) Reliability calculation method of the thermal error model of a machine tool based on deep neural network and the monte carlo method
CN113311803B (en) On-orbit spacecraft flywheel fault detection method based on kernel principal component analysis
CN106980910B (en) Medium-and-long-term power load measuring and calculating system and method
CN113188794B (en) Gearbox fault diagnosis method and device based on improved PSO-BP neural network
CN110244224B (en) Load parameter identification method for misalignment fault of wind driven generator rotor
Han et al. Fault detection of pneumatic control valves based on canonical variate analysis
CN114356897A (en) Electromechanical actuator fault diagnosis method based on data fusion
CN115409067A (en) Method for predicting residual life of numerical control machine tool assembly
CN114897103A (en) Industrial process fault diagnosis method based on neighbor component loss optimization multi-scale convolutional neural network
Yang et al. Similarity-based information fusion grey model for remaining useful life prediction of aircraft engines
Wang et al. An improved triplet network for electromechanical actuator fault diagnosis based on similarity strategy
CN113344099A (en) Mechanical equipment degradation point identification method and system based on variational self-encoder
Meng et al. Health condition identification of rolling element bearing based on gradient of features matrix and MDDCs-MRSVD
Tobon-Mejia et al. CNC machine tool health assessment using dynamic bayesian networks
Cempel Decomposition of the symptom observation matrix and grey forecasting in vibration condition monitoring of machines
Gao et al. An optimized clustering approach using simulated annealing algorithm with HMM coordination for rolling elements bearings’ diagnosis
CN111475986A (en) L STM-AON-based gear residual life prediction method
CN112651444A (en) Self-learning-based non-stationary process anomaly detection method
Han et al. A new condition monitoring and fault diagnosis system of induction motors using artificial intelligence algorithms
Jiang et al. Condition monitoring of rolling element bearing based on Phase-PCA

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