Method for estimating the remaining service life of the equipment under test

Abstract

A method for estimating the remaining useful life of a test equipment has a preliminary phase comprising the following steps: collecting test observations (step 10) and generating at least one signed test time series (step S x ); dividing the test time series to obtain severity categories corresponding to the aging stages of the test equipment equipment (step 14); performing initial learning of the diagnostic model on the test equipment equipment (step 45); performing second learning of the signature prediction model (step 51); the working stage includes the following steps: collecting observations when working on the test equipment and using the prediction model to generate an extrapolated time series; using the diagnostic model to classify the extrapolated time series and derive the remaining useful life of the test equipment.

Description

[0001] The present invention relates to the field of health monitoring processes. Background Art

[0002] The process of health monitoring, often referred to as PHM (for Prognostic and Health Management) but less commonly known as condition-based maintenance, predictive maintenance, or predictive maintenance, is used to predict the degradation of equipment or systems, particularly in aviation.

[0003] The health monitoring process thus makes it possible to avoid breakdowns, optimize service life, predict and plan maintenance operations (and possible disassembly) and thus reduce the costs of equipment or system maintenance operations.

[0004] The health monitoring process includes five main steps: data acquisition, data processing (or preconditioning), signature extraction, diagnosis and prognosis.

[0005] The data collection step includes collecting data generated by sensors that measure parameters representing the performance of the equipment or system. For example, data collection is performed via testing or downloading to the equipment or system while the equipment or system is running.

[0006] The data preconditioning step involves filtering the data to remove measurement noise from the data.

[0007] The signature extraction step includes processing the data to extract one or more signatures that are indicative of the state of the equipment or system, and one or more signatures that may be indicative of a degradation in the operation of the equipment or system.

[0008] The diagnostic step consists in determining the current health state of the equipment or system based on a set of fault categories (failure modes), each fault category being associated with a severity level.

[0009] The prognostic step involves estimating the remaining useful life (RUL) of the equipment or system.

[0010] Purpose of the Invention

[0011] The object of the present invention is to reliably and accurately estimate the remaining useful life of a piece of equipment under test. Summary of the Invention

[0012] To achieve this goal, a method for estimating the remaining useful life of the equipment under test is proposed. The method has a preliminary stage, which includes the following steps:

[0013] Collect test observations from test equipment similar to the equipment under test;

[0014] generating at least one signed test time series based on the test observation;

[0015] Dividing the test time series to obtain severity categories corresponding to aging stages of the test equipment;

[0016] Conduct first learning of diagnostic models on test equipment;

[0017] performing a second learning of the signature prediction model;

[0018] And the working phase includes the following steps:

[0019] Collect observations while working on the equipment under test;

[0020] Using the predictive model to generate an extrapolated time series representing the evolution of the signature on the test equipment;

[0021] This extrapolated time series is classified using a diagnostic model;

[0022] Calculating the membership function of the extrapolated time series to severity categories;

[0023] Based on the membership function, the remaining service life of the tested equipment is derived.

[0024] The estimation method according to the invention thus implements, during a preliminary phase, for example carried out in a laboratory, a first learning of a diagnostic model and a second learning of a predictive model for at least one signature derived from observations made on a test equipment device.

[0025] Subsequently, during the operational phase, the equipment under test is observed in real time, and at the "current" time t, the evolution of the signature on the equipment under test is extrapolated due to the predictive model and subsequently classified due to the diagnostic model. The severity class membership function then provides a reliable and accurate estimate of the remaining useful life of the equipment under test, based on the behavior of test equipment devices similar to the equipment under test under similar stresses.

[0026] An estimation method such as that just described is also proposed, in which the estimated remaining useful life is such that:

[0027]

[0028] wherein q is the estimated end-of-life time and t is the current time, and wherein at the estimated end-of-life time q, the last severity category membership function becomes greater than the penultimate severity category membership function.

[0029] An estimation method such as the one just described has also been proposed, in which a temporal consistency criterion is used to define several severity classes resulting from the partitioning.

[0030] Additionally, an estimation method such as the one just described is proposed, in which the temporal consistency criterion is evaluated using a first metric defined as follows:

[0031]

[0032] where x k (i) is the i-th test observation collected for test rig k, and where ct is the temporal consistency between two consecutive test observations and is defined using the following formula:

[0033]

[0034] where y k (i) is the observation x k (i) severity category.

[0035] An estimation method such as the one just described is also proposed, in which a second metric (which is an energy metric) is evaluated to select a partitioning algorithm.

[0036] We also propose an estimation method as described above, where the energy metric is defined by:

[0037]

[0038] where Ω i is the test observation set of category i, is the test observation, g i is the centroid of the i-th severity category, d is the distance measure and is defined by:

[0039]

[0040] Among them A d is a positive semidefinite matrix defining the distance metric, and where the symbol “|||| d ” denotes the norm corresponding to the distance d.

[0041] An estimation method such as the one just described is also proposed, in which the partitioning algorithm is the fuzzy C-means algorithm.

[0042] An estimation method such as the one just described has also been proposed, in which the prediction model is trained using a multi-step iterative method.

[0043] An estimation method such as the one just described is also proposed, in which the prediction model is constructed based on nonlinear SVR regression.

[0044] In addition, an estimation method such as that described above is proposed, in which the first learning of the diagnosis model uses the fuzzy SVM-FA classification algorithm.

[0045] The invention will be better understood upon reading the following description of specific non-limiting embodiments of the invention.

[0046] BRIEF DESCRIPTION OF THE DRAWINGS

[0047] With reference to the accompanying drawings, in which:

[0048] Figure 1 shows a roller bearing mounted on a shaft rotated by an electric motor;

[0049] Figure 2 The various working steps of the diagnosis and prognosis process according to the present invention are shown;

[0050] Figure 3 The various working steps of the diagnosis and prognosis process according to the present invention are shown;

[0051] Figure 4 shows the partitioning of the time series of signatures;

[0052] Figure 5 is a plot of the global temporal consistency curve as a function of the number of categories for the four partitioning algorithms;

[0053] Figure 6 is a graph showing curves of energy criteria as a function of the number of classes for four partitioning algorithms;

[0054] Figure 7 It is a plot of the energy criterion versus the number of categories for the C-fuzzy mean partitioning algorithm;

[0055] Figure 8 is a graph showing vibration signal power versus time for four test rig devices, the curves being divided according to four severity categories;

[0056] Figure 9 is a graph showing vibration signal power versus time for four test rig devices, the curves being divided according to four severity categories;

[0057] Figure 10 is a graph showing vibration signal power versus time for four test rig devices, the curves being divided according to four severity categories;

[0058] Figure 11 is a diagram showing an example of discrete classification;

[0059] Figure 12 is a graph showing classification results and severity levels;

[0060] Figure 13 The steps of the learning algorithm for the prediction model are shown;

[0061] Figure 14 is a graph including a curve of signal power versus time for a test equipment device, the curve being obtained using a predictive model;

[0062] Figure 15 is a graph including a plot of relative mean absolute error versus time for predicting the power of a vibration signal on one of the test rig devices;

[0063] Figure 16 is a graph including a plot of relative mean absolute error versus time for predicting the power of a vibration signal on one of the test rig devices;

[0064] Figure 17 The steps of the prediction and classification algorithms are shown;

[0065] Figure 18 is a graph showing the accuracy and precision performance of the remaining useful life prediction algorithm;

[0066] Figure 19 This is a graph comparing the accuracy of two remaining useful life prediction algorithms;

[0067] Figure 20 is a graph showing the accuracy and precision performance of the remaining useful life prediction algorithm;

[0068] Figure 21 is a graph showing estimated severity class membership curves for the equipment under test;

[0069] Figure 22 is a graph including remaining useful life estimation curves from two different times for the equipment under test.

[0070] Detailed description of the invention

[0071] The present invention is implemented in the context of health monitoring and relates to a method for estimating the remaining useful life of equipment under test.

[0072] Therefore, the present invention relates more specifically to the "prognostic" stage of the health monitoring process.

[0073] Here, by “prognosis” and in accordance with ISO 13381, we mean: “an estimate of the operational duration until failure and the risk of the presence or subsequent occurrence of one or more failure modes”. By adopting a system-oriented interpretation that takes into account the interactions between system components, a component is considered to have reached the end of its lifetime when its damaged state causes or accelerates the damage of other components.

[0074] The term "equipment under test" should be interpreted broadly. The equipment under test may be any equipment, for example an LRU (line replaceable unit) such as an actuator, or a system comprising several equipment pieces interacting with each other, or an electrical or mechanical assembly of any complexity.

[0075] The equipment under test here is embodied in an aircraft, but the present invention is not limited to such applications.

[0076] In the following, it is assumed that the equipment under test has not been updated and has not undergone any maintenance actions to restore it to its previous healthy state. Therefore, its remaining useful life (RUL) decreases linearly with time t:

[0077] RUL=t EOL -t,

[0078] Among them, t EOL It is the end of the life of the tested equipment.

[0079] Depending on the application, time can be considered in different ways.

[0080] The time since commissioning can be calculated. The accumulated operating time can also be used. For electromechanical actuators, this corresponds to the phase when the actuator is powered on. In the context of aviation systems, where component traceability is very strict, both types of time information are available. The time considered can also be the number of operating cycles, considering a deterministic speed / load profile.

[0081] Here, the accumulated working time is used.

[0082] It is also considered that the equipment under test degrades only during operation. Indeed, for actuators of aviation systems, and in particular for actuators integrated into single-aisle aircraft, the waiting phases on the apron or in the hangar are short compared to the operating time.

[0083] The estimation process consists of several steps. These steps are performed in a preliminary phase and subsequently in a working phase. The preliminary phase is performed on test equipment similar to the equipment under test and is performed at least partially in a laboratory or design office. The working phase is performed while the equipment under test is operating.

[0084] Each step will be described theoretically and then illustrated with an application example.

[0085] refer to Figure 1 ,The example application is taken from a public database, namely the benchmark bearing aging database, ,from tests conducted at the Intelligent Maintenance Systems (IMS) Center at the University of Cincinnati.

[0086] This public data was provided by NASA (Lee J, Qiu H, Yu G, Lin J, Rexnord Technical Services). Bearing Data Set [Internet], IMS, University of Cincinnati, NASA Ames Prognostics Data Repository (http: / / ti.arc.nasa.gov / project / prognostic-data-repository), NASA Ames Research Center, Moffett Field, CA; 2007. They are available at: https: / / ti.arc.nasa.gov / tech / dash / pcoe / prognostic-data-repository / , Qiu et al., 2006; Lee et al., 2007).

[0087] Four roller bearings 1a, 1b, 1c, 1d are placed on a single shaft 2 which is rotated at 2000 rpm by a motor 3. A spring mechanism applies a constant radial load of 6000 pounds or 26698 Newtons.

[0088] refer to Figure 2 and Figure 3 The preliminary phase first comprises a step 10 of collecting test observations, which are here performed via testing on a test rig. In the example application, the test rig and the test rig are roller bearings 1. Note that they could just as well be ball bearings.

[0089] In n u Test equipment on the m u Observation.

[0090] In this application example, the observation is performed with the aid of an accelerometer 4 integrated in the bearing 1 .

[0091] The preliminary stage of the estimation process then comprises a step 11 of preconditioning the data.

[0092] The preliminary phase then comprises a signature extraction step 12 .

[0093] The vector of fault signatures is denoted by X:

[0094] X = [x1, x2, ..., xp], where p is the number of fault signatures. A signature is a quantity that is sensitive to degradation of the test equipment device (and the equipment under test).

[0095] The signatures used here are the power of the vibration signal and the RMS value of the vibration signal.

[0096] The signature is stored in the database 13 .

[0097] Then generate a time series signature. X The signature time series corresponds to the evolution of the signature over time.

[0098] The estimation process then includes an automatic partitioning step 14 of the data to create coherent categories. These categories are also called failure mode, failure degree or severity categories. The categories correspond to the aging stages of the test equipment equipment.

[0099] The goal of the division is to u m on the test equipment u The different severity defects in the database are identified in the observations. Therefore, it is necessary to u *n u Each observation corresponds to a vector with a shape of p parameters. Therefore, the test equipment device is represented by p time series (one time series for each shape vector signature). The resulting severity class should be consistent over time on all test equipment devices.

[0100] For each test equipment, each category will be represented by the start time t start and end time t stop The severity categories are arranged from least severe to most severe. Therefore, the t stop will be the next category of t start .

[0101] exist Figure 4 In the example in , four categories c1, c2, c3 and c4 are shown.

[0102] First, try to define the optimal number of categories.

[0103] The partitioning of time series classically includes three categories of problems: partitioning of the entire time series into clusters, subsequence clustering, and time point clustering.

[0104] Partitioning the entire time series involves treating each sequence as an object and forming groups of the entire series. This type of partitioning is equivalent to classical partitioning, but uses integer time series instead of shape vectors. Subsequence partitioning aims to identify subsequences that repeat within a single time series. Partitioning the time series points involves assigning categories based on temporal proximity and values ​​within a single series.

[0105] None of these categories corresponds exactly to the problem of the present invention. Partitioning the entire time series is not suitable, because the goal is to find a coherent sequence. Subsequence partitioning is also not suitable, because it works from a single sequence. The same applies to partitioning time points. The problem here is to deal with n related to each other. u *p sequences.

[0106] Therefore, a temporal consistency criterion for partitioning is defined and used to select the number of categories.

[0107] Here is a definition of temporal consistency of partitions. It is based on three conditions:

[0108] - The number of categories is the same for different test equipment devices;

[0109] - Each category is continuous in time: there are no "jumps" in time;

[0110] - For different test equipment devices, these categories follow each other in the same order.

[0111] The following first metric is proposed to evaluate the temporal consistency of observations from a test rig, which are sorted in increasing time:

[0112]

[0113] where x k (i) is the i-th observation measured for test equipment device k, and where ct is the temporal consistency between two consecutive observations and is defined using the following formula:

[0114]

[0115] where y k (i) is the observation x k (i) Category.

[0116] For a given test equipment device k, the resulting partitioning is for CT k =c-1 is perfect agreement, where c is the number of categories.

[0117] On the contrary, if the division is completely inconsistent, then CT = m u -1, where m u is the number of observations. To facilitate the explanation of the first measure and its use in comparisons, an equivalent version is defined using values ​​in [0,1]:

[0118]

[0119] For a completely time-consistent partition, we have:

[0120] CTN=0.

[0121] For completely inconsistent partitions, there are:

[0122] CTN=1.

[0123] For different test equipment, there are different number of categories:

[0124] CTN<0.

[0125] Now, the overall timing consistency across all test rig devices can be calculated using the following formula:

[0126]

[0127] The “|” is the absolute value symbol.

[0128] The first metric verifies the first two conditions, namely the number of identical classes of devices per test rig and temporal continuity, but not the third. However, the third condition can be verified by plotting the evolution of the signature of the partitioned observations over time.

[0129] A partitioning algorithm should also be selected.

[0130] For a given number of classes, if several partitioning algorithms are completely consistent, a second metric is evaluated to select the partitioning algorithm. The second metric is an energy metric that must be minimized to have the most compact classes possible. The second metric is calculated using the following formula:

[0131]

[0132] where Ω i is the set of observations of category i, is the observation, g i is the centroid of the i-th class, d is the distance measure and is defined by the following formula:

[0133]

[0134] Among them A d is a positive semidefinite matrix (chosen by the user) defining the distance metric d, and where the symbol “|||| d ” denotes the norm corresponding to the distance d.

[0135] Here, the Euclidean distance is used:

[0136] A d =I, the identity matrix.

[0137] Four algorithms classically used for unsupervised segmentation are evaluated: K-means, K-means, C-means fuzzy, and C-fuzzy means probabilistic. Figure 5 The four algorithms process the database for numbers of categories ranging from 2 to 10. For each number of categories and each algorithm, the normalized global temporal consistency is calculated: curve 16 corresponds to the K-means algorithm, curve 17 corresponds to the K-means algorithm, curve 18 corresponds to the C-means fuzzy algorithm, and curve 19 corresponds to the C-means fuzzy probability algorithm.

[0138] Note that for c > 6, all partitions are inconsistent regardless of the algorithm chosen. Note subsequently that only the fuzzy C-means algorithm gives consistent results between c = 4 and c = 6.

[0139] refer to Figure 6 To confirm this result, the second energy metric is evaluated for values ​​of c up to 6. Curve 20 corresponds to the K-means algorithm, curve 21 corresponds to the K-means algorithm, and curve 22 corresponds to the C-means fuzzy algorithm.

[0140] The curve corresponding to the C-means fuzzy probability algorithm is not included in this figure because it reaches values ​​that are too high. Clearly, regardless of the number of classes in the interval considered, the fuzzy C-means algorithm has the lowest energy. Therefore, this algorithm is selected.

[0141] Note that the C-fuzzy mean algorithm was developed as an improvement of the K-fuzzy algorithm. It is assumed that an observation can belong to several categories at the same time, but to different degrees.

[0142] refer to Figure 7 , we cannot choose too many categories here, as this would mean a significant increase in computational cost. Therefore, as with parameter selection, the choice of the number of categories will be made at the "elbow" of the criterion evolution. Therefore, here, we will consider cases with three, four, and five categories.

[0143] refer to Figures 8 to 10 , shows the results for c=4.

[0144] Figure 8 The signature of interest is the strength of the vibration signal. Severity category 1 corresponding to curve segment 25, severity category 2 corresponding to curve segment 26, severity category 3 corresponding to curve segment 27, and severity category 4 corresponding to curve segment 28 are distinguished for the four test rigs.

[0145] Figure 9The signature of interest is the strength of the vibration signal. Severity category 1 corresponding to curve segment 30, severity category 2 corresponding to curve segment 31, severity category 3 corresponding to curve segment 32, and severity category 4 corresponding to curve segment 33 are distinguished for the four test rigs.

[0146] It can be seen that for this number of categories, temporal consistency is maintained and the energy criterion converges to a reasonable level (reduced by two-thirds). The first category corresponds to the run-in phase. The end-of-life phase is well isolated. The living environment is divided into two categories.

[0147] refer to Figure 10 ,We can see that, for example, when the number of categories is c=8, using the same partitioning algorithm for selection will result in completely inconsistent category divisions.

[0148] Figure 10 The signature of interest is the intensity of the vibration signal. For the four test equipment devices, severity categories n°1, corresponding to curve section 35, severity category n°2, corresponding to curve section 36, severity category n°3, corresponding to curve section 37, severity category n°4, corresponding to curve section 38, severity category n°5, corresponding to curve section 39, severity category n°6, corresponding to curve section 40, severity category n°7, corresponding to curve section 41, and severity category n°8, corresponding to curve section 42, were distinguished.

[0149] Then, the preliminary stage of the estimation process according to the invention comprises a diagnosis step 44 (cf. Figure 2 and Figure 3 ), which first includes performing a first learning 45 of the severity diagnosis model.

[0150] The diagnosis consists in assigning each working point to the closest class for a given distance type. Therefore, a classification model must be trained to perform this task using data available in a dedicated database, said data coming from test equipment.

[0151] use Figure 11 The classification method shown in will be possible.

[0152] In this approach, the similarity between a new observation and classes c1, c2, c3 is measured, where class c1 corresponds to no defect, class c2 corresponds to the first defect, and class c3 corresponds to the second defect.

[0153] This classification is discrete and does not take into account the degree of membership in a given category.

[0154] However, using the fuzzy membership function, another classification degree is selected to evaluate the severity degree of each failure mode.

[0155] The first learning of the severity diagnosis model is based on a fuzzy SVM (space vector machine) classification algorithm with multi-class membership functions.

[0156] The diagnostic model is then implemented (step 46) and validated (step 47). The severity classification results are given in Figure 12 Visible in.

[0157] Then, the preliminary stage of the estimation process according to the invention comprises a prognostic step 50 (cf. Figure 2 and Figure 3 ), which first includes performing a first learning 51 of a severity prediction model.

[0158] Prognosis involves extrapolating the signature into a limited range which will be called the prediction range ( Figure 2 The extrapolated signature is classified in a diagnostic sense and its severity is assessed at a predefined prediction range. The remaining useful life can be derived from this extrapolation.

[0159] First, the phase of building a prediction model (also called a regression model) is described.

[0160] The learning method of the forecast model was chosen based on the so-called "multi-step" strategy, where forecasts are made at several time steps.

[0161] These predictions are based on a series of observations at past moments. They represent themselves in the form:

[0162]

[0163] in:

[0164] X t-n0→t is a sequence of n0 observations (up to t) of signature values. These observations may relate to a single signature, so X t-n0→t The coefficients of , or these observations may involve multiple signatures, so the coefficients will have where p is the number of signatures;

[0165] fp is the prediction model;

[0166] Θ is the parameter of the prediction model;

[0167] - is the estimated sequence up to the prediction horizon hp.

[0168] The number of observations and the forecast horizon are typically set by the user.

[0169] If all coefficients It is called univariate extrapolation. If all coefficients This is called multivariate extrapolation. The second case allows taking into account the correlations between signatures.

[0170] The multi-step iterative approach involves learning a model to make predictions at time t+1. Then, iteratively running the same model in series to make predictions at time t+hp.

[0171] The advantage of this approach is that it's easy to implement. The disadvantage is that errors accumulate as the predictions progress. Within the constraints of aviation, it's necessary to consider the size factor: the memory required to store each model. Given these constraints, this iterative approach seems like a good first approach.

[0172] Figure 13 The algorithm of this method is summarized in .

[0173] Now let's describe the linear regression model itself.

[0174] The goal of temporal regression is to determine the parameter vector Θ, which defines the fp prediction model. This model (also known as a learning model) is based on nonlinear SVR (Support Vector Machine for Regression) regression. The latter aims to provide a prediction model based on a set of m labeled data {(x1, y1), ... To predict the output model.

[0175] Therefore, the algorithm must find the function This function minimizes the error on the learning set while also generalizing well. The regression problem can be formulated using a support vector machine. This is based on an SVM with flexible margins. For nonlinear cases, the problem is reformulated using the kernel technique. The goal is to maximize the following function:

[0176]

[0177] Under the following constraints:

[0178]

[0179] 0≤α i ≤C

[0180] 0≤α i * ≤C

[0181] where K is a kernel function that can be chosen such that:

[0182] K(x,x')=tanh(ak.x.x'+b k),

[0183] And where (αi, αi*) are the primary and dual Lagrangian coefficients. Then, the prediction function is expressed as:

[0184]

[0185] The constant b can be obtained from the KKT (Karush-Kuhn-Tucker) condition applied to one of the support points (i.e., the points corresponding to the following):

[0186]

[0187] in:

[0188]

[0189] The prediction model is then implemented (step 52) and verified (step 53) on all signatures.

[0190] The above description applies to Figure 1 The roller or ball bearing 1 is shown.

[0191] SVR is through Statistics and Machine Learning Toolbox in TM (Statistics and Machine Learning Toolbox TM Three prediction start times were tested: a start time of t = sp1 = 1 day (where no observations from the test equipment were used for learning), a start time of t = sp2 = 10.45 days, which corresponds to the end of the run-in phase, and a start time of t = sp3 = 20.65 days, which corresponds to the start of the second-to-last severity category.

[0192] The prediction results are evaluated in terms of mean absolute percentage error (MAPE), which is the dual of cumulative relative accuracy.

[0193] The predicted results are comparable to those of the test equipment.

[0194] Figures 14 to 16 Test results for bearing 1 a are shown.

[0195] Figure 14 Prediction of the power evolution of vibration signals is involved.

[0196] Curve 60 corresponds to the starting point at t = sp1, curve 61 corresponds to the starting point at t = sp2, and curve 62 corresponds to the starting point at t = sp3. Curve 63 corresponds to the target, ie the actual development observed.

[0197] Figure 15 The relative mean absolute error (MAPE) of the predictions for the vibration signal strength is shown.

[0198] Curve 64 corresponds to the starting point at t=sp1, curve 65 corresponds to the starting point at t=sp2, and curve 66 corresponds to the starting point at t=sp3.

[0199] Figure 16 The relative mean absolute error (MAPE) of the predictions of the vibration signal strength is shown.

[0200] Curve 67 corresponds to the starting point at t=sp1, curve 68 corresponds to the starting point at t=sp2, and curve 69 corresponds to the starting point at t=sp3.

[0201] The error increases over time, consistent with the use of an iterative approach. Generally speaking, the greater the number of observations used for the test equipment, the smaller the error. However, the error increases dramatically in the final prediction. This is due to an exponential overestimation of the signature. However, this occurs after the unit's lifetime has expired (i.e., 28 days).

[0202] refer to Figure 17 , we will now describe the testing phase that was performed in order to experimentally verify the estimation process according to the present invention.

[0203] As previously mentioned, observations from four test rig devices are available. By using these observations and comparing the obtained results with the remaining useful lives actually observed on these test rig devices, the prediction results according to the present invention, and thus the effectiveness of the estimation process, are evaluated. One of the test rig devices then becomes the device under test.

[0204] The MDLC classification model as well as the MDLR prediction model are considered to have been trained.

[0205] The working phase includes the steps of collecting observations while working on the equipment under test. Generate a signature sequence.

[0206] Use the predictive model to generate an extrapolated time series representing the evolution of the signature on the test equipment.

[0207] Get the predicted signature from t+1 to t+hp.

[0208] Classify the extrapolated time series using a diagnostic model.

[0209] Membership functions are calculated for the time series extrapolated to each severity category. Based on these membership functions, the remaining useful life of the equipment under test is derived.

[0210] Membership function u j (j from 1 to c) indicates the assigned fault category and its severity level.

[0211] The prediction range is when the algorithm reaches performance α e (accuracy) and β e (accuracy) time t αeβe and end of life time t EOL The intervals between them are as follows:

[0212] PH=t EOL -t αeβe .

[0213] refer to Figure 18 , the criterion can be expressed in terms of the mean of the RUL distribution or the distribution itself. Then note that:

[0214] RUL(k): time t k The actual RUL at the location;

[0215] Time t k The estimated RUL at

[0216] -π: RUL distribution;

[0217] -β e : the minimum acceptable confidence (probability, likelihood, etc.) about the center of gravity of the distribution π;

[0218] The actual RUL value at time t k à±α e The distribution of .

[0219] The goal here is simply to provide trends in advance so that maintenance actions can be predicted.

[0220] Therefore, the constant α e It will be defined as a ±15% interval around the initial RUL.

[0221] The relative accuracy is achieved without any additional information. In contrast, the prognostic range and performance require further information.

[0222] Since RUL is not given as a distribution here, it is not necessary to define a confidence level β e Therefore, the prognostic range will be calculated as follows:

[0223] PH=t EOL -t αe ,in:

[0224] t αe is the RUL estimated inspection time until end of life:

[0225]

[0226] The best range is obtained when the algorithm is always in the required performance range, and the worst range is obtained when the algorithm is never in the required performance range. In this case, α e is the uncertainty interval of the initial RUL, RUL(0).

[0227] Figure 19 and Figure 20 An example illustrating the performance of the prognostic range. Figure 19 According to the criterion α e Two algorithms are compared. Figure 20 According to the criterion α e and β e Evaluation of algorithms from distributions.

[0228] RUL (Remaining Useful Life) corresponds to the remaining lifetime from the current time t. The classification model then calculates the membership u of the severity class j (j ranges from 1 to c.) From this, the RUL is derived as follows.

[0229] The estimated RUL is equal to the difference between the current (predicted) time and the estimated end-of-life time q. Therefore, for 1≤q≤hp:

[0230]

[0231] wherein q is the estimated end-of-life time and t is the current time, and wherein at the estimated end-of-life time q, the last severity category membership function becomes greater than the penultimate severity category membership function.

[0232] Finally, the membership of class c is higher than that of class c-1, and therefore:

[0233]

[0234] As an example, Figure 21 The graph in gives the membership estimate for bearing 1d based on the prediction at t=10.45 days (learning was performed on bearings 1a, 1b and 31c).

[0235] Curve 70 corresponds to the u1 class 1 membership function, curve 71 corresponds to the u2 class 2 membership function, curve 72 corresponds to the u3 class 3 membership function, and curve 73 corresponds to the u4 class 4 membership function.

[0236] The limit value 74 corresponds to the beginning of category 4 (ie the last category) and therefore to the end of life time.

[0237] Therefore, the end of life is estimated to be:

[0238] q = 27.23 days, and the remaining lifespan is:

[0239] =qt=27.23-10.45=16.78 days.

[0240] In addition to the prediction of t=sp2=10.45 days, the remaining service life of the 1d roller was estimated by two other prediction models t=sp1=7.5 days and t=sp3=20.65 days.

[0241] Figure 22 Convergence to the correct EOL value is shown approaching the EOL within a starting tolerance of + / - 15%, which converges to 0% at the actual EOL.

[0242] Curve 75 is the actual value of the remaining service life. Curve 76 corresponds to the lower tolerance, while curve 77 corresponds to the upper tolerance. Curve 78 corresponds to the prediction t=sp1, curve 79 corresponds to the prediction t=sp2, and curve 80 corresponds to the prediction t=sp3.

[0243] It can be seen that all three estimates are very close to the actual RUL, which shows that the estimation process according to the present invention is very accurate.

[0244] Of course, the invention is not limited to the described embodiments but covers any alternative falling within the ambit of the invention such as defined by the claims.

Claims

1. A method for estimating the remaining useful life of a test equipment, comprising a preliminary stage and a working stage, wherein the preliminary stage comprises the following steps: Collecting test observations of test equipment similar to the equipment under test; generating at least one signed test time sequence based on the test observations; Dividing the test time series to obtain severity categories corresponding to aging stages of the test equipment devices; performing a first learning of a diagnostic model on the test equipment device; performing a second learning of a prediction model for the signature; And the working phase includes the following steps: collecting the observations while operating on the equipment under test; generating, using the predictive model, an extrapolated time series representing the evolution of the signature on the test equipment; classifying the extrapolated time series using the diagnostic model; calculating a membership function of the extrapolated time series to the severity category; According to the membership function, the remaining service life of the test equipment is derived, The estimated remaining useful life is such that: wherein q is the estimated end-of-life time and t is the current time, and wherein at the estimated end-of-life time q, the last severity category membership function becomes greater than the penultimate severity category membership function.

2. The method according to claim 1, characterized in that The temporal consistency criterion is used to define the several severity categories generated by the partitioning.

3. The method according to claim 2, characterized in that The temporal consistency criterion is evaluated using a first metric defined as follows: where x k (i) is the i-th test observation collected for the test equipment k, m u is the number of observations, and where ct is the temporal consistency between two consecutive test observations and is defined using the following formula: where y k (i) is the observation x k (i) severity category.

4. The method according to any one of the preceding claims, characterized in that A second metric, which is an energy metric, is evaluated to select a partitioning algorithm.

5. The method according to claim 4, characterized in that The energy metric is defined as: where Ω i is the test observation set of category i, is the test observation, g i is the centroid of the i-th severity category, d is a distance measure and is defined by: Among them A d is a positive semidefinite matrix defining a distance metric, and where the symbol "|||| d ” denotes the norm corresponding to the distance d.

6. The method according to claim 4, characterized in that The partitioning algorithm is the C-fuzzy mean algorithm.

7. The method according to claim 1, characterized in that The prediction model is trained using a multi-step iterative method.

8. The method according to claim 1, characterized in that The prediction model was constructed based on nonlinear SVR regression.

9. The method according to claim 1, characterized in that The first learning of the diagnostic model uses the fuzzy SVM-FA classification algorithm.

Images