CN109800487B - Rail transit rolling bearing service life prediction method based on fuzzy security domain - Google Patents

Rail transit rolling bearing service life prediction method based on fuzzy security domain Download PDF

Info

Publication number
CN109800487B
CN109800487B CN201910001130.XA CN201910001130A CN109800487B CN 109800487 B CN109800487 B CN 109800487B CN 201910001130 A CN201910001130 A CN 201910001130A CN 109800487 B CN109800487 B CN 109800487B
Authority
CN
China
Prior art keywords
algorithm
time
domain
clustering
fuzzy
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
CN201910001130.XA
Other languages
Chinese (zh)
Other versions
CN109800487A (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 CN201910001130.XA priority Critical patent/CN109800487B/en
Publication of CN109800487A publication Critical patent/CN109800487A/en
Application granted granted Critical
Publication of CN109800487B publication Critical patent/CN109800487B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

The invention provides a life prediction method based on a fuzzy security domain, which extracts a full-life operation characteristic vector, performs segmentation based on a multivariate fuzzy segmentation algorithm, performs sample time normalization by adopting a dynamic time normalization algorithm to finish sample data preprocessing, divides the states of components by adopting the fuzzy security domain algorithm, and finally establishes the life of time-varying Markov model prediction equipment in the current state.

Description

Rail transit rolling bearing service life prediction method based on fuzzy security domain
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 above object, the reliability analysis and lifetime prediction method based on the fuzzy security domain of the present invention 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;
let the sample sequence be
Figure GDA0002779157230000011
Where N is the data length of time series X, XiFor multiple samples, xi＝{xij1, 2., M }, where M is the dimension of the sample; one segment of X may represent a time-continuous series of sample points za,za+1,...,zbAnd is denoted 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 GDA0002779157230000021
Wherein a is1＝1,bc＝m,ak＝bk-1+1, presence of segment boundary s1＜s2＜…scSo that Sk(sk-1+1,sk)；tiIs a time sequence;
defining a sample sequence segmentation objective function as
Figure GDA0002779157230000022
Figure GDA0002779157230000023
By minimizing the objective function, the region boundary a is calculatedk,bk K 1, 2.., c, obtaining an optimal segmentation position;
Figure GDA0002779157230000024
wherein the content of the first and second substances,
Figure GDA0002779157230000025
is the cluster center of the kth security sub-domain,
Figure GDA0002779157230000026
is xiDistance to the center of the cluster, βk(ti) Is xiA membership function pertaining to the kth security sub-domain;
and obtaining the optimal region division by minimizing the sum of squares of weighted distances between the data points and the center of the clustering prototype by using a multivariate mixed Gaussian function as a sequence fitting function of the clustering prototype:
Figure GDA0002779157230000027
wherein, muk,iRepresenting sample data points
Figure GDA0002779157230000028
A degree of membership in the kth security sub-domain, k ═ 1, 2., c; m ═ 1, infinity) is the weighting index that determines the fuzzy clusters that are generated; d2(zi,ηk) Representing a distance between the sample data point and the cluster prototype; etakA clustering prototype function for the kth security sub-domain;
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 of the sample data point belonging to c safety sub-domains
Figure GDA0002779157230000029
Wherein the content of the first and second substances,
Figure GDA00027791572300000210
p(ηk) It is the unconditional clustering probability that,
Figure GDA00027791572300000211
ηkfor the cluster prototype function of the kth safety sub-domain, ηk＝{p(ηk),vk,Fk|k＝1,2,...,c}；
And is
Figure GDA0002779157230000031
Wherein alpha iskIs the initial probability of clustering;
Figure GDA0002779157230000032
the distance between the time variable of the ith sample data point and the center of the clustering time variable;
Figure GDA0002779157230000033
is the distance of the feature variable in the ith sample data point from the center of the clustered feature,
Figure GDA0002779157230000034
for the cluster center of the kth safety sub-domain in the feature space, r is the feature variable distance norm AkThe rank of (d); a. thek＝FkIn which F iskIs a fuzzy covariance matrix of a multivariate gaussian distribution,
Figure GDA0002779157230000035
let covariance matrix FkWith q non-zero eigenvalues (in descending order) λ k,l1,2, q and their corresponding feature vectors,
is provided with
Figure GDA0002779157230000036
For the characteristic variable x of sample dataiReducing to q dimension by PCA algorithm to obtain
Figure GDA0002779157230000037
Wherein
Figure GDA0002779157230000038
According to the PPCA algorithm,
Figure GDA0002779157230000039
wherein R iskIs any one of the q x q orthogonal transformation matrices,
Figure GDA00027791572300000310
(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 GDA00027791572300000311
The remaining life of the battery.
Preferably, the step (2) further comprises the following steps:
the fuzzy security domain automatic partitioning algorithm is converted into an optimization problem:
an objective function:
Figure GDA0002779157230000041
constraint conditions are as follows:
Figure GDA0002779157230000042
solving by adopting an alternative optimization AO algorithm and a Picard iterative algorithm:
the method comprises the following basic steps:
initialization: giving the segment number of the sample sequence X and the feature vector space dimension reserved by the PCA algorithm, selecting a proper stop condition of more than 0, and initializing
Figure GDA0002779157230000043
And (3) cyclic calculation:
computing clustering prototype ηkThe parameters of (2):
clustering initial probability:
Figure GDA0002779157230000044
clustering center:
Figure GDA0002779157230000045
wherein
Figure GDA0002779157230000046
And (3) updating the weight:
Figure GDA0002779157230000047
and (3) updating the variance:
Figure GDA0002779157230000048
distance norm:
Figure GDA0002779157230000049
calculating the clustering prototype center and standard deviation parameter of the sample sequence time parameter:
Figure GDA00027791572300000410
clustering and merging:
Figure GDA0002779157230000051
wherein S iso,SpFor two adjacent security sub-domains, 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 GDA0002779157230000052
wherein the content of the first and second substances,
Figure GDA0002779157230000053
is So,SpA feature vector clustering center;
the fuzzy decision algorithm is adopted to measure the clustering similarity of each safety subdomain in the whole, and the whole similarity matrix of the decision process is H ═ Ho,p|1≤o≤c,1≤p≤c}，
Figure GDA0002779157230000054
When h is generatedo,o+1Above a set threshold, the safety sub-field So,So+1And merging until the end of the algorithm is reached.
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 warping by adopting a dynamic time warping algorithm to finish sample data preprocessing.
The non-extensive entropy Tsallis entropy is an expression of entropy containing an index q proposed by Constanino Tsallis in Brazil physicist in 1988, is used for describing the chaos 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 GDA0002779157230000061
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 GDA0002779157230000062
Where 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 parts, denoted as ScX＝{Sk(ak,bk) L 1 is less than or equal to k is less than or equal to c, wherein a1＝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 GDA0002779157230000063
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 GDA0002779157230000064
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 GDA0002779157230000065
Wherein the content of the first and second substances,
Figure GDA0002779157230000066
is the cluster center of the kth security sub-domain;
Figure GDA0002779157230000067
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 GDA0002779157230000071
μi,krepresenting sample data points
Figure GDA0002779157230000072
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)i|η) Representing the probability density function that the sample data point belongs to the c security sub-domains.
Figure GDA0002779157230000073
Wherein the content of the first and second substances,
Figure GDA0002779157230000074
p(ηk) It is the unconditional clustering probability that,
Figure GDA0002779157230000075
η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 GDA0002779157230000076
wherein alpha iskIs the initial probability of clustering;
Figure GDA0002779157230000081
the distance between the time variable of the ith sample data point and the center of the clustering time variable;
Figure GDA0002779157230000082
is the distance of the feature variable in the ith sample data point from the center of the clustered feature,
Figure GDA0002779157230000083
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 GDA0002779157230000084
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, there are
Figure GDA0002779157230000085
For the characteristic variable x of sample dataiReducing to q dimension by PCA algorithm to obtain
Figure GDA0002779157230000086
Wherein
Figure GDA0002779157230000087
According to the PPCA (Probabilistic principal component analysis) algorithm,
Figure GDA0002779157230000088
wherein R iskIs any one of the q x q orthogonal transformation matrices,
Figure GDA0002779157230000089
so far, the fuzzy security domain automatic partitioning algorithm is converted into an optimization problem:
an objective function:
Figure GDA0002779157230000091
constraint conditions are as follows:
Figure GDA0002779157230000092
the problem can be solved by adopting an alternating optimization AO (optimization) algorithm and a Picard iteration algorithm.
The method comprises the following basic steps:
initialization:
giving the segment number of the sample sequence X and the feature vector space dimension reserved by the PCA algorithm, selecting a proper stop condition of more than 0, and initializing
Figure GDA0002779157230000093
And (3) cyclic calculation: for the
Computing clustering prototype ηkThe parameters of (2):
clustering initial probability:
Figure GDA0002779157230000094
clustering center:
Figure GDA0002779157230000095
wherein
Figure GDA0002779157230000096
Figure GDA0002779157230000097
And (3) updating the weight:
Figure GDA0002779157230000098
and (3) updating the variance:
Figure GDA0002779157230000099
distance norm:
Figure GDA00027791572300000910
calculating the clustering prototype center and standard deviation parameter of the sample sequence time parameter:
Figure GDA00027791572300000911
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 GDA0002779157230000101
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 GDA0002779157230000102
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 GDA0002779157230000103
When h is generatedo,o+1Above a set threshold, the safety sub-field So,So+1And merging until the end of the algorithm 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 stay time of the current state i as stti remaining life.
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 the initial state probability distribution 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 matrix
Figure GDA0002779157230000111
Wherein
Figure GDA0002779157230000112
Indicating that the residence time of the device in state i at time t is
Figure GDA0002779157230000113
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 GDA0002779157230000114
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 GDA0002779157230000115
When t is 1, the initial forward probability is:
Figure GDA0002779157230000116
Figure GDA0002779157230000117
Figure GDA0002779157230000118
when T is 2,3, …, T, the forward probability recurrence formula is:
Figure GDA0002779157230000119
Figure GDA00027791572300001110
Figure GDA00027791572300001111
(2) state at time ti has a residence time of
Figure GDA00027791572300001112
The state sequence(s) generated in the following T-T timet+1,st+2,...,sT) Has a backward probability of
Figure GDA00027791572300001113
When T is T, the initial backward probability formula is:
Figure GDA00027791572300001114
Figure GDA00027791572300001115
Figure GDA00027791572300001116
when T is more than 0 and less than T, the backward probability recurrence formula is as follows:
Figure GDA0002779157230000121
Figure GDA0002779157230000122
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 GDA0002779157230000123
(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 GDA0002779157230000124
Figure GDA0002779157230000125
ζ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 GDA0002779157230000126
setting etat(i, st) is the probability of the system staying at the state i for the time st at the moment t,
Figure GDA0002779157230000127
by
Figure GDA0002779157230000128
In a clear view of the above, it is known that,
Figure GDA0002779157230000131
Figure GDA0002779157230000132
represents the desired number of transitions from state i to state j while the system remains in state i for a time st
Figure GDA0002779157230000133
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 GDA0002779157230000134
Assuming that the probability of residence time in state i in model M obeys the mean value μiVariance is
Figure GDA0002779157230000135
The distribution of the gaussian component of (a) is,
Figure GDA0002779157230000136
Figure GDA0002779157230000137
(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) is F' (t), and the reliability function is R (t), F (t) + R (t) is 1, and λ (t) is F (t)/R (t).
When the total number of samples is M, M (t) represents the number of samples with faults occurring before time t, Δ M (t) represents the number of samples with faults occurring between time intervals (t, t + Δ t), then,
Figure GDA0002779157230000138
to pair
Figure GDA0002779157230000139
There is R (t) > 0, so λ (t) can be considered as an estimate of the probability of a fault condition in the (t, t + Δ t) interval.
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 GDA00027791572300001310
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
Figure GDA0002779157230000141
Figure GDA0002779157230000142
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 remaining life of the device at time t, entering state i and at the time of stay st in state i, is therefore:
Figure GDA0002779157230000143
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 result for Bearing1_1 is [0,360], [361,1800], [1801,2740],
[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, which are consistent with actual use, and the validity 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 GDA0002779157230000161
Figure GDA0002779157230000171
(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 stay 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 GDA0002779157230000172
TABLE 3 mean residence time of bearings in each safety domain
Figure GDA0002779157230000173
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 GDA0002779157230000174
Figure GDA0002779157230000181
TABLE 5
Figure GDA0002779157230000182
TABLE 6
Figure GDA0002779157230000183
TABLE 7
Figure GDA0002779157230000184
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 (2)

1. A rail transit rolling bearing service life prediction method based on a fuzzy security domain is characterized by comprising the following steps:
(1) extracting a whole-life running vibration characteristic vector and a temperature characteristic vector of a rolling bearing of the rail transit, 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;
let the sample sequence be
Figure FDA0002746389720000011
Where N is the data length of time series X, XiFor multiple samples, xi＝{xij1, 2., M }, where M is the dimension of the sample; one segment of X may represent a time-continuous series of sample points za,za+1,...,zbAnd is denoted 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 FDA0002746389720000012
Wherein a is1＝1,bc＝m,ak＝bk-1+1, presence of segment boundary s1＜s2＜…scSo that Sk(sk-1+1,sk)；tiIs a time sequence;
defining a sample sequence segmentation objective function as
Figure FDA0002746389720000013
Figure FDA0002746389720000014
By minimizing the objective function, the region boundary a is calculatedk,bkK 1, 2.., c, obtaining an optimal segmentation position;
Figure FDA0002746389720000015
wherein the content of the first and second substances,
Figure FDA0002746389720000016
is the cluster center of the kth security sub-domain,
Figure FDA0002746389720000017
is xiDistance to the center of the cluster, βk(ti) Is xiA membership function pertaining to the kth security sub-domain;
and obtaining the optimal region division by minimizing the sum of squares of weighted distances between the data points and the center of the clustering prototype by using a multivariate mixed Gaussian function as a sequence fitting function of the clustering prototype:
Figure FDA0002746389720000018
wherein, muk,iRepresenting sample data points
Figure FDA0002746389720000019
A degree of membership in the kth security sub-domain, k ═ 1, 2., c; m ═ 1, infinity) is the weighting index that determines the fuzzy clusters that are generated; d2(zi,ηk) Representing a distance between the sample data point and the cluster prototype; etakA clustering prototype function for the kth security sub-domain;
assuming that the sample sequence obeys the expectation of vkThe covariance matrix is FkMultivariate Gaussian distribution of (1), p (z)iη) represents the probability density function of the sample data point under the c safety sub-domains
Figure FDA0002746389720000021
Wherein the content of the first and second substances,
Figure FDA0002746389720000022
p(ηk) It is the unconditional clustering probability that,
Figure FDA0002746389720000023
ηkfor the cluster prototype function of the kth safety sub-domain, ηk＝{p(ηk),vk,Fk|k＝1,2,...,c}；
And is
Figure FDA0002746389720000024
Wherein alpha iskIs the initial probability of clustering;
Figure FDA0002746389720000025
the distance between the time variable of the ith sample data point and the center of the clustering time variable;
Figure FDA0002746389720000026
is the distance of the feature variable in the ith sample data point from the center of the clustered feature,
Figure FDA0002746389720000027
for the cluster center of the kth safety sub-domain in the feature space, r is the feature variable distance norm AkThe rank of (d); a. thek＝FkIn which F iskIs a fuzzy covariance matrix of a multivariate gaussian distribution,
Figure FDA0002746389720000028
let covariance matrix FkWith q non-zero eigenvalues lambdak,l1,2, q and their corresponding feature vectors,
is provided with
Figure FDA0002746389720000029
Uk＝[uk,1,uk,2,...,uk,q]；
For the characteristic variable x of sample dataiReducing to q dimension by PCA algorithm to obtain
Figure FDA0002746389720000031
Wherein
Figure FDA0002746389720000032
According to the PPCA algorithm,
Figure FDA0002746389720000033
wherein R iskIs any one of the q x q orthogonal transformation matrices,
Figure FDA0002746389720000034
(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 FDA0002746389720000035
The remaining life of the battery.
2. The method of claim 1, wherein step (2) further comprises the following:
the fuzzy security domain automatic partitioning algorithm is converted into an optimization problem:
an objective function:
Figure FDA0002746389720000036
constraint conditions are as follows:
Figure FDA0002746389720000037
solving by adopting an alternative optimization AO algorithm and a Picard iterative algorithm:
the method comprises the following basic steps:
initialization: giving the segment number of the sample sequence X and the feature vector space dimension reserved by the PCA algorithm, selecting a proper stop condition of more than 0, and initializing Wk,
Figure FDA0002746389720000038
μk,i；
And (3) cyclic calculation:
computing clustering prototype ηkThe parameters of (2):
clustering initial probability:
Figure FDA0002746389720000039
clustering center:
Figure FDA00027463897200000310
wherein
Figure FDA00027463897200000311
And (3) updating the weight:
Figure FDA0002746389720000041
and (3) updating the variance:
Figure FDA0002746389720000042
distance norm:
Figure FDA0002746389720000043
calculating the clustering prototype center and standard deviation parameter of the sample sequence time parameter:
Figure FDA0002746389720000044
clustering and merging:
Figure FDA0002746389720000045
wherein S iso,SpFor two adjacent security sub-domains, 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 FDA0002746389720000046
wherein v iso x,vp xIs So,SpA feature vector clustering center;
the fuzzy decision algorithm is adopted to measure the clustering similarity of each safety subdomain in the whole, and the whole similarity matrix of the decision process is H ═ Ho,p|1≤o≤c,1≤p≤c}，
Figure FDA0002746389720000047
When h is generatedo,o+1Above a set threshold, the safety sub-field So,So+1And merging until the end of the algorithm is reached.
CN201910001130.XA 2019-01-02 2019-01-02 Rail transit rolling bearing service life prediction method based on fuzzy security domain Active CN109800487B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910001130.XA CN109800487B (en) 2019-01-02 2019-01-02 Rail transit rolling bearing service life prediction method based on fuzzy security domain

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910001130.XA CN109800487B (en) 2019-01-02 2019-01-02 Rail transit rolling bearing service life prediction method based on fuzzy security domain

Publications (2)

Publication Number Publication Date
CN109800487A CN109800487A (en) 2019-05-24
CN109800487B true CN109800487B (en) 2020-12-29

Family

ID=66558315

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910001130.XA Active CN109800487B (en) 2019-01-02 2019-01-02 Rail transit rolling bearing service life prediction method based on fuzzy security domain

Country Status (1)

Country Link
CN (1) CN109800487B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112165161A (en) * 2020-02-11 2021-01-01 吴龙圣 Intelligent power grid monitoring method and system based on Internet of things
CN112504676B (en) * 2020-12-24 2022-04-01 温州大学 Rolling bearing performance degradation analysis method and device
CN112947290B (en) * 2021-05-16 2021-08-20 北京赛博联物科技有限公司 Edge cloud cooperation-based equipment state monitoring method and system and storage medium

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010011918A2 (en) * 2008-07-24 2010-01-28 University Of Cincinnati Methods for prognosing mechanical systems
CN102196420A (en) * 2011-06-02 2011-09-21 河海大学常州校区 Secure clustering routing management method for wireless sensor network
CN104502103A (en) * 2014-12-07 2015-04-08 北京工业大学 Bearing fault diagnosis method based on fuzzy support vector machine
CN104820750A (en) * 2015-05-11 2015-08-05 广西大学 Structure reliability dynamic response surface method based on discriminant analysis
CN105022912A (en) * 2015-05-28 2015-11-04 北京交通大学 Rolling bearing fault prediction method based on wavelet principal component analysis

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010011918A2 (en) * 2008-07-24 2010-01-28 University Of Cincinnati Methods for prognosing mechanical systems
CN102196420A (en) * 2011-06-02 2011-09-21 河海大学常州校区 Secure clustering routing management method for wireless sensor network
CN104502103A (en) * 2014-12-07 2015-04-08 北京工业大学 Bearing fault diagnosis method based on fuzzy support vector machine
CN104820750A (en) * 2015-05-11 2015-08-05 广西大学 Structure reliability dynamic response surface method based on discriminant analysis
CN105022912A (en) * 2015-05-28 2015-11-04 北京交通大学 Rolling bearing fault prediction method based on wavelet principal component analysis

Also Published As

Publication number Publication date
CN109800487A (en) 2019-05-24

Similar Documents

Publication Publication Date Title
CN109800487B (en) Rail transit rolling bearing service life prediction method based on fuzzy security domain
CN109740255B (en) Life prediction method based on time-varying Markov process
CN112699913B (en) Method and device for diagnosing abnormal relationship of household transformer in transformer area
CN106484949B (en) Momenttum wheel fail-safe analysis and method for predicting residual useful life based on degraded data
CN112629863B (en) Bearing fault diagnosis method for dynamic joint distribution alignment network under variable working conditions
CN111220387B (en) Vehicle bearing residual life prediction method based on multi-feature-quantity correlation vector machine
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
CN109800537A (en) A kind of machine tool thermal error model reliability degree calculation method based on deep neural network and Monte Carlo method
CN106980910B (en) Medium-and-long-term power load measuring and calculating system and method
Benmoussa et al. Remaining useful life estimation without needing for prior knowledge of the degradation features
CN110210126B (en) LSTMPP-based gear residual life prediction 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
CN113554148A (en) BiLSTM voltage deviation prediction method based on Bayesian optimization
Sun et al. Composite-graph-based sparse subspace clustering for machine fault diagnosis
CN113311803A (en) On-orbit spacecraft flywheel fault detection method based on kernel principal component analysis
CN112650146A (en) Fault diagnosis optimization method, system and equipment of numerical control machine tool under multiple working conditions
CN115409067A (en) Method for predicting residual life of numerical control machine tool assembly
CN110504878B (en) Soft measurement method for rotor speed and displacement of bearingless permanent magnet synchronous motor
CN114356897A (en) Electromechanical actuator fault diagnosis method based on data fusion
Yang et al. Similarity-based information fusion grey model for remaining useful life prediction of aircraft engines
CN117540285A (en) Bearing running state evaluation method based on energy entropy and regression type support vector machine
CN111475986A (en) L STM-AON-based gear residual life prediction method
CN111475987A (en) Gear residual life prediction method based ON SAE and ON-L STM

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