Multi-Task Data Imputation for Time-Series Forecasting in Turbomachinery Health Prognostics

Abstract

Time-series forecasting is the core of the prognostics and health management (PHM) of turbomachinery. However, missing data often occurs due to several reasons, such as the failure of sensors. These partially missing and irregular data greatly affect the quality of time-series modeling and prediction as most time-series models assume that the data are sampled uniformly over time. Meanwhile, the training process of models requires a large number of samples and time. Due to various reasons, it is difficult to obtain a significant amount of monitoring data, and the PHM of turbomachinery has high timeliness and accuracy requirements. To fix these problems, we propose a multi-task Gaussian process (MTGP)-based data imputation method that leverages knowledge transfer across multiple sensors and even equipment. Thereafter, we adopt long short-term memory (LSTM) neural networks to build time-series forecasting models based on the imputed data. In addition, the model integrates the methods of denoising and dimensionality reduction. The superiority of this integrated time-series forecasting framework, termed MT-LSTM, has been verified in various data imputation scenarios of a synthetic dataset and a real turbomachinery case.

1. Introduction

Turbomachinery, such as steam turbines, turbocompressors, and gas turbines, is the key equipment type in the electric power, metallurgical, petrochemical, and aviation industries supporting the national economic lifeline. With the development of technology, turbomachinery is heading toward high levels of speed, integration, precision, and intelligence. High availability and reliability standards must be satisfied [1,2]. However, due to their complex structures, special operating environments, and long service time, the failure rate of turbine units is high, which creates serious damage. Therefore, the prognostics and health management (PHM) of turbomachinery, which involves health monitoring, feature extraction, fault diagnosis, and fault forecasting, has emerged as a crucial issue that must be resolved [3,4].

The main research object of PHM is multivariate time-series data. The acquired original data are distinguished by high accuracy, wide coverage, variable patterns, different correlations, and rapid evolution due to the widespread use of inexpensive sensors and the growth of the Internet of Things. Tools for robust multivariate time-series modeling and analysis are urgently needed. The two main kinds of existing approaches are data-driven methodology and physics-model-based methodologies. The physics-model-based approaches must provide a thorough mathematical model to represent the system’s physical properties and fault mechanisms [5,6], whereas, building a physical model of complicated mechanical systems is challenging, and the approaches rely on domain knowledge, which limits the viability and adaptability of these approaches. The challenge of fault diagnosis in complicated mechanical systems is made easier by data-driven approaches. Hence, data-driven methods have gained widespread attention. Data-driven models have experienced two stages [7]. The first is traditional intelligent models based on machine learning (ML). This paradigm extracts statistical features of signals in the time domain, frequency domain, or time–frequency domain and feeds them into ML models such as support vector machines (SVM) [8], k-nearest algorithm [9], or naive Bayes [10] to identify the status of equipment. In traditional ML methods, many efforts have aimed at the manual design of feature extraction algorithms as traditional models cannot extract representative features from raw signals. Such feature-designing processes require abundant professional knowledge and experience, which necessitate a significant amount of human labor, which limits the application of these kinds of methods. The second is modern intelligent models based on deep learning (DL). In contrast, this end-to-end paradigm directly feeds the raw signals (time domain) or pre-processed signals (frequency domain and time–frequency domain) into the deep neural network model to learn the mappings between data and status [11].

Deep learning (DL) is the most popular embranchment of machine learning and has created significant progress in a wide range of areas [12]. Usually, DL-based methods construct a multi-layer structure to extract deeper and more essential features. Recently, various DL-based methods, for example, restricted Boltzmann machines (RBM) [13,14], convolutional neural networks (CNN) [15,16], recurrent neural networks (RNN) [17], auto-encoders (AE) [18], long short-term memory (LSTM) networks [19], or their variants [20,21] have been widely investigated. They can adaptively obtain effective features from time-domain, frequency-domain, and time–frequency-domain data. Then, the system’s operating state can be identified by inputting the features into the regressor or the classifier. However, the training of DL models requires many samples, which is difficult to achieve in practice due to their scarcity. Moreover, the training process requires a lot of time as each layer needs a huge number of training epochs [22,23].

To resolve the first drawback, Lim et al. [24] proposed an expedient assessment framework of LED reliability. Multi-output Gaussian process regression (MOGPR) was performed on sparse reference data to establish an accepfigure threshold for the degradation signature. The framework achieved a reduction in the validation time from 6000 h to 100 h. Jia et al. [25] presented a method named deep normalized convolutional neural network (DNCNN) with normalized layers and a weight normalization strategy. This approach effectively increases the fault diagnosis accuracy with a limited number of fault samples. Ensemble deep auto-encoders (EDAEs), a framework that includes several auto-encoders with various properties, was proposed by Shao et al. [26]. The accuracy of this approach can reach 93.5% when there are 12 health statuses of bearings and 100 fault samples. Based on the generative adversarial network (GAN), Mao et al. [27] proposed a fault diagnosis approach. The fault diagnosis model’s generalization ability was improved by generating synthetic samples for the minority fault class. However, deep neural networks form the foundation for the above methods. We know this requires significant training time and epochs, which impedes the model’s ability to update quickly.

To resolve the second drawback, Guo et al. [23] proposed a CNN–LSTM deep learning model-based method for fault diagnosis on the hydraulic system of the blanket transfer device Mover in the Chinese Fusion Engineering Test Reactor (CFETR). In the fault diagnosis experiment, the model had the highest accuracy of 98.56% on the test set and good efficiency in computation time compared to the other three models. Wang et al. [22] obtained a stable distribution of activation values in training by putting batch normalization (BN) in each layer of SAEs. This method can reduce the training time and epochs. Qin et al. [28] proposed a new activation function, Isigmoid, which combines the traditional Sigmoid with ReLU. The results showed that the approach found a solution to the vanishing gradient in the backpropagation process of the deep network and has advantages in convergence. Zhuang et al. [29] presented a new method called multi-scale deep convolutional neural network (MS-DCNN). This intelligent method extracts characteristics of various sizes in parallel using convolution kernels of different sizes, which has a better learning effect, stronger robustness, and reduced training time. Unfortunately, none of the approaches mentioned above can be used for few-shot learning or missing data.

To overcome the above disadvantages, we design a multi-task time-series data analysis model for the health management of turbomachinery. To account for the scarce-data scenario in the PHM of turbomachinery, we performed data imputation of the target task by leveraging knowledge from auxiliary data (tasks) through the modeling of a multi-task Gaussian process (MTGP), which is a non-parametric multi-task model. Thereafter, we employed long-short term memory (LSTM) to build the sequence model of the target data upon the imputed time-series data to perform a prediction in the future. Meanwhile, we integrated denoising and dimensionality reduction techniques. In comparison to solely modeling the scarce sequence data of the target task, the proposed multi-task LSTM (MT-LSTM) provides more accurate and efficient modeling and predictions, which thus benefits downstream fault diagnosis, remaining useful life estimation, etc.

The remainder of this paper is structured as follows. In Section 2, we expound on the theory and ingredients of the proposed multi-task time-series data analysis model. Section 3 conducts empirical research to verify the superiority of the proposed model. Finally, Section 4 summarizes the conclusions.

2. Methods

For the time-series forecasting of the performance of turbomachinery, missing data is commonly encountered due to the failure of sensors or human error, and it greatly affects the quality of time-series modeling and prediction. Therefore, it is important to fill in the missed time-series data before forecasting. The turbomachinery operation datasets, meanwhile, often have the following characteristics: time series with noise, non-significant periodicity, and high and correlated dimensions (multi-task). Hence, we need to integrate multiple methods, including data denoising, dimensionality reduction, normalization, imputation, and prediction. The flowchart in Figure 1 presents the process for time-series forecasting using datasets collected from turbomachinery. Suppose that we collected scarce time-series data (target task) regarding the interested performance indicator of turbomachinery, and some related and sufficient time-series datasets (auxiliary tasks) collected from other performance indicators are available. The proposed MT-LSTM first performs data denoising using discrete wavelet packet transform (DWPT), followed by a reduction in dimensionality via probabilistic principal component analysis (PPCA). Thereafter, the MT-LSTM conducts the data imputation of the target task through MTGP. Finally, the LSTM model is built upon the imputed data to predict of the target task in the future. The following sections elaborate on the key ingredients in the proposed MT-LSTM model.

2.1. Data Denoising

In the complex working environment of turbomachinery, sensor time-series data usually comprise a wide range of defects, for example, noise, error, and incoherence. Hence, it is necessary to remove the noise or uncertainty characteristics of data before modeling. In past decades, scholars have proposed many methods to remove noise from sensor data. However, few studies have evaluated multiple defects simultaneously. In this paper, we propose Bayesian wavelet denoising and use multiscale wavelet analysis to analyze the contexts in the signal.

Specifically, we adopt discrete wavelet packet transform (DWPT) to denoise and multiscale analyze the original time series. In DWPT, a given one-dimensional time-series signal, $f\left(t_i\right)\left(i=1, 2,\dots ,n\right)$, with $n$ historical data points, is simultaneously decomposed into a sequence of scaling coefficients, $s_j\left(k\right)$, and wavelet coefficients, $w_j\left(k\right)$. Then, we can use the inverse wavelet transform and the DWPT coefficients to reconstruct the time series as follows:

$$\overline{f}(t_i)={\displaystyle \sum _{k\in Z}{\displaystyle \sum _{j\in Z}[s_j(k){\phi }_{j,k}(t_i)+w_j(k){\psi }_{j,k}(t_i)]}}$$

(1)

where $\overline{f}\left(t_i\right)$ is the reconstructed time series; $k$ and $j$ are the time and frequency indices, respectively; ${\phi }_{j,k}\left(t_i\right)$ is the scaling function; ${\psi }_{j,k}\left(t_i\right)$ is the wavelet function; and the double summation means that the wavelet subspace and scaling subspace are simultaneously decomposed into second-level subspaces, which reveal the time and frequency resolution of the signal.

For the benefit of universality, we set $d_{jk}$ to the $k$-th decomposed coefficients at the $j$-th level ($j$ = $j_0$, …, $J$ − 1; $k$ = 0, 1, …, $2^j$ − 1). According to the Bayesian rule, the posterior distribution of the cleaned coefficients, ${\overline{d}}_{jk}$, can be expressed as

$${\overline{d}}_{jk}|{\gamma }_{jk},{\widehat{d}}_{jk},{\sigma }_j^2~N({\gamma }_{jk}\frac{{\tau }_j^2}{{\sigma }_j^2+{\tau }_j^2}{\widehat{d}}_{jk},{\gamma }_{jk}\frac{{\sigma }_j^2{\tau }_j^2}{{\sigma }_j^2+{\tau }_j^2})$$

(2)

where the symbol ~ represents a distribution; $N\left(a,b\right)$ denotes a normal distribution, where $a$ is the mean and $b$ is the variance; ${\widehat{d}}_{jk}$ denotes the decomposed coefficients distribution, ${\widehat{d}}_{jk}|{\overline{d}}_{jk}, {\sigma }_j^2~N\left({\overline{d}}_{jk}, {\sigma }_j^2\right)$; ${\gamma }_{jk}$ is a binary random variable with independent Bernoulli distribution ${\pi }_j$, i.e., P(${\gamma }_{jk}$ = 1) = 1 − P(${\gamma }_{jk }$ = 0) = ${\pi }_j$, which decides whether ${\overline{d}}_{jk}$ are zero (${\gamma }_{jk }$ = 0) or nonzero (${\gamma }_{jk}$ = 1); and the variance ${\tau }_j^2$ denotes the magnitude of ${\overline{d}}_{jk}$ at the $j$-th decomposition level. For all the coefficients in the $j$-th level, ${\tau }_j^2$ and ${\pi }_j$ may be assigned the same values; ${\sigma }_j$ is the standard deviation, which is estimated from the wavelet coefficients of the $j$-th level DWPT decomposition by dividing the median of the wavelet coefficients by a factor.

Then, the Bayesian hypothesis test is applied to threshold the decomposed coefficients by determining whether to accept the null hypothesis, $H_0$: ${\overline{d}}_{jk}$ = 0. The thresholding rule is defined as follows:

$${\tilde{d}}_{jk}={\widehat{d}}_{jk}I({\eta }_{jk}<1)$$

(3)

where ${\eta }_{jk}$ is the posterior odds ratio of ${\gamma }_{jk }$ = 0 versus ${\gamma }_{jk }$ = 1. I(.) denotes an indicator function. When ${\eta }_{jk}$ < 1, I(.) is equal to unity. Thus, $H_0$ is rejected and ${\widehat{d}}_{jk}$ estimates the coefficient ${\overline{d}}_{jk}$. If not, I(.) is equal to zero, and ${\widehat{d}}_{jk}$ is removed. The cleaned time series $\overline{f}\left(t_i\right)$ based on Equation (1) is then reconstructed by using the acquired coefficients ${\tilde{d}}_{jk}$ from Equation (3).

In summary, we can decompose the time-series data into various levels of wavelet coefficients by the DWPT analysis. Then, a Bayesian hypothesis test is carried out for each level of wavelet coefficients to eliminate possible defects. To avoid over-denoising, the ratio of posterior odds Bayesian method offers a simple way to determine whether the decomposed coefficients have any defects. Meanwhile, the DWTP method provides more coefficients denoting additional subtle details of the time-series signal, thus avoiding under-denoising.

2.2. Dimensionality Reduction

In the analysis of multi-variable time-series data of turbomachinery, due to the correlation between variables and the inherent uncertainty of data, too many variable dimensions increase the calculation assumption and the analysis complexity. After cleaning the multivariate time-series data, the PPCA method is mainly used to (1) reduce data dimensionality, (2) decouple the multivariate correlation, and (3) filter data uncertainty. The method plays an important role in improving the accuracy of subsequent data analysis and reducing the calculation time.

Let $y_i=\tilde{g}\left(t_i\right)\in R^q$ represent a $q$-dimensional real number and let $\mathrm{Y} ={[y_1,\dots ,y_n]}^{\mathrm{T}}$ represent the $n$ × $q$ data matrix: there are $q$ variables, and each variable contains $n$ cleaned time-series data points, $\tilde{g}\left(t_i\right)$. Let $\mathsf{\Phi } ={[{\theta }_1,\dots ,{\theta }_n]}^{\mathrm{T}}$ represent the $n$ × $d$ data matrix and ${\theta }_i\in R^d\left(d \le q\right)$ represent $d$ unobservable latent variables (factors), each contains $n$ corresponding positions in the latent space. All ${\theta }_i$ are usually defined as independently distributed Gaussian variables obeying zero mean and unit variance, i.e., ${\theta }_i~\mathrm{N}\left(0,\mathrm{I}\right)$. The latent variable model combines the measurable variable $y_i$ in the relevant data matrix $\mathrm{Y}$ with the corresponding uncorrelated latent variable matrix $\mathsf{\Phi }$. Then, the form of the Gaussian distribution is

$$y_i|\left({\theta }_i, \mathrm{W}, \mathsf{\psi }\right) ~ N\left(\mathrm{W}{\theta }_i+\mu , \mathsf{\psi }\right), i=1, 2,\dots ,n$$

(4)

where $\mathsf{\psi }={\sigma }^2\mathrm{I}$ is the isotropic noise covariance, the $q$ × $d$ weight matrix $\mathrm{W}$ represents the relation between $y_i$ and ${\theta }_i$, and $q$ mean values that are acquired from $\mathrm{Y}$ form the parameter vector $\mu$, i.e., $\mu =\left(1/n\right){\sum }_{i=1}^n{y}_i$. In the lower dimensional latent space, the unadjusted data ${\mathsf{\Phi }}^{*}={[{\theta }_1^{*},\dots ,{\theta }_n^{*}]}^{\mathrm{T}}$ can be expressed as

$${\mathsf{\Phi }}^{*}={ \mathrm{M}}_{\mathrm{ML}}^{-1}{\mathrm{W}}_{\mathrm{ML}}^{\mathrm{T}}\mathrm{Y}$$

(5)

where ${\mathrm{M}}_{\mathrm{ML}}={\sigma }_{ML}^2\mathrm{I} +{ \mathrm{W}}_{\mathrm{ML}}^{\mathrm{T}}\mathrm{W}$, ${\sigma }_{ML}^2$, and ${\mathrm{W}}_{\mathrm{ML}}$ are the maximum likelihood of ${\sigma }^2$ and $\mathrm{W}$. Equation (5) has the mean of ${\mathrm{M}}_{\mathrm{ML}}^{-1}{\mathrm{W}}_{\mathrm{ML}}^{\mathrm{T}}\mu$. The $\mathrm{d}$-dimensional principal components (PC) matrix ${\mathsf{\Phi }}^{*}$ is used in the following data-driven modeling.

2.3. MTGP Model

In the data analysis of turbomachinery, the lack of data is often caused by the damage of sensor equipment, the insufficient maintenance of equipment, the power failure of equipment, the aging fault of equipment hardware and software, the communication interruption of equipment, and new equipment recently put into operation. Data imputation and augmentation are crucial to the performance of the deep learning model for time-series data, which can effectively improve the size and quality of the training data. We explore a method for time-series data imputation using the multi-task GP model.

We assume a given dataset T = $\{\mathrm{x}$, $\mathrm{y}\}$, $\mathrm{x} =\left\{x_i\right| i = 1, \dots , n\}$ is specified as the index set, such as the observed time, and $\mathrm{y}=\left\{y_i\right| i = 1, \dots , n\}$ is the corresponding value of the observed training data. The purpose is to learn a regression model: $y = f\left(x\right) + \epsilon$, where $f\left(x\right)$ is a latent function and $\epsilon \sim N\left(0, {\sigma }^2\right)$ is a noise term. Given the test points ${\mathrm{x}}^{*}=\left\{x_i^{*}\right| i = 1, \dots , p\}$, we can predict the uncharted test observations ${\mathrm{y}}^{*}=\left\{y_i^{*}\right| i = 1, \dots , p\}$.

Gaussian processes (GP) are useful for time-series modeling [30] because they can accommodate irregular and uncertain observations, and provide probabilistic predictions. Gaussian process models suppose that $f\left(x\right)$ can be expressed as a probability distribution of the following functions:

$$\mathbf{y}=f(\mathbf{x})~GP(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime }))$$

(6)

where $m\left(\mathrm{x}\right)$ is the mean function and $k\left(\mathrm{x}, \mathrm{x}{}^{\prime }\right)$ is a covariance function describing the coupling between the inputs. The prior knowledge of the function behavior, which we want to model, can be encoded by modifying the covariance function. In the covariance of the vector $\mathrm{x}\in {\mathrm{R}}^n$, the size of the covariance matrix $\mathrm{K}\left(\mathrm{x}, \mathrm{x}\right)$ is $n\times n$, and the covariance function $k\left(x_i, x_j\right)$ provides the element ${\mathrm{K}}_{ij}$. The covariance function usually satisfies Mercer’s theorem; that is, $\mathrm{K}\left(\mathrm{x}, \mathrm{x}\right)$ must be symmetric and positive semidefinite. The affine transformation of the basic covariance function can be used to construct the complex covariance function.

For a given set T, we can obtain the predictions at the test indices ${\mathrm{x}}^{*}$ through calculating the conditional distribution $p({\mathrm{y}}^{*}|{\mathrm{x}}^{*}, \mathrm{x}, \mathrm{y})$, where $p({\mathrm{y}}^{*}|{\mathrm{x}}^{*}, \mathrm{x}, \mathrm{y})$ is a Gaussian distribution:

$$p({\mathbf{y}}^{\ast }|{\mathbf{x}}^{\ast },\mathbf{x},\mathbf{y})~N(m({\mathbf{y}}^{\ast }),{\sigma }^2({\mathbf{y}}^{\ast }))$$

(7)

where $m\left({\mathrm{y}}^{*}\right)$ is the mean and ${\sigma }^2\left({\mathrm{y}}^{*}\right)$ is the variance. Without the loss of generality, it is usually assumed that the mean function $m\left(\mathrm{x}\right)$ is zero. Thus, $m\left({\mathrm{y}}^{*}\right)$ and ${\sigma }^2\left({\mathrm{y}}^{*}\right)$ can be expressed as

$$m({\mathbf{y}}^{\ast })=\mathbf{K}{(\mathbf{x},{\mathbf{x}}^{\ast })}^{\mathrm{T}}\mathbf{K}{(\mathbf{x},\mathbf{x})}^{-1}\mathbf{y}$$

(8)

$${\sigma }^2({\mathbf{y}}^{\ast })=\mathbf{K}({\mathbf{x}}^{\ast },{\mathbf{x}}^{\ast })-\mathbf{K}{(\mathbf{x},{\mathbf{x}}^{\ast })}^{\mathrm{T}}\mathbf{K}{(\mathbf{x},\mathbf{x})}^{-1}\mathbf{K}(\mathbf{x},{\mathbf{x}}^{\ast })$$

(9)

Furthermore, we assume that $\mathrm{X} =\left\{x_i^j\right|j = 1, \dots , m, i = 1, \dots , n^j \}$ are the training indices, and $\mathrm{Y} =\left\{y_i^j\right| j = 1, \dots , m, i = 1, \dots , n^j \}$ are the observations for $m$ tasks, where task $j$ has $n^j$ training points. In an actual situation, we may encounter task-specific with $n^j$ training indices where the $n^j$ number is far lower than the other tasks. This is a challenge for the problem.

The problem of modeling $m$ tasks simultaneously (e.g., multiple time-series) drives the extension of the MTGP model, where all models have the same index set $\mathrm{x}$. Training an STGP (single-task GP) model for each task separately is a straightforward method, as shown in Figure 2a. For the modeling of $m$ tasks via MTGP, to reveal the affiliation of index $x_i^j$ and observation $y_i^j$ to task $j$, we have to add a label $l^j$ to the model as an additional input with $l^j=j$, as illustrated in Figure 2b.

In MTGP, two independent covariance functions, $k_c$ and $k_t$, can be assumed as

$$k_{MTGP}(x,x{}^{\prime },l,l{}^{\prime })=k_c(l,l{}^{\prime })\times k_t(x,x{}^{\prime })$$

(10)

where $k_c$ denotes the correlation between tasks and only depends on the labels $l$ and $k_t$ denotes the temporal covariance functions within a task and only depends on the indices $x$. Let $n^j = n$ for $j = 1, \dots , m$; the covariance matrix ${\mathrm{K}}_{\mathrm{MTGP}}$ can be expressed as

$${\mathbf{K}}_{MTGP}(\mathbf{X},\mathbf{l},{\mathit{\theta }}_c,{\mathit{\theta }}_t)={\mathbf{K}}_c(\mathbf{l},{\mathit{\theta }}_c)\otimes {\mathbf{K}}_t(\mathbf{X},{\mathit{\theta }}_t)$$

(11)

where $l = \{j | j = 1, \dots , m\}$; ⊗ is the Kronecker product; ${\theta }_c$ and ${\theta }_t$ are vectors, including hyperparameters for ${\mathrm{K}}_c$ and ${\mathrm{K}}_t$; ${\mathrm{K}}_c$ has a size of $m\times m$; and ${\mathrm{K}}_t$ has a size of $n\times n$. Thus, the matrix, ${\mathrm{K}}_{MTGP}$ has a size of $mn\times mn$. We named ${\mathrm{K}}_c$ as the correlation matrix. In practice, to adapt the model to the usage of task-specific training data, we can relax the assumption $n^j = n$ for $j = 1, \dots , m$.

This model has some useful features: we may have task-specific training indices $n^j$ (i.e., training data may be observed at task-specific times); when fitting the covariance function in (11), the correlation within the tasks can be learned automatically; the tasks have analogous temporal characteristics and hyperparameters ${\theta }_t$ in the framework assumption.

In summary, as MTGP is not a parametric approach, the number of parameters can grow as more observation data are collected. The advantage of the MTGP model over other comparable related methods, such as support vector regression (SVR), is that it may easily express prior knowledge of function behavior (e.g., periodicity or smoothness). Given multiple time series, the goal is to increase the overall modeling accuracy by employing a single MTGP model rather than several STGPs. MTGPs have been used in financial time series [31], environmental sensor networks [32], and the classification of ECG signals [33]. It is important to note that MTGPs can include nonuniformly sampled time series in a single model without further downsampling or interpolation. In addition, it is easy to include prior knowledge between time series, such as time shifts or assumptions about similar temporal. Thus, the method has many applications, such as prediction, correlation analysis, the estimation of time shifts between signals, and the compensation of missing data.

2.4. LSTM Model

It is very important to predict the operation trend for turbomachinery reliability. Recurrent neural networks (RNN) can predict trends quickly and accurately. However, the traditional RNN has the problem of gradient disappearance or explosion and cannot solve long-term dependence, which affects its practical application [34]. To solve this problem, we established a prediction model based on LSTM (long short-term memory). The LSTM is an improvement on RNN [35], which is composed of multiple connected recurrent units recursively. Three gate structures, i.e., forgetting gate, input gate, and output gate, constitute the recurrent unit of LSTM [23]. The advantages of LSTM are mainly embodied in the following aspects: (1) the three gating structures of LSTM can learn long-term memory information and provide a solution to long-term dependence; and (2) the sigmoid function and tank function are combined to form the activation function in LSTM, which keeps the gradient almost unchanged during the derivation of backpropagation, avoiding the disappearance or explosion of the gradient, and greatly accelerating the convergence speed of the model.

The LSTM neural network has a variety of input forms. Here, we consider the univariate time-series data. We suppose that there are $\mathrm{n}$ discrete time points with a uniform distribution, and $y=\{y_t | t =1, \dots , \mathrm{n}\}$ is a multi-dimensional time variable at each observed time point. These random variables serve as the input of the LSTM neural network.

According to Figure 3, at time $t$, LSTM receives three inputs: the network’s input $y_t$ in the present, the previous time output $h_{t-1}$, and the previous time cell state $C_{t-1}$. At the same time, LSTM has two outputs: output $h_t$ and cell state $C_t$. We can notice that $y_t$, $h_t$, and $C_t$ are all vectors. LSTM has three gate-like structures that filter through the information [36,37].

The forgetting gate controls how much information of the previous time cell state $C_{t-1}$ is passed to $C_t$ at the present time:

$$f_t=\sigma (W_f[h_{t-1},y_t]+b_f)$$

(12)

where $\sigma \left(·\right)$ is the sigmoid activation function, $W_f$ is the weight matrix of the forgetting gate, and $b_f$ is the bias.

The input gate controls how much information of input $y_t$ is passed to $C_t$ at the present time (i.e., it is used to update the cell state):

$$i_t=\sigma (W_i[h_{t-1},y_t]+b_i)$$

(13)

$${\tilde{C}}_t=\tanh (W_c[h_{t-1},y_t]+b_c)$$

(14)

$$C_t=f_t\ast C_{t-1}+i_t\ast {\tilde{C}}_t$$

(15)

where ${\mathrm{i}}_{\mathrm{t}}$ is the input gate, $W_i$ is the weight matrix of the input gate, and $b_i$ is the bias. ${\tilde{C}}_t$ is the input cell state, $W_c$ is the weight matrix of the cell state, and $b_c$ is the bias of the cell state.

The output gate determines how much information of $C_t$ is passed to the present output $h_t$:

$$o_t=\sigma (W_o[h_{t-1},y_t]+b_o)$$

(16)

$$h_t=o_t\ast \tanh (C_t)$$

(17)

where $o_t$ is the output gate, $W_o$ is the weight matrix of the output gate, and $b_o$ is the bias. The final output of the LSTM is determined by the output gate and the cell state.

3. Experiments and Results

In the following sections, we first illustrate the performance of the proposed MT-LSTM on an optical marker (OM) dataset. Thereafter, the superiority of MT-LSTM is further verified on a real turbine operation dataset.

3.1. OM Dataset

We first investigated whether the proposed MT-LSTM could augment scarce time-series data among multiple tasks and improve prediction accuracy. We used multivariate biomedical data from the literature [38]. These data consisted of three optical markers (OM) time series that were located at the midline of the body: the thorax (OM1), the lower breastbone (OM2), and near the navel (OM3). The three-dimensional positions of all sensors were simplified as their first principal component. In addition, the fourth task corresponded to a respiration belt (RB) placed around the torso near OM2. The monitoring time and sampling frequency (2.6 Hz) were the same for all tasks. Each task included 200 data points, as shown in Figure 4. We investigated two situations—a normal situation and an extreme situation—to prove the method’s validity. The normal situation included a few missing data (named S1) and different sampling frequencies (S2). We named the extreme situation S3. It is important to note that since the OM dataset was clean and low-dimensional, we did not conduct data pre-processing such as data denoising and dimensionality reduction (Figure 1).

In this paper, we used some indicators to evaluate prediction accuracy. The predicted test observations are denoted as ${\widehat{y}}^{*}=\left\{{\widehat{y}}_i^{*}\right| i = 1, \dots , n\}$. We used the mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE), defined as

$$MAE({\widehat{\mathit{y}}}^{\ast },{\mathit{y}}^{\ast })=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|{\widehat{y}}^{\ast }-y^{\ast }\right|}$$

(18)

$$RMSE({\widehat{\mathit{y}}}^{\ast },{\mathit{y}}^{\ast })=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\widehat{y}}^{\ast }-y^{\ast })}^2}}$$

(19)

$$MAPE({\widehat{\mathit{y}}}^{\ast },{\mathit{y}}^{\ast })=\frac{100\%}{n}{\displaystyle \sum _{i=1}^n\left|\frac{{\widehat{y}}^{*}-y^{*}}{y^{*}}\right|}$$

(20)

Lower MAE, RMSE, and MAPE indicate higher prediction accuracy.

3.1.1. Multi-Task Data Imputation

(1)

S1—Few data missing

For the S1 case, we assumed that a few data were missing for each of the four tasks (OM1-3, RB). Specifically, each task picked the first 150 data points and had 30 data points removed at different time points. Figure 5 and

Table 1. Error results of MTGP and STGP for S1.

Error of MTGP	Task 1	Task 2	Task 3	Task 4
MAE	0.194	0.087	0.160	0.240
RMSE	0.255	0.120	0.193	0.315
MAPE (%)	17.059	70.374	20.725	15.138
Error of STGP	Task 1	Task 2	Task 3	Task 4
MAE	1.626	0.726	0.191	1.762
RMSE	2.041	0.893	0.234	2.246
MAPE (%)	97.723	254.216	29.908	88.490

show the imputation and error results of each task using MTGP and STGP. Obviously, the MTGP model fit well with the original data of each task by transferring knowledge across tasks. The results were very accurate, and the related MAE and RMSE were close to 0. In contrast, for the results of STGP, it was impossible to observe the trend of the tasks, and the corresponding MAPE was much larger than that of MTGP.

(2)

S2—Different Sampling Frequencies

For the S2 case, we assumed that the sampling frequency of the tasks was different. In practice, the different sampling frequencies of various sensors led to different data volumes of signals, which introduced challenges to the subsequent modeling. We set the sampling frequency of Task 3 to be half that of the other tasks, i.e., Task 3 had only 75 data points. Tasks 1–2 and Task 4 had 30 data points removed (same as S1). As shown in Figure 6 and

Table 2. Error results of MTGP and STGP for S2.

Error of MTGP	Task 1	Task 2	Task 4
MAE	0.163	0.101	0.583
RMSE	0.224	0.136	0.720
MAPE (%)	15.127	78.622	68.081
Error of STGP	Task 1	Task 2	Task 4
MAE	1.745	0.317	1.840
RMSE	2.172	0.395	2.257
MAPE (%)	103.502	127.008	85.353

, the MTGP was found to effectively model these tasks with different sampling frequencies and have accurate predictions.

3.1.2. Time-Series Forecasting via Multi-Task Data Imputation

For the S3 case, we considered an extreme case—some tasks had very scarce data at finite early time points. For Task 1, we kept the first 13 points and removed the remaining 137 data points (91.3% data missing), which eliminated the periodicity of the data itself. The other tasks were kept unchanged. We predicted that the few-shot scenario would undoubtedly introduce great challenges to data augmentation and prediction.

First, we performed data imputation for Task 1 with STGP. As shown in Figure 7, the predictions quickly failed. Next, we used the MTGP model and observed significantly improved predictions (Figure 7). This was also verified in

Table 3. Error results of STGP and MTGP for S3.

Error	STGP	MTGP
MAE	1.84	0.494
RMSE	2.128	0.648
MAPE (%)	170.474	68.578

, which showed that in comparison to STGP, the RMSE of MTGP decreased sharply from 2.128 to 0.648. Consequently, the prediction results of the MT-LSTM in Figure 8 and

Table 4. Error results of MT-LSTM, ST-LSTM, and LSTM for S3.

Error	ST-LSTM	MT-LSTM	LSTM on Task 1
MAE	1.853	0.508	3.254
RMSE	2.187	0.71	3.704
MAPE (%)	156.440	37.002	392.669

were also more accurate predictions in comparison to the ST-LSTM based on STGP. Moreover, when the incomplete data (13 data points) were directly used for LSTM training, the predictions were extremely inaccurate (see Figure 8 and Table 4).

Hence, the validity of the proposed MT-LSTM was preliminarily verified on the OM dataset.

3.2. Real-World Datasets

The performance of our method was verified on a set of actual turbine operating data. In a large centrifugal compressor manufacturing enterprise, more than 10 compressors suffered accidents during operation in several years, many of which were new units recently put into operation, resulting in huge losses. Meanwhile, some accidents were accompanied by sensor failures. To analyze these accidents, we studied one of the compressor units. After the new compressor was brought into service, the vibration condition and the amplitude were stable at about 20 μm, which was much lower than the alarm value. However, the vibration value rose and exceeded the alarm value after more than 10 days of operation. It is obvious that the scarcity of available data poses significant challenges to this research. The original time-series data were obtained from various sensors positioned on the centrifugal compressor, and the deployment of the sensors is shown in Figure 9. The collected data contained various operating parameters: for example, inlet and exhaust pressure, exhaust flow, bearing temperature, and axis vibration. The sensors measured the operating data every hour, and a total of 11 performance variables were collected in the same time period in this paper. Hence, there were 11 tasks, each with only 200 data points. Obviously, these data contain noise.

Table 5. Details of tasks acquired from the centrifugal compressor.

Variable	Task	Minimum	Maximum
Inlet pressure	Task 1	0.099	0.103
Inlet temperature	Task 2	1.6	28.8
Exhaust pressure	Task 3	0.5	0.514
Exhaust temperature	Task 4	80	91.4
Exhaust flow	Task 5	90,086	99,872
Inlet guide vane	Task 6	30.2	63
Bearing 1’s temperature	Task 7	52.4	55.8
Bearing 2’s temperature	Task 8	75.7	78.4
Axis temperature	Task 9	73.3	79.3
Axis vibration	Task 10	19.4	26.3
Axis displacement	Task 11	−0.304	−0.266

shows the details of the tasks.

If we consider an extreme case such as S3 in the OM dataset, some tasks had very scarce data collected at the early time points. Specifically, each task picked the first 150 data points. For Task 10, we kept the first 20 data points and removed the remaining 130 data points (86.7% data missing). The other tasks were unchanged.

3.2.1. Results without Data Pre-Processing

The MT-LSTM method we propose in this paper was applied to this dataset. First, we used the MTGP model directly without data pre-processing. The calculation procedure was the same as that for the S3 of the OM dataset. Figure 10 and

Table 6. Error results of STGP and MTGP for the extreme situation.

Error	STGP	MTGP
MAE	0.597	0.534
RMSE	0.774	0.722
MAPE (%)	2.678	2.391

show the results. It was found that the imputation results of STGP failed quickly, and the prediction variance was high. As expected, the predictions of MTGP improved significantly over STGP. Consequently, the prediction of MT-LSTM achieved the best results (Figure 11 and

Table 7. Error results of MT-LSTM, ST-LSTM, and LSTM for the extreme situation.

Error	ST-LSTM	MT-LSTM	LSTM on Task 10
MAE	0.642	0.544	3.275
RMSE	0.791	0.686	3.648
MAPE (%)	6.418	5.646	13.994

).

3.2.2. Results with Data Pre-Processing

To investigate the impact of data pre-processing, we first normalized each task before MTGP modeling. The aim was to lessen the influence of various data value scales as normalization can compare multiple time-series variables at the same time. It prevents large-valued variables from overcontrolling small-valued variables, which improves accuracy. As shown in Figure 12 and Figure 13,

Table 8. Error results of STGP and MTGP (after normalization).

Error	STGP	MTGP
MAE	0.052	0.032
RMSE	0.066	0.042
MAPE (%)	31.938	19.518

and

Table 9. Error results of MT-LSTM, ST-LSTM, and LSTM (after normalization).

Error	ST-LSTM	MT-LSTM	LSTM on Task 10
MAE	0.056	0.036	0.222
RMSE	0.067	0.044	0.247
MAPE (%)	77.799	73.499	84.734

, the MTGP and MT-LSTM results were improved compared to those without data pre-processing.

Thereafter, we further performed data pre-processing, such as data denoising and dimensionality reduction; the first five principal components are shown in Figure 14. We found that partial principal components accounted for more information in the original data: the first principal component comprised 73.1%, the second principal component comprised 12.8%, the third principal component comprised 4%, the fourth principal component comprised 3.2%, and the fifth principal component comprised 2.8%. We normalized the first and second principal components (totally accounting for 85.9% of the information in the original data) and integrated them into the MT-LSTM model. The results are shown in Figure 15 and Figure 16,

Table 10. Error results of STGP and MTGP (with data pre-processing).

Error	STGP	MTGP
MAE	0.415	0.172
RMSE	0.554	0.218
MAPE (%)	34.557	20.218

and

Table 11. Error results of MT-LSTM, ST-LSTM, and LSTM (with data pre-processing).

Error	ST-LSTM	MT-LSTM	LSTM on Task 10
MAE	0.441	0.183	1.188
RMSE	0.566	0.228	1.424
MAPE (%)	70.872	67.051	99.277

. We found that the accuracy was significantly high, and the calculation time was about 5 min.

According to the experimental results, it could be found that the LSTM was poor in the scarce-data scenario. For the proposed MT-LSTM model, we performed data imputation of the target task by leveraging knowledge from auxiliary tasks through the modeling of an MTGP and then employed the LSTM to build the sequence model of the target data upon the imputed time-series data to perform a prediction in the future. The knowledge transfer across tasks helps the proposed MT-LSTM perform significantly better than the vanilla LSTM.

4. Conclusions and Discussion

Health management for turbomachinery is a challenging task. Deep learning has gradually entered this industry over the past few years. However, in practical application, the training process of deep learning models requires a large number of samples and time. Due to various reasons, it is difficult to obtain a significant amount of monitoring data, and the health management of turbomachinery has high timeliness and accuracy requirements. To account for the scarce-data scenario in the PHM of turbomachinery, this paper proposes integrating multi-task Gaussian processes, long short-term memory, Bayesian wavelet denoising, and probabilistic principal component analysis to achieve rapid automatic data imputation and time-series modeling to alleviate the need for big data as well as improve predictions with scarce data. Meanwhile, the framework can solve the optimization cleaning of the original multivariate time-series monitoring data, conduct the possible correlation and uncertainty in the multivariate data, and reduce the reliance on professional knowledge and experience. The experimental results show that the proposed method can accurately and efficiently reconstruct and predict with scarce time-series data via multi-task modeling. The superiority of our method is especially obvious for real centrifugal compressor operation data. When there are only 10~20 data points, the proposed method still has high prediction accuracy. This is significant to the prognostics and health management (PHM) of turbomachinery.

Due to the high precision requirement of turbomachinery operation, a slight error may lead to serious accidents; thus, the model’s accuracy should be improved. The main links of the model—data denoising, dimensionality reduction, and MT-LSTM—are processed separately, which is inconvenient for practical application. At present, the model has only been verified by an actual compressor, which cannot effectively investigate its robustness and universality. In fact, large engineering applications usually have hundreds of different types of turbines, and the working conditions are very complex. Therefore, in the future, we will carry out the robustness and universality verification of the model under multiple scenarios, such as multisource data fusion, different turbomachinery types, and different working conditions. In addition, we would like to continue development work in visualization and human–computer interaction in the future.