Estimation of machinery's remaining useful life in the presence of non-Gaussian noise by using a robust extended Kalman filte

Estimation of machinery’s remaining useful life in the presence of non-Gaussian noise by using a robust extended Kalman filter Hamid Shiri ∗a, Pawel Zimroza, Agnieszka Wyłomańskab, Radosław Zimroza,a aFaculty of Geoengineering, Mining and Geology, Wroclaw University of Science and Technology, Na Grobli 15, 50-421 Wroclaw, Poland bFaculty of Pure and Applied Mathematics, Hugo Steinhaus Center, Wroclaw University of Science and Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland Abstract Estimation of the remaining useful life (RUL) of industrial machinery is essential for condition-based maintenance (CBM). While numerous papers have explored this issues, challenges arise as machinery often works in non-stationary conditions, particularly in harsh environments (like mining machines, wind turbines, helicopters, etc). The data collected from such environments are affected by non-Gaussian noise, posing difficulties for traditional approaches to non-linear state estimation or prediction. The widely used extended Kalman filter (EKF) suffers from the non-Gaussian noise effect due to its recursive minimum L2-norm filtering. To address these issues, we propose a robust EKF based on the maximum correntropy criterion. This method effectively estimates the RUL of the time-varying degradation process in the presence of nonGaussian noise, also enabling confidence interval computation for uncertainty management. The efficiency of our approach was confirmed through application to simulated and benchmark data sets, outperforming Kalman filter-based methods for both simulated and real-world scenarios. Keywords: prognostics, remaining useful life (RUL), non-linear degradation, extended Kalman filter, robust methods, non-Gaussian noise. Nomenclature CBM Condition-based maintenance CI confidence interval CKF Cubature Kalman filter CM Condition monitoring EKF Extended Kalman Filter EOL End of the life FPT First prediction time HI Health index IMS Intelligent maintenance system data set KF Kalman Filter MAE Mean absolute erro MCEKF Maximum correntropy-based Extended Kalman Filter MCSKF maximum correntropy criterion switching Kalman filter ML Machine learning MSE Mean square error PDF Probability density function PF Particle filter PHM Prognostics and health management ∗Corresponding author. ∗∗ Corresponding author, [email protected] Email address: [email protected] (Hamid Shiri ∗) Preprint submitted to Measurement November 12, 2025 QKF Quadratic Kalman filter RMS Root mean square RUL Remaining useful life SMC Sequential Monte Carlo TRB tapered roller bearing UKF Unscented Kalman filter 1. Introduction Condition-based maintenance (CBM) revolutionises maintenance strategies by continuous monitoring of equipment and systems in real time, allowing precise maintenance decisions based on their current condition instead of rigid schedules. By leveraging advanced sensors and data analysis, the CBM rapidly identifies deviations from normal operation, ensuring optimal asset health assessment. This enables organisations to minimise downtime, extend the life of equipment, and reduce maintenance costs, leading to unparalleled operational efficiency in diverse industries. In the context of CBM, the role of RUL or prognostics takes centre stage, providing critical insights into an asset’s operational longevity based on its current state and performance trends. RUL serves as an early warning system for potential failures, empowering proactive and cost-effective maintenance planning that further reduces downtime and maximises asset lifespan. This data-driven methodology equips organizations to make informed maintenance choices, strategically prioritize vital components, and refine maintenance resource allocation, ultimately bolstering operational efficiency and curbing maintenance expenditures. Prognostic approaches can be categorised into three main groups: the data-driven approach (including the machine learning (ML) based approach [1,2], and the statistical model-based approach [3–5] ), physics-based approach [6,7], and hybrid approach [8,9]. Physics-based model approaches try to explain the degradation process by taking advantage of the physics of this process based on damage and fracture mechanics [10]. This family of approaches could provide accurate results. However, the exact physics knowledge of complex systems is not always available or too expensive to extract [11]. Therefore, in real applications, the use of physics-based model approaches is subject to significant restrictions. The data-driven approach aims to build a model of the degradation process by using historical data, particularly valuable when examining intricate engineering systems such as wind turbines, aircraft, and mining machines. Within this approach, there are two prominent subclasses: the machine learning-based approach and the statistical model-based approach, each offering distinct advantages and prerequisites. The machine learning-based approach serves as a powerful tool for modelling, segmenting, and predicting data from complex systems with degradation processes that defy easy correlation with physics or statistics. Nevertheless, these methods demand a substantial volume of historical degradation data for effective training, a resource that often proves elusive for many practical equipment scenarios. In contrast, statistical model-based approaches require comparatively less degradation data and do not require advanced knowledge of equipment mechanics. Moreover, they exhibit considerable potential to accommodate degradation uncertainties. When opting for a statistical model-based approach, the selection of an appropriate model that accurately characterises the degradation process becomes crucial. In summary, while both machine learning and statistical model-based approaches are integral to data-driven methodologies, the former thrives in complex, data-rich environments but demands more historical data, whereas the latter excels in scenarios with limited data and uncertainty considerations, prioritising the choice of an effective degradation model. The hybrid approach aims to benefit from the advantages of both physics-based and data-driven approaches by combining them. More detailed information on the hybrid approach can be found in the literature [8,9]. Kalman filter (KF) is one of the well-known methods to perform state estimation for engineering purposes, which has been used regularly for diagnostics and prognostics [12–15]. However, most of the real cases, especially in the area of prognostics, are non-linear, which creates a limitation for using KF. To solve the system’s non-linear state estimation problem, many extensions of KF, such as EKF [16], quadratic KF (QKF) [17], unscented KF (UKF) [18], cubature KF (CKF) [19], have been proposed in the literature. EKF is computationally efficient and simpler to implement rather than other mentioned filters, however, might not accurately handle with highly non-linear systems. The fundamental idea of EKF is to linearise the system using the first-order term of the Taylor series expansion of non-linear functions. 2 KF family usually uses second-order statistics and can perform optimal solutions when the noise is Gaussian distributed. However, most of machines in industries work in harsh environments [20]. Due to this, the data collected from these machines are usually affected by noise with non-Gaussian characteristics, which can pose a great challenge for the classic KF family approach. We should clarify that in this paper we consider heavy-tailed noise, which is a specific type of non-Gaussian noise characterised by the presence of extreme or outlying values in a data set, often leading to a distribution with thicker tails than the Gaussian (normal) distribution. Non-Gaussian noise, unlike the more typical Gaussian noise, is frequently encountered in practical applications, such as monitoring the health of mining equipment and wind turbines. For instance, when assessing the condition of mining machinery dealing with the screening of falling ore, the noise in the data often exhibits non-Gaussian characteristics due to irregular patterns and unexpected behaviors during the sieving process. Similarly, in the context of wind turbines, turbulence in the surrounding wind flow can introduce non-Gaussian noise into the collected data. Additionally, in scenarios involving electrical devices or electromagnetic fields, non-Gaussian noise may arise from the unpredictable behavior of electromagnetic interference or fluctuations in power supply. By acknowledging the existence of non-Gaussian noise originating from sources like falling ore in sieving screens, wind turbulence, and electromagnetic fields, we can better adapt our diagnostic and prognostic methods to accommodate these unique statistical properties. Recently, many methods have been developed to handle the effect of non-Gaussian noise in classic KF. A part of these approaches is developed based on considering the non-Gaussian distribution for noise, which is fulfilled for distributions like alpha-stable distributions, Student’s t distribution, etc [21,22]. Nevertheless, finding the analytic solutions for these distributions is a hard task and a high computational cost is also required. Another approach uses the Gaussian mixture distribution to approximate non-Gaussian noise [23–26]. The high computational cost can be mentioned as the main weakness of this family of methods. The third family of approaches that attempts to address a non-Gaussian noise problem is sequential Monte Carlo (SMC) sampling, which can approximate any distribution [27]. The particle filter (PF) can be the best example of this family of methods [24,28]. As with previous methods, the high computational cost is the most significant disadvantage of this family of approaches. The computational cost is a key parameter in the selection of methods in the subject of prognostics and RUL estimation, mainly because it is needed to process a large amount of data. To treat the high computational cost and the non-Gaussian noise effect, some papers have proposed a replacement for the mean-square error (MSE) criterion (inherent to the standard KF). For example, based on the definition of corr entropy which is a similarity measure (for introduction to topic, see [29]), some research used the correntropy criterion to drive KF [12,30–33]. Since the correntropy criterion uses a higher-order statistic of the noise instead of second-order statistics, the performance of state estimation will increase in comparison to standard KF in the presence of noise with non-Gaussian characteristics. In [12,34], Shiri et al. introduced a novel approach that utilises the maximum correntropy switching Kalman filter (MCSKF) for state detection when encountering non-Gaussian noise in a system. Their findings highlight the promise of MCSKF for online state detection, although their research primarily focusses on machine state detection without addressing predictive applications. It is important to note that the methods proposed in their work rely on a linear KF framework, which may not be universally applicable due to its linear assumptions in many practical scenarios. In this study, we aim to enhance the accuracy of RUL predictions and offer probabilistic outcomes to tackle the challenges posed by dynamic, non-linear degradation processes combined with noise exhibiting non-Gaussian characteristics. To achieve this, we introduce a robust EKF that is based on the maximum correntropy criteria, i.e. the maximum correntropy EKF (MCEKF). This approach seeks to address the complexities of the RUL estimation problem by considering both the uncertainties in the degradation process and the non-Gaussian characteristics of the noise. The main contribution of this paper can be summarized as follows: •a robust version of EKF is driven based on maximum correntropy criteria to deal with non-linear degradation processes combined with noise with non-Gaussian characteristics; •comparative simulations with the other filters are executed by Monte Carlo analysis in the presence of Gaussian and non-Gaussian noise to indicate the dominance of MCEKF; 3 •the proposed MCEKF is carried out on three benchmark real data sets to estimate RUL and verify the efficacy of the approach in real environments. The rest of the paper is organised as follows. In Section 2, the lifetime curve models are examined based on the literature. In Section 3, the theory for the critical parts of the processing methods is described. Then, in Section 4, the results for the simulated data are shown. Ultimately, the results for three benchmark data sets obtained using the proposed approach are demonstrated in Section 5with an expression of all intermediate steps. Conclusions are formed in Section 6. 2. Lifetime curve model In the PHM community, there are few statistical models that are used to describe the degradation process. These models are usually made up of the deterministic part (which is trying to qualify the global trend of the degradation process) and the random part (which is trying to consider the uncertainty of the degradation process). They are selected based on the behaviour of the degradation process with different deterministic and random trends [20]. In this work, we selected a time-varying model with three different regimes. The first regime has a constant deterministic trend, which refers to the healthy state of the system (the machine is working stable). The second regime has a linear growing trend which is assigned to the degradation of the machine. The last regime has an exponential trend, which is attributed to the critical state of the degradation machine. Moreover, the FPT starts from the last regime in this model. Also, as discussed in the literature [20], modelling of the random part is very important to reflect the uncertainty. In some models it is assumed that the variance of the noise is constant and the distribution of the noise is Gaussian [35–37], see Fig. 1. However, Zulawinski et al. [20] demonstrated that the random part may be non-homogeneous (has a time-varying scale), and the noise distribution is far from the Gaussian distribution. Therefore, we used a model with non-stationary characteristics and non-Gaussian noise; see Fig. 2. The scale (variance) of the random part changes during the degradation process based on the regime in which we are in. Figure 1: The long-term degradation model with three stages. Figure 2: The proposed long-term degradation model with three stages. It is important to note that our study operates under the assumption of a known FPT or first point of the critical stage, from which our prediction methodology commences. Due to this fact, there is no valuable information on the healthy stage for prediction, and starting the prediction from the degradation stage will lead to overor underestimation of RUL. It is worth highlighting that the accurate identification of the FPT is a crucial issue in the context of the estimation of the RUL, bearing a significant effect on the eventual results. However, we must clarify that the exploration of FPT detection lies beyond the scope of this current work, which is centred on the prediction task. Consequently, we leverage the findings derived from our previous research endeavours, where the dedicated focus was directed toward FPT detection within real-world data sets. For those interested in an in-depth analysis of this facet, we kindly direct them to our previous comprehensive publications that explore extensively this subject [2,12,38,39]. In the present study, we build upon the results of our prior work [12], which introduced a robust methodology 4 based on the switching maximum correntropy Kalman filter (SMCKF) to address the challenges of threshold problem (selecting the threshold for warning and alarm not for EOL) and online diagnostics, especially in the presence of non-Gaussian noise, utilising condition monitoring (CM) data. This approach leverages a set of dynamic system models to elucidate various stages of degradation, employing robust Bayesian estimation techniques. Due to its foundation in dynamic behaviour modelling, this approach obviates the need for a fixed threshold in the diagnostic process. Therefore, we can characterise our current paper as an extension of our previous work, with a specific focus on prediction tasks and addressing nonlinearity. 3. Methodology This section outlines the methodology grounded in a degradation process for the estimation of RUL. The RUL estimation workflow, based on our proposed approach, is illustrated in Fig. 3. The data obtained, such as sensor-recorded vibration data, serves as a means to monitor the machine’s health condition. Through signal processing techniques, HI is constructed from the measured signals, effectively representing the machine’s health status. Depending on the evolving degradation patterns of HI, the life of the machinery is categorised into three distinct stages: health, degradation, and the critical stage, where the FPT is identified as the initial entry into the critical stage. In the critical stage, we defined a state-space degradation model, and its parameters are estimated using the MCEKF and updated in real time as new HI values become available. An EOL threshold is defined, enabling the calculation of the PDF of the estimated RUL based on the constructed degradation model and the EOL threshold. Furthermore, the details of the state space degradation model and the theoretical basis of this process, which includes the derivation of the EKF and MCEKF, are described in the following subsections. Start Raw vibration data. Use features as health index. Detect FPT. Estimate parameters of the state space degradation model by using MCEKF and simulate trajectories by using estimated parameters. Set EOL threshold. Calculate PDF of RUL. End Figure 3: Complete procedure for derive distribution of RUL 3.1. RUL RUL is a critical metric in the field of machinery prognostics, representing the duration for which a machine is expected to operate efficiently before the need for repair or replacement. Fig. 4illustrates the concept of RUL. Our methodology focusses on training a predictive model using historical data exclusively. Subsequently, we employed this model to forecast future values of the degradation process, accompanied 5 by CI. The estimated degradation curve is then examined to identify the point in time at which it initially crosses a predetermined EOL threshold, denoted tEOL. The predicted RUL is calculated as the difference between this failure time tEOL and the time of the last observation in the training data set (t0), such that the predicted RUL =tEOL −t0. Figure 4: Complete procedure for derive distribution of RUL. 3.2. State space degradation model As confirmed by the findings in [20], the exponential model consistently demonstrates a better global regression performance compared to polynomial and other conventional models. Consequently, we have opted to adopt the ensuing exponential model as the degradation pattern designated for our analysis of the third regime of the degradation process (critical stage): yt=aexp(bt) + c, (1) where observation ytcorresponds to the HI, tis the discrete time and a, b, c represent the parameters of the model that need to be initialized with real data. Recognising the intricate and time sensitive nature of the degradation process, this study postulates that the model parameters at, bt, ct,are subject to temporal variations. This assumption aligns closely with the dynamics of the degradation process in the real world. Let the dynamic state xtbe defined as a column vector xt= [at, bt, ct]T. Furthermore, recognising the impact of external noise and inherent uncertainty, we introduce the subsequent state space model: (xt=xt−1+qt yt=atexp(btt) + ct+mt ,(2) where, the process noise is denoted as qt, presumed to be uncorrelated zero-mean white noise. This noise factor is characterised by the covariance matrix Qt. Besides, we assume mtto be a scalar noise, not necessarily Gaussian, to suit heavy tailed noise case inherent in the real data (her it has no longer the interpretation of the measurement noise). The noise mtis characterized with a time changing scale parameter Mt. Keeping an analogy with the degradation model from Eq. (5), we map the relationship between HI and the parameters featured in Eq. (2) denoting it with the non-linear functions fand has follows: f(xt−1) = xt−1,(3) h(xt) = atexp(btt) + ct.(4) 6 Furthermore, we set the initial state x0and the initial values for Q0and M0as well as initial value P0|0 for the covariance matrix Pt|twhich is described in the next subsection. In the following two subsections we describe EKF and introduce MCEKF considering more general notation. Please note, that in general approach ytwould be replaced with a vector yt, the scalar function hwith a vector function hand the scale parameter Mtwould be replaced with the matrix Mtwhich refers to the covariance matrix in the Gaussian case. 3.3. Extended Kalman filter The EKF works based on the use of Taylor series to transform non-linear filtering problems into linear forms. The discrete non-linear state-space representation of the model can be presented as follows. xt=f(xt−1) + qt, yt=h(xt) + mt,(5) where xtis the unknown dynamic state at time tand ytis its observation, respectively – both of them are in general vectors of selected size. Furthermore, qtand mtrepresent processes and measurement noise that are considered Gaussian distributed: qt∼ N(0,Qt),mt∼ N(0,Mt). Furthermore, fand hare non-linear state transition function and measurement function, respectively. Let us denote by Pt|t−1the prior estimate covariance matrix: Pt|t−1=E[(xt−ˆxt|t−1)(xt−ˆxt|t−1)T](6) and with Pt|tthe posterior estimate covariance matrix: Pt|t=E[(xt−ˆxt|t)(xt−ˆxt|t)T].(7) The hat operator here means the estimation procedure made by KF conditioned on the specific time (normally on the current step tor the previous step t−1). The EKF can be summarised as follows with the two steps iterated subsequently, namely the prediction step: ˆxt|t−1=f(ˆxt−1|t−1), ˆ Pt|t−1=Atˆ Pt−1|t−1AT t+Qt,(8) and the update step: ˆxt|t=ˆxt|t−1+Kt(yt−h(ˆxt|t−1)), ˆ Pt|t= (I−KtHt)ˆ Pt|t−1(I−KtHt)T+KtMtKT t.(9) where Tmeans matrix transposition, Iis identity matrix of the proper size. The term Ktin the above equations is called Kalman gain and for the standard EKF its optimal value is given by the formula: Kt=ˆ Pt|t−1HT t(Htˆ Pt|t−1HT t+Mt)−1,(10) where Atand Htare extracted by calculating the Jacobian matrices fand h, respectively, as follows: At=∂f ∂x|ˆxt−1(11) Ht=∂h ∂x|xt(12) 3.4. Maximum correntropy Extended Kalman filter In the following, for the model expressed with Eq. (5) another form of an objective function is introduced. Basis for this modification is the maximum correntropy criterion. In this work, we use Gaussian kernel function, and thus the correntropy of two random vectors Xand Yis defined as following [29]: C(X,Y) = E[G(∥X−Y∥)].(13) 7 where Gσ(u) = exp(−u2 2σ2)and σis the size of the kernel. The maximum correntropy criterion states that the distance (measured with Euclidean norm) between two random vectors is minimised when the correntropy, defined as the expected value of an arbitrary kernel function, is maximised. The benefit of using correntropy is that, unlike the squared error cost function, in its Maclaurin series expansion it has terms of higher order than second and is insensitive to outliers (as was described in [31]). Based on this criterion, there are different possible formulations of the cost function for MCEKF as a discrete sum (sample mean) of kernel function values. In [40] it is summed through all eigenvalues of both Mtand Pt|t−1. We, on the other hand, define it as a sum of two values of the kernel function, similar to [31] and our previous work [12]: JMC =1 2Gσ(∥yt−h(xt)∥M−1 t) + 1 2Gσ( xt−f(ˆxt−1|t−1) P−1 t|t−1 ).(14) Note that it takes the two standardised noise terms of the same model that EKF does in the previous section, but it combines them in a different way. Using this cost function, compared to using a quadratic cost function from the standard KF, makes the estimation more robust in the presence of non-Gaussian noise, which is the case when the real observation process does not necessarily fulfil the assumptions of the KF model (see [31]). The formulation of MCEKF can be derived by finding xtthat maximizes the objective function JMC given by Eq. (14). The resulting prediction and update steps are driven by the same equations as for the standard EKF; see Eq. (8), (9) and (11), (12) but the optimal Kalman gain for the MCEKF is given with the formula: Kt=ˆ Pt|t−1λtHT t(Htˆ Pt|t−1λtHT t+Mt)−1,(15) where λt=Gσ yt−h(ˆxt|t−1) M−1 t.(16) 3.5. Discussion of the hyperparameters setting In the state estimation procedure using MCEKF, when a significant outlier appears, the term yt− h(ˆxt|t−1)in Eq. (9) (usually called innovation or prefit residual) diverges, but Ktcontrols the divergence of the estimator ˆxt|t. See [30–33] for more details on the procedure of driving equations and stability. Also, we should clarify that, in the MCEKF, the σparameter (kernel size) controls how sensitive the filter is to differences between predicted and measured states. A larger σmakes the filter more tolerant to deviations, while a smaller σmakes it more sensitive to small differences, which could help to reject outliers but might also be more affected by noise. Therefore, selecting a proper value of σis crucial for balancing noise robustness and sensitivity to data deviations. It often requires some experimentation or domain knowledge to determine the best σvalue for a given application. In addition, the measurement error, denoted as mt, is determined by calculating the variance of stationary measurements acquired when the machinery is in a healthy operational state. It is important to note that the specific constant value of Mt=Mmay vary among individual cases. The assumption of stationary measurements from well-functioning machinery, as applied in this study, is rigorously validated by analysing data from multiple cases, thus confirming its practical applicability. As discussed in [41], the process error, represented as qt, encapsulates the uncertainty associated with the filter’s ability to model real-world dynamics. A lower process error enhances the filter’s fitting of condition monitoring (CM) data, making the model more responsive to delicate changes, although potentially making it susceptible to overfitting. In our approach, the process error, denoted as qt=q, is fine-tuned for the MCEKF model based on historical data from similar defect cases and is assumed to be consistent across cases. In this study, an arbitrary covariance matrix is employed, denoted P0|0. With the establishment of these matrices and the definition of the initial parameters, the MCEKF model for the prognostic evaluation of bearing degradation is formulated. 8 3.6. EOL The EOL in machinery prognostics is a critical parameter, and selecting the right threshold is paramount for effective maintenance and operation. It marks the point when a machine is expected to become unreliable, which makes it crucial for scheduling preventive maintenance, cost reduction, safety, resource optimisation, reliability, and environmental impact. The key to success lies in data-driven approaches, such as predictive maintenance techniques and thorough data analysis, which enable organisations to accurately identify and set appropriate EOL thresholds for their machinery. It is worth noting that choosing an appropriate EOL value for real data sets presents a considerable challenge, and this has been the subject of extensive research over the past decade [37,42]. The complexity of this task endures. In our study, we adopted a random selection approach for the EOL value. This choice comes from the flexibility of our proposed method, which is not restricted to predefined EOL values. It can be easily computed for a range of EOL values, offering versatility in its application. Additionally, in Appendix Appendix B, pseudocodes detailing our proposed methodology are provided for reference. 4. Simulation and results 4.1. Model of the degradation curve In the following section according to [20], the HI data is generated. The simulation is restricted only to the third regime starting from the FPT as was argued in Section 2. As was discussed there, the HI is constructed from two parts: HI(t) = D(t) + R(t),(17) where D(t)and R(t)are as deterministic and random parts of the HI, respectively. The deterministic part in Eq. (17) is described as follows: D(t) = a1exp(b1t) + c1,(18) where a1,b1, and c1are constant values that are used to build the exponential function that corresponds to the critical stage. In the simulation part, we consider Gaussian and Student’s t distributions for the random component of the degradation process (for more details about Student’s t distribution, see Appendix A). The time-varying random component corresponding to R(t)is generated in the following way: R(t) = SC(t)˜ R(t),(19) where SC(t)corresponds to the scale of the random part (variance) which is constructed as follows: SC(t) = a2exp(b2t),(20) where values a2,b2are constant. The ˜ R(t)function represents a series of independent identically distributed random variables (iid). For the Gaussian distribution case, we consider ˜ R(t)∼ N(0,1), and for the Student’s t case, we have ˜ R(t)∼t(ν) where νis the parameter called the number of degrees of freedom. 4.2. Results for simulated signal In this subsection, we apply our proposed methodology to simulated data based on Subsection 4.1. To show the performance of the proposed method and make a comparison with other filters from the KF family, we applied the classic EKF and UKF as well-known modifications of the classical KF for non-linear state estimation. For more details on EKF and UKF we encourage the reader to read [43,44]. 9 sets. The first and second data sets originate from the FEMTO and IMS benchmarks, which are renowned for their extensive use as a benchmark in various studies. In particular, these data sets present evidence of a Gaussian distribution. The last data set under scrutiny pertains to wind turbine bearing data and notably encompasses solid components that do not stick to the Gaussian distribution. The existence of such non-Gaussian behaviour introduces an additional layer of complexity to our investigation, warranting an in-depth analysis of the model performance in such scenarios. Also, note that for each data set, all initial parameters of the filters were selected as the same values. 5.1. Data sets description 5.1.1. FEMTO data set This data set was made publicly available during the IEEE International conference on PHM in 2012 and was generously provided by the Franche-Comté Electronics Mechanics Thermal Science and Optics–Sciences and Technologies Institute [45]. Comprising 17 run-to-failure data sets of rolling element bearings acquired from the PRONOSTIA platform, see Fig. 18, this data set was utilised to perform accelerated degradation tests, effectively simulating several years’ worth of bearing wear in a few hours. This was achieved by applying a high-level radial force that exceeded the maximum dynamic load of the bearings. Throughout the tests, the rotational speed of the bearings remained constant. Data acquisition was facilitated using two accelerometers and a thermocouple, capturing vibration signals and bearing temperatures. Bearing failure was deemed to have occurred when the amplitude of the vibration signal exceeded 20 g. Notably, this data set represents natural-occurring bearing degradation, devoid of any pre-seeded faults or artificial interventions. Compounding the complexity of the data set is the limited availability of only two training units for each operational condition. Furthermore, diverse fault patterns and lifetimes characterise different units, even when subjected to the same conditions. This heterogeneity between the training and testing units exacerbates the challenge of predicting RUL. Vibration signals within the data set exhibit a relatively low frequency resolution, as each sample spans a time length of 0.1 s, resulting in a frequency resolution of 10 Hz. Consequently, conventional fault diagnostic methods that rely on detailed frequency analysis are not suitable for this data set. Furthermore, the data set exemplifies the propensity of incipient faults in one component to propagate to other components through frequent contact, resulting in the simultaneous occurrence of various fault patterns, as visually represented in Fig. 19. This distinctive behaviour underscores the suitability of the data set for addressing multi-fault prognostic challenges. This data set is commonly used for health index construction [46–53], health stage evaluation [2,54–58], and predicting RUL [59–64]. Figure 18: FEMTO test rigs [45]. Figure 19: Different faults of rolling element bearings [45]. FEMTO case study In this paper, we apply our proposed methodology to estimate the RUL of the bearing labeled as 1−1 within the FEMTO data set. The data in this data set have been gathered under specific operating conditions, precisely at 1800 rpm and with an applied load of 4000 N, as shown in Fig. 20. 16 For each set of vibration data, we calculate the RMS, which is subsequently used as the health index for our analysis. Our choice of this specific case study was predicated on the compelling characteristics that it exhibits, rendering it an exemplary illustration of a three-stage degradation model, marked by predominantly monotonic behaviour. Furthermore, as indicated in [65], the random component of this case study closely approximates a Gaussian distribution. This aspect lends substantial support to our objective of demonstrating the efficacy of our proposed model in a real world scenario where the underlying distribution aligns with a Gaussian distribution. Figure 20: Raw vibration data for case study bearing 1−1. 5.1.2. IMS data set The data set under consideration was provided by the centre for intelligent maintenance systems (IMS) at the University of Cincinnati [66], and it has been made available through the NASA prognostic data repository of NASA [67]. This data set comprises three distinct subsets of bearing degradation tests. Each test involved the installation of four Rexnord ZA-2115 double-row bearings on a common shaft, see Fig. 21. Accelerometers were placed on bearing housings to capture vibration signals and an oil circulation system was engineered to ensure bearing lubrication. For debris collection, a magnetic plug was integrated into the oil feedback pipe. The execution of the test was designed to be adaptive, with an electrical switch that triggers cessation once the accumulation of debris exceeded a predetermined threshold. Following the conclusion of these tests, the bearings underwent a detailed examination and the observed fault patterns were meticulously recorded. Vibration signals exhibit high frequency resolution, making them well suited for fault diagnosis through frequency analysis techniques. According to the data instructions [67], each file consists of 20,480 data samples sampled at 20 kHz. Consequently, operators can derive valuable frequency-domain features to monitor the degradation progression of specific bearing components, such as rollers, outer race, and inner race. At the conclusion of the tests, detailed fault patterns for each bearing were observed and recorded, as depicted in Fig. 22 offering a rich source of information for researchers to explore the relationships between different fault patterns and their corresponding degradation trends. Moreover, note that this data set has been extensively used in various publications such as segmentation [68,69], RUL prediction [70–74], and condition monitoring [66]. 17 Figure 21: IMS test rig [66]. Figure 22: Different faults of rolling element bearings [66]. IMS case study In this part, we employ the IMS data set, which comprises three distinct subsets. For our study, we specifically focus on Subset 3 and select the vibration of bearing number 3 as our case study, see Fig. 23. Subset 3 captures data over a recording duration spanning from April 4, 2004, at 09:27:46 to April 4, 2004, at 19:01:57. This data set encompasses vibration data, which were collected from four channels. The channel arrangement is as follows: bearing1-channel1; bearing2-channel2; bearing3-channel3; bearing4-channel4. Data are recorded at a 10-minute interval in ASCII file format. The experiment concludes with an outer race failure event on bearing number 3, see Fig. 22. Our choice of this particular case within the IMS data set was underpinned by its distinctive features, most notably the remarkably brief transition between the degradation and critical stages, which might even warrant classification as a two-stage model rather than a multistage model. This specific case also provides an additional scenario where the noise distribution closely approximates a Gaussian distribution, further enhancing its relevance for our study. Figure 23: Raw vibration data for IMS subset 3 bearing 3. 18 5.1.3. Wind turbine data set The wind turbine data set comprises sensor data collected from the high-speed bearing shaft of a 2.2 MW wind turbine, visually depicted in Fig. 24. This data set encompasses measurements of the bearing’s inner race energy, sampled at 10-minute intervals over a period of approximately 50 days. In particular, the bearing under investigation shows an inner race fault and is identified as SKF 32222 J2 tapered roller bearings. This particular tapered roller bearing (TRB) boasts a 200 mm outer diameter, a bore of 110 mm, and a total length of 56 mm, featuring 20 rolling elements set at a 16 °taper angle and weighing approximately 20 pounds. The bearing in question is supported by two pillow blocks and is equipped with a load cell to measure the force applied to the bearing. The wind turbine operates at variable speeds ranging from 2 to 100 Hz, with a load cell capable of measuring loads up to 1000 pounds, although most tests were conducted at 150 or 300 pounds. To mitigate the risk of catastrophic gearbox failure, the highest applied test load was maintained at 50% of the rated power. Further comprehensive information on this data set can be found in [75]. Remarkably, this data set, which we refer to as the wind turbine data set, has received considerable attention in recent years as a topic of interest for prognostics of RUL in various research endeavours, exemplified by works such as [75–78]. Figure 24: Wind turbine test rigs [75]. Wind turbine data set case study For the application of our proposed approach to the wind turbine data set, the inner race energy is chosen judiciously as the representative health index. For more details on the procedure of extraction of health indexes, please see [79]. In particular, as shown in Fig. 37, this data set is characterised by a conspicuous presence of outliers. Consequently, it is prudent to categorise this data set as possessing a trend punctuated by non-Gaussian noise. Furthermore, the health index exhibits fluctuations, which can be attributed to dynamic factors such as load variations or even phenomena such as self-healing. This complex behaviour, observed within real-world scenarios, presents a big challenge when it comes to the tasks of segmentation and RUL prediction. The presence of such intricate dynamics underscores the complexity inherent in accurately characterising and forecasting RUL in practical applications. 19 5.2. Results 5.2.1. Results for FEMTO data set In this part, the proposed method is used to estimate the RUL of the The FEMTO data set. Panel (a) in Fig. 25 shows the complete degradation curves for this case study. Also, the small window shows the area after FPT, which is used to estimate RUL. It should be noted that several papers focused on detecting FPT in the literature, including our previous papers [2,12,38,39], in which we used different methods. In this article, we utilise the results of [12] that are shown in panel (b) in Fig. 25. Therefore, here we focus on the estimation of the RUL that is starting from this point that is already identified. Figure 25: (a) Health index (RMS) of FEMTO data set bearing1−1(b) FPT detected that is provided in this [12]. Fig. 26, shows the results of applying the proposed approach to the FEMTO data set. Top panel illustrates the degradation curves and the level of EOL. The bottom panel shows the estimated results using EKF, UKF, and the proposed MCEKF. The black dashed line and the purple area show the real RUL value and the ±20% accuracy bound for predicting RUL. As expected, when the number of data is small, the accuracy of the estimated RUL by using all filters is not acceptable. Based on bottom panel, EKF, UKF, and MCEKF first overestimated RUL. After a while, around t= 30, it can be seen that all methods are reaching the accuracy bound and remain approximately inside until the end of life. In addition, the selected parameters for MCEKF are presented in Table 1and the parameters of EKF and UKF are selected as the same value as MCEKF. 20 Figure 26: Final results of estimated RUL for FEMTO data set, top panel: health index plus EOL, bottom panel: predicted RUL results performed by EKF, UKF, and proposed MCEKF and 20% accuracy bound. The PDF showing the estimated RUL values from the MCEKF technique is shown in Fig. 27. Evidently, in the initial stages, where the data set contains an insufficient number of data points, the PDF of the estimated RUL does not yield favourable results. However, a discernible improvement becomes apparent as the data set size increases, allowing the PDF to progressively enclose the RUL distribution. In particular, beyond the temporal threshold of t= 30, the mean of the PDF closely aligns with the actual value of the RUL. Furthermore, the 90% CI related to the estimated RUL, derived from 100 simulated trajectories using the MCEKF model, is presented in Fig. 28. As can be seen; after t= 30, the RUL values are within the scope of the 90% CI. Figure 27: PDF of estimated RUL by using MCEKF for FEMTO data set. Figure 28: The 90% CI of estimated RUL by using MCEKF for FEMTO data set. To illustrate the sensitivity of our proposed methodology to the hyperparameter σ(kernel bandwidth for the MCEKF), we conducted an analysis employing different ranges of σ. The results related to the estimated RUL are shown in Fig. 29. In particular, in this scenario, the selection of a higher value of σyields more accurate results compared to opting for a lower value σ. To reinforce this observation, we present the MAE computed between the real RUL and the estimated RUL using our proposed method in Fig. 30, further affirming the findings from the preceding Fig. 29. 21 Figure 29: Predicted RUL results performed by proposed MCEKF with different value σand 20% accuracy bound. Figure 30: MAE between the predicted RUL performed by proposed MCEKF with different value σand real RUL. Table 1: Selected parameters value for MCEKF in case of FEMTO data set. MCEKF parameters std of Qtstd of Mtσ FEMTO data set [0.00001,0.00001, 0.00001] 0.0038 0.87 5.2.2. Results for IMS data set In this part, the proposed method is used to estimate the RUL of IMS data set. Panel (a) in Fig. 31 shows the complete degradation curves for this case study. Also, the small window shows the area after FPT, which is used to estimate RUL. It should be noted that several papers focused on the detection of FPT in the literature, including our previous papers [2,12,38,39], in which we used different methods to identify this point. So in this paper we used the results of [12] that is shown in panel (b) in Fig. 31. Therefore, this paper only focused on the estimation of RUL from this point. Figure 31: (a) Health index (RMS) of IMS data set subset 3 bearing3(b) FPT detected by using a method that is introduced in this [12]. The results of applying all of the methods on the IMS data set are presented in Fig. 32. Top panel illustrates the degradation curves and EOL level. Bottom panel shows the estimated results using EKF, UKF, and proposed MCEKF. The black dash line and the purple area show the real RUL value ±20% 22 accuracy bound for predicting RUL. As expected, when the number of data is small, the accuracy of the estimated RUL by using all filters is not acceptable. Based on bottom panel, EKF, UKF, and MCEKF overestimated RUL at first. After a while, around t= 30, it can be seen that MCEKF is coming to the accuracy bound and remaining approximately in the accuracy bound until the end of the life, while we can see that UKF and EKF could not provide accurate results. Furthermore, we provide the details of the parameters chosen for the MCEKF in Table. 2. It is important to note that, for the fair comparative analysis, the parameters of the EKF and the UKF have been set to the same values as those used in the MCEKF. Figure 32: Final results of estimated RUL for IMS data set, top panel: health index plus EOL, bottom panel: predicted RUL results performed by EKF, UKF, and proposed MCEKF and 20% accuracy bound. The PDF of the estimated RUL values obtained using the MCEKF technique is shown in Fig. 33. Evidently, in the initial stages, where the data set contains an insufficient number of data points, the PDF of the estimated RUL does not yield favourable results. However, a discernible improvement becomes apparent as the data set size increases, allowing the PDF to progressively enclose the real RUL. In particular, beyond the temporal threshold of t = 30, the mean of the PDF aligns closely with the actual RUL value. Furthermore, the 90% CI related to the estimated RUL, derived from 100 simulated trajectories using the MCEKF model, is presented in Fig. 34. As it can be seen; after t= 30, the RUL values are in the scope of the 90% CI. Figure 33: PDF of estimated RUL by using MCEKF for IMS data set. Figure 34: The 90% CI of estimated RUL by using MCEKF for IMS data set. 23 To illustrate the sensitivity of our proposed methodology to the hyperparameter σ(kernel bandwidth for the MCEKF), we conducted an analysis employing different ranges of σ. The results pertaining to the estimated RUL are shown in Fig. 35. In particular, in this scenario, the selection of a σvalue between 0.25−0.4 yields more accurate results compared to opting for a lower or higher value σ. To reinforce this observation, we present the MAE calculated between the real RUL and the estimated RUL using our proposed method in Fig. 36, further affirming the findings of the preceding Fig. 35. Figure 35: Predicted RUL results performed by proposed MCEKF with different value σand 20% accuracy bound. Figure 36: MAE between the predicted RUL performed by proposed MCEKF with different value σand real RUL. Table 2: Selected parameters value for MCEKF in case of IMS data set. MCEKF parameters std of Qtstd of Mtσ IMS data set [0.00001,0.00001, 0.00001] 0.0038 0.27 5.2.3. Results for wind turbine data set In this part, the proposed method is used to estimate the RUL of the wind turbine data set. Panel (a) in Fig. 37 shows the complete degradation curves for this case study. Also, the small window shows the area after FPT, which is used to estimate RUL. It should be noted that several papers focused on detecting FPT in the literature for this data set, including our previous papers [2,12,38,39], in which we used different methods. In this article, we utilise the results of [12] that are shown in panel (b) in Fig. 37. Therefore, here we focus on the estimation of the RUL that is starting from this point 24 Figure 37: (a) Wind turbine health index (inner race energy), (b) FPT detected that is provided in this [12]. Fig. 38, shows the results of applying all the methods to the wind turbine data set. Top panel illustrates the degradation curves and the level of EOL. Bottom panel shows the estimated results using EKF, UKF, and the proposed MCEKF. The black dash line and the purple area show the real RUL value ±20% accuracy bound to predict RUL. As expected, for scenarios where the data set is limited, the RUL estimates by all filtering methodologies exhibit sub-optimal accuracy. Analysing bottom panel reveals initial results for both EKF, UKF, and MCEKF to overestimate RUL. In particular, a noticeable convergence towards the accuracy bond can be discerned around t= 180. However, it is imperative to acknowledge that due to the non-linear trajectory of the health index, achieving precise results remains a challenge. Subsequently, the results of all methodologies as time progresses are shown. As we can seen; after t= 300, the mean RUL estimated by MCEKF closely aligns with the accuracy bound, while the other results produced by EKF and UKF exhibit less accurate results. It should be noted that around t= 310 and t= 440, the RUL estimates derived from EKF and UKF appear to be susceptible to the influence of non-Gaussian noise. In contrast, the MCEKF, as anticipated, reduced the impact of such non-Gaussian noise. In this study, we have carefully selected the parameters for the MCEKF applicable to the case study at hand. These parameter choices are detailed in Table. 3. Importantly, we would like to highlight that, for consistent comparison, the measurement noise and process noise parameters have been set to the same values for the EKF, the UKF, and the MCEKF. 25 spectral kurtosis-derived indices and svr, Applied Acoustics 120 (2017) 1–8. [77] L. Saidi, E. Bechhoefer, J. B. Ali, M. Benbouzid, Wind turbine high-speed shaft bearing degradation analysis for run-tofailure testing using spectral kurtosis, in: 2015 16th International Conference on Sciences and Techniques of Automatic Control and Computer Engineering (STA), IEEE, 2015, pp. 267–272. [78] J. B. Ali, L. Saidi, S. Harrath, E. Bechhoefer, M. Benbouzid, Online automatic diagnosis of wind turbine bearings progressive degradations under real experimental conditions based on unsupervised machine learning, Applied Acoustics 132 (2018) 167–181. [79] E. Bechhoefer, R. Schlanbusch, Generalized prognostics algorithm using kalman smoother, IFAC-PapersOnLine 48 (21) (2015) 97–104. Appendix A. Student’s t distribution The PDF of the (unscaled) Student’s t distribution is the following: f(x) = Γ( ν+1 ν) ν 2 1 √νπ 1 (1+ x2 ν) ν+1 2 ,(A.1) where νis called the number degrees of freedom and Γ(·)is the gamma function. It should be noted that: the mean of the Student’s t distribution is defined if ν > 1as µ= 0; the variance of the Student’s t distribution is defined if ν > 0as σ=ν ν−2. Appendix B. Pseudocodes for proposed method Subsequently, to enhance clarity and facilitate understanding of the RUL estimation process by the MCEKF, we propose the following pseudocodes: Algorithm 1 RUL Estimation with MCEKF 1: Input: Observations ytfor discrete time t 2: Output: Predicted RUL, Probability 3: Initialize State and Covariance Matrix: 4: x0←Initial state vector [a0, b0, c0]T 5: P0←Initial covariance matrix 6: Initialize Process and Measurement Noise Covariance and maximum correntropy bandwidth: 7: Qt←Process noise covariance matrix 8: Mt←Time-changing scale parameter for measurement noise 9: σ←Bandwidth for MCEKF 10: for tin observations do 11: MCEKF Update: 12: xt,Pt←MCEKF-Update(xt−1,Pt−1,yt,Qt,Mt,σ) 13: Prediction of Next Observation: 14: yt←h(xt) 15: Prediction of RUL: 16: tEOL ←Time when degradation curve crosses EOL threshold 17: t0←Time of the last observation in the training data set 18: Calculate Probability of Predicted RUL: 19: P(Predicted RUL)←CalculateProbability(xt,Pt, tEOL, t0) 20: return Predicted RUL, P(Predicted RUL) 21: end for 32 Algorithm 2 MCEKF 1: Initialize state estimate ˆ x0and covariance matrix P0 2: for each time step tdo 3: Prediction: 4: ˆ xt|t−1=f(ˆ xt−1)▷Predicted state 5: At=∂f ∂x(ˆ xt−1)▷Jacobian of process function 6: Pt|t−1=AtPt−1AT t+Qt▷Predicted covariance 7: Update: 8: Ht=∂h ∂x(ˆ xt|t−1)▷Jacobian of measurement function 9: zt=yt−h(ˆ xt|t−1)▷Measurement residual 10: λt=Gσ∥zt∥M−1 t▷Scale factor 11: Kt=ˆ Pt|t−1λtHT tHtˆ Pt|t−1λtHT t+Mt−1 ▷Kalman gain 12: ˆ xt=ˆ xt|t−1+Ktzt▷Updated state estimate 13: Pt= (I−KtHt)ˆ Pt|t−1▷Updated covariance 14: end for 33