Genetic Algorithm-based Series-Parallel Mode NARX Model (GABPM-NARX) for Remaining Useful Life Estimation

1. Introduction

Both in production and service industries, prognostics, and health management (PHM) of engineering systems is very important for safe, reliable, and cost-efficient operation. Failure of the equipment during operation may have serious consequences. In terms of safety, there may be serious injuries, loss of life, and damage to the environment. In terms of financial losses, companies incur corrective maintenance costs (which are usually much higher than preventive maintenance costs) and opportunity costs due to the stoppage of production or services. Prognostics deals with estimating the time to failure of a system which is called as remaining useful life (RUL) in the literature [1]. Estimation of an engineering system’s RUL enables a company to plan the maintenance workforce and spare part supply resulting in more effective and efficient health management.

Prognostics models for RUL estimation are classified as physics-based, data-driven, and hybrid models [2]. Physics-based models express the system degradation by mathematical equations exploiting the system knowledge [3]. This approach requires expertise in system structure and dynamics, and it could be very difficult to construct mathematical equations representing the system degradation. Also, these models are system-specific so they cannot be generalized for other systems. Data-driven models do not require expertise in the system domain, and they can be generalized for other systems. These approaches utilize sensor technology to collect data that reflects the system deterioration. In the last two decades data-driven methods have been used extensively in PHM [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. The downside of the data-driven approaches is that they require a lot of data. In the case of limited data, hybrid models can be used for estimation of RUL. Hybrid models usually utilize data driven modelling for estimating the system state based on the sensor measurements (measurement model) then system degradation is expressed by physics-based models to estimate the RUL (degradation model) [21,22]. For the review of the topic see [22,23,24]. In this manuscript for a sufficient data case, a new data-driven methodology is proposed for the estimation of RUL utilizing Nonlinear Auto Regressive neural network with eXogenous inputs (NARX). NARX is a special kind of recurrent neural network. The proposed methodology is applied to NASA C-MAPSS data which is a benchmark dataset from the NASA data repository [25]. Although the data belongs to turbofan engines, the developed methodology can be generalized for any system.

Artificial Neural Networks (ANNs) are one of the data-based prognostics commonly used in literature. ANNs are computational algorithms inspired by the complex functionality of the human brain, where interconnected neurons process information in parallel [26]. These models do not use traditional statistical methods, but they are built as massively parallel networks, and trained to predict outputs [27]. ANNs are used effectively in modeling nonlinear relations and dynamical systems [28]. ANN is a very successful instrument in modeling the relationship between inputs and outputs in complex systems where the obvious exact relationship between inputs and outputs in the data is unknown [29]. With its ability to learn effectively from supervised and unsupervised input data, it can be applied to various fields such as modeling, time series analysis, pattern recognition, signal processing, and control [30,31]. ANNs can be divided into two categories: feedforward neural networks (FNNs) and recurrent neural networks (RNNs) [32]. The difference between them is that an RNN has a feedback mechanism between neurons, while an FNN does not. There would be one or more feedback connections between RNN neurons. The feedback mechanism can be quite diverse, from one neuron feeding back to another, a neuron feeding back to more than one neuron, more than one neuron feeding back to a neuron, or even a neuron feeding back to itself. The state of a network at any given moment depends on its previous state, which is achieved by feedback [33]. RNNs make predictions with dynamic modeling by iteratively simulating the dynamic behavior underlying nonlinear time series with an autonomous system with multi-step long-term predictions [34]. This structure is very useful for the analysis of dynamic systems. Since an engineering system and the sensors that monitor its state have a very dynamic and complex relation, we propose a new RUL prediction methodology based on the NARX model. The contributions of this paper can be summarized as follows:

• The predictive power of the NARX for solving dynamic and complex problems is recognized and proposed as a promising new model for the estimation of RUL.
• The hidden layer size and the delay size are two important parameters for the NARX architecture. In literature, as of our knowledge, these parameters are determined by trial. We utilized the genetic algorithm (GA) to determine the optimal values for hidden layer size and delay size, so the proposed model is called the genetic algorithm-based series-parallel mode NARX (GABPM-NARX) model.
• The developed GA fitness function for parameter optimization is independent of the model form so that it could be used for optimizing the parameters of other data-driven RUL estimation models in the literature.
• The proposed GABPM-NARX framework exhibits superior prediction accuracy in Remaining Useful Life (RUL) estimation compared to existing methodologies in literature. Owing to its robust performance and generalization capabilities, this method is poised to play a significant role in future RUL prediction research and applications. Its distinctive architecture empowers the model to extract enriched representations from complex degradation patterns, positioning GABPM-NARX as a promising benchmark in prognostics and health management.

To decrease the negative effect of the sensor noise on the performance of the NARX model, median filtering is utilized. The proposed median filtering approach also can be utilized in the data preparation stage for other existing models in literature.

The remainder of this paper is organized as follows. Section 2 provides a review of the literature. For fast access to concise information, the most crucial aspects of literature are summarized in a table. Section 3 provides a brief review of the NARX model and GA as preliminary. A detailed explanation of the proposed methodology is given in Section 4. A numerical experiment of the proposed methodology run on the C-MAPPS data set and given in Section 5. Finally, Section 6 provides conclusions.

2. Related Work

This section provides a concise review of the current literature on RUL estimation that utilizes the NASA C-MAPSS turbofan engine data. We decided to give the literature review in a table format since it gives the most important information in a concise form and makes the comparison of literature easy (Table 1). The papers in Table 1 are compared in terms of the data preprocessing steps, the framework used in modeling, and the metrics used for performance evaluation. In addition, the papers in the literature were analyzed in terms of keywords by years in Figure 1. The network density diagram shows that the most frequently used keywords in the papers are remaining useful life, deep learning, prognostics and health management, prognostics, C-MAPSS, predictive maintenance, prediction, Long Short-Term Memory (LSTM), turbofan engine, convolutional neural network (CNN), and machine learning.

3. Preliminary

3.1. NARX Model

The Nonlinear Auto Regressive neural network with eXogenous inputs (NARX) model is a type of RNN that can capture the nonlinear and dynamic relationships between input and output variables. In time series it is natural to have long-term dependencies. Auto Regressive part of the NARX allows the output depend not only on the previous outputs but also on the current and past inputs. Due to this feature, NARX is successful in predicting time series [6,40,41]. The NARX model can be used also for system identification [42], control [43] and signal processing [44] applications. It can handle complex data with high accuracy and robustness. The ability of the NARX model to maintain predictive or actual time series of historical data affects its performance [45]. NARX converges faster than other neural network models, generalizes well, and shows effective learning [46]. In addition, the NARX model belongs to the model class that is computationally equivalent to Turing machines [47].

The basic element of any kind of ANN is a neuron. The structure of a neuron is shown in Figure 2. Here the $k^{th}$ neuron has $n$ input arcs, and each arch has a weight $w_{ki}$ associated with it. Weighted input signals $X_1, X_2,\dots ,X_n$ are combined and a bias $b_k$ added to the sum to obtain intermediate output $v_k={\sum }_{i=1}^n{w}_{ki}{X}_i+b_k$ . The final output $Yk$ is obtained by applying an activation function $\phi (.)$ to${ v}_k$. The NARX model is formed by the interconnected layers of neurons in a special recurrent architecture.

For a general NARX model, input and output relation can be defined by the formula in Equation (1):

$${\widehat{y}}_t=f(y_{t-1}^{\text{'}},y_{t-2}^{\text{'}},\dots ,y_{t-n_y}^{\text{'}},x_{t-1},x_{t-2},\dots ,x_{t-n_x})$$

(1)

While $y_{t-1}^{\text{'}},y_{t-2}^{\text{'}},\dots ,y_{t-n_y}^{\text{'}}$ represents either real outputs or network outputs depending on the architecture, $x_{t-1},x_{t-2},\dots ,x_{t-n_x}$ represents network inputs. Time delay parameters expressed by $n_y$ and $n_x$ respectively indicate the number of past outputs and past inputs to be used for feedback, and $t$ denotes the discrete time step. Thus, the feedback and input layers are time-delayed and contain historical data [47]. The network structure of a general NARX model is shown in Figure 3. The model comprises an input layer, a hidden layer, and an output layer.

NARX models have two different architectures: series-parallel (also called as open loop) (Figure 4a) and parallel (also called as closed loop) architecture (Figure 4b). In the series-parallel architecture, the real output values are given as additional inputs whereas in the parallel architecture predicted output values are given as additional inputs. Formulation of series-parallel and parallel architecture is given in Equations (2) and (3) respectively:

$${\widehat{y}}_{t+1}=f(y_t,y_{t-1},y_{t-2},\dots ,y_{t-n_y},x_t,x_{t-1},x_{t-2},\dots ,x_{t-n_x})$$

(2)

$${\widehat{y}}_{t+1}=f({\widehat{y}}_t,{\widehat{y}}_{t-1},{\widehat{y}}_{t-2},\dots ,{\widehat{y}}_{t-n_y,}{x}_t,x_{t-1},x_{t-2},\dots ,x_{t-n_x})$$

(3)

In Equation (2) $y_t$ is the real output value at time $t$and in Equation (3) $\widehat{y}\left(t+1\right)$ is the prediction of the output for time $t+1$ at time $t$. $x(.)$is input values, and $f(.)$ is the mapping function. As it can be seen from Equation (2), in the series-parallel architecture ${\widehat{y}}_{t+1}$ is predicted using past real output values with time delay $n_y$ and input values with time delay $n_x$. As opposed to in parallel architecture ${\widehat{y}}_{t+1}$ is predicted using past predicted output values with time delay $n_y$ and input values with time delay $n_x$ as shown in Equation (3). The series-parallel and parallel architectures are shown in Figure 4.

The advantage of the series-parallel architecture is that the accuracy of the prediction increases by the use of the actual output values also as inputs. For example, if the problem is to estimate next year's wheat yield, past years' yield values can be used. RUL estimation appropriate NARX architecture is series-parallel architecture.

3.2. Population-Based Metaheuristics Method: Genetic Algorithm (GA)

The genetic algorithm (GA) is a population-based metaheuristic search algorithm inspired by evolutionary theory [48,49]. It focuses on the selection of suitable individuals by natural selection from a population of characters (a set of potential solutions). This process is based on transferring the good genes of the current generation to the future generations so that the chromosomes of future individuals are shaped accordingly. GA iteratively constructs the chromosome structures of future individuals from the good genes in the chromosome structures of ancestors. It is a good method for solving difficult combinatorial optimization problems. In this way, GA optimizes the model expressed around a problem-specific fitness function. GA consists of a population with a character set, control parameters, fitness function, crossover, mutation, selection, and coding operations [50]. The features of GA use a solution space in the search, provide probabilistic transitions, work with a parameter set, and focus on the return of the fitness function [51].

GA uses crossover and mutation functions to generate appropriate chromosome structures that optimize the fitness function. Crossover is the function to generate daughter solutions, that is, new chromosomes, by randomly selecting a crossover point from two parent solutions. Crossover GA is also one of the most important functions for generating suitable new chromosomes. The mutation function, on the other hand, ensures that the algorithm does not get stuck in the local optimal solution, so that genetic diversity is preserved, and the valuable information contained in the chromosomes is not prematurely lost (Grossmann & Sargent, 1978).

4. Proposed methodology for RUL estimation

4.1. Flowchart of the Proposed Methodology for RUL Estimation

The flowchart of the proposed methodology for RUL estimation is shown in Figure 5. First, feature selection is performed based on the raw data set to identify the most relevant features that give information about the degradation of the system. For this purpose, graphical and correlation analyses are utilized. The selected features are normalized using the z-score normalization method. Then, the RUL column is rectified with the piece-wise linear function. Finally, the noise in the data is removed by the median filter method. These steps are explained in detail and exemplified in the numerical example section. The processed training set is transferred to the NARX model. The parallel architecture is selected, and the GA is used to determine the number of hidden layers in the model architecture. Finally, the trained NARX model is used for the estimation of RUL values.

4.2. Genetic Algorithm-Based Series-Parallel Mode NARX Model (GABPM-NARX)

Mathematical formulation of the parallel mode NARX model is developed hierarchically starting from the neuron which is the basic element of the architecture. Let us assume that the architecture has m hidden layer, and each layer has s neurons. The formulation of $j^{th}$ neuron of $l^{th}$ hidden layer for time $t$ is given in Equation (4).

$$n_{jt}^l=\gamma ({\beta }_{j0}+\sum _{p=1}^n{\beta }_{jp}{x}_{pt}+{\rho }_{jp}(x_{p\left(t-1\right)}+\dots +x_{p\left(t-n_x\right)})+{\rho }_{n+1}({\widehat{y}}_{\left(t-1\right)}+\dots +{\widehat{y}}_{(t-n_y)})j=1,\dots ,s$$

(4)

$$\gamma \left(a\right)={\displaystyle \frac{1}{1+e^{-a}}}$$

(5)

In Equation (4), $\gamma$ is the logistic activation function and its formulation is given in Equation (5). $j=1,\dots ,s$ and $p=\mathrm{1,2},\dots ,n$ represents the indexes for neuron and exogenous variable sizes respectively. ${\beta }_{j0}$ is the input bias. $x_{pt}$ is the value of $p^{th}$ exogenous variable at time $t$ and ${\beta }_{jp},$ are the input layer weight of it at neuron $j$. $x_{p(t-1)},\dots , x_{p(t-n_x)}$ represents the time delayed value of exogenous variable and ${\widehat{y}}_{(t-1)},\dots ,{\widehat{y}}_{(t-n_y)}$ represents time delayed estimated output. ${\rho }_{jp}$ and ${\rho }_{j(n+1)}$ represent the input layer weights for $p^{th}$time delayed exogenous variable and time delayed estimated output, respectively.

A hidden layer is composed of s neurons. The formulation of $l^{th}$ hidden layer at time $t$ is given in Equation (6).

$$h_t^l=\sum _{j=1}^s{\delta }_j{n}_{jt}^l, l=1,\dots ,m{h}_t^l=\sum _{j=1}^s{\delta }_j\gamma \left({\beta }_{j0}+\sum _{p=1}^n{\beta }_{jp}{x}_{pt}+{\rho }_{jp}{x}_{j(t-n_x)}+{\rho }_{n+1}{\widehat{y}}_{(t-n_y)}\right)l=1,\dots ,m$$

(6)

Here ${\delta }_j$ is the weight of $j^{th}$ neurons. Finally, the output layer combines all the hidden layers and it is formulation given in Equation (7).

$${\widehat{y}}_t=\theta \left[{\delta }_0+\sum _{l=1}^m{{\omega }_{l}h}_t^l\right]$$

(7)

Here $\theta$ is the activation function for the output, which is $\theta \left[x\right]=x$. The output bias is ${\delta }_0$ and the output layer weights are ${\omega }_l$. Inserting Equations (4) and (6) in Equation (7), we obtained Equation (8) which is the expression of the parallel mode NARX model.

$${\widehat{y}}_t=\theta \left[{\delta }_0+\sum _{l=1}^m{\omega }_l\sum _{j=1}^s{\delta }_j\gamma \left({\beta }_{j0}+\sum _{p=1}^n{\beta }_{jp}{x}_{pt}+{\rho }_{jp}{x}_{j(t-n_x)}+{\rho }_{n+1}{\widehat{y}}_{(t-n_y)}\right)\right]$$

(8)

Hidden layers learn the relationship between inputs and output. Therefore, the size of hidden layers affects the performance of the NARX model. One of the other important parameters is the delay size that defines the auto-regressive part. In previous studies that utilized NARX model, the hidden layer and the delay size were determined manually by trial and error [45,52]. In this study, GA is used for determining the optimal/near-optimal values of these parameters.

To perform the optimization by GA, a special fitness function is constructed. The fitness function of GA is called as $RMSE\_(sample\_(NARX\_(fitness\_func ) ) )$, and its calculation formula is expressed in Equation (9).

$$RMS{E}_{sampl{e}_{NAR{X}_{fitnes{s}_{func}}}}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left({RUL}_i^{estimated\_sample}-{RUL}_i^{actual\_sample}\right)}^2}$$

(9)

Here, ${RUL}_i^{estimated\_sample}$and ${RUL}_i^{actual\_sample}$ are the predicted RUL and actual RUL values of the $i^{th}$ row of the training dataset, respectively. $N$is the total number of samples in the training data. ${RUL}_i^{estimated\_sample}$ is found utilizing Levenberg-Marquardt learning algorithm which is the popular algorithm in literature [53,54].

Once the GABPM-NARX model is constructed the training data is fed into the model to learn the relation between inputs and output. Then trained model can be used for estimating RUL values.

5. Experimental Study

5.1. C-MAPSS Dataset

The dataset provides the run-to-failure trajectories of a fleet of turbofan engines. It is generated utilizing the Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) model developed at NASA [25], [55]. The dataset includes four different subsets (FD001, FD002, FD003, and FD004) that are simulated under different combinations of operating conditions and failure modes. Table 2 provides detailed information on the data subsets. The data subsets consist of multivariate time series that record the information in the sensor channels to be able to characterize the failure evolution. Each subset is divided into training and test sets. The engines always run normally at the beginning of the time series and develop a fault at some point during the series. In the training set, the error grows until it becomes a system error. In the test set, the time series ends sometime before the system failure. The goal is to estimate the RUL on the C-MAPSS test set. In addition, a vector of actual RUL values is provided for the test data. The data is presented in 26 columns. In the dataset column 1 gives the engine ID number, column 2 gives the running cycle, column 3-5 gives operational settings, and column 6-26 gives the 21 different sensor measurements. Each row of the data set is a snapshot of the data obtained during a single processing cycle.

A diagram of the turbofan engine is shown in Figure 6. Key components include a fan, low-pressure compressor (LPC), high-pressure compressor (HPC), fuel oil, high-pressure turbine (HPT), low-pressure turbine (LPT), and nozzle. Sensors monitor the health of the key components of the turbofan engine. The description of the 21 sensors is shown in Table 3.

We used the FD001 data set to show the effectiveness of the proposed methodology. It is the most used data set in the literature [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20].

5.2. Data Preparation

Data preparation steps in order are feature selection, normalization, RUL label rectification, and noise-cleaning. Details of these steps are given below.

5.2.1. Feature Selection-Graphical and Correlation Analysis

The measurement data of some sensors do not change at all during the life of the turbofan engine or change around fixed values, so they do give none or very little information about the health state of the engine. Including irrelevant sensors in the study would increase variation and the noise thereby decreasing the prediction accuracy. These sensors should be removed by feature selection methods. In this study, feature selection is performed by graphical and correlation analysis.

Figure 7 gives examples (due to space limitation) of time series graphs showing the different behaviors of the sensor measurements. Measurements of sensors 1, 5, 10, 16, 18, and 19 stay constant through the whole engine life as in Figure 7a. Measurements of sensor 6 change between 21.6 and 21.61 without any trend and is given in Figure 7b. Measurements of sensors 2-4, 7-10, 11-15, 17, and 20-21 show a trend as in Figure 7c.

Pearson correlation coefficient analysis is a filtering method for determining exclusions from a dataset based on the degree of correlation. Pearson correlation coefficient measures the linear relation between two variables and is expressed in Equation (10). As $r$ gets closer to +1 positive linear relation between two variables increases, as $r$gets closer to -1 negative relation between two variables increases, and as r get closer to 0 relation between two variables decreases. Table 4 gives Pearson correlation coefficients between RUL and each sensor. Table 4 supports the results of the graphical analysis. The correlation between RUL and constant sensors is zero. The correlation between RUL and Sensor 6 is -0.11. Correlation between RUL and sensors with a trend ranging from -0.78 to 0.75. As a result of two analyses, sensors 2,3,4,7,8,9,11,12,13,14,15,17,20, and 21 are selected as features.

The age of an engine (the running cycle) also provides important information for the estimation of the RUL. The Pearson correlation coefficient between RUL and the running cycle is calculated as -0.746 indicating that the running cycle should be added to the feature set. As a result, the feature set included 14 selected sensor values and the running cycle.

$$r={\displaystyle \frac{{\sum }_{t=1}^n(X_t-\overline{X})(Y_t-\overline{Y})}{\sqrt{{\sum }_{t=1}^n{(X_t-\overline{X})}^2}\sqrt{{\sum }_{t=1}^n{(Y_t-\overline{Y})}^2}}}$$

(10)

5.2.2. Z-Score Normalization

The scale of the measurement values of the sensors are different from each other, and this will negatively impact the model training process. Therefore, the sensor data must be normalized. There are many ways to normalize variables, in this study z-score normalization, the most popular technique, was utilized [56,57]. The formula of z-score normalization is shown in Equation (11) [58]. Here $x_{ij}$ is the $c$ measurement of $i^{th}$ sensor, ${\mu }_i$ is the mean of $i^{th}$ sensor measurements, ${\sigma }_i$ is the corresponding standard deviation, and $x_{ij}^{\text{'}}$ is the normalized value of the $j^{th}$ measurement of $i^{th}$ sensor.

$$x_{ij}^{\text{'}}={\displaystyle \frac{x_i-{\mu }_i}{{\sigma }_i}}$$

(11)

5.2.3. Piece-Wise Linear Function

RUL labels for the training dataset start with maximum life of the engines and decrease by one for each flight cycle until the failure occurs. When failure occurs, the RUL label becomes zero. Accordingly, the shape of RUL is a time-decreasing linear function. But in practice, the deterioration of the system is almost constant at the beginning, and it starts decreasing with the first signs of engine degradation. Reflecting on this observation, a piece-wise linear RUL function has been proposed in [59] and adopted in the literature [15,60,61]. RUL label rectification formula is expressed in Equation (12). Here ${RUL}_{early}$ represents the constant RUL value in early life. In this study ${RUL}_{early}$ is set to 125 cycle as suggested in [62]. Based on this ${RUL}_{early}$ value, graph of the piece-wise linear RUL function is shown in Figure 8.

$$RUL=\left\{\begin{array}{c}{RUL}_{early}, if RUL\ge {RUL}_{early}\\ RUL, if RUL<{RUL}_{early}\end{array}\right.$$

(12)

5.2.4. Noise Cleaning with Median Filter

Noise in the sensor measurements decreases the model performance. Median filtering, developed by [63], is a useful non-linear signal processing technique for regulating noise in data. It is a popular filtering method in image processing [64,65,66]. In image processing, the median filter is applied to the data in sliding windows of odd size. By filtering the center of the window is replaced by the median of the window. Figure 9 shows the a median filter and a mean (smoothing) filter applied to discrete step function, a ramp function, a pulse function, and a triangle function with a window of five pixels [67]. It can be seen from the Figure 9 that median filter does not affect the step function and ramp functions. For the pulse functions, pulses whose periods are less than one-half of the window width, are suppressed. For the triangle function it flattens the peak of the triangle.

In this study, median filtering is used to regulate the noise in FD001 dataset. Window size is taken as a whole data size for not losing important information. As an example, Figure 10 shows the signal distribution for Sensor 3 and Sensor 15 with and without median filter. It can be seen that median filter smooths the noise without changing the shape of the signal trajectories.

5.3. Experimental Setup

Let us summarize the preparation of FD001 dataset, given in section 5.2. In the first step sensors 2,3,4,7,8,9,11,12,13,14,15,17,20, and 21 are selected using graphical and correlation analysis. Based on its high correlation the running cycle is added to the feature set. Time series plots of the selected sensors are given in Figure 11 In the second step normalization is performed applying the z-score method to the sensors. The running cycle is not normalized. In the third step RUL labels are rectified by setting ${RUL}_{early}$ as 125 flight cycle. In the fourth step median filtering is utilized to clean the noise in the sensor data. Although there is a separate training data set, to evaluate the internal performance of the model testing data set is divided into 5% validation and 5% testing in the fifth step.

Time delay size is determined by examining the sensor values. The sensors values are approximately similar in the 4 cycle-time windows, as shown in Figure 12. From another point of view, when examined within the framework of maintenance logic, a time period of about 4 cycles can inform the malfunction.

Hidden layer size and delay size determined by the GA using the fitness function given in Equation (9). The GA parameters are expressed in Table 5.

As stated in Section 4.2, Levenberg-Marquardt learning algorithm is used for the estimation of RUL_i^{estimated_sample} values in Equation (9). Figure 13 shows the fitness function value of the GA. In the initial generations, there was a rapid improvement in fitness function, and in the next generations, appropriate value was obtained with minimal improvements. Optimum value of hidden layer size determined as 10 and delay size determined as 4.

Next parallel NARX model is constructed. Model parameters of the NARX Model (time delay size, size of hidden layer, type of learning algorithm, and error formula) are given in Table 6.

Using 10 hidden layers, the proposed GABPM-NARX model is constructed, and its architecture is given in Figure 14.

Once the GABPM-NARX model is constructed, it is trained using three different learning algorithms, the Levenberg-Marquardt Method, Scaled Conjugate Gradient Method, and Bayesian Regulation Method. The best results were obtained when the Levenberg-Marquardt Method was used, and only the result of this model is given below.

5.4. Evaluation of Training Process

The training data set can be represented in matrix format as follows. There are 20631 rows of training data belonging to 100 engines. Each row represents one flight cycle. The values of 14 sensors are selected as features. Let ${\tau }_i$ be the lifetime of the system $i$, $i=1, 2,\dots , 100$. The maximum lifetime is 362 flight cycles belonging the engine number 69 and the minimum lifetime is 128 flight cycles belonging to engine 39. Sensor measurements of system $i$ during its lifetime ${\tau }_i$ can be represented by the matrix $X_i{=[x}_1^i x_2^i\dots .x_m^i]$, $X_i\in R^{15\times {\tau }_i}$ . Here $x_j^i$ is the $1\times {\tau }_i$ vector of $j^{th}$ sensor measurements for $i^{th}$ system during its lifetime ${\tau }_i$. So, the whole data set is a matrix with dimension 15$\times$

20631. The dataset was partitioned into training, validation, and test subsets, with respective allocations of 60%, 20%, and 20%.

Using this data, the GABPM-NARX model is trained using the Levenberg-Marquardt algorithm to learn the relation between inputs and output. The regression determination coefficient of the simple linear regression model between output values $\left({\widehat{y}}_i\right)$ and target values $y_i$ can be used as the first metric to evaluate the success of the model training process. As is well known, R can take values in the interval [0,1]. The closer the value of R to 1, the better the training process is. The regression results are shown in Figure 15. The R-value is approximately 0.9973, indicating that the training process was successful.

As a second metric for the evaluation of the training process, the Mean squared error (MSE) is utilized. MSE is the mean of the square differences between output and target ($MSE={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{\left({y_i-\widehat{y}}_i\right)}^2)$. Low MSE values indicate successful training. Progression of MSE as a function of training iteration (epoch) g is given in Figure 16. Model performance was assessed based on the Mean Squared Error (MSE) across training, validation, and test sets over 31 epochs, as presented in Figure 16. The training curve exhibited a consistent decline in MSE, reflecting effective learning and convergence. The validation curve showed moderate fluctuation and reached its lowest value (24.2886) at epoch 25, indicating optimal generalization capacity. The test curve closely followed the validation trend and remained stable throughout, supporting the model’s robustness on unseen data. The slight deviation observed beyond epoch 25 suggests early signs of overfitting, reinforcing the suitability of epoch 25 as the preferred checkpoint for model deployment. Overall, these results demonstrate adequate convergence and generalization, confirming the reliability of the model for predictive applications.

The error, $\left({y_i-\widehat{y}}_i\right),$ an error histogram was generated to analyze the distribution of prediction errors at a sampling frequency of 20 Bins (Figure 17). The x-axis represents discrete error intervals, while the y-axis indicates the number of occurrences per interval. A prominent peak is observed at the zero-error bin, exceeding 2.0 × 10⁶ instances. This pronounced concentration suggests that the majority of predictions fall precisely on the expected values, indicative of a highly accurate modeling process. The surrounding error bins exhibit negligible frequencies, further reinforcing the model’s precision and minimal bias. The histogram’s narrow distribution and central peak support the hypothesis that the system's residuals are tightly clustered around zero, confirming consistent and reliable performance under the given operating conditions.

Residual diagnostics were performed to assess the adequacy of the time series model. As depicted in Figure 18, the autocorrelation function (ACF) of the residuals revealed a significant spike at lag 0, indicating strong instantaneous dependence. However, autocorrelations at other lags were negligible and remained within the 95% confidence bounds, suggesting the absence of serial correlation. This pattern implies that the residuals exhibit characteristics consistent with white noise, supporting the model’s effectiveness in capturing temporal structure.

Model training dynamics were evaluated over 31 epochs using gradient magnitude, adaptive mu parameter, and validation failure count as diagnostic indicators (Figure 19). The gradient curve exhibited a high initial value with a declining trend, stabilizing toward later epochs (final value: 310.181), which indicates effective convergence behavior. The mu parameter showed notable fluctuations, with a pronounced drop around epoch 15 and stabilizing at a minimal value of 0.01 by epoch 31—suggesting progressive adjustment of the learning rate. Validation failure count remained minimal during the initial training phase but escalated in the final third of the epochs, reaching 6 failures at epoch 31. This late-stage rise may indicate the onset of overfitting or instability, warranting further hyperparameter tuning or early stopping mechanisms. Collectively, these metrics provide critical insight into the model’s convergence and generalization capacity.

5.5. Evaluation of the Test Process

After the training process, FD001 test data is fed into the GABPM-NARX model. The root-mean-square error (RMSE) and Score measurements were used to compare the performance of the method with other methods in the literature.

5.5.1. The Root-Mean-Square Error (RMSE)

The RMSE, which is widely used as a metric to evaluate regression models, expresses the deviation between the estimated value and the actual value. The RMSE gives the same penalty value for the under estimation and over estimation. The calculation formula is given in Equation (13). N is the total number of samples in the test data.

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left({RUL}_i^{estimated}-{RUL}_i^{actual}\right)}^2}$$

(13)

5.5.2. Score

The score measurement function calculates the penalty in case of under and over estimation by different formulas. As it can be seen from Equation (14), Score measurement penalizes over estimation more [25]. The reason is that the overestimating RUL can cause serious results.

$$Score=\sum _{i=1}^N{SF}_i,{SF}_i=\left\{\begin{array}{c}{e}^{-({RUL}_i^{estimated}-{RUL}_i^{actual})/13}-1,{(RUL}_i^{estimated}-{RUL}_i^{actual})<0\\ e^{({RUL}_i^{estimated}-{RUL}_i^{actual})/10}-1,{(RUL}_i^{estimated}-{RUL}_i^{actual})\ge 0\end{array}\right.$$

(14)

Metric sensitivity to error magnitude was evaluated by plotting Score and Root Mean Squared Error (RMSE) functions against varying error values (Figure 20). The Score Function exhibits a linear V-shaped profile, symmetric around the origin, indicating a direct proportionality between the metric and the absolute error. In contrast, the RMSE Function follows a parabolic trajectory, with its minimum centered at zero error, reflecting quadratic sensitivity to error magnitude. The linear behavior of the Score Function emphasizes uniform penalization regardless of error direction, whereas the RMSE’s curvature amplifies larger deviations. This comparative analysis highlights the differential weighting mechanisms embedded within each metric, offering insight into their respective implications for model optimization and error evaluation.

5.5.3. Results

The RUL estimation results of 8 randomly selected engines are shown in Figure 21. It can be seen that GABPM-NARX model successfully catches the trend of actual rule values. To intuitively validate the effectiveness of the proposed GABPM-NARX framework at different regression complexities, we present comparisons of the estimated and actual RUL values at different degrees for different engines on a sample basis (Figure 21). The GABPM-NARX model estimates consistently show better agreement with the actual RUL curves. This performance can be attributed to the spatio-temporal feature extraction and fusion mechanism, which utilises temporal degradation models through both sensor-level correlations via self-observation and levels appropriately determined by NARX in GA. These example predictions demonstrate the superior generalisation capacity and robustness of the proposed GABPM-NARX across different data complexities and model configurations.

$RMSE$ and $score$ values of the GABPM-NARX model are 3.0331 and 12.5216 respectively. These metrics are calculated for the end of the observation period for all engines. Comparison of these values with the literature is given in Table 1. In terms of $RMSE$, the GABPM-NARX model was better than literature.

6. Discussion

The proposed methodology, which combines median filtering, NARX modeling and genetic algorithm-based parameter tuning, represents an innovative pipeline for the prediction of remaining useful life (RUL). The median filter was found to be effective in reducing noise in the sensor data, an essential pre-processing step that is often undervalued in existing prediction studies. This is consistent with previous work that emphasizes the importance of robust data cleaning to improve signal integrity and model performance.

The use of the NARX model takes advantage of the temporal dynamics inherent in degradation trends, as time series data typically exhibits memory-dependent behavior in RUL estimation. In contrast to traditional feedforward networks, NARX structures are better suited to capture such dependencies. Our results support the growing body of literature suggesting that recurrent architectures, especially those with external feedback loops, are useful in prognostic health management (PHM) systems.

While previous research often determines the architecture of the model empirically, our third contribution — a genetic algorithm-based search for the optimal hidden layer size and time delay — introduces a more systematic and scalable optimization strategy. This optimization approach is not only model-independent but also generalizable, so that it could benefit a larger number of data-driven forecasting systems beyond NARX.

Although the GABPM-NARX model does not outperform all benchmark methods on the FD001-C-MAPSS dataset, its performance confirms its suitability as a predictive model and highlights several opportunities for refinement. These include exploring ensemble-based schemes, hybridizing physics-informed constraints with data-driven models, and incorporating transfer learning strategies to account for dataset shifts under different operating conditions or machine types.

Future research could also address adaptive or real-time tuning of NARX configurations as degradation progresses to maintain prediction accuracy in transient environments. Furthermore, cross-validation across multiple C-MAPSS subsets would strengthen the claim of generalizability and pave the way for use in real-world scenarios.

Overall, this study highlights the potential of integrating evolutionary optimization with dynamic neural networks and creates a foundation for more intelligent and adaptive predictive modeling systems.

7. Conclusions

This paper presents a methodology for predicting RUL based on NARX models and GA. As the first contribution, we proposed to utilize median filtering to clean the noise from data in the data preparation step. Since long-term dependencies are naturally observed in time series, we proposed to utilize the NARX model for RUL prediction. In literature hidden layers and delay size of the NARX model are determined experimentally. As our third contribution, we find the optimal hidden layer and delay size for the NARX model utilizing GA. The fitness function developed for the GA is not dependent on the model form so it can be used for parameter optimization of other data-based RUL estimation models.

The performance of the proposed methodology is shown on the FD001-C-MAPPS dataset. Although the performance of the GABPM-NARX model is not the best in the literature, this study shows that the NARX model is promising in prognostic studies. Future research will be conducted to improve the prediction accuracy of the proposed GABPM-NARX model.