Fault Detection based Multiple Local Manifold Learning and Its Application to Blast Furnace Ironmaking Process

1. Introduction

Industrial processes are becoming more complex and their hazards to the environment are receiving increasing attention. Process safety is a non-negligible component of industrial processes, which comprises several steps such as hazard identification and analysis [1]. In particular, the identification and analysis of hazards is a key step in the prevention and mitigation of major process accidents.

In industrial processes, timely and accurately identifying abnormal operating conditions can prevent major accidents and improve operational efficiency, thus achieving compliance with environmental and safety regulations. Dynamic process monitoring for hazard/fault identification is already a trend in the future development of process safety and risk management [2]. Real-time monitoring of process operations to ensure safety measures is an essential step in the modern process industry [3].

Fault detection plays a pivotal role in guaranteeing operation safety and reducing downtime in complex industrial processes [4,5]. Broadly speaking, fault detection techniques can be categorized into three classes, model-based, knowledge-based, and data-driven based methods [6]. Model-based methods rely on the mathematical model. However, the mathematical model is often difficult or time-consuming to establish for complex industrial processes such as the blast furnace ironmaking process. For knowledge-based methods, the model is built from expert knowledge or qualitative information, which limits its applications for complex industrial processes. Conversely, only the measured process variables are required for data-driven fault detection methods. Thus, data-driven techniques are more suitable and efficient for the fault detection of complex industrial processes [7]. With the advance of the sensor, communication, and computing technologies, large amounts of data are collected in modern industrial processes. Under such circumstances, data-driven fault detection has gained an explosive amount of attention in recent years both from academia and industry [8].

To handle the highly correlated high-dimensional process data, multivariate analysis (MVA) has been widely employed in industrial processes [9]. In MVA, the process behavior is modeled by transforming the high-dimensional data into a lower-dimensional space. The features are extracted for establishing monitoring statistics. Among the MVA based fault detection methods, principal component analysis (PCA) has gained widespread popularity in process monitoring and fault diagnosis in the past decades [10]. In PCA, process data is projected into a lower-dimensional space to preserve the significant variability information as much as possible. Due to its efficiency and simplicity, PCA has been successfully applied in a large number of industrial processes [11]. Despite this, PCA is regarded as a kind of globality-based linear dimensionality reduction technique. However, the process data mostly lie on or close to a low-dimensional manifold. Compared to globality-based methods, manifold learning is an approach to nonlinear dimensionality reduction through discovering the manifold structure of data. In manifold learning, the input data is assumed to be sampled from a low-dimensional manifold. Representative manifold learning methods include Isomap [12], Locally Linear Embedding (LLE) [13], Laplacian eigenmaps (LE) [14], Local Tangent Space Alignment (LTSA) [15], Locality Preserving Projections (LPP) [16], Neighborhood preserving embedding (NPE) [17], and Hessian eigenmaps (or called Hessian LLE) [18].

In NPE, each data point is represented as a linear combination of the neighboring data points. Then, an optimal embedding is found to preserve the neighborhood structure in the dimensionality reduced space [17]. Chen et al. [19] applied eigenvalue decomposition and generalized eigenvalue decomposition to solve the unstable problem caused by singularity problem in NPE, and developed an NPE-based incipient fault detection method for small-scale cyber-physical systems. Since the NPE method can preserve the local manifold structure of different modes, Song et al. [20] performed NPE on the time-lagged variables for multimode dynamic process monitoring. LPP is designed to find the optimal linear approximations to the eigenfunctions of the Laplace Beltrami operator on the manifold through the nearest neighbor search in the low-dimensional space [16]. Duan et al. [21] employed LPP to preserve the local structure of process data, then adopted least squares support vector machine to predict the key-performance-indicator. Zhang et al. [22] combined LPP and PCA to preserve both global and local structures of the data set and developed fault detection and identification method by utilizing the extracted features. LLE attempts to discover nonlinear structure in high-dimensional data by exploiting the local symmetries of linear reconstructions [18]. Wu et al. learned structure information by LLE and incorporated the extracted local information into canonical correlation analysis (CCA) for quality-relevant nonlinear process monitoring [23]. Li and Zhang implemented the supervised locally linear embedding projection method for bearing fault diagnosis and illustrated its validity using the experimental data [24]. Different manifold learning methods focus on uncovering the manifold structure with different criteria. It relies on the knowledge and experience of experts for their own purposes. Therefore, only partial information from the underlying manifold is learned by each existing local manifold learning method. To take advantages of different manifold learning methods for better discovering the underlying manifold structure, Xing et al. [25] provided a common framework to synthesize the partial information extracted from different local manifold learning methods under local tangent coordinates.

Motivated by the above discussions, a novel data-driven fault detection based on fused local manifold learning (FLML) is proposed in this paper. In the proposed FLML, the partial information on the geometric structure of the underlying manifold is firstly extracted from LE, LLE, and Hessian Locally Linear Embedding (HLLE) methods, respectively. A novel objective function is formulated to fuse the extracted partial information. On the basis of the optimization results, FLML can learn the geometric information from different local methods. The geometric structure of the underlying manifold is more thoroughly explored by the proposed FLML, compared to LE, LLE, and HLLE. Like the PCA-based fault detection method, two monitoring statistics including Hotelling’s $T^2$ and Q statistics are established. The effectiveness and advantages of the proposed FLML-based fault detection are illustrated by an industrial Tennessee Eastman process benchmark and a real blast furnace ironmaking process.

The rest of this paper is organized as follows. Section II briefly introduces the ideas of LE, LLE, and HLLE. Section III illustrates the proposed FLML method and its application in fault detection in detail. In Section IV, the proposed FLML-based fault detection approach is verified through an industrial Tennessee Eastman (TE) process benchmark and a real blast furnace ironmaking process. Finally, Section V provides the conclusion.

2. Brief Review of LE, LLE and HLLE

2.1. LE

LE is a well-known manifold learning method to extract local structure features in the original sample space. The main idea behind LE is that the corresponding projections of neighboring points in the low-dimensional space should be close if the neighboring points in the high-dimensional space are close [14].

Given a set $\mathit{X}={[{\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_n]}^T\in {\mathbb{R}}^{n\times D}$ of n sample points in ${\mathbb{R}}^D$, where D is the number of variables. The goal of LE is to find a set of points $\mathit{Y}={[{\mathit{y}}_1,{\mathit{y}}_2,...,{\mathit{y}}_n]}^T\in {\mathbb{R}}^{n\times d}$ ($d\ll D$) in the low-dimensional space to have the neighbor relations between sample points $\mathit{X}$ in the high-dimensional space.

In LE, the first step is to construct the adjacency graph. k-nearest neighbors ($KNN$) is a widely used neighbors selection strategy, due to its simplicity. To model the neighborhood relations between sample points, an adjacency matrix $\mathit{S}\in {\mathbb{R}}^{n\times n}$ is employed, where each element $S_{ij}$ represents the neighborhood relations between ${\mathit{x}}_i$ and its neighbor ${\mathit{x}}_j$. Using the Gaussian heat kernel, the adjacency matrix $S_{ij}$ can be formed as follows,

$$S_{ij}=\left\{\begin{array}{ccc}\hfill & {\mathit{e}}^{-\parallel {\mathit{x}}_i-{\mathit{x}}_j{\parallel }^2/2{\delta }^2}\hfill & ,{\mathit{x}}_i\mathrm{and}{\mathit{x}}_j\in \mathrm{neighbors}\hfill \\ \hfill & 0\hfill & ,\mathrm{otherwise}\hfill \end{array}\right.$$

(1)

where $\delta$ is the kernel width.

The objective function of LE can be cast as,

$$\begin{array}{ccc}\underset{{\mathit{y}}_i,{\mathit{y}}_j}{min}\hfill & J_{LE}=\sum _{i,j=1}^n{\parallel {\mathit{y}}_i-{\mathit{y}}_j\parallel }^2{S}_{ij}={\mathit{Y}}^T{\mathcal{L}}_e\mathit{Y}\hfill & \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathit{Y}}^T\mathit{D}\mathit{Y}=\mathit{I}\hfill \end{array}$$

(2)

where $\mathit{D}({\mathit{D}}_{ii}={\sum }_{j=1}^n{S}_{ij})$ is a diagonal matrix, and ${\mathcal{L}}_e({\mathcal{L}}_e=\mathit{D}-\mathit{S})$ is defined as a Laplacian matrix.

The optimization problem (2) is equivalent to generalized eigenvalue problem,

$${\mathcal{L}}_e\mathit{Y}=\lambda \mathit{D}\mathit{Y}$$

(3)

It can be readily solved through eigenvalue decomposition.

2.2. LLE

In LLE, the key assumption is that each data point and its neighbors are lied on or closed to a locally linear patch. A sample can be represented as the linear combination of multiple samples from its neighborhood [13]. The optimization problem of LLE can be formulated as,

$$\begin{array}{ccc}\underset{{\theta }_{ij}}{min}\hfill & J_{LLE}=\sum _{i=1}^n{\parallel {\mathit{x}}_i-\sum _{j=1}^n{\theta }_{ij}{\mathit{x}}_j\parallel }^2\hfill & \\ \mathrm{s}.\mathrm{t}.\hfill & \sum _{j=1}^n{\theta }_{ij}=1\hfill \end{array}$$

(4)

where ${\theta }_{ij}$ is the weight coefficient and only when ${\mathit{x}}_i$ and ${\mathit{x}}_j$ are neighbors the corresponding ${\theta }_{ij}$ has a value, otherwise it is 0. Similarly, the neighbors selection can be determined by $KNN$ method.

Denoted the low-dimensional features as $\mathsf{\Phi }={[{\mathit{\varphi }}_1^T,{\mathit{\varphi }}_2^T,...,{\mathit{\varphi }}_n^T]}^T$, the minimization problem (4) can be formulated as,

$$\begin{array}{cc}\underset{w_{ij}}{min}\hfill & J_{LLE}=\sum _{i=1}^n{\parallel {\mathit{\varphi }}_i-\sum _{j=1}^n{\theta }_{ij}{\mathit{\varphi }}_j\parallel }^2=tr\left({\mathsf{\Phi }}^T{\mathcal{L}}_l\mathsf{\Phi }\right)\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathsf{\Phi }}^T\mathsf{\Phi }=\mathit{I}\hfill \end{array}$$

(5)

where $\mathit{I}$ is the identity matrix, $\mathsf{\Theta }={\left[{\theta }_{ij}\right]}_{i=1,...,n,j=1,...,n}$ is the weight coefficient matrix and ${\mathcal{L}}_l={(\mathit{I}-\mathsf{\Theta })}^T(\mathit{I}-\mathsf{\Theta })$.

Similarly, the optimization problem (5) can be solved through eigenvalue decomposition.

2.3. HLLE

HLLE is regarded as a variant of LLE. Assumed that the low-dimensional data representation is locally isometric to an open and connected subset, the idea behind HLLE is to minimize the curviness of the high-dimensional manifold while embedding it into a low-dimensional space [18].

Assumed that the set $\mathit{X}={[{\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_n]}^T\in {\mathbb{R}}^{n\times D}$ is located on a smooth manifold $\mathcal{M}\subset {\mathbb{R}}^D$ with an intrinsic dimension $d\ll D$, PCA is firstly performed on $\mathit{X}$ to obtain d eigenvectors $\mathit{V}=[{\mathit{v}}_1,{\mathit{v}}_2,...,{\mathit{v}}_d]\in {\mathbb{R}}^{D\times d}$.

It is assumed that the number of neighbors of the sample ${\mathit{x}}_i$ is k. Then, the local tangent coordinate ${\mathit{u}}_j,j=1,...,k$ of the sample ${\mathit{x}}_i$ can be calculated by projecting the local neighborhood into the tangent subspace,

$${\mathit{u}}_j^i=\langle \mathit{V},{\mathit{x}}_j-{\mathit{x}}_i\rangle ={\mathit{V}}^T({\mathit{x}}_j-{\mathit{x}}_i)$$

(6)

Then, the Hessian matrix ${\mathit{H}}^i$ containing local information will be calculated by the projection of neighbors ${\mathit{x}}_j$ in the tangent coordinate ${\mathit{u}}_j$. Through defining ${\mathit{U}}^i=[{\mathit{U}}_1^i,{\mathit{U}}_2^i,...,{\mathit{U}}_k^i]$ constructing by the ${\mathit{u}}_j^i$[26], ${\mathit{H}}^i$ can be constructed by the last ${\left({\mathit{U}}^i\right)}^{†}$ where ${(.)}^{†}$ is the pseudo-inverse symbol. Therefore, the local objective of HLLE can be estimated with,

$$\begin{array}{ccc}\underset{{\mathit{f}}^i}{min}{J}_{HLLE}\hfill & =\frac{1}{N}\sum _{i=1}^N{\mathit{f}}^{iT}{\mathit{H}}^{iT}{\mathit{H}}^i{\mathit{f}}^i\hfill & \hfill \\ \mathrm{s}.\mathrm{t}.{\mathit{f}}^{iT}{\mathit{f}}^i=1\hfill \end{array}$$

(7)

where ${\mathit{f}}^i={[f\left({\mathit{x}}_1\right),f\left({\mathit{x}}_2\right),...,f\left({\mathit{x}}_k\right)]}^T$ is the local projection and $f\left({\mathit{x}}_j\right),j=1,...,k$ is the smooth function used to estimate the local neighborhood information at a fixed point ${\mathit{x}}_i$. Then, extending the local neighborhood information to the all samples, we can get the global projection $\mathit{f}={[{\mathit{f}}^1,{\mathit{f}}^2,...,{\mathit{f}}^N]}^T$. The optimization objective of HLLE can be rewritten as,

$$\begin{array}{ccc}\underset{\mathit{f}}{min}{J}_{HLLE}\hfill & =\frac{1}{N}\sum _{i=1}^N{\mathit{f}}^T{\mathit{S}}^{iT}{\mathit{H}}^{iT}{\mathit{H}}^i{\mathit{S}}^i\mathit{f}\hfill & \hfill \\ & =tr\left({\mathit{f}}^T{\mathcal{L}}_h\mathit{f}\right)\hfill \\ \mathrm{s}.\mathrm{t}.{\mathit{f}}^T\mathit{f}=1\hfill \end{array}$$

(8)

where ${\mathcal{L}}_h=\frac{1}{N}{\sum }_{i=1}^N{\mathit{S}}^{iT}{\mathit{H}}^{iT}{\mathit{H}}^i{\mathit{S}}^i$ is the local objective and the neighboring selection matrix ${\mathit{S}}^i$ can convert local projection to global projection ${\mathit{f}}^i={\mathit{S}}^i\mathit{f}$ and each entry $S_{ij}^i$ of ${\mathit{S}}^i$ can be obtained as,

$$S_{ij}^i=\left\{\begin{array}{ccc}\hfill & 1,\hfill & {\mathit{x}}_i\mathrm{and}{\mathit{x}}_j\in \mathrm{neighbors}\hfill \\ \hfill & 0,\hfill & \mathrm{Otherwise}\hfill \end{array}\right.$$

(9)

For solving the above optimization problem (8), eigenvalue decomposition is also utilized. Details of HLLE can be found in [18].

3. Proposed Method

3.1. FLML: Fused Local Manifold Learning

From the previous section, it can be found that LE, LLE, and HLLE can explore the local geometric structure information from different perspectives. To synthesize these local information, Xing et al. [25] fused the local information obtained by multiple manifold learning methods including LE, LLE, HLLE and LTSA by reformulating the different local manifold algorithms under the local tangent coordinate system to reveal the underlying manifold of the dataset. Similar to [25], the fused local objective ${\mathcal{F}}_f$ is defined as

$${\mathcal{F}}_f=c_1{\mathcal{L}}_e+c_2{\mathcal{L}}_l+(1-c_1-c_2){\mathcal{L}}_h$$

(10)

to integrate LE, LLE, and HLLE. Here, $c_1+c_2\le 1,c_1,c_2\in [0,1]$ are the fusion coefficients.

Remark 1. In [25], the local objectives are optimized to determine the fusion coefficients and global embedding coordinates simultaneously. In [27], the selection of c is determined by employing the alternating optimization method which iterative updates c and $\mathit{Y}$ in an alternating fashion. For simplicity, only the global embedding coordinates are obtained in this study. The fusion coefficients are considered fixed values.

Similar to LE, LLE, and HLLE, the proposed FLML method is a kind of linear projection method. To facilitate the online fault detection, an explicit linear mapping from the original space to the low-dimensional space is provided. Thus, the main goal of FLML is to seek a transformation matrix that maps the high-dimensional data to low-dimensional data. Supposed that $\mathit{w}\in {\mathbb{R}}^{D\times 1}$ is the transformation vector from $\mathit{X}$ to $\mathit{y}$, therefore, the projections in the low-dimensional space can be represented as $\mathit{y}={\mathit{w}}^T\mathit{X}$.

To synthesize the local geometric structure information and impose constraints to prevent multiple solutions, the objective function of FLML is formulated as follows,

$$\begin{array}{ccc}\underset{\mathit{w}}{min}\hfill & J_{FLML}={\mathit{w}}^T{\mathit{X}}^T{\mathcal{F}}_f\mathit{X}\mathit{w}\hfill & \\ \mathrm{s}.\mathrm{t}.\hfill & {\mathit{w}}^T{\mathit{X}}^T\mathit{X}\mathit{w}=1\hfill \end{array}$$

(11)

Remark 2. As shown in (10) and (11), it is noticed that the corresponding local structure will be extracted from LE, LLE and HLLE while $c_1$ and $c_2$ are selected to be $(c_1=1,c_2=0)$, $(c_1=0,c_2=1)$ and $(c_1=0,c_2=0)$, respectively.

To solve the optimization problem (11), we use the technique of Lagrange multipliers as follows,

$$J_{\mathrm{FLML}}={\mathit{w}}^T{\mathit{X}}^T{\mathcal{F}}_f\mathit{X}\mathit{w}+\lambda /2(1-{\mathit{w}}^T{\mathit{X}}^T\mathit{X}\mathit{w})$$

(12)

where $\lambda$ is the Lagrange multiplier.

While $\frac{\partial J_{\mathrm{FLML}}}{\partial \mathit{w}}=0$, it results in

$${\mathit{X}}^T{\mathcal{F}}_f\mathit{X}\mathit{w}=\lambda {\mathit{X}}^T\mathit{X}\mathit{w}$$

(13)

Hence, we can use generalized eigenvalue decomposition to obtain the transformation vector $\mathit{w}$ from (13). Finally, the transformation matrix $\mathit{W}\in {\mathbb{R}}^{D\times d}$ can be assembled by the eigenvectors corresponding to the smallest d eigenvalues derived from the result of generalized eigenvalue decomposition.

3.2. FLML based Fault Detection

Generally, data-driven fault detection methods contain two steps, offline modeling and online monitoring. In the offline modeling step, the process data $\mathit{X}\in {\mathbb{R}}^{n\times D}$ are collected under normal operating condition for training. Using the training data $\mathit{X}$, an FLML model is established. Subsequently, the low-dimensional data $\mathit{Y}\in {\mathbb{R}}^{n\times d}$ is obtained. With the transformation matrix $\mathit{W}$, the relation between $\mathit{X}$ and $\mathit{Y}$ can be represented,

$$\mathit{Y}=\mathit{X}\mathit{W}$$

(14)

The residual matrix $\mathit{E}\in {\mathbb{R}}^{n\times D}$ is,

$$\mathit{E}=\mathit{X}-\mathit{Y}{\mathit{W}}^T=\mathit{X}(\mathit{I}-\mathit{W}{\mathit{W}}^T)$$

(15)

In the online monitoring step, the low-dimensional data point ${\mathit{y}}_{new}$ and residual ${\mathit{e}}_{new}$ of a new standardized sample ${\mathit{x}}_{new}$ is obtained,

$$\left\{\begin{array}{cc}\hfill & {\mathit{y}}_{new}={\mathit{W}}^T{\mathit{x}}_{new}\hfill \\ \hfill & {\mathit{e}}_{new}=(\mathit{I}-\mathit{W}{\mathit{W}}^T){\mathit{x}}_{new}\hfill \end{array}\right.$$

(16)

The Hotelling’s $T^2$ and squared prediction error (SPE) statistics (also called Q statistic) are often used for monitoring. For $T^2$ statistic, the Mahalanobis distance is used to evaluate the variations of $\mathit{y}$. The Hotelling’s $T^2$ statistic is defined as,

$$T^2={\mathit{y}}_{new}^T{\mathsf{\Lambda }}^{-1}{\mathit{y}}_{new}$$

(17)

where $\mathsf{\Lambda }$ is the covariance matrix of $\mathit{Y}$ which is extracted from the normal operating condition data $\mathit{X}$.

For Q statistic, the Euclidean distance is adopted to evaluate the magnitude of vector in the residual space as follows,

$$Q={\mathit{e}}_{new}^T{\mathit{e}}_{new}$$

(18)

Under the assumption that all operating parameters and prediction errors have a Gaussian distribution, the upper control limits (UCLs) of Hotelling’s $T^2$ and Q statistics are determined by means of F distribution and ${\chi }^2$ distribution, respectively. Thus, with a level of significance $\alpha$, the UCLs are calculated as follows,

$$\left\{\begin{array}{ccc}\hfill & T_{UCL}^2\hfill & \sim \frac{d(n-1)}{n-d}F(d,n-d,\alpha )\hfill \\ \hfill & Q_{UCL}\hfill & \sim g{\chi }_{h,\alpha }^2\hfill \end{array}\right.$$

(19)

where $g={\sigma }^2/2\mu$ and $h=2{\mu }^2/{\sigma }^2$, $\mu$ is the sample mean and ${\sigma }^2$ is the sample variance of the statistic Q. g and h are calculated using normal operating condition data.

In FLML, the parameters such as the number of the neighbors k, the bandwidth $\delta$ of Gaussian heat kernel, the fusion weights $c_1$ and $c_2$, and the number of latent variables d are determined as follows,

• k: The finding of neighborhood relation is related to the selection of k [23]. To balance the computation complexity and generalization capability, we choose k in the range $[1:10]$ with the smallest mean false alarm rate (FAR). The definition of FAR can be found in next section.
• $\delta$: If the bandwidth of the Gaussian heat kernel function is too small, the kernel will be sensitive to noise. A large bandwidth may create an overly smooth mapping [28]. Empirically, the bandwidth $\delta$ is chosen as $\delta =bm{\sigma }^2$ where m is the size of variables, b is a constant, and ${\sigma }^2$ represents the variance of the data, which is 1 as the original data is normalized [22]. In the case studies, $b=50$ are selected.
• $c_1$ and $c_2$: It is noticeable that the hyper-parameters $c_1$ and $c_2$ have a important influence on the performance of the proposed FLML method. However, it is a challenging work to choose a set of optimal hyper-parameters. As a traditional way of performing hyper-parameter optimization, the grid search method is employed. For this purpose, a finite set of $c_1\in \{0,0.05,...,0.95\}$ and $c_2\in \{0,0.05,...,0.95\}$ are explored by minimizing the mean FAR.
• d: Similar to NPE-based and LPP-based methods, the number of latent variables d is selected by searching for eigenvalues similar to the smallest non-zero eigenvalue.

In summary, the procedure of the FLML based fault detection is described as follows,

Offline training:

• Step 1. Standardize the collected dataset $\mathit{X}$ under normal operating condition.
• Step 2. Extract the local structure information by calculating ${\mathcal{L}}_e$, ${\mathcal{L}}_l$ and ${\mathcal{L}}_h$ from (2), (5), and (8).
• Step 3. Establish the fused local structure $\mathcal{F}$ defined in (10).
• Step 4. Perform eigenvalue decomposition as (13) to solve the optimization problem (11) and retain the eigenvectors corresponding to the smallest d eigenvalues as the transformation matrix $\mathit{W}$.
• Step 5. Compute the UCLs $T_{UCL}^2$ and $Q_{UCL}$ using (19).

Online monitoring:

• Step 1. Standardize a new collected sample ${\mathit{x}}_{new}$.
• Step 2. Obtain the score ${\mathit{y}}_{new}$ and residual ${\mathit{e}}_{new}$ from (16).
• Step 3. Compute the online statistics $T^2$ and Q according to (17) and (18), respectively.
• Step 4. If any of $T^2$ and Q statistics violate the confident limit, a fault is declared.

4. Case Studies

In this section, the proposed FLML-based fault detection method is verified by conducting experiments on the Tennessee Eastman process and a real-world blast furnace ironmaking process.

4.1. Tennessee Eastman Process

The TE process is a well-known and widely used industrial benchmark for comparing the performance of process monitoring and control [29]. In the TE process, there are five main units including a separator, a compressor, a reactor, a vapour/liquid separator, a stripper and a condenser, and eight components A-H. It also has 12 manipulated variables (XMV1-11) and 41 measured variables (XMEAS1-41). Among these variables, XMEAS(23-41) are the composition of A-H measured with 6 min sampling intervals in different positions. Other variables are collected with 3 min sampling intervals. There are 21 preprogrammed faults in TE process. Of 21 faults, due to the absence of observable change in the mean and standard deviation between their corresponding faulty and normal operation, faults 3, 9 and 15 are very difficult to detect [30]. In this study, faults 3,9 and 15 are ignored. The descriptions of the faults considered in this study are listed in Table 1.

A widely used dataset of TE process can be found in http://web.mit.edu/braatzgroup/links.html. In this study, we also adopt this dataset. Specifically, 22 measured variables XMEAS(1:22) and 11 manipulated variables XMV(1:11) are selected as $\mathbf{x}\in {\mathbb{R}}^{33}$. 500 samples collected under the normal operating condition (IDV(0)) are used as a training dataset. For each faulty dataset, it has 960 samples in total where the fault is injected from the 160^th sample.

To assess the fault detection performance, three indices including Missed detection rate (MDR), Detection delay (DD) and false detection rate (FAR) are used. DD is defined as the time interval from the start of fault to the detection time, which is expressed as the first time of 5 consecutive rises. MDR and FAR with a $99.9\%$ confidence level can be calculated as follows,

$$\begin{array}{c}\mathrm{MDR}(\%)=\frac{\mathrm{No}.\mathrm{of}\mathrm{samples}(J\le J_{UCL}|\mathrm{Fault})}{\mathrm{total}\mathrm{samples}\left(\mathrm{Fault}\right)}\times 100\%\hfill \\ \mathrm{FAR}(\%)=\frac{\mathrm{No}.\mathrm{of}\mathrm{samples}(J>J_{UCL}|\mathrm{Normal})}{\mathrm{total}\mathrm{samples}\left(\mathrm{Normal}\right)}\times 100\%\hfill \end{array}$$

where $J\in \{T^2,Q\}$ and $J_{UCL}\in \{T_{UCL}^2,Q_{UCL}\}$.

For comparative study, several typical fault detection methods including PCA, NPE, LPP, Principal Component Pursuit (PCP)[31], kernel PCA (KPCA)[32], mixed KPCA (MKPCA) and structured joint sparse PCA (SJSPCA)[33] are employed. For PCA, $d=19$ is set according to the $95\%$ cumulative percentage of variance (CPV). For the NPE, $k=5$ and $d=19$ is selected. For LPP, $k=5$, $\delta =1650$, and $d=19$ are selected. For PCP, $d=19$ is selected. For KPCA, the kernel width $C=1000$ and $d=22$ is selected. For MKPCA, the kernel widths $C_1$, $C_2$ are set as 500 and 1000 respectively, and $d=28$ is selected. For SJSPCA, $k=5$, ${\lambda }_1=600$, ${\lambda }_2=50$, and $d=19$ are selected. As displayed in Figure 1, $k=5$ is selected for the TE case. Similarly, As shown in Figure 2, $c_1=0.25$ and $c_2=0.25$ are chosen. For FLML, $\delta =1650$, $c_1=0.25$, $c_2=0.25$, and $d=19$ are selected.

Table 2 and Table 3 list the MDRs and FARs of 18 faults. It can be observed that the NPE and LPP can provide superior performance over PCA. The utilization of local structure information can enhance fault detection performance. On the other hand, FLML $T^2$ statistic offers the lowest MDRs among all the comparative methods. The average of MDRs of FLML $T^2$ statistic is $11.49\%$. Compared to NPE and LPP methods, the average of MDRs of FLML $T^2$ statistic is increasingly reduced. Furthermore, the FARs of FLML $T^2$ statistic reach the same level of NPE and LPP. The average of FARs of FLML $T^2$ statistic is $0.45\%$. The DD index represents the sensitivity of monitoring statistics. A smaller DD means the monitoring statistic can detect the fault earlier. Table 4 lists the DDs for all methods, where the DD is indicated in the unit of the hour. As displayed in Table 4, FLML $T^2$ statistic can derive the smallest DD among the all methods.

To further illustrate the superiority of the proposed FLML method, the monitoring results are depicted in Figure 3 and Figure 4 by PCA, NPE, LPP, PCP, KPCA, MKPCA, SJSPCA, and FLML methods for IDV (10) and IDV (19), respectively. Fault 10 is designed by adding a random disturbance on the C feed temperature (Stream 4). Compared to the step type fault scenarios such as fault 1, fault 2, fault 10 is more difficult to detect. As shown in Figure 3, it is observed that FLML $T^2$ can detect most of the faulty samples. In contrast, other monitoring statistics such as PCA $T^2$, NPE $T^2$, PCP $T^2$, KPCA $T^2$, MKPCA $T^2$, SJSPCA $T^2$, and LPP $T^2$ can not effectively detect the fault 10, since most of the monitoring statistics are below the corresponding UCLs. Fault 19 is an unknown faulty type. The monitoring results of fault 19 is displayed in Figure 4. By comparison, it can be observed that PCA, NPE, LPP, PCP, KPCA, MKPCA, and SJSPCA fail to detect the occurrence of fault 19. The online statistics are almost lower than the UCLs after 161^th sample. A promising result can be obtained using FLML $T^2$ as shown in Figure 4(h).

4.2. Blast Furnace Ironmaking Process

In this section, the proposed method is verified through a real blast furnace (BF) ironmaking process at a steel company in South China. In the steel manufacturing process, the blast furnace ironmaking process plays a vital role. The blast furnace ironmaking process is considered one of the most complex industrial processes. The basic units of the blast furnace ironmaking process are depicted in Figure 5. As it can be seen from Figure 5, the ironmaking process can be mainly divided into several sub-systems including the charging system, gas processing system, hot air system, the pulverized coal injection system, the iron system, and the BF body. The inner structure of the blast furnace is vertical. In the BF, the iron ore and coke are fed from the top along the vertical direction. The ${1000}^{\circ }C$ hot air and coal powder are blown into the furnace from the bottom. Through complex chemical reactions, the molten iron and slag are generated and accumulated in the hearth. In a periodical way, the molten iron and slag are discharged from the bottom of the furnace through the tap hole. As a byproduct, the flux gas escapes from the top of the furnace.

In the ironmaking process, it often suffers from abnormal furnace conditions, due to the effects of unreasonable daily operation and various disturbances. If these abnormal furnace conditions are not timely detected, the product quality will be degraded, and even the safety of the plant would not be ensured. Thus, the effective detection of abnormalities becomes an indispensable component of the operation of a blast furnace. In this study, we consider the detection of the channeling fault. The channeling fault may be caused by several reasons such as the low-quality coke pulverized coal, the inappropriate adjustment of air volume. In a channeling accident, the high-temperature furnace gas passes in the path of least resistance at a high velocity. The furnace gas increases the heat load at the wall and top of the furnace, resulting in possible equipment damage such as burning of the bag dust catcher [34].

Figure 5. Diagram of Blast Furnace Ironmaking Process [35].

Figure 5. Diagram of Blast Furnace Ironmaking Process [35].

In this case, a data set was collected from 2013 12/21 to 2014 01/05. The sampling interval is 10s. As demonstrated in [36], we select 7 variables that are the most sensitive to faults. These variables are listed in Table 5. A total of 1000 samples are collected from 12/21 00:20 to 20-12/28 09:40 under the normal operating conditions as the training dataset. According to the accident report of the operation personnel, from 2013 12/30 23: 44 to 2013 12/31 05: 18, the ironworks accident report has recorded the occurrence of the channeling. To facilitate the verification of fault detection performance, a testing dataset are generated, where 1200 samples were collected where the fault occurred from the 200^th sample.

For PCA, $d=3$ is selected. For NPE, the neighbor parameter k and number of principle components d are selected to be 5 and 3, respectively. For LPP, $k=5$, $d=3$ and $\delta =7\ast 50$ are determined. For PCP, $d=4$ is selected. For KPCA, $C=200$ and $d=3$. For MKPCA, $C_1=150$, $C_2=200$ and $d=3$. For SJSPCA, $k=5$, ${\lambda }_1=600$, ${\lambda }_2=50$, and $d=3$. For fair comparison, $k=5$, $d=3$, $\delta =7\ast 50$ are selected for FLML. The result of grid search for selecting $c_1$ and $c_2$ is plotted in Figure 6. On the basis of grid search results, $c_1=0.3$ and $c_2=0.3$ are determined.

The monitoring results are plotted in Figure 7. It can be observed that NPE and LPP $T^2$ statistics can detect more faulty samples than PCA $T^2$ statistic. But, the improvement is limited as shown in Figure 7. The reason may be that NPE and LPP only extract partial local structure information. It is also noticed that there are more false alarms for PCA Q and PCP Q statistics. Compared to the other methods, the proposed FLML method can provide much better performance as shown in Figure 7(h). Only a few faulty samples are missed. Table 6 lists the MDRs, FARs, and DDs. As shown in Table 6, the MDR and FAR of FLML $T^2$ are $2.80\%$ and $0.00\%$, respectively. The DD of FLML $T^2$ is $0.17$ mins. The channeling condition can be timely and accurately detected by FLML $T^2$. Among these comparative methods, the proposed FLML method can offer the best monitoring performance.

5. Conclusions

In this paper, a novel data-driven fault detection technique based on fusing local manifold learning methods is proposed for complex industrial processes. The proposed method aims to explore a more comprehensive local structure by synthesizing partial information learned from LE, LLE, and HLLE. With the exploit of local geometric structure, the process data is projected into a lower-dimensional space. Hotelling’s $T^2$ and Q statistics are employed for fault detection. Case studies on the widely used TE process benchmark and a real blast furnace ironmaking process show the superior monitoring performance of the proposed methods, by comparison with other related methods. Future work will investigate the application of the proposed method on fault identification and classification.