#### Abstract

In order to improve the accuracy and calculation efficiency of aeroengine rotor vibration reliability analysis, a time-varying rotor vibration reliability analysis method under the aeroengine operating state is proposed. Aiming at the highly nonlinear and strong coupling of factors affecting the reliability of aeroengine rotor vibration, an intelligent neural network modeling framework (short form-INNMF) is proposed. The proposed method is based on DEA, with QAR information as the analysis data, and four factors including engine working state, fuel/oil working state, aircraft flight state, and external conditions are considered to analyse the rotor vibration reliability. INNMF is based on the artificial neural network (ANN) algorithm through improved particle swarm optimization (PSO) algorithm and Bayesian Regularization (BR) optimization. Through the analysis of the rotor vibration reliability of the B737-800 aircraft during a flight mission from Beijing to Urumqi, the time-varying rotor vibration reliability was obtained, which verified the effectiveness and feasibility of the method. The comparison of INNMF, random forest (RF), and ANN shows that INNMF improves analysis accuracy and calculation efficiency. The proposed method and framework can provide useful references for aeroengine rotor vibration analysis, special treatment, maintenance, and design.

#### 1. Introduction

As the core of the aircraft, aeroengines operate under extreme conditions such as high speeds, high temperatures, and variable loads. Their reliability is the most important factor supporting the performance and survivability of the aircraft [1–3]. In recent years, the rapid development of the world aviation industry has put forward new and higher requirements for the reliability and sustainable airworthiness of the aeroengine. Countries with developed aviation industries have integrated engine performance condition monitoring systems and fault diagnosis systems. Among them, aeroengine rotor vibration value prediction and vibration trend analysis are considered to be one of the important concerns. Therefore, in order to avoid safety accidents caused by aeroengine vibration failures, it is necessary for researchers to explore methods for evaluating the reliability of aeroengine rotor vibration during operation. The evaluation results can provide the basis for scientifically formulating engine maintenance plans, extending the service life of engines, and ensuring aircraft flight safety.

Rotation and transmission components are the key components of engine energy conversion and power output, as well as important components that affect engine safety [4]. Rotor vibration always exists when the aeroengine is working, and abnormal vibration is usually considered as a precursor to failure. On the contrary, vibration has been proven to be one of the most reliable and sensitive technical methods for fault diagnosis of rotating and transmission components. The relative maturity of vibration testing methods makes fault diagnosis methods based on vibration analysis effective for popularization and application [5]. Some scholars have studied the vibration characteristics of the aeroengine through experiments and numerical methods, such as Minh Hai and Bonello [6, 7] proposed a method suitable for solving the dynamic response of the high-dimensional complex rotor system in the time and frequency domain, established the complex finite element model of double rotor and three rotor engines, and studied the vibration response characteristics of the aeroengine under multifrequency excitation; to monitor the rolling bearing operating status with casings in real time efficiently and accurately, Yan-Ting [8] proposed a fusion method based on *n*-dimensional characteristic parameters’ distance (*n*-DCPD) for rolling bearing fault diagnosis with two types of signals including vibration signal and acoustic emission signals; Pennacchi et al. [9] proposed a self-adaptive rotor system imbalance recognition method, which can effectively avoid the problem of weight coefficient selection when using the least square method to recognize imbalance; Williams [10] carried out the modeling of the blade-casing rubbing fault, established a detailed model of the wear of the inner bushing of the casing in the model, and verified the applicability of the model; Li et al. [11] established the dynamic model of the drum-disk-shaft structure and studied the effects of support stiffness, coupling stiffness, and rotating speed on vibration behaviours of the system on the basis of considering the elastic support, flexible connection, and the Coriolis and centrifugal effects on the drum due to the rotation. However, in the practical engineering application, the actual state of each engine is very different, and the performance degradation rate is also different in the use process, so the prediction accuracy of the model may not reach the ideal level. In addition, the existing aeroengine vibration assessment work is not sufficient in the comprehensive evaluation of the operating conditions, operation status, and rotor vibration reliability of the research object. With the development of the aviation industry, the traditional methods cannot meet the needs of the industry, so the ability of rotor vibration reliability analysis during the operation of aeroengines needs to be improved urgently.

The International Civil Aviation Organization (ICAO) emphasized in the safety management system (SMS) that the identification of hazards should be positive and forward-looking, and it is necessary to actively identify hazards that have not yet occurred [12]. Federal Aviation Administration (FAA) proposed the safety “monitoring/data analysis” (MSAD) process [13]. The development of advanced sensor technology, data acquisition, and transmission technology has created conditions for the reliability evaluation of rotor vibration during aeroengine operation based on operating data [14, 15]. It is generally believed that the aircraft’s quick access recorder (QAR) data is used for daily operation monitoring [16, 17]. QAR records a large amount of data, a complete range of parameters and high recording frequency. In addition, the wireless QAR transmission technology developed in recent years has been gradually applied. The analysis and application of QAR data has been paid more and more attention by researchers. Operational vibration reliability analysis not only considers the characteristics of flight missions and the operating characteristics of aeroengines but also the impact of the equipment’s own health, external environment, system operating conditions, and system operating behaviours. In this study, based on the extracted QAR data, combined with engine operating state, fuel/oil operating state, aircraft flight state, and aircraft operation environment, the rotor vibration reliability of the aeroengine during operation is analysed.

QAR data collected during aeroengine operation are characterized by many parameters, large quantity, time-varying, strong coupling, and nonlinearity, which lead to low analysis accuracy and long calculation time [18, 19]. The artificial neural network (ANN) algorithm is considered as an intelligent learning method, which has strong nonlinear mapping ability and strong robustness. In view of the advantages of high computational efficiency and high accuracy, ANN has been widely used in data mining and pattern recognition [20, 21]. In the field of aviation, Dong applies the deep neural network to aircraft parameter identification to detect and characterize aircraft icing [22]; Omar Alkhamisi and Mehmood used the integration of the machine learning algorithm and deep learning algorithm to improve risk prediction in the aviation system [23]; Zhang and Mahadevan trained two different types of deep learning models from different angles to predict flight path and proposed a risk prediction method based on the deep learning long-term short-term memory structure recurrent neural network [24]. It can be seen that ANN has strong nonlinear mapping ability and simplification ability and can learn historical data to quantitatively predict the trend of selected parameters. However, in the process of fitting time-varying and highly nonlinear functions, over fitting and local optimization problems often occur in the training process, which affects the prediction accuracy and limits its further application in operational vibration reliability analysis.

The purpose of this study is to propose a time-varying rotor vibration reliability analysis method based on Data Envelopment Analysis (DEA), which considers four factors: aeroengine operating status, fuel/oil operating status, aircraft operating status, and aircraft operating environment. An intelligent neural network modeling framework (INNMF) is developed to analyse the relationship between features and predict the vibration reliability. Hereinto, the improved particle swarm optimization (PSO) algorithm is used to search the initial weights and thresholds of ANN, and the final weights and thresholds are trained based on the training performance function of the Bayesian Regularization (BR) algorithm. The feasibility and effectiveness of the proposed method and framework are verified by the vibration safety analysis of the aeroengine running in a specific aircraft mission.

#### 2. Reliability Analysis Method of Rotor Vibration

##### 2.1. Feature Selection

The reliability analysis of aeroengine rotor vibration during operation needs to consider four kinds of factors: aeroengine operating status, fuel/oil operating status, aircraft operating status, and aircraft operating environment. 21 specific characteristics of the four factors are selected, as shown in Table 1.

##### 2.2. Evaluating Indicator

Combined with the time-space relationship of aircraft operation, the DEA method is used to analyse the operating status characteristics and vibration reliability of the aeroengine. As a classical nonparametric model, DEA was proposed by A. Charnes and W. W. Cooper in 1978. Based on the concept of “relative efficiency,” it is a systematic analysis method to evaluate the performance of decision-making units (DMU) according to the data of multi-index input and output. The advantage of DEA is that it does not need to know the correlation between goals in advance and can achieve the objective evaluation of the performance level, eliminating the error caused by human factors [25].

Assuming that **X**_{j} and **Y**_{j} are the input vector and output vector of the DMU, respectively, use (**X**, **Y**) to represent the entire activity of this DMU. The set **T** formed by the activities of *m* types of inputs and *s* types of outputs of *n* decision-making units is called the production possible set, which is expressed as follows:

Assign appropriate weights to each input and output, and assume that the input and output weight vectors are as follows:

The efficiency index refers to the ratio of output to input when the input is and the output is *u*^{T}*y*_{j} under the weight coefficient **v** and **u**. The efficiency index of the decision-making unit DMU_{j} (1 ≤ *j* ≤ *n*) can be expressed as

The efficiency index is used as a constraint to constitute a fractional programming problem as follows:

Through Charnes–Cooper transformation, it can be expressed as

Convert the fractional programming model into a linear programming model, as shown in the following equation:

The input and output data of *n* decision-making units to be evaluated are substituted into equation (6) to get the performance evaluation index of each evaluation unit, and the result is **B**_{i} (*i* = 1, 2, …, *n*):where **R** is used as the reliability evaluation index of aeroengine rotor vibration. Among them, **R** is obtained by min-max standardization of each performance evaluation index and **C** (as shown in equation (8)).

#### 3. Basic Theory of INNMF

##### 3.1. ANN

Backpropagation artificial neural network (BP-ANN) can effectively fit the input vector **x**_{i} and output the response **y**_{i} (**x**_{i}), where **x**_{i} = [*x*_{1}, *x*_{2}, … *x*_{i}] (*i* = 1, 2, …, *n*). Because of its arbitrary shape and strong adaptive ability, it can accurately fit the complex functional relationship between random variables and response variables. Therefore, the output response can be obtained by constructing the BP-ANN model without solving a large number of dynamic equations, which greatly reduces the amount of calculation and improves the calculation speed. The topology model of BP-ANN is shown in Figure 1. The number of neurons in the input layer *n*_{i} and the number of neurons in the output layer *n*_{o} are determined by the number of input vectors and the number of output responses, respectively. The number of neurons in the hidden layer *n*_{h} can be obtained by the following equation:where *r* is the empirical number during [0, 5].

The random variable **x** (**x** ∈ **R**^{n}) and a group of dynamic responses **Y**_{i} (**x**_{i}) (**y** ∈ **R**) are fitted by the ANN regression function *f* (*x*). The function is defined aswhere **W** is the weight and threshold vector, and the problem of function fitting is transformed into the problem of finding the best weight and threshold. The training performance function of ANN is expressed aswhere *e*(·) is the training error function.

Assuming that the weight and threshold vectors after the *k*th iteration are **W**^{k}, the training error function *e* (**W**) can be estimated by Taylor series expansion of **W**^{k} position, such aswhere ** J**(·) is the Jacobian matrix and

*o*(·) is infinitesimal of the higher order.

Since **W** changes very little between the two iterations, **W**^{k+1} − **W**^{k} is a very small value. Therefore, equation (12) can be translated into

Then, the training function can be written aswhere Δ(**W**^{k}) is the variation of **W** between the *k*th and (*k* + 1)th iterations.

According to the three-layer BP-ANN model, the structure of the training performance function is defined aswhere **W**_{ij} is the connection weight between the input layer node *i* and hidden layer neuron *j*, **b**_{j} is the *j*th threshold of the hidden layer, **W**_{jk} is the connection weight between the hidden layer node *j* and output layer node *k*, **B**_{k} is the output layer threshold, *f*_{1}(·) is the transfer function of the hidden layer, *f*_{2}(·) is the transfer function of the output layer, *n* is the number of input layer nodes, and *m* is the number of hidden layer nodes.

##### 3.2. INNMF Modeling

For the training problem of the training performance function in aeroengine rotor vibration reliability analysis, high nonlinearity and strong coupling lead to high complexity of the function, which always makes it difficult to obtain accurate weights and thresholds in training. In order to solve the above problems, INNMF is developed to improve the training performance function from two aspects: (1) PSO algorithm is improved to find the initial optimal weight and threshold and avoid premature convergence. (2) BR algorithm is used to obtain the final optimal weight and threshold through network training, so as to achieve better generalization ability and acceptable prediction accuracy.

###### 3.2.1. Initial Weight and Threshold

The optimization of initial weight and threshold of ANN is very important to improve the prediction accuracy of ANN. As an important search algorithm based on particle swarm optimization, PSO has excellent search ability and low requirements for parameters, and the process of PSO is easy to realize [26, 27]. In this study, the PSO algorithm is used to optimize ANN to improve the accuracy. However, the fixed inertia weight and learning factor of the conventional PSO algorithm will fall into the blind search and obtain the local optimal solution, which greatly affects the search efficiency and accuracy of the global initial optimal weight and threshold in the solution space.

In order to solve the above problems, this section relies on the improved PSO algorithm to optimize ANN to improve the recognition accuracy. The dynamic inertia weight and dynamic learning factor with the number of iterations are proposed. The purpose is to complete the dynamic search of particle swarm optimization and balance the global search ability and local search ability, so as to obtain a better optimal solution set.

The basic idea of the improved PSO algorithm is as follows. (1) A group of particles are initialized in space, each particle is a potential solution, and the string of weights and thresholds of the neural network are used as the position of particles. (2) Taking the training error function as the fitness function, all particles follow the current optimal particle to search in the solution space, and update the individual position by tracking the individual extreme value and group extreme value. (3) The optimal particle is selected to update the individual extreme value and the position of the group extreme value until the optimal solution is found, that is, the initial optimal weight and threshold of the neural network.

The update of particle position and velocity is shown inwhere *i* is the *i*th particle, *k* is the current number of iterations, *K* is the maximum number of iterations, **V**_{id} is the current particle velocity, **X**_{id} is the current particle position, **P**_{id} is the current individual extremum, **P**_{gd} is the current population extremum, *r*_{1} and *r*_{2} are random numbers between [0, 1], *iw*(*k*) is the dynamic nonlinear inertia weight, *iw*_{0} is the initial inertia weight, *iw*_{K} is the inertia weight obtained under the maximum number of iterations, *c*_{1}(*k*) and *c*_{2}(*k*) are the dynamic nonlinear individual learning factor and the dynamic nonlinear global learning factor in the *k*th iteration, respectively, *c*_{1o} and *c*_{2o} are the initial individual learning factor and the initial population learning factor, respectively, and *c*_{1K} and *c*_{2K} are individual learning factors and group learning factors in the *k*th iteration, respectively. The inertia weight **w** reflects the degree that the particle’s current velocity inherits the previous velocity. Larger inertia weight is beneficial to global search, while smaller inertia weight is more beneficial to local search.

###### 3.2.2. Trained by BR

BR theory shows that smaller weights and thresholds have less over fitting and faster convergence and have good generalization performance for input variables outside the training set [28]. By adding prior conditions, the solution space is reduced and the possibility of finding wrong solutions is reduced. Studies have shown that smaller weight and threshold are beneficial to the smoothness and simplicity of the training performance function. On the one hand, the BR algorithm can reduce the training error by reducing the network weight to improve the training performance. On the other hand, the over fitting of BP-ANN is avoided, so the calculation accuracy of BP-ANN is improved. The training performance function **E** based on BR is introduced as follows:where *k*_{1} and *k*_{2} are proportional coefficients, **w**_{j} is the network weight, *ε* is the expected output error function, **W** is the weight threshold vector of each layer of the network, *K* is the number of iterations, **Z** is the Jacobian matrix of *ε*, and *λ* is the iteration variable.

##### 3.3. Strength of INNMF

Based on the above analysis, INNMF is expected to improve the efficiency and accuracy of operational safety analysis by integrating ANN and improving PSO and BR algorithm. The main reasons are as follows:(1)ANN has strong fitting ability and can deal with the nonlinear relationship between input random variables and output responses. The computational task of solving complex equations is reduced and the computational efficiency is improved.(2)The dynamic inertia weight and dynamic learning factor in the improved PSO algorithm can reduce the search cost and improve the search accuracy. It effectively avoids the problem of poor initial training value setting and easy to produce the local optimal solution in INNMF.(3)The dynamic performance function based on the BR algorithm reduces the training error and limits the complexity of the intelligent neural network model.(4)The proposed INNMF is expected to improve the calculation efficiency and accuracy of rotor vibration reliability analysis during operation.

##### 3.4. Analysis Process of Rotor Vibration Reliability

According to the characteristics of QAR data and the requirements of aeroengine rotor vibration reliability analysis, the key features are extracted. The DEA method is used to evaluate the vibration margin *R* of the aeroengine. The INNMF is used to fit and predict, and the contribution of each characteristic to the operation vibration margin *R* is analysed. The basic idea is as follows:(1)According to the need of aeroengine rotor vibration reliability analysis, the features of the aeroengine are determined by combining the comprehensive modal theory with aeroengine operation principle(2)Acquire QAR data, extract input variables, and preprocess data(3)Combined with the DEA method, the vibration reliability *r* of the aeroengine during operation is evaluated comprehensively(4)According to the vibration reliability analysis, a large number of samples are extracted as training samples(5)The number of nodes in each layer of the neural network is defined, and the neural network model is established(6)The improved PSO is used to search the initial optimal weight and threshold(7)The network is trained by the BR algorithm to establish INNMF

The basic flow chart of rotor vibration reliability analysis of the aeroengine with INNMF is shown in Figure 2.

#### 4. Case Study

This study takes B737-800 aircraft as an example and takes QAR data of a flight from Beijing to Urumqi of an aviation company as an example to analyse the rotor vibration reliability of its aeroengine during operation. The QAR data containing 21 features of four factors is extracted. If the feature value of a line is missing, the previous value fills in the missing value. The B737-800 selected in the calculation example is equipped with two CFM56-7B engines, which have excellent performance, abundant sensors, and good safety. As mentioned in Section 2, DEA is used to evaluate and analyse the reliability of rotor vibration during aeroengine operation. Then, INNMF is used to fit the data to study the importance of each feature. According to the calculation results, the vibration reliability law of the aeroengine rotor during operation is summarized.

##### 4.1. Evaluating the Reliability of Rotor Vibration *R*

Assuming that the rotor vibration margin *R* represents the rotor vibration reliability during the operation of the aeroengine, and the set **D** represents the values of the extracted 21 features, and the relationship between *R* and **D** can be expressed as

According to the rotor vibration reliability method during the operation of the aeroengine proposed in Section 2, combined with the extracted QAR data, the time-varying rotor vibration margin *R* during the operation of the aeroengine of the case is analysed. The results are shown in Table 2.

Table 2 shows the intercepted QAR data and *R* during take-off, cruise, and landing phases. Time 9:00:01–9:00:04 is intercepted from the take-off phase. As the aircraft climbs, the barometric altitude of aircraft (BALT) and computed airspeed (CAS) continue to increase. At this time, the fuel flow (FF) is large, the exhaust gas temperature (EGT) is high, the rotor keeps running at a high speed, and *R* is around 0.7. Time 10:00:01–10:00:04 is intercepted from the cruise phase. At this time, the aircraft is running smoothly, the BALT remains at 36,096 feet, the EGT and FF are lower than the take-off phase, and the value of *R* is close to 0.9. Time 12:00:01–12:00:04 is intercepted from the landing phase. With the decrease of BALT, FF, low-pressure compressor rotor speed (*N*1), high-pressure compressor rotor speed (*N*2), and EGT decrease obviously. However, the aircraft operates frequently, and the engine runs continuously for more than 3 hours. Compared with the cruise phase, *R* is smaller, around 0.74.

Figure 3 visually shows the distribution and variation of rotor vibration margin *R* during the whole flight process of the case flight. It can be seen that *R* of the take-off phase and landing phase is lower than that of the cruise phase, and the fluctuation of the take-off phase and landing phase is larger. In the take-off phase, *R* drops sharply to less than 0.1 when the aircraft climbs with high thrust after leaving the ground and then quickly recovers to around 0.8. This is due to the fact that FF, *N*1, and *N*2 are very large in high thrust climb, the engine is close to full load operation, and the rotor vibrates violently, resulting in low *R*. In the cruise phase, *R* is high, mostly between 0.8 and 1, but there are some fluctuations. This is because in the process of operation, the changes of external temperature and wind speed will have a certain impact on the engine operation. In the landing stage, the value of *R* is low on average, widely distributed and fluctuated. Because this process is gradient descent and environmental factors are constantly changing, *R* is not stable.

##### 4.2. Establishing the IESMF Mathematical Model

According to the calculation results in Table 2, 3000 groups of data are extracted by the random sampling method according to the time interval. Among them, 2400 are used as training samples to create INNMF, and 600 are used as test samples for INNMF verification.

To make INNMF suitable for these training samples, according to the characteristics of input variables and output response, the ANN model uses a 21-8-1 three-layer structure. Choose “tansig” as the transfer function from the input layer to the hidden layer; “purelin” as the transfer function from the hidden layer to the output layer; “trainbr” as the training function. The number of particles *N* = 40 and the particle dimension = 185 are selected for the model. After 100 iterations, Figure 4 shows the optimal fitness value curves.

The network weight threshold optimized by improved PSO is assigned to the network model, and the INNMF is obtained after training by the BR algorithm. The weights and thresholds are shown in the following equation:

##### 4.3. Regression Prediction

The regression prediction of *R* in the process of aeroengine operation is carried out by using INNMF, and the results are shown in Figure 5. Figure 6 shows the prediction error of INNMF.

Figure 7 shows the importance of each feature to the vibration margin *R*. The five characteristics that have great influence on *R* are TAT, BALT, wind speed (WSFMC), angle of attack (AOAL and AOAR), and Greenwich Mean Time (GMT). The aircraft relies on aerodynamics to achieve flight and control. The flight management computer can adjust different aircraft attitude and flight speed to make the aircraft in a safe envelope. The direct and effective way to control the speed is to control the thrust of the engine. In the whole process of operation, the altitude and speed are constantly changing, and the external environment is also changing, so the flight time is closely related to the environmental factors (TAT and WSFMC) and BALT. In addition, the change of the angle of attack represents the adjustment of the aircraft operating status, which is usually accompanied by the adjustment of engine performance, so the change of the angle of attack will affect *R* to a certain extent. On the contrary, the characteristics that have the least influence on the value of the rotor vibration margin *R* during the operation of the aeroengine are the oil quantity (EGOILQ1 and EGOILQ2). In summary, combined with the actual operation of the aircraft and expert experience, the above characteristic importance analysis is reasonable. To sum up, the above analysis results are in line with the actual operation of the aircraft and the experience of experts, which is reasonable and feasible.

##### 4.4. Validation of INNMF

To verify the applicability and performance of the INNMF proposed in this research, random forest (RF) and ANN are used to fit and predict the rotor vibration margin *R*, respectively. RF is a machine learning algorithm based on the decision tree proposed by American scholar Leo Breiman [29]. As an integrated learning method, it has the advantages of superior performance and simple structure. The calculation time and accuracy are shown in Table 3.

As shown in Table 3, as the sample size increases, the calculation time of INNMF is less than that of RF and ANN. As the simulation time increases, the computational efficiency of INNMF is higher than that of ANN. The results show that INNMF, which combines ANN and improved PSO and uses BR training, has high computational efficiency in the analysis of rotor vibration reliability during aeroengine operation. The reason is that (i) INNMF can quickly fit and train uncertain parameters. (ii) The neural network improved by PSO can quickly obtain accurate initial weights and thresholds, thus saving a lot of time and improving the efficiency of analysis.

Regarding the calculation accuracy of INNMF, almost all sampling points are fitted by it with low training error, showing the strong generalization ability and nonlinear fitting ability of INNMF. As shown in Table 3, compared with RF and ANN, INNMF is more precise. The reasons are as follows: (i) INNMF has strong nonlinear mapping ability and generalization ability in the training process. The obtained more accurate mathematical model provides a guarantee for the calculation accuracy. (ii) Combined with the improved PSO algorithm, the search accuracy of global initial optimal weights and thresholds is improved. (iii) The BR algorithm reduces the training error by reducing the weight of the network.

To sum up, the INNMF guarantees the calculation accuracy and greatly improves the calculation efficiency. This method has good adaptability to the reliability analysis of aeroengine rotor vibration during operation, and the comprehensive performance is satisfactory.

#### 5. Conclusions

(1)The proposed method for analysing the reliability of rotor vibration during aeroengine operation, which considers 4 types of factors and 21 characteristics, can effectively calculate the time-varying*R.*Compared with the traditional method, it can reflect the reliability of aeroengine short-term vibration, which has the significance of theoretical exploration.(2)Based on the 21 proposed features, combined with the calculation results of

*R*, the importance of each feature’s influence on

*R*during the operation of the aeroengine is obtained. The six most important features are TAT, BALT, WSFMC, AOAL, AOAR, and GMT.(3)By comparing with traditional methods, the results show that the INNMF proposed in this research has high efficiency and high precision for reliability analysis of rotor vibration during engine operation. As the number of sample size increases, the advantages of INNMF are more obvious.(4)The work of this paper enriches the theory of rotor vibration reliability analysis and provides a certain reference for the operation monitoring, maintenance, and optimization design of aeroengines.

#### Data Availability

The data used to support the findings of this study are included within this article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 51875465). The authors are grateful for the support.