Turbofan Engine Health Prediction Model Based on ESO-BP Neural Network

Abstract

Predicting the health index of turbofan engines is critical in reducing downtime and ensuring aircraft safety. This study introduces the elite snake optimizer-back propagation (ESO-BP) model to address the challenges of low accuracy and poor stability in predicting the health index of turbofan engines through neural networks. Firstly, the snake optimizer (SO) was improved into the elite snake optimizer (ESO) through an elite-guided strategy and a reverse learning mechanism. The performance improvement was validated using benchmark functions. Additionally, feature importance was introduced as a feature selection method. Finally, the optimization results of the ESO were employed to set the initial weights and biases of the BP neural network, preventing convergence to local optima. The prediction performance of the ESO-BP model was validated using the C-MAPSS datasets. The ESO-BP model was compared with the CNN, RNN, LSTM, baseline BP, and unimproved SO-BP models. The results demonstrated that the ESO-BP model has a superior accuracy with an impressive R-squared ($R^2$) value of 0.931 and a root mean square error ($RMSE$) of 0.060 on the FD001 sub-dataset. Furthermore, the ESO-BP model exhibited lower standard deviations of evaluation metrics on 100 trials. According to the study, ESO-BP demonstrated a greater prediction accuracy and stability when compared to commonly used models such as CNN, RNN, LSTM, and BP.

1. Introduction

Turbofan engines are one of the core components of an aircraft [1]. Predicting the remaining useful life (RUL) of the engine and implementing prognostics and health management (PHM) are critical for formulating a targeted maintenance strategy [2]. The turbofan engine health prediction model helps to reduce the risk of unplanned downtime and maximize the operational lifespan [3].

The prediction model for the RUL of turbofan engines is developed based on either their physical structure model or the degradation trends among various monitoring sensors. There are two types of existing methods for estimating RUL: those based on physical models [4,5] and data-driven approaches [6,7,8,9]. Data-driven methods analyze data and identify complex data relationships, achieving accurate predictions of changes in the turbofan engine’s state [10]. Data-driven predictions for the RUL of turbofan engines are influenced by various factors, including the number of operating cycles, operating parameters, and sensor data.

The back-propagation (BP) neural network, which was proposed by Rumelhart et al. [11] in 1986, is a feed-forward neural network trained with the error back-propagation algorithm. It is currently acknowledged as one of the most extensively employed neural network models, celebrated for its simple principles, high computational speed, and robust regression prediction capabilities [12]. However, BP neural networks frequently experience slow convergence and a tendency to become stuck in local optima due to the random selection of initial weights and thresholds [13]. Currently, there is a surge in the use of neural networks coupled with intelligent optimization algorithms. In 2022, Yu et al. [14] introduced a chaotic krill herd algorithm optimized BP neural network, achieving a 3.65% reduction in prediction error for kidney bean yield. Ji and Ding [12] enhanced BP neural networks by utilizing the improved Sparrow Search Algorithm, resulting in an R² value of 0.99 for sea surface temperature retrieval. Lai et al. [15] introduced the IA-PSO-BP model, which achieved an R² value of 0.97 for ground pressure prediction.

In the field of intelligent optimization algorithms, Hashim et al. [16] presented a novel algorithm named SO in 2022. However, this algorithm has persistent drawbacks, such as the random movement of individual positions, which may lead to local optima and poor stability. In 2023, Deng and Liu [17] introduced the snow ablation optimizer (SAO), which incorporates elite sets, demonstrating both rapid convergence and a high convergence accuracy.

This study presents an ESO-BP model for predicting the health of turbofan engines. The goal is to enhance the prediction accuracy and stability of the BP neural network in the context of predicting turbofan engine health. The subsequent sections are organized as follows: Section 2 introduces the elite snake optimizer (ESO). Section 3 elaborates on the proposed model architecture. Section 4 covers the experimental specifics, results, and analysis. Section 5 discusses the findings. Finally, Section 6 concludes the study.

2. Elite Snake Optimizer

Intelligent optimization algorithms, also known as optimization techniques or metaheuristic algorithms, are computational methods designed to iteratively explore solution spaces with the aim of discovering the optimal solution to a problem [18]. These algorithms are often inspired by natural processes or phenomena and strive to efficiently navigate complex and high-dimensional spaces in search of the best possible solution [18].

2.1. Snake Optimizer

The SO draws inspiration from the behavioral traits of snakes, specifically, their fighting, foraging, mating, and reproduction behaviors. The algorithm utilizes a group of snake agents to explore the search space and aims to achieve a global optimum based on the fitness values of individual agents. The algorithm also takes into account environmental factors that influence different snake behaviors, such as temperature and food quantity. The equations for updating the environmental factors are provided below [16].

$$Temp=e^{-\frac{t}{T}}$$

(1)

$$Q=C_1{e}^{\frac{t-T}{T}}$$

(2)

Here, $Temp$ represents the temperature, $t$ signifies the current iteration, $T$ denotes the maximum iterations, $Q$ stands for the food quantity, and $C_1$ is a constant value.

Unlike other algorithms, the SO takes a step further in emulating biological behavior. The SO directs the population individuals through unique positional transitions based on the current iterative phase during each iteration of the algorithm [19]. While there are bio-inspired intelligent optimization algorithms such as the sparrow search algorithm (SSA) [20], beluga whale optimization (BWO) [21], and dung beetle optimizer (DBO) [22], which only mimic the predation paths of organisms in their position updating, the SO improves its perceptual capabilities of the search space by simulating the physiological behaviors of snakes in response to environmental changes. This enhancement accelerates the convergence of SO and improves its precision in achieving convergence [19].

2.2. Elite-Guided Strategy

In intelligent optimization algorithms, elite-guided strategies involve selecting or generating a set of individuals with a higher fitness during the population initialization phase. These individuals subsequently influence the movement of others in the following iterations, enhancing the search efficiency and convergence speed of the swarm [23].

Various intelligent optimization algorithms use different elite-guided strategies. For example, in particle swarm optimization (PSO), each particle updates its velocity and position by considering its individual historical best position and the global best position [24]. In ant colony optimization (ACO), ants release pheromones while searching for a path, and subsequent ants select paths based on the concentration of pheromones [25]. The whale optimization algorithm (WOA) designates the current optimal solution as the leader for the population in each iteration, and the rest adjust their positions based on the leader’s location [26].

Elite-guided strategies aim to direct individuals towards regions with a superior fitness, accelerating the convergence process [27]. However, this may result in premature convergence, trapping the algorithm in local optima [27]. Therefore, designing an appropriate elite-guided strategy is a crucial aspect of research on intelligent optimization algorithms.

During the population initialization stage of ESO, an elite-guided strategy is employed. Elite subsets are selected from both the male and female snake populations, each comprising four individuals. The subsets consist of the top three individuals ranked by the fitness of their current position, along with a new individual calculated as the mean coordinate of all individuals in the current population.

Figure 1 illustrates that Elite 1, Elite 2, and Elite 3 expedite the convergence speed of SO in the search space, since they are the three agents closest to the optimal solution. Elite 4 is the new individual calculated as the mean coordinate of all individuals in the current population, which emphasizes the diversity of the population and helps SO to avoid local optima. Equation (3) is used to calculate Elite 4. This strategy aims to achieve a balance between convergence speed and accuracy, thereby enhancing the global search efficiency of SO.

$$x_{elite4,i}=\frac{{\sum }_{j=1}^n{x}_{j,i}}{n}$$

(3)

Here, $x_{elite4,i}$ is the $i$-th dimension value of Elite 4, $x_{i,j}$ is the $i$-th dimension value of $j$-th individual, and $n$ is the population size.

2.3. Reverse Learning Mechanism

Reverse learning mechanisms are utilized as an enhancement method for intelligent optimization algorithms. The idea is to use symmetric positions within the search space, such as point symmetry, axial symmetry, planar symmetry, or rotational symmetry, to generate reverse populations [28]. The reverse population complements the original solution, maximizes the utilization of the search space, and enhances the algorithm’s ability to overcome local optima [29]. The aim of this method is to improve the algorithm’s global search capability, thereby enhancing the comprehensiveness and effectiveness of the optimization process [29].

During the position update phase of ESO, a reverse learning mechanism is incorporated. This involves assigning opposite values to the positions of the new population within the limits of each dimension, resulting in a mirrored population. The reverse learning mechanism expands the existing population while enhancing the algorithm’s capability to overcome local optima. The calculation for the position of the reverse individual is presented in Equation (4) [28].

$$x_i^{*}=l_u+l_l-x_i$$

(4)

In the Equation (4), $x_i^{*}$ represents the value of the reverse individual for the $i$-th dimensional coordinate, $l_u$ is the upper boundary of the current search space, $l_l$ is the lower boundary of the current search space, $x_i$ is the value of the original individual for the $i$-th dimensional coordinate.

2.4. Benchmark Test Results

The SO has been enhanced into the elite snake optimizer through the improvements mentioned above. To evaluate the convergence accuracy and stability of the ESO, 30 trials were performed using the 30-dimensional benchmark functions F7, F18, and F30 from the 2017 IEEE Congress on Evolutionary Computation (CEC2017) [30]. Figure 2a–c show the corresponding three-dimensional surface plots with a contour plot underneath.

The CEC2017 test function set is a widely recognized and challenging tool for validating algorithm performance [30]. All functions in the test set are rotated and shifted, which increases the difficulty of algorithmic optimization search. The set consists of 29 test functions, all of which are rotated and shifted to increase the difficulty of the algorithmic optimization search. It is worth noting that the F2 function was officially removed from the original function set due to its instability.

Functions F1 and F3 have a single peak and no local minima, making them suitable for testing the algorithm’s convergence ability. Functions F4–F10 have multiple peaks with local extrema, making them suitable for testing the algorithm’s ability to escape local optima [30]. Functions F11–F20 are hybrids, composed of three or more CEC 2017 baseline functions that have been rotated or shifted. Each subfunction is assigned a weight. Functions F21–F30 are compositions of at least three hybrid functions or CEC 2017 baseline functions that have been rotated or shifted. Each subfunction has an additional bias value and weight, making these combined functions even more challenging to optimize.

The ESO was compared to eight mainstream intelligent optimization algorithms: differential evolution algorithm (DE) [31], PSO, grey wolf optimization (GWO) [32], chicken swarm optimization (CSO) [33], SSA, BWO, DBO, and SO. These algorithms were aimed toward obtaining the global minimum of fitness. The best fitness values of the optimization algorithm were recorded for each trial. The results, as shown in Figure 2b,d,f, demonstrate that ESO has lower means, maximums, and minimums, which indicate a better convergence accuracy. The smaller interquartile range represents the better stability of the ESO algorithm.

3. Model Construction

This study presents an ESO-BP neural network model for predicting the health index of turbofan engines. The model utilizes the optimization results of ESO to initialize the weights and biases of the BP neural network, which prevents it from converging to local optima and improves both prediction accuracy and stability. Figure 3 displays the complete flowchart of ESO-BP.

3.1. Health Index Degradation

The C-MAPSS dataset defines the remaining useful life of a turbofan engine as the number of remaining operating cycles until a fault occurs. Current research primarily centers on predicting the remaining operating cycles of turbofan engines [6,7,8,9,34,35,36]. However, the remaining operating cycles are influenced by various external and internal factors, resulting in significant errors and a low coefficient of determination between the predicted and actual values. Most prediction modeling studies that focus on the number of remaining operating cycles struggle to capture the degradation of an engine with a substantial RUL greater than 100 cycles.

This study differs from traditional $RUL$ prediction by utilizing a remaining health index ($RHI$) evaluation criterion for turbofan engines. The $RHI$ is constructed based on $RUL$, and the health of each brand-new engine is initially recorded as ‘1’. The RHI decreases by a certain amount for each operating cycle until the engine stops due to a malfunction, at which point, the $RHI$ becomes ‘0’ [37]. In simpler terms, the $RHI$ is the ratio between the current number of operating cycles remaining for the engine and the maximum number of operating cycles for this engine ($Ma{x}_c$), as calculated by Equation (5).

$$RHI=\frac{RUL}{Ma{x}_c}$$

(5)

Figure 4 illustrates a visualization that compares the $RHI$ degradation function with existing segmented $RUL$ degradation functions.

3.2. Introduction to Datasets

This study uses the publicly available turbofan condition monitoring data from the Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) dataset provided by the National Aeronautics and Space Administration (NASA) to predict the remaining health index of turbofan engines [38]. The C-MAPSS dataset comprises four subsets, each containing various operating and faulty working conditions. Each subset is then divided into a training set and test set, as shown in

Table 1. C-MAPSS dataset information.

Sub-Dataset	FD001		FD002		FD003		FD004
	Training Set	Test Set	Training Set	Test Set	Training Set	Test Set	Training Set	Test Set
Operating condition	1	1	6	6	1	1	6	6
Fault mode	1	1	1	1	2	2	2	2
Number of engines	100	100	260	259	100	100	249	248
Number of operating cycles	20,631	13,096	53,759	33,991	24,720	16,596	61,249	41,214
Mean of engine operating cycles	206.31	130.96	206.77	131.24	247.20	165.96	245.98	166.19
Standard error of engine operating cycles	0.32	0.47	0.20	0.34	0.55	0.67	0.30	0.45
Maximum of engine operating cycles	362	303	378	367	525	475	543	486
Minimum of engine operating cycles	128	31	128	21	145	38	128	19

.

The training set records various state parameters of the turbofan engine. These variables include the number of operating cycles and 24 pieces of sensor information. Data are recorded as time series throughout its complete lifespan cycle, from the normal state to the fault state. In contrast, in the test set randomly stops collecting the cycle parameters, meaning that the complete life cycle is not recorded for each engine.

These sensor data include information about the engine’s operating parameters, pressure, and temperature for 24 lots of sensor information, as shown in

Table 2. Feature importance of 4 sub-datasets. The underlined numbers represent the redundant features in the corresponding sub-datasets.

Feature Description	Units	Symbol	$\mathbf{Feature} \mathbf{Importance} {\mathit{C}}_{\mathit{i}}$
			FD001	FD002	FD003	FD004
Number of operating cycles	-	CYCLE	5.574	10.141	5.386	9.868
Height	km	H	0.090	0.867	0.047	0.231
Mach	-	MA	0.086	0.343	0.002	0.493
Throttle resolver angle	°	TRA	0.000	0.141	0.000	0.146
Total temperature at fan inlet	°R	T2	0.000	0.193	0.000	0.224
Total temperature at LPC outlet	°R	T24	1.625	0.710	0.800	1.855
Total temperature at HPC outlet	°R	T30	1.134	1.525	1.191	0.301
Total temperature at LPT outlet	°R	T50	1.285	0.866	0.891	1.297
Pressure at fan inlet	psia	P2	0.000	0.146	0.000	0.118
Total pressure in bypass-duct	psia	P15	0.017	0.291	1.973	0.655
Total pressure at HPC outlet	psia	P30	1.334	0.698	0.966	0.812
Physical fan speed	rpm	NF	1.649	0.766	0.825	1.085
Physical core speed	rpm	NC	1.137	0.645	1.011	0.928
Engine pressure ratio (P50/P2)	-	EPR	0.000	0.149	0.264	0.518
Static pressure at HPC outlet	psia	PS30	0.971	1.569	1.004	1.026
Ratio of fuel flow to Ps30	pps/psi	PHI	1.357	0.497	1.586	0.886
Corrected fan speed	rpm	NRF	1.261	1.378	0.813	2.271
Corrected core speed	rpm	NRC	0.799	0.959	1.290	1.124
Bypass ratio	-	BPR	1.371	1.084	1.158	1.835
Burner fuel-air ratio	-	FARB	0.000	0.677	0.000	0.355
Bleed enthalpy	-	HT_BLEED	0.829	0.925	0.317	0.594
Demanded fan speed	rpm	NF_DMD	0.000	0.259	0.000	0.180
Demanded corrected fan speed	rpm	PCNFR_DMD	0.000	0.025	0.000	0.101
HPT coolant bleed	lbm/s	W31	0.752	0.554	1.652	0.382
LPT coolant bleed	lbm/s	W32	1.636	0.816	1.548	1.009

, which can be used to train and test fault diagnosis and prediction models. This dataset is widely used in the fields of machine learning and data mining, which provides valuable data support for aero-engine health management.

3.3. Feature Selection

In the context of predicting the health of turbofan engines, there are several redundant features. Appropriate feature selection methods can enhance model performance, reduce computational costs, and provide an improved model interpretability.

This study employs a feature selection method based on a random forest. In the random forest model, each decision tree is built by conducting bootstrap sampling from the original dataset [39]. The corresponding out-of-bag data are used to evaluate the prediction error of the $j$-th tree in the random forest. In order to assess the importance of the $i$-th feature, the values of $i$-th feature are randomly permuted and the out-of-bag data prediction error is reevaluated. The importance of the $i$-th feature is determined by calculating the ratio of the average and standard deviation of out-of-bag data errors for all trees in the entire random forest, both before and after disrupting the $i$-th feature. Equation (6) is used to calculate the standard deviation ${\sigma }_i$, and Equation (7) is used to calculate the feature importance $C_i$.

$${\sigma }_i=\sqrt{\frac{{\sum }_{j=1}^{n_t}{\left(er{r}_{i,j}-er{r}_{0,j}\right)}^2}{n_t}}$$

(6)

$$C_i=\frac{{\sum }_{j=1}^{n_t}\left(er{r}_{i,j}-er{r}_{0,j}\right)}{n_t·{\sigma }_i}$$

(7)

Here, ${\sigma }_i$ represents the standard deviation of the errors, $er{r}_{i,j}$ represents the prediction error after disrupting the $i$-th feature, $er{r}_{0,j}$ represents the prediction error before disrupting the $i$-th feature, $n_t$ is the number of trees, and $C_i$ is the importance of the $i$-th feature. In this section, we present the feature importance analysis for FD001, FD002, FD003, and FD004. The results are summarized in Table 2.

Table 2 shows that the FD001 subset has seven redundant features, and the FD003 subset has six redundant features. Redundant features are those with ‘0’ feature importance in Table 2. Upon examining of the raw data, it is evident that the values of these features remain constant during the operating cycles of turbofan engines. To reduce computational costs, these features are excluded from both the comparative and proposed models in this study.

3.4. Data Normalization and Correlation Analysis

When predicting the health index of turbofan engines, there is a significant difference in the order of magnitude between operational parameters and sensor parameters. To train neural networks, each feature needs to be weighted and summed. Without normalizing the original data, some features may carry excessively large weights, while others may have excessively small weights. This can complicate the training of neural networks and affect the predictive performance of the model.

Equation (8) expresses the Min–Max scaling calculation used to normalize both the operational parameters and sensor parameters in the original dataset. The goal is to mitigate the dimensional impact among different features.

$$X_{i-norm}=\frac{X_i-X_{min}}{X_{max}-X_{min}}$$

(8)

Here, $X_{min}$ represents the minimum of the feature, $X_{max}$ represents the maximum of the feature, $X_i$ is the original value, and $X_{i-norm}$ is the normalized value.

The Pearson linear correlation coefficient is the most commonly used linear correlation coefficient [40]. We calculated the Pearson correlation coefficient in order to analyze the correlation between features. Equation (9) shows the calculation of the Pearson correlation coefficient.

$$\mathrm{r}=\frac{{\sum }_{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(X_i-\overline{X}\right)}^2{\sum }_{j=1}^n{\left(Y_j-\overline{Y}\right)}^2}}$$

(9)

Here, $n$ is the sample size of the feature, $X_i$ represents the $i$-th sample value of the feature, $\overline{X}$ represents the mean value of the feature, $Y_i$ is the $i$-th value of the number of operating cycles, and $\overline{Y}$ is the mean value of the number of operating cycles.

For the No. 100 engine in the training set of the FD001, we selected the 15 features with feature importance larger than 0.1 (${\mathrm{C}}_{\mathrm{i}}>0.1$) from its records. We then normalized the original data using Min–Max normalization to obtain the corresponding normalized values. The results of the normalization are presented in Figure 5.

From Figure 5, it seems that 10 out of the 15 features exhibit an increasing trend with the number of engine cycles, while the remaining 5 features show a decreasing trend. Analyzing these changes in degradation can enhance our understanding of the relationship between the engine’s operational state and the features, providing valuable validation for feature selection result. The Pearson correlation coefficients between the number of operating cycles and the other 14 features are presented in

Table 3. Pearson correlation coefficients between number of operating cycles and all the 25 features of the NO. 100 engine in the training set of FD001. The underlined values present the redundant features.

Feature	Pearson Correlation Coefficient	Feature	Pearson Correlation Coefficient
CYCLE	1	EPR	0
H	−0.037	PS30	0.840
MA	0.013	PHI	−0.824
TRA	0	NRF	0.789
T2	0	NRC	−0.057
T24	0.692	BPR	0.739
T30	0.672	FARB	0
T50	0.809	HT_BLEED	0.682
P2	0	NF_DMD	0
P15	0.121	PCNFR_DMD	0
P30	−0.787	W31	−0.743
NF	0.780	W32	−0.714
NC	0.459

. All the 10 features in Figure 5a show positive Pearson correlation coefficients in Table 3, while the 5 features in Figure 5b present negative values in Table 3. In this way, the values in Table 3 justify the observed result of normalization.

3.5. BP Neural Network

The BP neural network model’s performance is significantly impacted by the number of hidden layers, number of neurons in the hidden layer, initial weights, and biases. The weights indicate the strength of the connections between neurons, while the biases are additional learnable parameters for neurons that introduce shifts and offsets to the neuron inputs, adjusting the input to the activation function [14].

This study consistently used a single hidden layer for the BP neural networks to control the variable. The number of neurons in the hidden layer was determined using the empirical formula, as shown in Equation (10) [15].

$$h_n=\sqrt{i_n+o_n}+\alpha$$

(10)

In the equation, $h_n$ represents the number of neurons in the neural network’s hidden layer, $i_n$ represents the number of neurons in the input layer, $o_n$ represents the number of neurons in the output layer, and $\alpha$ is a positive integer with a value ranging from 1 to 10.

The optimal value of $\alpha$ for finding the number of neurons in the hidden layer was found using grid search. The initial weights and biases of the BP neural network were optimized using ESO. The optimization problem for the weights and biases was formulated in a high-dimensional space, and the fitness function was based on the root mean square error of the BP neural network on training sets [13]. Figure 6 illustrates the optimization process. Equations (11) and (12) are used to calculate the dimensions and fitness values [14].

$$d=\left(i_n+1\right)\ast h_n+\left(h_n+1\right)\ast o_n$$

(11)

$$f_i=\frac{{\sum }_{i=1}^M{\left(t_i-{\widehat{t}}_i\right)}^2}{M}$$

(12)

Here, $d$ signifies the search space dimension of ESO, $f_i$ is the fitness value of the $i$-th iteration of ESO, $M$ indicates the number of samples in the training set, $t_i$ denotes the actual value of the $i$-th sample, and ${\widehat{t}}_i$ represents the model’s predicted value for the $i$-th sample.

3.6. Parameters of ESO

After determining the number of neurons in the hidden layer using grid search, the dimensionality of the ESO search space is defined based on the model structure. Subsequently, ESO optimizes the initial values of the weights and biases in the BP neural network by establishing the upper and lower bounds of the solution space and specifying the population size [13].

The population size is a crucial parameter in intelligent optimization algorithms. A large population size may increase the computation time required to execute the algorithm, potentially resulting in wasted computational resources. Conversely, a small population size may decrease population diversity, hindering the effective exploration of the solution space and potentially leading to a slower convergence. In addition, an excessively small population may become trapped in local optima, which can cause the global optimum to be missed. Therefore, the parameters of ESO are configured based on the values presented in

Table 4. Parameter values of ESO.

Parameter Name	Value
Limit of the solution space	[−2, 2]
Population size	10
Iterations	100

.

3.7. ESO Optimization Process

The overall framework of the proposed ESO-BP prediction model is illustrated in Figure 7. The implementation of ESO optimization involves the following specific steps:

1.

Employ grid search to determine the optimal number of neurons in the hidden layer and establish the optimal structure for the BP neural network model.

2.

Initialize the parameters for ESO, define the dimensions and value ranges of the search space based on the optimal model structure, and specify the population size, maximum iteration count, and fitness function using Equations (10)–(12).

3.

Randomly generate the initial population, compute the fitness value for each individual, and identify the elite sets for both the male and female populations.

4.

In the early iterations, the populations enter the exploration phase due to the quantity of food being lower than 0.25. In this phase, individuals in both the male and female populations randomly select an elite guider and move to positions near it. The exploration phase is presented in Equation (13) [16].

$$X_i\left(t+1\right)=X_{elite}\left(t\right)\pm 0.05\times \exp \left(\frac{-f_{elite}}{f_i}\right)\times \left(\left(l_u-l_l\right)\times rand+l_l\right)$$

(13)

Here, $X_i$ is the $i$-th individual position and $t$ is the current iteration count. $X_{elite}$ is the position of a randomly selected elite guider in an elite set. $f_{elite}$ is the fitness of the elite guider. $f_i$ is the fitness of the $i$-th individual. $l_u$ and $l_l$ are the upper and lower bounds of the problem, respectively.

5.

In the later stages of iteration, the populations enter the exploitation phase when the food quantity exceeds 0.25. In this phase, the behavior of the populations is influenced by the temperature. If the temperature is greater than the threshold, both male and female individuals move toward the population’s optimal position (food). If the temperature is lower than the threshold, there is a 40% probability that male and female populations engage in conflict, and a 60% probability that they engage in mating. The snake fight stage is defined in Equations (14) and (15). The snake mating stage is defined in Equations (16) and (17). The calculation for updating the worst individual is defined in Equation (18).

$$X_{i,m}\left(t+1\right)=X_{i,m}\left(t\right)\pm 2\times \exp \left(\frac{-f_{best,f}}{f_{i,m}}\right)\times rand\times \left(Q\times X_{best,f}-X_{i,m}\left(t\right)\right)$$

(14)

$$X_{i,f}\left(t+1\right)=X_{i,f}\left(t\right)\pm 2\times \exp \left(\frac{-f_{best,m}}{f_{i,f}}\right)\times rand\times \left(Q\times X_{best,m}-X_{i,f}\left(t\right)\right)$$

(15)

$$X_{i,m}\left(t+1\right)=X_{i,m}\left(t\right)\pm 2\times \exp \left(\frac{-f_{i,j}}{f_{i,m}}\right)\times rand\times \left(Q\times X_{i,f}\left(t\right)-X_{i,m}\left(t\right)\right)$$

(16)

$$X_{i,f}\left(t+1\right)=X_{i,f}\left(t\right)\pm 2\times \exp \left(\frac{-f_{i,m}}{f_{i,f}}\right)\times rand\times \left(Q\times X_{i,m}\left(t\right)-X_{i,f}\left(t\right)\right)$$

(17)

$$X_{worst}=l_l+rand\times \left(l_u-l_l\right)$$

(18)

In Equations (14)–(18), $X_{i,m}$ and $X_{i,f}$ are the positions of the $i$-th individual in the male and female populations, $X_{best,m}$ and $X_{best,f}$ are the positions of the optimal individuals in the male and female populations, $f_{best,m}$ and $f_{best,f}$ are the fitness values of the optimal individuals in the male and female populations, $f_{i,m}$ and $f_{i,f}$ are the fitness values of current male or female individuals, and $Q$ is the food quantity defined in Equation (2) [16].

6.

Generate reverse populations using Equation (4). If there is a new individual with better fitness, the ESO algorithm will replace the original agent with the new individual.

7.

Continuously adjust the position of the snake agent until reaching the maximum iteration count, obtaining the optimal spatial position.

8.

Based on the coordinates of the optimal solution, set the weights from the input layer to the hidden layer, the biases for the hidden layer, the weights from the hidden layer to the output layer, and the biases for the output layer for the BP neural network.

The pseudo-code of ESO is defined as Algorithm 1:

Algorithm 1 The pseudo-code of ESO

Stage 1. Initialization Initialized Problem Setting (d, $l_l$, $l_u$, n, T)

Generate the ESO population $S_n\left(n=1,2,\dots ,N\right)$.

Calculate the fitness ($f_i$) of each individual of population.

Divide the ESO population $S_n$ into two equal groups. Generate the Elite-Guided set.

Stage 2. ESO Iteration   for $i=1:T$     Calculate Temp using Equation (1).     Calculate Q using Equation (2).     if $Q>1$ then       $Q=1$     end if     Calculate Elite 4 using Equation (3).     Select Elite 1, Elite2, and Elite 3 according to the fitness.     if $Q<0.25$ then      Perform exploration using Equation (13).     else if $Temp>0.6$ then      Perform exploitation, all individuals move to food.     else      if $rand>0.6$ then      Male population fight with female population using Equations (14) and (15).      else       Snakes mate using Equations (16) and (17).       Update the worst individual using Equation (18).      end if     end if     Generate reverse populations using Equation (4).     Calculate fitness of reverse individuals $f_{i-r}$.     if $f_{i-r}<f_i$      Change the best solution to reverse individual.     end if    end for Return optimal solution.

4. Simulation Experiments and Result Analysis

4.1. Description to Hyperparameters

To achieve a better forecasting performance, it is essential to systematically optimize the relevant hyperparameters of the models. These hyperparameters include the number of hidden layer nodes, learning rate, and the size of the training batch. To ensure the control of the variables in later comparison experiments among different sub-datasets and deep learning models, we used empirical values of 1 for the number of hidden layers, 1000 for the max iteration epochs, 128 for the batch size, and 0.01 for the initial learning rate of BP neural networks. To achieve a higher accuracy, the number of hidden nodes can be tuned sequentially in the range given by Equation (11), while keeping the other parameters unchanged. Figure 8 presents the grid search result of the number of hidden nodes. According to Figure 8, the number of nodes in the hidden layer of BP neural networks was taken as 13.

4.2. Model Evaluation Metrics

To evaluate the accuracy of the prediction model proposed in this study, three evaluation metrics were used. The root mean square error was used to measure the difference between the predicted and actual values. This metric indicates the average deviation between predicted and actual values and is commonly used in regression tasks. A smaller RMSE value indicates a greater accuracy, and its unit is the same as that of the original data. The $RMSE$ is calculated using Equation (19).

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left({\widehat{y}}_i-y_i\right)}^2}{n}}$$

(19)

The coefficient of determination, which represents the model’s explanatory power for the dependent variable, has a range of values from 0 to 1. A higher value indicates a better fit of the model to the data. The $R^2$ is calculated using Equation (20).

$$R^2=1-\frac{{\sum }_{i=1}^n{\left({\widehat{y}}_i-y_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(20)

It is better for the predicted health index of the turbofan engine to be slightly underestimated rather than overestimated [35]. If the predicted health index cannot reach zero before a failure occurs, it poses a serious safety hazard to the aircraft. To address this issue, this study employs a scoring function ($Score$) proposed by the International Conference on Prognostics and Health Management (PHM08) Data Challenge [38]. The scoring function penalizes situations where the model’s predicted values are larger more severely. The score is calculated by Equations (21) and (22).

$$S_i=\left\{\begin{array}{c}{e}^{-\frac{{\widehat{y}}_i-y_i}{13}}-1,\hspace{1em}\hspace{1em}{\widehat{y}}_i<y_i\\ e^{\frac{{\widehat{y}}_i-y_i}{10}}-1,\hspace{1em}\hspace{1em}{\widehat{y}}_i\ge y_i\end{array}\right.$$

(21)

$$Score=\sum _{i=1}^n{S}_i$$

(22)

From Equation (19) to Equation (22), ${\widehat{y}}_i$ indicates the predicted value of the $i$-th sample, while $y_i$ denotes its real value and $n$ is the total number of samples.

4.3. Optimization Results

To evaluate the optimization ability of various intelligent optimization algorithms on BP neural networks, we optimized the initial weights and biases of the networks using the eight comparison algorithms mentioned above. The population size, number of iterations, and limits of the solution space for all eight algorithms were controlled to be the same as ESO above. The fitness function of the optimization problem was selected as the RMSE of the BP neural network on the training set. The convergence curves of the nine algorithms are depicted in Figure 9.

In Figure 9, all nine algorithms improved the accuracy of BP neural networks. However, BWO and SSA performed as the worst and the second worst, respectively, in the simulation test of the benchmark functions F18 and F30. Correspondingly, BWO performed the worst in optimizing the initial weights and biases in FD003 and FD004, while SSA performed the second worst in FD001 and FD003. Furthermore, ESO demonstrated a faster and more precise convergence than SO in all four sub-datasets. The performance indicators in simulation and practice were consistent with each other. These findings lead to three conclusions: 1. Using intelligent optimization to initialize the weights and biases of the BP neural network is an effective method for improving model prediction accuracy; 2. The more suitable the weights and biases found by initialization methods, the more accurate the prediction performed by BP neural networks; and 3. ESO, enhanced by the elite-guided strategy and reverse learning mechanism, has both theoretical and practical performance improvements.

Figure 9 shows the convergence curves of DE, PSO, GWO, SSA, BWO, DBO, CSO, SO, and ESO in optimizing the fitting of BP neural networks to four different sub-datasets. After 100 iterations, ESO achieved a 7.356% reduction in the RMSE of the BP neural network on the FD001 sub-dataset. This result demonstrates that ESO can enhance the prediction accuracy of the BP neural network for the RHI prediction problem in turbofan engines within a relatively small number of iterations. Compared to the other eight mainstream optimization algorithms, ESO converged to smaller fitness values. This indicates that ESO is more suitable for optimizing the BP neural network on turbofan engine RHI prediction.

4.4. Prediction Results

The ESO-BP model was used to predict the health index of engines in subsets FD001, FD002, FD003, and FD004 of C-MAPSS. In FD001, the predicted and actual health indexes of the last operating cycle for all engines in the test set were compared. Figure 10 shows slight errors in the RHI prediction of each engine’s last cycle. This indicates that the proposed method can accurately predict the RHI of turbofan engines. Additionally, four engines were randomly selected from the test set to compare the predicted and actual values of the health degradation curves. Figure 11 shows the fitting ability of the degradation trend for four engines in FD001. The results indicate that ESO-BP demonstrates a high accuracy in multivariate input and univariate output prediction.

4.5. Comparative Experiments

To evaluate the prediction accuracy and stability of ESO-BP, we compared its performance with that of convolutional neural network (CNN), recurrent neural network (RNN), long short-term memory (LSTM), and BP neural networks trained using two different functions. The CNN, RNN, and LSTM models employed in this study used the Adam optimizer [41] for training. BP(SCG) and ESO-BP(SCG) were trained using the scaled conjugate gradient (SCG) backpropagation algorithm [42]. BP(LM) and ESO-BP(LM) were trained using the Levenberg–Marquardt (LM) backpropagation algorithm [43].

Table 5. Comparison of different prediction models for 100 trials. The bold number represents the best model.

Dataset	Model	RMSE		R2		Score
		Avg	Std	Avg	Std	Avg	Std
FD001	CNN	6.641 × 10−2	4.346 × 10−3	9.163 × 10−1	1.158 × 10−2	5.728 × 10+1	4.813 × 100
	RNN	6.773 × 10−2	2.841 × 10−3	9.132 × 10−1	7.560 × 10−3	5.764 × 10+1	3.445 × 100
	LSTM	6.253 × 10−2	1.722 × 10−3	9.261 × 10−1	4.219 × 10−3	5.152 × 10+1	2.342 × 100
	BP(SCG)	6.337 × 10−2	1.346 × 10−3	9.241 × 10−1	3.249 × 10−3	5.175 × 10+1	1.812 × 100
	ESO-BP(SCG)	6.223 × 10−2	7.070 × 10−4	9.268 × 10−1	1.664 × 10−3	5.021 × 10+1	8.930 × 10−1
	BP(LM)	6.070 × 10−2	4.840 × 10−4	9.304 × 10−1	1.115 × 10−3	4.891 × 10+1	5.164 × 10−1
	ESO-BP(LM)	6.040 × 10−2	4.160 × 10−4	9.311 × 10−1	9.580 × 10−4	4.859 × 10+1	4.111 × 10−1
FD002	CNN	8.233 × 10−2	6.792 × 10−3	8.839 × 10−1	2.008 × 10−2	1.854 × 10+2	1.569 × 10+1
	RNN	8.645 × 10−2	7.401 × 10−3	8.719 × 10−1	2.294 × 10−2	1.976 × 10+2	1.629 × 10+1
	LSTM	7.543 × 10−2	6.836 × 10−3	9.024 × 10−1	1.938 × 10−2	1.675 × 10+2	1.461 × 10+1
	BP(SCG)	7.427 × 10−2	4.604 × 10−3	9.058 × 10−1	1.183 × 10−2	1.656 × 10+2	1.389 × 10+1
	ESO-BP(SCG)	7.195 × 10−2	3.950 × 10−3	9.116 × 10−1	9.829 × 10−3	1.583 × 10+2	1.283 × 10+1
	BP(LM)	6.528 × 10−2	2.730 × 10−4	9.275 × 10−1	6.080 × 10−4	1.380 × 10+2	7.994 × 10−1
	ESO-BP(LM)	6.515 × 10−2	1.880 × 10−4	9.278 × 10−1	4.160 × 10−4	1.374 × 10+2	4.245 × 10−1
FD003	CNN	7.155 × 10−2	4.230 × 10−3	9.098 × 10−1	1.164 × 10−2	7.954 × 10+1	6.028 × 100
	RNN	7.795 × 10−2	3.068 × 10−3	8.932 × 10−1	8.565 × 10−3	8.724 × 10+1	3.304 × 100
	LSTM	7.133 × 10−2	1.807 × 10−3	9.106 × 10−1	4.604 × 10−3	7.732 × 10+1	1.962 × 100
	BP(SCG)	6.892 × 10−2	1.986 × 10−3	9.165 × 10−1	4.976 × 10−3	7.478 × 10+1	3.088 × 100
	ESO-BP(SCG)	6.729 × 10−2	7.630 × 10−4	9.205 × 10−1	1.812 × 10−3	7.269 × 10+1	1.384 × 100
	BP(LM)	6.779 × 10−2	6.860 × 10−4	9.193 × 10−1	1.640 × 10−3	7.187 × 10+1	7.749 × 10−1
	ESO-BP(LM)	6.744 × 10−2	4.090 × 10−4	9.202 × 10−1	9.690 × 10−4	7.146 × 10+1	5.216 × 10−1
FD004	CNN	9.433 × 10−2	5.466 × 10−3	8.447 × 10−1	1.851 × 10−2	2.600 × 10+2	1.525 × 10+1
	RNN	9.982 × 10−2	8.187 × 10−3	8.256 × 10−1	2.961 × 10−2	2.756 × 10+2	1.847 × 10+1
	LSTM	9.214 × 10−2	4.700 × 10−3	8.520 × 10−1	1.556 × 10−2	2.518 × 10+2	1.156 × 10+1
	BP(SCG)	8.889 × 10−2	3.833 × 10−3	8.623 × 10−1	1.210 × 10−2	2.462 × 10+2	1.314 × 10+1
	ESO-BP(SCG)	8.615 × 10−2	2.283 × 10−3	8.708 × 10−1	6.897 × 10−3	2.365 × 10+2	7.276 × 100
	BP(LM)	7.546 × 10−2	6.580 × 10−4	9.010 × 10−1	1.744 × 10−3	1.973 × 10+2	1.875 × 100
	ESO-BP(LM)	7.529 × 10−2	4.020 × 10−4	9.014 × 10−1	1.054 × 10−3	1.970 × 10+2	1.345 × 100

displays the averages (Avg) and standard deviations (Std) of the evaluation metrics after 100 model training turns. The results indicate that ESO-BP consistently outperformed mainstream models in all the evaluation metrics, proving the effectiveness of the ESO-BP neural network model in predicting the RHI of turbofan engines.

The results demonstrate that ESO-BP(LM) attained the highest prediction stability, showcasing the smallest standard deviation across all evaluation metrics. This suggests that combining the LM algorithm with improved initial values assists neural networks in avoiding local optima, leading to increased stability. In terms of prediction accuracy, ESO-BP(LM) outperformed in FD001, FD002, and FD004, showcasing a superior predictive capability across various operational conditions. The R² and RMSE of ESO-BP(SCG) on FD003 were slightly better than those of ESO-BP(LM), indicating that predicting the RHI of turbofan engines under multi-fault conditions poses a particular challenge. ESO improved the performance of BP neural networks with both SCG and LM, indicating that using ESO to initialize BP neural networks is a beneficial approach. The BP neural networks in the following study were trained with the LM algorithm for a better performance, according to the results in Table 5.

To illustrate the deviation between the predicted and real RHI values for each operating cycle of all test set engines in FD001, the real RHI value serves as the horizontal axis, and the predicted RHI value serves as the vertical axis in Figure 12. The deviation between predicted and real values for the different models is shown in Figure 12. The majority of the scatter points from the four models align closely with the reference line. This figure illustrates that the ESO-BP(LM) neural network model exhibits fewer anomalous scatter points and achieves the best prediction results. Figure 13 shows that the predicted RHI degradation curves of the turbofan engines using the proposed method are closer to the real degradation curves compared to other models.

4.6. Ablation Experiment

An ablation experiment was conducted to assess the impact of each enhancement approach on the BP model. The experiment compared the baseline BP(LM) with SO-BP(LM) and ESO-BP(LM).

Table 6. Averages of evaluation metrics for the ablation experiment. The bold number represents the best model.

Dataset	Model	RMSE	R2	Score
FD001	BP	6.070 × 10−2	9.304 × 10−1	4.891 × 10+1
	SO-BP	6.065 × 10−2	9.305 × 10−1	4.859 × 10+1
	ESO-BP	6.040 × 10−2	9.311 × 10−1	4.859 × 10+1
FD002	BP	6.528 × 10−2	9.275 × 10−1	1.380 × 10+2
	SO-BP	6.524 × 10−2	9.276 × 10−1	1.380 × 10+2
	ESO-BP	6.515 × 10−2	9.278 × 10−1	1.374 × 10+2
FD003	BP	6.779 × 10−2	9.193 × 10−1	7.187 × 10+1
	SO-BP	6.761 × 10−2	9.197 × 10−1	7.182 × 10+1
	ESO-BP	6.744 × 10−2	9.202 × 10−1	7.146 × 10+1
FD004	BP	7.546 × 10−2	9.010 × 10−1	1.973 × 10+2
	SO-BP	7.610 × 10−2	8.993 × 10−1	1.994 × 10+2
	ESO-BP	7.529 × 10−2	9.014 × 10−1	1.970 × 10+2

presents the average evaluation metrics of 100 trials for BP(LM), SO-BP(LM), and ESO-BP(LM).

The evaluation metrics derived from 100 trials demonstrate that SO optimized BP neural networks on the FD001, FD002, and FD003 sub-datasets of the C-MAPSS. ESO further enhanced the optimization effect and showed improvement on the FD004 subset, where SO-BP performed poorly. In order to visualize the improvements of SO and ESO on BP neural networks, Figure 14 shows the RMSEs of BP(LM), SO-BP(LM), and ESO-BP(LM).

5. Discussion

The objective of training a neural network is to fine-tune weights and biases to find a parameter set that aligns the predicted values with the actual values as closely as possible [44]. Typically, this process involves finding the global optimal solution, which requires determining the combination of parameters that minimizes the loss function across the entire parameter space. The exceptional prediction performance of the ESO-BP model is attributed to three primary factors: optimizing the initial weights and biases through ESO, employing the RHI metric for turbofan engine degradation, and selecting LM as the training function for BP models.

1.

Initializing neural network weights was demonstrated to accelerate convergence, prevent entrapment in local optima, simplify network architecture, identify feature importance, and improve prediction performance [45]. In this study, adjustments to the initial weights and biases were applied before training to improve the convergence accuracy and stability of the BP neural networks. ESO exhibited a superior convergence compared to SO in benchmark functions. Similarly, in the ablation experiment, ESO-BP outperformed SO-BP. This suggests that providing better initial weights and biases results in higher accuracy and stability in BP neural networks. ESO determined the initial weights and biases that minimized the RMSE after its iteration process, preventing the BP neural network from easily becoming stuck in local optima and experiencing stability issues due to random parameter selection.

2.

In contrast to RUL prediction, Jiang et al. [46] also developed a health index prediction model for turbofan engines in 2023. According to their study, the average accuracy of 10 trials for hybrid methods was significantly lower than that of single methods. Several methods were combined to create multiple new models, including empirical mode decomposition (EMD), variational mode decomposition (VMD), scale-adaptive attention mechanism (SAA), dynamic step size-based fruit fly optimization algorithm (DSSFOA), bidirectional long short-term memory network (BiLSTM), support vector regression (SVR), and LSTM. The average accuracies of the FD002 sub-dataset are compared in

Table 7. Comparison of RMSE for hybrid methods [46] and ESO-BP model.

Model	Number of Trial	RMSE
EMD-SAA-BiLSTM	10	0.101
VMD-SAA-BiLSTM	10	0.082
VMD-DSSFOA-SVR	10	0.156
VMD-DSSFOA-LSTM	10	0.093
VMD-DSSFOA-SAA-LSTM	10	0.086
ESO-BP	100	0.065

.

3.

The Levenberg–Marquardt algorithm is particularly effective in training BP neural networks for regression problems [47]. The LM algorithm provides numerical solutions for nonlinear minimization. It combines the strengths of the Gauss–Newton algorithm and the gradient descent method by dynamically adjusting parameters during execution to address the limitations of both methods. When the gradient drops rapidly, the LM algorithm behaves more like the Gauss–Newton algorithm. Otherwise, it behaves more like gradient descent [43].

6. Conclusions

This study improved the convergence performance of an existing optimization algorithm through an elite-guided strategy and a reverse learning mechanism. Simulation experiments were conducted using three benchmark test functions to compare the elite snake optimizer (ESO) with eight mainstream intelligent optimization algorithms. The results indicated that ESO exhibited a higher convergence accuracy and improved stability. Furthermore, ESO was used to initialize the weights and biases of the BP neural network, resulting in the creation of the ESO-BP model. The ESO-BP model was compared with traditional CNN, RNN, LSTM, and a baseline BP neural network in terms of prediction accuracy and stability using the C-MAPSS dataset. The findings showed that the ESO-BP model enhanced the prediction accuracy and stability of BP across all the sub-datasets of C-MAPSS. In conclusion, the ESO-BP model has theoretical significance in improving intelligent optimization algorithms and practical value in predicting the health index of turbofan engines.

In the future, intelligent optimization algorithms can be employed to initialize the parameters for other neural networks and applied across various domains in machine learning. The primary concern is the unacceptable time cost required for intelligent optimization algorithms to search for the global optimum in high-dimensional spaces, attributed to the large number of weights and biases. Accordingly, we will develop effective pruning methods to eliminate specific dimensions of the connection weights, reducing time costs and streamlining the networks. This is a promising research direction.