An Iterative Reduced KPCA Hidden Markov Model for Gas Turbine Performance Fault Diagnosis

To improve gas-path performance fault pattern recognition for aircraft engines, a new data-driven diagnostic method based on hidden Markov model (HMM) is proposed. A redundant sensor somewhat interferes with fault diagnostic results of the HMM, and it also increases the computational burden. The contribution of this paper is to develop an iterative reduced kernel principal component analysis (IRKPCA) algorithm to extract fault features from original high-dimension observation without large additional calculation load and combine it with the HMM for engine gas-path fault diagnosis. The optimal kernel features are obtained by iterative sequential forward selection of the IRKPCA, and the features with lower dimensions are contracted through a trade-off between the fault information and modeling data scale in reduced kernel space. The similarity degree is designed to simplify the HMM modeling data using fault kernel features. Test results show that the proposed methodology brings a significant improvement in diagnostic confidence and computational efforts in the applications of a turbofan engine fault diagnosis during its steady and dynamic process.


Introduction
The gas turbine engine is the power source of aircraft, and its reliability directly affects aircraft safety and performance [1].The engine is an easy-fault piece of machinery since it has complex structure and runs in harsh operating conditions.During the course of an engine's life, various physical failures might happen, such as corrosion, erosion, fouling, and foreign object damage [2,3].These failures lead to gas-path performance degradation, either gradually or abruptly, which is recognized as engine gas-path fault and is greatly harmful to flight safety [4,5].For the purpose of enhancing operating reliability and reducing maintenance costs of aircraft propulsion systems, engine gas-path fault diagnosis technology has attracted interest.
Generally speaking, gas turbine fault diagnosis approaches are divided into model-based and data-driven ones [6,7].Variants of Kalman filters are the typical model-based diagnostic methods for gas turbine engines [8,9].It requires a reliable engine model, which relies on physical characteristics and aero-thermodynamic theory.The engine modeling uncertainties and operation non-linearity during the transient process negatively affects the performance of model-based technologies [10].Engine-to-engine variation makes it difficult for the general engine model to represent every individual engine.In addition, the observability of general Kalman filters relates to the number of sensor measurements and their types, and it cannot recognize fault patterns as available sensors less than health parameters.
The data-driven approach is another important method of fault diagnosis for complex non-linear systems, especially in the rich data circumstance [11][12][13].It depends on the collected data of fault modes but not an accurate mathematical model of a gas turbine engine.It is particularly suitable to a complex strong non-linear system, and it is not limited to the sensor measurement number.In the data-driven field, much attention has been paid to neural network approaches, which are theoretically developed from empirical risk minimization with simple mathematical expressions [8].Bartolini applied neural networks to micro gas turbines, and then Amozegar developed an ensemble of dynamic neural network identifiers for engine fault detection and isolation [14,15].The desirable topological structure of a neural network is usually selected by experience, and the diagnostic performance by a neural network easily fluctuates with stochastic measurement noise.
Different from the neural network approaches, the hidden Markov model (HMM) is a classic data-driven fault diagnosis for non-linear stochastic systems [16].The HMM has rigorous theoretical deduction and definite model structure [17].HMM statistical characteristics in modeling and classification makes it outperform fault pattern recognition for a mechanical system with clear randomness and uncertainties [18].However, it is noted that the HMM's computational cost increases dramatically when the measurement dimension increases.The sensor data from various operating cross-sections for engine gas-path fault diagnosis is usually complex and physically correlative in transient process in the flight envelope.Consequently, it is necessary to extract fault features from raw measurement sequences to decrease test data dimensions and simplify the HMM structure.
Principal component analysis (PCA) is introduced into the HMM to extract fault features to reduce computational effort.PCA is achieved by projecting the linear data matrix onto an uncorrelated subspace with less information loss, but its performance is reduced in the plant with strong non-linearity.Kernel-PCA (KPCA) is developed to overcome shortcomings of conventional PCA-to-linear issues [19,20].Provided a kernel matrix (an n × n matrix where n is the number of the dataset) to map an original dataset into the feature space, the KPCA computational complexity is O(n3) in the principal component feature extraction.Taouali proposed a novel RKPCA to improve the sparsity capability [21], but it is difficult to balance the useful information capacity and data scale to reach the appropriate sample number.
To improve fault diagnostic confidence and computational efforts, a novel fault diagnostic approach is proposed using the combination of iterative reduced KPCA and HMM for engine gas-path fault diagnosis.In this paper, the IRKPCA is developed with a sample-reduction mechanism, which is designed to decline the redundant information in the initial observation sequences for gas-path fault feature extraction.The similarity degree and forward kernel inverse are employed to simplify the sample data, and then the kernel fault feature by the IRKPCA is used by the HMM to perform gas-path fault pattern recognition.The systematical tests are carried out to evaluate fault diagnosis performance of the proposed methodology, and it runs on a two-spool turbofan engine simulation in the steady and transient process at various cycle numbers during its life.The results indicate the superiority of the IRKPCA-HMM, and it supports our viewpoints.
The remainder of paper is organized as follows.The IRKPCA is developed from the basic KPCA, and the comparisons of the involved KPCAs are followed by feature extraction performance using benchmark datasets in Section 2. In Section 3, the IRKPCA-HMM is presented by the combination of IRKPCA and HMM to simplify the fault diagnostic model using reduced-kernel fault features.Simulation and analysis are given on a turbofan engine in dynamics for gas-path fault diagnosis in Section 4. Section 5 draws a conclusion and discusses future research directions.

KPCA
As the non-linear extension of PCA, KPCA has been applied to various engineering fields and exhibited satisfactory performance in complicated non-linear systems [19,22].The kernel transformation of the KPCA maps the measured data space to higher-dimension feature space.Given a set of N mean-variance scaled training data, and L is the original dimension of training measurement.The covariance matrix in the feature space F is where transformation matrix φ( ) : R L → R H is a non-linear map from input vector into H-dimensional feature space F. Non-zero eigenvectors γ of the covariance matrix C is calculated by eigenvalue decomposition where λ is the eigenvalues of C in feature space.Combine two expressions in Equation ( 2) and multiply the kernel matrix φ(x k ) in both sides, and we can obtain Then, a kernel matrix K = {k ij } with the dimension of N × N is introduced to achieve the inner product in non-linear feature space A radial basis kernel k(x i , and σ is a dispersion factor.The kernel matrix K is centralized before kernel transformation where I is an identity matrix with the size of N × N.Then, Equation ( 3) is specified as The kernel principal components are selected as the p largest positive eigenvalues λ 1 ≥ λ 2 ≥ • • • ≥ λ p > 0, and their eigenvectors 1 , . . ., p form the kernel principal component set.The projection parameter of a new data x new onto the eigenvectors in the principal component set is called the KPCA-transformed feature variable.The i-th feature z i for x new is expressed From Equation (7), the kernel matrix is computed to obtain the feature vector for each sample in training dataset.The burden of computation and storage would be dramatically serious for drawing Energies 2018, 11, 1807 4 of 21 fault features as the increase of kernel matrix calculation number and the sample data accumulation.For these purposes, an iterative reduced KPCA (IRKPCA) is developed from the KPCA algorithm.

IRKPCA
Given the selected samples set from the training dataset X s {x i }, i = 1, • • • , N s , N s means the selected sample number.The base vectors in feature subspace Φ = {φ(x 1 ), . . . ,φ(x s )} related to the selected samples are sufficient to express the whole transformed data in feature space.Thus, the mapping estimate of any sample vector x i can be represented by a linear combination of these base vectors where a i is the coefficient vector on the feature subspace basis denoted as A new objective function to find these coefficients is designed, and it includes two minimization elements: the estimated mapping errors and coefficient vector norm.The objective function is expressed as where ρ is the regularization parameter to weight off the estimation error and coefficient scales.
The former part indicates that the modeling-variable mapping is as close as possible to the real mapping, and the latter one compacts the scales.Let (∂J i /∂a i ) = 0, we have (10) where I s is an identity matrix of the size of s × s.Combining Equations ( 9) and (10), the solution of minimization J i can be rewritten (11) where K ss represents the kernel matrix of the selected vectors, and K si is the dot product between x i and the selected base vectors, K si = [k(x 1 ,x i ), . . ., k(x s ,x i )] T .Since the calculation procedure of objective function is treated to each sample in the training dataset, the overall objective function is defined by Equation ( 12) can be rewritten as The sample reduction procedure is an iterative sequential forward process.The sample chooses preforms from the training dataset at each step of the procedure, which are combined with the prior selected ones to generate the optimal objective function of Equation (13).The procedure stops as the mean of overall objective is no longer increasing clearly when the sample is added.The mean residual is given by Energies 2018, 11, 1807 To reduce the number of training sample, the similarity degree is used.The similarity degree between a new sample x i and the selected samples Kernel matrix inverse is calculated when a new training sample is selected during the iteration.It is very time-consuming to repeatedly calculate the K ss inverse.The forward kernel inverse strategy is used to compute the (s + 1)-th sample given base vectors X s and R s = (K ss + ρI s ) −1 Provided an invertible matrix A, and matrices D, V, U, we have Since the matrix R s has been obtained at the s-th iteration, R s can be directly acquired by applying Equation (17) to Equation ( 16) for the matrix inversion where , and the kernel inverse number becomes 1 as a new training sample is considered.The reduction scheme of IRKPCA training sample is developed by simplifying training data number for fault feature extraction.In the IRKPCA, the sample similarity degree is checked in advance, and the sample is such that the similarity index larger than its threshold is deleted.The forward kernel inverse strategy is employed to produce the matrix R s+1 with respect to new training data.

Feature Extraction Performance Test
The benchmark datasets from UCI Machine Learning [23] are used to validate feature extraction performance of the basic KPCA and proposed IRKPCA.The six datasets include Iris, Wine, Glass, E.coli, Balance and Vehicle.The feature number (#Feature), the number of classes (#Class) and sample number (#SN) of each dataset are listed in Table 1.The comparisons of basic KPCA and IRKPCA are investigated from reduced features, discriminative power, sparsity, and reduced time in Table 2.The reduced feature number depicts the principal component size after feature extraction by the KPCA.Discriminative power contains the mean correct classification ratio and standard errors using the reduced features in the transformed space.Sparsity is expressed by 2-norm/1-norm (L 2 /L 1 ), and the larger L 2 /L 1 value is more in favor of classification [24,25].The reduced-dimension time gives the computational efforts of dataset transformation from the original measurement space to feature space.The regularization parameter ρ and kernel parameter σ are yielded separately from the bound of [2 -3 , 2 3 ] and [2 0 , 2 10 ] by the interval 2 times.From Table 2, we can see that both KPCA algorithms reduce the feature number, and the number by IRKPCA is less than that by KPCA in the most examined datasets.In terms of the discriminate power index, the IRKPCA produces more than 2.70% and 2.89% of mean correct ratio compared to KPCA in the E.coli and Balance datasets, respectively.The differences of mean correct ratios by the two KPCAs in the rest of the datasets are no more than 1%.The indices L 2 /L 1 by IRKPCA is slightly larger than those by KPCA except in the Glass dataset, so the former produces less information loss than the learning feature.The KPCA generates one magnitude larger reduced time than IRKPCA.The reduced time of KPCA increases faster with the sample size, and the indices are 1.2387 s by KPCA and 0.0645 s by IRKPCA in the Vehicle dataset.The IRKPCA consumes less computational time and owns a lower feature number without discriminative power reduce compared to KPCA.Hence, it implies that the IRKPCA might be a better candidate for gas-path fault feature extraction.Figure 1 shows the extraction effect on benchmark datasets of Iris and Wine by the proposed IRKPCA method.The first three dimensions of benchmark datasets are presented in the form of scatter plots in Figure 1, and the points with the same color belongs to the same cluster.From Figure 1 we can see that the distance of different color points increases after the dimension transformation by the IRKPCA.The processed datasets are more easily distinguished than the original dataset in scatter plots, and it supports further classification.

Hidden Markov Model
Hidden Markov model is extended from Markov chains to yield inferential statistical information on a state sequence [26,27].It comprises a finite hidden state number, and each state relates to an observation at a time point.There are two kinds of probabilities of every hidden state: transition probability and observation probability.HMM have double stochastic processes: the stochastic transition from one hidden state to another one with transition probability and the stochastic output observation generated from every hidden state with observation probability.
The state sequence of HMM is unobservable, but it can be estimated from the observation sequence.Let hidden state sets = { , , ⋯ , } and observation vector = { , , ⋯ , } with ∈ { , , ⋯ , }, where M is the number of hidden states, N is the number of observation per state and T is the length of observation sequence.The variable ( ∈ ) indicates the hidden state, and ( ∈ ) the observation at time t.In general, the compact definition λ = (π, A, B) is used to specify an HMM model, and it includes three components: state transition probability matrix A, observation probability matrix B, and initial state distribution π.
An initial state distribution π illustrates the probability of the starting hidden state si where πi = Pr(q0 = si), A transition probability A contains the element aij { } ( )

Hidden Markov Model
Hidden Markov model is extended from Markov chains to yield inferential statistical information on a state sequence [26,27].It comprises a finite hidden state number, and each state relates to an observation at a time point.There are two kinds of probabilities of every hidden state: transition probability and observation probability.HMM have double stochastic processes: the stochastic transition from one hidden state to another one with transition probability and the stochastic output observation generated from every hidden state with observation probability.
The state sequence of HMM is unobservable, but it can be estimated from the observation sequence.Let hidden state sets where M is the number of hidden states, N is the number of observation per state and T is the length of observation sequence.The variable q t (q t ∈ S) indicates the hidden state, and o t (o t ∈ O) the observation at time t.In general, the compact definition λ = (π, A, B) is used to specify an HMM model, and it includes three components: state transition probability matrix A, observation probability matrix B, and initial state distribution π.
An initial state distribution π illustrates the probability of the starting hidden state s i (t = 0) where Energies 2018, 11, 1807 where a ij is the probability of transition from one hidden state s i at time t to another hidden state s j at time t + 1, and it follows a ij = Pr(q t+1 = s j |q t = s i ), 1 ≤ i, j ≤ M. The observations are emitted from each state according to the conditional probability, b i (v j ) = Pr(o t = v j |q t = s i ), which is represented as a matrix B{b i (o j )} as well.The elements of B evidently satisfy the constraints: The Baum-Welch algorithm is used to obtain the topological parameters of HMM, where an iterative calculation is implemented to maximize the observed sequence probability Pr(X|λ) from an initial λ.The observed sequence X T×H = [x 1 , x 2 , . . ., x T ] T is divided into two sections: X t×H = [x 1 , x 2 , . . . ,x t ] T with t measurement data number, and X (T−t)×H = [x t+1 , x t+2 , . . . ,x T ] T with T − t.H is the dimension of measurement vector.Two probabilities corresponding to the forward variable and backward variable are computed.The former is denoted by the probability of observing the partial measurement data X t×H ending in the state s j , and the latter by the probability of observing the remaining data X (T−t)×H as the state s i .The forward variable α t (j) is where α 1 where β T (i) = 1, 1 ≤ i ≤ M. The detailed Baum-Welch algorithm to obtain the probability of sequence O can be referred to in paper [28].The scalar quantization is to generate the partition vector and the codebook vector according to the quantization distortion of the input before training HMM [29].
The observation sequence of HMM is formed by measurement data from various sensors, which are equipped along a gas path to depict engine operation.As mentioned earlier, redundant information of various operating data will increase fault diagnosis difficulty.The larger dimension of measurement vector x i leads to more complex HMM topology [5,19].It would not only increase the computational time of forward-backward procedure in Baum-Welch algorithm but also negatively impact on classification accuracy.

IRKPCA-HMM Algorithm
Gas-path fault feature extraction based on the IRKPCA is combined with HMM for fault diagnosis in this section.The reduced observation sequence X* T×R = [x * 1 , x * 2 , . . ., x * T ] T with R feature dimensions (R < H) is generated from the original sequence X T×H .Thus, the re-estimation process of HMM with the new sequence X* T×R is as follows Energies 2018, 11, 1807 where a (k) t is k-th forward variable, and t+1 is k-th backward variable.In IRKPCA-HMM test phase, Equation ( 26) is specified as The log-likelihood probability (LL) is a fault index defined by the logarithm calculation to P(X * |λ) in Equation (27).The larger LL indicates the greater consistency of observation sequence to the HMM.
The IRKPCA-HMM achieves gas-path fault diagnosis in the lower dimension feature space from the original measurement space, and it decreases the calculation effort not only in the training process but also in testing process.The IRKPCA-HMM algorithm schematic is presented for fault diagnosis of aero engines in Figure 2.
Energies 2018, 11, x FOR PEER REVIEW 9 of 20 The log-likelihood probability (LL) is a fault index defined by the logarithm calculation to ( *| ) P X λ in Equation (27).The larger LL indicates the greater consistency of observation sequence to the HMM.
The IRKPCA-HMM achieves gas-path fault diagnosis in the lower dimension feature space from the original measurement space, and it decreases the calculation effort not only in the training process but also in testing process.The IRKPCA-HMM algorithm schematic is presented for fault diagnosis of aero engines in Figure 2. The proposed IRKPCA-HMM procedure is as follows: (1) Initialization.
Set initial kernel parameters, regulation parameter and information loss threshold.Select the number of hidden states and observation states, and topological structure of HMM.
(2.2.) Select the first optimal sample.Calculate the objective function Equation (13) for each sample, and the sample x1 with the maximum objective value is added into the sample dataset Xs, (2.3.)Compute the sample similarity degree.The similarity degree is computed between the selected sample in Xs and rest samples xi in t s s − − X X X .If sample similarity degree is larger than similarity threshold, then  The proposed IRKPCA-HMM procedure is as follows: (1) Initialization.Set initial kernel parameters, regulation parameter and information loss threshold.Select the number of hidden states and observation states, and topological structure of HMM.(2) Search optimal samples by iterations.
(2.1.)The original training sample set X 0 , and let the dataset Xs = ∅, Xs = ∅, Xt = X 0 .(2.2.) Select the first optimal sample.Calculate the objective function Equation (13) for each sample, and the sample x 1 with the maximum objective value is added into the sample dataset Xs, Xs = Xs ∪ {x 1 }.(2.3.)Compute the sample similarity degree.The similarity degree is computed between the selected sample in Xs and rest samples x i in Xt − Xs − Xs.If sample similarity degree is larger than similarity threshold, then Xs = Xs ∪ {x i }, i = I + 1, and continue until the rest samples are all checked.
(2.4.)Calculate the objective Equation ( 18) for each sample and select the sample of maximum objective value x i .Compute the mean objective residual as Equation ( 14) for the selected optimal sample x i .If the sample mean objective residual is larger than the difference threshold, go to 1.3, else continue.
(3) Principal component computation in kernel space.
(3.1.)Calculate the centralized kernel matrix K with selected samples.
(5.3.)Form the HMM libraries for gas-path fault modes.

IRKPCA-HMM Based Engine Gas-Path Fault Diagnosis
The proposed IRKPCA-HMM approach to gas-path fault diagnosis is tested on a virtual two-spool turbofan engine developed by the component-level engine model [30,31].The examined turbofan engine is mainly composed of inlet, fan, compressor, bypass, combustor, high-pressure turbine (HPT) and low-pressure turbine (LPT), mixer and nozzle, and it is illustrated in Figure 3.The inlet supplies airflow into the fan, and then the air is divided to two streams: one flowing into the compressor and the other passing through the bypass.Air leaving the compressor moves to the combustor, where fuel is injected and burns to produce hot gas to drive the turbines.The fan and compressor are driven by the LPT and HPT, respectively.Gas from LPT and air from bypass mix in the mixer, and then leaves the engine through the nozzle.Closed-loop control strategy of spool speed is applied to aero engine with safety protection [32].The engine station numbers in Figure 3  (5) HMM modeling with reduced observation.
(5.3.)Form the HMM libraries for gas-path fault modes.

IRKPCA-HMM Based Engine Gas-Path Fault Diagnosis
The proposed IRKPCA-HMM approach to gas-path fault diagnosis is tested on a virtual twospool turbofan engine developed by the component-level engine model [30,31].The examined turbofan engine is mainly composed of inlet, fan, compressor, bypass, combustor, high-pressure turbine (HPT) and low-pressure turbine (LPT), mixer and nozzle, and it is illustrated in Figure 3.The inlet supplies airflow into the fan, and then the air is divided to two streams: one flowing into the compressor and the other passing through the bypass.Air leaving the compressor moves to the combustor, where fuel is injected and burns to produce hot gas to drive the turbines.The fan and compressor are driven by the LPT and HPT, respectively.Gas from LPT and air from bypass mix in the mixer, and then leaves the engine through the nozzle.Closed-loop control strategy of spool speed is applied to aero engine with safety protection [32].The engine station numbers in Figure 3   The data are generated from the numerical engine model [33,34] to evaluate the involved methods in the steady behavior of the maximum power operation and transient behavior including acceleration and deceleration.The involved engine parameters are reported in Table 3 The data are generated from the numerical engine model [33,34] to evaluate the involved methods in the steady behavior of the maximum power operation and transient behavior including acceleration and deceleration.The involved engine parameters are reported in Table 3.The control variables include fuel flow W f and Nozzle area A 8 , which define the operating point of the engine.The health parameters are unmeasurable and represent engine gas-path health, containing indicators of fan efficiency SE 1 , fan flow SW 1 , compressor efficiency SE 2 , compressor flow SW 2 , HPT efficiency SE 3 , HPT flow SW 3 , LPT efficiency SE 4 , and LPT flow SW 4 .The available measurements are used to calculate health parameters, and they are low-pressure spool speed N L , high-pressure spool speed N H , compressor inlet temperature T 22 , compressor inlet pressure P 22 , compressor outlet pressure P 3 , compressor outlet temperature T 3 , LPT inlet temperature T 43 , LPT inlet pressure P 43 , LPT outlet pressure P 5 and mixing chamber inlet temperature T 6 [35].The maximum power point on the ground is defined as corrected percentage of high-pressure spool speed N Hcor = 100%, and corresponds to the corrected normalized values of engine control variables: fuel flow W f = 100%, nozzle area A 8 = 47%.The measurement noise follows time-uncorrelated zero-mean Gaussian noise, and the magnitude of these noises can be referred to in paper [36].Both gradual and abrupt performance deterioration causes health parameter variations.The health parameters resulting from performance gradual degradation is long term, and all health parameters synchronously diverge from their nominal quantities with the cycle number increase.It starts from a healthy engine (all health parameters at their nominal values) at initial cycle number CN = 0, and with the linearly deviation at the end of cycle number CN = 6000.The first factory overhaul occurs at one quarter of the engine's whole lifetime, and three cycle number points before this overhaul, including CN = 0, CN = 807 and CN = 1558, are addressed in this paper.Table 4 shows health parameter deviations under gradual degradation with regard to cycle numbers.The health parameters move suddenly from their nominal values in gas-path fault scenarios, and the shift quantities of each fault case are given in Table 5.There are thirteen operating scenarios in total at one cycle number, including twelve fault cases and a no-derivation case.Sensor malfunction, such as bias or drift, is not considered in this study.The historical measured data sample is used offline to build up gas-path fault HMM libraries of the engine by the proposed methodology in training stage.Every gas-path fault case relates to one HMM fault library, and they are independent each other.The IRKPCA-HMM libraries are ζ = {ζ 1 , ζ 2 , . . . ,ζ K } as the count of fault case equals to K. The available engine measurements in Table 3 are recorded online in sequence, and IRKPCA-HMM runs in the left-right type [37].Figure 4 shows a gas-path fault diagnosis framework based on IRKPCA-HMM, and the processed data sequence is fed into HMM libraries and each LL of HMM will be calculated.3 are recorded online in sequence, and IRKPCA-HMM runs in the left-right type [37].Figure 4 shows a gas-path fault diagnosis framework based on IRKPCA-HMM, and the processed data sequence is fed into HMM libraries and each LL of HMM will be calculated.The optimal kernel samples are obtained from an online-sensed sequence by IRKPCA in the test stage.The probabilities related to gas-path fault libraries are calculated from reduced observation by the Baum-Welch algorithm.The index LL is used to recognize gas-path fault mode from the observation, and it belongs to the fault library that owns the largest LL [18].The examined algorithms including HMM, KPCA-HMM and IRKPCA-HMM are performed on a Windows 10 PC with CPU i5-2450 M @2.50 GHZ (Intel, Santa Clara, CA, USA) and 8 GB RAM using MATLAB R2012b software (The MathWorks, Inc, Natick, MA, USA).The Monte-Carlo simulation is conducted, and the performance indices are from ten tries.The correct diagnostic ratio Acc and its standard deviation Std are separated, to assess gas-path fault diagnostic accuracy and stability: where Nc is sample number of correct recognition, Nt is total sample number in one fault scenario, and U is the count of fault scenarios.The optimal kernel samples are obtained from an online-sensed sequence by IRKPCA in the test stage.The probabilities related to gas-path fault libraries are calculated from reduced observation by the Baum-Welch algorithm.The index LL is used to recognize gas-path fault mode from the observation, and it belongs to the fault library that owns the largest LL [18].The examined algorithms including HMM, KPCA-HMM and IRKPCA-HMM are performed on a Windows 10 PC with CPU i5-2450 M @2.50 GHZ (Intel, Santa Clara, CA, USA) and 8 GB RAM using MATLAB R2012b software (The MathWorks, Inc., Natick, MA, USA).The Monte-Carlo simulation is conducted, and the performance indices are from ten tries.The correct diagnostic ratio Acc and its standard deviation Std are separated, to assess gas-path fault diagnostic accuracy and stability:

Fault Diagnosis in Steady Process
where Nc is sample number of correct recognition, Nt is total sample number in one fault scenario, and U is the count of fault scenarios.

Fault Diagnosis in Steady Process
To evaluate fault diagnosis capability of the examined algorithms in steady process, tests are conducted in engine gas-path fault scenarios mixed with gradual degradation at full power on the ground.Gas-path faults are separately injected into the nominal performance deterioration at CN = 0, CN = 807, CN = 1558 in Table 4.The health parameters deviate with the cycle number accumulated over time, and they have an abrupt shift with the constant bias related to every gas-path fault scenario in the steady process.The available measurements are recorded as Table 3, and the hidden state number of HMMs are obtained by searching from 2 to 8 with unit interval.After several tries, the IRKPCA-HMM iteration stop conditions in training process are as follows: iterative step exceeds 100 or convergence error (LL difference between current step and last step) is below 0.01.The similarity degree is 0.9, and the training data for stochastic modeling is the observation sequence with the length of 100 sampling points.There are 1300 samples in total used as training data.Figure 5 gives the effect of gas-path fault feature extraction by IRKPCA in steady behavior.Similar to that of benchmark dataset, the first three dimensions of aero engine dataset are presented in the form of scatter plots, where points with the same color belong to the same cluster from Figure 5.We have a more distinct version of thirteen engine fault patterns in steady process after feature extraction by IRKPCA.
The test data of ten gas-path fault scenarios are different from their training observation sequences.Table 6 shows maximum log-likelihood probability LL* and correct recognized number Nc by the HMMs at CN = 1558 in the steady process.The IRKPCA-HMM produces the least absolute LL* except in the case of XI and the largest Nc except in cases VIII and XI among the involved HMMs.It implies that IRKPCA-HMM is superior to HMM and KPCA-HMM with confidence and correct ratio regards at CN = 1558 in the steady process.The performance comparisons of HMMs are presented at various cycle numbers in the steady behavior in Table 6, and it shows average performance indices of HMMs in 13 scenarios.Similar to that of benchmark dataset, the first three dimensions of aero engine dataset are presented in the form of scatter plots, where points with the same color belong to the same cluster from Figure 5.We have a more distinct version of thirteen engine fault patterns in steady process after feature extraction by IRKPCA.
The test data of ten gas-path fault scenarios are different from their training observation sequences.Table 6 shows maximum log-likelihood probability LL* and correct recognized number Nc by the HMMs at CN = 1558 in the steady process.The IRKPCA-HMM produces the least absolute LL* except in the case of XI and the largest Nc except in cases VIII and XI among the involved HMMs.It implies that IRKPCA-HMM is superior to HMM and KPCA-HMM with confidence and correct ratio regards at CN = 1558 in the steady process.The performance comparisons of HMMs are presented at various cycle numbers in the steady behavior in Table 6, and it shows average performance indices of HMMs in 13 scenarios.The fault feature number and optimal hidden state number of basic HMM are the same at three cycle numbers, and they are larger than those of KPCA-HMM and IRKPCA-HMM.From Table 6, KPCA-HMM and IRKPCA-HMM have much simpler topological structure than the basic HMM due to less reduced feature number and hidden state number.The quantities of confusion matrix and transition matrix of IRKPCA-HMM in fault mode 1 are shown in Table 7.The fault diagnostic accuracy indices of Acc, Std, and execution time t test by three HMMs are discussed in Table 8.The Acc of KPCA-HMM and IRKPCA-HMM are clearly larger than that of HMM, while Std are smaller than that of HMM at three cycle numbers.The larger Acc represents better confidence of fault diagnosis result, and the less Std illustrates better stability.Hence, both of KPCA-HMM and IRKPCA-HMM have better fault diagnostic confidence and stability than basic HMM, and KPCA-HMM is a little worse than IRKPCA-HMM.When it comes to executing time t test , KPCA-HMM and IRKPCA-HMM consume less time compared to basic HMM due to the former two having more simplified topology.It is also found that IRKPCA-HMM has almost half the executing time of KPCA-HMM.The performance index Acc decreases a bit with cycle number accumulation over time, while computational time of three HMMs are hardly changed at three cycle numbers.The IRKPCA-HMM produces the least Std in all cases in Table 8.It implies that IRKPCA-HMM has more outstanding diagnostic accuracy, stability and computational efforts compared to the HMM and KPCA-HMM.Hence, it is a satisfactory method of gas-path fault diagnosis in the steady behavior of turbofan engine.In addition, the tests are simulated in the steady process of high-altitude operation (H = 5000, Ma = 1, W f = 100%, A 8 = 52%), and the results are as shown in Table 9.As seen from Table 9, the performance indices of HMM, KPCA-HMM and IRKPCA-HMM are similar to those in the steady process of ground operation.

Fault Diagnosis in Transient Process
The transient test is performed including acceleration and deceleration in the flight envelope to further reveal the performance of proposed methodology for gas-path fault diagnosis.The engine starts from the idle (W f = 68%, A 8 = 100%), and gradually increases to full power (W f = 100%, A 8 = 47%).After dwelling 0.5 s it moves sharply back to the idle, and the whole operation lasts 9.5 s on the ground.The variations of control variables are shown in the transient behavior in Figure 6.The deviation quantity of the combination of gradual degradation and abrupt degradation are added into health parameters, and the simulations run at three cycle numbers.
starts from the idle (Wf = 68%, A8 = 100%), and gradually increases to full power (Wf = 100%, A8 = 47%).After dwelling 0.5 s it moves sharply back to the idle, and the whole operation lasts 9.5 s on the ground.The variations of control variables are shown in the transient behavior in Figure 6.The deviation quantity of the combination of gradual degradation and abrupt degradation are added into health parameters, and the simulations run at three cycle numbers.The sampling rate of 10-dimension measurements is 0.1 s, and the length of observation sequence for training is 95.The hidden state number and the iteration stop parameters in the dynamics are set as the same as those in Section 4.1.There are 1235 samples in total used as training data for 13 fault scenarios, and average training time of these fault scenarios by HMM, KPCA-HMM and IRKPCA-HMM are 23.27 s, 19.53 s and 16.89 s, respectively.We can find that the training computational efforts of IRKPCA-HMM are the least among the examined algorithms.The scale of test data for each gas-path fault scenario is 10 observation sequences, and the length of every sequence is the same as the training one.The indices of LL* and Nc by the examined HMMs at CN = 0 in the transient process of ground operation is given in Table 10.The sampling rate of 10-dimension measurements is 0.1 s, and the length of observation sequence for training is 95.The hidden state number and the iteration stop parameters in the dynamics are set as the same as those in Section 4.1.There are 1235 samples in total used as training data for 13 fault scenarios, and average training time of these fault scenarios by HMM, KPCA-HMM and IRKPCA-HMM are 23.27 s, 19.53 s and 16.89 s, respectively.We can find that the training computational efforts of IRKPCA-HMM are the least among the examined algorithms.The scale of test data for each gas-path fault scenario is 10 observation sequences, and the length of every sequence is the same as the training one.The indices of LL* and Nc by the examined HMMs at CN = 0 in the transient process of ground operation is given in Table 10.From Table 10, the indices of LL* and Nc by IRKPCA-HMM are the largest ones in the most fault cases as engine experiences from idle to full power and then back to idle.The topological parameters of three HMMs and average performance indices of all fault scenarios at CN = 0, CN = 807, CN = 1558 are reported in Table 11.The fault feature number and optimal hidden state number of basic HMM are clearly larger than the rest HMMs, and it means that the topologies of the latter two HMMs are simplified.This is positive for reducing the computational time of gas-path fault diagnosis.Both KPCA-HMM and IRKPCA-HMM produce the similar performance indices of Acc and Std, which outperforms basic HMM.The confidence and stability of fault diagnostics are improved by fault feature extraction of KPCAs.When the performance index t test is concerned, KPCA-HMM is obviously different from IRKPCA-HMM and no longer better than basic HMM.The feature extraction dominates computational time of gas-path fault diagnosis in the transient process.The IRKPCA-HMM consumes less time for feature extraction due to the forward kernel inverse and sample simplification scheme, and it is the best way of weighting off the diagnostic accuracy and computational efforts.
Furthermore, transient performance tests of the proposed methodology are in the flight envelope, and control variables fuel flow Wf and nozzle area A8 change along flight operation H and Ma.The engine starts from the ground point H = 0, Ma = 0, climbs to high altitude (H = 5000, Ma = 1), and the whole operation lasts 10 s.It runs at the full power operation using closed-loop control strategy of spool speed.Figure 7 shows the change route of four input variables during transient process in the flight envelope.The deviation quantities related to each fault scenario are initially added into health parameters as well as that of the ground operation.From Table 10, the indices of LL* and Nc by IRKPCA-HMM are the largest ones in the most fault cases as engine experiences from idle to full power and then back to idle.The topological parameters of three HMMs and average performance indices of all fault scenarios at CN = 0, CN = 807, CN = 1558 are reported in Table 11.The fault feature number and optimal hidden state number of basic HMM are clearly larger than the rest HMMs, and it means that the topologies of the latter two HMMs are simplified.This is positive for reducing the computational time of gas-path fault diagnosis.Both KPCA-HMM and IRKPCA-HMM produce the similar performance indices of Acc and Std, which outperforms basic HMM.The confidence and stability of fault diagnostics are improved by fault feature extraction of KPCAs.When the performance index ttest is concerned, KPCA-HMM is obviously different from IRKPCA-HMM and no longer better than basic HMM.The feature extraction dominates computational time of gas-path fault diagnosis in the transient process.The IRKPCA-HMM consumes less time for feature extraction due to the forward kernel inverse and sample simplification scheme, and it is the best way of weighting off the diagnostic accuracy and computational efforts.
Furthermore, transient performance tests of the proposed methodology are implemented in the flight envelope, and control variables fuel flow Wf and nozzle area A8 change along flight operation H and Ma.The engine starts from the ground point H = 0, Ma = 0, climbs to high altitude (H = 5000, Ma = 1), and the whole operation lasts 10 s.It runs at the full power operation using closed-loop control strategy of spool speed.Figure 7 shows the change route of four input variables during transient process in the flight envelope.The deviation quantities related to each fault scenario are initially added into health parameters as well as that of the ground operation.12.
Table 12.LL* and Nc by three HMMs during dynamic process in the flight envelope at CN = 0.   12. From Table 12, the indices of LL* and Nc by IRKPCA-HMM at CN = 0 are the largest ones in the most fault cases in the flight envelope.The topological parameters of three HMMs and average performance indices in all fault scenarios at CN = 0, CN = 807, CN = 1558 are reported in Table 12.The fault feature number and optimal hidden state number of basic HMM are clearly larger than the others.It indicates that the topological structure of the HMMs is clearly simplified after feature extraction by the KPCA and IRKPCA, and it is positive for reducing computational time of gas-path fault diagnosis.

Conclusions
This paper develops a systematic approach to fault feature extraction and pattern recognition, which leads to an improved data-driven fault diagnosis method.The novelty of this methodology lies in the development of IRKPCA and HMM in combination to facilitate gas-path fault diagnosis for turbofan engines.The reduced samples from IRKPCA in feature space decrease the measurement dimension while the principal information of fault feature is retained.The IRKPCA is evaluated using general benchmark datasets, and the results reveal that IRKPCA is superior to plain KPCA regarding discriminative power, sparsity and reduced dimension time.The simplified observation sequence by IRKPCA is utilized by HMM to develop an IRKPCA-HMM algorithm.The goal of this methodology is to increase gas-path fault diagnostic accuracy and relieve computational effort both in steady and transient behaviors.The proposed methodology is evaluated in the scenarios of gas-path abrupt fault mixed with gradual degradation in the flight envelope, and test data are generated from a dual-spool turbofan engine model.The stochastic diagnostic modeling framework is presented and numerically assessed by several performance indices.The advantage of the proposed methodology is that it does not only produce more reliable results but also consumes less computational efforts of fault diagnosis.
This research establishes a new direction in data-driven fault diagnosis by proposing IRKPCA-HMM technique that is specifically beneficial to gas-path stochastic fault diagnosis for turbofan engine applications.The methodology developed in this study is not only limited to turbofan engine, but also extended to other types of gas turbine engine.There are some important topics for further study related to this work.First, further studies can be carried out to investigate various kernels used to map the measurement space to feature space.Second, extensions of the cases that have more than one gas-path abrupt fault, added to gradual degradation and the tests of semi-physical hardware in the loop, are worthy of future exploration.

Figure 1 .
Figure 1.Extraction effect on Iris dataset and wine dataset by the IRKPCA.

Figure 1 .
Figure 1.Extraction effect on Iris dataset and wine dataset by the IRKPCA.

( 2 .
1.) The original training sample set X0, and let the dataset 0 and continue until the rest samples are all checked.(2.4.)Calculate the objective Equation (18) for each sample and select the sample of maximum objective value xi.Compute the mean objective residual as Equation (14) for the selected optimal sample xi.If the sample mean objective residual is larger than the difference threshold, go to 1.3, else continue.
are as follows: inlet exit marked by 2, compressor inlet by 22, compressor exit by 3, HPT entrance by 43, LPT entrance by 5, and LPT exit by 6.

Figure 4 .
Figure 4. Gas-path fault diagnosis framework for turbofan engine based on IRKPCA-HMM.

Figure 4 .
Figure 4. Gas-path fault diagnosis framework for turbofan engine based on IRKPCA-HMM.

Energies 2018 ,
11, x FOR PEER REVIEW 13 of 20 (a) Original data (b) Data with feature extraction

Figure 5 .
Figure 5.The effect comparison of gas-path fault feature extraction using IRKPCA.

Figure 5 .
Figure 5.The effect comparison of gas-path fault feature extraction using IRKPCA.

Figure 6 .
Figure 6.Variations of fuel flow and nozzle area of the engine in the transient process.

Figure 6 .
Figure 6.Variations of fuel flow and nozzle area of the engine in the transient process.

Figure 7 .
Figure 7.The change route of engine input variables in the transient process in flight envelope.The observation sequence length of training data is 100, and there are 1300 samples in total used for 13 fault cases.The average training time of fault cases by HMM, KPCA-HMM and IRKPCA-HMM are 26.75 s, 22.43 s and 17.75 s, respectively.The training time of the IRKPCA-HMM is the least among the examined algorithms.The indices of LL* and Nc by the examined HMMs at CN = 0 in the transient process is given in Table12.

Figure 7 .
Figure 7.The change route of engine input variables in the transient process in flight envelope.
length of training data is 100, and there are 1300 samples in total used for 13 fault cases.The average training time of fault cases by HMM, KPCA-HMM and IRKPCA-HMM are 26.75 s, 22.43 s and 17.75 s, respectively.The training time of the IRKPCA-HMM is the least among the examined algorithms.The indices of LL* and Nc by the examined HMMs at CN = 0 in the transient process is given in Table

Table 1 .
Specification of these benchmark datasets.

Table 2 .
Performance comparisons of basic KPCA and IRKPCA on the benchmark datasets.

Table 3 .
Turbofan engine control variables, health parameters, and measured variables.
3 -Compressor outlet pressure SW 3 -HPT flow indicator T 3 -Compressor outlet temperature SE 4 -LPT efficiency indicator T 43 -Low pressure turbine inlet temperature SW 4 -LPT flow indicator P 43 -Low pressure turbine inlet pressure P 5 -Low pressure turbine outlet pressure T 6 -Mixer inlet temperature

Table 4 .
Gradual degradation of turbofan engine performance over time.

Table 5 .
Turbofan engine gas-path performance fault cases.

Table 6 .
LL* and Nc by three HMMs in the steady process at CN = 1558.

Table 6 .
LL* and Nc by three HMMs in the steady process at CN = 1558.

Table 8 .
Fault diagnosis comparisons of the HMMs in the ground steady process.

Table 9 .
Comparisons of fault diagnosis methods in high altitude steady operation.

Table 10 .
LL* and Nc by three HMMs at CN = 0 in the transient process of ground operation.

Table 10 .
LL* and Nc by three HMMs at CN = 0 in the transient process of ground operation.

Table 11 .
Fault diagnosis result comparisons in the transient process of ground operation over time.

Table 11 .
Fault diagnosis result comparisons in the transient process of ground operation over time.

Table 12 .
LL* and Nc by three HMMs during dynamic process in the flight envelope at CN = 0.