Cointegration Testing Method for Monitoring Nonstationary Processes

Ind. Eng. Chem. Res. 2009, 48, 3533–3543

3533

Cointegration Testing Method for Monitoring Nonstationary Processes Qian Chen,*,† Uwe Kruger,‡ and Andrew Y. T. Leung§ College of Aerospace Engineering, Nanjing UniVersity of Aeronautics and Astronautics, Nanjing, P. R. China, School of Electroncs, Electrical Engineering and Computer Science, Queens UniVersity, Belfast, U.K., and Department of Building and Construction, City UniVersity of Hong Kong, Hong Kong, P. R. China

This paper introduces cointegration testing method for nonstationary process monitoring, which yields a longrun dynamic equilibrium relationship for nonstationary process systems. The process variables are examined, and then a cointegration model of the tested nonstationary variables is identified. The residual sequence of the cointegration model describes the dynamic equilibrium errors of the nonstationary process system and can be further analyzed for condition monitoring and fault detection purposes. The autocorrelated residual sequence is filtered with AR model first, then compensated to keep the fault signatures from being distorted by the filtering process. An application case study to an industrial distillation unit with a nonstatioanry process shows that a tidy cointegration model can describe the dynamic equilibruim state of the unit and correctly detect abnormal behavior of the process. 1. Introduction The monitoring of complex processes in the chemical and general manufacturing industries is of ever increasing importance and attention often results from growing environmental concerns. Monitoring is, however, a difficult task, as most industrial processes are nonstationary in nature. Such behavior may be caused by seasonal changes, processes that involve emptying and filling cycles, throughput changes, the presence of unmeasured disturbances, and operator interventions, etc. In the literature, a variety of paradigms have been proposed for process monitoring, among which are signal-based techniques,1-4 model-based techniques,5-8 rule-based techniques,9-12and more recently knowledge-based techniques.13-16 For complex systems, however, statistically based techniques collectively referred to as multivariate statistical process control (MSPC) have gained attention because of their conceptual simplicity. MSPC is a datadriven approach and mainly based on principal component analysis or partial least-squares or their extensions.17-23 Problems with MSPC technology arise if the recorded data violate the assumptions imposed on the derived monitoring statistics, for example, data correlation,24,25 time varying behavior,26,27 and nonlinear behavior.28,29 A problem that has only sporadically been addressed, however, is that of nonstationary behavior. For example, Wang et al.30 showed how to adapt the mean and variance of process variables and update the confidence limits of MSPC monitoring model for mildly nonstationary processes. This approach, however, may also accommodate incipient fault conditions, which is undesirable. Box et al.31 and del Castillo32 showed that ARIMA models can describe a nonstationary behavior. For process monitoring, Berthouex and Box33 applied such ARIMA models to predict the performance of a wastewater treatment plant. However, if the identified ARIMA models are inverted to produce white noise sequences, the problem of forecast recovery emerges.34,35 Also, the data-differencing process can cause the loss of dynamic information in the data that leads to poor predictive model performance. To the best of the authors knowledge, a reliable * To whom correspondence should be addressed. E-mail: Q.Chen@ nuaa.edu.cn. Tel.: +86-25-84893221. † Nanjing University of Aeronautics and Astronautics. ‡ Queens University. § City University of Hong Kong.

framework for monitoring nonstationary processes has not yet been reported in the literature. This paper proposes the use of the cointegration testing method to construct nonstationary process models that produce stationary residual sequences which can be further exploited in a statistical process control framework. If a nonstationary variable becomes stationary after being differenced, it is said to be an integrated variable. Engle and Granger36 introduced the concept of cointegration of a set of integrated nonstationary variables to describe the relation bewteen the variables. More precisely, they showed that it is possible for a linear combination of a set of integrated nonstationary variables to be stationary if the nonstationary variables are integrated of the same order and share common trends. In this case the variables are said to be cointegrated. Over the last two decades, economists and statisticians have applied cointegration testing methods to analyze volatile financial or economical systems, which are typical candidates to represent nonstationary processes. A linear combination of a set of nonstationary processes is stationary only if those nonstationary processes have a long-term balanced relationship. In recent years, cointegration testing methods were developed to model the relationship between energy consumption and GDP.37-39 This paper shows that the cointegration testing methodology can be developed into a tool for condition monitoring and fault diagnosis of nonstationary dynamic engineering systems, such as those often encountered in the chemical industry. On the basis of a commonly used model for linear dynamic systems including feedback control and the presence of measured and unmeasured disturbances, it is shown that the assumptions imposed on economic systems can also be imposed on nonstationary engineering systems. This follows from two fundamental facts: (i) the nonstationary variables of a dynamic process system must be correlated to each other and (ii) they are governed by internal physical and/or chemical mechanisms. The cross-correlation shows deterministic and stochastic common trends, and the physical and chemical mechanisms must maintain a long-run dynamic equilibrium relation between the nonstationary variables to carry the system through the designed operating condition. These facts provide the basis for the existence of cointegration relationships, implying that cointegrated variables maintain a dynamic equilibrium state of the system. Since

10.1021/ie801611s CCC: $40.75  2009 American Chemical Society Published on Web 02/10/2009

3534 Ind. Eng. Chem. Res., Vol. 48, No. 7, 2009

cointegrated nonstationary variables produce a long-term equilibrium, the residuals of the cointegration model are stochastic and stationary, and are, in fact, the dynamic equilibrium errors of the nonstationary system. The application of the cointegration testing method to recorded data from a distillation process which possesses nonstationary process behavior, demonstrates that it can be utilized for process monitoring. Since the cointegration does not require differencing the data, as in the case for the ARIMA model, the undesirable impact of forecast recovery is circumvented. 2. Analysis of Existing Methods The description of nonstationary systems has been studied extensively in Box et al.31 and to a lesser extent in del Castillo.32 Since a nonstationary sequence can be obtained by integrating a stationary sequence, autoregressive integrated moving average (ARIMA) models can represent nonstationary sequences. For process monitoring, Berthouex and Box39 applied such ARIMA models to forecast the performance of a wastewater treatment plant. If an ARIMA model is inverted, the nonstationary sequence should be the input to the inverse model and the stationary and serially uncorrelated sequence will be the output. Such an inverse model may be considered for statistical process monitoring. This, however, would also imply that any fault signature would pass through the inverse ARIMA model, and that this inverse model includes a procedure of differencing the nonstationary sequence. Therefore, if a fault condition takes the shape of a step function, the differencing procedure will produce an impulse response function, which implies that the fault signature will decay over a very short period of time. This phenomenon is referred to as forecast recoVery.33,34 The dynamic MSPC approach44–46 can be employed to model a dynamic process system, but the system must be stationary. If the system is nonstationary, the dynamic MSPC approach is not applicable in principle. Wang et al.30 proposed an MSPCbased approach to monitor mildly nonstationary processes; as well as updating the sample mean and variance of each variable, the impact of nonstationary behavior was accommodated by constantly altering the confidence limits or associated univariate monitoring statistics. Although two simulation examples and the analysis of recorded data from a distillation unit showed that fault conditions could be detected, this approach may be insensitive to incipiently developing faults. This deficiency is a consequence of adapting the sample mean and variance, the underlying MSPC monitoring model and the confidence limits. In addition, there is no practically useful criterion to judge whether a process is mildly or strongly nonstationary; the conclusions can be very subjective. Hence, the existing frameworks to monitor nonstationary processes have not been very reliable and practical in the literature so far. This work attempts to show that the cointegration testing method, designed to model nonstationary behavior in economic systems, could be a better choice for its simplicity and efficiency. 3. Preliminaries of Cointegration After a brief outline of the cointegration testing method in subsection 3.1, subsection 3.2 gives a discussion of how to determine a cointegration model. Subsection 3.3 then demonstrates that the cointegration testing method can be applied to nonstationary process monitoring. More information concerning the cointegration testing method is given in the Appendix.

3.1. The Basis of the Cointegration Testing Method. Cointegrated variables share common trends. For simplicity, consider the case of two stochastic variables: z1 ) w1 + φ1

(1a)

z2 ) w2 + φ2

(1b)

Each wi represents a nonstationary sequence (e.g., a random walk) and each φi a stationary sequence (white noise), respectively. If these variables are cointegrated, there exists a set of nonzero parameters β1 and β2 such that β1z1 + β2z2 is a stationary stocastic sequence, that is, β1z1 + β2z2 ) β1(w1 + φ1) + β2(w2 + φ2) ) (β1w1 + β2w2) + (β1φ1 + β2φ2)

(2)

β1z1 + β2z2 is stationary if and only if β1w1 + β2w2 vanishes, that is β1w1 + β2w2 ) 0, or β1w1 ) -β2w2, which implies that w1 and w2 are nonstationary sequences that are equal up to a scaling factor -β2/β1. The parameters βi must be such that they purge the trends from the linear combination. Therefore, the two stochastic variables z1 and z2 must share common trends if cointegrated. The common trends can be deterministic or stochastic. In practice, the time-series will potentially be much more complicated. They can accommodate a constant shift and/or a deterministic trend. For example, the AR(2) models of the sequences are40 z(t) ) F1z(t - 1) + F2z(t - 2) + φ(t)

(3a)

z(t) ) R + F1z(t - 1) + F2z(t - 2) + φ(t)

(3b)

z(t) ) R + µt + F1z(t - 1) + F2z(t - 2) + φ(t)

(3c)

The regressor in eq 3a is the simplest one which only involves stochastic trends. Other two types of trending regressors affect the sampling distribution of the coefficient estimators of F1 and F2. Therefore, the three different regressors will have different coefficient asymptotic distributions which are the key for cointegration tests and for estimators of cointegrating vectors (see Dickey-Fuller test in Appendix for details). If a nonstationary process variable z becomes stationary after differencing d times, it is said to be integrated of order d, denoted as z ≈ I(d).31 Engle and Granger36 considered a set of d-integrated nonstationary variables, zT ) (z1 z2 · · · zN) ∈ RN, that hold a long-run dynamic equilibrium if and only if β1z1 + β2z2 + · · · + βnzN ) ξ

(4)

where ξ ∈ R is the equilibrium residual sequence and βT ) (β1 β2 · · · βN) ∈ RN is a cointegrating Vector. This is the basis of a cointegration model. If the equilibrium is to describe the long-run relationship between the nonstationary variables, then the residual sequence ξ must be stationary, that is, integrated of order 0, ξ ≈ I(0). Therefore, cointegration modeling is a process of identifying coefficient vectors β such that the residual sequence ξ is a stationary sequence. Then the variables are said to be cointegrated variables. Definition. The variable set z is said to be cointegrated of order (d, b), denoted by (z1, z2, · · · , zN) ≈ CI(d,b) if and only if36(1) all variables in z are integrated of the same order d; (2) there exists a parameter vector β such that βTz is integrated of order I(d-b), d g b > 0, and b is the reduction in integration order through the linear combination of the variables in z; and (3) the variables are cointegrated variables and β is referred as

Ind. Eng. Chem. Res., Vol. 48, No. 7, 2009 3535

Figure 1. Reference sequences of nonstationary variables. Table 1. ADF Test Results for Reference Data of Debutaniser Process variable tray 14 temp column overhead pressure tray 2 temp reflux vessel level butane product flow reflux flow reflux flow (difference) percentage C3 in C4 percentage C5 in C4 fresh feed temp butane product flow butane product flow (differ.) reboiler vessel level reboiler Steam Flow percentage C4 in C5 percentage C4 in C5 (differ.) tray 31 temp fresh feed flow reboiler outlet temp reboiler outlet temp (difference)

AR order ADF statistic 2 6 8 9 11 7 (6) 1 1 9 4 (3) 12 11 15 (14) 8 10 19 (18)

-4.5270 -42.2729 -4.0302 -4.9493 -6.6441 -2.0933 -51.3552 -5.0843 -9.5888 -4.8332 -2.3245 -46.4042 -31.2877 -7.5211 -3.0378 -34.0679 -4.7066 -6.0250 -2.6136 -29.4889

Table 2. Johanson H0 Test to Determine Number of Cointegrating Vectors β

test result

H0 hypothesis

trace test statistic

critical value (1%)

stationary stationary stationary stationary stationary nonstationary stationary stationary stationary stationary nonstationary stationary stationary stationary nonstationary stationary stationary stationary nonstationary stationary

re0 re1 re2 re3

124.6432 40.4071 10.9969 2.8182

54.46 35.65 20.04 6.65

a cointegrating vector if b ) d, that is βTz produces a stationary stochastic sequence. In an engineering system, the measured variables must share stochastic and/or deterministic common trends as discussed in the introduction. According to the above definition, the parameters of the cointegrating vector must be determined such that the common trends within the variables are purged. The residual sequence ξ should be a stationary sequence if the cointegration model represents the long-run equilibrium of an engineering system in normal condition. The first and second-order moments of the residuals of the cointegration model remain constant for as long as the process operates under normal condition. In the case that an unmeasured disturbance is introduced into the system for a short period of time, the first and second-order statistics of the residual sequence are also altered but will recover soon after the disturbance ceases. If, however, the dynamic equilibrium is distorted and does not recover, that could indicate that a permanent change

has occurred in the system. The same analysis applies to the presence of a deterministic fault condition, for example, a sensor bias or the performance deterioration of a processing unit. Having given the above discussion, it is reasonable to consider a cointegration model of a set of nonstationary variables as a process monitoring model. Compared with dynamic MSPC and ARIMA models, a cointegration model is conceptually simpler and easy to handle. It should be noted that if z includes a total of N variables, each of which is integrated of order d and the N variables are coinetegrated of order d-b, there may be as many as N - 1 linearly independent cointegrating vectors, which means there are potentially N - 1 cointegration models to describe the longterm dynamic equilibruim state. In theory, these models are equally valid for representing the system equilibrium relationships. However, multiple cointegrating vectors are evaluated by solving a generalized eigenvalue problem numerically as detailed below. As the eigenvectors are associated with the corresponding eigenvalues, the eigenvector corresponding to the largest eigenvalue normally possesses the best precision which establishes most accurate equilibrium relationship. 3.2. Cointegration Modeling Process. The determination of a cointegration model involves the construction of a vector error correction model (VECM), shown in eq A5. The following five steps comprise the procedure of establishing a VECM model: (1) the ordinary least squares (OLS) approach is used to identify an autoregressive (AR) model for each variable in z; (2) a unit root test is carried out using the augmented Dickey-Fuller (ADF)41 method (Appendix A) for each AR model to examine whether the examined variables are integrated of order I(1); (3)


Figure 2. Residual sequences of cointegration model.

Figure 3. Autocorrelation function for residual sequences of cointegration model. Table 3. First Two Cointegrating Vectors Including Offset Term term

reflux flow

butane flow

C4 in C5 concentration

reboiler temp

offset

1st vector 2nd vector

-63.246 -0.005

63.235 19.891

-0.944 2.295

127.804 -267.21

-509.825 1326.563

OLS is used to identify two vector autoregression (VAR) models,42,43 that is, ∆z(t) and z(t) are regressed on ∆z(t-k), k ) 1, 2,..., p, to obtain two residual vectors utilized to form an eigenvalue problem (Appendix B) required for the next step; (4) the eigenvalue problem is solved to obtain the eigenvectors of eq A10a, which are the cointegrating Vectors βj; and, (5) the dominant eigenvector is selected (which is associated with the

largest eigenvalue) to be the cointegrating vector according to the Johanson test (Appendix B). 4. A Process Monitoring Strategy To implement the proposed nonstationary process monitoring scheme, a cointegration model is identified on a training data


Figure 4. Recorded sequences of nonstationary variables describing fresh feed drops.

Figure 5. Residuals sequences of cointegration model describing fresh feed drops.

set and then tested on a validation data set. Those data sets should be sampled when the system is known to be under normal operating conditions. Afterward, the model can be applied to online operating data to monitor the process system by testing the stationarity of the equilibrium residuals with the ADF test method. If the equilibrium residual sequence stays stationary, the process should be considered normal, otherwise abnormal. Also, some of the statistics of the residual sequence can be extracted and verified against the normal criteria, for example,

probability distribution density limits. If the statistic values are consistent with those criteria representing normal operating condition, the process should be regarded to be in statistical control. However, if the values of the statistics violate the criteria, the process could well be out of statistical control, which indicates the presence of a fault condition. In a process system under feedback control, the variables under closed-loop control are assumed to be stationary, while the remaining variables can be nonstationary depending on the system behavior and the input


Figure 6. Hotellings T2 statistics for residuals of cointegration model describing impact of fresh feed drops.

simplify the monitoring charts and therefore reduce the computational work load for the subsequent monitoring process;

Table 4. Parameters of AR Model to Remove Auto Correlation from Residuals lagged term

1

2

3

4

5

value

0.3193

0.2418

0.1780

0.1221

0.0922

variable properties. After identifying a set of nonstationary variables using the ADF test, a process model, that is, a cointegration relationship, can be established between the nonstationary variable set by determining the cointegrating vectors. The steps for establishing a cointegration model are summarized in subsection 4.2 and more details of cointegration testing are given in the Appendix. The main existing statistical-based monitoring techniques include dynamic MSPC and ARIMA models. The dynamic MSPC methods are designed to deal with stationary dynamic process, and would therefore produce false alarms as a result of the nonstationary behavior.30 The ARIMA model can be utilized to describe nonstationary behavior of a multivariate process system,31-33 however, it directly applies differenced data to the model, which, in turn, is associated with the loss of fault signatures and forecast recovery problem.34,35 In comparison with the existing techniques, the cointegration testing method has the following advantages: (A) The cointegration approach does not apply differenced data to the model, and hence the dynamics and fault features embedded in the data remain intact. Although the differencing process is used when determining the matrices of the eigenvalue equation (eq A11) for the cointegrating vectors, the model quality would not be affected for the differenced data are not directly applied to the model. Therefore, the cointegration model will not suffer from the same problems as the ARIMA model like forecast recovery; (B) The cointegration model only produces one residual sequence compared to dynamic MSPC or ARIMA models, which may produce a few residual sequences. This would

(C) A cointegration model is only a simple linear combination of a set of nonstationary variables. It does not involve a complex model structural selection process, which is the case for all other modeling techniques including the ARIMA and dynamic MSPC models; It should be noted, however, that the residual sequences are normally autocorrelated as a result of controller feedback and the presence of unmeasured disturbances which may be autocorrelated. In an MSPC framework, Kruger et al.24 and Xie et al.25 showed that autocorrelation has the potential to significantly affect the number of type I and II errors, which is undesirable. Removing autocorrelation from the residual sequences can therefore circumvent these unwanted effects. However, Dyer et al.47 showed that the application of inverse time-series filters for autocorrelation-removing may alter the magnitude and signature of fault conditions. To address this issue, Lieftucht et al.48 introduced compensation terms to remove the fault information from the original and residual sequences first. The approach then puts the fault features back into the filtered residual sequence. For step-type faults, a consistent estimation of the fault signature can be obtained, whereas general deterministic faults can be approximated by a series of small step-type faults. For the context of a cointegration model, the next subsection discusses the construction of residuals that can be employed in a statistical process control framework. Then the determination of the compensation term for diagnosing general fault conditions is discussed. 4.1. Autocorrelation in Residual Sequences. Assuming that the process variables are auto- and cross-correlated, the following assumption must be imposed on the residual sequence ξ.


Lemma. In the presence of controller feedback and unmeasured disturbances which are autocorrelated, the residual sequence ξ must, in general, also be assumed to be autocorrelated. Proof. Giving a set of N nonstationary sequences of K observations, stored in vectors z1, z2,..., zN, that are autocorrelated within a window (1 τ0j], that is, ∑iK)-τ τ zj(i) zj(i + τ) > 0, 1 e τ e τ0j the autocorrelation function of the residual sequence ξ ) ∑jN) 1 βjzj is given by K-τ

Rξξ(τ) ) K-τ

)

(

N

)

N

i)τ N

j j

j)1 K-τ

+2

∑ ∑

+2

(

2 j j

j

j j

) ( ) ∑ ∑ (∑ i)τ

∑ ∑ β z (i) z (i + τ) i)τ N

)

(

K-τ

∑ ξ(i) ξ(i + τ))∑ ∑ β z (i)∑ β z (i + τ)

j)1 K-τ β2j zj(i) zj(i + τ) j)1 i)τ

N

j)1 N

∑∑∑ i)τ N

βjβk

j)1 k)j+1

N

)

2 j zjzj(τ) + 2

zj(i) zk(i + τ)

i)τ

N

∑β R j)1

) )

βjβkzj(i) zk(i + τ) j)1 k)j+1 N K+τ

N

∑∑

j)1 k)j+1

βjβkRzjzk (5)

where Rzjzj(τ) and Rzjzk(τ) are the autocorrelation and crosscorrelation functions of the nonstationary variables. From eq 5, only if the nonstationary variables are not auto- and crosscorrelated, then ξ is not autocorrelated. However, this assumption cannot hold in general because of the intervention of the feedback controller and the autocorrelated unmeasured disturbances. Therefore, there must exist

∑ a [ξ(t - i) - je(t - i)] i

(9)

i)1

where eˆ(t) is the filtered residual sequence after compensation, m ) t - t0 is the duration in which the step-type fault occurred at instance t0. For a general deterministic fault sequence, Lieftucht et al.48 proposed an approximation of the sequence by a series of step-type faults. More precisely, the fault sequence is segmented into small disjoint windows and if each of the samples in the kth window have been processed, the compensation term, je(t) is set to zero, m is set to one, and the first sample of the (k + 1)th window is investigated. Although successfully applied to a simulation example and an industrial process, this approach can have poor performance for the first few samples in each of the disjoint windows. This can be overcome by employing a moving window approach instead, where the update for the compensation term in eq 9 becomes je(t) ) je(t - 1) +

1 [e(t) - e(t - m0)] m0

(10)

in which m0 is the variable window size. The next section presents an application study on measured data from an industrial distillation unit in which a number of variables show nonstationary behavior. 5. Industrial Application Case Study

Rξξ(τ) > 0 Since the autocorrelation can be described by an AR sequence, an inverse AR filter can be used to produce an uncorrelated sequence. Let e(t) be a stationary and serially uncorrelated sequence, the stationary autocorrelated residual sequence ξ(t) can be determined as follows: 1 ξ(t) ) e(t) A(z-1)

(6)

Here, A(z-1) ) 1 - (a1z-1 + a2z-2 + · · · + aqz-q), where q is the order of the denominator polynomial. Inverting the above AR filter produces e(z-1) by filtering ξ(z-1), that is: e(t) ) A(z-1)ξ(t) or q

e(t) ) ξ(t) -

q

eˆ(t) ) ξ(t) -

∑ a ξ(t - i) i

(7)

i)1

The filtered residual sequence e(t) can now be utilized to construct univariate monitoring statistics for online process monitoring. 4.2. Compensation Terms. Re-examining eq 7 by considering the presence of a deterministic fault sequence f(t) gives which implies that the filtered residual e(t) is the sum of the

original sequence e0(t) and the filtered fault sequence ef(t). For a step-type fault, Lieftucht et al.48 showed that the ef(t) can considerably depart from f(t) and proposed the following compensation scheme to extract an estimate of f(t) from ef(t). The compensation term is defined as

5.1. The Distillation Unit. This distillation process purifies a mix of hydrocarbons, mainly butane (C4), pentane (C5) and impurities of propane (C3) that enter the unit as fresh feed. The separation is achieved with 33 trays within the main distillation tower. The temperatures along this main unit are measured at trays 2 and 31 and between trays 13 and 14. Further temperature measurements include the outlet stream temperature of the reboiler unit entering the bottom of the tower and the fresh feed temperature. The measured flow rates are the fresh feed flow, the reflux flow, the reboiler steam flow and the butane product flow. For safety reasons, the column overhead pressure, the reflux, and reboiler vessel level are monitored. To monitor the product quality of the distillation process, the proportions of the C3 in C4 and C5 in C4 concentrations of the top draw are measured and analyzed. Furthermore, to ensure an economic production, the C4 in C5 concentration of the bottom draw is also recorded. The data were recorded at a sampling interval of 30 s and included approximately 8000 samples that described normal operation (reference data) and a set including around 4500 samples showing the impact of fresh feed fluctuation. The top and bottom concentrations are not under feedback control and the fresh feed flow appears to be varying. Thus, any abnormal drop or increase in the feed or feed composition could result in undesired alterations of these concentrations, which occasionally arise, as demonstrated in this study. Therefore, it is important to monitor the process to detect any unsatisfied product quality caused by the excess fluctuation in fresh feed. The task of this section is to establish a cointegration model for the distillation system to monitor its performance. To build a cointegration model, certain number of variables should be selected first by ADF test as well as a priori knowledge about the distillation unit. The statistic behavior of the filtered and compensated residual sequence of the cointegration model is screened with a monitoring chart to detect any abnormal operating conditions.


5.2. Determining the Cointegration Model. To appropriately rescale the recorded variables, it is a common practice in economics to apply a natural logarithm to each nonstationary sequence. For obtaining a cointegration model, the rescaled variables describing the reference data set, shown in Figure 1, were then analyzed following the steps outlined in section 4.2, which first involved the construction of AR models for each of the 16 recorded variables. Table 1 shows the order of the AR models, estimated using the AIC and SC criteria49,50 for each process variable. Next, the ADF test was applied to determine (i) which of the variables are nonstationary and (ii) if the variables are integrated of order 1, that is, I(1). Table 1 includes these results and also indicates that the variables z6 (reflux flow), z10 (butane product flow), z13 (percentage of C4 in C5) and z16 (reboiler temperature) are nonstationary. The critical value was obtained with respect to a confidence of 99%, which produced a critical value of -3.43; that means that each variable associated with a larger value is nonstationary. After differencing the nonstationary sequences and applying the same test, each of the four differenced variables is stationary confirming that they are I(1) sequences. The VAR and VECM models were then built using OLS to obtain the residual vectors that subsequently formed the basis for constructing the residual covariance matrix, described in eq A10a. The order of the VAR model was determined to be 12 by applying the criteria discussed in ref 51 and ref 52. This was followed by evaluating the eigenvectors and eigenvalues of eq A10a, determining how many cointegrating vectors can be extracted and which one is the most suitable cointegrating vector. Table 2 shows the results of the null hypothesis associated with the Johansen test for a confidence of 99%, highlighting the existence of two cointegrating vectors. Table 3 gives the values of both these cointegrating vectors. As discussed in subsection 4.2, step 5, the first cointegrating vector was selected, and the cointegration model used for process monitoring was therefore: ξ ) -63.2z6 + 63.2z10 - 0.944z13 + 127.8z16 - 509.8 or in a normalized form, ξ ) -z6 + z10 - 0.015z13 + 2.02z16 - 8.06

(11)

The upper plot in Figure 2 shows the residual sequence of the above cointegration models and shows that it possesses a significant autocorrelation, also noticeable from the associated autocorrelation function (upper plot in Figure 3). On the basis of the BIC criterion as proposed in Hannan53 and Pauler,54 the best model structure should be five time lags for the AR model to remove autocorrelation. By the recommendation in Kruger et al.,24 the cost function used for the BIC criterion included the mismatch between an uncorrelated sequence and the estimated autocorrelation function. Table 4 shows the parameters of the identified AR model. The lower plot in Figure 3 confirms that the residual sequence of an inverse AR(5) filter only produces negligible autocorrelation. The lower plot in Figure 2 shows the filtered residual sequence of the cointegration model. 5.3. Analyzing Abnormal Condition Caused by Fresh Feed Drops. Some drops in fresh feed arose just after 1800 and 2900 samples into the second data set. The plant operators counteracted the presence of the first drop by slightly reducing the reflux flow rate; however, the second drop, where the feed reduced from an average value of 25 T/h to 17 T/h over a period of 15 h, remained unnoticed. After a significant increase in impurities within the butane product flow, the plant operators

introduced another reduction in the reflux flow rate to reduce the considerable variation of the material and energy balance, resultant from the second drop. This example therefore shows that it is essential to detect such drops as early as possible to initiate appropriate countermeasures for preventing produce quality from being compromised. Figure 4 shows the recorded sequences of the four nonstationary variables. While the first drop had an insignificant effect upon the operation of the unit, the second prolonged drop produced substantial fluctuations of the material as well as the enthalpy balance, noticeable by the butane product flow (material balance) and the reboiler temperature (energy balance). Assuming that any change in the material and energy balance within the unit leads to an alteration in the variable interrelationships governed by the cointegration model, the top plot in Figure 5 shows that both events could be detected. However, after the removal of autocorrelation, both events could only be sporadically detected. This results from the impact of the inverse AR filter upon the fault signature, discussed in subsection 5.2. Applying the compensation term defined in eq 10 for m0 ) 15 removed the undesired effect and extracted the fault signature (bottom plot in Figure 5) from the filtered residuals of the cointegration model. Comparing the top and bottom plot in Figure 5, it should be noted that the fault signature accurately described both fresh feed drops as well as the disturbance in the material and energy balance of the unit from the second drop. In contrast, the fault signature negligibly departs from zero for the first 1800 samples and between the 2100th and 2900th sample. The above analysis is an explanation of the cause of the changes based on the production diary and domain knowledge, which does not suggest the cointegration model can detect the causes of the changes. To isolate the causes, further diagnosis is needed which will be dealt with in later work. Although the residual sequences in Figure 5 can give some information about the abnormal conditions, a more informative monitoring chart should be constructed. A Hotellings T2 statistic with confidence limits is a simple but adequate option. Figure 6 indicates that the Hotellings T2 statistic of the compensated and filtered residual sequence of the cointegration model could clearly detect the presence and impacts of both drops, while the filtered statistic, as expected, only produced occasional violations that did not reveal the correct fault signature. The bottom plot in Figure 6 indicates that an almost complete removal of the fault signature from the compensated residual sequence could be achieved. The “fault-free” residual sequence e0 can be obtained by e0 ) eˆ - je, where eˆ and je, are determined by eq 9. The violations between 3350 and 4400 emerged as a result of the sharp spikes in the butane product flow (Figure 4). 6. Conclusions This paper has introduced the cointegration paradigm as a condition monitoring strategy for complex nonstationary processes systems. In comparison with existing work, which mainly focuses on the use of ARIMA models, the cointegration testing method can overcome the problem of forecast recovery. On the basis of a linear time-invariant model structure for a process under feedback control, the paper has shown that the variables are cointegrated and that the equilibrium sequences (residuals of the cointegration model) can be exploited in a statistical process control framework. This has been demonstrated by the application of the developed monitoring scheme to the measured data from an industrial distillation unit, for which sequences of fresh feed drops (faults) were recorded.


Since the process variables are autocorrelated and crosscorrelated and the process plant is under closed-loop control, the stationary residuals of the cointegration model are also autocorrelated. To address this issue, autocorrelation filtering and fault compensation measures are performed on the residual sequence of the cointegration model to prevent the production of significant type I and II errors. Moreover, to extract the correct fault signature from the filtered residual sequence of the cointegration model, a recently proposed compensation scheme has been utilized. Applying this monitoring approach to the industrial data, a fresh feed drop could be detected and the impacts of these events upon the material and energy balance of the unit were correctly identified. This work has demonstrated that the cointegration testing method can be a useful methodology for engineering system condition monitoring and fault diagnosis. Acknowledgment The authors gratefully acknowledge that this work is jointly supported by Natural Science Foundation of China (grant number: 10772080) and Research Grant Council of Hong Kong (grant number 1161/05).

R ) 0 to determine the type of the deterministic trends for choosing a right ADF criterion table. Johanson Test. Johansen42 and Juselius43 introduced a cointegration test based on vector auto-regression (VAR), which is applied to a set of nonstationary variables and integrated of order I(1). In general, a cointegration model of N nonstationary variables is given by the following linear combination: (A3) ξ(t) ) βTz(t) Consider the VAR(p) model representing a multivariate AR time-series structure of order p: p

z(t) )

∑ Π z(t - k) + φ(t)

(A4)

k

k)1

where the Πk are parameter matrices of the multivariate AR structure and φ(t) is the vector of white noise components. A vector error correction model (VECM) can be obtained by subtracting z(t - 1) from both side of eq A4: p-1

∆z(t) ) Πz(t - 1) +

∑ Φ ∆z(t - k) + φ(t)

(A5)

k

k)1

where p

Appendix A Unit Root Test. The unit root test is employed in association with the tests including the Dickey-Fuller (DF)55 and Augmented Dickey Fuller (ADF) test.41 The DF or ADF test is the starting point for cointegration modeling process. [ Dickey-Fuller Test. For simplicity consider the nonstationary time series in eq 3a, that is z(t) ) F1z(t - 1) + F2z(t 2) + φ(t) ) (F1 + F2)z(t - 1) - F2∆z(t - 1) + φ(t), where φ(t) is assumed to be a normally distributed white noise sequence, φ(t) ≈ Ν{0,σ2}. The parameters F1 and F2 can be estimated with an ordinary least square (OLS) method. Assuming that |F2| < 1 the hypothesis H0: F1 + F2 ) 1 can be tested. If H0 is accepted then z(t) has a unit root and eq A6a becomes ∆z(t) ) -F2∆z(t - 1) + φ(t)

(A1)

Since |F2| < 1, ∆z(t) is a stationary zero-mean sequence, this implies that z(t) is an I(1) nonstationary sequence. If H0 is rejected, one carries on testing the AR(2) model for z1(t) ) ∆z(t). If H0 is now accepted z(t) is an I(2) nonstationary sequence. If H0 is rejected again, one repeats the test for d ) 3, etc. until H0 is accepted. To test the hypothesis H0, that is F1 + F2 is significantly close to 1, the t-statistic, tF ) (Fˆ 1 + Fˆ 2)/(se(Fˆ 1 + Fˆ 2)) can be employed, which follows the Dickey-Fuller distribution. Note that Fˆ 1 and Fˆ 2 are the OLS estimates of F1 and F2, respectively, and se( · ) is the standard deviation of the estimated parameters. On the basis of the above equation, the t-statistic is evaluated and compared against 1%, 5%, or 10% of the level values in the DF table. Three different DF tables should be used for the three models of time series with different determinitic common trends. [ Augmented Dickey-Fuller Test. In the DF test, because it cannot be guaranteed that φt is a white noise, the parameter estimation is biased. In 1979 and 1980 Dickey and Fuller extended the DF test method to the so-called augmented DF (ADF) test for the more general case. Consider the presence of a deterministic trend in the AR(2) model of eq 3c, ∆zt ) R + µt - F2∆zt-1 + (F1 + F2)zt-1 + φt

(A2)

after estimating the parameters of the model, the null hypothesis of the ADF test inlcudes H0: F1 + F2 ) 1, H0: µ ) 0 and H0:

Π ) -In +

∑Π

k

(A6a)

k)1

p

Φk ) -

∑

Π j,

k ) 1, 2, · · · , p - 1

(A6b)

j)k+1

Π can be decomposed into a product of three matrices, that is, Π(1) ) U(1) M(1) V(1) in which M(1) is a diagonal matrix and has rank r, Π ) -Π(1) also has rank r. Let β denote an n × r matrix whose columns form a basis for the row space of Π, so that every row of Π can be written as a linear combination of the rows of βT; thus one can write Π ) δβT, where δ is an n × r matrix with full column rank. Equation A5 then becomes p-1

∆z(t) ) δβTz(t - 1) +

∑ Φ ∆z(t - k) + φ(t) k

(A7)

k)1

The residual vector ξ(t - 1) can now be obtained as follows: p-1

ξ(t - 1) ) (δTδ)-1δT(∆z(t) -

∑ Φ ∆z(t - k) - φ(t)) k

k)1

(A8) Equation A8 indicates that the elements of ξ(t - 1) are integrated of I(0), that is, stationary sequences. Thus, the linear combination of eq A3 is I(0), and the column vectors of β are cointegrating Vectors. Assuming that the VAR model has k < n unit roots, the rank of Π is r ) n - k. The null hypotheses H0 is therefore: H0 : rank(Π) e r with Π ) δβT T

(A9)

β z(t) could be interpreted as the “equilibrium” of a dynamic system, ξ(t - 1) as the vector of “equilibrium errors” and eq A7 describes the self-correcting mechanism of the system from the “non-equilibrium” state. The estimates for δ, β, and Π are subject to some restrictions so that δ and β can be separately identified. The matrix Π can be found by a nonlinear least-squares regression of ∆z(i). The parameter matrices Φ1, Φ2,..., Φp-1 can then be obtained using


OLS regression of ∆z(t) - δβTz(t - 1) on the lagged differences. To this end, one carries out the regression of ∆z(t) on the lagged differences giving the residuals eˆ0, and then the regression of z(t - 1) on the lagged differences giving the residuals eˆ1. It can be shown that another regression model eˆ0t ˆ eˆ1t + εt can be obtained based on the theorem of uniqueness )Π of vector orthogonal projection. Johansen56 showed that testing ˆ is equivalent to testing for the number of for the rank of Π canonical correlations between eˆ0 and eˆ1 that are different from zero using the following test statistics: n

λtr(r) ) -T

∑ ln(1 - λˆ ) i

(A10a)

i)r+1

λmax(r, r + 1) ) -T ln(1 - λˆ r+1)

(A10b)

- 1 S01 with where the λˆ ks are the eigenvalues of the matrix S10S00 respect to the matrix S11, ordered in descending order 1 > λˆ 1 > · · · > λˆ n > 0, where T is the sample size of the variables, and Sij ) 1/T∑ eîteˆTjt, i,j ) 0,1. The null hypothesis of the trace statistic in eq A10a tests that the number of cointegrating vectors is less than or equal to r. Note that the further the eigenvalues are from zero the more negative ln(1 - λˆ i) becomes. Likewise, the maximim eigenvalue in eq A10b tests a null hypothesis of r cointegrating vectors against the specific alternative of r + 1. As above, if λˆ r + 1 is close to zero, the statistic will be insignificant. Further, if the null hypothesis is accepted, the r cointegrating vectors contained in β* can be estimated as the first r columns ˆ ) [vˆ1 vˆ2 · · · vˆN]: of matrix V

ˆ i ) 0, (λiS11 - S10S-1 00 S01)v

i ) 1, 2, · · · , n

(A11)

subject to vîTS11vˆj ) δij. If β has been obtained, the estimates of δ, Πk, and φ(t) in eq A5) can be computed by inserting βˆ in their corresponding OLS formulation. It should be noted that if only the cointegration model is required rather than the VECM model, then the eigenvectors of the eq A11 complete the identification procedure. Literature Cited (1) Bardon, O.; Sidahmed, M. Early Detection of Leakages in The Exhaust and Discharge Systems of Reciprocating Machines by Vibration Analysis. Mech. Syst. Signal Process. 1996, 8, 551. (2) Chen, Y. D.; Du, R.; Qu, L. S. Fault Features of Large Rotating Machinery and Diagnosis Using Sensor Fusion. J. Sound Vib. 1995, 188, 227. (3) Kim, K.; Parlos, A. G. 2003, Reducing The Impact of False Alarms in Induction Motor Fault Diagnosis. J. Dyn. Syst., Meas. Control, Trans. ASME 2003, 125, 80. (4) Hu, N.; Chen, M.; Wen, X. The Application of Stochastic Resonance Theory for Early Detecting Rub Impact Fault of Rotor Systems. Mech. Syst. Signal Process. 2003, 17, 883. (5) Leyval, L.; Gentil, S.; Feray-Beaumont, S. Modal Based Causal Reasoning for Process Supervision. Automatica 1994, 30, 1295. (6) Mylaraswamy, D.; Kavuri, S.; Venkatasubramanian, V. A Framework for Automated DeVelopment of Causal Models for Fault Diagnosis. Proceedings of the Annual AIChE Meeting, 1994, San Francisco, CA. (7) Frank, P. M. Analytical and Qualitative Model-Based Fault DiagnosissA Survey and Some New Results. Eur. J. Control 1996, 6. (8) Rajaraman, S.; Hahn, J.; Mannan, M. S. A Methodology for Fault Detection, Isolation and Identification for Nonlinear Processes with Parametric Uncertainties. Ind. Eng. Chem. Res. 2004, 43, 6744. (9) Kramer, M. A.; Palowitch, B. L. A Rule-Based Approach to Fault Diagnosis Using the Signed Directed Graph. AIChE J. 1987, 33, 1067. (10) Iserman, R. Fault Diagnosis of Machines via Parameter Estimation and Knowledge Processing: A Tutorial Paper: Fault Detection, Supervision and Safety of Technical Processes. Automatica 1993, 29, 815.

(11) Shin, R.; Lee, L. S. Use of Fuzzy Cause-Effect Digraph for Resolution Fault Diagnosis of Process Plants -I. Fuzzy Case-Effect Digraph. Ind. Eng. Chem. Res. 1995, 34, 1688. (12) Upadhyaya, B. R.; Zhao, K.; Lu, B. Fault Monitoring of Nuclear Power Plant Sensor and Field Devices. Prog. Nucl. Energy 2003, 43, 337. (13) Lehane, M.; Dube, F.; Halasz, M.; Orchard, R.; Wylie, R.; Zaluski, M. Integrated Diagnosis System (ids) for Aircraft Fleet Maintenance. Proceedings of the 5th National Conference on Artificial Intelligence (AAAI 98), 1998, Menlo Park, CA. (14) Ming, R.; Haibin, Y.; Heming, Y. Integrated Distribution Intelligent Systems Architecture for Incidents Monitoring and Diagnosis. Comput. Ind. 1998, 37, 143. (15) Qing, Z.; Zhihan, X. Design of A Novel Knowledge-Based Fault Detection and Isolation Scheme. IEEE Trans. Syst., Man Cybernet., Part B 2004, 34, 1089. (16) Shing, C. T.; Chee, P. L. Application of An Adaptive Neural Network with Symbolic Rule Extraction to Fault Detection and Diagnosis in A Power Generation Plant. IEEE Trans. Energy ConVers. 2004, 19, 369. (17) MacGregor, J. F.; Kourti, T. Statistical Process Control of Multivariate Processes. Control Eng. Pract. 1995, 3, 403. (18) Nomikos, P.; MacGregor, J. F. Multivariate SPC Charts for Monitoring Batch Processes. Technometrics 1995, 37, 41. (19) Wise, B. M.; Gallagher, N. B. The Process Chemometrics Approach to Process Monitoring and Fault Detection. J. Process Control 1996, 6, 329. (20) Dunia, R.; Qin, S. J. Subspace Approach to Multidimensional Fault Identification and Reconstruction. AIChE J. 1998, 44, 1813. ¨ nday, C. Statistical Process and Controller Performance (21) Çinar, A.; U Monitoring: A Tutorial on Current Methods and Future Directions. Proc. Am. Control Conf. 1999, 2625. (22) Kruger, U.; Chen, Q.; Sandoz, D. J.; McFarlane, R. C. Extended PLS Approach for Enhanced Condition Monitoring of Industrial Processes. AIChE J. 2001, 47, 2076. (23) Kourti, T. Application of Latent Variable Methods to Process Control and Multivariate Statistical Process Control in Industry. Int.l J. Adapt. Control Signal Process. 2005, 19, 213. (24) Kruger, U.; Zhou, Y.; Irwin, G. W. Improved Principal Component Monitoring of Large-Scale Processes. J. Process Control 2005, 14, 879. (25) Xie, L.; Kruger, U.; Lieftucht, D.; Littler, T.; Chen, Q.; Wang, S.Q. Statistical Monitoring of Dynamic Multivariate Processes -Part 1. Modelling Autocorrelation and Cross-Correlation. Ind. Eng. Chem. Res. 2006, 45, 1659. (26) Li, W.; Yue, H.; Valle-Cervantes, S.; Qin, S. J. Recursive PCA for Adaptive Process Monitoring. J. Process Control 2000, 10, 471. (27) Wang, X.; Kruger, U.; Irwin, G. W. Process Monitoring Approach Using Fast Moving Window PCA. Ind. Eng. Chem. Res. 2005, 44, 5691. (28) Jia, F.; Martin, E. B.; Morris, A. J. Non-Linear Principal Component Analysis with Application to Process Fault Detection. Comput. Chem. Eng. 1998, 22, 851. (29) Shao, R.; Jia, F.; Martin, E. B.; Morris, A. J. Wavelets and NonLinear Principal Components Analysis for Process Monitoring. Control Eng. Pract. 1999, 7, 865. (30) Wang, X.; Kruger, U.; Lennox, B. Recursive Partial Least Squares Algorithms for Monitoring Complex Industrial Processes. Control Eng. Pract. 2003, 11, 613. (31) Box, G. E. P.; Jenkins, G. M.; Reinsel, G. C. Time Series Analysis -Forecasting and Control, 3rd ed; Prentice-Hall: Englewood Cliffs, NJ, 1994. (32) Castillo, E. Del Statistical Process Adjustment for Quality Control; John Wiley & Sons: New York, 2002. (33) Berthouex, P. M.; Box, G. E. Time Series Models for Forecasting Wastewater Treatment Plant Performance. Water Res. 1996, 30, 1865. (34) Superville, D. W.; Adams, B. M. An Evaluation of Forecast-Based Quality Control Schemes. Commun. Stat.: Simul. Comput. 1994, 23, 645. (35) Apley, D. W.; Jianjun, S. The GLRT for Statistical Process Control of Autocorrelated Processes. IEEE Trans. 1999, 31, 1123. (36) Engle, R. F.; Granger, C. W. J. Cointegration and Error-Correction: Representation, Estimation and Testing. Econometrica 1987, 55, 251. (37) Ramanathan, R. Short- and Long-Run Elasticities of Gasoline Demand in India: An Empirical Analysis Using Cointegration Techniques. Energy Econom. 1999, 21, 321. (38) Han, Z.-Y.; Fan, Y.; Wei, Y.-M. Study on The Cointegration and Causality Between GDP and Energy Consumption in China. Int. J. Global Energy 2004, 22, 225. (39) Stern, D. I. Multivariate Cointegration Analysis of The Role of Energy in The US Macroeconomy. Energy Econom. 2000, 22, 267. (40) Waterson, M. W. Vector Auto-Regressions and Cointegrations In Handbook of Economics, version 4; Engle, R. F., McFadden, D. L., Eds.; Elsevier Science, BV: New York, 1994.

Ind. Eng. Chem. Res., Vol. 48, No. 7, 2009 3543 (41) Dickey, D. A.; Fuller, W. A. Distributions of The Estimators for Auto-Regressive Time Series with A Unit Root. J. Am. Stat. Assoc. 1979, 74, 427. (42) Johansen, S. Statistical Analysis of Cointegrating Vectors. J. Econom., Dyn. Control 1988, 12, 231. (43) Johansen, S.; Juselius, K. Maximum Likelihood Estimation and Inference on Cointegration with Applications to The Demand for Money. Oxford Bull. Econom. Stat. 1990, 52, 21. (44) MacGregor, J. F.; Marlin, T. E.; Kresta, J.; Skagerberg, B. MultiVariate Statistical Methods in Process Analysis and Control; AIChE Symposium of the 4th International Conference on Chemical Process Control; AIChE: New York, 1991; Publications P-67. (45) Treasure, R. J.; Kruger, U.; Cooper, J. E. Dynamic Multivariate Statistical Process Control Using Subspace Identification. J. Process Control 2004, 14, 279. (46) Li, P.; Treasure, R. J.; Kruger, U. Dynamic Principal Component Analysis Using Subspace Model Identification; Lecture Notes in Computer Science; Springer-Verlag: New York; 2005; Vol. 3644, pp 727. (47) Dyer, J. N.; Adams, B., M.; Conerly, M. D. A Simulation Study of The Impact of Forecast Recovery for Control Charts Applied to ARMA Processes. J. Mod. Appl. Stat. Methods 2002, 1, 343. (48) Lieftucht, D.; Kruger, U.; Xie, L.; Littler, T.; Chen, Q.; Wang, S.Q. Statistical Monitoring of Dynamic Multivariate Processes — Part 2. Identifying Fault Magnitude and Signature. Ind. Eng. Chem. Res. 2006, 45, 1677.

(49) Akaike, H. A New Look at The Statistical Model Identification. IEEE Trans. Automat. Control 1974, 19, 716. (50) Schwarz, G. Estimating The Dimension of A Model. Ann. Stat. 1978, 6, 461. (51) Lutkepohl, H. Introduction to Multiple Time Series Analysis; Springer-Verlag: Berlin, 1991. (52) Paulsen, J. Order Determination of Multivariate Autoregressive Time Series with Unit Roots. J. Time Ser. Anal. 1984, 5, 115. (53) Hannan, J. The Estimation of The Order of An ARMA Process. Ann. Stat. 1980, 8, 1071. (54) Pauler, D. K. The Schwarz Criterion and Related Methods for Normal Linear Models. Biometrika 1998, 85, 13. (55) Dickey, D. A.; Fuller, W. A. Likelihood Ratio Statistics for Autoregressive Time Series with a Unit Root. Econometrica 1981, 49, 1057. (56) Johnsen, S. Likelihood-Based Inference in Cointegrated Vector AutoregressiVe Models; Oxford University Press: Oxford, 1995.

ReceiVed for reView October 23, 2008 ReVised manuscript receiVed December 31, 2008 Accepted January 9, 2009 IE801611S

Cointegration Testing Method for Monitoring Nonstationary Processes

Recommend Documents