CN109740255A - Life-span prediction method based on time-varying markoff process - Google Patents
Life-span prediction method based on time-varying markoff process Download PDFInfo
- Publication number
- CN109740255A CN109740255A CN201910001112.1A CN201910001112A CN109740255A CN 109740255 A CN109740255 A CN 109740255A CN 201910001112 A CN201910001112 A CN 201910001112A CN 109740255 A CN109740255 A CN 109740255A
- Authority
- CN
- China
- Prior art keywords
- state
- time
- probability
- life
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 33
- 230000008569 process Effects 0.000 title claims abstract description 19
- 230000007704 transition Effects 0.000 claims description 17
- 230000011218 segmentation Effects 0.000 claims description 16
- 238000009826 distribution Methods 0.000 claims description 14
- 239000013598 vector Substances 0.000 claims description 10
- 238000005315 distribution function Methods 0.000 claims description 4
- 230000014759 maintenance of location Effects 0.000 claims description 3
- BULVZWIRKLYCBC-UHFFFAOYSA-N phorate Chemical compound CCOP(=S)(OCC)SCSCC BULVZWIRKLYCBC-UHFFFAOYSA-N 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 2
- 230000006870 function Effects 0.000 description 24
- 238000006731 degradation reaction Methods 0.000 description 20
- 230000015556 catabolic process Effects 0.000 description 18
- 239000011159 matrix material Substances 0.000 description 9
- 238000000513 principal component analysis Methods 0.000 description 7
- 108020001568 subdomains Proteins 0.000 description 7
- 238000005096 rolling process Methods 0.000 description 6
- 238000005070 sampling Methods 0.000 description 6
- 238000012360 testing method Methods 0.000 description 5
- 230000001133 acceleration Effects 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000000638 solvent extraction Methods 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 description 2
- 229910052782 aluminium Inorganic materials 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- BASFCYQUMIYNBI-UHFFFAOYSA-N platinum Chemical group [Pt] BASFCYQUMIYNBI-UHFFFAOYSA-N 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000005299 abrasion Methods 0.000 description 1
- 238000004026 adhesive bonding Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000000739 chaotic effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000013210 evaluation model Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000002955 isolation Methods 0.000 description 1
- 230000007787 long-term memory Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 229910052697 platinum Inorganic materials 0.000 description 1
- 229920000642 polymer Polymers 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 229940036051 sojourn Drugs 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
Landscapes
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The present invention provides a kind of life-span prediction method based on time-varying markoff process, this method is by establishing time-varying Markov model, initial forward probability is calculated using improved Forward-backward algorithm, the condition probability formula of backward probability and status switch, and revaluation is carried out to model parameter and obtains equipment in the expectation of state residence time according to equipment current operating conditions are had determined that, and equipment is solved in the state duration of each state, calculate the remaining life of current state residence time.
Description
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 time-varying MarkovThe model calculates an initial forward probability, a backward probability and a conditional probability formula of a state sequence by utilizing 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 asThe 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 stay time of the system in the state i at the time T is T*Probability of (2)
When t is 1, the initial forward probability is:
when T is 2,3, …, T, the forward probability recurrence formula is:
2) the residence time of state i at time t isThe state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
When T is T, the initial backward probability formula is:
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
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
3) Parameter estimation
Obtaining a parameter re-estimation formula in a time-varying semi-Markov process model based on a forward-backward probability formula
ζ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.
γ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,
let ηt(i, st) is the probability of the system staying at the state i for the time st at the moment t,
byIn a clear view of the above, it is known that,
represents the desired number of transitions from state i to state j while the system remains in state i for a time stRepresents 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
The probability of residence time in state i in model M obeys the mean value μiVariance isThe distribution of the gaussian component of (a) is,
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,
to pairWith 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,
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
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:
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:
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 beWhere N is the data length of time series X. If xiFor multiple samples, xi＝{xijI 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 asWherein 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 asIt 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.
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.
Wherein,is the cluster center of the kth security sub-domain;
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:
μi,krepresenting sample data pointsA 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 the distance between the sample data point and the cluster prototype ηkIs 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)iL η) represents the probability density function that the sample data point belongs to the c security sub-domains.
Wherein,
p(ηk) It is the unconditional clustering probability that,
ηkas a function of the cluster prototype for the kth security 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:
wherein, αkIs the initial probability of clustering;
the distance between the time variable of the ith sample data point and the center of the clustering time variable;
is the distance of the feature variable in the ith sample data point from the center of the clustered feature,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,
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
For the characteristic variable x of sample dataiReducing to q dimension by PCA algorithm to obtainWherein
According to the ppca (probabilistic principal component analysis) algorithm,wherein R iskIs any one of the q x q orthogonal transformation matrices,
so far, the fuzzy security domain automatic partitioning algorithm is converted into an optimization problem:
an objective function:
constraint conditions are as follows:
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
And (3) cyclic calculation: for the
Compute cluster prototypes ηkThe parameters of (2):
clustering initial probability:
clustering center:
wherein
And (3) updating the weight:
and (3) updating the variance:
distance norm:
calculating the clustering prototype center and standard deviation parameter of the sample sequence time parameter:
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,
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:
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},
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 current running state i of the equipment according to the determined current running state i of the equipmentThe expectation of the residence time is solved, and the state residence time of the equipment in the states i +1, i +2, …, N is solved; calculating the dwell time of the current state i asThe 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. Although the conventional HSMM model improves the Markov chain of the HMM model, that is, the state dwell time obeys general distribution, but does not consider the influence of the state dwell time on the state transition probability, so that a time-varying transition state is introduced into the Markov chain model to obtain a state transition matrix (sojourn time based state transition probabilities) related to the dwell time, and a reliability evaluation model for implementing the time-varying Markov process is recorded as M ═ a, ST, where pi is an initial state probability distribution, and pi ═ a, STi},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 matrixWhereinIndicating that the residence time of the device in state i at time t isThe 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 satisfiesST 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)
When t is 1, the initial forward probability is:
when T is 2,3, …, T, the forward probability recurrence formula is:
(2) the residence time of state i at time t isThe state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
When T is T, the initial backward probability formula is:
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
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
(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
ζ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,
let ηt(i, st) is the probability of the system staying at the state i for the time st at the moment t,
byIn a clear view of the above, it is known that,
represents the desired number of transitions from state i to state j while the system remains in state i for a time stIndicating that the system is transitioning from state i when the system stays in state i for a time stThe desired number of times of emergence, therefore aijThe formula for reevaluation of (st) is
Assuming that the probability of residence time in state i in model M obeys the mean value μiVariance isThe distribution of the gaussian component of (a) is,
(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,
to pairWith 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 device is denoted by L, RUL (t) denotes that the device is operating normally at time t, fromt to the remaining life expectancy of the system failing,
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
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。
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:
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
(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
TABLE 3 mean residence time of bearings in each safety domain
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
TABLE 5
TABLE 6
TABLE 7
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 time-varying Markov model, calculating initial forward probability by improved forward-backward algorithm, and calculatingConditional probability formulas for probability and state sequences; 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 asThe 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 stay time of the system in the state i at the time T is T*Probability of (2)
When t is 1, the initial forward probability is:
when T is 2,3, …, T, the forward probability recurrence formula is:
2) the residence time of state i at time t isThe state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
When T is T, the initial backward probability formula is:
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
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
3) Parameter estimation
Obtaining a parameter re-estimation formula in a time-varying semi-Markov process model based on a forward-backward probability formula
ζ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.
γ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,
let ηt(i, st) is the probability of the system staying at the state i for the time st at the moment t,
byIn a clear view of the above, it is known that,
represents the desired number of transitions from state i to state j while the system remains in state i for a time stIndicates the desired number of transitions from state i when the system stays in state i for a time st, sinceA is aijThe formula for reevaluation of (st) is
The probability of residence time in state i in model M obeys the mean value μiVariance isThe distribution of the gaussian component of (a) is,
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,
to pairWith 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,
the known device 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
The method comprises the steps that equipment enters a state i at a time t, and when the stay time of the state i is st, the effective remaining stay time of the state i is operated; ruliFor the desired dwell time of the device in state i, ruli＝μi；
The device enters state i at time t, and the remaining lifetime when the dwell time in state i is st is:
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 true CN109740255A (en) | 2019-05-10 |
| CN109740255B 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) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111260504A (en) * | 2020-02-11 | 2020-06-09 | 吴龙圣 | Intelligent power grid monitoring method and system and intelligent power grid controller |
| TWI706149B (en) * | 2019-12-04 | 2020-10-01 | 財團法人資訊工業策進會 | Apparatus and method for generating a motor diagnosis model |
| CN112016866A (en) * | 2019-05-31 | 2020-12-01 | 北京京东尚科信息技术有限公司 | Order data processing method and device, electronic equipment and readable medium |
Citations (2)
| 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 |
-
2019
- 2019-01-02 CN CN201910001112.1A patent/CN109740255B/en active Active
Patent Citations (2)
| 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 |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112016866A (en) * | 2019-05-31 | 2020-12-01 | 北京京东尚科信息技术有限公司 | Order data processing method and device, electronic equipment and readable medium |
| 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 |
| CN111260504A (en) * | 2020-02-11 | 2020-06-09 | 吴龙圣 | Intelligent power grid monitoring method and system and intelligent power grid controller |
| CN111260504B (en) * | 2020-02-11 | 2020-11-17 | 南京瀚元科技有限公司 | Intelligent power grid monitoring method and system and intelligent power grid controller |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109740255B (en) | 2021-02-09 |
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 | |
| CN110866314B (en) | Method for predicting residual life of rotating machinery of multilayer bidirectional gate control circulation unit network | |
| Shen et al. | A new intermediate-domain SVM-based transfer model for rolling bearing RUL prediction | |
| CN109800537A (en) | A kind of machine tool thermal error model reliability degree calculation method based on deep neural network and Monte Carlo method | |
| CN111220387B (en) | Vehicle bearing residual life prediction method based on multi-feature-quantity correlation vector machine | |
| CN106980910B (en) | Medium-and-long-term power load measuring and calculating system and method | |
| Han et al. | Fault detection of pneumatic control valves based on canonical variate analysis | |
| CN110175425B (en) | Prediction method of residual life of gear based on MMALSTM | |
| US20210055706A1 (en) | Condition monitoring system | |
| CN113311803A (en) | On-orbit spacecraft flywheel fault detection method based on kernel principal component analysis | |
| CN113188794A (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 | |
| CN114356897A (en) | Electromechanical actuator fault diagnosis method based on data fusion | |
| CN115409067A (en) | Method for predicting residual life of numerical control machine tool assembly | |
| CN112149953B (en) | Electromechanical equipment operation safety assessment method based on multimode linkage and multistage cooperation | |
| CN117609679A (en) | Electric push rod fault detection method based on multi-source data | |
| CN113627090A (en) | Method and system for adjusting axis of hydraulic turbine unit and electronic equipment | |
| CN117540285A (en) | Bearing running state evaluation method based on energy entropy and regression type support vector machine | |
| Yang et al. | Fault diagnosis system of induction motors using feature extraction, feature selection and classification algorithm | |
| Liu et al. | DLVR-NWP: a novel data-driven bearing degradation model for RUL estimation | |
| Wang et al. | An improved triplet network for electromechanical actuator fault diagnosis based on similarity strategy | |
| CN110504878B (en) | Soft measurement method for rotor speed and displacement of bearingless permanent magnet synchronous motor | |
| CN108959787A (en) | Consider the thermal deformation prediction technique and system of the macro dual drive system of actual condition |
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 |