CN117367800A - Multi-feature index bearing residual life prediction method based on random correlation - Google Patents

Abstract

The invention belongs to the technical field of bearing life prediction, and particularly discloses a multi-feature index bearing residual life prediction method based on random correlation, which comprises the following steps of: extracting characteristic signals in the multi-sensor by using methods such as time domain, time domain and the like, selecting characteristics according to evaluation indexes such as monotonicity, trend, robustness and the like, and analyzing various relations such as correlation among the multi-characteristics, correlation among the multi-characteristics and the multi-indexes, correlation among the multi-indexes and the like; and then constructing a plurality of performance indexes of the characterization equipment, taking correlation among the indexes and different degradation modes into consideration, constructing a multi-index degradation state space model, defining three failure modes according to first time definition, deducing a corresponding residual life probability density function, and providing important basis for predicting system faults.

Description

Multi-feature index bearing residual life prediction method based on random correlation

Technical Field

The invention belongs to the technical field of bearing life prediction, and particularly relates to a multi-feature index bearing residual life prediction method considering random correlation.

Background

Degradation such as abrasion, fatigue, cracks and the like can occur in long-term operation of industrial equipment, so that the operation performance of the equipment is reduced and even fails. For the state of degradation of the device, it is mostly not directly detectable, but indirectly monitored by various sensing devices such as vibration, temperature and pressure sensors. Meanwhile, some sensor signals cannot directly reflect the degradation state of the device, and only the degradation state of the device can be reflected by extracting the characteristics of the sensing information.

At present, the multi-index constructed based on multi-characteristics exists in the following ways: the first is that the relationship between multiple features and multiple indexes is one-to-one, which corresponds to the construction of a single index by a single feature. For a data-driven multi-index model, a machine learning algorithm and other intelligent technologies are utilized to represent the degradation behavior of the system, and the system has good nonlinear mapping capability, and as the structure of the methods is a black box, the correlation among multiple indexes cannot be intuitively embodied. And the analytical expression of the multivariate degradation level joint distribution and the system reliability can not be obtained when the degradation rate interaction relation is complex based on the degradation rate correlation modeling type. The method has the advantages that the relation between characteristic signals and indexes is one-to-one correspondence, only a single index constructed by a single characteristic is considered, and the reflected information is possibly incomplete.

The second common method is to construct an HI index by fusing a plurality of features, wherein the method only fuses the plurality of features, and the fused features have redundant features, so that the physical meaning is difficult to be clarified, and the correlation among the plurality of features cannot be embodied. The two methods cannot reflect the complex relationship between the multi-feature data and the multi-performance index. The current multi-index modeling method constructed by using multiple features needs to consider multiple correlation relationships between the multiple features and multiple indexes, which may improve the RUL (remaining service life) prediction performance of the system.

Furthermore, current multi-index degradation models are often assumed to be single degradation models, which limits the versatility of the method. In a practical system, different multi-index systems may have different degradation modes, and thus modeling of multiple different degradation modes of the multi-index system needs to be considered. The distribution of RUL predictions for the system is different in different failure modes. Thus, different failure modes have a significant impact on the RUL prediction of the system.

Disclosure of Invention

In order to solve the technical problems in the prior art, the invention provides a multi-feature index bearing residual life prediction method based on random correlation.

In order to achieve the above purpose, the technical scheme adopted by the invention is as follows: the multi-feature index bearing residual life prediction method based on random correlation comprises the following specific steps:

step one, constructing a state space model of multi-feature degradation, selecting monotonicity, robustness and trend of multi-feature data through the state space model, and analyzing the relation between the multi-feature data and the multi-index data;

step two, considering the correlation among multiple features, the correlation among multiple features and multiple indexes, the correlation among multiple indexes and the degradation characteristic of a reflection system in a state space model, considering multiple degradation modes in the state space model, and carrying out degradation state and parameter joint estimation by adopting a maximum expectation algorithm and a Kalman filtering method;

and step three, defining a redundancy failure mode, a fusion failure mode and a competition failure mode among multiple indexes, and establishing a residual life prediction model by using a state space model.

In step one, feature selection is performed using monotonicity, robustness, and trending:

1) Monotonicity characterizes the consistency of system performance degradation, and the monotonicity Mon (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For characteristic sequences, ε (Y _i -Y _i-1 ) Is Y _i -Y _i-1 Is a unit step function of [0,1 ]]Within the range;is epsilon (Y) _i -Y _i-1 )-ε(Y _i-1 -Y _i ) From the summation of i=2 to n,is->Mon (Y) is the monotonicity of the feature sequence Y;

2) The trend index reflects the linear correlation degree between the characteristic sequence of performance degradation and the service life, and the trend Corr (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For the feature sequence, n is the total number of monitoring in the whole performance degradation process, t= (T ₁ ,t ₂ ...t _n ) For the corresponding sequence of monitoring instants, the value is 0,1]Within the range;is Y _i t _i From the summation of i=1 to n,is->Corr (Y) is the trend of the characteristic sequence Y,for->Performing open square calculation;

3) The robustness index reflects the robustness of the system performance degradation characteristic to external interference, and the robustness Rob (Y) of the characteristic sequence Y is defined as:

wherein,the value range of the trend sequence of the performance degradation characteristic is [0,1 ]]，For->Obtaining an exponential function value->For->The function divides the summation of i=1 to n by n, rob (Y) is the robustness of the feature sequence Y, and based on characteristics including monotonicity, trending and robustness, the choice of multiple features is converted into a weighted summation, the calculation formula λ is:

λ＝[φ ₁ Mon(Y)+φ ₂ Tre(Y)+φ ₃ Rob(Y)] (4)

wherein lambda is E [0,1 ]]，φ _i (i=1, 2, 3) is the relative importance of the eigenvalue Y on three selection criteria, λ is a positive correlation relationship with monotonicity, correlation and robustness.

In the second step, a state space model is adopted to carry out degradation modeling on a multi-index system constructed by multiple features, extracted features are used for directly expressing degradation indexes, and the features of sensing signals reflect the degradation state of equipment;

the correlation between the multiple indexes is reflected by a correlation matrix of the state space model, and when the system is observable, the system adopts a KF algorithm (Kalman filtering) to perform state estimation.

In step three, the degradation state distribution X of the multi-index system is obtained _k+l The method comprises the following steps:

wherein,for matrices φ (k+l, j+1) and ζ _j+1 The sum of the products from j=k to k+l-1 gives further X _k+l Is->Sum of variances P _k+l ：

Let, S (l) _k )＝X(l _k +t _k )-X(t _k ) S (0) =0, i.eWherein S (l) _k ) From t for the degenerate state X _k To t _k+l Degradation delta over a period of time;

1) Failure of competition

Judging that the system fails and the corresponding residual life L exists in a plurality of degradation indexes as long as one index exceeds the corresponding failure threshold ₁ The method comprises the following steps:

L ₁ ＝inf{l＞0:X ₁ (k+l)≥w ₁ ,or...,orX _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (7)

wherein X is _i In (k+l), k+l is time t _k+l I is the number of indexes, w _i A fault threshold value preset for the index i, wherein inf { · } is the lower bound of the function;

residual life distribution function of system in contention failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )＜w ₁ ,orX ₂ (l _k +t _k )＜w ₂ ,...,orX _n (l _k +t _k )＜w _n |X _k ,θ _k Is at X } _k And theta _k In the degradation state of n indexes under the condition, one or more indexes exceed the corresponding failure threshold value X ₁ (l _k +t _k )＜w ₁ ,orX ₂ (l _k +t _k )＜w ₂ ,...,orX _n (l _k +t _k )＜w _n Is a function of the probability of (1),at the current time t as index i _k Is used to determine the amount of degradation of (1),for n index degradation states at t _k A joint probability density function of the time increment S,is->N-fold integration of the function;

2) Failure of redundancy

When the multiple degradation indexes exceed the corresponding failure thresholds, judging that the system fails and the corresponding residual life L ₂ The method comprises the following steps:

L ₂ ＝inf{l＞0:X ₁ (k+l)≥w ₁ ,...,and X _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (9)

residual life distribution function of system in redundancy failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n |X _k ,θ _k Is at X } _k And theta _k The degradation states of n indexes all exceed the corresponding failure threshold value X under the condition ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n Probability of (2);

3) Failure of fusion

When all degradation indexes do not exceed the corresponding threshold value, but the weighted summation of a plurality of degradation indexes exceeds the total failure threshold value, judging that the system fails and the corresponding residual life L ₃ The method comprises the following steps:

L ₃ ＝inf{l＞0:m ₁ X ₁ (k+l)+...+m _n X _n (k+l)≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n |X _k ,θ _k } (11)

wherein W is a fusion weight, w=λ ₁ w ₁ +..λ _i w _i +...λ _n w _n ，λ _i As the weight coefficient corresponding to index i, the residual life distribution function of the system in the fusion failure modeThe method comprises the following steps:

P{m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n |X _k ,θ _k is at X } _k And theta _k Under the condition that all indexes do not exceed the corresponding threshold value, but the degradation state weighted summation of n indexes exceeds the total failure threshold value m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n The corresponding probabilities;

then, the distribution function obtained by the formulas (8), (10) and (12) is usedObtaining RUL probability density function of the system through a finite difference formula:

wherein,for the purpose of distributing the function->Seeking a derivative;

obtaining the current time t according to the full probability formula _k Considering residual life probability density functions of implicit state systems in different failure modesThe method comprises the following steps:

wherein, psi (X) _k |θ _k ,Y _1:k ) Probability density function for implicit state of system, i.e Is state X _k Mean and variance of (a) and obtained by a filtering algorithm, < + >>For->And psi (X) _k |θ _k ,Y _1:k ) Integration at- ≡to factafter the product,to be at theta _k ,X _k ,Y _1:k Residual life under conditions l _k Probability density functionIs a mathematical expectation of (a).

Compared with the prior art, the invention has the following specific beneficial effects: in the method, various relations among a plurality of indexes and the characteristics are considered, an RUL prediction framework from the plurality of characteristics to the plurality of indexes is constructed, in the framework, the mapping relation between the characteristics and the indexes is analyzed, the plurality of indexes are adopted to reflect the degradation state of a system, a plurality of performance indexes are constructed according to the corresponding mapping relation to reflect the degradation condition of a component by utilizing a plurality of characteristic signals, the RUL prediction and degradation modeling are carried out on the component by considering different failure modes and degradation modes of the plurality of indexes, and unknown parameters in the model are calculated and updated by combining an EM, RTS and KF method. The model provided by the invention can intuitively reflect the correlation among indexes, provide RUL probability distribution, well predict the degradation trend of industrial equipment, improve the RUL prediction performance, and provide important basis for predicting the failure of the system by reflecting different degradation degrees of the system by different indexes.

Drawings

Fig. 1 is an overall flow chart of the present invention.

FIG. 2 is a diagram illustrating three correlation modeling schemes of multiple features and multiple indicators according to the present invention.

Fig. 3 is a structural diagram of the experimental platform.

Fig. 4 is a vibration signal monitoring diagram of the bearings 1 to 4, fig. 4 (a) is a monitoring diagram of a horizontal vibration signal, and fig. 4 (b) is a monitoring diagram of a vertical vibration signal.

Fig. 5 is a three-layer wavelet packet decomposed sub-band energy-to-bar graph of a vibration signal of the bearing1-4, fig. 5 (a) is a three-layer wavelet packet decomposed sub-band energy-to-bar graph of a horizontal vibration signal, and fig. 5 (b) is a three-layer wavelet packet decomposed sub-band energy-to-bar graph of a vertical vibration signal.

Fig. 6 is a schematic diagram of time domain signal feature extraction, fig. 6 (a) is a schematic diagram of horizontal vibration signal feature extraction, and fig. 6 (b) is a schematic diagram of vertical vibration signal feature extraction.

Fig. 7 is a schematic diagram of time-frequency domain signal feature extraction.

Fig. 8 is a performance index degradation state estimation diagram of the bearings 1 to 4, fig. 8 (a) is a horizontal performance index degradation state estimation diagram, and fig. 8 (b) is a vertical performance index degradation state estimation diagram.

Fig. 9 is a multi-index remaining life prediction method constructed by different methods, fig. 9 (a) is a multi-index remaining life prediction method constructed by an M1 method, fig. 9 (b) is a multi-index remaining life prediction method constructed by an M2 method, fig. 9 (c) is a multi-index remaining life prediction method constructed by an M3 method, and fig. 9 (d) is a multi-index remaining life prediction method constructed by an M4 method.

Fig. 10 is a schematic diagram showing a comparison of residual life of different prediction methods, fig. 10 (a) is a schematic diagram showing a comparison of true values and average values, and fig. 10 (b) is a schematic diagram showing a comparison of MSE values.

FIG. 11 is a graph showing the comparison of the relative errors between the predicted and actual RULs according to the present invention.

Fig. 12 is a schematic view of time domain signal feature extraction, fig. 12 (a) is a schematic view of horizontal vibration signal feature extraction, and fig. 12 (b) is a schematic view of vertical vibration signal feature extraction.

Fig. 13 is a schematic diagram of time-frequency domain signal feature extraction.

Fig. 14 is a schematic view of two performance index degradation state estimation of the bearing2-2, fig. 14 (a) is a schematic view of performance index degradation state estimation in the horizontal direction, and fig. 14 (b) is a schematic view of performance index degradation state estimation in the vertical direction.

Fig. 15 is a residual life prediction result map of the bearing2-2, fig. 15 (a) is a residual life prediction result map of the method M0, fig. 15 (b) is a residual life prediction result map of the method M1, fig. 15 (c) is a residual life prediction result map of the method M2, and fig. 15 (d) is a residual life prediction result map of the method M3.

FIG. 16 is a graph showing the comparison of the actual and predicted values of the remaining life of the bearing 2-2.

Fig. 17 is an error comparison chart of different prediction methods, fig. 17 (a) is an RMSE value comparison chart of different prediction methods, fig. 17 (b) is an MRE value comparison chart of different prediction methods, fig. 17 (c) is a SCORE value comparison chart of different prediction methods, and fig. 17 (d) is an MSE value comparison chart of different prediction methods.

Detailed Description

In order to make the technical problems, technical schemes and beneficial effects to be solved more clear, the invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.

As shown in fig. 1 and 2, the method for predicting the residual life of the multi-feature index bearing based on random correlation uses the 3 indexes to perform feature selection.

1) Monotonicity characterizes the consistency of system performance degradation, and the monotonicity Mon (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For characteristic sequences, ε (Y _i -Y _i-1 ) Is Y _i -Y _i-1 Is a unit step function of [0,1 ]]Within the range;is epsilon (Y) _i -Y _i-1 )-ε(Y _i-1 -Y _i ) From the summation of i=2 to n,is->Mon (Y) is the monotonicity of the feature sequence Y, which is a value of [0,1]In the range, the closer the monotonicity is to 1, the better the monotonicity trend of the characteristic is, and the performance degradation modeling and the residual life prediction can be better performed.

2) The trend index reflects the linear correlation degree between the characteristic sequence of performance degradation and the service life (time), and represents the change trend of the performance degradation curve along with the time to a certain extent, wherein the trend Corr (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For the feature sequence, n represents the total number of monitors throughout the performance degradation process,is Y _i t _i Summation from i=1 to n, +.>Is->Corr (Y) is the trend of the characteristic sequence Y, ++>For->Open square calculation was performed, t= (T ₁ ,t ₂ ...t _n ) For the corresponding sequence of monitoring instants, the value is 0,1]Within the range; the closer to 1 is the higher the correlation of the sequence of features with the run time, and the better the degradation of the device can be characterized for the features.

3) The robustness index characterizes the robustness of the system performance degradation characteristic to external interference, and the robustness Rob (Y) of the characteristic sequence Y is defined as:

wherein,the trend sequence for the corresponding performance degradation characteristic can be obtained by a smoothing technology, and the value range is still 0,1]The smoother the characteristic over time, the greater the robustness index, the higher the accuracy of the performance degradation modeling and rest life prediction, and +.>For->Obtaining an exponential function value->To pair(s)The function is divided by n from the sum of i=1 to n.

Based on characteristics including monotonicity, trending, and robustness, the choice of multiple features can be seen as a weighted sum, calculated as:

λ＝[φ ₁ Mon(Y)+φ ₂ Tre(Y)+φ ₃ Rob(Y)] (4)

wherein lambda is E [0,1 ]]；φ _i (i=1, 2, 3) is the eigenvalue Y among three choicesThe relative importance of the index is given by engineering experience. The relation of lambda and monotonicity, correlation and robustness is positive, and when the lambda value of a certain feature is higher, the feature is shown to have better comprehensive performance, and the degradation process of equipment performance can be reflected better.

Degradation modeling and parameter estimation:

because the state space model can represent various correlations between the multi-feature and the multi-index, the multi-index system constructed by the multi-feature is subjected to degradation modeling by adopting the state space model, Y is set as a feature vector of a sensing signal capable of reflecting equipment degradation, and X is set as an index vector reflecting degradation.

Y＝F(X)+v (5)

Wherein v is the comprehensive reflection of errors caused by sensing noise and feature extraction; f (X) is a functional relation between a sensing signal characteristic vector and an index vector reflecting degradation, in general, a specific form of F (X) is difficult to express in an analytic mode, and can be modeled by adopting a machine learning method, such as a neural network, bayesian estimation and the like, but the machine learning method needs to obtain a large number of data samples. In practical problems, the degradation index cannot be directly measured and can only be expressed by the feature vector of the sensing signal, if modeling is performed by using the model (5), the method is theoretically feasible, but the machine learning method is difficult to apply because a large amount of (X, Y) sample data cannot be obtained. If the extracted features are used to express the degradation indicator directly, this is present in many cases, where y=x+v, where C is the measurement coefficient, may also be y=cx+v. Due to the complex relationship between the plurality of features and the plurality of indicators. Therefore, the actual situation is more reflected by y=cx+v.

For example, there are multiple correlations between multi-features and multi-metrics, which may be one-to-one, many-to-one, and many-to-many correspondence. Thus, in constructing the multiple performance indicators, it is desirable to consider that a certain feature or set of features of the sensing signal reflects the degradation state of the device, and thus, it is possible to obtain:

wherein,for the characteristic value extracted for the sensor signal i reflecting the degradation of the device at time k,To reflect the index vector of system degradation, +.>C for comprehensively reflecting errors brought by sensor noise and feature extraction at k moment i _ip (p=1. N. n is characterized by i and the corresponding measurement coefficients between the indices p. When i=1..m, it is possible to obtain:

in actual engineering, degradation processes such as abrasion of aircraft components, simulated light emitting diodes and bearings can be modeled by using a Wiener process. Therefore, for the degradation process of the multi-index system, which is assumed to be a wiener process, the degradation process of the index is considered to depend on the degradation state of the index, and is related to the degradation states influenced by other indexes, so that the following can be obtained:

wherein a is _ji (i, j=1..n) is a correlation coefficient between the index i and the index j;as the diffusion coefficient of index i, Δt _k-1 As a nonlinear function, B _k -B _k-1 For the increment of brownian motion at time k, when i=1..n, it is available:

wherein a is _ji (i, j=1..n) is a correlation coefficient between the index i and the index j.

Because the different sensor data are the system responses of the same component degradation process and are interfered by the same external environment, certain correlation exists between the observation errors, and the correlation exists between the extracted multiple features, and the correlation is embodied in the covariance of the measurement noise. The process is carried out by the steps of,from the formulae (7) and (9):

wherein, xi _k Noise being a state of degradation of the system, i.e. ζ _k ～MVN(0,Q _k ) And corresponding varianceτ is the sampling interval, ω _k For measuring noise, i.e. omega _k ～MVN(0,R _k ) And corresponding variance->A _k-1 Is a correlation coefficient matrix, i.e.)>C is the measurement coefficient matrix, i.e.)>

Correlation among multiple indexes passes through a correlation matrix A of a state space model _k-1 To reflect it. From the perspective of the correlation between the degradation states, the modeling framework can intuitively characterize the correlation between the indices. For multi-index correlation modeling of multi-feature construction, a certain data sample needs to be obtained to perform parameter estimation and RUL prediction. Thus, under constructionThe multi-characteristic in the mode process is used as a measured value Y, and the degradation state X of the multi-index is estimated by a filtering algorithm. For a random linear system, the system can be used to measure Y _k And taking into account the characteristics of noise acting on the system, the state of degradation X of the system can be estimated by constructing a corresponding dynamic system _k . When the system is observable, the system can perform state estimation by adopting a KF algorithm. The conditions that the judging system can observe are as follows:

1) When Deltat _k-1 =τ, at this timeWhen->The system is observable. Wherein rank (·) is rank for the matrix.

2) When Deltat _k-1 As a nonlinear function, as shown in table 1:

TABLE 1 different nonlinear functions

For a state space model containing nonlinear degradation functions, it is necessary to construct an observable M when the system is determined to be observable _ds The (0, k) matrix is:

wherein phi is the transfer matrix and phi (j, 0) =a _j-1 φ(j-1,0)＝A _j-1 A _j-2 .......A ₁ A ₀ ；R ^-1 To invert the matrix R. Phi (j, 0) ^T Is the transpose of the matrix phi (j, 0),phi (j, 0) ^T C(j) ^T R ^-1 (j) C (j) Φ (j, 0) is the sum of j=0 to k.

When M _ds (0,k)Is positive matrix, i.e. M _ds Determinant of (0, k) is greater than 0, |M _ds If (0, k) | is not equal to 0, the system is observable, and the measurement data set Y is arranged at the moment k _0:k ＝{Y ₀ ,Y ₁ ,...,Y _k When the system is able to observe, use the measurement Y _k Estimating state X _k ；

1. Due to the statistical properties of the initial state:

where E {. Cndot. Is mathematical expectation and cov {. Cndot. Is variance.

Y can then be obtained ₀ Is the mathematical expectation E { Y } ₀ Covariance function v _Y (k ₀ ) And X ₀ Is a cross covariance matrix sigma of (2) _XY (k ₀ ) The method comprises the following steps:

wherein C is ^T Is the transpose of the measurement coefficient matrix C. cov { Y ₀ ,Y ₀ Y is } is ₀ And Y is equal to ₀ Covariance function between cov { X ] ₀ ,Y ₀ X is } is ₀ And Y is equal to ₀ Covariance function between.

Thus, time Y is 0 ₀ For X ₀ Estimate of (2)The method comprises the following steps:

K ₀ ＝P ₀ C ^T [CP ₀ C ^T +R ₀ ] ^-1 (15)

wherein K is ₀ A filtering gain of 0 time instant. E { X ₀ |Y ₀ Is at Y } ₀ Find X under the value ₀ Is a mathematical expectation of (a). And for X ₀ Covariance P of estimation error _0|0 The method comprises the following steps:

wherein,is->And->Covariance function between cov { X ] ₀ ,X ₀ |Y ₀ Is at Y } ₀ Find X under the value ₀ And X is ₀ Covariance function between.

2. At time k-1, it is assumed that measurement data Z has been obtained _k-1 ＝Y _0:k-1 ＝{Y ₀ ,Y ₁ ,...,Y _k-1 }。

Based on Z _k-1 For state X _k-1 Estimate of (2)Covariance P of estimation error _k-1|k-1 The method comprises the following steps of:

wherein cov { X ] _k-1 ,X _k-1 |Z _k-1 Is at Z } _k-1 Find X under the value _k-1 And X is _k-1 Covariance function between.For matrix->And->Is a mathematical expectation after transposed product.

Using information Z _k-1 For state X _k Estimate of (2)Namely:

wherein E { X } _k |Z _k-1 Is at Z } _k-1 Find X under the value _k Corresponding prediction errorThe method comprises the following steps:

corresponding estimation errorAnd->Covariance matrix P between _k|k-1 The method comprises the following steps:

reuse Z _k-1 For Y _k Forecast to obtain

Wherein E { Y } _k |Z _k-1 Is at Z } _k-1 Under the value, find Y _k Corresponding prediction errorThe method comprises the following steps:

corresponding estimation errorAnd->The synergetic difference matrix sigma _Y (k|k-1) and +.>And X is _k Sigma between _XY (k|k-1) is:

wherein cov { Y ] _k ,Y _k |Z _k-1 Is at Z } _k-1 Under the value, find Y _k And Y is equal to _k A covariance function between the two values,for matrix->And->Mathematical expectation after transposed product of (a)，Is X _k And->Covariance function between.

For sigma _XY (k|k-1), further obtainable:

∑ _XY (k|k-1)＝cov{X _k ,Y _k |Z _k-1 } (25)

wherein cov { X ] _k ,Y _k |Z _k-1 Is at Z } _k-1 Find X under the value _k And Y is equal to _k Covariance function between.

3. At time k, a new measurement Y is again obtained _k Reuse { Z ] _k-1 ,Y _k Z, i.e _k For X _k Performing state estimation, based on the nature of the innovation sequence, equivalent toFor X _k Make an estimation, then stateful estimation +.>The updated formula of (2) is:

wherein,is at->Under the condition of finding X _k K, K _k For filtering at time kWave gain, assumed estimation error->The method comprises the following steps:

then the error is estimatedAnd->Covariance matrix P of (2) _k|k The method comprises the following steps: />

Wherein,to at Z _k-1 And->Find X under the value _k And X is _k Covariance function between.

It is finally possible to obtain the product,using θ as an unknown parameter in the degradation model, i.e., θ= [ A, C, Q, G, X _0|0 ,P _0|0 ]The invention adopts an EM algorithm to solve the problem of maximum likelihood estimation of parameters when unknown states exist. And aiming at unknown parameters in the degradation model, on the basis of the online updating degradation state process of the KF algorithm, carrying out self-adaptive estimation on the unknown parameters by using a maximum Expectation (EM) algorithm and a smoothing filter (RTS) algorithm.

Wherein,for iteration of parameter estimation at step l+1 at time k, lnf (X _k ,Y _1:k The method comprises the steps of carrying out a first treatment on the surface of the θ) is a logarithmic function of the joint probability density of the system degradation state and the measured data, argmax (·) is a maximum value of the function.To at Y _1:k And a log likelihood function lnf under θ (X _k ,Y _1:k The method comprises the steps of carrying out a first treatment on the surface of the θ) is provided.

1) Expected calculation of log likelihood function:

wherein logf (X) _i |X _i-1 θ) is in X _i-1 And X under the condition of theta _i Log likelihood function of (1), logf (Y) _i |X _i θ) is in X _i And Y under the condition of theta _i Is a function of the log-likelihood of (a),is logf (X) _i |X _i-1 ,θ)+logf(Y _i |X _i θ) function from i=1 to k.

The expectation is calculated: v=e [ lnf (X _k ,Y _1:k ；θ)]

Wherein ln|P _0|0 I is P _0|0 The matrix tr (·) is the trace used for the matrix,is a matrix (X) _j -A _j-1 X _j-1 ) ^T Q ^-1 (X _j -A _j-1 X _j-1 ) Summation from j=1 to k.

The method further comprises the following steps:

in order to calculate the conditional expectation of the log-likelihood function, the expectation is found using the RTS smoothing filter algorithm:

wherein,is the mean value of backward filtering, P _j|k Variance of backward filtering, P _j,j-1|k Covariance matrix for backward filtering, and M _j Is the gain of the backward filtering. Since the hidden variables in equation (33) are many, it cannot be directly maximized. For this purpose, the invention adopts RTS smoothing algorithm to process the formula (33) to obtain:

2) And (3) deriving:

wherein,parameter estimation value for step l+1 +.> For->The function is maximized.

Let the partial differentiation of the v function be 0 to estimate, and obtain:

5. residual life prediction

For complex multi-feature and multi-index systems, three failure modes, namely competition failure, redundancy failure and fusion failure, exist, and the corresponding prediction results of the residual life of the system under different failure modes are different. Before solving for the system RUL in three failure modes, the degradation state distribution X of the multi-index system is firstly solved _k+l ：

Wherein,for matrices φ (k+l, j+1) and ζ _j+1 Summation from j=k to k+l-1 after the product.

Further obtain X _k+l Is not limited by the desire of (a)Sum of variances P _k+l ：

Let, S (l) _k )＝X(l _k +t _k )-X(t _k ) S (0) =0, i.eWherein S (l) _k ) From t for the degenerate state X _k To t _k+l Degradation increases over a period of time.

1) Failure of contention, failure of a system being ascertained as soon as one of the multiple degradation indicators exceeds a corresponding failure threshold, and corresponding remaining life L ₁ Is that

L ₁ ＝inf{l＞0:X ₁ (k+l)≥w ₁ ,or...,orX _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (43)

Wherein X is _i In (k+l), k+l is time t _k+l I is the number of indexes, w _i A fault threshold preset for index i. inf {.cndot. } is the lower bound of the function.

Residual life distribution function of system in contention failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )＜w ₁ ,orX ₂ (l _k +t _k )＜w ₂ ,...,orX _n (l _k +t _k )＜w _n |X _k ,θ _k Is at X } _k And theta _k One or more indexes in the degradation states of n indexes under the condition exceed the corresponding failure threshold value X ₁ (l _k +t _k )＜w ₁ ,orX ₂ (l _k +t _k )＜w ₂ ,...,orX _n (l _k +t _k )＜w _n Is a probability of (2).At the current time t as index i _k Is used to determine the amount of degradation of (1),for n index degradation states at t _k A joint probability density function of the time increment S,is->N-fold integration of the function.

2) Redundancy failure: failure of system and corresponding remaining life L ₂ Is that;

L ₂ ＝inf{l＞0:X ₁ (k+l)≥w ₁ ,...,and X _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (45)

residual life distribution function of system in redundancy failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n |X _k ,θ _k Is at X } _k And theta _k The degradation states of n indexes all exceed the corresponding failure threshold value X under the condition ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n Is a probability of (2).

3) Fusion failure, system failure and corresponding residual life of L ₃ 。

L ₃ ＝inf{l＞0:m ₁ X ₁ (k+l)+...+m _n X _n (k+l)≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n |X _k ,θ _k }

(47)

Wherein W is a fusion weight, w=λ ₁ w ₁ +..λ _i w _i +...λ _n w _n ；λ _i The weight coefficient corresponding to the index i.

Residual life distribution function of system under fusion failure mode according to first time definitionThe method comprises the following steps:

wherein P { m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n |X _k ,θ _k Is at X } _k And theta _k Under the condition that all indexes do not exceed the corresponding threshold value, but the degradation state weighted summation of n indexes exceeds the total failure threshold value m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)＜w ₁ ....and X _n (k+l)＜w _n The corresponding probabilities.

Then, the obtained values are obtained by using the formulas (44), (46) and (48)Obtaining RUL probability density function of the system through a finite difference formula:

wherein,for the purpose of distributing the function->And (5) deriving.

Since the degradation state of the equipment is implicit, the current time t can be obtained according to a full probability formula _k RUL probability density function of system in different failure modesThe method comprises the following steps: />

Wherein, psi (X) _k |θ _k ,Y _1:k ) Probability density function for implicit state of system, i.e Is state X _k And is obtained by a filtering algorithm.For->And psi (X) _k |θ _k ,Y _1:k ) Integration at- ≡to factafter the product.To be at theta _k ,X _k ,Y _1:k Residual life under conditions l _k Probability density functionIs a mathematical expectation of (a).

In order to further verify the validity of the residual life prediction method provided by the invention, the bearing test data on the PRONOSTIA test platform shown in figure 3 is utilized for verification analysis. Table 2 lists training data sets and test data sets for bearings under three conditions of load and rotational speed, which can be used to train and verify the RUL prognosis model. The measured data has a large error due to interference of noise and other factors in the bearing degradation process, as shown in fig. 4. Firstly, carrying out feature extraction on horizontal and vertical vibration signals of the bearings 1-4 by adopting a time domain and time frequency domain statistical method to form a candidate feature set, and carrying out feature selection through evaluation indexes such as monotonicity, trend, robustness and the like. And then constructing two performance indexes of horizontal vibration and vertical vibration representing the health state of the test device by using the selected characteristics, and selecting data of a degradation stage as test data. The invention adopts the degradation model in the formula (10) to describe the implicit degradation track of the bearing and carries out residual life prediction. Wherein fig. 5 shows the three-layer wavelet packet decomposition subband energy ratio of the bearing1-4 vibration signal. The time to failure of the bearings 1-4 for the first time was 2278 minutes.

Table 2 bearing operating conditions

Feature extraction is performed by using statistics such as time domain feature maximum value, average value, peak-to-peak value, variance, standard deviation, mean square amplitude, skewness, pulse factor, peak factor and the like, as shown in fig. 6; for the time-frequency domain signals, 8 vibration signals with the largest frequency band energy in the bearings 1-4 are selected for feature extraction, see fig. 7; and selecting the extracted characteristics through three evaluation indexes such as monotonicity, trend, robustness and the like, wherein the weight of the characteristic evaluation index is phi ₁ ＝0.3,φ ₂ ＝0.3,φ ₃ ＝0.4。

The features are ranked according to the composite index score, and features with composite indices greater than 0.5 will be retained. Finally, the eigenvalues of each group are selected as shown in table 3. Therefore, the system can observe, so that 11 characteristics are utilized, and the degradation state estimation is carried out by adopting the method provided by the invention, so that the performance index 1 and the performance index 2 of the bearing1-4 are constructed. The actual degradation data of the performance index constructed by the multiple eigenvalues in the vibration signals of the bearings 1-4 are shown in fig. 8. As seen from Table 4, the index 1 degradation function ln (b ₃ t _k )-ln(b ₃ t _k-1 ) The corresponding AIC value is the smallest and the index 2 degradation functionThe corresponding AIC value is the smallest. Δt of the two-index degradation model of the experiment _k-1 Respectively adopts ln (b) ₃ t _k )-ln(b ₃ t _k-1 ) And->The performance degradation failure thresholds of the bearings 1-4 index 1 and index 2 are 0.98 and 1.

Table 3 evaluation results of different characteristic data of vibration signals in case

TABLE 4 selection of four different degradation functions for multiple indicators

For bearings 1-4, it is assumed that only one of the performance degradations for index 1 and index 2 exceeds the corresponding threshold, and the system fails. The weights are 0.5 and 0.5, respectively. After obtaining the estimated values of the model parameters at each moment, the probability density function of the residual life of the bearing at different moments can be obtained by combining the formula (50), as shown in fig. 9. Taking 300 to 360 monitoring moments as an example, the comparison of the predicted value of the remaining life with the actual value of the three methods is shown in fig. 10. The health index method constructed by only RMS features is marked as an M0 method; the method of fusing a plurality of characteristics into one health index is denoted as an M1 method; the method of various correlations between the various features and the various indices (multi-index single degradation model) proposed by the method is denoted as the M2 method. The method presented herein (multi-index multiple degradation model) is denoted as the M3 method. In addition, four different prediction methods were compared herein using the four evaluation indicators in table 4. The smaller the (root mean square error) RMSE, (average relative error) MRE and (mean square error) MSE are, the better the fitting effect for the model is; the scoring function (Score) is mainly used to describe the relation between the predicted value and the true value overestimation and underestimation. Generally, the smaller the Score value, the more accurate the prediction result.

It can be seen from fig. 10 that for bearings 1-4, the M3 method is close to the true value update with a small relative error. The MSE of the M3 method is smaller compared to the M0, M1 and M2 methods. The analysis of the reason is that the M0 method only considers the mean square amplitude characteristic signal, which cannot fully reflect the degradation information of the bearings 1-4, and can lead to inaccurate RUL prediction. For the M1 method, although a plurality of characteristics are considered, the method only uses simple fusion to construct a composite HI index, and does not consider the correlation between the characteristics, so that the RUL prediction of the system is inaccurate. Although the M2 method builds multiple indices based on multiple features, the degradation model of the multiple indices is considered the same. It can also be seen from Table 5 that as the monitoring time increases, the predicted RUL for the M3 method is closer to the true value, with the MSE value being minimal compared to the M0, M1 and M2 methods. In addition, as can be seen from fig. 11, the relative error of the method of the present invention is smaller than that of methods M0, M1, M2, etc., so that the prediction of the method of the present invention can obtain higher prediction accuracy. In summary, for complex multi-index system RUL prediction, the three methods M0, M1 and M2 are not applicable. Therefore, the degradation model provided by the invention can more accurately reflect the degradation trend and the health state of the equipment in the online prediction framework.

TABLE 5 prediction of RUL and MSE for Bearing1-4 at different monitoring times

Feature extraction is performed by using statistics such as time domain feature maximum value, average value, peak-to-peak value, variance, standard deviation, mean square amplitude, skewness, pulse factor, peak factor and the like, see fig. 12; for the time-frequency domain signals, 8 vibration signals with the largest frequency band energy in the bearing2-2 are selected for feature extraction, see fig. 13. Wherein the weight of the characteristic evaluation index is phi ₁ ＝0.3,φ ₂ ＝0.3,φ ₃ ＝0.4。

Firstly, the proposed features are subjected to feature selection through indexes such as monotonicity, tendencies, robustness and the like, then the features are ranked according to the comprehensive index score, the features with the comprehensive index being more than 0.5 are reserved, and finally 7 feature data are selected. Normalizing the 7 selected characteristic data, and obtaining the M through calculation _ds (0,k)| _2×2 Not equal to 0. Therefore, the system can observe, 7 characteristics are utilized, and the construction of the Bearing2-2 performance degradation index is carried out according to the model provided by the invention, and the construction result is shown in fig. 14. The performance degradation failure thresholds of the bearing2-2 index 1 and the index 2 are 1 and 0.97. Based on the structural characteristics of the bearing system, it is assumed that in the contention failure mode, index 1 and index 2 occur only if one exceeds the threshold systemFailure occurs.

By AIC information criterion evaluation, delta t of two index degradation model of the experiment _k-1 Respectively adoptAndfor the bearing component 2-2, the last measurement is typically set to its failure threshold value [4, 18]. In order to verify the effectiveness and feasibility of the multi-index degradation modeling and residual life prediction method provided by the invention, the multi-index degradation modeling and residual life prediction method is compared with different multi-index prediction methods in a competition failure mode. First, degradation state estimation and unknown parameter estimation are performed, then, residual life prediction is performed on the system by using the estimated degradation state and parameters, and the residual life prediction results of the four methods are shown in fig. 15 by taking 300 th to 440 th data as an example.

As can be seen from fig. 16, as the monitoring time increases, the mean value of the remaining life predicted based on the M3 method gradually approaches the true value as compared with the M0, M1, and M2 methods. As shown in fig. 17, the prediction error of the M3 method among the four evaluation indexes is smaller than that of the other three methods.

By comparing with different advanced methods, the proposed model can intuitively reflect the correlation between indexes, provide RUL probability distribution, well predict the degradation trend of industrial equipment and improve the RUL prediction performance. In addition, the present invention contemplates correlation between various features and indicators. Because different indexes reflect different severity of system degradation, it is significant to construct multiple degradation indexes using multiple features.

The foregoing description of the preferred embodiment of the invention is not intended to be limiting, but rather is intended to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.

Claims (4)

1. The multi-feature index bearing residual life prediction method based on random correlation is characterized by comprising the following specific steps of:

step one, constructing a state space model of multi-feature degradation, selecting monotonicity, robustness and trend of multi-feature data through the state space model, and analyzing the relation between the multi-feature data and the multi-index data;

step two, considering the correlation among multiple features, the correlation among multiple features and multiple indexes, the correlation among multiple indexes and the degradation characteristic of a reflection system in a state space model, considering multiple degradation modes in the state space model, and carrying out degradation state and parameter joint estimation by adopting a maximum expectation algorithm and a Kalman filtering method;

and step three, defining a redundancy failure mode, a fusion failure mode and a competition failure mode among multiple indexes, and establishing a residual life prediction model by using a state space model.

2. The method for predicting residual life of a multi-feature index bearing based on random correlation as claimed in claim 1, wherein in step one, feature selection is performed by adopting monotonicity, robustness and trend:

1) Monotonicity characterizes the consistency of system performance degradation, and the monotonicity Mon (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For characteristic sequences, ε (Y _i -Y _i-1 ) Is Y _i -Y _i-1 Is a unit step function of [0,1 ]]Within the range;is epsilon (Y) _i -Y _i-1 )-ε(Y _i-1 -Y _i ) From the summation of i=2 to n,is->Mon (Y) is the monotonicity of the feature sequence Y;

2) The trend index reflects the linear correlation degree between the characteristic sequence of performance degradation and the service life, and the trend Corr (Y) of the characteristic sequence Y is defined as:

wherein y= (Y) ₁ ,Y ₂ ..Y _n ) For the feature sequence, n is the total number of monitoring in the whole performance degradation process, t= (T ₁ ,t ₂ ...t _n ) For the corresponding sequence of monitoring instants, the value is 0,1]Within the range;is Y _i t _i From the summation of i=1 to n,is->Corr (Y) is the trend of the characteristic sequence Y,for->Performing open square calculation;

3) The robustness index reflects the robustness of the system performance degradation characteristic to external interference, and the robustness Rob (Y) of the characteristic sequence Y is defined as:

wherein,the value range of the trend sequence of the performance degradation characteristic is [0,1 ]]，To pair(s)Obtaining an exponential function value->For->The function divides the summation of i=1 to n by n, rob (Y) is the robustness of the feature sequence Y, and based on characteristics including monotonicity, trending and robustness, the choice of multiple features is converted into a weighted summation, the calculation formula λ is:

λ＝[φ ₁ Mon(Y)+φ ₂ Tre(Y)+φ ₃ Rob(Y)] (4)

wherein lambda is E [0,1 ]]，φ _i (i=1, 2, 3) is the relative importance of the eigenvalue Y on three selection criteria, λ is a positive correlation relationship with monotonicity, correlation and robustness.

3. The method for predicting the residual life of the multi-feature index bearing based on the random correlation according to claim 2, wherein in the second step, a state space model is adopted to conduct degradation modeling on a multi-feature constructed multi-index system, extracted features are used for directly expressing degradation indexes, and the features of sensing signals reflect degradation states of equipment;

the correlation among the multiple indexes is reflected by a correlation matrix of the state space model, and when the system is observable, the system adopts a KF algorithm to perform state estimation.

4. The method for predicting residual life of multi-feature-index bearing based on random correlation as recited in claim 3, wherein in step three, a degradation state distribution X of the multi-index system is obtained _k+l The method comprises the following steps:

wherein,for matrices φ (k+l, j+1) and ζ _j+1 The sum of the products from j=k to k+l-1 gives further X _k+l Is->Sum of variances P _k+l ：

Let, S (l) _k )＝X(l _k +t _k )-X(t _k ) S (0) =0, i.eWherein S (l) _k ) From t for the degenerate state X _k To t _k+l Degradation delta over a period of time;

1) Failure of competition

Judging that the system fails and the corresponding residual life L exists in a plurality of degradation indexes as long as one index exceeds the corresponding failure threshold ₁ The method comprises the following steps:

L ₁ ＝inf{l>0:X ₁ (k+l)≥w ₁ ,or...,orX _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (7)

wherein X is _i In (k+l), k+l is time t _k+l I isNumber of indexes, w _i A fault threshold value preset for the index i, wherein inf { · } is the lower bound of the function;

residual life distribution function of system in contention failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )<w ₁ ,orX ₂ (l _k +t _k )<w ₂ ,...,orX _n (l _k +t _k )<w _n |X _k ,θ _k Is at X } _k And theta _k In the degradation state of n indexes under the condition, one or more indexes exceed the corresponding failure threshold value X ₁ (l _k +t _k )<w ₁ ,orX ₂ (l _k +t _k )<w ₂ ,...,orX _n (l _k +t _k )<w _n Is a function of the probability of (1),at the current time t as index i _k Is used to determine the amount of degradation of (1),for n index degradation states at t _k A joint probability density function of the time increment S,is->N-fold integration of the function;

2) Failure of redundancy

When the multiple degradation indicators each exceed a corresponding failure threshold,then conclude that the system is malfunctioning and the corresponding remaining lifetime L ₂ The method comprises the following steps:

L ₂ ＝inf{l>0:X ₁ (k+l)≥w ₁ ,...,and X _i (k+l)≥w _i ,or...,i＝1,2,...,n|X _k ,θ _k } (9)

residual life distribution function of system in redundancy failure mode according to first time definitionThe method comprises the following steps:

wherein P { X ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n |X _k ,θ _k Is at X } _k And theta _k The degradation states of n indexes all exceed the corresponding failure threshold value X under the condition ₁ (l _k +t _k )≥w ₁ ,and X ₂ (l _k +t _k )≥w ₂ ,...,and X _n (l _k +t _k )≥w _n Probability of (2);

3) Failure of fusion

When all degradation indexes do not exceed the corresponding threshold value, but the weighted summation of a plurality of degradation indexes exceeds the total failure threshold value, judging that the system fails and the corresponding residual life L ₃ The method comprises the following steps:

L ₃ ＝inf{l>0:m ₁ X ₁ (k+l)+...+m _n X _n (k+l)≥Wand X ₁ (k+l)<w ₁ ....and X _n (k+l)<w _n |X _k ,θ _k } (11)

wherein W is a fusion weight, w=λ ₁ w ₁ +..λ _i w _i +...λ _n w _n ，λ _i Is a weight coefficient corresponding to index i, is in a fusion failure modeResidual life distribution function of systemThe method comprises the following steps:

P{m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)<w ₁ ....and X _n (k+l)<w _n |X _k ,θ _k is at X } _k And theta _k Under the condition that all indexes do not exceed the corresponding threshold value, but the degradation state weighted summation of n indexes exceeds the total failure threshold value m ₁ X ₁ (l _k +t _k )+...+m _n X _n (l _k +t _k )≥Wand X ₁ (k+l)<w ₁ ....and X _n (k+l)<w _n The corresponding probabilities;

then, the distribution function obtained by the formulas (8), (10) and (12) is usedObtaining RUL probability density function of the system through a finite difference formula:

wherein,for the purpose of distributing the function->Seeking a derivative;

obtaining the current time t according to the full probability formula _k Considering the remaining life of an implicit state system under different failure modesRate Density functionThe method comprises the following steps:

wherein, psi (X) _k |θ _k ,Y _1:k ) Probability density function for implicit state of system, i.e Is state X _k Mean and variance of (a) and obtained by a filtering algorithm, < + >>For->And psi (X) _k |θ _k ,Y _1:k ) Integration at- ≡to factafter the product,to be at theta _k ,X _k ,Y _1:k Residual life under conditions l _k Probability density functionIs a mathematical expectation of (a).