A new process monitoring method based on noisy time structure independent component analysis☆

Lianfang Cai,Xuemin Tian*

College of Information and Control Engineering,China University of Petroleum,Qingdao 266580,China

Keywords:Process monitoring Independent component analysis Measurement noises Kurtosis Mixing matrix Contribution plot Sensitivity analysis

ABSTRACT Conventional process monitoring method based on fast independent component analysis(FastICA)cannot take the ubiquitous measurement noises into account and may exhibit degraded monitoring performance under the adverse effects of the measurement noises.In this paper,a new process monitoring approach based on noisy time structure ICA(NoisyTSICA)is proposed to solve such problem.A NoisyTSICA algorithm which can consider the measurement noises explicitly is firstly developed to estimate the mixing matrix and extract the independent components(ICs).Subsequently,a monitoring statistic is built to detect process faults on the basis of the recursive kurtosis estimations of the dominant ICs.Lastly,a contribution plot for the monitoring statistic is constructed to identify the fault variables based on the sensitivity analysis.Simulation studies on the continuous stirred tank reactor system demonstrate that the proposed Noisy TSICA-based monitoring method outperforms the conventional Fast ICA-based monitoring method.

1.Introduction

As modern industrial processes become more and more complex,process monitoring plays a key role in ensuring the process safety,product quality and economic profit.With a large number of process variables measured in industrial processes,multivariate statistical process monitoring(MSPM)techniques[1-4]which are a popular kind of data-driven methods are developing rapidly for the process fault detection and identification.The MSPM methods transform highdimensional measured data into low-dimensional space to extract meaningful information for detecting and identifying abnormal situations of industrial processes.Among the numerous MSPM methods,principal component analysis(PCA)is a classical approach.It extracts uncorrelated latent variables called principle components(PCs)which capture the most variances of the original variables and has many extended versions for fault detection and identification[5-8].However,PCA can only consider up to the second-order variance-covariance statistic and cannot make use of higher-order statistical information contained in the process data[9],which may lead to inadequate feature extraction and degraded monitoring performance.Furthermore,PCA makes a strict assumption that the extracted PCs follow multivariate Gaussian distributions to determine the control limits of the monitoring statistics.In most cases,the essential latent variables of industrial processes obey non-Gaussian distributions[10],and thus,there is a great possibility that this assumption cannot be in accordance with realistic industrial situations.

Recently,a newly emerging MSPM approach called independent component analysis(ICA)is attracting more and more attention from both the academic researchers and process engineers.Different from PCA,ICA can further utilize the higher-order negentropy statistic[11]or the second-order time-delayed covariance statistic[12]to recover mutually independent non-Gaussian latent variables called independent components(ICs)from the original measured variables and can be applied to deal with non-Gaussian processes which are more practical in the real-world manufacturing environment[13].Kano et al.[14]proposed an ICA-based univariate statistical process monitoring(USPM)method and demonstrated its superiority over the PCA-based MSPM.But the number of ICs'monitoring charts is comparatively large,which would raise the burden on the process monitoring.To solve such problem,Lee et al.[15]constructed three Mahalanobis-type monitoring statistics based on the extracted ICs and thus developed an ICA-based MSPM method.On the basis of this work,various improved versions of ICA were proposed by taking different process characteristics into account.Stefatos and Hamza[16]proposed a dynamic independent component analysis(DICA)method for capturing process dynamic pattern.Odiowei and Cao[17]combined canonical variate analysis(CVA)with ICA and developed a state-space ICA method for dynamic process monitoring.Tian et al.[18]proposed a multiway kernel independent component analysis(KICA)method for monitoring nonlinear batch process.Considering the extraction of non-Gaussian feature and the preservation of the local neighborhood structure simultaneously,Cai et al.[19]integrated KICA with locality preserving projections(LPP)for nonlinear process monitoring.To accountfor the process multimodal characteristic,Zhang et al.[20]developed a modified KICA based process monitoring method by the introduction of the Kronecker product.Rashid and Yu[21]proposed a hidden Markov model based adaptive ICA approach for monitoring processes with multimodality.In these studies,the well-known fast ICA algorithm(FastICA)[22]was adopted as a promising feature extraction technique.However,these Fast ICA-based process monitoring methods commonly used noise-free ICA model and thus could not take the influence of measurement noises into consideration.In fact,measurement noises always exist in industrial processes and thus inevitably contaminate measured data[23,24].Under the adverse effects of measurement noises,the monitoring performances of these Fast ICA-based monitoring methods may decline drastically.Consequently,developing an improved ICA-based monitoring method which can eliminate or attenuate the effects of measurement noises is extremely important for enhancing the process monitoring performance.To the best of our knowledge,currently,there is one monitoring method based on probabilistic ICA(PICA)[24]that can consider the measurement noises explicitly and can obtain the ICs without noises from the normal offline data.Nevertheless,the ICs calculated online are still corrupted by noises and thus cannot effectively reflect the real process information.

Motivated by the above analysis,a new process monitoring method based on noisy time structure ICA(NoisyTSICA)is proposed in this paper.The process data are described with a noisy ICA model and conducted robust prewhitening with the covariance matrix of the measurement noises to obtain the whitened data.The NoisyTSICA objective function is then constructed based on the time structure of the whitened data and optimized by the gradient descent algorithm to estimate the mixing matrix and extract the ICs.Furthermore,a monitoring statistic is built to conduct fault detection by using the recursive kurtosis estimations of the dominant ICs.Lastly,a contribution plot for the constructed monitoring statistic is established to identify the fault variables by applying sensitivity analysis.A simple example of estimating the mixing matrix and a process monitoring example of the continuous stirred tank reactor are used to demonstrate the effectiveness of the proposed monitoring method.

2.Fast ICA-based Monitoring Method

The conventional monitoring methods based on ICA usually adopt the noise-free ICA model as follows

where x=[x1,x2,…,xm]T∈Rm×1denotes the vector of the zero-mean measured variables,s=[s1,s2,…,sm]T∈ Rm×1is the vector of the zeromean ICs and A∈Rm×mis the unknown mixing matrix.

ICA tries to estimate both A and s only from x.Equivalently,the objective of ICA can be defined as follows:to find a de-mixing matrix W∈Rm×mwhich can make the elements of the estimated vector ?∈Rm×1given by

be as independent of each other as possible.

Usually,the measured variables need to be whitened firstly.The whitened variables z∈Rm×1are obtained by

where Q=diag(λ1,λ2,…,λm)?1/2[β1,β2,…,βm]T∈ Rm×mdenotes the whitening matrix,λi,i=1,2,…,m are the eigenvalues of the measured variables' covariance matrix Cx=E(xxT)∈Rm×mand satisfy the condition λ1≥ λ2≥ … ≥ λm,diag(λ1,λ2,…,λm)denotes the m × m diagonal matrix with λ1,λ2,…,λmas its diagonal elements,βi,i=1,2,…,m are the eigenvectors of Cxand E(?)denotes the expectation operator.

Then,the whitened variables satisfy the condition

where Im∈Rm×mis the identity matrix.

For the mathematical convenience,all the ICs can be assumed to have the unit variance without loss of generality.Then Eq.(2)is reformulated as

where W=UQ,U=[u1,u2,…,um]T∈Rm×mis an orthogonal matrix due to the reason that E(??T)=E(UzzTUT)=U E(zzT)UT=UUT=Im.

Thus,the problem of estimating the ordinary matrix W is converted to a simpler problem of estimating the orthogonal matrix U.To calculate U,the FastICA algorithm based on the maximum negentropy criterion[22]is widely used in the conventional ICA-based monitoring methods[13-18].The optimization objective of FastICA is defined as follows

where uT∈ R1×mis a row vector of the orthogonal matrix U,υ is a Gaussian variable with zero mean and unit variance,and G(?)is a nonquadratic function and can be chosen as G(uTz)=? exp(?(uTz)2/2).The specific details of FastICA can be found in reference[22].

Once the orthogonal matrix U is obtained,the ICs are estimated by Eq.(5)and arranged in the descending order according to their non-Gaussian degrees measured by the negentropy statistic[18].The row vectors of U are also ordered correspondingly.The mixing matrix is then estimated by ? =(UQ)?1.To conduct fault detection,the two monitoring statistics are defined as follows[9,11,15-18,24]

where x(t)denotes the sample value of x at the sample time t,Udis composed of the first d row vectors of U,?d∈ Rd×1is the vector of the d dominant ICs,d is the number of the dominant ICs,I2(t)is used to monitor the systematic part of the process variation and SPE(t)is used to monitor the residual part of the process variation.

Once a fault is detected by Eqs.(7)-(8),the variable contributions of x(t)to I2(t)and SPE(t)are calculated to identify the fault variables by using the following equations,respectively[15]

3.NoisyTSICA-based Monitoring Method

It is a well accepted fact that the measured data are usually contaminated by the measurement noises with different intensities.Under the negative influences of the measurement noises,the feature extracted by the conventional FastICA-based monitoring method cannot reliably reveal the process information,which may result in the unsatisfactory monitoring performance.In this section,a NoisyTSICA algorithm,which can explicitly consider the measurement noises,is developed,and then a NoisyTSICA-based monitoring method to improve the process monitoring performance is proposed.

3.1.NoisyTSICA algorithm

Different from the noise-free ICA model of Eq.(1),the following noisy ICA model is adopted to develop the NoisyTSICA algorithm

where x=[x1,x2,…,xm]T∈ Rm×1denote the measured variables corrupted by the measurement noise variables ε=[ε1,ε2,…,εm]T∈ Rm×1.The noise variables ε and the ICs s satisfy the conditions:(1)the elements of s are zero-mean and mutually independent;(2)the elements of ε are zero-mean and uncorrelated Gaussian white noise variables;(3)the elements ofεand the elements of s are mutually independent.

The objective of the NoisyTSICA algorithm is to estimate the mixing matrix A and extract the ICs from x.Denote Ω as the covariance matrix of the noise variables.Then,the measured variables' covariance matrix Cxcan be written as

The robust prewhitening on the measured variables to obtain the whitened variables y∈Rm×1is conducted by the following equation

Substituting Eq.(11)into Eq.(13),the following relationship can be obtained

where VT=(Cx?Ω)?1/2A=[v1,v2,…,vm]∈Rm×mis the orthogonal matrix because of the following reason

Based on Eq.(14),the weighted sum of the whitened variables'different time-delayed covariance matrices can be calculated as

where γij,i=1,2,…,lag ? 1,j=i+1,i+2,…,lag are the weighted coefficients,lag is the predefined maximum time lags and~Cyis called the time structure of the whitened variables y in this paper.

Combining Eq.(14)with Eq.(16),the following expression can be derived based on the conditions that ε and s satisfy

where Λs∈ Rm×mdenotes the diagonal matrix withk=1,2,…,m as its diagonal elements,andis called the time structure of the k th IC in this paper.

From Eq.(17),the following result can be obtained by using the orthogonal matrices VTand V

Based on Eq.(18),the objective function of the NoisyTSICA can be constructed as follows

It is obvious that the i th optimal solution vimakes Eq.(19)be equal to the i th minimum value among the diagonal elements of Λs.Thus,Eq.(19)can make use of the time structure of the whitened variables y to seek the column vectors of the orthogonal matrix VTone by one.

The gradient descent algorithm can be adopted to optimize Eq.(19)as

where vi(k)denotes the value of viat the k th iteration and θ is the step size parameter.

Upon obtaining the matrix VT,the mixing matrix can then be estimated as

Thus,the NoisyTSICA algorithm is developed to estimate the mixing matrix A.The procedure of the NoisyTSICA algorithm is summarized as follows.

(1)Conduct the robust prewhitening on the measured variables x by Eq.(13).

(2)Calculate the time structure of the whitened variables y by Eq.(16).

(3)Let i=1.

(4)Set an initial value with the unit norm for vi.

(5)Update viby using Eq.(20).

(6)Do the orthogonalization vi=vjand normalization vi=vi/‖vi‖2.

(7)If vidoes not converge,go back to 5);else,save the vector vi.

(8)Let i=i+1.If i≤m,go back to 4);else,proceed to the next step.

(9)Estimate the mixing matrix by Eq.(21).

After obtaining the estimation of the mixing matrix,its inverse matrixis substituted into Eq.(11)to extract the ICs as

3.2.Fault detection strategy

From Eq.(22),it is found that the extracted ICs~s are contaminated with the Gaussian noise variables.If the monitoring statistic is simply built up based onaccording to Eq.(7),the Gaussian noise variables~ε will also be included in the monitoring statistic and still impose the adverse effects on the process monitoring performance.

Based on the above analysis,a fourth-order statistic is introduced,which is the kurtosis statistic to eliminate the influences of the Gaussian noise variables by taking advantage of the fact that the kurtosis statistic of a zero-mean Gaussian variable is zero[25].The kurtosis of the ith extracted ICcan be calculated as

From Eq.(23),it can be easily seen that the kurtosis of the extracted ICis the same as the kurtosis of the real IC siand the effect of the Gaussian noise variableon the kurtosis statistic is thus removed.In light of this result,the kurtosis statistics of the extracted ICs are used to construct a monitoring statistic for fault detection.Thus,the following recursive expression[26]is firstly taken to estimate the kurtosis of the ith extracted ICat the sample time t

Then, the extracted ICs are arranged in the descending order according to their non-Gaussian degrees which can be measured by the kurtosis absolute values of the extracted ICs in Eq.(23)and the first c extracted ICs are chosen as the dominant ICs.The row vectors of the matrixare then arranged according to the same order.Without loss of generality,are denoted as the dominant ICs.Finally,in order to conduct process fault detection,a monitoring statistic can be constructed as follows by using the recursive kurtosis estimations of the c dominant ICs

After the construction of the monitoring statistic D2,the confidence limit needs to be determined to judge whether the process is in control or not.Theδconfidence limit for the built monitoring statistic can be determined as follows [27]. Split the normal measured data into two parts:the training data with N1samples and the validating data with N2samples.Based on the training data,the mixing matrix,the ICs and the covariance matrix Φ of kurt(|t)are estimated.Based on the validating data,the monitoring statistic valuesare calculated.Round N2(1 ? δ)towards the nearest integer r.The rth highest value ofis adopted as the confidence limit of D2statistic.

The fault detection strategy based on the developed NoisyTSICA algorithm can be divided into two stages:the offline modeling stage and the online monitoring stage.The detailed steps can be demonstrated as follows.

Stage I:the offline modeling stage.

(1)Divide the normal measured data into two parts: the training data and the validating data.

(2)Based on the training data,estimate the mixing matrix and extract the ICs by the NoisyTSICA algorithm.Calculate the covariance matrix Φ of kurt(|t).

(3)Based on the validating data,extract the ICs by Eq.(22),calculate the monitoring statistic values by Eq.(25)and determine the δ confidence limit.

Stage II:the online monitoring stage.

(1)Obtain the current measured data,extract the current ICs by Eq.(22)and calculate the current monitoring statistic value by Eq.(25).

(2)Determine whether the current monitoring statistic value exceeds the confidence limit,and give an alarm if some fault is detected.

In the above procedure,the training data,validating data and current measured data are all normalized with the means and variances of the measured variables in the training data.The recursive kurtosis estimations of the dominant ICs respectively in the training data,validating data and current measured data are all normalized with the means and variances of the recursive kurtosis estimations of the dominant ICs in the training data.

3.3.New contribution plot for fault identification

Once a fault is detected by the constructed monitoring statistic,the next important step of process monitoring is the fault identification,which is to identify the fault variables for locating the root cause of the fault.However,the D2monitoring chart can only indicate the deviation from the normal operating conditions and lacks the ability of pointing out the fault variables causing the process out of control.

In the present study,the contribution plot methods display their effectiveness in identifying the fault variables.The contribution of each process variable to the monitoring statistic can be calculated and the variable with the greatest contribution value usually indicates the fault source.In this section,based on the sensitivity analysis[28]which calculates the rates of change in the system output variables resulting from the small perturbations in the problematic parameters,a new contribution plot for the constructed monitoring statistic D2to identify the fault variables is built up,which defines the measured variables' contributions to the D2statistic as the following vector form

where “°”denotes the Hadamard product of two vectors.

Substituting Eqs.(22)and(25)into Eq.(26),CD2(t)can be further derived as

Based on Eq.(27),the variable contributions of the measured data x(t)for the monitoring statistic value D2(t)can be calculated and the measured variables which have relatively larger contribution values can be identified as the fault variables.

4.Simulation Studies

In this section,the developed NoisyTSICA algorithm and the FastICA algorithm are firstly tested by a simple example of estimating the mixing matrix.Then the proposed NoisyTSICA-based monitoring method and the conventional FastICA-based monitoring method are evaluated by a continuous stirred tank reactor(CSTR)system.

4.1.Simple example of estimating the mixing matrix

Consider three ICs s=[s1,s2,s3]Tthat have the following distributions

where fs=1000 is the sample frequency,t is the sample number ranging from 0 to 1999,the function rem(t,23)denotes t?t/23」*23 andt/23」rounds t/23 to the nearest integer towards zero.These ICs are linearly mixed as Eq.(11)with a randomly generated mixingto get the measured signals.The covariance matrix of the Gaussian noise variables is Ω=diag(0.05,0.10,0.15).The NoisyTSICA and FastICA algorithms are respectively adopted to estimate the mixing matrix from the measured signals corrupted with the Gaussian noises.For NoisyTSICA,the parameter lag in Eq.(16)is set to 100 and all the weighted coefficients are simply set to 1.For comprehensive comparison of the two algorithm performances,one hundred independent experiments are conducted and the Gaussian noises are generated randomly in each experiment.The following performance index of the permutation error(PE)[29]which is widely used for evaluating ICA algorithms is adopted for quantitative comparison

where Gijis the element at the position of the i th row and the j th column of the matrix G.For the FastICA algorithm,the estimation of the mixing matrix is ? =(UQ)?1and thus G can be calculated as G= ??1A=UQA.For the NoisyTSICA algorithm,the estimation of the mixing matrix is ? =(Cx? Ω)1/2VTand thus G can be calculated as.A smaller PE suggests that the corresponding algorithm has a higher accuracy for estimating the mixing matrix.

The PE values of the two algorithms in the 100 independent experiments are shown in Fig.1.Itis easily observed that most of the PE values of the NoisyTSICA algorithm are much lower than the PE values of the FastICA algorithm.Table 1 lists the maximum,minimum and average PE values of the two algorithms.From the Table 1,it is obvious that all the three statistical indexes for the PE values of the NoisyTSICA algorithm are lower than those of the FastICA algorithm.From the above comparisons,it can be concluded that the developed NoisyTSICA can effectively attenuate the adverse effects of the noises on the mixing matrix estimation and has better estimating performance for the mixing matrix and the ICs than FastICA.Thus,it is reasonable to develop the NoisyTSICA-based monitoring method to conduct the fault detection and identification.

Fig.1.PE values comparison of the two algorithms.

Table 1 The maximum,minimum and average PE values of the two algorithms

4.2.Process monitoring in the CSTR system

The CSTR system is widely applied for testing the performance of process monitoring methods[30,31].A schematic diagram of the CSTR system with the cascade control system is shown in Fig.2.In the CSTR system,the reactant A flows into the reactor and then a first order irreversible reaction A→B happens.The product B flows out as an outlet stream.Heat from the exothermic reaction is removed by the cooling flow of the jacket.The liquid level of the reactor is controlled by manipulating the outlet flow.The temperature of the reactor is controlled by manipulating the coolant flow.

Based on the mass,energy and component balances,the dynamic model of this CSTR system can be constructed as

Fig.2.CSTR system with cascade control.

The model parameters in Eqs.(30)-(33)are described in Table 2.The more specific details about the CSTR system can be found in reference by Johannesmeyer et al.[31].

The process data of normal operating condition and faulty conditions are generated by simulating the CSTR process.Ten processvariables are measured and the Gaussian noises are added to all the measurements in the simulation procedure.Table 3 lists the measured variables,the variances of the corresponding noise variables and the noise intensity which is simply defined as the ratio percentage of thenoise variable's variance over the measured variable's variance.The simulation data are sampled every 2 s.The normal operating data with 2000 samples are recorded and divided into two parts:the training data with N1=1000 samples for estimating the mixing matrix and the validating data with N2=1000 samples for determining the control limit.Nine fault patterns are simulated and the measured data of each fault pattern are recorded with 900 samples.For each fault pattern,the fault is introduced at the 190th sample.The simulated fault patterns are shown in Table 4,which include the process parameters change and sensor bias.

Table 2 Parameters of the CSTR system

Table 3 Measured variable,the variance of the noise variable and the noise intensity

The conventional FastICA-based monitoring method and the proposed NoisyTSICA-based monitoring method are applied for monitoring the CSTR system and the process monitoring performance of the two methods is compared.The δ=99%confidence limit is adopted as the alarming threshold.The fault detection performance is evaluated from the fault detection rate which is defined as the percentage of the fault samples whose monitoring statistic values exceed the confidence limit in all the fault samples and the fault detection time.In order to decrease the false alarm,a fault is detected only when six consecutive monitoring statistic values exceed the confidence limit and the fault detection time is then defined as the first sample at which the confidence limit is exceeded[32].In the conventional FastICA-based monitoring method,the number of the dominant ICs is selected as d=5 so that the negentropy cumulative sum of the dominant ICs is above 90%of the negentropy cumulative sum of the extracted ICs.In the NoisyTSICA-based monitoring method,the number of the dominant ICs is selected as c=7 so that the cumulative sum of the dominant ICs'kurtosis absolute values is also above 90%of the cumulative sum of all the extracted ICs'kurtosis absolute values.Both the choices are on the basis of the commonly used cumulative percent variance(CPV)criterion[11,33].The parameter settings for Eq.(16)are the same as those in Section 4.1.The learning rate μ in Eq.(24)is set as 0.197 and how to choose an appropriate learning rate is discussed in Appendix A.

The monitoring charts of the fault 6 are firstly illustrated in Figs.3-4 for method comparison.In order to facilitate the comparison,in each monitoring chart,the monitoring statistic values are divided by their corresponding confidence limit so that the confidence limit is equal to 1.The monitoring statistic values are plotted as the solid line and the confidence limit is plotted as the dashed line.From the D2monitoring chart of the NoisyTSICA-based monitoring method in Fig.4,when the sample number is greater than 350,almost all the monitoring statistic values exceed the confidence limit apparently.By contrast,from the two monitoring charts of the FastICA-based monitoring method in Fig.3,only when the sample number is greater than 380 and 400 respectively,can the statistic values exceed their corresponding confidence limits.The monitoring chart comparison shows that the D2monitoring chart of the NoisyTSICA-based monitoring method can give a much earlier and more obvious fault indication for the fault 6 than both the two monitoring charts of the FastICA-based monitoring method.

The monitoring charts of the fault8 are also illustrated in Figs.5-6 to show the effectiveness of the proposed monitoring method.It is obvious that the D2monitoring chart of the NoisyTSICA-based monitoring method in Fig.6 reacts much more quickly and sharply to the occurring fault than both the two monitoring charts of the FastICA-based monitoring method in Fig.5.Specifically,in the D2monitoring chart,almost all the monitoring statistic values exceed the confidence limit after about the 260th sample,whereas in both the two monitoring charts of the FastICA-based monitoring method,most monitoring statistic values from the 260th sample to the 350th sample are still below or around the corresponding confidence limits and thus cannot give a fault indication.The above comparison of the monitoring charts also demonstrates that the NoisyTSICA-based monitoring method achieves much earlier fault detection time and higher fault detection rate than the FastICA-based monitoring method under the fault pattern 8.

Next,the fault detection performance of the FastICA-based monitoring method and NoisyTSICA-based monitoring method on all the nine fault scenarios of the CSTR system is investigated.The fault detection times and fault detection rates of the two monitoring methods are tabulated in Tables 5 and 6,respectively.From Table 5,it can be observed that both the monitoring methods can detect the occurring fault immediately for the step-change fault patterns 1-4.Compared with the stepchange faults,the ramp-change faults are more difficult to detect,which can be confirmed by the fault detection times for the ramp-change fault patterns 5-9.For the challenging problem of detecting the ramp-change faults,the NoisyTSICA-based monitoring method performs much better than the FastICA-based monitoring method,which can be examined from Table 5.Taking the fault 5 as an example,the D2monitoring statistic of the NoisyTSICA-based monitoring method detects the fault at the 289th sample,while the I2and SPE monitoring statistics of the FastICA-based monitoring method only detect the fault at the 321st and the 340th sample,respectively.The comparison results in Table 5 clearly illustrate the superior ability of the NoisyTSICA-based monitoring method in shortening the fault detection delay for the challenging ramp-change fault patterns 5-9.From Table 6,it is obvious that both the monitoring methods achieve 100%or near to 100%fault detection rates for the step-change fault patterns 1-3,which suggests that almost all the fault samples are successfully detected by both the monitoring methods.However,for the step-change fault pattern 4 and the rampchange fault patterns 5-9,the NoisyTSICA-based monitoring method achieves much higher fault detection rates than the FastICA-based monitoring method,reflecting the more excellent and more reliable fault detection ability of the NoisyTSICA-based monitoring method.The average fault detection rates of the two FastICA-based monitoring statistics for the fault patterns 4-9 are 0.7573 and 0.7209,respectively,whereas the average fault detection rate of the NoisyTSICA-based monitoring statistic for the fault patterns 4-9 is up to 0.8178.According to the above analysis,it can be concluded that the fault detection performance of the NoisyTSICA-based monitoring method outperforms that of the FastICA-based monitoring method especially for the challenging ramp-change faults with slow change of the process parameters.

After a fault is detected,fault identification is the next important task to identify the fault variables and help eliminate the influence of fault.The fault1 is firstly taken as an example to test the faultidentification performance of the FastICA-based monitoring method and NoisyTSICA-based monitoring method.From Table 5,it is found that both the monitoring methods detect the fault 1 immediately after the fault occurs at the 190th sample.The fault identification results of the two monitoring methods at the 190th sample are shown in Figs.7-8.From Fig.8,it is indicated that the No.1 process variable has the dominant contribution to the D2statistic and thus the NoisyTSICA-based monitoring method identifies the first process variable as the fault variable.On the other hand,the FastICA-based monitoring method obtains two very different fault identification results from the two contribution plots in Fig.7 for the I2and SPE statistics,respectively.In particular,the contribution plotfor the I2statistic identifies the No.5,and No.7-No.10 process variables as the fault variables,whereas the contribution plot for the SPE statistic identifies the No.1 process variable as the fault variable,which leads to an inconsistent and confusing identification result.Based on the process knowledge described in Tables 3-4,it is known that the fault variable for the fault 1 is the No.1 process variable.Therefore,the NoisyTSICA-based monitoring method can provide a clearer and more definite indication of the fault variable than the FastICA-based monitoring method under the fault pattern 1.

Fig.3.FastICA-based monitoring charts of the fault 6.

Fig.4.NoisyTSICA-based monitoring chart of the fault 6.

Table 4 Fault patterns of the CSTR system

Fig.5.FastICA-based monitoring charts of the fault 8.

Fig.6.NoisyTSICA-based monitoring chart of the fault 8.

Another faultidentification case is the fault3.From Table 5,this fault is also detected quickly at the 190th sample by both the monitoring methods.The identification results of this fault are demonstrated in Figs.9-10.From Fig.9,the contribution plot for the I2statistic identifies the No.9 process variable as the most probable fault variable,whereas the contribution plot for the SPE statistic identifies the No.7 process variable as the most likely fault variable.By contrast,the contribution plot for the D2statistic identifies the No.5 process variable as the one andonly fault variable explicitly.From Tables 3 and 4,it is shown that the fault 3 is caused by the measurement bias of the No.5 process variable.Thus,the NoisyTSICA-based monitoring method gives a meaningfuland right indication of the fault variable,whereas the FastICA-based monitoring method cannot achieve this point.

Table 5 Comparison of the fault detection time(sample number)

Table 6 Comparison of the fault detection rate

Fig.7.Identification results of fault 1 using FastICA-based monitoring method.

Fig.8.Identification result of fault 1 using NoisyTSICA-based monitoring method.

Fig.9.Identification results of the fault 3 using the FastICA-based monitoring method.

5.Conclusions

Fig.10.Identification result of the fault3 using the NoisyTSICA-based monitoring method.

This paper proposed a Noisy TSICA-based monitoring method.The contributions of this paper are in two aspects.Firstly,a NoisyTSICA algorithm is developed,which can explicitly take the effects of the measurement noises into account to estimate the mixing matrix and extract the ICs.Secondly,an effective process monitoring method on the basis of NoisyTSICA is proposed,which includes the fault detection by constructing the D2monitoring statistic based on the recursive kurtosis estimations of the dominant ICs and the fault identification by establishing a contribution plot for the D2monitoring statistic based on the sensitivity analysis.Simulation results on the example of estimating the mixing matrix and the CSTR process monitoring example demonstrate that Noisy TSICA estimates the mixing matrix more accurately than FastICA,and the proposed NoisyTSICA-based monitoring method performs much better than the conventional FastICA-based monitoring method in terms of the fault detection efficiency and fault identification accuracy.

There are two important issues that should be noted.Firstly,the covariance matrix of the measurement noises is assumed to be available based on the process prior knowledge.How to carry out an appropriate estimation of the matrix only from the measured data will be included in our future work.Secondly,the measurement noises are assumed to have Gaussian distributions in this paper.In our future study,we will further focus on developing an effective monitoring method that needs no requirement for the distributions of measurement noises.

Nomenclature

A mixing matrix

? estimation of the mixing matrix

Cxcovariance matrix of the measured variables

CI2contribution plot for the I2monitoring statistic

CSPE contribution plot for the SPE monitoring statistic

CD2contribution plot for the D2monitoring statistic

c number of the dominant ICs in the NoisyTSICA-based monitoring method

D2monitoring statistic builtin the NoisyTSICA-based monitoring method

d number of the dominant ICs in the FastICA-based monitoring method

fssample frequency

G non-quadratic function used in FastICA

Imidentity matrix

I2monitoring statistic of FastICA-based monitoring method

J1objective function of FastICA

J2objective function of NoisyTSICA

kurt(si) kurtosis of the i th IC

kurt(si|t)recursive kurtosis estimation of the i th IC

kurt(sc|t)vector of the recursive kurtosis estimations of the dominant ICs

m number of measured variables

N1number of samples in the training dataset

N2number of the samples in the validating dataset

Q whitening matrix used in FastICA

SPE monitoring statistic of FastICA-based monitoring method

s original ICs

t sample number

U orthogonal matrix calculated by FastICA

uTrow vector of U

V orthogonal matrix calculated by NoisyTSICA

vTrow vector of V

W de-mixing matrix

x measured variables

y whitened variables in NoisyTSICA

z whitened variables in FastICA

δ confidence limit

ε Gaussian noise variables

θ step size

μ learning rate in the recursive kurtosis estimations of the ICs

υ zero-mean and unit-variance Gaussian variable

Ω covariance matrix of Gaussian noise variables

Appendix.The empirical selection of the learning rate μ

Currently,there is no method for determining the optimal learning rate μ.This parameter can only be selected by trial and error.Here,an empirical method to choose an appropriate value for this parameter is developed.The basic idea of this heuristic method can be explained as follows.

In the off-line modeling stage,the measured data under the normal operating condition are adopted to constitute the training dataset and the validating dataset.For the reason that the samples whose monitoring statistic values exceed the δ confidence limit are referred to as fault samples,a value is chosen for the learning rateμto make the monitoring statistic values of the validating samples be below the confidence limit as much as possible.With this strategy,the region of normal operating condition may be better described.

For the constructed monitoring statistic D2,an index is defined as Eq.(A1)to measure the difference between the δ confidence limit and the monitoring statistic valuesfor the validating data

where R1denotes the number of the monitoring statistic values betweenand?σ,R2denotes the number of the monitoring statistic values less than? σ,denotes the δ confidence limit of the D2monitoring statistic and σ is a predefined constant satisfying the condition σ ＜

A smaller value of η indicates that more monitoring statistic values are far below the confidence limit and thus the region of normal operating condition is better depicted.Based on the fact that the index η is directly affected by the learning rate μ,the index η can be minimized to determine an appropriate value for the parameterμ.Then,the following search procedure can be adopted to choose the learning rate μ.

(1)Set the search range for the learning rate from μminto μmaxthat covers all the possible choices of the learning rate.Choose the required δ and δ0values which satisfy the condition δ0＜ δ.

(2)Start with μ = μmin.

(3)Calculate the monitoring statistic valuesfor the validating data.

(4)Determine the δ confidence limitfor the D2monitoring statistic.

If μ = μmin,the δ0confidence limit Dlim,δ02for the D2monitoring statistic is also determined and set σ =Dlim,δ2? Dlim,δ2.

0

(5)Calculate the index η according to Eq.(A1).

(6)Set μ = μ + μmin.If μ ≤ μmax,go back to step 3);otherwise,stop the procedure.

The value in the search range μminto μmaxthat approximately minimizes the index η may be selected as the appropriate learning rate.

Thus,the learning rate for monitoring the CSTR system can be chosen as follows:The search range for μ is set from μmin=0.001 to μmax=0.5 and the values of δ and δ0are set to 0.99 and 0.96,respectively.The index η as the function of the learning rate μ is demonstrated in Fig.A1.It can be observed from Fig.A1 that,a very smallμmay lead to a large number of normal-operating samples whose monitoring statistic values fluctuate around the confidence limit.This may result in a high false alarm rate for on-line fault detection.Indeed,it is confirmed from our simulation experience that a too small μ may result in a large false alarm rate.On the other hand,it is also suggested from our empirical experience that a too large μ may decrease the sensitivity of the algorithm to the occurring fault and thus raises the missed alarm rate.What is more,Fig.A1 points out that a value from 0.197 to 0.237 can be chosen for the learning rate μ.The results of Fig.A1 together with our empirical experience provide a meaningful suggestion that the value of0.197 is an appropriate choice for μ which can offer a good compromise between the false alarm rate and the missed alarm rate.

Fig.A1.Relationship between the index η and the learning rate μ for the CSTR system.